Luminescence (version 0.3.4)

calc_FiniteMixture: Apply the finite mixture model (FMM) after Galbraith (2005) to a given De distribution

Description

This function fits a k-component mixture to a De distribution with differing known standard errors. Parameters (doses and mixing proportions) are estimated by maximum likelihood assuming that the log dose estimates are from a mixture of normal distributions.

Usage

calc_FiniteMixture(input.data, sigmab, n.components, sample.id = "unknown sample", 
    n.iterations = 200, grain.probability = FALSE, main = "Finite Mixture Model", 
    dose.scale, pdf.weight = TRUE, pdf.sigma = "sigmab", pdf.colors = "gray", 
    pdf.scale, plot.proportions = TRUE)

Arguments

input.data

RLum.Results or data.frame (required): for data.frame: two columns with De (input.data[,1]) and De error (values[,2])

sigmab

numeric (required): spread in De values given as a fraction (e.g. 0.2). This value represents the expected overdispersion in the data should the sample be well-bleached (Cunningham & Wallinga 2012

n.components

numeric (required): number of components to be fitted. If a vector is provided (e.g. c(2:8)) the finite mixtures for 2, 3 ... 8 components are calculated and a plot and a statistica

sample.id

character (with default): sample id

n.iterations

numeric (with default): number of iterations for maximum log likelihood estimates

grain.probability

logical (with default): prints the estimated probabilities of which component each grain is in

main

character (with default): plot main title

dose.scale

numeric: manually set the scaling of the y-axis of the first plot with a vector in the form of c(min,max)

pdf.weight

logical (with default): weight the probability density functions by the components proportion (applies only when a vector is provided for n.components)

pdf.sigma

character (with default): if "sigmab" the components normal distributions are plotted with a common standard deviation (i.e. sigmab) as assumed by the FFM. Alternatively,

pdf.colors

character (with default): color coding of the components in the the plot. Possible options are "gray", "colors" and "none"

pdf.scale

numeric: manually set the max density value for proper scaling of the x-axis of the first plot

plot.proportions

logical (with default): plot barplot showing the proportions of components

Value

• Returns a terminal output. In addition a list is returned containing the following elements:
• mle.matrix
  matrix covariance matrix of maximum likelihood estimates.
• grain.probability
  matrix with estimated probabilities of which component each grain is in.
• meta
  data.frame containing model parameters (sample.id, sigmab, n.components, llik, bic).
• components
  data.frame containing fitted components.
• single.comp
  data.frame containing log likelihood and BIC for a single component.
• If a vector for n.components is provided (e.g. c(2:8)), mle.matrix, grain.probability and meta are lists containing matrices of the results for each iteration of the model. The output should be accessed using the function get_RLum.Results

Function version

0.32 (2014-05-09 17:14:59)

Details

This model uses the maximum likelihood and Bayesian Information Criterion (BIC) approaches. Indications of overfitting are: - increasing BIC - repeated dose estimates - covariance matrix not positive definite - covariance matrix produces NaNs - convergence problems Plot If a vector (c(k.min:k.max)) is provided for n.components a plot is generated showing the the k components equivalent doses as normal distributions. By default pdf.weight is set to FALSE, so that the area under each normal distribution is always 1. If TRUE, the probability density functions are weighted by the components proportion for each iteration of k components, so the sum of areas of each component equals 1. While the density values are on the same scale when no weights are used, the y-axis are individually scaled if the probability density are weighted by the components proportion. The standard deviation (sigma) of the normal distributions is by default determined by a common sigmab (see pdf.sigma). For pdf.sigma = "se" the standard error of each component is taken instead. The stacked barplot shows the proportion of each component (in per cent) calculated by the FFM. The last plot shows the achieved BIC scores and maximum log-likelihood estimates for each iteration of k.

References

Galbraith, R.F. & Green, P.F., 1990. Estimating the component ages in a finite mixture. Nuclear Tracks and Radiation Measurements, 17, pp. 197-206. Galbraith, R.F. & Laslett, G.M., 1993. Statistical models for mixed fission track ages. Nuclear Tracks Radiation Measurements, 4, pp. 459-470. Galbraith, R.F. & Roberts, R.G., 2012. Statistical aspects of equivalent dose and error calculation and display in OSL dating: An overview and some recommendations. Quaternary Geochronology, 11, pp. 1-27. Roberts, R.G., Galbraith, R.F., Yoshida, H., Laslett, G.M. & Olley, J.M., 2000. Distinguishing dose populations in sediment mixtures: a test of single-grain optical dating procedures using mixtures of laboratory-dosed quartz. Radiation Measurements, 32, pp. 459-465. Galbraith, R.F., 2005. Statistics for Fission Track Analysis, Chapman & Hall/CRC, Boca Raton. Further reading Arnold, L.J. & Roberts, R.G., 2009. Stochastic modelling of multi-grain equivalent dose (De) distributions: Implications for OSL dating of sediment mixtures. Quaternary Geochronology, 4, pp. 204-230. Cunningham, A.C. & Wallinga, J., 2012. Realizing the potential of fluvial archives using robust OSL chronologies. Quaternary Geochronology, 12, pp. 98-106. Rodnight, H., Duller, G.A.T., Wintle, A.G. & Tooth, S., 2006. Assessing the reproducibility and accuracy of optical dating of fluvial deposits. Quaternary Geochronology, 1, pp. 109-120. Rodnight, H. 2008. How many equivalent dose values are needed to obtain a reproducible distribution?. Ancient TL, 26, pp. 3-10.

See Also

calc_CentralDose, calc_CommonDose, calc_FuchsLang2001, calc_MinDose3, calc_MinDose4

Examples

## load example data
data(ExampleData.DeValues, envir = environment())

## (1) apply the finite mixture model
## NOTE: the data set is not suitable for the finite mixture model,
## which is why a very small sigmab is necessary
calc_FiniteMixture(ExampleData.DeValues,
                   sigmab = 0.08, n.components = 2,
                   grain.probability = TRUE)

## (2) repeat the finite mixture model for 2, 3 and 4 maximum number of fitted
## components and save results
## NOTE: The following example is computationally intensive. Please un-comment
## the following lines to make the example work.
#res<- calc_FiniteMixture(ExampleData.DeValues,
#                         sigmab = 0.01, n.components = c(2:4),
#                         pdf.weight = TRUE, dose.scale = c(2200,4500))

## show structure of the results
#res

## show the results on equivalent dose, standard error and proportion of
## fitted components
#get_RLum.Results(object=res, data.object="components")