Merge pull request #260 from rbchan/gdistsamp_ZIP

Add ZIP support to gdistsamp, gmultmix, and gpcount

12 files changed, 268 insertions, 61 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index ba0d28e..b658bba 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: unmarked
-Version: 1.3.2
-Date: 2023-07-06
+Version: 1.3.2.9001
+Date: 2023-07-27
 Type: Package
 Title: Models for Data from Unmarked Animals
 Authors@R: c(
diff --git a/R/gdistsamp.R b/R/gdistsamp.R
index b728d01..c43076a 100644
--- a/R/gdistsamp.R
+++ b/R/gdistsamp.R
@@ -2,7 +2,7 @@
 gdistsamp <- function(lambdaformula, phiformula, pformula, data,
     keyfun=c("halfnorm", "exp", "hazard", "uniform"),
     output=c("abund", "density"), unitsOut=c("ha", "kmsq"),
-    mixture=c('P', 'NB'), K, starts, method = "BFGS", se = TRUE, engine=c("C","R"),
+    mixture=c("P", "NB", 'ZIP'), K, starts, method = "BFGS", se = TRUE, engine=c("C","R"),
     rel.tol=1e-4, threads=1, ...)
 {
 if(!is(data, "unmarkedFrameGDS"))
@@ -111,11 +111,13 @@ else {
 if(identical(mixture, "NB")) {
     nOP <- 1
     nbPar <- "alpha"
-    }
-else {
+} else if(identical(mixture, "ZIP")) {
+    nOP <- 1
+    nbPar <- "psi"
+} else {
     nOP <- 0
     nbPar <- character(0)
-    }
+}
 
 nLP <- ncol(Xlam)
 nP  <- nLP + nPP + nDP + nSP + nOP
@@ -163,8 +165,9 @@ halfnorm = {
 
         switch(mixture,
             P = f <- sapply(k, function(x) dpois(x, lambda)),
-            NB = f <- sapply(k, function(x) dnbinom(x, mu=lambda,
-                size=exp(pars[nP]))))
+            NB = f <- sapply(k, function(x) dnbinom(x, mu=lambda, size=exp(pars[nP]))),
+            ZIP = f <- sapply(k, function(x) dzip(rep(x, length(lambda)), lambda=lambda, psi=plogis(pars[nP])))
+        )
         for(i in 1:M) {
             mn <- matrix(0, lk, T)
             for(t in 1:T) {
@@ -231,8 +234,9 @@ exp = {
 
         switch(mixture,
             P = f <- sapply(k, function(x) dpois(x, lambda)),
-            NB = f <- sapply(k, function(x) dnbinom(x, mu=lambda,
-                size=exp(pars[nP]))))
+            NB = f <- sapply(k, function(x) dnbinom(x, mu=lambda, size=exp(pars[nP]))),
+            ZIP = f <- sapply(k, function(x) dzip(rep(x, length(lambda)), lambda=lambda, psi=plogis(pars[nP])))
+        )
         for(i in 1:M) {
             mn <- matrix(0, lk, T)
             for(t in 1:T) {
@@ -301,8 +305,9 @@ hazard = {
 
         switch(mixture,
             P = f <- sapply(k, function(x) dpois(x, lambda)),
-            NB = f <- sapply(k, function(x) dnbinom(x, mu=lambda,
-                size=exp(pars[nP]))))
+            NB = f <- sapply(k, function(x) dnbinom(x, mu=lambda, size=exp(pars[nP]))),
+            ZIP = f <- sapply(k, function(x) dzip(rep(x, length(lambda)), lambda=lambda, psi=plogis(pars[nP])))
+        )
         for(i in 1:M) {
             mn <- matrix(0, lk, T)
             for(t in 1:T) {
@@ -364,8 +369,9 @@ uniform = {
         p <- 1
         switch(mixture,
             P = f <- sapply(k, function(x) dpois(x, lambda)),
-            NB = f <- sapply(k, function(x) dnbinom(x, mu=lambda,
-                size=exp(pars[nP]))))
+            NB = f <- sapply(k, function(x) dnbinom(x, mu=lambda, size=exp(pars[nP]))),
+            ZIP = f <- sapply(k, function(x) dzip(rep(x, length(lambda)), lambda=lambda, psi=plogis(pars[nP])))
+        )
         for(i in 1:M) {
             mn <- matrix(0, lk, T)
             for(t in 1:T) {
@@ -398,7 +404,7 @@ if(engine =="C"){
   if(output!='density'){
     A <- rep(1, M)
   }
-  mixture_code <- switch(mixture, P={1}, NB={2})
+  mixture_code <- switch(mixture, P={1}, NB={2}, ZIP={3})
   n_param <- c(nLP, nPP, nDP, nSP, nOP)
   Kmin <- apply(yt, 1, max, na.rm=TRUE)
 
@@ -453,6 +459,13 @@ if(identical(mixture, "NB"))
         covMat = as.matrix(covMat[nP, nP]), invlink = "exp",
         invlinkGrad = "exp")
 
+if(identical(mixture,"ZIP")) {
+    estimateList@estimates$psi <- unmarkedEstimate(name="Zero-inflation",
+        short.name = "psi", estimates = ests[nP],
+        covMat=as.matrix(covMat[nP, nP]), invlink = "logistic",
+        invlinkGrad = "logistic.grad")
+}
+
 umfit <- new("unmarkedFitGDS", fitType = "gdistsamp",
     call = match.call(), formula = form, formlist = formlist,
     data = data, estimates = estimateList, sitesRemoved = D$removed.sites,
diff --git a/R/gmultmix.R b/R/gmultmix.R
index 6dabe9b..a6da9d5 100644
--- a/R/gmultmix.R
+++ b/R/gmultmix.R
@@ -1,7 +1,7 @@
 
 # data will need to be an unmarkedMultFrame
 gmultmix <- function(lambdaformula, phiformula, pformula, data,
-    mixture=c('P', 'NB'), K, starts, method = "BFGS", se = TRUE,
+    mixture=c("P", "NB", "ZIP"), K, starts, method = "BFGS", se = TRUE,
     engine=c("C","R"), threads=1, ...)
 {
 if(!is(data, "unmarkedFrameGMM"))
@@ -59,7 +59,7 @@ if(T==1) {
     phiPars <- colnames(Xphi)
     }
 nDP <- ncol(Xdet)
-nP <- nLP + nPP + nDP + (mixture=='NB')
+nP <- nLP + nPP + nDP + (mixture%in%c('NB','ZIP'))
 if(!missing(starts) && length(starts) != nP)
     stop(paste("The number of starting values should be", nP))
 
@@ -102,9 +102,10 @@ nll_R <- function(pars) {
     cp[,,R+1] <- 1 - apply(cp[,,1:R,drop=FALSE], 1:2, sum, na.rm=TRUE)
 
     switch(mixture,
-        P = f <- sapply(k, function(x) dpois(x, lambda)),
-        NB = f <- sapply(k, function(x) dnbinom(x, mu=lambda,
-            size=exp(pars[nP]))))
+      P = f <- sapply(k, function(x) dpois(x, lambda)),
+      NB = f <- sapply(k, function(x) dnbinom(x, mu=lambda, size=exp(pars[nP]))),
+      ZIP = f <- sapply(k, function(x) dzip(rep(x, length(lambda)), lambda=lambda, psi=plogis(pars[nP])))
+    )
     g <- matrix(as.numeric(NA), M, lk)
     for(i in 1:M) {
         A <- matrix(0, lk, T)
@@ -134,8 +135,8 @@ if(engine=="R"){
   kmytC <- kmyt
   kmytC[which(is.na(kmyt))] <- 0
 
-  mixture_code <- switch(mixture, P={1}, NB={2})
-  n_param <- c(nLP, nPP, nDP, mixture=="NB")
+  mixture_code <- switch(mixture, P={1}, NB={2}, ZIP={3})
+  n_param <- c(nLP, nPP, nDP, mixture%in%c("NB","ZIP"))
   Kmin <- apply(yt, 1, max, na.rm=TRUE)
 
   nll <- function(params) {
@@ -157,8 +158,7 @@ covMat <- invertHessian(fm, nP, se)
 ests <- fm$par
 fmAIC <- 2 * fm$value + 2 * nP
 
-if(identical(mixture,"NB")) nbParm <- "alpha"
-	else nbParm <- character(0)
+nbParm <- switch(mixture, P={character(0)}, NB={"alpha"}, ZIP={"psi"})
 
 names(ests) <- c(lamPars, phiPars, detPars, nbParm)
 
@@ -192,6 +192,13 @@ if(identical(mixture,"NB"))
         covMat = as.matrix(covMat[nP, nP]), invlink = "exp",
         invlinkGrad = "exp")
 
+if(identical(mixture,"ZIP")) {
+    estimateList@estimates$psi <- unmarkedEstimate(name="Zero-inflation",
+        short.name = "psi", estimates = ests[nP],
+        covMat=as.matrix(covMat[nP, nP]), invlink = "logistic",
+        invlinkGrad = "logistic.grad")
+}
+
 umfit <- new("unmarkedFitGMM", fitType = "gmn",
     call = match.call(), formula = form, formlist = formlist,
     data = data, estimates = estimateList, sitesRemoved = D$removed.sites,
diff --git a/R/gpcount.R b/R/gpcount.R
index 2b7afd6..af4dade 100644
--- a/R/gpcount.R
+++ b/R/gpcount.R
@@ -1,15 +1,13 @@
 
 # data will need to be an unmarkedMultFrame
 gpcount <- function(lambdaformula, phiformula, pformula, data,
-    mixture=c('P', 'NB'), K, starts, method = "BFGS", se = TRUE,
+    mixture=c('P', 'NB', 'ZIP'), K, starts, method = "BFGS", se = TRUE,
     engine=c('C', 'R'), threads=1, ...)
 {
 if(!is(data, "unmarkedFrameGPC"))
     stop("Data is not of class unmarkedFrameGPC.")
 mixture <- match.arg(mixture)
 engine <- match.arg(engine)
-if(identical(mixture, "ZIP") & identical(engine, "R"))
-    stop("ZIP mixture not available when 'engine=R'")
 
 formlist <- list(lambdaformula = lambdaformula, phiformula = phiformula,
     pformula = pformula)
@@ -49,7 +47,7 @@ nLP <- ncol(Xlam)
 nPP <- ncol(Xphi)
 phiPars <- colnames(Xphi)
 nDP <- ncol(Xdet)
-nP <- nLP + nPP + nDP + (mixture=='NB')
+nP <- nLP + nPP + nDP + (mixture%in%c('NB','ZIP'))
 if(!missing(starts) && length(starts) != nP)
     stop("There should be", nP, "starting values, not", length(starts))
 
@@ -66,8 +64,9 @@ nll <- function(pars) {
     for(i in 1:I) {
         f <- switch(mixture,
             P = dpois(M, lam[i], log=TRUE),
-            NB = dnbinom(M, mu=lam[i], size=exp(pars[nP]), log=TRUE))
-#            ZIP = dzip())
+            NB = dnbinom(M, mu=lam[i], size=exp(pars[nP]), log=TRUE),
+            ZIP = log(dzip(M, lambda=lam[i], psi=plogis(pars[nP])))
+        )
         ghi <- rep(0, lM)
         for(t in 1:T) {
             gh <- matrix(-Inf, lM, lM)
@@ -118,8 +117,7 @@ covMat <- invertHessian(fm, nP, se)
 ests <- fm$par
 fmAIC <- 2 * fm$value + 2 * nP
 
-if(identical(mixture,"NB")) nbParm <- "alpha"
-	else nbParm <- character(0)
+nbParm <- switch(mixture, P={character(0)}, NB={"alpha"}, ZIP={"psi"})
 
 names(ests) <- c(lamPars, phiPars, detPars, nbParm)
 
@@ -150,6 +148,13 @@ if(identical(mixture,"NB"))
         covMat = as.matrix(covMat[nP, nP]), invlink = "exp",
         invlinkGrad = "exp")
 
+if(identical(mixture,"ZIP")) {
+    estimateList@estimates$psi <- unmarkedEstimate(name="Zero-inflation",
+        short.name = "psi", estimates = ests[nP],
+        covMat=as.matrix(covMat[nP, nP]), invlink = "logistic",
+        invlinkGrad = "logistic.grad")
+}
+
 umfit <- new("unmarkedFitGPC", fitType = "gpcount",
     call = match.call(), formula = form, formlist = formlist,
     data = data, estimates = estimateList, sitesRemoved = D$removed.sites,
diff --git a/R/ranef.R b/R/ranef.R
index 669dc1a..f99efaf 100644
--- a/R/ranef.R
+++ b/R/ranef.R
@@ -445,12 +445,15 @@ setMethod("ranef", "unmarkedFitGMMorGDS",
     post <- array(0, c(nSites, K+1, 1))
     colnames(post) <- M
     mix <- object@mixture
-    if(identical(mix, "NB"))
+    if(identical(mix, "NB")){
         alpha <- exp(coef(object, type="alpha"))
+    } else if(identical(mix, "ZIP")){
+        psi <- plogis(coef(object, type="psi"))
+    }
     for(i in 1:nSites) {
         switch(mix,
                P  = f <- dpois(M, lambda[i]),
-               # FIXME: Add ZIP
+               ZIP = f <- dzip(M, lambda[i], psi),
                NB = f <- dnbinom(M, mu=lambda[i], size=alpha))
         g <- rep(1, K+1) # outside t loop
         for(t in 1:T) {
@@ -528,11 +531,15 @@ setMethod("ranef", "unmarkedFitGPC",
     mix <- object@mixture
     if(identical(mix, "NB"))
         alpha <- exp(coef(object, type="alpha"))
+    if(identical(mix, "ZIP"))
+        psi <- plogis(coef(object, type="psi"))
+
     for(i in 1:R) {
         switch(mix,
                P  = f <- dpois(M, lambda[i]),
-               # FIXME: Add ZIP
-               NB = f <- dnbinom(M, mu=lambda[i], size=alpha))
+               NB = f <- dnbinom(M, mu=lambda[i], size=alpha),
+               ZIP = f <- dzip(M, lambda=lambda[i], psi=psi)
+        )
         ghi <- rep(0, lM)
         for(t in 1:T) {
             gh <- matrix(-Inf, lM, lM)
diff --git a/R/unmarkedFit.R b/R/unmarkedFit.R
index d0fac62..95c832c 100644
--- a/R/unmarkedFit.R
+++ b/R/unmarkedFit.R
@@ -826,6 +826,11 @@ setMethod("fitted", "unmarkedFitGMM",
     T <- data@numPrimary
     J <- ncol(y) / T
     lambda <- drop(exp(Xlam %*% coef(object, 'lambda') + Xlam.offset))
+    if(identical(object@mixture, "ZIP")) {
+        psi <- plogis(coef(object, type="psi"))
+        lambda <- (1-psi)*lambda
+    }
+
     if(T==1)
         phi <- 1
     else
@@ -2869,7 +2874,12 @@ setMethod("simulate", "unmarkedFitGMM",
         switch(mixture,
             P = M <- rpois(n=n, lambda=lam),
             NB = M <- rnbinom(n=n, mu=lam,
-                size=exp(coef(object, type="alpha"))))
+                size=exp(coef(object, type="alpha"))),
+            ZIP = {
+              psi <- plogis(coef(object['psi']))
+              M <- rzip(n, lambda=lam, psi=psi)
+            }
+        )
 
         N <- rbinom(n*T, size=M, prob=phi.mat)
         # bug fix 3/16/2010
@@ -2935,10 +2945,15 @@ setMethod("simulate", "unmarkedFitGPC",
     simList <- list()
     for(s in 1:nsim) {
         switch(mixture,
-               P = M <- rpois(n=R, lambda=lam),
+            P = M <- rpois(n=R, lambda=lam),
                #               FIXME: Add ZIP
-               NB = M <- rnbinom(n=R, mu=lam,
-               size=exp(coef(object, type="alpha"))))
+            NB = M <- rnbinom(n=R, mu=lam,
+               size=exp(coef(object, type="alpha"))),   
+            ZIP = {
+              psi <- plogis(coef(object['psi']))
+              M <- rzip(R, lambda=lam, psi=psi)
+            }
+        )
 
         N <- rbinom(R*T, size=M, prob=phi.mat)
         N <- matrix(N, nrow=R, ncol=T, byrow=FALSE)
@@ -3028,7 +3043,12 @@ setMethod("simulate", "unmarkedFitGDS",
                 NB = {
                     alpha <- exp(coef(object, type="alpha"))
                     Ns <- rnbinom(1, mu=lambda[i], size=alpha)
-                    })
+                    },
+                ZIP = {
+                    psi <- plogis(coef(object['psi']))
+                    Ns <- rzip(1, lambda[i], psi)
+                    }
+            )
             for(t in 1:T) {
                 N <- rbinom(1, Ns, phi[i,t])
                 cp.it <- cpa[i,,t]
diff --git a/man/gdistsamp.Rd b/man/gdistsamp.Rd
index 323f5cc..1e15cac 100644
--- a/man/gdistsamp.Rd
+++ b/man/gdistsamp.Rd
@@ -11,7 +11,7 @@ to be modeled using the negative binomial distribution.
 \usage{
 gdistsamp(lambdaformula, phiformula, pformula, data, keyfun =
 c("halfnorm", "exp", "hazard", "uniform"), output = c("abund",
-"density"), unitsOut = c("ha", "kmsq"), mixture = c("P", "NB"), K,
+"density"), unitsOut = c("ha", "kmsq"), mixture = c("P", "NB", "ZIP"), K,
 starts, method = "BFGS", se = TRUE, engine=c("C","R"), rel.tol=1e-4, threads=1, ...)
 }
 \arguments{
@@ -39,8 +39,8 @@ starts, method = "BFGS", se = TRUE, engine=c("C","R"), rel.tol=1e-4, threads=1,
           kilometers, respectively.
 }
   \item{mixture}{
-      Either "P" or "NB" for the Poisson and negative binomial models of
-          abundance.
+      Either "P", "NB", or "ZIP" for the Poisson, negative binomial, or
+      zero-inflated Poisson models of abundance.
 }
   \item{K}{
       An integer value specifying the upper bound used in the integration.
diff --git a/man/gmultmix.Rd b/man/gmultmix.Rd
index b54701c..40f6007 100644
--- a/man/gmultmix.Rd
+++ b/man/gmultmix.Rd
@@ -9,7 +9,7 @@ The three model parameters are abundance, availability, and detection
 probability.
 }
 \usage{
-gmultmix(lambdaformula, phiformula, pformula, data, mixture = c("P", "NB"), K, 
+gmultmix(lambdaformula, phiformula, pformula, data, mixture = c("P", "NB", "ZIP"), K, 
          starts, method = "BFGS", se = TRUE, engine=c("C","R"), threads=1, ...)
 }
 \arguments{
@@ -18,8 +18,8 @@ gmultmix(lambdaformula, phiformula, pformula, data, mixture = c("P", "NB"), K,
     \item{phiformula}{RHS formula describing availability covariates}
     \item{pformula}{RHS formula describing detection covariates}
     \item{data}{An object of class unmarkedFrameGMM}
-    \item{mixture}{Either "P" or "NB" for Poisson and Negative Binomial mixing
-        distributions.}
+    \item{mixture}{Either "P", "NB", or "ZIP" for the Poisson, negative binomial, or
+      zero-inflated Poisson models of abundance}
     \item{K}{The upper bound of integration}
     \item{starts}{Starting values}
     \item{method}{Optimization method used by \code{\link{optim}}}
@@ -35,16 +35,16 @@ gmultmix(lambdaformula, phiformula, pformula, data, mixture = c("P", "NB"), K,
     }
 \details{
 The latent transect-level super-population abundance distribution
-\eqn{f(M | \mathbf{\theta})}{f(M | theta)} can be set as either a
-Poisson or
-a negative binomial random variable, depending on the setting of the
-\code{mixture} argument. \code{mixture = "P"} or \code{mixture = "NB"}
-select
-the Poisson or negative binomial distribution respectively.  The mean of
-\eqn{M_i} is \eqn{\lambda_i}{lambda_i}.  If \eqn{M_i \sim NB}{M_i ~ NB},
-then an
-additional parameter, \eqn{\alpha}{alpha}, describes dispersion (lower
-\eqn{\alpha}{alpha} implies higher variance).
+\eqn{f(M | \mathbf{\theta})}{f(M | theta)} can be set as a
+Poisson, negative binomial, or zero-inflated Poisson random variable, 
+depending on the setting of the \code{mixture} argument. 
+\code{mixture = "P"}, \code{mixture = "NB"}, and \code{mixture = "ZIP"}
+select the Poisson, negative binomial, and zero-inflated Poisson distributions 
+respectively.  The mean of \eqn{M_i} is \eqn{\lambda_i}{lambda_i}.  
+If \eqn{M_i \sim NB}{M_i ~ NB}, then an additional parameter, \eqn{\alpha}{alpha}, 
+describes dispersion (lower \eqn{\alpha}{alpha} implies higher variance). If 
+\eqn{M_i \sim ZIP}{M_i ~ ZIP}, then an additional zero-inflation parameter
+\eqn{\psi}{psi} is estimated.
 
 The number of individuals available for detection at time j
 is a modeled as binomial:
diff --git a/man/gpcount.Rd b/man/gpcount.Rd
index f847baa..b7e3b1a 100644
--- a/man/gpcount.Rd
+++ b/man/gpcount.Rd
@@ -10,7 +10,7 @@ population size, availability, and detection probability.
 }
 \usage{
 gpcount(lambdaformula, phiformula, pformula, data,
-mixture = c("P", "NB"), K, starts, method = "BFGS", se = TRUE,
+mixture = c("P", "NB", "ZIP"), K, starts, method = "BFGS", se = TRUE,
 engine = c("C", "R"), threads=1, ...)
 }
 \arguments{
@@ -27,7 +27,8 @@ engine = c("C", "R"), threads=1, ...)
     An object of class unmarkedFrameGPC
 }
   \item{mixture}{
-    Either "P" or "NB" for Poisson and negative binomial distributions
+    Either "P", "NB", or "ZIP" for Poisson, negative binomial, or 
+    zero-inflated Poisson distributions
 }
   \item{K}{
     The maximum possible value of M, the super-population size.
@@ -58,11 +59,13 @@ engine = c("C", "R"), threads=1, ...)
 \details{
   The latent transect-level super-population abundance distribution
   \eqn{f(M | \mathbf{\theta})}{f(M | theta)} can be set as either a
-  Poisson or a negative binomial random variable, depending on the
+  Poisson, negative binomial, or zero-inflated Poisson random variable, depending on the
   setting of the \code{mixture} argument. The expected value of
   \eqn{M_i} is \eqn{\lambda_i}{lambda_i}.  If \eqn{M_i \sim NB}{M_i ~ NB},
   then an additional parameter, \eqn{\alpha}{alpha}, describes
-  dispersion (lower \eqn{\alpha}{alpha} implies higher variance).
+  dispersion (lower \eqn{\alpha}{alpha} implies higher variance). If 
+  \eqn{M_i \sim ZIP}{M_i ~ ZIP}, then an additional zero-inflation parameter
+  \eqn{\psi}{psi} is estimated.
 
   The number of individuals available for detection at time j
   is a modeled as binomial:
diff --git a/tests/testthat/test_gdistsamp.R b/tests/testthat/test_gdistsamp.R
index f215d90..c993013 100644
--- a/tests/testthat/test_gdistsamp.R
+++ b/tests/testthat/test_gdistsamp.R
@@ -644,3 +644,68 @@ test_that("gdistsamp simulate method works",{
 
 })
 
+test_that("gdistsamp with ZIP mixture works",{
+    #Line
+    set.seed(343)
+    #R <- 500 # for accuracy checks
+    #T <- 10
+    R <- 50
+    T <- 5
+    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)
+    x <- rnorm(R)
+    lambda <- exp(1.3 + beta[1]*x)
+    psi <- 0.3
+    phi <- plogis(as.matrix(0.4 + beta[2]*covs))
+    sigma <- exp(3)
+    J <- length(breaks)-1
+    y <- array(0, c(R, J, T))
+    M <- numeric(R)
+    for(i in 1:R) {
+        M[i] <- unmarked:::rzip(1, lambda[i], psi=psi) # Individuals within the 1-ha strip
+        for(t in 1:T) {
+            # Distances from point
+            d <- runif(M[i], 0, strip.width)
+            # Detection process
+            if(length(d)) {
+                cp <- phi[i,t]*exp(-d^2 / (2 * sigma^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
+
+    umf <- unmarkedFrameGDS(y = y, survey="line", unitsIn="m",
+                            siteCovs=data.frame(par1=x),
+                            yearlySiteCovs=list(par2=covs),
+                            dist.breaks=breaks,
+                            tlength=rep(transect.length, R), numPrimary=T)
+
+    # R and C give same result
+    fm_R <- gdistsamp(~par1, ~par2, ~1, umf, mixture="ZIP", output="abund", se=FALSE, engine="R",
+                      control=list(maxit=1))
+    fm_C <- gdistsamp(~par1, ~par2, ~1, umf, mixture="ZIP", output="abund", se=FALSE, engine="C",
+                      control=list(maxit=1))
+    expect_equal(coef(fm_R), coef(fm_C))
+    
+    # Fit model
+    fm_C <- gdistsamp(~par1, ~par2, ~1, umf, mixture="ZIP", output="abund")
+    
+    expect_equivalent(coef(fm_C), c(2.4142,0.3379,-1.2809,0.25916,2.91411,-1.12557), tol=1e-4)
+    
+    # Check ZIP-specific methods
+    ft <- fitted(fm_C)
+    r <- ranef(fm_C)
+    b <- bup(r)
+    #plot(M, b)
+    #abline(a=0,b=1)
+    s <- simulate(fm_C)
+
+})
diff --git a/tests/testthat/test_gmultmix.R b/tests/testthat/test_gmultmix.R
index 8389b66..5caafa6 100644
--- a/tests/testthat/test_gmultmix.R
+++ b/tests/testthat/test_gmultmix.R
@@ -287,3 +287,43 @@ test_that("getP works when there is only 1 site", {
   expect_equal(dim(gp), c(1,4))
 
 })
+
+test_that("gmultmix ZIP model works",{
+  # Simulate independent double observer data
+  nSites <- 50
+  lambda <- 10
+  psi <- 0.3
+  p1 <- 0.5
+  p2 <- 0.3
+  cp <- c(p1*(1-p2), p2*(1-p1), p1*p2)
+  set.seed(9023)
+  N <- unmarked:::rzip(nSites, lambda, psi)
+  y <- matrix(NA, nSites, 3)
+  for(i in 1:nSites) {
+    y[i,] <- rmultinom(1, N[i], c(cp, 1-sum(cp)))[1:3]
+  }
+
+  # Make unmarkedFrame
+  observer <- matrix(c('A','B'), nSites, 2, byrow=TRUE)
+  expect_warning(umf <- unmarkedFrameGMM(y=y, obsCovs=list(observer=observer),
+         type="double",numPrimary=1))
+
+  # Compare R and C engines
+  fmR <- gmultmix(~1, ~1, ~observer, umf, mixture="ZIP", engine="R", se=FALSE,
+                  control=list(maxit=1))
+  fmC <- gmultmix(~1, ~1, ~observer, umf, mixture="ZIP", engine="C", se=FALSE,
+                  control=list(maxit=1))
+  expect_equal(coef(fmR), coef(fmC))
+
+  #  Fit model full
+  fm <- gmultmix(~1,~1,~observer, umf, mixture="ZIP")
+  expect_equivalent(coef(fm), c(2.2289,0.1858,-0.9285,-0.9501), tol=1e-4)
+
+  # Check methods
+  ft <- fitted(fm)
+  r <- ranef(fm)
+  b <- bup(r)
+  #plot(N, b)
+  s <- simulate(fm)
+
+})
diff --git a/tests/testthat/test_gpcount.R b/tests/testthat/test_gpcount.R
index 022c2ce..7f23a76 100644
--- a/tests/testthat/test_gpcount.R
+++ b/tests/testthat/test_gpcount.R
@@ -117,3 +117,50 @@ test_that("gpcount R and C++ engines give same results",{
   fmR <- gpcount(~x, ~yr, ~o1, data = umf, K=23, engine="R", control=list(maxit=1))
   expect_equal(coef(fm), coef(fmR))
 })
+
+test_that("gpcount ZIP mixture works", {
+  
+  set.seed(123)
+  M <- 100
+  J <- 5
+  T <- 3
+  lam <- 3
+  psi <- 0.3
+  p <- 0.5
+  phi <- 0.7
+
+  y <- array(NA, c(M, J, T))
+
+  N <- unmarked:::rzip(M, lambda=lam, psi=psi)
+
+  for (i in 1:M){
+    for (t in 1:T){
+      n <- rbinom(1, N[i], phi)
+      for (j in 1:J){
+         y[i,j,t] <- rbinom(1, n, p)
+      }
+    }
+  }
+
+  ywide <- cbind(y[,,1], y[,,2], y[,,3])
+  umf <- unmarkedFrameGPC(y=ywide, numPrimary=T)
+  
+  # check R and C engines match
+  fitC <- gpcount(~1, ~1, ~1, umf, mixture="ZIP", K=10, engine="C",
+                se=FALSE, control=list(maxit=1))
+  fitR <- gpcount(~1, ~1, ~1, umf, mixture="ZIP", K=10, engine="R",
+                se=FALSE, control=list(maxit=1))
+  expect_equal(coef(fitC), coef(fitR))
+  
+  # Properly fit model
+  fit <- gpcount(~1, ~1, ~1, umf, mixture="ZIP", K=10)
+  expect_equivalent(coef(fit), c(1.02437, 0.85104, -0.019588, -1.16139), tol=1e-4)
+
+  # Check methods
+  ft <- fitted(fit)
+  r <- ranef(fit)
+  b <- bup(r)
+  #plot(N, b)
+  #abline(a=0, b=1)
+  s <- simulate(fit)
+})