nmixTTD {unmarked}

Fit N-mixture Time-to-detection Models

Description

Fit N-mixture models with time-to-detection data.

Usage

nmixTTD(stateformula= ~1, detformula = ~1, data, K=100,
    mixture = c("P","NB"), ttdDist = c("exp", "weibull"), starts, method="BFGS", 
    se=TRUE, engine = c("C", "R"), threads = 1, ...)

Arguments

stateformula

Right-hand sided formula for the abundance at each site.

detformula

Right-hand sided formula for mean time-to-detection.

data

unmarkedFrameOccuTTD object that supplies the data (see unmarkedFrameOccuTTD). Note that only single-season models are supported by nmixTTD.

K

The upper summation index used to numerically integrate out the latent abundance. This should be set high enough so that it does not affect the parameter estimates. Computation time will increase with K.

mixture

String specifying mixture distribution: "P" for Poisson or "NB" for negative binomial.

ttdDist

Distribution to use for time-to-detection; either "exp" for the exponential, or "weibull" for the Weibull, which adds an additional shape parameter k.

starts

optionally, initial values for parameters in the optimization.

method

Optimization method used by optim.

se

logical specifying whether or not to compute standard errors.

engine

Either "C" or "R" to use fast C++ code or native R code 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 model extends time-to-detection (TTD) occupancy models to estimate site abundance using data from single or repeated visits. Latent abundance can be modeled as Poisson (mixture="P") or negative binomial (mixture="NB"). Time-to-detection can be modeled as an exponential (ttdDist="exp") or Weibull (ttdDist="weibull") random variable with rate parameter \lambda and, for the Weibull, an additional shape parameter k. Note that occuTTD puts covariates on \lambda and not 1/\lambda, i.e., the expected time between events.

Assuming that there are N independent individuals at a site, and all individuals have the same individual detection rate, the expected detection rate across all individuals \lambda is equal to the the individual-level detection rate r multipled by the number of individuals present N.

In the case where there are no detections before the maximum sample time at a site (surveyLength) is reached, we are not sure if the site has N=0 or if we just didn't wait long enough for a detection. We therefore must censor (C the exponential or Weibull distribution at the maximum survey length, Tmax. Thus, assuming true abundance at site i is N_i, and an exponential distribution for the TTD y_i (parameterized with the rate), then:

y_i \sim Exponential(r_i * N_i) C(Tmax)

Note that when N_i = 0, the exponential rate lambda = 0 and the scale is therefore 1 / 0 = Inf, and thus the value will be censored at Tmax.

Because in unmarked values of NA are typically used to indicate missing values that were a result of the sampling structure (e.g., lost data), we indicate a censored y_i in nmixTTD instead by setting y_i = Tmax_i in the y matrix provided to unmarkedFrameOccuTTD. You can provide either a single value of Tmax to the surveyLength argument of unmarkedFrameOccuTTD, or provide a matrix, potentially with a unique value of Tmax for each value of y. Note that in the latter case the value of y that will be interpreted by nmixTTD as a censored observation (i.e., Tmax) will differ between observations!

Value

unmarkedFitNmixTTD object describing model fit.

Author(s)

Ken Kellner [e-mail]

References

Strebel, N., Fiss, C., Kellner, K. F., Larkin, J. L., Kery, M., & Cohen, J (2021). Estimating abundance based on time-to-detection data. Methods in Ecology and Evolution 12: 909-920.

See Also

unmarked, unmarkedFrameOccuTTD

Examples

## Not run: 

# Simulate data
M = 1000 # Number of sites
nrep <- 3 # Number of visits per site
Tmax = 5 # Max duration of a visit
alpha1 = -1 # Covariate on rate
beta1 = 1 # Covariate on density
mu.lambda = 1 # Rate at alpha1 = 0
mu.dens = 1 # Density at beta1 = 0

covDet <- matrix(rnorm(M*nrep),nrow = M,ncol = nrep) #Detection covariate
covDens <- rnorm(M) #Abundance/density covariate
dens <- exp(log(mu.dens) + beta1 * covDens)
sum(N <- rpois(M, dens)) # Realized density per site
lambda <- exp(log(mu.lambda) + alpha1 * covDet) # per-individual detection rate
ttd <- NULL
for(i in 1:nrep) {
  ttd <- cbind(ttd,rexp(M, N*lambda[,i]))  # Simulate time to first detection per visit
}
ttd[N == 0,] <- 5 # Not observed where N = 0; ttd set to Tmax
ttd[ttd >= Tmax] <- 5 # Crop at Tmax

#Build unmarked frame
umf <- unmarkedFrameOccuTTD(y = ttd, surveyLength=5,
                            siteCovs = data.frame(covDens=covDens),
                            obsCovs = data.frame(covDet=as.vector(t(covDet))))

#Fit model
fit <- nmixTTD(~covDens, ~covDet, data=umf, K=max(N)+10)

#Compare to truth
cbind(coef(fit), c(log(mu.dens), beta1, log(mu.lambda), alpha1))

#Predict abundance/density values
head(predict(fit, type='state'))


## End(Not run)