diff options
| author | Ken Kellner <ken@kenkellner.com> | 2019-09-23 21:14:49 -0400 |
|---|---|---|
| committer | Ken Kellner <ken@kenkellner.com> | 2019-09-23 21:14:49 -0400 |
| commit | eb36fb97e2e993e7e6b9895c9d48ed85e4ea7e31 (patch) | |
| tree | 328a44ef637937ec9092a12b2f156812bc55bbb6 | |
| parent | 83ed0eb41d6eabf3a4c2adab2e01c0f72366d620 (diff) | |
Handle missing values in nmix
| -rw-r--r-- | examples/example_nmix.jl | 2 | ||||
| -rw-r--r-- | src/nmix.jl | 6 | ||||
| -rw-r--r-- | test/test_nmix.jl | 14 |
3 files changed, 14 insertions, 8 deletions
diff --git a/examples/example_nmix.jl b/examples/example_nmix.jl index e4fb2f9..1a31209 100644 --- a/examples/example_nmix.jl +++ b/examples/example_nmix.jl @@ -24,7 +24,7 @@ fit_all = nmix(allsub(λ_formula), allsub(p_formula), umd); #Model selection table fit_all -#Missing values (broken right now) +#Missing values yna = Array{Union{Int,Missing}}(deepcopy(umd.y)) yna[1,:] = fill(missing, 5) yna[2,1] = missing diff --git a/src/nmix.jl b/src/nmix.jl index 087bf2d..33551ab 100644 --- a/src/nmix.jl +++ b/src/nmix.jl @@ -25,7 +25,7 @@ function nmix(λ_formula::FormulaTerm, p_formula::FormulaTerm, y = data.y N, J = size(y) - K = isnothing(K) ? maximum(y) + 20 : K + K = isnothing(K) ? maximum(skipmissing(y)) + 20 : K function loglik(β::Array) @@ -44,7 +44,7 @@ function nmix(λ_formula::FormulaTerm, p_formula::FormulaTerm, fg = zeros(eltype(β),K+1) for k = 0:K - if maximum(y[i,:]) > k + if maximum(skipmissing(y[i,:])) > k continue end fk = pdf(Poisson(λ[i]), k) @@ -168,7 +168,7 @@ function residuals(fit::Nmix) end function chisq(fit::Nmix) - return sum(residuals(fit) .^ 2 ./ fitted(fit)) + return sum(skipmissing(residuals(fit) .^ 2 ./ fitted(fit))) end """ diff --git a/test/test_nmix.jl b/test/test_nmix.jl index 3be70ee..c20fce1 100644 --- a/test/test_nmix.jl +++ b/test/test_nmix.jl @@ -5,6 +5,13 @@ umd1 = simulate(Nmix, @formula(ψ~elev+forest), @formula(p~wind+precip), (10,3), [0,0,0,0,0,0]) fit1 = nmix(@formula(ψ~elev), @formula(p~wind), umd1) + + yna = Array{Union{Int,Missing}}(deepcopy(umd1.y)) + yna[1,:] = fill(missing,3) + yna[2,1] = missing + yna[3,[1,3]] = fill(missing,2) + umd2 = UmData(yna, umd1.site_covs, umd1.obs_covs) + fit2 = nmix(@formula(ψ~elev), @formula(p~wind), umd2) #Test fitted ft = Unmarked.fitted(fit1) @@ -21,14 +28,13 @@ #Test chi-square calc @test isequal(round(Unmarked.chisq(fit1),digits=4), 23.0500) #with missing values - #@test isequal(round(Unmarked.mb_chisq(fit2),digits=4), 3.6675) + @test isequal(round(Unmarked.chisq(fit2),digits=4), 19.3308) #Test parametric bootstrap Random.seed!(123) gf = gof(fit1) @test isequal(length(gf.tstar), 30) @test isequal(mean(gf.t0 .< gf.tstar), 0.7000) - #gf2 = gof(fit2) - #@test isequal(mean(gf2.t0 .< gf2.tstar), 0.86) - + gf2 = gof(fit2) + @test isequal(mean(gf2.t0 .< gf2.tstar), 0.6000) end |