Add ZIP to gdistsamp, not tested

4 files changed, 31 insertions, 8 deletions

diff --git a/R/gdistsamp.R b/R/gdistsamp.R
index b728d01..54ff6fb 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"))
@@ -21,6 +21,10 @@ survey <- data@survey
 unitsIn <- data@unitsIn
 mixture <- match.arg(mixture)
 
+if(mixture == "ZIP" & engine=="R"){
+  stop("R engine does not support ZIP", call.=FALSE)
+}
+
 formlist <- list(lambdaformula = lambdaformula, phiformula = phiformula,
     pformula = pformula)
 check_no_support(formlist)
@@ -112,6 +116,10 @@ if(identical(mixture, "NB")) {
     nOP <- 1
     nbPar <- "alpha"
     }
+if(identical(mixture, "ZIP")) {
+    nOP <- 1
+    nbPar <- "psi"
+    }
 else {
     nOP <- 0
     nbPar <- character(0)
@@ -398,7 +406,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 +461,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/ranef.R b/R/ranef.R
index 9dca184..350ed6d 100644
--- a/R/ranef.R
+++ b/R/ranef.R
@@ -434,12 +434,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) {
diff --git a/R/unmarkedFit.R b/R/unmarkedFit.R
index d0fac62..d259732 100644
--- a/R/unmarkedFit.R
+++ b/R/unmarkedFit.R
@@ -3028,7 +3028,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.