Averages heteroscedastic data either using the ordinary weighted mean, or using a random effects model with two sources of variance. Computes the MSWD of a normal fit without overdispersion. Implements a modified Chauvenet criterion to detect and reject outliers. Only propagates the systematic uncertainty associated with decay constants and calibration factors after computing the weighted mean isotopic composition. Does not propagate the uncertainty of any initial daughter correction, because this is neither a purely random or purely systematic uncertainty.
weightedmean(x, ...)# S3 method for default
weightedmean(
x,
from = NA,
to = NA,
random.effects = FALSE,
detect.outliers = TRUE,
plot = TRUE,
levels = NA,
clabel = "",
rect.col = c("#00FF0080", "#FF000080"),
outlier.col = "#00FFFF80",
sigdig = 2,
oerr = 3,
ranked = FALSE,
hide = NULL,
omit = NULL,
omit.col = NA,
...
)
# S3 method for other
weightedmean(
x,
from = NA,
to = NA,
random.effects = FALSE,
detect.outliers = TRUE,
plot = TRUE,
levels = NA,
clabel = "",
rect.col = c("#00FF0080", "#FF000080"),
outlier.col = "#00FFFF80",
sigdig = 2,
oerr = 3,
ranked = FALSE,
hide = NULL,
omit = NULL,
omit.col = NA,
...
)
# S3 method for UPb
weightedmean(
x,
random.effects = FALSE,
detect.outliers = TRUE,
plot = TRUE,
from = NA,
to = NA,
levels = NA,
clabel = "",
rect.col = c("#00FF0080", "#FF000080"),
outlier.col = "#00FFFF80",
sigdig = 2,
type = 4,
cutoff.76 = 1100,
oerr = 3,
cutoff.disc = discfilter(),
exterr = FALSE,
ranked = FALSE,
common.Pb = 0,
hide = NULL,
omit = NULL,
omit.col = NA,
...
)
# S3 method for PbPb
weightedmean(
x,
random.effects = FALSE,
detect.outliers = TRUE,
plot = TRUE,
from = NA,
to = NA,
levels = NA,
clabel = "",
rect.col = c("#00FF0080", "#FF000080"),
outlier.col = "#00FFFF80",
sigdig = 2,
oerr = 3,
exterr = FALSE,
common.Pb = 2,
ranked = FALSE,
hide = NULL,
omit = NULL,
omit.col = NA,
...
)
# S3 method for ThU
weightedmean(
x,
random.effects = FALSE,
detect.outliers = TRUE,
plot = TRUE,
from = NA,
to = NA,
levels = NA,
clabel = "",
rect.col = c("#00FF0080", "#FF000080"),
outlier.col = "#00FFFF80",
sigdig = 2,
oerr = 3,
ranked = FALSE,
Th0i = 0,
hide = NULL,
omit = NULL,
omit.col = NA,
...
)
# S3 method for ArAr
weightedmean(
x,
random.effects = FALSE,
detect.outliers = TRUE,
plot = TRUE,
from = NA,
to = NA,
levels = NA,
clabel = "",
rect.col = c("#00FF0080", "#FF000080"),
outlier.col = "#00FFFF80",
sigdig = 2,
oerr = 3,
exterr = FALSE,
ranked = FALSE,
i2i = FALSE,
hide = NULL,
omit = NULL,
omit.col = NA,
...
)
# S3 method for KCa
weightedmean(
x,
random.effects = FALSE,
detect.outliers = TRUE,
plot = TRUE,
from = NA,
to = NA,
levels = NA,
clabel = "",
rect.col = c("#00FF0080", "#FF000080"),
outlier.col = "#00FFFF80",
sigdig = 2,
oerr = 3,
exterr = FALSE,
ranked = FALSE,
i2i = FALSE,
hide = NULL,
omit = NULL,
omit.col = NA,
...
)
# S3 method for ThPb
weightedmean(
x,
random.effects = FALSE,
detect.outliers = TRUE,
plot = TRUE,
from = NA,
to = NA,
levels = NA,
clabel = "",
rect.col = c("#00FF0080", "#FF000080"),
outlier.col = "#00FFFF80",
sigdig = 2,
oerr = 3,
exterr = FALSE,
ranked = FALSE,
i2i = TRUE,
hide = NULL,
omit = NULL,
omit.col = NA,
...
)
# S3 method for ReOs
weightedmean(
x,
random.effects = FALSE,
detect.outliers = TRUE,
plot = TRUE,
from = NA,
to = NA,
levels = NA,
clabel = "",
rect.col = c("#00FF0080", "#FF000080"),
outlier.col = "#00FFFF80",
sigdig = 2,
oerr = 3,
exterr = FALSE,
ranked = FALSE,
i2i = TRUE,
hide = NULL,
omit = NULL,
omit.col = NA,
...
)
# S3 method for SmNd
weightedmean(
x,
random.effects = FALSE,
detect.outliers = TRUE,
plot = TRUE,
from = NA,
to = NA,
levels = NA,
clabel = "",
rect.col = c("#00FF0080", "#FF000080"),
outlier.col = "#00FFFF80",
sigdig = 2,
oerr = 3,
exterr = FALSE,
ranked = FALSE,
i2i = TRUE,
hide = NULL,
omit = NULL,
omit.col = NA,
...
)
# S3 method for RbSr
weightedmean(
x,
random.effects = FALSE,
detect.outliers = TRUE,
plot = TRUE,
from = NA,
to = NA,
levels = NA,
clabel = "",
rect.col = c("#00FF0080", "#FF000080"),
outlier.col = "#00FFFF80",
sigdig = 2,
oerr = 3,
exterr = FALSE,
i2i = TRUE,
ranked = FALSE,
hide = NULL,
omit = NULL,
omit.col = NA,
...
)
# S3 method for LuHf
weightedmean(
x,
random.effects = FALSE,
detect.outliers = TRUE,
plot = TRUE,
from = NA,
to = NA,
levels = NA,
clabel = "",
rect.col = c("#00FF0080", "#FF000080"),
outlier.col = "#00FFFF80",
sigdig = 2,
oerr = 3,
exterr = FALSE,
i2i = TRUE,
ranked = FALSE,
hide = NULL,
omit = NULL,
omit.col = NA,
...
)
# S3 method for UThHe
weightedmean(
x,
random.effects = FALSE,
detect.outliers = TRUE,
plot = TRUE,
from = NA,
to = NA,
levels = NA,
clabel = "",
rect.col = c("#00FF0080", "#FF000080"),
outlier.col = "#00FFFF80",
sigdig = 2,
oerr = 3,
ranked = FALSE,
hide = NULL,
omit = NULL,
omit.col = NA,
...
)
# S3 method for fissiontracks
weightedmean(
x,
random.effects = FALSE,
detect.outliers = TRUE,
plot = TRUE,
from = NA,
to = NA,
levels = NA,
clabel = "",
rect.col = c("#00FF0080", "#FF000080"),
outlier.col = "#00FFFF80",
sigdig = 2,
oerr = 3,
exterr = FALSE,
ranked = FALSE,
hide = NULL,
omit = NULL,
omit.col = NA,
...
)
Returns a list with the following items:
a two or three element vector with:
t: the weighted mean. An asterisk is added to the plot title
if the initial daughter correction is based on an isochron
regression, to mark the circularity of using an isochron to compute
a weighted mean.
s[t]: the standard error of the weighted mean, excluding the
uncertainty of the initial daughter correction. This is because
this uncertainty is neither purely random nor purely systematic.
a two-element vector with the (over)dispersion and its standard error.
the Mean Square of the Weighted Deviates (a.k.a. `reduced Chi-square' statistic)
the number of degrees of freedom of the Chi-square test for homogeneity (\(df=n-1\), where \(n\) is the number of samples).
the p-value of a Chi-square test with \(df\) degrees of freedom, testing the null hypothesis that the underlying population is not overdispersed.
vector of logical flags indicating which steps are included into the weighted mean calculation
list of plot parameters for the weighted mean
diagram, including mean (the mean value), ci (a grey
rectangle with the (1 s.e., 2 s.e. or 100[1-\(\alpha\)]%,
depending on the value of oerr) confidence interval ignoring
systematic errors), ci.exterr (a grey rectangle with the
confidence interval including systematic errors), dash1 and
dash2 (lines marking the confidence interval augmented by
\(\sqrt{mswd}\) overdispersion if random.effects=FALSE),
and marking the confidence limits of a normal distribution whose
standard deviation equals the overdispersion parameter if
random.effects=TRUE).
a two column matrix of values (first column) and their
standard errors (second column) OR an object of class
UPb, PbPb, ThPb, ArAr, KCa,
ReOs, SmNd, RbSr, LuHf, ThU,
fissiontracks or UThHe
optional arguments
minimum y-axis limit. Setting from=NA scales the
plot automatically.
maximum y-axis limit. Setting to=NA scales the
plot automatically.
if TRUE, computes the weighted mean
using a random effects model with two parameters: the mean and
the dispersion. This is akin to a `model-3' isochron
regression.
if FALSE, attributes any excess dispersion to an
underestimation of the analytical uncertainties. This akin to a
`model-1' isochron regression.
logical flag indicating whether outliers should be detected and rejected using Chauvenet's Criterion.
logical flag indicating whether the function should produce graphical output or return numerical values to the user.
a vector with additional values to be displayed as different background colours of the plot symbols.
label of the colour legend
Fill colour for the measurements or age estimates. This can
either be a single colour or multiple colours to form a colour
ramp (to be used if levels!=NA):
a single colour: rgb(0,1,0,0.5), '#FF000080',
'white', etc.;
multiple colours: c(rbg(1,0,0,0.5),
rgb(0,1,0,0.5)), c('#FF000080','#00FF0080'),
c('blue','red'), c('blue','yellow','red'), etc.;
a colour palette: rainbow(n=100),
topo.colors(n=100,alpha=0.5), etc.; or
a reversed palette: rev(topo.colors(n=100,alpha=0.5)),
etc.
For empty boxes, set rect.col=NA
if detect.outliers=TRUE, the outliers are
given a different colour.
the number of significant digits of the numerical values reported in the title of the graphical output.
indicates whether the analytical uncertainties of the output are reported in the plot title 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').
plot the aliquots in order of increasing age?
vector with indices of aliquots that should be removed from the weighted mean plot.
vector with indices of aliquots that should be plotted but omitted from the weighted mean calculation.
colour that should be used for the omitted aliquots.
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).
the age (in Ma) below which the
\(^{206}\)Pb/\(^{238}\)U age and above which the
\(^{207}\)Pb/\(^{206}\)Pb age is used. This parameter is
only used if type=4.
discordance cutoff filter. This is an object of
class discfilter
propagate decay constant uncertainties?
common lead correction:
0: none
1: use the Pb-composition stored in
settings('iratio','Pb207Pb206') (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','Pb206Pb208') and
settings('iratio','Pb207Pb208') (if x has class
UPb and x$format=7 or 8).
2: remove the common Pb by projecting the data along an
inverse isochron. Note: choosing this option introduces a degree of
circularity in the weighted age calculation. In this case the
weighted mean plot just serves as a way to visualise the residuals
of the data around the isochron, and one should be careful not to
over-interpret the numerical output.
3: use the Stacey-Kramers two-stage model to infer the
initial Pb-composition (only applicable if x has class
UPb)
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).
`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',...).
Note that choosing this option introduces a degree of circularity in the weighted age calculation. In this case the weighted mean plot just serves as a way to visualise the residuals of the data around the isochron, and one should be careful not to over-interpret the numerical output.
Let \(\{t_1, ..., t_n\}\) be a set of n age estimates
determined on different aliquots of the same sample, and let
\(\{s[t_1], ..., s[t_n]\}\) be their analytical
uncertainties. IsoplotR then calculates the weighted mean of
these data using one of two methods:
The ordinary error-weighted mean:
\(\mu = \sum(t_i/s[t_i]^2)/\sum(1/s[t_i]^2)\)
A random effects model with two sources of variance:
\(\log[t_i] \sim N(\log[\mu], \sigma^2 = (s[t_i]/t_i)^2 + \omega^2 )\)
where \(\mu\) is the mean, \(\sigma^2\) is the total variance and \(\omega\) is the 'overdispersion'. This equation can be solved for \(\mu\) and \(\omega\) by the method of maximum likelihood.
IsoplotR uses a modified version of Chauvenet's criterion for outlier detection:
Compute the error-weighted mean (\(\mu\)) of the \(n\) age determinations \(t_i\) using their analytical uncertainties \(s[t_i]\)
For each \(t_i\), compute the probability \(p_i\) that that \(|t-\mu|>|t_i-\mu|\) for \(t \sim N(\mu, s[t_i]^2 MSWD)\) (ordinary weighted mean) or \(\log[t] \sim N(\log[\mu],s[t_i]^2+\omega^2)\) (random effects model)
Let \(p_j \equiv \min(p_1, ..., p_n)\). If \(p_j<0.05/n\), then reject the j\(^{th}\) date, reduce \(n\) by one (i.e., \(n \rightarrow n-1\)) and repeat steps 1 through 3 until the surviving dates pass the third step.
If the analytical uncertainties are small compared to the scatter between the dates (i.e. if \(\omega \gg s[t]\) for all \(i\)), then this generalised algorithm reduces to the conventional Chauvenet criterion. If the analytical uncertainties are large and the data do not exhibit any overdispersion, then the heuristic outlier detection method is equivalent to Ludwig (2003)'s `2-sigma' method.
The uncertainty budget of the weighted mean does not include the uncertainty of the initial daughter correction (if any). This uncertainty is neither a purely systematic nor a purely random uncertainty and cannot easily be propagated with conventional geochronological data processing algorithms. This caveat is especially pertinent to chronometers whose initial daughter composition is determined by isochron regression. You may note that the uncertainties of the weighted mean are usually much smaller than those of the isochron. In this case the isochron errors are more meaningful, and the weighted mean plot should just be used to inspect the residuals of the data around the isochron.
central
ages <- c(251.9,251.59,251.47,251.35,251.1,251.04,250.79,250.73,251.22,228.43)
errs <- c(0.28,0.28,0.63,0.34,0.28,0.63,0.28,0.4,0.28,0.33)
weightedmean(cbind(ages,errs))
attach(examples)
weightedmean(LudwigMean)
Run the code above in your browser using DataLab