F.huggins.estim(capture, recapture=NULL, histories, remove=FALSE, cap.init, recap.init, nhat.v.meth=1, df=NA, link="logit", control=mra.control())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'. 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 probabilties and both depend on the
model in capture. Note that NULL models are specified without the ~.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 matricies in
capture), it is replicated to have length NX.capture. This parameter does not usually need to be specified. recapture. This parameter does not usually need to be specified. 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 == 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" estimated="" coefficients.="" otherwise,="" if="" F.cjs.estim
for a plot of the link functions)
mra.control() for a description of the individual control parameters.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:
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.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.df*(df+1))
/ (NAN - df - 1)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.exit.code is in message.cov.code.p.hat.n.hat. Computed using method
specified in nhat.v.meth.n.hat. n.hat.n.hat.
Currently, confidence level cannot be changed
from 95%. n.hat. Currently,
this is 1 only.histories matrix.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 $B$'s and
$G$'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.
Amstrup, S. C., T. L. McDonald, and B. F. J. Manly (editors). 2005. Handbook of Capture-Recapture Analysis. Princeton University Press.
print.hug,
F.cjs.estim
# 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)
Run the code above in your browser using DataLab