Bchron (version 4.7.3)

Bchronology: Runs the Compound Poisson-Gamma chronology model of Haslett and Parnell (2008)

Description

Fits a non-parametric chronology model to age/position data according to the Compound Poisson-Gamma model defined by Haslett and Parnell (2008) <DOI:10.1111/j.1467-9876.2008.00623.x>. This version uses a slightly modified Markov chain Monte Carlo fitting algorithm which aims to converge quicker and requires fewer iterations. It also a slightly modified procedure for identifying outliers

Usage

Bchronology(
  ages,
  ageSds,
  positions,
  positionThicknesses = rep(0, length(ages)),
  calCurves = rep("intcal20", length(ages)),
  ids = NULL,
  outlierProbs = rep(0.01, length(ages)),
  predictPositions = seq(min(positions), max(positions), length = 100),
  pathToCalCurves = system.file("data", package = "Bchron"),
  jitterPositions = FALSE,
  iterations = 10000,
  burn = 2000,
  thin = 8,
  extractDate = 1950 - as.numeric(format(Sys.time(), "%Y")),
  maxExtrap = 1000,
  thetaMhSd = 0.5,
  muMhSd = 0.1,
  psiMhSd = 0.1,
  ageScaleVal = 1000,
  positionNormalise = TRUE
)

Arguments

ages

A vector of ages (most likely 14C)

ageSds

A vector of 1-sigma values for the ages given above

positions

Position values (e.g. depths) for each age

positionThicknesses

Thickness values for each of the positions. The thickness value should be the full thickness value of the slice. By default set to zero.

calCurves

A vector of values containing either 'intcal20', 'shcal20', 'marine20', or 'normal' (older calibration curves are also supported, e.g. intcal13). Should be the same length the number of ages supplied. Non-standard calibration curves can be used provided they are supplied in the same format as those previously mentioned and are placed in the same directory, or created via createCalCurve. Normal indicates a normally-distributed (non-14C) age.

ids

ID names for each age

outlierProbs

A vector of prior outlier probabilities, one for each age. Defaults to 0.01

predictPositions

A vector of positions (e.g. depths) at which predicted age values are required. Defaults to a sequence of length 100 from the top position to the bottom position

pathToCalCurves

File path to where the calibration curves are located. Defaults to the system directory where the 3 standard calibration curves are stored.

jitterPositions

Whether to jigger the positions at startup or not. Default is FALSE but if there are lots of dates at similar depths this may resolve some initialisation problems

iterations

The number of iterations to run the procedure for

burn

The number of starting iterations to discard

thin

The step size for every iteration to keep beyond the burnin

extractDate

The top age of the core. Used for extrapolation purposes so that no extrapolated ages go beyond the top age of the core. Defaults to the current year

maxExtrap

The maximum number of extrapolations to perform before giving up and setting the predicted ages to NA. Useful for when large amounts of extrapolation are required, i.e. some of the predictPositions are a long way from the dated positions

thetaMhSd

The Metropolis-Hastings standard deviation for the age parameters

muMhSd

The Metropolis-Hastings standard deviation for the Compound Poisson-Gamma mean

psiMhSd

The Metropolis-Hastings standard deviation for the Compound Poisson-Gamma scale

ageScaleVal

A scale value for the ages. Bchronology works best when the ages are scaled to be approximately between 0 and 100. The default value is thus 1000 for ages given in years.

positionNormalise

Whether to normalise the position values. Bchronology works best when the positions are normalise to be between 0 and 1 The default value is TRUE

Value

A list of class BchronologyRun which include elements:

theta

The posterior estimated values of the ages

phi

The posterior estimated outlier values (1=outlier, 2=not outlier). The means of this parameter give the posterior estimated outlier probabilities

mu

The posterior values of the Compound Poisson-Gamma mean

psi

The posterior values of the Compound Poisson-Gamma scale

thetaPredict

The posterior estimated ages for each of the values in predictPosition

predictPositions

The positions at which estimated ages were required

calAges

The calibrated ages as output from BchronCalibrate

inputVals

All of the input values to the Bchronology run

Details

The Bchronology function fits a compound Poisson-Gamma distribution to the increments between the dated levels. This involves a stochastic linear interpolation step where the age gaps are Gamma distributed, and the position gaps are Exponential. Radiocarbon and non-radiocarbon dates (including outliers) are updated within the function also by MCMC.

References

Haslett, J., and Parnell, A. C. (2008). A simple monotone process with application to radiocarbon-dated depth chronologies. Journal of the Royal Statistical Society, Series C, 57, 399-418. DOI:10.1111/j.1467-9876.2008.00623.x Parnell, A. C., Haslett, J., Allen, J. R. M., Buck, C. E., and Huntley, B. (2008). A flexible approach to assessing synchroneity of past events using Bayesian reconstructions of sedimentation history. Quaternary Science Reviews, 27(19-20), 1872-1885. DOI:10.1016/j.quascirev.2008.07.009

See Also

BchronCalibrate, BchronRSL, BchronDensity, BchronDensityFast

Examples

# NOT RUN {
# Data from Glendalough
data(Glendalough)

# Run in Bchronology - all but first age uses intcal20
GlenOut = Bchronology(ages=Glendalough$ages,ageSds=Glendalough$ageSds,
                      calCurves=Glendalough$calCurves,positions=Glendalough$position,
                      positionThicknesses=Glendalough$thickness,ids=Glendalough$id,
                      predictPositions=seq(0,1500,by=10))

# Summarise it a few different ways
summary(GlenOut) # Default is for quantiles of ages at predictPosition values
summary(GlenOut, type='convergence') # Check model convergence
summary(GlenOut, type='outliers') # Look at outlier probabilities

# Predict for some new positions
predictAges = predict(GlenOut, newPositions = c(150,725,1500), newPositionThicknesses=c(5,0,20))

# Plot the output
plot(GlenOut,main="Glendalough",xlab='Age (cal years BP)',ylab='Depth (cm)',las=1)
# }
# NOT RUN {
# }