Models the data as a Normal distribution with two sources of variance. Estimates the mean and `overdispersion' using the method of Maximum Likelihood. Computes the MSWD of a Normal fit without overdispersion. Implements a modified Chauvenet Criterion to detect and reject outliers. Only propagates the analytical uncertainty associated with decay constants and \(\zeta\) and J-factors after computing the weighted mean isotopic composition.
weightedmean(x, ...)# S3 method for default
weightedmean(x, detect.outliers = TRUE, plot = TRUE,
rect.col = rgb(0, 1, 0, 0.5), outlier.col = rgb(0, 1, 1, 0.5),
sigdig = 2, alpha = 0.05, ...)
# S3 method for UPb
weightedmean(x, detect.outliers = TRUE, plot = TRUE,
rect.col = rgb(0, 1, 0, 0.5), outlier.col = rgb(0, 1, 1, 0.5),
sigdig = 2, type = 4, cutoff.76 = 1100, cutoff.disc = c(-15, 5),
alpha = 0.05, exterr = TRUE, common.Pb = 0, ...)
# S3 method for PbPb
weightedmean(x, detect.outliers = TRUE, plot = TRUE,
rect.col = rgb(0, 1, 0, 0.5), outlier.col = rgb(0, 1, 1, 0.5),
sigdig = 2, alpha = 0.05, exterr = TRUE, common.Pb = 1, ...)
# S3 method for ThU
weightedmean(x, detect.outliers = TRUE, plot = TRUE,
rect.col = rgb(0, 1, 0, 0.5), outlier.col = rgb(0, 1, 1, 0.5),
sigdig = 2, alpha = 0.05, i2i = TRUE, ...)
# S3 method for ArAr
weightedmean(x, detect.outliers = TRUE, plot = TRUE,
rect.col = rgb(0, 1, 0, 0.5), outlier.col = rgb(0, 1, 1, 0.5),
sigdig = 2, alpha = 0.05, exterr = TRUE, i2i = FALSE, ...)
# S3 method for ReOs
weightedmean(x, detect.outliers = TRUE, plot = TRUE,
rect.col = rgb(0, 1, 0, 0.5), outlier.col = rgb(0, 1, 1, 0.5),
sigdig = 2, alpha = 0.05, exterr = TRUE, i2i = TRUE, ...)
# S3 method for SmNd
weightedmean(x, detect.outliers = TRUE, plot = TRUE,
rect.col = rgb(0, 1, 0, 0.5), outlier.col = rgb(0, 1, 1, 0.5),
sigdig = 2, alpha = 0.05, exterr = TRUE, i2i = TRUE, ...)
# S3 method for RbSr
weightedmean(x, detect.outliers = TRUE, plot = TRUE,
rect.col = rgb(0, 1, 0, 0.5), outlier.col = rgb(0, 1, 1, 0.5),
sigdig = 2, alpha = 0.05, exterr = TRUE, i2i = TRUE, ...)
# S3 method for LuHf
weightedmean(x, detect.outliers = TRUE, plot = TRUE,
rect.col = rgb(0, 1, 0, 0.5), outlier.col = rgb(0, 1, 1, 0.5),
sigdig = 2, alpha = 0.05, exterr = TRUE, i2i = TRUE, ...)
# S3 method for UThHe
weightedmean(x, detect.outliers = TRUE, plot = TRUE,
rect.col = rgb(0, 1, 0, 0.5), outlier.col = rgb(0, 1, 1, 0.5),
sigdig = 2, alpha = 0.05, ...)
# S3 method for fissiontracks
weightedmean(x, detect.outliers = TRUE, plot = TRUE,
rect.col = rgb(0, 1, 0, 0.5), outlier.col = rgb(0, 1, 1, 0.5),
sigdig = 2, alpha = 0.05, exterr = TRUE, ...)
a two column matrix of values (first column) and their
standard errors (second column) OR an object of class
UPb, PbPb, ArAr, ReOs, SmNd,
RbSr, LuHf, ThU, fissiontracks or
UThHe
optional arguments
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.
the fill colour of the rectangles used to show the measurements or age estimates.
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.
the confidence limits of the error bars/rectangles.
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), or the (Wetherill) concordia age (type=5)
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.
two element vector with the maximum and minimum
percentage discordance allowed between the
\(^{207}\)Pb/\(^{235}\)U and \(^{206}\)Pb/\(^{238}\)U
age (if \(^{206}\)Pb/\(^{238}\)U < cutoff.76) or
between the \(^{206}\)Pb/\(^{238}\)U and
\(^{207}\)Pb/\(^{206}\)Pb age (if
\(^{206}\)Pb/\(^{238}\)U > cutoff.76). Set
cutoff.disc=NA if you do not want to use this filter.
propagate decay constant uncertainties?
apply a common lead correction using one of three methods:
1: use the isochron intercept as the initial Pb-composition
2: use the Stacey-Kramer two-stage model to infer the initial
Pb-composition
3: use the Pb-composition stored in
settings('iratio','Pb206Pb204') and
settings('iratio','Pb207Pb204')
`isochron to intercept': calculates the initial (aka
`inherited', `excess', or `common')
\(^{40}\)Ar/\(^{36}\)Ar, \(^{207}\)Pb/\(^{204}\)Pb,
\(^{87}\)Sr/\(^{86}\)Sr, \(^{143}\)Nd/\(^{144}\)Nd,
\(^{187}\)Os/\(^{188}\)Os or \(^{176}\)Hf/\(^{177}\)Hf
ratio from an isochron fit. Setting i2i to FALSE
uses the default values stored in
settings('iratio',...). When applied to data of class
ThU, setting i2i to TRUE applies a
detrital Th-correction.
Returns a list with the following items:
a three element vector with:
x: the weighted mean
s[x]: the estimated analytical uncertainty of x
ci[x]: the studentised \(100(1-\alpha)\%\) confidence
interval for x.
a two element vector with the (over)dispersion and its corresponding \(100(1-\alpha)\%\) confidence interval.
the degrees of freedom for the Chi-square test (\(n-2\))
the \(100(1-\alpha/2)\) percentile of a
t-distribution with df+1 degrees of freedom
the Mean Square of the Weighted Deviates (a.k.a. `reduced Chi-square' statistic)
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
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 assuming a Normal distribution with two sources of
variance:
\(t_i \sim N(\mu, \sigma^2 = s[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(0,\sqrt{s[t_i]^2+\omega^2) }\)
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 analtyical 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.
Ludwig, K. R. User's manual for Isoplot 3.00: a geochronological toolkit for Microsoft Excel. Berkeley Geochronology Center Special Publication, 2003.
# NOT RUN {
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))
data(examples)
weightedmean(examples$LudwigMean)
# }
Run the code above in your browser using DataLab