weightedmean: Calculate the weighted mean age

Description Usage Arguments Details Value See Also Examples

View source: R/weightedmean.R

Description

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.

Usage

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
weightedmean(x, ...)

## Default S3 method:
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,
  alpha = 0.05,
  ranked = FALSE,
  hide = NULL,
  omit = NULL,
  omit.col = NA,
  ...
)

## S3 method for class '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,
  alpha = 0.05,
  cutoff.disc = discfilter(),
  exterr = TRUE,
  ranked = FALSE,
  common.Pb = 0,
  hide = NULL,
  omit = NULL,
  omit.col = NA,
  ...
)

## S3 method for class '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,
  alpha = 0.05,
  exterr = TRUE,
  common.Pb = 2,
  ranked = FALSE,
  hide = NULL,
  omit = NULL,
  omit.col = NA,
  ...
)

## S3 method for class '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,
  alpha = 0.05,
  ranked = FALSE,
  i2i = TRUE,
  detritus = 0,
  hide = NULL,
  omit = NULL,
  omit.col = NA,
  ...
)

## S3 method for class '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,
  alpha = 0.05,
  exterr = TRUE,
  ranked = FALSE,
  i2i = FALSE,
  hide = NULL,
  omit = NULL,
  omit.col = NA,
  ...
)

## S3 method for class '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,
  alpha = 0.05,
  exterr = TRUE,
  ranked = FALSE,
  i2i = FALSE,
  hide = NULL,
  omit = NULL,
  omit.col = NA,
  ...
)

## S3 method for class '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,
  alpha = 0.05,
  exterr = TRUE,
  ranked = FALSE,
  i2i = TRUE,
  hide = NULL,
  omit = NULL,
  omit.col = NA,
  ...
)

## S3 method for class '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,
  alpha = 0.05,
  exterr = TRUE,
  ranked = FALSE,
  i2i = TRUE,
  hide = NULL,
  omit = NULL,
  omit.col = NA,
  ...
)

## S3 method for class '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,
  alpha = 0.05,
  exterr = TRUE,
  ranked = FALSE,
  i2i = TRUE,
  hide = NULL,
  omit = NULL,
  omit.col = NA,
  ...
)

## S3 method for class '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,
  alpha = 0.05,
  exterr = TRUE,
  i2i = TRUE,
  ranked = FALSE,
  hide = NULL,
  omit = NULL,
  omit.col = NA,
  ...
)

## S3 method for class '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,
  alpha = 0.05,
  exterr = TRUE,
  i2i = TRUE,
  ranked = FALSE,
  hide = NULL,
  omit = NULL,
  omit.col = NA,
  ...
)

## S3 method for class '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,
  alpha = 0.05,
  ranked = FALSE,
  hide = NULL,
  omit = NULL,
  omit.col = NA,
  ...
)

## S3 method for class '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,
  alpha = 0.05,
  exterr = TRUE,
  ranked = FALSE,
  hide = NULL,
  omit = NULL,
  omit.col = NA,
  ...
)

Arguments

x

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

from

minimum y-axis limit. Setting from=NA scales the plot automatically.

to

maximum y-axis limit. Setting to=NA scales the plot automatically.

random.effects

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.

detect.outliers

logical flag indicating whether outliers should be detected and rejected using Chauvenet's Criterion.

plot

logical flag indicating whether the function should produce graphical output or return numerical values to the user.

levels

a vector with additional values to be displayed as different background colours of the plot symbols.

clabel

label of the colour legend

rect.col

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

outlier.col

if detect.outliers=TRUE, the outliers are given a different colour.

sigdig

the number of significant digits of the numerical values reported in the title of the graphical output.

alpha

the confidence limits of the error bars/rectangles.

ranked

plot the aliquots in order of increasing age?

hide

vector with indices of aliquots that should be removed from the weighted mean plot.

omit

vector with indices of aliquots that should be plotted but omitted from the weighted mean calculation.

omit.col

colour that should be used for the omitted aliquots.

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 age 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

exterr

propagate decay constant uncertainties?

common.Pb

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','Pb208Pb206') and settings('iratio','Pb208Pb207') (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)

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',...).

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.

detritus

detrital ^{230}Th correction (only applicable when x$format=1 or 2).

0: no correction

1: project the data along an isochron fit

2: correct the data using an assumed initial ^{230}Th/^{232}Th-ratio for the detritus.

3: correct the data using the measured present day ^{230}Th/^{238}U, ^{232}Th/^{238}U and ^{234}U/^{238}U-ratios in the detritus.

Details

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:

  1. The ordinary error-weighted mean:

    μ = ∑(t_i/s[t_i]^2)/∑(1/s[t_i]^2)

  2. A random effects model with two sources of variance:

    \log[t_i] \sim N(\log[μ], σ^2 = (s[t_i]/t_i)^2 + ω^2 )

    where μ is the mean, σ^2 is the total variance and ω is the 'overdispersion'. This equation can be solved for μ and ω by the method of maximum likelihood.

IsoplotR uses a modified version of Chauvenet's criterion for outlier detection:

  1. Compute the error-weighted mean (μ) of the n age determinations t_i using their analytical uncertainties s[t_i]

  2. For each t_i, compute the probability p_i that that |t-μ|>|t_i-μ| for t \sim N(μ, s[t_i]^2 MSWD) (ordinary weighted mean) or \log[t] \sim N(\log[μ],s[t_i]^2+ω^2) (random effects model)

  3. 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 ω \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.

Value

Returns a list with the following items:

mean

a 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.

ci[t]: the 100(1-α)\% confidence interval for t

disp

a three-element vector with the (over)dispersion and the lower and upper half-widths of its 100(1-α)\% confidence interval.

mswd

the Mean Square of the Weighted Deviates (a.k.a. ‘reduced Chi-square’ statistic)

df

the number of degrees of freedom of the Chi-square test for homogeneity (df=n-1, where n is the number of samples).

p.value

the p-value of a Chi-square test with df degrees of freedom, testing the null hypothesis that the underlying population is not overdispersed.

valid

vector of logical flags indicating which steps are included into the weighted mean calculation

plotpar

list of plot parameters for the weighted mean diagram, including mean (the mean value), ci (a grey rectangle with the 100[1-α]% confidence interval ignoring systematic errors), ci.exterr (a grey rectangle with the 100[1-α]% confidence interval including systematic errors), dash1 and dash2 (lines marking the 100[1-α]% confidence interval augmented by √{mswd} overdispersion if random.effects=FALSE), and marking the 100[1-α]% confidence limits of a normal distribution whose standard deviation equals the overdispersion parameter if random.effects=TRUE).

See Also

central

Examples

1
2
3
4
5
6
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)

IsoplotR documentation built on July 10, 2021, 1:06 a.m.