pcount: Fit the N-mixture model of Royle (2004)

View source: R/pcount.R

pcountR Documentation

Fit the N-mixture model of Royle (2004)

Description

Fit the N-mixture model of Royle (2004)

Usage

pcount(formula, data, K, mixture=c("P", "NB", "ZIP"),
    starts, method="BFGS", se=TRUE, engine=c("C", "R", "TMB"), threads=1, ...)

Arguments

formula

Double right-hand side formula describing covariates of detection and abundance, in that order

data

an unmarkedFramePCount object supplying data to the model.

K

Integer upper index of integration for N-mixture. This should be set high enough so that it does not affect the parameter estimates. Note that computation time will increase with K.

mixture

character specifying mixture: "P", "NB", or "ZIP".

starts

vector of starting values

method

Optimization method used by optim.

se

logical specifying whether or not to compute standard errors.

engine

Either "C", "R", or "TMB" to use fast C++ code, native R code, or TMB (required for random effects) during the optimization.

threads

Set the number of threads to use for optimization in C++, if OpenMP is available on your system. Increasing the number of threads may speed up optimization in some cases by running the likelihood calculation in parallel. If threads=1 (the default), OpenMP is disabled.

...

Additional arguments to optim, such as lower and upper bounds

Details

This function fits N-mixture model of Royle (2004) to spatially replicated count data.

See unmarkedFramePCount for a description of how to format data for pcount.

This function fits the latent N-mixture model for point count data (Royle 2004, Kery et al 2005).

The latent abundance distribution, f(N | \mathbf{\theta}) can be set as a Poisson, negative binomial, or zero-inflated Poisson random variable, depending on the setting of the mixture argument, mixture = "P", mixture = "NB", mixture = "ZIP" respectively. For the first two distributions, the mean of N_i is \lambda_i. If N_i \sim NB, then an additional parameter, \alpha, describes dispersion (lower \alpha implies higher variance). For the ZIP distribution, the mean is \lambda_i(1-\psi), where psi is the zero-inflation parameter.

The detection process is modeled as binomial: y_{ij} \sim Binomial(N_i, p_{ij}).

Covariates of \lambda_i use the log link and covariates of p_{ij} use the logit link.

Value

unmarkedFit object describing the model fit.

Author(s)

Ian Fiske and Richard Chandler

References

Royle, J. A. (2004) N-Mixture Models for Estimating Population Size from Spatially Replicated Counts. Biometrics 60, pp. 108–105.

Kery, M., Royle, J. A., and Schmid, H. (2005) Modeling Avaian Abundance from Replicated Counts Using Binomial Mixture Models. Ecological Applications 15(4), pp. 1450–1461.

Johnson, N.L, A.W. Kemp, and S. Kotz. (2005) Univariate Discrete Distributions, 3rd ed. Wiley.

See Also

unmarkedFramePCount, pcountOpen, ranef, parboot

Examples


## Not run: 

# Simulate data
set.seed(35)
nSites <- 100
nVisits <- 3
x <- rnorm(nSites)               # a covariate
beta0 <- 0
beta1 <- 1
lambda <- exp(beta0 + beta1*x)   # expected counts at each site
N <- rpois(nSites, lambda)       # latent abundance
y <- matrix(NA, nSites, nVisits)
p <- c(0.3, 0.6, 0.8)            # detection prob for each visit
for(j in 1:nVisits) {
  y[,j] <- rbinom(nSites, N, p[j])
  }

# Organize data
visitMat <- matrix(as.character(1:nVisits), nSites, nVisits, byrow=TRUE)

umf <- unmarkedFramePCount(y=y, siteCovs=data.frame(x=x),
    obsCovs=list(visit=visitMat))
summary(umf)

# Fit a model
fm1 <- pcount(~visit-1 ~ x, umf, K=50)
fm1

plogis(coef(fm1, type="det")) # Should be close to p


# Empirical Bayes estimation of random effects
(fm1re <- ranef(fm1))
plot(fm1re, subset=site %in% 1:25, xlim=c(-1,40))
sum(bup(fm1re))         # Estimated population size
sum(N)                  # Actual population size


# Real data
data(mallard)
mallardUMF <- unmarkedFramePCount(mallard.y, siteCovs = mallard.site,
obsCovs = mallard.obs)
(fm.mallard <- pcount(~ ivel+ date + I(date^2) ~ length + elev + forest, mallardUMF, K=30))
(fm.mallard.nb <- pcount(~ date + I(date^2) ~ length + elev, mixture = "NB", mallardUMF, K=30))


## End(Not run)


unmarked documentation built on July 9, 2023, 5:44 p.m.