peakfit: Finite mixture modelling of geochronological datasets

In IsoplotR: Statistical Toolbox for Radiometric Geochronology

Finite mixture modelling of geochronological datasets

Description

Implements the discrete mixture modelling algorithms of Galbraith and Laslett (1993) and applies them to fission track and other geochronological datasets.

Usage

peakfit(x, ...)

## Default S3 method:
peakfit(x, k = "auto", sigdig = 2, oerr = 3, log = TRUE, np = 3, ...)

## S3 method for class 'fissiontracks'
peakfit(x, k = 1, exterr = TRUE, sigdig = 2, log = TRUE, oerr = 3, np = 3, ...)

## S3 method for class 'UPb'
peakfit(
  x,
  k = 1,
  type = 4,
  cutoff.76 = 1100,
  cutoff.disc = discfilter(),
  common.Pb = 0,
  exterr = TRUE,
  sigdig = 2,
  log = TRUE,
  oerr = 3,
  np = 3,
  ...
)

## S3 method for class 'PbPb'
peakfit(
  x,
  k = 1,
  exterr = TRUE,
  sigdig = 2,
  log = TRUE,
  common.Pb = 0,
  oerr = 3,
  np = 3,
  ...
)

## S3 method for class 'ArAr'
peakfit(
  x,
  k = 1,
  exterr = TRUE,
  sigdig = 2,
  log = TRUE,
  i2i = FALSE,
  oerr = 3,
  np = 3,
  ...
)

## S3 method for class 'ThPb'
peakfit(
  x,
  k = 1,
  exterr = TRUE,
  sigdig = 2,
  log = TRUE,
  i2i = FALSE,
  oerr = 3,
  np = 3,
  ...
)

## S3 method for class 'KCa'
peakfit(
  x,
  k = 1,
  exterr = TRUE,
  sigdig = 2,
  log = TRUE,
  i2i = FALSE,
  oerr = 3,
  np = 3,
  ...
)

## S3 method for class 'ReOs'
peakfit(
  x,
  k = 1,
  exterr = TRUE,
  sigdig = 2,
  log = TRUE,
  i2i = TRUE,
  oerr = 3,
  np = 3,
  ...
)

## S3 method for class 'SmNd'
peakfit(
  x,
  k = 1,
  exterr = TRUE,
  sigdig = 2,
  log = TRUE,
  i2i = TRUE,
  oerr = 3,
  np = 3,
  ...
)

## S3 method for class 'RbSr'
peakfit(
  x,
  k = 1,
  exterr = TRUE,
  sigdig = 2,
  log = TRUE,
  i2i = TRUE,
  oerr = 3,
  np = 3,
  ...
)

## S3 method for class 'LuHf'
peakfit(
  x,
  k = 1,
  exterr = TRUE,
  sigdig = 2,
  log = TRUE,
  i2i = TRUE,
  oerr = 3,
  np = 3,
  ...
)

## S3 method for class 'ThU'
peakfit(
  x,
  k = 1,
  exterr = FALSE,
  sigdig = 2,
  log = TRUE,
  oerr = 3,
  Th0i = 0,
  np = 3,
  ...
)

## S3 method for class 'UThHe'
peakfit(x, k = 1, sigdig = 2, log = TRUE, oerr = 3, np = 3, ...)

Arguments

x	either an [nx2] matrix with measurements and their standard errors, or an object of class fissiontracks, UPb, PbPb, ThPb, ArAr, KCa, ReOs, SmNd, RbSr, LuHf, ThU or UThHe
...	optional arguments (not used)
k	the number of discrete age components to be sought. Setting this parameter to 'auto' automatically selects the optimal number of components (up to a maximum of 5) using the Bayes Information Criterion (BIC).
sigdig	number of significant digits to be used for any legend in which the peak fitting results are to be displayed.
oerr	indicates whether the analytical uncertainties of the output are reported in the plot legend as: 1: 1\sigma absolute uncertainties. 2: 2\sigma absolute uncertainties. 3: absolute (1-\alpha)% confidence intervals, where \alpha equales the value that is stored in settings('alpha'). 4: 1\sigma relative uncertainties (\%). 5: 2\sigma relative uncertainties (\%). 6: relative (1-\alpha)% confidence intervals, where \alpha equales the value that is stored in settings('alpha').
log	take the logs of the data before applying the mixture model?
np	number of parameters for the minimum age model. Must be either 3 or 4.
exterr	propagate the external sources of uncertainty into the component age errors?
type	scalar indicating whether to plot the ^{207}Pb/^{235}U age (type=1), the ^{206}Pb/^{238}U age (type=2), the ^{207}Pb/^{206}Pb age (type=3), the ^{207}Pb/^{206}Pb-^{206}Pb/^{238}U age (type=4), the concordia age (type=5), or the ^{208}Pb/^{232}Th age (type=6).
cutoff.76	the age (in Ma) below which the ^{206}Pb/^{238}U and above which the ^{207}Pb/^{206}Pb age is used. This parameter is only used if type=4.
cutoff.disc	discordance cutoff filter. This is an object of class discfilter.
common.Pb	common lead correction: 0:none 1: use the Pb-composition stored in settings('iratio','Pb206Pb204') (if x has class UPb and x$format<4); settings('iratio','Pb206Pb204') and settings('iratio','Pb207Pb204') (if x has class PbPb or x has class UPb and 3<x$format<7); or settings('iratio','Pb208Pb206') and settings('iratio','Pb208Pb207') (if x has class UPb and x$format=7 or 8). 2: use the isochron intercept as the initial Pb-composition 3: use the Stacey-Kramers two-stage model to infer the initial Pb-composition (only applicable if x has class UPb)
i2i	‘isochron to intercept’: calculates the initial (aka ‘inherited’, ‘excess’, or ‘common’) ^{40}Ar/^{36}Ar, ^{40}Ca/^{44}Ca, ^{207}Pb/^{204}Pb, ^{87}Sr/^{86}Sr, ^{143}Nd/^{144}Nd, ^{187}Os/^{188}Os, ^{230}Th/^{232}Th, ^{176}Hf/^{177}Hf or ^{204}Pb/^{208}Pb ratio from an isochron fit. Setting i2i to FALSE uses the default values stored in settings('iratio',...).
Th0i	initial ^{230}Th correction. 0: no correction 1: project the data along an isochron fit 2: if x$format is 1 or 2, correct the data using the measured present day ^{230}Th/^{238}U, ^{232}Th/^{238}U and ^{234}U/^{238}U activity ratios in the detritus. If x$format is 3 or 4, correct the data using the measured ^{238}U/^{232}Th activity ratio of the whole rock, as stored in x by the read.data() function. 3: correct the data using an assumed initial ^{230}Th/^{232}Th-ratio for the detritus (only relevant if x$format is 1 or 2).

Details

Consider a dataset of n dates \{t_1, t_2, ..., t_n\} with analytical uncertainties \{s[t_1], s[t_2], ..., s[t_n]\}. Define z_i = \log(t_i) and s[z_i] = s[t_i]/t_i. Suppose that these n values are derived from a mixture of k>2 populations with means \{\mu_1,...,\mu_k\}. Such a discrete mixture may be mathematically described by P(z_i|\mu,\omega) = \sum_{j=1}^k \pi_j N(z_i | \mu_j, s[z_j]^2 ) where \pi_j is the proportion of the population that belongs to the j^{th} component, and \pi_k=1-\sum_{j=1}^{k-1}\pi_j. This equation can be solved by the method of maximum likelihood (Galbraith and Laslett, 1993). IsoplotR implements the Bayes Information Criterion (BIC) as a means of automatically choosing k. This option should be used with caution, as the number of peaks steadily rises with sample size (n). If one is mainly interested in the youngest age component, then it is more productive to use an alternative parameterisation, in which all grains are assumed to come from one of two components, whereby the first component is a single discrete age peak (\exp(m), say) and the second component is a continuous distribution (as descibed by the central age model), but truncated at this discrete value. IsoplotR uses a normal likelihood function (section 6.11 of Galbraith, 2005) for the minimum age model. This may result in some inaccuracy for young and/or uranium-poor fission track samples.

Value

Returns a list with the following items:

peaks

a 2 x k matrix with the following rows:

t: the ages of the k peaks

s[t]: the standard errors of t

props

a 2 x k matrix with the following rows:

p: the proportions of the k peaks

s[p]: the standard errors of p

L

the log-likelihood of the fit

legend

a vector of text expressions to be used in a figure legend

References

Galbraith, R.F. and Laslett, G.M., 1993. Statistical models for mixed fission track ages. Nuclear Tracks and Radiation Measurements, 21(4), pp.459-470.

Galbraith, R.F. 2005, Statistics for fission track analysis. Chapman and Hall/CRC, 229p.

See Also

radialplot, central

Examples

attach(examples)
peakfit(FT1,k=2)

peakfit(LudwigMixture,k='min')

IsoplotR documentation built on Nov. 10, 2023, 9:08 a.m.