diff options
author | Ken Kellner <ken@kenkellner.com> | 2023-04-05 18:13:13 -0400 |
---|---|---|
committer | Ken Kellner <ken@kenkellner.com> | 2023-04-05 18:13:13 -0400 |
commit | d5caacc1c20d8dbc253dfa8d31c26af3119fa1bf (patch) | |
tree | 3c45b43433ec06dd7ecf55ea78df94963eefd531 | |
parent | 845a8e74b5e5eca078bd9e1f5581c8ae94ee266e (diff) |
Fix multinomPois ranef bug when NAs in removal sampling
-rw-r--r-- | DESCRIPTION | 4 | ||||
-rw-r--r-- | R/ranef.R | 23 | ||||
-rw-r--r-- | tests/testthat/test_multinomPois.R | 11 | ||||
-rw-r--r-- | vignettes/occuMulti.Rmd | 2 |
4 files changed, 28 insertions, 12 deletions
diff --git a/DESCRIPTION b/DESCRIPTION index c8511af..de0ed2f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: unmarked -Version: 1.2.5.9012 -Date: 2023-03-08 +Version: 1.2.5.9013 +Date: 2023-04-05 Type: Package Title: Models for Data from Unmarked Animals Authors@R: c( @@ -249,24 +249,35 @@ setMethod("ranef", "unmarkedFitMPois", lam <- predict(object, type="state")[,1] R <- length(lam) cp <- getP(object) - cp <- cbind(cp, 1-rowSums(cp)) + cp <- cbind(cp, 1-rowSums(cp, na.rm=TRUE)) N <- 0:K - post <- array(0, c(R, K+1, 1)) + post <- array(NA, c(R, K+1, 1)) colnames(post) <- N for(i in 1:R) { f <- dpois(N, lam[i]) g <- rep(1, K+1) - if(any(is.na(y[i,])) | any(is.na(cp[i,]))) + yi <- y[i,] + cpi <- cp[i,] + if(any(is.na(y[i,])) | any(is.na(cp[i,]))){ + # This only handles cases when all NAs are at the end of the count vector + # Not when they are interspersed. I don't know how to handle that situation + if(object@data@piFun == "removalPiFun"){ + warning("NAs in counts and/or covariates for site ",i, ". Keeping only counts before the first NA", call.=FALSE) + notNA <- min(which(is.na(y[i,]))[1], which(is.na(cp[i,]))[1], na.rm=TRUE) - 1 + yi <- yi[c(1:notNA)] + cpi <- cpi[c(1:notNA, length(cpi))] + } else { next + } + } for(k in 1:(K+1)) { - yi <- y[i,] ydot <- N[k] - sum(yi) if(ydot<0) { g[k] <- 0 next } - yi <- c(yi, ydot) - g[k] <- g[k] * dmultinom(yi, size=N[k], prob=cp[i,]) + yik <- c(yi, ydot) + g[k] <- g[k] * dmultinom(yik, size=N[k], prob=cpi) } fudge <- f*g post[i,,1] <- fudge / sum(fudge) diff --git a/tests/testthat/test_multinomPois.R b/tests/testthat/test_multinomPois.R index 6b4b693..1103839 100644 --- a/tests/testthat/test_multinomPois.R +++ b/tests/testthat/test_multinomPois.R @@ -101,9 +101,14 @@ test_that("multinomPois can fit a removal model",{ res <- residuals(m2_C) expect_equal(dim(res), dim(umf1@y)) - expect_warning(r <- ranef(m2_C)) - expect_equal(dim(r@post), c(3,56,1)) - expect_equal(bup(r), c(10.794,0.000,2.655), tol=1e-4) + expect_warning(r <- ranef(m2_C, K=50)) + expect_equal(dim(r@post), c(3,51,1)) + expect_equal(bup(r), c(10.794,6.9317,2.655), tol=1e-4) + + umf2 <- unmarkedFrameMPois(y=y, siteCovs=data.frame(x1=rnorm(5)), type="removal") + m4 <- multinomPois(~1~x1, umf2) + r <- ranef(m4, K=30) + expect_equal(dim(r@post), c(5,31,1)) expect_warning(s <- simulate(m2_C, 2, na.rm=FALSE)) expect_equal(length(s), 2) diff --git a/vignettes/occuMulti.Rmd b/vignettes/occuMulti.Rmd index 25eae85..59f602e 100644 --- a/vignettes/occuMulti.Rmd +++ b/vignettes/occuMulti.Rmd @@ -64,7 +64,7 @@ The data are a simplified version of the data used in @Rota2016, with the data c The dataset is included with `unmarked` and is called `MesoCarnivores`. First, we need to load in the dataset, which is a list with several components. -```{r, eval=FALSE} +```{r} library(unmarked) data(MesoCarnivores) names(MesoCarnivores) |