diff options
author | Ken Kellner <kenkellner@users.noreply.github.com> | 2023-07-27 15:23:31 -0400 |
---|---|---|
committer | GitHub <noreply@github.com> | 2023-07-27 15:23:31 -0400 |
commit | 516906afa4166978133c61541fa32f80ebda4416 (patch) | |
tree | 2aa38ffb6f9e78cdd4d0164f19589bd9336c3983 | |
parent | cd4783a0ee6aa9e17c7fe6e9901e780b4bf9407c (diff) | |
parent | ae793d82bf05b88826bdad4ccb1f685fc7c68def (diff) |
Merge pull request #260 from rbchan/gdistsamp_ZIP
Add ZIP support to gdistsamp, gmultmix, and gpcount
-rw-r--r-- | DESCRIPTION | 4 | ||||
-rw-r--r-- | R/gdistsamp.R | 39 | ||||
-rw-r--r-- | R/gmultmix.R | 25 | ||||
-rw-r--r-- | R/gpcount.R | 21 | ||||
-rw-r--r-- | R/ranef.R | 15 | ||||
-rw-r--r-- | R/unmarkedFit.R | 30 | ||||
-rw-r--r-- | man/gdistsamp.Rd | 6 | ||||
-rw-r--r-- | man/gmultmix.Rd | 26 | ||||
-rw-r--r-- | man/gpcount.Rd | 11 | ||||
-rw-r--r-- | tests/testthat/test_gdistsamp.R | 65 | ||||
-rw-r--r-- | tests/testthat/test_gmultmix.R | 40 | ||||
-rw-r--r-- | tests/testthat/test_gpcount.R | 47 |
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, @@ -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) +}) |