Support in R ll function, update fitted method and add tests

3 files changed, 82 insertions, 12 deletions

diff --git a/R/gdistsamp.R b/R/gdistsamp.R
index 54ff6fb..2c09bbf 100644
--- a/R/gdistsamp.R
+++ b/R/gdistsamp.R
@@ -21,10 +21,6 @@ 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)
@@ -171,8 +167,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) {
@@ -239,8 +236,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) {
@@ -309,8 +307,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) {
@@ -372,8 +371,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) {
diff --git a/R/unmarkedFit.R b/R/unmarkedFit.R
index d259732..a856288 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
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)
+
+})