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.
basta(object, ...)# S3 method for default
basta(object, studyStart, studyEnd, minAge = 0, model = "GO",
shape = "simple", covarsStruct = "fused", niter = 11000,
burnin = 1001, thinning = 20, recaptTrans = studyStart,
nsim = 1, parallel = FALSE, ncpus = 2, lifeTable = TRUE,
updateJumps = TRUE, ...)
A 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).
Basic deviance information criterion (DIC) calculations to be used for model selection (Spiegelhalter et al. 2002, Celeux et al. 2006). 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”.
If 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”.
If 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”.
If 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).
If 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.
Model specifications inidicating the model, the shape, the covariate structure and the minimum age that were specified by the user.
If 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.
If requested, summary matrix of estimated times of birth including median and upper and lower 95% predictive intervals.
If requested, summary matrix of estimated times of birth including median and upper and lower 95% predictive intervals.
If requested, summary matrix of estimated ages at death including median and upper and lower 95% predictive intervals.
If requested or called by functions plot or plot.basta median and 95% predictive intervals for the estimated mortality rates.
If requested or called by functions plot or plot.basta median and 95% predictive intervals for the estimated survival probability.
If requested or called by functions plot or plot.basta thinned sequences, including burn-in, of parameter traces used to plot outputs.
If requested and specified in the argument lifeTable, a cohort life table calculated from the estimated ages at death.
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, \dots, n_t-1\) represent the observation window of \(n_t\) years. This is followed (optionally) by columns for categorical and continuous covariates (see details).
The first year of the study.
The last year of the study.
Age at which the analysis should start (see details)
The underlying mortality model to be used. "EX" = exponential,"GO" = Gompertz, "WE" = Weibull and "LO" = logistic (see details).
The overall shape of the model. Values are: simple = no extra parameters added; Makeham = a constant parameter is added to the mortality; and bathtub = a Gompertz declining mortality for early ages and a constant parameter are added to the mortality model (see details).
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 hazards structure; “prop.haz”, which puts all covariates under a proportional hazards structure; and “all.in.mort” puts all covariates as a multilevel function of the mortality parameters (see details).
The total number of MCMC steps.
The number of iterations for the burn in (see details).
The number of skipped MCMC steps to minimize serial autocorrelation (see details).
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, so the recaptTrans is simply defined as a single element vector of the first year of the observation period (e.g. c(1985)). If recapture probabilities were known to change at year say, 1990, the RPTP should be defined as c(1985, 1990).
A numerical value for the number of simulations to be run.
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 is ran in series.
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 pkgsnowfall is not installed, the simulations are run in series.
A logical argument indicating whether or not to produce life tables. If TRUE, a cohort life table is calculated using function MakeLifeTable.
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)
Fernando Colchero colchero@imada.sdu.dk, Owen R. Jones jones@biology.sdu.dk, and Maren Rebke maren.rebke@avitec-research.de
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 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 functions, namely (a) Exponential (“EX”; Cox and Oakes 1974), with constant mortality with age specified as
$$\mu_b(x | b) = b,$$
where \(b > 0\), (b) the Gompertz mortality function (“GO”; Gompertz 1925, Pletcher 1999), calculated as
$$\mu_b(x | b) = exp(b_0 + b_1 x),$$
where \(-\infty < b_0, b_1 < \infty\), (c) the Weibull mortality model (“WE”; Pinder III et al. 1978) calculated as
$$\mu_b(x | b) = b_0 b_1^{b_0} x^{b_0 -1},$$
where \(b_0, b_1 > 0\), and (d) the logistic mortality model (“LO”; Pletcher 1999), calculated as
$$\mu_b(x | b) = exp(b_0 + b_1 x) / (1 + b_2 exp(b_0)/b_1 (exp(b_1 x)-1)),$$ where \(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, 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)$$,
where \(-\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 model with constant mortality of the form \(\mu_0(x | b) = b\), the parameter is equal to \(b = b_0 + b_1 z_1 + \dots, b_k z_k\), 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”.
Another additional argument is returnAges, which outputs a value estAges, which is a matrix with all estimated ages after the burnin and thinning. This matrix also includes ages that were known, thus each row corresponds to the individual in the same row as the original dataset.
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.
Celeux, G., Forbes, F., Robert, C. P., and Titterington, D. M. (2006) Deviance information criteria for missing data models. Bayesian Analysis, 1(4), 651-673.
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.
summary.basta, print.basta, plot.basta to visualise summary outputs for objects of class “basta”.
CensusToCaptHist and MakeCovMat for raw data formatting.
## 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 = 100, 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