basta: Parametric Bayesian estimation of age-specific survival for left-truncated and right-censored capture-recapture/recovery data.

Description

This function performs multiple Markov Chain Monte Carlo (MCMC) simulations for the Bayesian estimation of age-specific mortality and survival trends when a large proportion of (or all) records have unknown times of birth and/or death. Survival parameters and unknown (i.e. latent) birth and death times are estimated, allowing the user to test a range of mortality patterns, and to test the effect of continuous and/or discrete covariates following Colchero and Clark's (2012) approach.

Usage

basta(object, ...)
## S3 method for class 'Additional arguments to be passed to function \code{basta':
\item{\dots (see \code{details}), covarsStruct = "fused", niter = 11000, 
    burnin = 1001, thinning = 20, recaptTrans = studyStart, 
    nsim = 1, parallel = FALSE, ncpus = 2, lifeTable = TRUE, 
    updateJumps = TRUE, ...)

Arguments

object

A data.frame to be used as an input data file for BaSTA. The first column is a vector of individual unique IDs, the second and third columns are birth and death years respectively. Columns 4-(nt-1) represent the observation window of nt years

studyStart

The first year of the study.

studyEnd

The last year of the study.

minAge

Age at which the analysis should start (see details)

model

The underlying mortality model to be used. "EX" = exponential, "GO"= Gompertz, "WE" = Weibull and "LO" = logistic (see details).

shape

The overall shape of the model. Values are: simple = no extra parameters added; Makeham = a constant parameter is added to the mortality rate; and bathtub = a Gompertz declining mortality for early ages and a constan

covarsStruct

Character string that indicates how covariates should be evaluated. The options are: fused, which defines all categorical variables as covariates for each mortality parameter and all continuous covariates under a proportional

niter

The total number of MCMC steps.

burnin

The number of iterations for the burn in (see details).

thinning

The number of skipped MCMC steps to minimize serial autocorrelation (see details).

recaptTrans

A vector (of length npi) defining the recapture probability transition times (RPTP). These are points (years) where the recapture probability is thought to change. The default setting is for the recapture probability to be constant throughout the study,

nsim

A numerical value for the number of simulations to be run.

parallel

A logical argument indicating whether the multiple simulations should be run in parallel or not. If TRUE, package snowfall is called and multiple simulations are run in parallel. If snowfall is not installed, the model i

ncpus

a numerical value that indicates the number of cpus to be used if parallel is TRUE and package snowfall is installed. The default is 2 cpus. If package pkg{snowfall} is not installed, the simulations are run in series.

lifeTable

A logical argument indicating whether or not to produce life tables. If TRUE, a cohort life table is calculated using function MakeLifeTable.

updateJumps

A logical argument indicating wheter to update jump standard deviations (adaptive independent Metropolis) until an update rate of 0.25 is achieved (see details).

...

Additional arguments to be passed to function basta (see details)

Value

coefficientsA matrix with estimated coefficients (i.e. mean values per parameter on the thinned sequences after burnin), which includes standard errors, upper and lower 95% credible intervals, update rates per parameter, serial autocorrelation on the thinned sequences and the potential scale reduction factor for convergence (see Convergence value below).
DICBasic deviance information criterion (DIC) calculations to be used for model selection (Spiegelhalter et al. 2002). Small differences between values should only be used a reference (see comments in Spiegelhalter et al. 2002). If all or some of the simulations failed, then the returned value is Not calculated.
convergenceIf requested, a matrix with convergence coefficients based on potential scale reduction as described by Gelman et al. (2004). If only one simulation was ran, then the returned value is Not calculated.
KullbackLeiblerIf called by summary, list with Kullback-Leibler discrepancy matrices between pair of parameters for categorical covariates (McCulloch 1989, Burnham and Anderson 2001) and McCulloch's (1989) calibration measure. If only one simulation was ran or if no convergence is reached, then the returned value is Not calculated.
paramsIf requested, a matrix with the thinned, converged parameter traces of all runs. This matrix is used to calculate quantiles for parameters, survival probability and mortality (see below).
settingsIf called by summary, this is a vector indicating the number of iterations for each MCMC, the burn in sequence, the thinning interval, and the number of simulations that were run.
modelSpecsModel specifications inidicating the model, the shape, the covariate structure and the minimum age that were specified by the user.
jumpPriorsIf requested or called by functions summary or summary.basta, a matrix with the jump standard deviations and prior means and standard deviation for the parameters.
birthQuantIf requested, summary matrix of estimated times of birth including median and upper and lower 95% predictive intervals.
deathQuantIf requested, summary matrix of estimated times of birth including median and upper and lower 95% predictive intervals.
agesQuantIf requested, summary matrix of estimated ages at death including median and upper and lower 95% predictive intervals.
mortQuantIf requested or called by functions plot or plot.basta median and 95% predictive intervals for the estimated mortality rates.
survQuantIf requested or called by functions plot or plot.basta median and 95% predictive intervals for the estimated survival probability.
parsForPlotIf requested or called by functions plot or plot.basta thinned sequences, including burn-in, of parameter traces used to plot outputs.
lifeTableIf requested and specified in the argument lifeTable, a cohort life table calculated from the estimated ages at death.

Details

To construct the input data object the function CensusToCaptHist can be used to build the capture-recapture matrix, while the covariate (design) matrix can be constructed with the MakeCovMat function.

basta uses parametric mortality functions to estimate age-specific mortality (survival) from capture-recapture/recovery data. The mortality rate function describes how the risk of mortality changes with age, and is defined as $\mu(x | \theta)$, where $x$ corresponds to age and $\theta$ is a vector of parameters to be estimated.

The model argument allows the user to choose between four basic mortality rate functions:

(a) Exponential (EX; Cox and Oakes 1974), with constant mortality rate:

$$\mu_b(x | b) = b$$ with $b > 0$.

(b) Gompertz (GO; Gompertz 1925, Pletcher 1999):

$$\mu_b(x | b) = exp(b_0 + b_1 x)$$ with $-\infty < b_0, b_1 < \infty$.

$$\mu_b(x | b) = b_0 b_1^(b_0) x^(b_0 -1)$$ with $b_0, b_1 > 0$.

(d) logistic (LO; Pletcher 1999):

$$\mu_b(x | b) = exp(b_0 + b_1 x) / (1 + b_2 exp(b_0)/b_1 (exp(b_1 x)-1))$$ with $b_0, b_1, b_2 > 0$.

The shape argument allows the user to extend these models in order to explore more complex mortality shapes. The default value is simple which leaves the model as defined above. With value Makeham, a constant is added to the mortality rate, making the model equal to $\mu_0(x | \theta)= \mu_b(x | b) + c$, where $\theta = [c, b]$. With value bathtub, concave shapes in mortality can be explored. This is achieved by adding a declining Gompertz term and a constant parameter to the basic mortality model, namely:

$$\mu_0(x | \theta) = exp(a_0 - a_1 x) + c + \mu_b(x | b)$$. with $-\infty < a_0 < \infty$, $a_1 \leq 0$ and $c \leq 0$.

To incorporate covariates into the inference process, the mortality model is further extended by including a proportional hazards structure, of the form:

$$\mu(x | \theta, \Gamma, Z_a, Z_c) = \mu_0(x | \theta, Z_a) exp(\Gamma Z_c)$$

where $\mu_0(x | \theta, Z_a)$ represents the mortality section as defined above, while the second term $exp(\Gamma Z_c)$ corresponds to the proportional hazards function. $Z_a$ and $Z_c$ are covariate (design) matrices for categorical and continuous covariates, respectively.

When covariates are included in the dataset, the basta function provides three different ways in which these can be evaluated by using argument covarsStruct:

1. fused will make the mortality parameters linear functions of all categorical covariates (analogous to a generalised linear model (GLM) structure) and will put all continuous covariates under a proportional hazards structure. Thus, for a simple exponential mortality rate of the form $\mu_0(x | b) = b$, the parameter is equal to $b = b_0 + b_1 z_1 + \dots$, where $[b_0, \dots, b_k]$ are paramters that link the mortality parameter $b$ with the categorical covariates $[z_1,\dots,z_k]$.

2. prop.haz will put all covariates under a proportional hazards structure irrespective of the type of variable.

3. all.in.mort will put all covariates as linear functions of the survival parameters as explained above. Since most models require the lower bounds for the mortality parameters to be equal to 0, the only model that can be used for this test is Gompertz with shape set to simple. In case these arguments are specified deferently, a warning message is printed noting that model will be forced to be GO and shape will be set to simple.

The burnin argument represents the number of steps at the begining of the MCMC run that is be discarded. This sequence commonly corresponds to the non-converged section of the MCMC sequence. Convergence and model selection measures are calculated from the remaining thinned parameter chains if multiple simulations are run, and all if all of them run to completion.

The thinning argument specifies the number of steps to be skipped in order to reduce serial autocorrelation. The thinned sequence, which only includes steps after burn in, is then used to calculate convergence statistics and model for selection.

The updateJumps argument specifies wether to run a simulation to find appropriate jump standard deviations for theta and gamma parameters. If argument nsim is set to 1, then the simulation runs with the update jumps routine active. If nsim is larger than 1, then an initial simulation is ran to find apropriate jumps before the main analysis is ran.

Additional arguments for priors, jumps and start values can be passed on the ...section. For instance, argument thetaStart can be specified as a vector defining the initial values for each parameter in the survival model. If this argument is not specified, a set of random parameters is generated for each simulation. Similarly, argument gammaStart can be specified for all parameters in the proportional hazards section of the model. Jump standard deviations (i.e. the standard error in the Metropolis step) can be specified with arguments thetaJumps and gammaJumps. As with thetaStart, default values are assigned if these arguments are not specified. Arguments thetaPriorMean, thetaPriorSd, gammaPriorMean and gammaPriorSd can be used to specify prior means and standard errors for each survival and proportional hazards parameters. Again, if not specified, default values are assigned.

The number of parameters in thetaStart, thetaJumps, thetaPriorMean and thetaPriorSd should be a vector or matrix for the parameters in the mortality function. The number of parameters will depend on the model chosen with model (see above). If the number of parameters specified does not match the number of parameters inherent to the model and shape selected, the function returns an error.

As described above, the number of parameters for gammaStart, gammaJumps, gammaPriorMean and gammaPriorSd arguments (i.e. section b), namely the proportional hazards section, will be a function of the number of continuous covariates if argument covarsStruct is fused, or to the total number of covariates when covarsStruct is prop.haz.

References

Burnham, K.P. and Anderson, D.R. (2001) Kullback-Leibler information as a basis for strong inference in ecological studies. Widlife Research, 28, 111-119.

Colchero, F. and J.S. Clark (2012) Bayesian inference on age-specific survival from capture-recapture data for censored and truncated data. Journal of Animal Ecology. 81, 139-149.

Colchero, F., O.R. Jones and M. Rebke. (2012) BaSTA: an R package for Bayesian estimation of age-specific survival from incomplete mark-recapture/recovery data with covariates. Method in Ecology and Evolution. 3, 466-470.

Cox, D. R., and Oakes D. (1984) Analysis of Survival Data. Chapman and Hall, London.

Gelman, A., Carlin, J.B., Stern, H.S. and Rubin, D.B. (2004) Bayesian data analysis. 2nd edn. Chapman & Hall/CRC, Boca Raton, Florida, USA.

Gompertz, B. (1825) On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies. Philosophical Transactions of the Royal Society of London, 115, 513-583.

King, R. and Brooks, S.P. (2002) Bayesian model discrimination for multiple strata capture-recapture data. Biometrika, 89, 785-806.

McCulloch, R.E. (1989) Local model influence. Journal of the American Statistical Association, 84, 473-478.

Pinder III, J.E., Wiener, J.G. and Smith, M.H. (1978) The Weibull distribution: a new method of summarizing survivorship data. Ecology, 59, 175-179.

Spiegelhalter, D.J., Best, N.G., Carlin, B.P. and van der Linde, A. (2002) Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society: Series B, 64, 583-639.

Examples

Run this code

## Load data:
data("sim1", package = "BaSTA")

## Check data consistency:
new.dat  <- DataCheck(sim1, studyStart = 51, 
            studyEnd = 70, autofix = rep(1,7))

## Run short version of BaSTA on the data:
out <- basta(sim1, studyStart = 51, studyEnd = 70, 
             niter = 200, burnin = 11, thinning = 10,
             updateJumps = FALSE)

## Print results:
summary(out, digits = 3)

## Plot traces for survival parameters:
plot(out)

## Plot traces for proportional hazards parameter:
plot(out, trace.name = "gamma")

## Plot survival and mortality curves:
plot(out, plot.trace = FALSE)

Run the code above in your browser using DataLab