weightedmean: Calculate the weighted mean age

Description Usage Arguments Details Value See Also Examples

View source: R/weightedmean.R

Description

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 ζ and J-factors after computing the weighted mean isotopic composition.

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
weightedmean(x, ...)

## Default S3 method:
weightedmean(x, from = NA, to = NA,
  random.effects = TRUE, 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 = TRUE,
  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, cutoff.disc = c(-15, 5), alpha = 0.05,
  exterr = TRUE, ranked = FALSE, common.Pb = 0, hide = NULL,
  omit = NULL, omit.col = NA, ...)

## S3 method for class 'PbPb'
weightedmean(x, random.effects = TRUE,
  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 = 1, ranked = FALSE, hide = NULL,
  omit = NULL, omit.col = NA, ...)

## S3 method for class 'ThU'
weightedmean(x, random.effects = TRUE,
  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, Th02 = c(0, 0),
  Th02U48 = c(0, 0, 1e+06, 0, 0, 0, 0, 0, 0), hide = NULL,
  omit = NULL, omit.col = NA, ...)

## S3 method for class 'ArAr'
weightedmean(x, random.effects = TRUE,
  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 = TRUE,
  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 'ReOs'
weightedmean(x, random.effects = TRUE,
  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 = TRUE,
  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 = TRUE,
  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 = TRUE,
  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 = TRUE,
  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 = TRUE,
  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, 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

a vector of two fill colours used to show the measurements or age estimates. If levels=NA, then only the first colour is used. If levels is a vector of numbers, then bg is used to construct a colour ramp.

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), or the (Wetherill) concordia age (type=5)

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

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.

exterr

propagate decay constant uncertainties?

common.Pb

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

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 or ^{176}Hf/^{177}Hf ratio from an isochron fit. Setting i2i to FALSE uses the default values stored in settings('iratio',...).

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.

Th02

2-element vector with the assumed initial ^{230}Th/^{232}Th-ratio of the detritus and its standard error. Only used if detritus==2

Th02U48

9-element vector with the measured composition of the detritus, containing X=0/8, sX, Y=2/8, sY, Z=4/8, sZ, rXY, rXZ, rYZ. Only used if isochron==FALSE and detritus==3

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 assuming a Normal distribution with two sources of variance:

t_i \sim N(μ, σ^2 = s[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(0,√{s[t_i]^2+ω^2) }

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

Value

Returns a list with the following items:

mean

a three element vector with:

x: the weighted mean

s[x]: the standard error of the weighted mean

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

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

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

data(examples)
weightedmean(examples$LudwigMean)

IsoplotR documentation built on Dec. 9, 2018, 1:04 a.m.