multmixOpen: Open population multinomial N-mixture model

In unmarked: Models for Data from Unmarked Animals

Open population multinomial N-mixture model

Description

Fit the model of Dail and Madsen (2011) and Hostetler and Chandler (2015) for designs involving repeated counts that yield multinomial outcomes. Possible data collection methods include repeated removal sampling and double observer sampling.

Usage

multmixOpen(lambdaformula, gammaformula, omegaformula, pformula,
    data, mixture=c("P", "NB", "ZIP"), K,
    dynamics=c("constant", "autoreg", "notrend", "trend", "ricker", "gompertz"),
    fix=c("none", "gamma", "omega"), immigration=FALSE, iotaformula = ~1,
    starts, method="BFGS", se=TRUE, ...)

Arguments

lambdaformula	Right-hand sided formula for initial abundance
gammaformula	Right-hand sided formula for recruitment rate (when dynamics is "constant", "autoreg", or "notrend") or population growth rate (when dynamics is "trend", "ricker", or "gompertz")
omegaformula	Right-hand sided formula for apparent survival probability (when dynamics is "constant", "autoreg", or "notrend") or equilibrium abundance (when dynamics is "ricker" or "gompertz")
pformula	A right-hand side formula describing the detection function covariates
data	An object of class unmarkedFrameMMO
mixture	String specifying mixture: "P", "NB", or "ZIP" for the Poisson, negative binomial, or zero-inflated Poisson distributions respectively
K	Integer defining upper bound of discrete integration. This should be higher than the maximum observed count and high enough that it does not affect the parameter estimates. However, the higher the value the slower the computation
dynamics	Character string describing the type of population dynamics. "constant" indicates that there is no relationship between omega and gamma. "autoreg" is an auto-regressive model in which recruitment is modeled as gamma*N[i,t-1]. "notrend" model gamma as lambda*(1-omega) such that there is no temporal trend. "trend" is a model for exponential growth, N[i,t] = N[i,t-1]*gamma, where gamma in this case is finite rate of increase (normally referred to as lambda). "ricker" and "gompertz" are models for density-dependent population growth. "ricker" is the Ricker-logistic model, N[i,t] = N[i,t-1]*exp(gamma*(1-N[i,t-1]/omega)), where gamma is the maximum instantaneous population growth rate (normally referred to as r) and omega is the equilibrium abundance (normally referred to as K). "gompertz" is a modified version of the Gompertz-logistic model, N[i,t] = N[i,t-1]*exp(gamma*(1-log(N[i,t-1]+1)/log(omega+1))), where the interpretations of gamma and omega are similar to in the Ricker model
fix	If "omega", omega is fixed at 1. If "gamma", gamma is fixed at 0
immigration	Logical specifying whether or not to include an immigration term (iota) in population dynamics
iotaformula	Right-hand sided formula for average number of immigrants to a site per time step
starts	Vector of starting values
method	Optimization method used by optim
se	Logical specifying whether or not to compute standard errors
...	Additional arguments to optim, such as lower and upper bounds

Details

These models generalize multinomial N-mixture models (Royle et al. 2004) by relaxing the closure assumption (Dail and Madsen 2011, Hostetler and Chandler 2015, Sollmann et al. 2015).

The models include two or three additional parameters: gamma, either the recruitment rate (births and immigrations), the finite rate of increase, or the maximum instantaneous rate of increase; omega, either the apparent survival rate (deaths and emigrations) or the equilibrium abundance (carrying capacity); and iota, the number of immigrants per site and year. Estimates of population size at each time period can be derived from these parameters, and thus so can trend estimates. Or, trend can be estimated directly using dynamics="trend".

When immigration is set to FALSE (the default), iota is not modeled. When immigration is set to TRUE and dynamics is set to "autoreg", the model will separately estimate birth rate (gamma) and number of immigrants (iota). When immigration is set to TRUE and dynamics is set to "trend", "ricker", or "gompertz", the model will separately estimate local contributions to population growth (gamma and omega) and number of immigrants (iota).

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.

For "constant", "autoreg", or "notrend" dynamics, the latent abundance state following the initial sampling period arises from a Markovian process in which survivors are modeled as S_{it} \sim Binomial(N_{it-1}, \omega_{it}), and recruits follow G_{it} \sim Poisson(\gamma_{it}). Alternative population dynamics can be specified using the dynamics and immigration arguments.

\lambda_i, \gamma_{it}, and \iota_{it} are modeled using the the log link. p_{ijt} is modeled using the logit link. \omega_{it} is either modeled using the logit link (for "constant", "autoreg", or "notrend" dynamics) or the log link (for "ricker" or "gompertz" dynamics). For "trend" dynamics, \omega_{it} is not modeled.

The detection process is modeled as multinomial: \mathbf{y_{it}} \sim Multinomial(N_{it}, \pi_{it}), where \pi_{ijt} is the multinomial cell probability for plot i at time t on occasion j.

Options for the detection process include equal-interval removal sampling ("removal"), double observer sampling ("double"), or dependent double-observer sampling ("depDouble"). This option is specified when setting up the data using unmarkedFrameMMO. Note that unlike the related functions multinomPois and gmultmix, custom functions for the detection process (i.e., piFuns) are not supported. To request additional options contact the author.

Value

An object of class unmarkedFitMMO

Warning

This function can be extremely slow, especially if there are covariates of gamma or omega. Consider testing the timing on a small subset of the data, perhaps with se=FALSE. Finding the lowest value of K that does not affect estimates will also help with speed.

Note

When gamma or omega are modeled using year-specific covariates, the covariate data for the final year will be ignored; however, they must be supplied.

If the time gap between primary periods is not constant, an M by T matrix of integers should be supplied to unmarkedFrameMMO using the primaryPeriod argument.

Secondary sampling periods are optional, but can greatly improve the precision of the estimates.

Author(s)

Ken Kellner [e-mail], Richard Chandler

References

Dail, D. and L. Madsen (2011) Models for Estimating Abundance from Repeated Counts of an Open Metapopulation. Biometrics. 67: 577-587.

Hostetler, J. A. and R. B. Chandler (2015) Improved State-space Models for Inference about Spatial and Temporal Variation in Abundance from Count Data. Ecology 96: 1713-1723.

Royle, J. A. (2004). Generalized estimators of avian abundance from count survey data. Animal Biodiversity and Conservation 27(1), 375-386.

See Also

multinomPois, gmultmix, unmarkedFrameMMO

Examples

   
  #Generate some data 
  set.seed(123)
  lambda=4; gamma=0.5; omega=0.8; p=0.5
  M <- 100; T <- 5
  y <- array(NA, c(M, 3, T))
  N <- matrix(NA, M, T)
  S <- G <- matrix(NA, M, T-1)

  for(i in 1:M) {
    N[i,1] <- rpois(1, lambda)
    y[i,1,1] <- rbinom(1, N[i,1], p)    # Observe some
    Nleft1 <- N[i,1] - y[i,1,1]         # Remove them
    y[i,2,1] <- rbinom(1, Nleft1, p)   # ...
    Nleft2 <- Nleft1 - y[i,2,1]
    y[i,3,1] <- rbinom(1, Nleft2, p)

    for(t in 1:(T-1)) {
      S[i,t] <- rbinom(1, N[i,t], omega)
      G[i,t] <- rpois(1, gamma)
      N[i,t+1] <- S[i,t] + G[i,t]
      y[i,1,t+1] <- rbinom(1, N[i,t+1], p)    # Observe some
      Nleft1 <- N[i,t+1] - y[i,1,t+1]         # Remove them
      y[i,2,t+1] <- rbinom(1, Nleft1, p)   # ...
      Nleft2 <- Nleft1 - y[i,2,t+1]
      y[i,3,t+1] <- rbinom(1, Nleft2, p)
    }
  }
  y=matrix(y, M)
  
  #Create some random covariate data
  sc <- data.frame(x1=rnorm(100))

  ## Not run: 
  #Create unmarked frame
  umf <- unmarkedFrameMMO(y=y, numPrimary=5, siteCovs=sc, type="removal")

  #Fit model
  (fit <- multmixOpen(~x1, ~1, ~1, ~1, K=30, data=umf))
  
  #Compare to truth
  cf <- coef(fit)
  data.frame(model=c(exp(cf[1]), cf[2], exp(cf[3]), plogis(cf[4]), plogis(cf[5])), 
             truth=c(lambda, 0, gamma, omega, p))

  #Predict
  head(predict(fit, type='lambda'))

  #Check fit with parametric bootstrap
  pb <- parboot(fit, nsims=15)
  plot(pb)

  # Empirical Bayes estimates of abundance for each site / year
  re <- ranef(fit)
  plot(re, layout=c(10,5), xlim=c(-1, 10))
  
## End(Not run)

unmarked documentation built on Sept. 11, 2024, 8:28 p.m.