Fit a hamstr Age-Depth Model — hamstr • hamstr

hamstr is used to fit an age-depth model to a set of age-control points. Ages should already be on the desired scale, e.g. calendar ages, and will not be calibrated. The function calibrate_14C_age can be used to calibrate radiocarbon dates prior to fitting a hamstr model.

Usage

hamstr(
  depth,
  obs_age,
  obs_err,
  min_age = 1950 - as.numeric(format(Sys.Date(), "%Y")),
  K_fine = NULL,
  K_factor = NULL,
  K,
  top_depth = NULL,
  bottom_depth = NULL,
  acc_mean_prior = NULL,
  acc_shape = 4,
  mem_mean = 0.5,
  mem_strength = 10,
  model_bioturbation = FALSE,
  n_ind = NULL,
  L_prior_mean = 7.5,
  L_prior_shape = 2.5,
  L_prior_sigma = NULL,
  model_displacement = FALSE,
  D_prior_mean = 2,
  D_prior_shape = 1.5,
  D_prior_sigma = NULL,
  model_hiatus = FALSE,
  H_top = NULL,
  H_bottom = NULL,
  H_max = NULL,
  sample_posterior = TRUE,
  hamstr_control = list(),
  stan_sampler_args = list()
)

Arguments

depth: depths of observed ages (age control points)
obs_age: observed age at each depth (age control points)
obs_err: error associated with each observed age (1 standard error)
min_age: the minimum age that the first modelled depth can be. Useful if extrapolating above the shallowest age control point to e.g. the surface. So set min_age to the year the core was collected. E.g. for a core collected in 1990, with ages in years BP this would be -40 (present = 1950 by convention). The default value is the current year in calendar age BP (e.g. -71 for 2021).
K_fine: the number of sections at the highest resolution of the model.
K_factor: the rate at which the thickness of the sections grows between subsequent levels.
K: deprecated, use K_fine and K_factor instead.
top_depth, bottom_depth: the top and bottom depths of the desired age-depth model. Must encompass the range of the data. Defaults to the shallowest and deepest data points.
acc_mean_prior: hyperparameter for the prior on the overall mean accumulation rate for the record. Units are obs_age / depth. E.g. if depth is in cm and age in years then the accumulation rate is in years/cm. The overall mean accumulation rate is given a weak half-normal prior with mean = 0, SD = 10 * acc_mean_prior. If left blank, acc_mean_prior is set to the mean accumulation rate estimated by fitting a robust linear model using rlm.
acc_shape: hyperparameter for the shape of the priors on accumulation rates. Defaults to 4.
mem_mean: hyperparameter; a parameter of the Beta prior distribution on "memory", i.e. the autocorrelation parameter in the underlying AR1 model. The prior on the correlation between layers is scaled according to the thickness of the sediment sections in the highest resolution hierarchical layer, delta_c, which is determined by the total length age-models and the parameter vector K. mem_mean sets the mean value for R (defaults to 0.5), while w = R^(delta_c)
mem_strength: hyperparameter: sets the strength of the memory prior, defaults to 10 as in Bacon >= 2.5.1
model_bioturbation: defaults to FALSE. If TRUE, additional uncertainty in the observed ages due to sediment mixing (bioturbation) is modelled via a latent variable process. The amount of additional uncertainty is a function of the mixing depth L, the sedimentation rate, and the number of particles (e.g. individual foraminifera) per measured date. See description for details.
n_ind: the number of individual particles (e.g. foraminifera) in each sample that was dated by e.g. radiocarbon dating. This can be a single value or a vector the same length as obs_age. Only used if model_bioturbation = TRUE.
L_prior_mean: mean of the gamma prior on mixing depth, defaults to 7.5.
L_prior_shape, L_prior_sigma: shape and standard deviation of the gamma prior on the mixing depth. Set only one of these, the other will be calculated. Defaults to shape = 2.5. If either the shape or sigma parameter is set to zero, the mixing depth is fixed at the value of L_prior_mean, rather than being sampled.
model_displacement: model additional error on observed ages that does not scale with the number of individual particles in a sample, for example due to incomplete mixing.
D_prior_mean: mean of the gamma prior on additional error on observed ages. Units are those of the depth variable, e.g. cm.
D_prior_shape, D_prior_sigma: shape and standard deviation of the gamma prior on the additional error on observed ages. Set only one of these, the other will be calculated. Defaults to shape = 1.5. If either the shape or sigma parameter is set to zero, the additional error if fixed at the value of D_prior_mean, rather than being sampled.
model_hiatus: optionally model a hiatus.
H_top, H_bottom: limits to the location of a hiatus. By default these are set to the top and bottom data points but can be set by the user
H_max: maximum length of the hiatus in age units
sample_posterior: if set to FALSE, hamstr skips sampling the model and returns only the data, model structure and prior parameters so that data and prior distributions can be plotted and checked prior to running a model. Defaults to TRUE
hamstr_control: additional arguments to hamstr useful for debugging or development. See hamstr_control for details.
stan_sampler_args: additional arguments to sampling passed as a named list. e.g. list(chains = 8, iter = 4000) to run 8 MCMC chains of 4000 iterations instead of the default 4 chains of 2000 iterations. See get_stan_sampler_args for details.

Value

Returns an object of class "hamstr_fit", which is a list composed of the output from the stan sampler .$fit, and the list of data passed to the sampler, .$data

Examples

if (FALSE) { # \dontrun{

fit <- hamstr(
  depth = MSB2K$depth,
  obs_age = MSB2K$age,
  obs_err = MSB2K$error)

plot(fit)

} # }