basta | R Documentation |
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, ...) ## Default S3 method: 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, ...)
object |
A |
studyStart |
The first year of the study. |
studyEnd |
The last year of the study. |
minAge |
Age at which the analysis should start (see |
model |
The underlying mortality model to be used. |
shape |
The overall shape of the model. Values are: |
covarsStruct |
Character string that indicates how covariates should be evaluated. The options are: “ |
niter |
The total number of MCMC steps. |
burnin |
The number of iterations for the burn in (see |
thinning |
The number of skipped MCMC steps to minimize serial autocorrelation (see |
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, so the |
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 |
ncpus |
a numerical value that indicates the number of cpus to be used if |
lifeTable |
A logical argument indicating whether or not to produce life tables. If |
updateJumps |
A logical argument indicating wheter to update jump standard deviations (adaptive independent Metropolis) until an update rate of 0.25 is achieved (see |
... |
Additional arguments to be passed to function |
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 μ(x | θ), where x corresponds to age and θ 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
μ_b(x | b) = b,
where b > 0, (b) the Gompertz mortality function (“GO
”; Gompertz 1925, Pletcher 1999), calculated as
μ_b(x | b) = exp(b_0 + b_1 x),
where -∞ < b_0, b_1 < ∞, (c) the Weibull mortality model (“WE
”; Pinder III et al. 1978) calculated as
μ_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
μ_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 μ_0(x | θ)= μ_b(x | b) + c, where θ = [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
μ_0(x | θ) = exp(a_0 - a_1 x) + c + μ_b(x | b)
,
where -∞ < a_0 < ∞, a_1 ≤q 0 and c ≤q 0.
To incorporate covariates into the inference process, the mortality model is further extended by including a proportional hazards structure, of the form
μ(x | θ, Γ, Z_a, Z_c) = μ_0(x | θ, Z_a) exp(Γ Z_c)
,
where μ_0(x | θ, Z_a) represents the mortality section as defined above, while the second term exp(Γ 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 μ_0(x | b) = b, the parameter is equal to b = b_0 + b_1 z_1 + …, b_k z_k, where [b_0, …, b_k] are paramters that link the mortality parameter b with the categorical covariates [z_1,…,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.
coefficients |
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 |
DIC |
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 “ |
convergence |
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 “ |
KullbackLeibler |
If called by |
params |
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). |
settings |
If called by |
modelSpecs |
Model specifications inidicating the |
jumpPriors |
If requested or called by functions |
birthQuant |
If requested, summary matrix of estimated times of birth including median and upper and lower 95% predictive intervals. |
deathQuant |
If requested, summary matrix of estimated times of birth including median and upper and lower 95% predictive intervals. |
agesQuant |
If requested, summary matrix of estimated ages at death including median and upper and lower 95% predictive intervals. |
mortQuant |
If requested or called by functions |
survQuant |
If requested or called by functions |
parsForPlot |
If requested or called by functions |
lifeTable |
If requested and specified in the argument |
Fernando Colchero colchero@imada.sdu.dk, Owen R. Jones jones@biology.sdu.dk, and Maren Rebke maren.rebke@avitec-research.de
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.