closedpCI {Rcapture}

Customization of a Loglinear Model and Confidence Interval for Abundance Estimation in Closed Population Capture-Recapture Experiments

Description

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.

Usage

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 class 'closedpCI'
print(x, ...)

## S3 method for class 'closedpCI'
boxplot(x, main="Boxplots of Pearson Residuals", ...)

## S3 method for class 'closedpCI'
plot(x, main="Scatterplot of Pearson Residuals", ...)

Arguments

X

The matrix of the observed capture histories (see Rcapture-package for a description of the accepted formats).

dfreq

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.

dtype

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).

t

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).

t0

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.

m

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.

mX

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.

h

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).

h.control

A list of elements to control the heterogeneous part of the model, if any.
For a Poisson or Gamma heterogeneous model:

theta

The 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:

neg

If 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:

initcoef

Initial values for the loglinear coefficients given to optim. The default values are taken from Darroch's model.

initsig

Initial value for the sigma parameter of the normal mixing distribution. The default is 0.2.

method

The method to be used by optim. The default is "BFGS".

mname

A character string specifying the name of the customized model. By default, it is derived from the arguments specifying the model.

alpha

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.

fmaxSupCL

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.

x.closedpCI

An object, produced by a closedpCI function, to produce a plot of the multinomial profile likelihood for N.

main

A main title for the plot

x

An object, produced by a closedpCI function, to print or to plot.

Details

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 ith 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.

Value

n

The number of captured units

t

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.

t0

For closedpCI.0 only: the value of the argument t0 used in the computations.

results

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).

bias

Not produced for normal heterogeneous models (h="Normal"): The asymptotic bias of the estimated population size.

fit

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.

fit.warn

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.

neg.eta

For Chao's lower bound models only: the position of the eta parameters set to zero in the loglinear parameter vector, if any.

CI

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:

0

no error or warning occured;

-1

an error occured while calculating the multinomial abundance estimation;

4

a warning occured while calculating the multinomial abundance estimation;

5

a warning occured while calculating the lower bound of the multinomial profile likelihood confidence interval;

6

a warning occured while calculating the upper bound of the multinomial profile likelihood confidence interval.

CI.err

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.

CI.warn

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.

alpha

1-the confidence level of the interval.

N.CI

Not produced for normal heterogeneous models (h="Normal"): The x-coordinates for plot.closedpCI.t.

loglik.CI

Not produced for normal heterogeneous models (h="Normal"): The y-coordinates for plot.closedpCI.t.

Note

For normal heterogeneous models, the closedpCI functions use optim from the stats package. Otherwise, models are fitted with glm and the code to compute the multinomial profile likelihood confidence interval calls the functions glm, optimize and uniroot, all from the stats package.

Author(s)

Louis-Paul Rivest [e-mail] and Sophie Baillargeon

References

Baillargeon, S. and Rivest, L.P. (2007) Rcapture: Loglinear models for capture-recapture in R. Journal of Statistical Software, 19(5), doi: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.

See Also

closedp, closedp.Mtb

Examples

# 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")