Description Usage Arguments Details Value Note Author(s) References See Also Examples
Estimates CormackJollySeber (CJS) capturerecapture models with individual, time, and individualtime varying covariates using the "regression" parametrization of Amstrup et al (2005, Ch 9). For live recaptures only. Losses on capture allowed. Uses a logistic link function to relate probability of capture and survival to external covariates.
1 2 3 
capture 
Formula specifying the capture probability model. Must be a formula object with
no response. I.e., "~" followed by the names of 2D arrays of covariates to fit in the capture model.
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 animals = number of rows in

survival 
Formula specifying the survival probability model. Must be a formula object with
no response. I.e., "~" followed by the names of 2D arrays of covariates to fit in the survival model.
For example: 'survival = ~ year + ageclass' where year and ageclass are matrices of size NAN X NS
containing year and ageclass covariate values. Number of matrices specified in the survival
model is assumed to be NY.
Time varying and individual varying vectors are fitted using 
histories 
A NAN X NS = (number of animals) X (number of capture occasions) matrix containing capture histories. Capture histories are comprised of 0's, 1',s and 2's. 0 in cell (i,j) means animal i was not captured on occasion j, 1 in cell (i,j) means animal i was captured on occasion j and released live back into the population, 2 in cell (i,j) means animal i was captured on occasion j and was not released back into the population (e.g., it died). Animals with '2' as the last nonzero entry of their history are considered 'censored'. Their lack of capture information is removed from the likelihood after the occasion with the 2. Rows of all zeros (i.e., no captures) are allowed in the history matrix, but do not affect coefficient or population size estimates. A warning is thrown if rows of all zeros exist. Capture and survival probabilities are computed for animals with all zero histories. In this way, it is possible to have the routine compute capture or survival estimates for combinations of covariates that do not exist in the data by associating the covariate combinations with histories that have all zero entries. 
cap.init 
(optional) Vector of initial values for coefficients in the capture model. One element
per covariate in 
sur.init 
(optional) Vector or initial values for coefficients in the survival model. One element
per covariate in 
group 
(optional) A vector of length NAN giving the (nonchanging) group membership of every
captured animal (e.g., sex). Group is used only for computing TEST 2 and TEST 3.
TEST 2 and TEST 3 are computed separately for each group. E.g., if group=sex,
TEST 2 and TEST 3 are computed for each sex. TEST 2 and TEST3 are used only
to estimate Chat. See 
nhat.v.meth 
Integer specifying method for computing variance estimates
of population size estimates. 
c.hat 
External (override) estimate of variance inflation factor ( 
df 
External (override) model degrees of freedom to use during estimation.
If 
intervals 
Time intervals. This is a vector of length 
conf 
Confidence level for the confidence intervals placed around estimates of population size. Default 95% confidence. 
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 for a plot of the link functions)

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 illbehaved models. See

This is the workhorse routine for estimating CJS models. It compiles all the covariate matrices, then calls a Fortran routine to maximize the CJS likelihood and perform goodnessoffit tests. HorvitzThompsontype population size estimates are also computed by default.
If control=mra.control(trace=1)
, a log file, named mra.log
, is written to the current directory. This file contains
additional details, such as individual Test 2 and Test 3 components, in a semifriendly
format. This file is overwritten each run. See help(mra.control)
for more details.
Model Specification: Both the capture
and survival
model can be
specified as any combination of 2d matrices (time and individual varying covariates),
1d time varying vectors, 1d individual
varying vectors, 1d time varying factors, and 1d individual varying factors.
Specification of time or individual varying effects uses the
tvar
(for 'time varying') and ivar
(for 'individual varying') functions.
These functions expand covariate vectors along the appropriate dimension to
be 2d matrices suitable for fitting in the model. ivar
expands
an individual varying vector to all occasions. tvar
expands a time
varying covariate to all individuals. To do the expansion, both tvar
and ivar
need to know the size of the 'other' dimension. Thus, tvar(x,100)
specifies a 2d matrix with size 100
by length(x)
.
ivar(x,100)
specifies a 2d matrix with size length(x)
by 100
.
For convenience, the 'other' dimension of time or individual varying covariates
can be specified as an attribute of the vector. Assuming x
is a NS vector and the 'nan' attribute of x
has been set as
attr(x,"nan") < NAN
, tvar(x,NAN)
and
tvar(x)
are equivalent. Same, but viseversa, for individual varying
covariates (i.e.,
assign the number of occasions using attr(x,"ns")<NS
). This saves
some typing in model specification.
Factors are allowed in ivar
and tvar
. When a factor is specified,
the contr.treatment
coding is used. By default, an intercept
is assumed and the first level of all
factors are dropped from the model (i.e., first levels are the reference levels,
the default R action). However, there are applications where more than one level
will need to be dropped, and the user has control over this via the drop.levels
argument to ivar
and tvar
. For example,
tvar(x,drop.levels=c(1,2))
drops the first 2 levels of factor x.
tvar(x,drop.levels=length(levels(x)))
does the SAS thing and drops
the last level of factor x
. If drop.levels
is outside the range
[1,length(levels(x))
] (e.g., negative or 0), no levels of the factor
are dropped. If no intercept is fitted in the model, this results in the
socalled cell means coding for factors.
Example model specifications: Assume 'age' is a NAN x NS 2d matrix of ages, 'effort' is a size NS 1d vector of efforts, and 'sex' is a size NAN 1d factor of sex designations ('M' and 'F').
capture= ~ 1 : constant effect over all individuals and time (intercept only model)
capture= ~ age : Intercept plus age
capture= ~ age + tvar(effort,NAN) : Intercept plus age plus effort
capture= ~ age + tvar(effort,NAN) + ivar(sex,NS) : Intercept plus age plus effort plus sex. Females (1st level) are the reference.
capture= ~ 1 + ivar(sex,NS,0) : sex as a factor, cell means coding
capture= ~ tvar(as.factor(1:ncol(histories)),nrow(histories),c(1,2)) : time varying effects
Values in 2d Matrix Covariates: Even though covariate matrices are required to be NAN x NS (same size as capture histories), there are not that many parameters. The first capture probability cannot be estimated in CJS models, and the NSth survival parameter does not exist. When a covariate matrix appears in the capture model, only values in columns 2:ncol(histories) are used. When a covariate matrix appears in the survival model, only values in columns 1:(ncol(histories)1) are used. See examples for demonstration.
An object (list) of class c("cjs","cr") with many components. Use print.cr
to print
it nicely. Use names(fit)
, where the call was fit < F.cr.estim(...)
,
to see names of all returned components. To see values of individual components,
issue commands like fit$s.hat, fit$se.s.hat, fit$n.hat, etc.
Components of the returned object are as follows:
histories 
The input capture history matrix. 
aux 
Auxiliary information about the fit, mostly stored input values. This is a list containing the following components:

loglik 
Maximized CJS likelihood value for the model 
deviance 
Model deviance = 2* 
aic 
AIC for the model = 
qaic 
QAIC (quasiAIC) = ( 
aicc 
AIC with small sample correction = AIC + (2* 
qaicc 
QAIC with small sample correction = QAIC + (2* 
vif 
Variance inflation factor used = estimate of c.hat = 
chisq.vif 
Composite Chisquare statistic from Test 2 and Test 3 used to compute 
vif.df 
Degrees of freedom for composite chisquare statistic from Test 2 and Test 3, based on pooling rules. 
parameters 
Vector of all coefficient estimates, NX capture probability coefficients first, then NY survival coefficients. This vector is length NX+NY regardless of estimated DF. 
se.param 
Standard error estimates for all coefficients. Length NX+NY. 
capcoef 
Vector of coefficients in the capture model. Length NX. 
se.capcoef 
Vector of standard errors for coefficients in capture model. Length NX. 
surcoef 
Vector of coefficients in the survival model. Length NY. 
se.surcoef 
Vector of standard errors for coefficients in survival model. Length NY. 
covariance 
Variancecovariance matrix for the estimated model coefficients. Size (NX+NY) X (NX+NY). 
p.hat 
Matrix of estimated capture probabilities computed from the model. One for each animal each occasion. Cell (i,j) is estimated capture probability for animal i during capture occasion j. Size NAN X NS. First column corresponding to first capture probability is NA because cannot estimate P1 in a CJS model. 
se.p.hat 
Matrix of standard errors for estimated capture probabilities. One for each animal each occasion. Size NAN X NS. First column is NA. 
s.hat 
Matrix of estimated survival probabilities computed from the model. One for each animal each occasion. Size NAN X NS. Cell (i,j) is estimated probability animal i survives from occasion j to j+1. There are only NS1 intervals between occasions. Last column corresponding to survival between occasion NS and NS+1 is NA. 
se.s.hat 
Matrix of standard errors for estimated survival probabilities. Size NAN X NS. Last column is NA. 
df 
The number of parameters assumed in the model. This value was used in the penalty term of AIC, AICc, QAIC, and QAICc.
This value is either the number of independent
parameters estimated from the rank of the variancecovariance matrix (by default),
or NX+NY, or a specified value, depending on the input value of DF parameter. See 
df.estimated 
The number of parameters estimated from the rank of the variancecovariance matrix. This
is stored so that 
control 
A list containing the input maximization and estimation control parameters. 
message 
A vector of strings interpreting various codes about the estimation.
The messages interpret, in this order, the codes for (1) maximization algorithm used,
(2) exit code from the maximization algorithm (interprets 
exit.code 
Exit code from the maximization routine. Interpretation for

cov.code 
A code indicating the method used to compute the covariance matrix. 
fn.evals 
The number of times the likelihood was evaluated prior to exit from the minimization routine.
If 
ex.time 
Execution time for the maximization routine, in minutes. This is returned for
2 reasons. First, this is useful for benchmarking. Second, in conjunction with

n.hat 
Vector of HorvitzThompson estimates of population size. The HorvitzThompson estimator of size is, N(ij) = sum( h(ij) / p(ij)) Length of 
se.n.hat 
Estimated standard errors for 
n.hat.lower 
Lower limit of 
n.hat.upper 
Upper limit of 
n.hat.conf 
Confidence level of intervals on 
nhat.v.meth 
Code for method used to compute variance of 
num.caught 
Vector of observed number of animals captured each occasion. Length NS. 
fitted 
Matrix of fitted values for the capture histories. Size NAN X NS.
Cell (i,j)
is expected value of capture indicator in cell (i,j) of 
residuals 
Matrix of Pearson residuals defined as, (h(ij)  Psi(ij))^2 / Psi(ij) ,
where Psi(ij) is the expected (or fitted) value for cell
(i,j) and h(ij) is the capture indicator for
animal i at occasion j. This matrix has size NAN X NS. See parts pertaining to the "overall test" in
documentation for 
resid.type 
String describing the type of residuals computed. Currently, only Pearson residuals are returned. 
MARK Users: Due to differences in the way MRA and MARK parameterize the sine link, coefficient estimates will differ between the two packages when this link is used to fit the same model in both packages. The fit (measured by deviance, AIC, etc.) will agree between the two packages. Capture and survival probability estimates will also agree between the two packages.
MARK does not contain a hazard rate link function.
Trent McDonald, WESTINC, [email protected]
Taylor, M. K., J. Laake, H. D. Cluff, M. Ramsay, and F. Messier. 2002. Managing the risk from hunting for the Viscount Melville Sound polar bear population. Ursus 13:185202.
Manly, B. F. J., L. L. McDonald, and T. L. McDonald. 1999. The robustness of markrecapture methods: a case study for the northern spotted owl. Journal of Agricultural, Biological, and Environmental Statistics 4:78101.
Huggins, R. M. 1989. On the statistical analysis of capture experiments. Biometrika 76:133140.
Amstrup, S. C., T. L. McDonald, and B. F. J. Manly (editors). 2005. Handbook of CaptureRecapture Analysis. Princeton University Press.
Peterson. 1986. Statistics and Probability Letters. p.227.
McDonald, T. L., and S. C. Amstrup. 2001. Estimation of population size using open capturerecapture models. Journal of Agricultural, Biological, and Environmental Statistics 6:206220.
tvar
, ivar
, print.cjs
, residuals.cjs
, plot.cjs
,
F.cjs.covars
, F.cjs.gof
, mra.control
, F.update.df
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36  ## Fit CJS model to dipper data, timevarying capture and survivals.
## Method 1 : using factors
data(dipper.histories)
ct < as.factor( paste("T",1:ncol(dipper.histories), sep=""))
attr(ct,"nan")<nrow(dipper.histories)
dipper.cjs < F.cjs.estim( ~tvar(ct,drop=c(1,2)), ~tvar(ct,drop=c(1,6,7)), dipper.histories )
## Method 2 : same thing using 2d matrices
xy < F.cjs.covars( nrow(dipper.histories), ncol(dipper.histories) )
# The following extracts 2D matrices of 0s and 1s
for(j in 1:ncol(dipper.histories)){ assign(paste("x",j,sep=""), xy$x[,,j]) }
dipper.cjs < F.cjs.estim( ~x3+x4+x5+x6+x7, ~x2+x3+x4+x5, dipper.histories )
## Values in the 1st column of capture covariates do not matter
x3.a < x3
x3.a[,1] < 999
dipper.cjs2 < F.cjs.estim( ~x3.a+x4+x5+x6+x7, ~x2+x3+x4+x5, dipper.histories )
# compare dipper.cjs2 to dipper.cjs
## Values in the last column of survival covariates do not matter
x3.a < x3
x3.a[,ncol(dipper.histories)] < 999
dipper.cjs2 < F.cjs.estim( ~x3+x4+x5+x6+x7, ~x2+x3.a+x4+x5, dipper.histories )
# compare dipper.cjs2 to dipper.cjs
## A plot to compare the link functions
sine.link < function(eta){ ifelse( eta < 4, 0, ifelse( eta > 4, 1, .5*(1+sin(eta*pi/8)))) }
eta < seq(5,5, length=40)
p1 < 1 / (1 + exp(eta))
p2 < sine.link(eta)
p3 < 1.0  exp( exp( eta ))
plot(eta, p1, type="l" )
lines(eta, p2, col="red" )
lines(eta, p3, col="blue" )
legend( "topleft", legend=c("logit", "sine", "hazard"), col=c("black", "red", "blue"), lty=1)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.