The closedpCI.t and closedpCI.0 functions fit a loglinear model specified by the user
and compute a confidence interval for the abundance estimation. For a normal heterogeneous model,
a log-transformed confidence interval (Chao 1987) is produced.
For any other model, the multinomial profile likelihood confidence interval (Cormack 1992) is produced.
The model is identified with the argument m or mX.
For heterogeneous models, the form of the heterogeneity is specified with the arguments
h and h.control. If h is given with mX, heterogeneity is added in mX.
These functions extend closedp.t and closedp.0 as
they broaden the range of models one can fit and they compute confidence intervals.
Unlike the closedp functions, it fits only one model at a time.
closedpCI.t(X, dfreq=FALSE, m=c("M0","Mt","Mh","Mth"), mX=NULL,
h=NULL, h.control=list(), mname=NULL, alpha=0.05,
fmaxSupCL=3, …)closedpCI.0(X, dfreq=FALSE, dtype=c("hist","nbcap"), t=NULL, t0=NULL,
m=c("M0","Mh"), mX=NULL, h=NULL, h.control=list(),
mname=NULL, alpha=0.05, fmaxSupCL=3, …)
plotCI(x.closedpCI, main="Profile Likelihood Confidence Interval", …)
# S3 method for closedpCI
print(x, …)
# S3 method for closedpCI
boxplot(x, main="Boxplots of Pearson Residuals", …)
# S3 method for closedpCI
plot(x, main="Scatterplot of Pearson Residuals", …)
The matrix of the observed capture histories (see Rcapture-package for a
description of the accepted formats).
A logical. By default FALSE, which means that X has one row per unit.
If TRUE, it indicates that the matrix X contains frequencies in its last column.
A characters string, either "hist" or "nbcap", to specify the type of data.
"hist", the default, means that X contains complete observed capture histories.
"nbcap" means that X contains numbers of captures (see Rcapture-package
for details on data formats).
Requested only if dtype="nbcap". A numeric specifying the total number of
capture occasions in the experiment. For closedpCI.0, the value t=Inf
is accepted. It indicates that captures occur in continuous time
(see Rcapture-package).
A numeric. Models are fitted considering only the frequencies of units captured 1 to t0 times.
By default, if t is not equal to Inf, t0=t. When t=Inf, the default value
of t0 is the larger number of captures observed.
A character string indicating the model to fit. For closedpCI.0 it can be either "M0"=M0 model
or "Mh"=Mh model. For closedpCI.t it can also be "Mt"=Mt model or "Mth"=Mth model.
The design matrix of the loglinear model. By default, the design matrix is built based
on the m argument. If a mX argument is given, it must be a matrix (or
an object that can be coerced to a matrix by as.matrix). The order of the capture
histories in the mX matrix must be as defined in the histpos.t or
histpos.0 function.
For the closedpCI.t function only, mX can also be an object of class "formula" or a hierarchical loglinear model name. The details of a formula specification are given under Details and the details of a hierarchical loglinear model name specification are given in the closedpMS.t documentation.
An intercept is always added to the model. Therefore, mX must not contain a column of ones if it is a matrix. Also, if the argument h is not NULL, one or more columns for heterogeneity are added to the design matrix.
A character string ("LB", "Chao", "Poisson", "Darroch", "Gamma" or "Normal") or a
numerical R function specifying the form of the column(s) for heterogeneity
in the design matrix.
"LB" and "Chao" represents Chao's lower bound model (the default),
"Poisson" represents the function \(f(k)=theta^k-1\),
"Darroch" represents the function \(f(k)=k^2/2\),
and "Gamma" represents the function \(f(k)=-\log(theta + k) + \log(theta)\)
where \(k\) is the number of captures and \(theta\) is a parameter specified with
the h.control element theta.
"Normal" represents heterogeneity
modeled with a normal mixing distribution, as presented in Rivest (2011).
If an R function is given, it must be the implementation of any convex
mathematical function \(f(k)\) (it must have only one argument).
A list of elements to control the heterogeneous part of the model, if any. For a Poisson or Gamma heterogeneous model:
thetaThe value of the parameter \(theta\) in \(f(k)=theta^k-1\) for the Poisson model (the default value is 2) and in \(f(k)=-\log(theta + k) + \log(theta)\) for the Gamma model (the default value is 3.5).
For a Chao's lower bound heterogeneous model:
negIf this option is set to TRUE (the default), negative eta parameters in
Chao's lower bound models are set to zero.
For a Normal heterogeneous model:
A character string specifying the name of the customized model. By default, it is derived from the arguments specifying the model.
A confidence interval with confidence level 1-alpha is constructed. The value of alpha
must be between 0 and 1; the default is 0.05.
A numeric (by default 3). The upper end point of the interval to be searched by
uniroot to find the upper bound of the multinomial profile likelihood
confidence interval (Cormack 1992) is defined by
fmaxSupCL\(\times \hat{N}\). In
this product, fmaxSupCL is a factor multiplying \(\hat{N}\), the population
size estimated by Poisson regression. If the upper bound obtained is equal
the upper end point fmaxSupCL\(\times \hat{N}\), the
element SupCL of the output value CI will begin with the symbol
\(>\). In that case, one could try to increase the value of fmaxSupCL
to find the upper bound of the multinomial profile likelihood confidence interval
instead of a minimum value for this upper bound.
Further arguments to be passed to glm, optim,
print.default, plot.default or boxplot.default.
An object, produced by a closedpCI function, to produce a plot
of the multinomial profile likelihood for \(N\).
A main title for the plot
An object, produced by a closedpCI function, to print or to plot.
The number of captured units
The total number of capture occasions in the data matrix X.
When captures occur in continuous time (input argument t=Inf),
this output value is the larger number of captures observed.
For closedpCI.0 only: the value of the argument t0 used in the computations.
A table containing, for the fitted model:
abundance:the estimated population size,
stderr:the standard error of the estimated population size,
deviance:the model's deviance,
df:the number of degrees of freedom,
AIC:the Akaike's information criterion,
BIC:the bayesian information criterion,
infoFit:a numerical code giving information about error or warnings
encountered when fitting the model (see Rcapture-package for details).
Not produced for normal heterogeneous models (h="Normal"): The asymptotic bias of the estimated population size.
The 'glm' object obtained from fitting the model except for
normal heterogeneous models (h="Normal"). These models are not fitted
with glm. For them, fit is a list with the following elements:
parameters:the matrix of parameters (loglinear coefficients + sigma parameter) estimates with their standard errors,
varcov:the estimated variance-covariance matrix of the estimated parameters,
y:the y vector used to fit the model,
fitted.values:the model fitted values,
initparam:the initial values for the parameters (loglinear coefficients +
sigma parameter) used by optim,
optim:the output produced by optim.
A vector of character strings. If the glm or optim
(for normal heterogeneous models) function generates
one or more warnings when fitting the model, a copy of these warnings are
stored in fit.warn. NULL if no warnings occured.
For Chao's lower bound models only: the position of the eta parameters set to zero in the loglinear parameter vector, if any.
Not produced for normal heterogeneous models (h="Normal"): A table containing
the abundance estimation and its multinomial profile likelihood confidence interval.
The last column of the table, named infoCI, contains a numerical code giving information about error or warnings encountered when calculating the confidence interval. Here is a description of the meaning of this numerical code:
no error or warning occured;
an error occured while calculating the multinomial abundance estimation;
a warning occured while calculating the multinomial abundance estimation;
a warning occured while calculating the lower bound of the multinomial profile likelihood confidence interval;
a warning occured while calculating the upper bound of the multinomial profile likelihood confidence interval.
Not produced for normal heterogeneous models (h="Normal"): If
an error occured while calculating the multinomial abundance, a copy of the error message is
stored in CI.err. NULL if no error occured.
Not produced for normal heterogeneous models (h="Normal"): If
one or more warnings occur while calculating the multinomial abundance estimation or the profile likelihood confidence interval, a copy of these warnings are
stored in CI.warn. NULL if no warnings occured.
1-the confidence level of the interval.
Not produced for normal heterogeneous models (h="Normal"):
The x-coordinates for plot.closedpCI.t.
Not produced for normal heterogeneous models (h="Normal"):
The y-coordinates for plot.closedpCI.t.
The closedpCI.t function fits models using the frequencies of the observable capture histories (vector of size \(2^t-1\)), whereas closedpCI.0 uses the number of units capture i times, for \(i=1,\ldots,t\) (vector of size \(t\)). Thus, closedpCI.0 can be used with data sets larger than those for closedpCI.t, but it cannot fit models with a temporal effect. See Rcapture-package for more details about the distinction between .t and .0 functions.
These functions do not work for closed population models featuring a behavioral effect, such as Mb and Mbh. The abundance estimation is calculated as the number of captured units plus the exponential of the Poisson regression intercept. However, models with a behavioral effect can by fitted with closedp.t (Mb and Mbh), closedp.Mtb and closedp.bc.
CHAO'S LOWER BOUND MODELS
Chao's lower bound models estimate a lower bound for the abundance. Rivest (2011)
presents a generalized loglinear model underlying this estimator. To test whether a
certain model for heterogeneity is adequate,
one can conduct a likelihood ratio test by subtracting the deviance of a Chao's lower
bound model to the deviance of the heterogeneous model under study. The two models should
have the same mX argument.
Under the null hypothesis of equivalence between the two models, the difference of deviances
follows a chi-square distribution with degrees of freedom equal to the difference between
the models' degrees of freedom.
A Chao's lower bound model contains \(t-2\) parameters, called
eta parameters, for the heterogeneity. These parameters should theoretically be greater
or equal to zero (see Rivest and Baillargeon (2007)). When the element neg of
the argument h.control is set to TRUE (the default), negative eta parameters are
set to zero (to do so, columns are removed from the design matrix of the model). Degrees
of freedom of Chao's lower bound model increase when eta parameters are set to zero.
ARGUMENT mX : FORMULA SPECIFICATION
For the closedpCI.t function, mX
can be an object of class "formula". The only accepted variables
in this formula are c1 to ct. The variable ci represents
a capture indicator (1 for a capture, 0 otherwise) for the \(i\)th capture occasions.
Also, the formula must not contain a response variable since
it is only used to construct the design matrix of the model.
For example, if t=3, the Mt model is fitted if
mX = ~ . or mX = ~ c1 + c2 + c3. The symbol . in this formula
is a shortcut for c1 + c2 + ... + ct. Formula mX arguments
facilitate the addition of interactions between capture occasions in the model. For
example, if t=3, the Mt model with an interaction between the first and
the second capture occasion is fitted if mX = ~ . + c1:c2.
See formula for more details of allowed formulae.
PLOT METHODS AND FUNCTIONS
The boxplot.closedpCI function produces a boxplot of the Pearson residuals of the customized model.
The plot.closedpCI function traces the scatterplot of the Pearson residuals in terms of \(f_i\)
(number of units captured i times) for the customized model.
The plotCI function produces a plot of the multinomial profile likelihood for \(N\).
The value of N maximizing the profile likelihood and the bounds of the confidence interval are identified.
Baillargeon, S. and Rivest, L.P. (2007) Rcapture: Loglinear models for capture-recapture in R. Journal of Statistical Software, 19(5), 10.18637/jss.v019.i05.
Chao, A. (1987) Estimating the population size for capture-recapture data with unequal catchability. Biometrics, 43(4), 783--791.
Cormack, R. M. (1992) Interval estimation for mark-recapture studies of closed populations. Biometrics, 48, 567--576.
Rivest, L.P. (2011) A lower bound model for multiple record systems estimation with heterogeneous catchability. The International Journal of Biostatistics, 7(1), Article 23.
Rivest, L.P. and Baillargeon, S. (2007) Applications and extensions of Chao's moment estimator for the size of a closed population. Biometrics, 63(4), 999--1006.
# NOT RUN {
# hare data set
CI <- closedpCI.t(hare, m = "Mth", h = "Poisson", h.control = list(theta = 2))
CI
plotCI(CI)
# HIV data set
mat <- histpos.t(4)
mX2 <- cbind(mat, mat[, 1] * mat[ ,2])
closedpCI.t(HIV, dfreq = TRUE, mX = mX2, mname = "Mt interaction 1,2")
# which can be obtained more conveniently with
closedpCI.t(HIV, dfreq = TRUE, mX = ~ . + c1:c2, mname = "Mt interaction 1,2")
# BBS2001 data set
CI0 <- closedpCI.0(BBS2001, dfreq = TRUE, dtype = "nbcap", t = 50, t0 = 20,
m = "Mh", h = "Gamma", h.control = list(theta = 3.5))
CI0
plot(CI0)
plotCI(CI0)
### As an alternative to a gamma model, one can fit a negative Poisson model.
### It is appropriate in experiments where very small capture probabilities
### are likely. It can lead to very large estimators of abundance.
# Third primary period of mvole data set
period3 <- mvole[, 11:15]
psi <- function(x) { 0.5^x - 1 }
closedpCI.t(period3, m = "Mh", h = psi)
### Example of normal heterogeneous models
### diabetes data of Bruno et al. (1994)
histpos <- histpos.t(4)
diabetes <- cbind(histpos, c(58,157,18,104,46,650,12,709,14,20,7,74,8,182,10))
# chosen interaction set I in Rivest (2011)
closedpCI.t(X = diabetes, dfreq = TRUE, mX = ~ . + c1:c3 + c2:c4 + c3:c4,
h = "Normal", mname = "Mth normal with I")
### Example of captures in continuous time
# Illegal immigrants data set
closedpCI.0(ill, dtype = "nbcap", dfreq = TRUE, t = Inf, m = "Mh", h = "LB")
# }
Run the code above in your browser using DataCamp Workspace