F.huggins.estim: F.huggins.estim - Estimation of Huggins closed population...

In tmcd82070/MRA: Mark-Recapture Analysis

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. For example: 'capture = ~ age + sex', where age and sex are matrices of size NAN X NS containing the age and sex covariate values. NAN = number of individuals = number of rows in histories matrix (see below). NS = number of samples = number of columns in histories matrix (see below). Number of matrices specified in the initial capture model is assumed to be NX. Specify intercept only model as 'capture = ~ 1'. Specify model without an intercept using 'capture = ~ -1 + x'.
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 the recapture model is NY. Total number of parameters this routine attempts to estimate is NX+NY. See df argument. If recapture=NULL, no recapture model (or the empty model) is estimated. In this case, recapture probabilities equal initial capture probabilities and both depend on the model in capture. Note that NULL models are specified without the ~.
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, 1 in cell (i,j) means individual i was captured on occasion j and released live back into the population. Because the population being sampled is assumed closed, deaths on capture (known removals) are not allowed. If deaths on capture occurred and an estimate of N at the beginning of the study is sought, remove the entire history, estimate N using this routine from the remaining histories, and add back the number of deleted histories.
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 in the model for initial capture also appear in the model for recaptures with the same coefficient values. See Details section. If remove is a vector, each entry specifies whether the corresponding effect in capture should be removed from the recapture model. If remove is shorter than NX (the number of matrices in capture), it is replicated to have length NX.
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. nhat.v.meth = 1 is a delta method estimator utilizing the derivative of P(ever captured) w.r.t. the capture parameters. This is the same estimator as used in program MARK.
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. If df <= 0, the number of parameters will be set to NX+NY = the number of estimated coefficients. Otherwise, if df > 0, the supplied value is used. Only AIC, QAIC, AICc, and QAICc are dependent on this value (in their penalty terms).
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 in help for F.cjs.estim for a plot of the link functions) The "logit" link is eta = log( p / (1 - p) ) with inverse p = 1 / (1 + exp(-eta)). The "sine" link is eta = 8*asin( 2*p - 1 ) / pi , which ranges from -4 to 4. The inverse "sine" link is p = 0.5*(1 + sin( eta*pi/8 )) for values of eta between -4 and 4. For values of eta < -4, p = 0. For values of eta > 4, p = 1. Scaling of the sine link was chosen to yield coefficients roughly the same magnitude as the logit link. The "hazard" link is eta = log( -log( 1 - p )), with inverse 1 - exp( -exp( eta )). The value of p from the inverse hazard link approaches 0 as eta decreases. For values of eta > 3, p = 1 for all intents and purposes.
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, changes to these variables can effect speed and performance for ill-behaved models. See mra.control() for a description of the individual control parameters.

Details

This routine compiles all the covariate matrices, then calls a Fortran routine to maximize the Huggins closed population likelihood. So-called heterogeneous 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[i,j] = B0 + B1*x1[i,j] + ... + Ba*xa[i,j]

and

c[i,j] = B0 + B1*x1[i,j] + ... + Ba*xa[i,j] + G0 + G1*z1[i,j] + ... + Gb*zb[i,j]

where the x's and z's are covariate values specified in capture and recapture, respectively, and the B's and G's are estimated coefficients. (For brevity, 'a' has been substituted 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[i,j] = B0 + B1*x1[i,j] + B2*x2[i,j]

and

c[i,j] = B0 + B2*x2[i,j] + G0 + G1*x1[i,j] + G2*x3[i,j].

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.

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

Author(s)

Trent McDonald, WEST-INC, [e-mail]

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
set.seed(3425)
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)

tmcd82070/MRA documentation built on May 31, 2019, 4:36 p.m.