mra (version 2.1)

F.huggins.estim: F.huggins.estim - Estimation of Huggins closed population capture-recapture model.

Description

Estimates Huggin's closed population capture-recapture models with individual, time, and individual-time varying covariates using the "regression" parameterization of Amstrup et al (2006, Ch 9). For live recaptures only. A logistic link function is used to relate probability of capture to external covariates.

Usage

F.huggins.estim(capture, recapture=NULL, histories, remove=FALSE, cap.init, recap.init,
    nhat.v.meth=1, df=NA, link="logit", control=mra.control())

Arguments

capture

Formula specifying covariates to included in the initial capture probability model. Must be a formula object without a response. Specify ~, followed by the names of 2-D arrays of covariates to relate to initial capture probability.

recapture

Formula specifying covariates to included in the recapture probability model, or NULL. Should be specified the same way as the capture model. For example: 'recapture = ~ behave + sex'. The number of covariates specified in

histories

A NAN X NS = (number of individuals) X (number of capture occasions) matrix containing capture histories. Capture histories are comprised of 0's and 1's only. 0 in cell (i,j) of this matrix means individual i was not captured on occasion j,

remove

A logical scalar, or vector of logical values, specifying which capture covariates to remove from the recapture model. By default (remove=FALSE), no capture covariates are removed, meaning all terms i

cap.init

(optional) Vector of initial values for coefficients in the initial capture model. One element per covariate in capture. This parameter does not usually need to be specified.

recap.init

(optional) Vector of initial values for coefficients in the recapture model. One element per covariate in recapture. This parameter does not usually need to be specified.

nhat.v.meth

Integer specifying method for computing variance estimate for the population size estimate. Currently, only nhat.v.meth = 1 is implemented. Plans are for nhat.v.meth = 2 to be a boot strap estimate of variance. <

df

External (override) model degrees of freedom to use during estimation. If df == NA, the number of parameters is estimated from the rank of the matrix of 2nd derivatives or Hessian, depending on cov.meth parameter.

link

The link function to be used. The link function converts linear predictors in the range (-infinity, infinity) to probabilities in the range (0,1). Valid values for the link function are "logit" (default), "sine", and "hazard". (see Examples

control

A list containing named control parameters for the minimization and estimation process. Control parameters include number of iterations, covariance estimation method, etc. Although the default values work in the vast majority of cases, change

Value

• An object (list) of class c("hug","cr") with many components. Use print.hug to print it nicely. Use names(fit), where the call was fit <- F.huggins.estim(...), to see names of all returned components. To see values of individual components, issue commands like fit$n.hat, fit$se.n.hat, etc. Components of the returned object are as follows:
• histories
  The input capture history matrix. Size NAN x NS
• aux
  Auxiliary information, mostly stored input values. This is a list containing: $call, $nan=number of individuals, $ns=number of samples, $nx=number of coefficients in the initial capture model, $ny=number of coefficients in recapture model, $cov.name=names of the covariates in both models (initial capture covariates first, then recapture covariates), $ic.name=name of capture history matrix, $mra.version=version number of this package, $R.version=R version used, $run.date=date the model was estimated.
• loglik
  Value of the Huggins log likelihood at it's maximum.
• deviance
  Model deviance = -2*loglik. This is relative deviance because all covariates are individual and time varying. When individual covariates are present, a saturated likelihood cannot be computed. Use this to compute deviance differences only.
• aic
  AIC for the model = deviance + 2*(df). df is either the estimated number of independent parameters (by default), or NX+NY, or a specified value, depending on the input value of df parameter.
• aicc
  AIC with small sample correction = AIC + (2*df*(df+1)) / (NAN - df - 1)
• capcoef
  Vector of estimated coefficients in the initial capture model. Length NX.
• se.capcoef
  Vector of estimated standard errors for coefficients in initial capture model. Length NX.
• recapcoef
  Vector of estimated coefficients in the recapture model. Length NY.
• se.surcoef
  Vector of standard errors for coefficients in recapture model. Length NY.
• covariance
  Variance-covariance matrix for the estimated model coefficients. Size (NX+NY) X (NX+NY).
• p.hat
  Matrix of estimated initial capture probabilities computed from the model. Size of this matrix is NAN x NS. Cell (i,j) is estimated probability of first capture for individual i during capture occasion j.
• se.p.hat
  Matrix of standard errors for estimated initial capture probabilities. Size NAN x NS.
• c.hat
  Matrix of estimated recapture probabilities computed from the model. Size NAN x NS. Cell (i,j) is estimated probability of capturing individual i during occasion j given that it was initially captured prior to j.
• se.c.hat
  Matrix of standard errors for estimated recapture probabilities. Size NAN X NS.
• df
  Number of estimable parameters in the model. df is either the estimated number of independent parameters (by default) based on rank of the variance matrix, or NX+NY, or a specified value, depending on the input value of df parameter.
• message
  A string indicating whether the maximization routine converged.
• exit.code
  Exit code from the maximization routine. Interpretation for exit.code is in message.
• cov.code
  A code indicating the method used to compute the covariance matrix.
• cov.meth
  String indicating method used to compute covariance matrix. Interprets cov.code.
• n.hat
  The Huggins estimate of population size. This estimate is sum( 1/ pstar(i) ), where pstar(i) is probability of observing individual i, which equals 1 - p.hat[i,1]*p.hat[i,2]* ... *p.hat[i,NS], where p.hat is the returned value of p.hat.
• se.n.hat
  Estimated standard error of n.hat. Computed using method specified in nhat.v.meth.
• n.hat.lower
  Lower limit of log based confidence interval for n.hat.
• n.hat.upper
  Upper limit of log based confidence interval for n.hat.
• n.hat.conf
  Confidence level for the interval on n.hat. Currently, confidence level cannot be changed from 95%.
• nhat.v.meth
  Code for method used to compute variance of n.hat. Currently, this is 1 only.
• num.caught
  Number of individuals ever captured = number of rows in the histories matrix.
• n.effective
  Effective sample size = number of observed individuals times number of occasions = NAN * NS

Details

This routine compiles all the covariate matrices, then calls a Fortran routine to maximize the Huggins closed population likelihood. So-called heterogenious models that utilize mixture distributions for probability of capture cannot be fitted via this routine. If remove=FALSE (default) the models for initial capture and subsequent recapture are, $$p_{ij} = \beta_0 + \beta_1 x_{ij1} + \ldots + \beta_a x_{ija}$$ and $$c_{ij} = \beta_0 + \beta_1 x_{ij1} + \ldots + \beta_a x_{ija} + \gamma_0 + \gamma_1 z_{ij1} + \ldots + \gamma_b z_{ijb}$$ where the x's and z's are covariate values specified in capture and recapture, respectively, and the $\beta$'s and $\gamma$'s are estimated coefficients. (For breavity, 'a' has been substitutied for NX, 'b' for NY.) In other words, by default all effects in the capture model also appear in the recapture model with the same estimated coefficients. This is done so that capture and recapture probabilities can be constrained to equal one another. If capture=~1 and recapture=NULL, capture and recapture probabilities are constant and equal to one another. If capture=~x1 and recapture=NULL, capture and recapture probabilities are equal, and both are the exact same function of covariate x1. A simple additive behavioral (trap happy or trap shy) effect is fitted by specifying an intercept-only model for recaptures, i.e., capture=~x1+x2+...+xp and recapture=~1. When a Huggin's model object is printed using the default print method (print.hug), a "C" (for "capture") appears next to coefficients in the recapture model that are also in the initial capture model. These coefficients are fixed in the recapture model. A "B" (for "behavioral") appears next to free coefficients in the recapture model that do not appear in the initial capture model. If remove is something other than FALSE, it is extended to have length NX, and if element i equals TRUE, the i-th effect in the capture model is removed from the recapture model. If remove=c(FALSE, TRUE, FALSE), capture=~x1+x2, and recapture=~x1+x3, the models for initial capture and subsequent recapture are, $$p_{ij} = \beta_0 + \beta_1 x_{ij1} + \beta_2 x_{ij2}$$ and $$c_{ij} = \beta_0 + \beta_2 x_{ij2} + \gamma_0 + \gamma_1 x_{ij1} + \gamma_2 x_{ij3}.$$ Note that x1 appears in the recapture equation, but with a different estimated coefficient. If remove=TRUE, all capture effects are removed from the recapture model and the models are completely separate. The ability to remove terms from the recapture model adds modeling flexibility. For example, if initial captures were hypothesized to depend on the variable sex, but recaptures were hypothesized to be constant (no sex effect), the arguments to fit this model would be capture=~sex, recapture=~1, and remove=TRUE. A pure time-varying model with different time effects in the initial and subsequent capture models can be fitted using capture=~tvar(1:ns,nan), recapture=~tvar(1:ns,nan), and remove=TRUE. In this case, the same model, but parameterized differently, can be fitted with remove=FALSE. See Details of help(F.cjs.estim) for ways that 2-d matrices, 1-d vectors, and 1-d factors can be specified in the capture and recapture models. If argument trace in a call to mra.control is set to something other than 0, a log file named mra.log is written to the current directory. See mra.control for actions associated with values of trace. CAREFUL: mra.log is overwritten each run. Values in 2-d Matrix Covariates: Even though covariate matrices are required to be NAN x NS (same size as capture histories), there are not that many recapture parameters. Recapture parameters for the first occasion are not defined. For all covariates in the recapture model, only values in columns 2:ncol(histories) are used. See examples for demonstration.

References

Huggins, R. M. 1989. On the statistical analysis of capture experiments. Biometrika 76:133-140. Amstrup, S. C., T. L. McDonald, and B. F. J. Manly (editors). 2005. Handbook of Capture-Recapture Analysis. Princeton University Press.

See Also

print.hug, F.cjs.estim

Examples

# Fake the data for these examples
ch.mat <- matrix( round(runif(30*5)), nrow=30, ncol=5)
ch.mat <- ch.mat[ apply(ch.mat,1,sum) > 0, ]  # no zero rows allowed
ct <- as.factor(1:ncol(ch.mat))
attr(ct,"nan") <- nrow(ch.mat)   # used to fit time varying factor
sex <- round(runif(nrow(ch.mat)))   # fake sex 
attr(sex,"ns") <- ncol(ch.mat)

# Models parallel to the 8 Otis et al. models.
# see Amstrup et al. (2005, p. 77)

# Constant model (model M(0)).
hug.0 <- F.huggins.estim( ~1, NULL, ch.mat )

# Time varying model (model M(t))
hug.t <- F.huggins.estim( ~tvar(ct), NULL, ch.mat)

# Additive Behavioral model (model M(b))
hug.b <- F.huggins.estim( ~1, ~1, ch.mat )

# Time and Behavioral model (model M(tb))
hug.tb <- F.huggins.estim( ~tvar(ct), ~1, ch.mat )

# Individual effects (model M(h))
hug.h <- F.huggins.estim( ~ivar(sex), NULL, ch.mat )

# Individual and Behavioral effects (model M(bh))
hug.bh <- F.huggins.estim( ~ivar(sex), ~1, ch.mat )

# Individual and time effects (model M(th))
hug.th <- F.huggins.estim( ~ivar(sex)+tvar(ct), NULL, ch.mat )

# Individual, time, and behavoral effects (model M(tbh))
hug.tbh <- F.huggins.estim( ~ivar(sex)+tvar(ct), ~1, ch.mat )

# Time varying initial captures, recaptures are constant and depend on sex.
hug.custom1 <- F.huggins.estim( ~tvar(ct), ~ivar(sex), ch.mat, remove=TRUE )

# Compare hug.custom1 to the following: Time varying initial captures with 
# time varying recaptures that depend on sex.
hug.custom2 <- F.huggins.estim( ~tvar(ct), ~ivar(sex), ch.mat, remove=FALSE )

# Values in first column of recapture covariates do not matter. 
# Below, mod.1 and mod.2 are identical.
mod.1 <- F.huggins.estim( ~tvar(ct), ~tvar( c( 0,1,2,3,4), nrow(ch.mat)), ch.mat, remove=TRUE)
mod.2 <- F.huggins.estim( ~tvar(ct), ~tvar( c(-9,1,2,3,4), nrow(ch.mat)), ch.mat, remove=TRUE)