Merge pull request #171 from kenkellner/simulate_multinom_NA

Fix simulate() to work with multinomial distributions when there are missing values

R/unmarkedFit.R

@@ -4233,7 +4233,8 @@ multinomOpenSim <- function(object, nsim, seed, na.rm){
       yst <- 1
         for(t in 1:T) {
           yend <- yst + J - 1
-          y.it <- as.integer(rmultinom(1, N[i,t], prob=cp[i,,t]))
+          #rmultinom2 in utils.R
+          y.it <- as.integer(rmultinom2(1, N[i,t], prob=cp[i,,t]))
           y.sim[i,yst:yend] <- y.it[1:J]
           yst <- yst + J
         }
@@ -4357,7 +4358,7 @@ setMethod("simulate", "unmarkedFitOccuFP",
               P[,1] <- Z*rbinom(M * J, 1, prob = (1-p)) + (1-Z)*rbinom(M * J, 1, prob = (1-fp))
               P[,2] <- (1-P[,1])*(1-Z) + (1-P[,1])*rbinom(M * J, 1, prob = (1-b))*Z
               P[,3] <- 1 - P[,1]-P[,2]
-              yvec <- sapply(1:(M*J),function(x) which(as.logical(rmultinom(1,1,P[x,])))-1)
+              yvec <- sapply(1:(M*J),function(x) which(as.logical(rmultinom2(1,1,P[x,])))-1)
               simList[[i]] <- matrix(yvec, M, J, byrow = TRUE)
             }
             return(simList)
@@ -4744,7 +4745,7 @@ setMethod("simulate", "unmarkedFitGMM",
                 pi.it <- cp.arr[i,t,]
                 na.it <- is.na(pi.it)
                 pi.it[na.it] <- 0
-                y.sim[i,,t] <- drop(rmultinom(1, N[i,t], pi.it))[1:J]
+                y.sim[i,,t] <- drop(rmultinom2(1, N[i,t], pi.it))[1:J]
                 y.sim[i,na.it[1:J],t] <- NA
             }
         }
@@ -4885,7 +4886,7 @@ setMethod("simulate", "unmarkedFitGDS",
     for(s in 1:nsim) {
         for(i in 1:M) {
             switch(mixture,
-                P = Ns <- rpois(1, lambda[1]),
+                P = Ns <- rpois(1, lambda[i]),
                 NB = {
                     alpha <- exp(coef(object, type="alpha"))
                     Ns <- rnbinom(1, mu=lambda[i], size=alpha)
@@ -4894,7 +4895,7 @@ setMethod("simulate", "unmarkedFitGDS",
                 N <- rbinom(1, Ns, phi[i,t])
                 cp.it <- cpa[i,,t]
                 cp.it[J+1] <- 1-sum(cp.it)
-                y.it <- as.integer(rmultinom(1, N, prob=cp.it))
+                y.it <- as.integer(rmultinom2(1, N, prob=cp.it))
                 ysim[i,,t] <- y.it[1:J]
                 }
             }

R/utils.R

@@ -838,3 +838,14 @@ getDistCP <- function(keyfun, param1, param2, survey, db, w, a, u){
     })
     cp * u
 }
+
+
+#Modified rmultinom for handling NAs
+rmultinom2 <- function(n, size, prob){
+  if(is.na(size)){
+    return(matrix(NA, length(prob), length(n)))
+  }
+  stats::rmultinom(n=n, size=size, prob=prob)
+
+}
+

inst/unitTests/runit.simulate.R

@@ -0,0 +1,52 @@
+test.simulate.GDS <- function(){
+
+    set.seed(343)
+    R <- 30 
+    T <- 3  
+    strip.width <- 50
+    transect.length <- 200 #Area != 1
+    breaks <- seq(0, 50, by=10)
+
+    covs <- as.data.frame(matrix(rnorm(R*T),ncol=T))
+    names(covs) <- paste0('par',1:3)
+
+    beta <- c(0.4,0.3,0.6)
+    lambda <- exp(1.3 + beta[1]*covs$par1) 
+    phi <- plogis(as.matrix(0.4 + beta[2]*covs)) 
+    sigma <- exp(as.matrix(3 + beta[3]*covs)) 
+    J <- length(breaks)-1
+    y <- array(0, c(R, J, T))
+    for(i in 1:R) {
+        M <- rpois(1, lambda[i]) # Individuals within the 1-ha strip
+        for(t in 1:T) {
+            # Distances from point
+            d <- runif(M, 0, strip.width)
+            # Detection process
+            if(length(d)) {
+                cp <- phi[i,t]*exp(-d^2 / (2 * sigma[i,t]^2)) # half-normal w/ g(0)<1
+                d <- d[rbinom(length(d), 1, cp) == 1]
+                y[i,,t] <- table(cut(d, breaks, include.lowest=TRUE))
+            }
+        }
+    }
+    y <- matrix(y, nrow=R) # convert array to matrix
+
+    covs$par1[2] <- NA
+    umf <- unmarkedFrameGDS(y = y, siteCovs=covs, yearlySiteCovs=covs,
+                            survey="line", unitsIn="m",
+                            dist.breaks=breaks,
+                            tlength=rep(transect.length, R), numPrimary=T)
+
+    fm <- gdistsamp(~par1, ~1, ~1, umf, se=FALSE, engine="C")
+
+    #This used to error due to rmultinom not accepting size=NA
+    s <- simulate(fm, nsim=2, na.rm=FALSE)
+    checkEqualsNumeric(length(s), 2)
+    checkEqualsNumeric(dim(s[[1]]), c(30,15))
+    checkTrue(!any(is.na(s[[1]][1,])))
+    checkTrue(all(is.na(s[[1]][2,])))
+
+    pb <- parboot(fm, nsim=3)
+    checkTrue(inherits(pb, "parboot"))
+
+}