diff options
author | Ken Kellner <ken@kenkellner.com> | 2023-01-31 13:01:48 -0500 |
---|---|---|
committer | Ken Kellner <ken@kenkellner.com> | 2023-01-31 13:05:03 -0500 |
commit | 52e222e804473ee9c42ac7a826b47ac667725c3d (patch) | |
tree | bcf3a54b5c37099c273ed3dac4b399d1f56a0700 | |
parent | 8bf6bf43ca45808a5839e091534c6b4a4b494741 (diff) |
Add ZIP to gdistsamp, not tested
-rw-r--r-- | R/gdistsamp.R | 19 | ||||
-rw-r--r-- | R/ranef.R | 7 | ||||
-rw-r--r-- | R/unmarkedFit.R | 7 | ||||
-rw-r--r-- | man/gdistsamp.Rd | 6 |
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, @@ -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. |