Package 'IsoplotR'

Calculate isotopic ages

Description

Calculates U-Pb, Pb-Pb, Th-Pb, Ar-Ar, K-Ca, Re-Os, Sm-Nd, Rb-Sr, Lu-Hf, U-Th-He, Th-U and fission track ages and propagates their analytical uncertainties. Includes options for single grain, isochron and concordia_ages.

Usage

age(x, ...)

## Default S3 method:
age(
  x,
  method = "U238-Pb206",
  oerr = 1,
  sigdig = NA,
  exterr = FALSE,
  J = c(NA, NA),
  zeta = c(NA, NA),
  rhoD = c(NA, NA),
  d = diseq(),
  ...
)

## S3 method for class 'UPb'
age(
  x,
  type = 1,
  exterr = FALSE,
  i = NULL,
  oerr = 1,
  sigdig = NA,
  common.Pb = 0,
  discordance = discfilter(),
  ...
)

## S3 method for class 'PbPb'
age(
  x,
  isochron = TRUE,
  common.Pb = 2,
  exterr = FALSE,
  i = NULL,
  oerr = 1,
  sigdig = NA,
  projerr = FALSE,
  ...
)

## S3 method for class 'ArAr'
age(
  x,
  isochron = FALSE,
  i2i = TRUE,
  exterr = FALSE,
  i = NULL,
  oerr = 1,
  sigdig = NA,
  projerr = FALSE,
  ...
)

## S3 method for class 'KCa'
age(
  x,
  isochron = FALSE,
  i2i = TRUE,
  exterr = FALSE,
  i = NULL,
  oerr = 1,
  sigdig = NA,
  projerr = FALSE,
  ...
)

## S3 method for class 'UThHe'
age(x, isochron = FALSE, central = FALSE, i = NULL, oerr = 1, sigdig = NA, ...)

## S3 method for class 'fissiontracks'
age(x, central = FALSE, i = NULL, oerr = 1, sigdig = NA, exterr = FALSE, ...)

## S3 method for class 'ThU'
age(
  x,
  isochron = FALSE,
  Th0i = 0,
  exterr = FALSE,
  i = NULL,
  oerr = 1,
  sigdig = NA,
  ...
)

## S3 method for class 'ThPb'
age(
  x,
  isochron = TRUE,
  i2i = TRUE,
  exterr = FALSE,
  i = NULL,
  oerr = 1,
  sigdig = NA,
  projerr = FALSE,
  ...
)

## S3 method for class 'ReOs'
age(
  x,
  isochron = TRUE,
  i2i = TRUE,
  exterr = FALSE,
  i = NULL,
  oerr = 1,
  sigdig = NA,
  projerr = FALSE,
  ...
)

## S3 method for class 'SmNd'
age(
  x,
  isochron = TRUE,
  i2i = TRUE,
  exterr = FALSE,
  i = NULL,
  oerr = 1,
  sigdig = NA,
  projerr = FALSE,
  ...
)

## S3 method for class 'RbSr'
age(
  x,
  isochron = TRUE,
  i2i = TRUE,
  exterr = FALSE,
  i = NULL,
  oerr = 1,
  sigdig = NA,
  projerr = FALSE,
  ...
)

## S3 method for class 'LuHf'
age(
  x,
  isochron = TRUE,
  i2i = TRUE,
  exterr = FALSE,
  i = NULL,
  oerr = 1,
  sigdig = NA,
  projerr = FALSE,
  ...
)

## S3 method for class 'PD'
age(
  x,
  nuclide,
  isochron = TRUE,
  i2i = TRUE,
  exterr = FALSE,
  i = NULL,
  oerr = 1,
  sigdig = NA,
  projerr = FALSE,
  ...
)

Arguments

x	can be: a scalar containing an isotopic ratio, a two element vector containing an isotopic ratio and its standard error, or the spontaneous and induced track densities Ns and Ni, a two element vector containing Ar40Ar39, s[Ar40Ar39], a two element vector containing K40Ca40 and s[K40Ca40], a six element vector containing U, s[U], Th, s[Th], He and s[He], an eight element vector containing U, s[U], Th, s[Th], He, s[He], Sm and s[Sm] a two element vector containing Sr87Rb87 and s[Sr87Rb87] a two element vector containing Os187Re187 and s[Os187Re187] a two element vector containing Nd143Sm147 and s[Nd144Sm147] a two element vector containing Hf176Lu176 and s[Hf176Lu176] a five element vector containing Th230U238, s[Th230/U238], U234U238, s[U234U238] and cov[Th230U238,U234U238] OR an object of class UPb, PbPb, ThPb, ArAr, KCa, ThU, RbSr, SmNd, ReOs, LuHf, UThHe or fissiontracks.
...	additional arguments
method	one of either 'U238-Pb206', 'U235-Pb207', 'Pb207-Pb206', 'Th232-Pb208', 'Ar-Ar', 'K-Ca', 'Th-U', 'Re-Os', 'Sm-Nd', 'Rb-Sr', 'Lu-Hf', 'U-Th-He' or 'fissiontracks'
oerr	indicates whether the analytical uncertainties of the output are reported 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'). (only used when isochron and central are FALSE)
sigdig	number of significant digits for the uncertainty estimate (only used if type=1, isochron=FALSE and central=FALSE).
exterr	propagate the external (decay constant and calibration factor) uncertainties?
J	two-element vector with the J-factor and its standard error.
zeta	two-element vector with the zeta-factor and its standard error.
rhoD	two-element vector with the track density of the dosimeter glass and its standard error.
d	an object of class diseq.
type	scalar flag indicating whether 1: each U-Pb analysis should be considered separately, 2: all the measurements should be combined to calculate a concordia_age, 3: a discordia_line should be fitted through all the U-Pb analyses using the maximum likelihood algorithm of Ludwig (1998), which assumes that the scatter of the data is solely due to the analytical uncertainties. 4: a discordia_line should be fitted ignoring the analytical uncertainties. 5: a discordia_line should be fitted using a modified maximum likelihood algorithm that accounts for overdispersion by adding a geological (co)variance term.
i	index of a particular aliquot
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','Pb206Pb208') and settings('iratio','Pb207Pb208') (if x has class UPb and x$format=7 or 8). 2: use the isochron intercept as the initial Pb-composition 3: use the Stacey-Kramer two-stage model to infer the initial Pb-composition
discordance	discordance calculator. This is an object of class discfilter, or a two element list containing: option: one of 1 or 't' (absolute age filter); 2 or 'r' (relative age filter); 3 or 'sk' (Stacey-Kramers common Pb filter); 4 or 'a' (perpendicular Aitchison distance); 5 or 'c' (concordia distance); 6 or 'p' (p-value of concordance); or NA (omit the discordance from the output) before: logical flag indicating whether the discordance should be calculated before (TRUE) or after (FALSE) the common-Pb correction.
isochron	logical flag indicating whether each analysis should be considered separately (isochron=FALSE) or an isochron age should be calculated from all analyses together (isochron=TRUE).
projerr	logical. If TRUE, propagates the uncertainty of the non-radiogenic isotope correction (the 'projection error') into the age uncertainty. Note that the resulting single grain age uncertainties may be strongly correlated with each other, but these error correlations are not reported in the output.
i2i	‘isochron to intercept’: calculates the initial (aka ‘inherited’, ‘excess’, or ‘common’) $^{40}$Ar/$^{36}$Ar, $^{40}$Ca/$^{44}$Ca, $^{87}$Sr/$^{86}$Sr, $^{143}$Nd/$^{144}$Nd, $^{187}$Os/$^{188}$Os, $^{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',...).
central	logical flag indicating whether each analysis should be considered separately (central=FALSE) or a central age should be calculated from all analyses together (central=TRUE).
Th0i	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).
nuclide	a text string corresponding to a valid entry for settings('lambda',...)

Value

1. if x is a scalar or a vector, returns the age using the geochronometer given by method and its standard error.

2. if x has class UPb and type=1, returns a table with the following columns: t.75, err[t.75], t.68, err[t.68], t.76, err[t.76], (t.82, err[t.82],) t.conc, err[t.conc], (disc) or err[p.conc],) containing the $^{207}$Pb/$^{235}$U-age and standard error, the $^{206}$Pb/$^{238}$U-age and standard error, the $^{207}$Pb/$^{206}$Pb-age and standard error, (the $^{208}$Pb/$^{232}$Th-age and standard error,) the single grain concordia_age and standard error, (and the % discordance or p-value for concordance,) respectively.

3. if x has class UPb and type=2, 3, 4 or 5, returns the output of the concordia function.

4. if x has class PbPb, ThPb, ArAr, KCa, RbSr, SmNd, ReOs, LuHf, ThU or UThHe and isochron=FALSE, returns a table of Pb-Pb, Th-Pb, Ar-Ar, K-Ca, Rb-Sr, Sm-Nd, Re-Os, Lu-Hf, Th-U or U-Th-He ages and their standard errors.

5. if x has class ThU and isochron=FALSE, returns a 5-column table with the Th-U ages, their standard errors, the initial $^{234}$U/$^{238}$U-ratios, their standard errors, and the correlation coefficient between the ages and the initial ratios.

6. if x has class PbPb, ThPb, ArAr, KCa, RbSr, SmNd, ReOs, LuHf, UThHe or ThU and isochron=TRUE, returns the output of the isochron function.

7. if x has class fissiontracks and central=FALSE, returns a table of fission track ages and standard errors.

8. if x has class fissiontracks or UThHe and central=TRUE, returns the output of the central function.

See Also

concordia, isochron, central

Examples

attach(examples)
tUPb <- age(UPb,type=1)
tconc <- age(UPb,type=2)
tdisc <- age(UPb,type=3)
tArAr <- age(ArAr)
tiso <- age(ArAr,isochron=TRUE,i2i=TRUE)
tcentral <- age(FT1,central=TRUE)

---

Predict isotopic ratios from ages

Description

Groups a set of functions that take one (or more) ages (and their uncertainties) as input and produces the U–Pb, Th–Pb, Pb–Pb, Ar–Ar, K–Ca, Rb–Sr, Sm–Nd, Lu–Hf, Re–Os, concordia or Stacey-Kramers ratios as output.

Usage

age2ratio(
  tt,
  st = 0,
  ratio = "Pb206U238",
  exterr = FALSE,
  d = diseq(),
  J,
  sJ = 0,
  bratio = 0.895
)

Arguments

tt	a scalar or (except when ratio = 'Wetherill', 'Tera-Wasserburg' or 'U-Th-Pb') vector of ages.
st	a scalar or (except when ratio = 'Wetherill', 'Tera-Wasserburg' or 'U-Th-Pb') vector with the standard error(s) of tt. Not used when ratio = 'Stacey-Kramers'.
ratio	one of 'Pb206U238', 'Pb207U235', 'U238Pb206', 'Pb207Pb206', 'Pb208Th232', 'Wetherill', 'Tera-Wasserburg', 'U-Th-Pb', 'Ar40Ar39', 'Ca40K40', 'Hf176Lu176', 'Sr87Rb87', 'Os187Re187', 'Nd143Sm147' or 'Stacey-Kramers'.
exterr	logical. If TRUE, propagates decay constant uncertainties into st. Not used when ratio = 'Stacey-Kramers'.
d	an object of class diseq.
J	the J-factor of the Ar–Ar system (only used if ratio is 'Ar40Ar39').
sJ	the standard error of J (only used if ratio is 'Ar40Ar39').
bratio	branching ratio of $^40$K

Value

If ratio is 'Pb207U235', 'U238Pb206', 'Pb207Pb206', 'Pb208Th232', 'Ar40Ar39', 'Ca40K40', 'Hf176Lu176', 'Sr87Rb87', 'Os187Re187' or 'Nd143Sm147': either a two-element vector or a two-column matrix with the predicted isotopic ratio(s) and its/their standard error(s).

If ratio is 'Wetherill', 'Tera-Wasserburg' or 'U-Th-Pb': a two-element list containing

x: the concordia ratios

cov: the covariance matrix of the concordia ratios

If ratio is 'Stacey-Kramers': a three-column matrix with predicted $^{206}$Pb/$^{204}$Pb, $^{207}$Pb/$^{204}$Pb and $^{208}$Pb/$^{204}$Pb ratios.

Examples

ratios <- c('Pb206U238','Pb207U235','U238Pb206','Pb207Pb206',
            'Pb208Th232','Wetherill','Tera-Wasserburg','U-Th-Pb',
            'Ar40Ar39','Ca40K40','Hf176Lu176','Sr87Rb87',
            'Os187Re187','Nd143Sm147','Stacey-Kramers')
for (ratio in ratios){
     r <- age2ratio(tt=1000,st=10,ratio=ratio,J=1,sJ=0.1)
     print(r)
}

---

Plot a (40Ar/39Ar) release spectrum

Description

Produces a plot of boxes whose widths correspond to the cumulative amount of $^{39}$Ar (or any other variable), and whose heights express the analytical uncertainties. Only propagates the analytical uncertainty associated with decay constants and J-factors after computing the plateau composition.

Usage

agespectrum(x, ...)

## Default S3 method:
agespectrum(
  x,
  oerr = 3,
  plateau = TRUE,
  random.effects = FALSE,
  levels = NA,
  clabel = "",
  plateau.col = c("#00FF0080", "#FF000080"),
  non.plateau.col = "#00FFFF80",
  sigdig = 2,
  line.col = "red",
  lwd = 2,
  xlab = "cumulative fraction",
  ylab = "X",
  hide = NULL,
  omit = NULL,
  ...
)

## S3 method for class 'other'
agespectrum(
  x,
  oerr = 3,
  plateau = TRUE,
  random.effects = FALSE,
  levels = NA,
  clabel = "",
  plateau.col = c("#00FF0080", "#FF000080"),
  non.plateau.col = "#00FFFF80",
  sigdig = 2,
  line.col = "red",
  lwd = 2,
  xlab = "cumulative fraction",
  ylab = "X",
  hide = NULL,
  omit = NULL,
  ...
)

## S3 method for class 'ArAr'
agespectrum(
  x,
  oerr = 3,
  plateau = TRUE,
  random.effects = FALSE,
  levels = NA,
  clabel = "",
  plateau.col = c("#00FF0080", "#FF000080"),
  non.plateau.col = "#00FFFF80",
  sigdig = 2,
  exterr = FALSE,
  line.col = "red",
  lwd = 2,
  i2i = FALSE,
  hide = NULL,
  omit = NULL,
  ...
)

Arguments

x	a three-column matrix whose first column gives the amount of $^{39}$Ar in each aliquot, and whose second and third columns give the age and its uncertainty. OR an object of class ArAr
...	optional parameters to the generic plot function
oerr	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').
plateau	logical flag indicating whether a plateau age should be calculated. If plateau=TRUE, the function computes the weighted mean of the largest succession of steps that pass the Chi-square test for age homogeneity. If TRUE, it returns a list with plateau parameters.
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.
levels	a vector with additional values to be displayed as different background colours of the plot symbols.
clabel	label of the colour legend
plateau.col	Fill colours of the rectangles used to mark the steps belonging to the age plateau. 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 plateau.col=NA
non.plateau.col	if plateau=TRUE, the steps that do NOT belong to the plateau are given a different colour.
sigdig	the number of significant digits of the numerical values reported in the title of the graphical output.
line.col	colour of the average age line
lwd	width of the average age line
xlab	x-axis label
ylab	y-axis label
hide	vector with indices of aliquots that should be removed from the plot.
omit	vector with indices of aliquots that should be plotted but omitted from age plateau calculation
exterr	propagate the external (decay constant and calibration factor) uncertainties?
i2i	‘isochron to intercept’: calculates the initial (aka ‘inherited’, ‘excess’, or ‘common’) $^{40}$Ar/$^{36}$Ar ratio from an isochron fit. Setting i2i to FALSE uses the default values stored in settings('iratio',...)

Details

IsoplotR defines the ‘plateau age’ as the weighted mean age (using a random effects model with two sources of dispersion) of the longest sequence (in terms of cumulative $^{39}$Ar content) of consecutive heating steps that pass the modified Chauvenet criterion (see weightedmean). Note that this definition is different (and simpler) than the one used by Isoplot (Ludwig, 2003). However, it is important to mention that all definitions of an age plateau are heuristic by nature and should not be used for quantitative inference. It is possible (and likely) that the plateau steps exhibit significant overdispersion. This overdispersion can be manually reduced by removing individual heating steps with the optional omit argument.

Value

If plateau=TRUE, returns a list containing the output of the weightedmean function, plus the following items:

fract

the fraction of $^{39}$Ar contained in the plateau

i

indices of the steps that are retained for the plateau age calculation

See Also

weightedmean

Examples

attach(examples)
par(mfrow=c(2,1))
agespectrum(ArAr)
# removing the first 6 steps yields the longest plateau
# that passes the chi-square test for homogeneity
agespectrum(ArAr,omit=1:6)

---

Plot continuous data as cumulative age distributions

Description

Plot a dataset as a Cumulative Age Distribution (CAD), also known as a ‘empirical cumulative distribution function’.

Usage

cad(x, ...)

## Default S3 method:
cad(
  x,
  pch = NA,
  verticals = TRUE,
  xlab = "age [Ma]",
  col = "black",
  hide = NULL,
  ...
)

## S3 method for class 'other'
cad(
  x,
  pch = NA,
  verticals = TRUE,
  xlab = "age [Ma]",
  col = "black",
  hide = NULL,
  ...
)

## S3 method for class 'detritals'
cad(
  x,
  pch = NA,
  verticals = TRUE,
  xlab = "age [Ma]",
  col = "rainbow",
  hide = NULL,
  ...
)

## S3 method for class 'UPb'
cad(
  x,
  pch = NA,
  verticals = TRUE,
  xlab = "age [Ma]",
  col = "black",
  type = 4,
  cutoff.76 = 1100,
  cutoff.disc = discfilter(),
  common.Pb = 0,
  hide = NULL,
  ...
)

## S3 method for class 'PbPb'
cad(
  x,
  pch = NA,
  verticals = TRUE,
  xlab = "age [Ma]",
  col = "black",
  common.Pb = 1,
  hide = NULL,
  ...
)

## S3 method for class 'ArAr'
cad(
  x,
  pch = NA,
  verticals = TRUE,
  xlab = "age [Ma]",
  col = "black",
  i2i = FALSE,
  hide = NULL,
  ...
)

## S3 method for class 'KCa'
cad(
  x,
  pch = NA,
  verticals = TRUE,
  xlab = "age [Ma]",
  col = "black",
  i2i = FALSE,
  hide = NULL,
  ...
)

## S3 method for class 'ThPb'
cad(
  x,
  pch = NA,
  verticals = TRUE,
  xlab = "age [Ma]",
  col = "black",
  i2i = TRUE,
  hide = NULL,
  ...
)

## S3 method for class 'ThU'
cad(
  x,
  pch = NA,
  verticals = TRUE,
  xlab = "age [ka]",
  col = "black",
  Th0i = 0,
  hide = NULL,
  ...
)

## S3 method for class 'ThPb'
cad(
  x,
  pch = NA,
  verticals = TRUE,
  xlab = "age [Ma]",
  col = "black",
  i2i = TRUE,
  hide = NULL,
  ...
)

## S3 method for class 'ReOs'
cad(
  x,
  pch = NA,
  verticals = TRUE,
  xlab = "age [Ma]",
  col = "black",
  i2i = TRUE,
  hide = NULL,
  ...
)

## S3 method for class 'SmNd'
cad(
  x,
  pch = NA,
  verticals = TRUE,
  xlab = "age [Ma]",
  col = "black",
  i2i = TRUE,
  hide = NULL,
  ...
)

## S3 method for class 'RbSr'
cad(
  x,
  pch = NA,
  verticals = TRUE,
  xlab = "age [Ma]",
  col = "black",
  i2i = TRUE,
  hide = NULL,
  ...
)

## S3 method for class 'LuHf'
cad(
  x,
  pch = NA,
  verticals = TRUE,
  xlab = "age [Ma]",
  col = "black",
  i2i = TRUE,
  hide = NULL,
  ...
)

## S3 method for class 'UThHe'
cad(
  x,
  pch = NA,
  verticals = TRUE,
  xlab = "age [Ma]",
  col = "black",
  hide = NULL,
  ...
)

## S3 method for class 'fissiontracks'
cad(
  x,
  pch = NA,
  verticals = TRUE,
  xlab = "age [Ma]",
  col = "black",
  hide = NULL,
  ...
)

Arguments

x	a numerical vector OR an object of class UPb, PbPb, ThPb, ArAr, KCa, UThHe, fissiontracks, ReOs, RbSr, SmNd, LuHf, ThU or detritals
...	optional arguments to the generic plot function
pch	plot character to mark the beginning of each CAD step
verticals	logical flag indicating if the horizontal lines of the CAD should be connected by vertical lines
xlab	x-axis label
col	if x has class detritals, the name of one of R's built-in colour palettes (e.g., 'heat.colors', 'terrain.colors', 'topo.colors', 'cm.colors'), OR a vector with the names or codes of two colours to use as the start and end of a colour ramp (e.g. col=c('yellow','blue')). For all other data formats, the name or code for a colour to give to a single sample dataset
hide	vector with indices of aliquots that should be removed from the plot.
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.
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','Pb206Pb208') and settings('iratio','Pb207Pb208') (if x has class UPb and x$format=7 or 8). 2: use the isochron intercept as the initial Pb-composition 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',...).
Th0i	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).

Details

Empirical cumulative distribution functions or cumulative age distributions are the most straightforward way to visualise the probability distribution of multiple dates. Suppose that we have a set of $n$ dates $t_i$. The CAD is a step function that sets out the rank order of the dates against their numerical value:

$CAD(t) = \sum_i 1(t<t_i)/n$

where 1($\ast$) = 1 if $\ast$ is true and 1($\ast$) = 0 if $\ast$ is false. CADs have two desirable properties (Vermeesch, 2007). First, they do not require any pre-treatment or smoothing of the data. This is not the case for histograms or kernel density estimates. Second, it is easy to superimpose several CADs on the same plot. This facilitates the intercomparison of multiple samples. The interpretation of CADs is straightforward but not very intuitive. The prominence of individual age components is proportional to the steepness of the CAD. This is different from probability density estimates such as histograms, in which such components stand out as peaks.

References

Vermeesch, P., 2007. Quantitative geomorphology of the White Mountains (California) using detrital apatite fission track thermochronology. Journal of Geophysical Research: Earth Surface, 112(F3).

See Also

kde, radialplot

Examples

attach(examples)
cad(DZ,verticals=FALSE,pch=20)

---

Fits random effects models to overdispersed datasets

Description

Computes the logratio mean composition of a continuous mixture of fission track or U-Th-He data and returns the corresponding age and fitting parameters. 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

central(x, ...)

## Default S3 method:
central(x, ...)

## S3 method for class 'UThHe'
central(x, compositional = FALSE, model = 1, ...)

## S3 method for class 'fissiontracks'
central(x, exterr = FALSE, ...)

Arguments

x	an object of class UThHe or fissiontracks, OR a 2-column matrix with (strictly positive) values and uncertainties
...	optional arguments
compositional	logical. If TRUE, calculates the 'barycentric' U-Th-He, age, i.e. the age corresponding to the weighted mean logratio composition.
model	only relevant if compositional is TRUE. If the scatter between the data points is solely caused by the analytical uncertainty, then the MSWD value should be approximately equal to one. There are three strategies to deal with the case where MSWD>1.choose one of the following statistical models: 1: assume that the analytical uncertainties have been underestimated by a factor $\sqrt{MSWD}$. 2: ignore the analytical uncertainties. 3: attribute any excess dispersion to the presence of geological uncertainty, which manifests itself as an added (co)variance term.
exterr	include the zeta or decay constant uncertainty into the error propagation for the central age?

Details

The central age assumes that the observed age distribution is the combination of two sources of scatter: analytical uncertainty and true geological dispersion.

1. For fission track data, the analytical uncertainty is assumed to obey Binomial counting statistics and the geological dispersion is assumed to follow a lognormal distribution.

2. For U-Th-He data, the U-Th-(Sm)-He compositions and uncertainties are assumed to follow a logistic normal distribution.

3. For all other data types, both the analytical uncertainties and the true ages are assumed to follow lognormal distributions.

The difference between the central age and the weighted mean age is usually small unless the data are imprecise and/or strongly overdispersed.

The uncertainty budget of the central age does not include the uncertainty of the initial daughter correction (if any), for the same reasons as discussed under the weightedmean function.

Value

If x has class UThHe and compositional is TRUE, returns a list containing the following items:

uvw

(if the input data table contains Sm) or uv (if it does not): the mean log[U/He], log[Th/He] (, and log[Sm/He]) composition.

covmat

the covariance matrix of uvw or uv.

mswd

the reduced Chi-square statistic of data concordance, i.e. $mswd=SS/df$, where $SS$ is the sum of squares of the log[U/He]-log[Th/He] compositions.

model

the fitting model.

df

the degrees of freedom ($2n-2$) of the fit (only reported if model=1).

p.value

the p-value of a Chi-square test with df degrees of freedom (only reported if model=1.)

age

a two- or three-element vector with:
t: the 'barycentric' age, i.e. the age corresponding to uvw.
s[t]: the standard error of t.
disp[t]: the standard error of t enhanced by a factor of $\sqrt{mswd}$ (only reported if model=1).

w

the geological overdispersion term. If model=3, this is a two-element vector with the standard deviation of the (assumedly) Normal dispersion and its standard error. w=0 if model<3.

OR, otherwise:

age

a two-element vector with the central age and its standard error.

disp

a two-element vector with the overdispersion (standard deviation) of the excess scatter, and its standard error.

mswd

the reduced Chi-square statistic of data concordance, i.e. $mswd=X^2/df$, where $X^2$ is a Chi-square statistic of the EDM data or ages

df

the degrees of freedom ($n-2$)

p.value

the p-value of a Chi-square test with df degrees of freedom

References

Galbraith, R.F. and Laslett, G.M., 1993. Statistical models for mixed fission track ages. Nuclear Tracks and Radiation Measurements, 21(4), pp.459-470.

Vermeesch, P., 2008. Three new ways to calculate average (U-Th)/He ages. Chemical Geology, 249(3), pp.339-347.

See Also

weightedmean, radialplot, helioplot

Examples

attach(examples)
print(central(UThHe)$age)

---

Confidence intervals

Description

Given a parameter estimate and its standard error, calculate the corresponding 1-sigma, 2-sigma or $100(1-\alpha)\%$ confidence interval, in absolute or relative units.

Usage

ci(x = 0, sx, oerr = 3, df = NULL, absolute = FALSE)

Arguments

x	scalar estimate
sx	scalar or vector with the standard error of x without and (optionally) with $\sqrt{MSWD}$ overdispersion multiplier.
oerr	indicates if the confidence intervals should be reported 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').
df	(optional) number of degrees of freedom. Only used if sx is a vector.
absolute	logical. Returns absolute uncertainties even if oerr is greater than 3. Used for some internal IsoplotR functions.

Details

Several of IsoplotR's plotting functions (including isochron, weightedmean, concordia, radialplot and helioplot) return lists of parameters and their standard errors. For ‘model-1’ fits, if the data pass a Chi-square test of homogeneity, then just one estimate for the standard error is reported. This estimate can be converted to a confidence interval by multiplication with the appropriate quantile of a Normal distribution. Datasets that fail the Chi-square test are said to be ‘overdispersed’ with respect to the analytical uncertainties. The simplest way (‘model-1’) to deal with overdispersion is to inflate the standard error with a $\sqrt{MSWD}$ premultiplier. In this case, IsoplotR returns two estimates of the standard error. To convert the second estimate to a confidence interval requires multiplication with the desired quantile of a t-distribution with the appropriate degrees of freedom.

Value

A scalar or vector of the same size as sx.

Examples

attach(examples)
fit <- isochron(PbPb,plot=FALSE,exterr=FALSE)
err <- ci(x=fit$age[1],sx=fit$age[-1],oerr=5,df=fit$df)
message('age=',signif(fit$age[1],4),'Ma, ',
        '2se=',signif(err[1],2),'%, ',
        '2se(with dispersion)=',signif(err[2],2),'%')

---

Geochronological data classes

Description

S3 classes to store geochronological data generated by read.data or diseq.

Usage

is.UPb(x)

is.PbPb(x)

is.ThPb(x)

is.ArAr(x)

is.KCa(x)

is.PD(x)

is.RbSr(x)

is.SmNd(x)

is.LuHf(x)

is.ReOs(x)

is.ThU(x)

is.UThHe(x)

is.fissiontracks(x)

is.detritals(x)

is.other(x)

is.diseq(x)

as.UPb(x, format = 3, ierr = 1, d = diseq())

as.PbPb(x, format = 1, ierr = 1)

as.ArAr(x, format = 3, ierr = 1)

as.ThPb(x, format = 1, ierr = 1)

as.KCa(x, format = 1, ierr = 1)

as.RbSr(x, format = 1, ierr = 1)

as.ReOs(x, format = 1, ierr = 1)

as.SmNd(x, format = 1, ierr = 1)

as.LuHf(x, format = 1, ierr = 1)

as.ThU(
  x,
  format = 1,
  ierr = 1,
  U8Th2 = 0,
  Th02i = c(0, 0),
  Th02U48 = c(0, 0, 1e+06, 0, 0, 0, 0, 0, 0)
)

as.UThHe(x, ierr = 1)

as.fissiontracks(x, format = 1, ierr = 1)

as.detritals(x)

as.other(x, format = NULL, ierr = 1)

Arguments

x	a data object returned by read.data or diseq.
format	data format. See read.data for details.
ierr	input error. See read.data for details.
d	an object of class diseq.
U8Th2	$^{238}$U/$^{232}$Th activity-ratio of the whole rock. Used to estimate the initial $^{230}$Th/$^{238}$U disequilibrium (for Th-U formats 3 and 4).
Th02i	2-element vector with the assumed initial $^{230}$Th/$^{232}$Th-ratio of the detritus (for Th-U formats 1 and 2) and its standard error.
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.

Details

IsoplotR uses the following S3 classes to store geochronological data: UPb, PbPb, ThPb, KCa, UThHe, fissiontracks, detritals and PD, where the latter is the parent class of the simple parent-daughter chronometers RbSr, SmNd, LuHf and ReOs. All these classes have overloaded versions of the generic length() function and `[` subsetting method.

Additional functions for each class include as.X(x), which converts the data table x to an object of class X; and is.X(x), which checks if x has class X.

• UPb: a list containing:

  x

  a matrix containing the isotopic measurements

  format

  a number between 1 and 8

  d

  an object of class diseq, i.e. the output of the diseq function

• ArAr: a list containing

  x

  a matrix containing the isotopic measurements

  J

  a two-element vector with the J-factor and its uncertainty

  format

  a number between 1 and 3

• ThU: a list containing

  x

  a matrix containing the isotopic measurements

  format

  a number between 1 and 4

  Th02

  a two element vector with the assumed initial $^{230}$Th/$^{232}$Th-ratio of Th-bearing detritus. Only aplicable to formats 1 and 2.

  Th02U48

  9-element vector with the measured composition of Th-bearing detritus

  U8Th2

  the measured $^{238}$U/$^{232}$Th activity ratio of the whole rock. Only applicable to formats 3 and 4

• PbPb, ThPb, KCa, PD, RbSr, SmNd, LuHf, or ReOs: a list containing

  x

  a matrix containing the isotopic measurements

  format

  a number between 1 and 3

• UThHe: a matrix of He, U, Th (and Sm) measurements

• fissiontracks: a list containing

  format

  a number between 1 and 3

  x

  a matrix of spontaneous and induced fission track counts (only included if format=1)

  rhoD

  the track density of the dosimeter glass, extracted from the input data (only included if format=1)

  zeta

  the zeta calibration constant extracted from the input data (only included if format<3)

  Ns

  a list containing the spontaneous fission track counts (only included if format>1)

  U

  a list of lists containing the U-concentration or U/Ca-ratio measurements for each of the analysed grains (only included if format>1)

  sU

  a list of lists containing the standard errors of the U-concentration or U/Ca-ratio measurements for each of the analysed grains (only include if format>1)

  spotSize

  the laser ablation spot size (only included if format>1)

• detritals: a list of named vectors, one for each detrital sample.

• diseq: is a class that contains the output of the diseq function, which stores initial disequilibrium data for U–Pb geochronology.

Value

is.X(x) returns a logical value.

as.X(x) returns an object of class X.

See Also

read.data diseq

Examples

attach(examples)
ns <- length(UPb)
concordia(UPb[-ns,])
if (is.PD(RbSr)) print('RbSr has class PD')

---

Concordia diagram

Description

Plots U-Pb data on Wetherill, Tera-Wasserburg or U-Th-Pb concordia diagrams, calculates concordia_ages and compositions, evaluates the equivalence of multiple ($^{206}$Pb/$^{238}$U-$^{207}$Pb/$^{235}$U, $^{207}$Pb/$^{206}$Pb-$^{206}$Pb/$^{238}$U, or $^{208}$Th/$^{232}$Th-$^{206}$Pb/$^{238}$U) compositions, computes the weighted mean isotopic composition and the corresponding concordia_age using the method of maximum likelihood, computes the MSWD of equivalence and concordance and their respective Chi-squared p-values. Performs linear regression and computes the upper and lower intercept ages (for Wetherill) or the lower intercept age and the $^{207}$Pb/$^{206}$Pb intercept (for Tera-Wasserburg), taking into account error correlations and decay constant uncertainties.

Usage

concordia(
  x = NULL,
  tlim = NULL,
  type = 1,
  show.numbers = FALSE,
  levels = NA,
  clabel = "",
  ellipse.fill = c("#00FF0080", "#FF000080"),
  ellipse.stroke = "black",
  concordia.col = "darksalmon",
  exterr = FALSE,
  show.age = 0,
  oerr = 3,
  sigdig = 2,
  common.Pb = 0,
  ticks = 5,
  anchor = 0,
  hide = NULL,
  omit = NULL,
  omit.fill = NA,
  omit.stroke = "grey",
  ...
)

Arguments

x	an object of class UPb
tlim	age limits of the concordia line
type	one of 1: Wetherill – ${}^{206}$Pb/${}^{238}$U vs. ${}^{207}$Pb/${}^{235}$U 2: Tera-Wasserburg – ${}^{207}$Pb/${}^{206}$Pb vs. ${}^{238}$U/${}^{206}$Pb 3: U-Th-Pb concordia – ${}^{208}$Pb/${}^{232}$Th vs. ${}^{206}$Pb/${}^{238}$U (only available if x$format=7 or 8)
show.numbers	logical flag (TRUE to show grain numbers)
levels	a vector with length(x) values to be displayed as different background colours within the error ellipses.
clabel	label for the colour legend (only used if levels is not NA).
ellipse.fill	Fill colour for the error ellipses. This can either be a single colour or multiple colours to form a colour ramp. Examples: 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 ellipses, set ellipse.fill=NA
ellipse.stroke	the stroke colour for the error ellipses. Follows the same formatting guidelines as ellipse.fill
concordia.col	colour of the concordia line
exterr	show decay constant uncertainties?
show.age	one of either: 0: plot the data without calculating an age 1: fit a concordia composition and age 2: fit a discordia line through the data using the maximum likelihood algorithm of Ludwig (1998), which assumes that the scatter of the data is solely due to the analytical uncertainties. In this case, IsoplotR will either calculate an upper and lower intercept age (for Wetherill concordia), or a lower intercept age and common ($^{207}$Pb/$^{206}$Pb)-ratio intercept (for Tera-Wasserburg). If mswd>0, then the analytical uncertainties are augmented by a factor $\sqrt{mswd}$. 3: fit a discordia line ignoring the analytical uncertainties 4: fit a discordia line using a modified maximum likelihood algorithm that includes accounts for any overdispersion by adding a geological (co)variance term.
oerr	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').
sigdig	number of significant digits for the concordia/discordia age
common.Pb	common lead projection: 0:none 1: use the Pb-composition stored in settings('iratio','Pb207Pb206') (if x$format<4); settings('iratio','Pb206Pb204') and settings('iratio','Pb207Pb204') (if 3<x$format<7); or settings('iratio','Pb208Pb206') and settings('iratio','Pb208Pb207') (if x$format>6). 2: use the isochron intercept as the initial Pb-composition. If show.age>1, the data are projected along the isochron line, but the isochron itself is based on the uncorrected data. 3: use the Stacey-Kramers two-stage model to infer the initial Pb-composition.
ticks	either a scalar indicating the desired number of age ticks to be placed along the concordia line, OR a vector of tick ages.
anchor	control parameters to fix the intercept age or common Pb composition of the isochron fit. This can be a scalar or a vector. If anchor[1]=0: do not anchor the isochron. If anchor[1]=1: fix the common Pb composition at the values stored in settings('iratio',...). If anchor[1]=2: force the isochron line to intersect the concordia line at an age equal to anchor[2]. If anchor[1]=3: anchor the non-radiogenic component to the Stacey-Kramers mantle evolution model.
hide	vector with indices of aliquots that should be removed from the concordia diagram
omit	vector with indices of aliquots that should be plotted but omitted from concordia or discordia age calculation
omit.fill	fill colour that should be used for the omitted aliquots.
omit.stroke	stroke colour that should be used for the omitted aliquots.
...	optional arguments passed on to scatterplot

Details

The concordia diagram is a graphical means of assessing the internal consistency of U-Pb data. It sets out the measured $^{206}$Pb/$^{238}$U- and $^{207}$Pb/$^{235}$U-ratios against each other (‘Wetherill’ diagram); or, equivalently, the $^{207}$Pb/$^{206}$Pb- and $^{206}$Pb/$^{238}$U-ratios (‘Tera-Wasserburg’ diagram). Alternatively, for data format 7 and 8, it is also possible to plot $^{208}$Pb/$^{232}$Th against the $^{206}$Pb/$^{238}$U. The space of concordant isotopic compositions is marked by a curve, the ‘concordia line’. Isotopic ratio measurements are shown as 100(1-alpha)% confidence ellipses. Concordant samples plot near to, or overlap with, the concordia line. They represent the pinnacle of geochronological robustness. Samples that plot away from the concordia line but are aligned along a linear trend form an isochron (or ‘discordia’ line) that can be used to infer the composition of the non-radiogenic (‘common’) lead or to constrain the timing of prior lead loss.

Value

If show.age=1, returns a list with the following items:

x

a named vector with the (weighted mean) U-Pb composition

cov

the covariance matrix of the (weighted mean) U-Pb composition

mswd

a vector with three items (equivalence, concordance and combined) containing the MSWD (Mean of the Squared Weighted Deviates, a.k.a the reduced Chi-squared statistic) of isotopic equivalence, age concordance and combined goodness of fit, respectively.

p.value

a vector with three items (equivalence, concordance and combined) containing the p-value of the Chi-square test for isotopic equivalence, age concordance and combined goodness of fit, respectively.

df

a three-element vector with the number of degrees of freedom used for the mswd calculation.

age

a two-or three-element vector with:
t: the concordia_age (in Ma)
s[t]: the standard error of t
disp[t]: the standard error of t augmented by $\sqrt{mswd}$ to account for any overdispersion.

If show.age=2, 3 or 4, returns a list with the following items:

model

the fitting model (=show.age-1).

par

a vector with the upper and lower intercept ages (if type=1) or the lower intercept age and common Pb intercept(s) (if type=2). If show.age=4, includes an overdispersion term as well.

cov

the covariance matrix of the elements in par.

logpar

the logarithm of par

logcov

the logarithm of cov

err

a matrix with on or two rows:

s: the standard errors of the parameter estimates

disp: the standard errors of the parameter estimates augmented by $\sqrt{mswd}$ to account for overdispersed datasets (only reported if show.age=2).

df

the degrees of freedom of the concordia fit (concordance + equivalence)

p.value

p-value of a Chi-square test for age homogeneity (only reported if type=3).

mswd

mean square of the weighted deviates – a goodness-of-fit measure. mswd > 1 indicates overdispersion w.r.t the analytical uncertainties (not reported if show.age=3).

n

the number of aliquots in the dataset

References

Ludwig, K.R., 1998. On the treatment of concordant uranium-lead ages. Geochimica et Cosmochimica Acta, 62(4), pp.665-676.

Examples

attach(examples)
concordia(UPb,show.age=2)

dev.new()
concordia(UPb,type=2,xlim=c(24.9,25.4),
          ylim=c(0.0508,0.0518),ticks=249:254,exterr=TRUE)

dev.new()
concordia(UPb,show.age=2,anchor=c(2,260))

---

Prepare geochronological data for York regression

Description

Takes geochronology data as input and produces a five-column table with the variables, their uncertainties and error correlations as output. These can subsequently be used for York regression.

Usage

data2york(x, ...)

## Default S3 method:
data2york(x, format = 1, ...)

## S3 method for class 'other'
data2york(x, ...)

## S3 method for class 'UPb'
data2york(x, option = 1, tt = 0, ...)

## S3 method for class 'ArAr'
data2york(x, inverse = TRUE, ...)

## S3 method for class 'ThPb'
data2york(x, inverse = FALSE, ...)

## S3 method for class 'KCa'
data2york(x, inverse = FALSE, ...)

## S3 method for class 'PbPb'
data2york(x, inverse = TRUE, ...)

## S3 method for class 'PD'
data2york(x, exterr = FALSE, inverse = FALSE, ...)

## S3 method for class 'UThHe'
data2york(x, ...)

## S3 method for class 'ThU'
data2york(x, type = 2, generic = TRUE, ...)

Arguments

x	a five or six column matrix OR an object of class UPb, PbPb, ThPb, ArAr, ThU, UThHe, or PD (which includes objects of class RbSr, SmNd, LuHf and ReOs), generated by the read.data(...) function
...	optional arguments
format	one of 1 or 2: X, s[X], Y, s[Y], rXY; where rXY is the error correlation between X and Y; or 3: X/Z, s[X/Z], Y/Z, s[Y/Z], X/Y, s[X/Y]; for which the error correlations are automatically computed from the redundancy of the three ratios.
option	one of 1: Wetherill concordia ratios: X=07/35, sX=s[07/35], Y=06/38, sY=s[06/38], rXY. 2: Tera-Wasserburg ratios: X=38/06, sX=s[38/06], Y=07/06, sY=s[07/06], rXY. 3: X=38/06, sX=s[38/06], Y=04/06, sY=s[04/06], rXY (only valid if x$format=4,5 or 6). 4: X=35/07, sX=s[35/07], Y=04/07, sY=s[04/07], rXY (only valid if x$format=4,5 or 6). 5: U-Th-Pb concordia ratios: X=06/38, sX=s[06/38], Y=08/32, sY=s[08/32], rXY (only valid if x$format=7,8). 6: X=38/06, sX=s[38/06], Y=08/06, sY=s[08/06], rXY (only valid if x$format=7,8). 7: X=35/07, sX=s[35/07], Y=08/07, sY=s[08/07], rXY (only valid if x$format=7,8). 8: X=32/08, sX=s[32/08], Y=06/08, sY=s[06/08], rXY (only valid if x$format=7,8). 9: X=32/08, sX=s[32/08], Y=07/08, sY=s[07/08], rXY (only valid if x$format=7,8).
tt	the age of the sample. This is only used if x$format=7 or 8, in order to calculate the inherited ${}^{208}$Pb/${}^{232}$Th ratio.
inverse	toggles between normal and inverse isochron ratios. data2york returns five columns X, s[X], Y, s[Y] and r[X,Y]. If inverse=TRUE, then X = ${}^{204}$Pb/${}^{206}$Pb and Y = ${}^{207}$Pb/${}^{206}$Pb (if x has class PbPb), or X = ${}^{232}$Th/${}^{208}$Pb and Y = ${}^{204}$Pb/${}^{208}$Pb (if x has class ThPb), or X = ${}^{39}$Ar/${}^{40}$Ar and Y = ${}^{36}$Ar/${}^{40}$Ar (if x has class ArAr), or X = ${}^{40}$K/${}^{40}$Ca and Y = ${}^{44}$Ca/${}^{40}$Ca (if x has class KCa), or X = ${}^{87}$Rb/${}^{87}$Sr and Y = ${}^{86}$Sr/${}^{87}$Sr (if x has class RbSr), or X = ${}^{147}$Sm/${}^{143}$Nd and Y = ${}^{144}$Nd/${}^{143}$Nd (if x has class SmNd), or X = ${}^{187}$Re/${}^{187}$Os and Y = ${}^{188}$Os/${}^{187}$Os (if x has class ReOs), or X = ${}^{176}$Lu/${}^{176}$Hf and Y = ${}^{177}$Hf/${}^{176}$Hf (if x has class LuHf). If inverse=FALSE, then X = ${}^{206}$Pb/${}^{204}$Pb and Y = ${}^{207}$Pb/${}^{204}$Pb (if x has class PbPb), or X = ${}^{232}$Th/${}^{204}$Pb and Y = ${}^{208}$Pb/${}^{204}$Pb (if x has class ThPb), or X = ${}^{39}$Ar/${}^{36}$Ar and Y = ${}^{40}$Ar/${}^{36}$Ar (if x has class ArAr), or X = ${}^{40}$K/${}^{44}$Ca and Y = ${}^{40}$Ca/${}^{44}$Ca (if x has class KCa), or X = ${}^{87}$Rb/${}^{86}$Sr and Y = ${}^{87}$Sr/${}^{86}$Sr (if x has class RbSr), or X = ${}^{147}$Sm/${}^{144}$Nd and Y = ${}^{143}$Nd/${}^{144}$Nd (if x has class SmNd), or X = ${}^{187}$Re/${}^{188}$Os and Y = ${}^{187}$Os/${}^{188}$Os (if x has class ReOs), or X = ${}^{176}$Lu/${}^{177}$Hf and Y = ${}^{176}$Hf/${}^{177}$Hf (if x has class LuHf).
exterr	If TRUE, propagates the external uncertainties (e.g. decay constants) into the output errors.
type	Return ‘Rosholt’ or ‘Osmond’ ratios? Rosholt (type=1) returns X=8/2, sX=s[8/2], Y=0/2, sY=s[0/2], rXY. Osmond (type=2) returns X=2/8, sX=s[2/8], Y=0/8, sY=s[0/8], rXY.
generic	If TRUE, uses the following column headers: X, sX, Y, sY, rXY. If FALSE and type=1, uses U238Th232, errU238Th232, Th230Th232, errTh230Th232, rXY or if FALSE and type=2, uses Th232U238, errTh232U238, Th230U238, errTh230U238, rXY.

Value

a five-column table that can be used as input for york-regression.

See Also

york

Examples

f <- system.file("RbSr1.csv",package="IsoplotR")
dat <- read.csv(f)
yorkdat <- data2york(dat)
fit <- york(yorkdat)

---

Set up a discordance filter

Description

Define a discordance cutoff to filter U–Pb data.

Usage

discfilter(option = 0, before = TRUE, cutoff)

Arguments

option	one of five options: 0: do not apply a discordance filter 1 or 't': the absolute age difference (Ma) between the $^{206}$Pb/$^{238}$U and $^{207}$Pb/$^{206}$Pb ages. 2 or 'r': the relative age difference (%) between the $^{206}$Pb/$^{238}$U and $^{207}$Pb/$^{206}$Pb ages. 3 or 'sk': percentage of common Pb measured along a mixing line connecting the measured composition and the Stacey-Kramers mantle composition in Tera-Wasserburg space. 4 or 'a': logratio distance (%) measured along a perpendicular line connecting Tera-Wasserburg concordia and the measured composition. 5 or 'c': logratio distance (%) measured along a line connecting the measured composition and the corresponding single grain concordia_age composition. Further details in Vermeesch (2021).
before	logical flag indicating whether the discordance filter should be applied before (TRUE) or after (FALSE) the common-Pb correction.
cutoff	a two-element vector with the minimum (negative) and maximum (positive) allowed discordance. Default values vary between the different options. To view them, enter discfilter(option) at the R command line.

Details

The most reliable U–Pb age constraints are obtained from (zircon) grains whose $^{206}$Pb/$^{238}$U and $^{207}$Pb/$^{206}$Pb ages are statistically indistinguishable from each other. U–Pb compositions that fulfil this requirements are called ‘concordant’. Those that violate it are called ‘discordant’. The discordance of the $^{206}$Pb/$^{238}$U and $^{207}$Pb/$^{206}$Pb systems can be defined in five different ways. By setting a cutoff for any of these criteria, U–Pb data can be filtered for data quality.

Value

a list with the input parameters. Default values for cutoff are

c(-48,140) if option=='t';

c(-5,15) if option=='r';

c(-0.36,0.96) if option=='sk';

c(-1.6,4.7) if option=='a'; and

c(-2,5.8) if option=='c'.

References

Vermeesch (2021) “On the treatment of discordant data in detrital zircon U–Pb geochronology”, Geochronology.

See Also

cad, kde, radialplot

Examples

dscf <- discfilter(option='c',before=TRUE,cutoff=c(-1,1))
weightedmean(x=examples$UPb,exterr=FALSE,sigdig=2,
             cutoff.disc=dscf,common.Pb=3)

---

Set up U-series disequilibrium correction for U-Pb geochronology

Description

The U-Pb method conventionally assumes initial secular equilibrium of all the intermediate daughters of the ${}^{238}$U-${}^{206}$Pb and ${}^{235}$U-${}^{207}$Pb decay chains. Violation of this assumption may produce inaccurate results. diseq sets up initial disequilibrium parameters that are subsequently passed on to the read.data function for incorporation in other functions.

Usage

diseq(
  U48 = list(x = 1, sx = 0, option = 0, m = 0, M = 20, x0 = 1, sd = 10),
  ThU = list(x = 1, sx = 0, option = 0, m = 0, M = 20, x0 = 1, sd = 10),
  RaU = list(x = 1, sx = 0, option = 0, m = 0, M = 20, x0 = 1, sd = 10),
  PaU = list(x = 1, sx = 0, option = 0, m = 0, M = 20, x0 = 1, sd = 10),
  buffer = 1e-05
)

Arguments

U48	a list containing seven items (x, sx, m, M, x0, sd and option) specifying the ${}^{234}$U/${}^{238}$U disequilibrium. If option=0, then x and sx are ignored and no disequilibrium correction is applied. If option=1, then x contains the initial ${}^{234}$U/${}^{238}$U ratio and sx its standard error. If option=2, then x contains the measured ${}^{234}$U/${}^{238}$U ratio and sx its standard error. m, M specify the minimum and maximum search limits of the initial ${}^{234}$U/${}^{238}$U activity ratio. x0 and sd specify the mean and standard deviation of the prior distribution for the the initial ${}^{234}$U/${}^{238}$U activity ratio.
ThU	a list containing seven items (x, sx, m, M, x0, sd and option) specifying the ${}^{230}$Th/${}^{238}$U disequilibrium. If option=0, then x and sx are ignored and no disequilibrium correction is applied. If option=1, then x contains the initial ${}^{230}$Th/${}^{238}$U ratio and sx its standard error. If option=2, then x contains the measured ${}^{230}$Th/${}^{238}$U ratio and sx its standard error. If option=3, then x contains the measured Th/U ratio of the magma (assumed or determined from the whole rock or volcanic glass) and sx its standard error. This only applies for Th-bearing U-Pb data formats 7 and 8. m, M, x0 and sd are analogous to the eponymous settings for ThU.
RaU	a list containing seven items (x, sx, m, M, x0, sd and option) specifying the ${}^{226}$Ra/${}^{238}$U disequilibrium. If option=0, then x and sx are ignored and no disequilibrium correction is applied. If option=1, then x contains the initial ${}^{226}$Ra/${}^{238}$U ratio and sx its standard error. m, M, x0 and sd are analogous to the eponymous settings for ThU.
PaU	a list containing seven items (x, sx, m, M, x0, sd and option) specifying the ${}^{231}$Pa/${}^{235}$U disequilibrium. If option=0, then x and sx are ignored and no disequilibrium correction is applied. If option=1, then x contains the initial ${}^{231}$Pa/${}^{235}$U ratio and sx its standard error. m, M, x0 and sd are analogous to the eponymous settings for ThU.
buffer	small amount of padding to avoid singularities in the prior distribution, which uses a logistic transformation: $y = \ln\left(\frac{x-m+buffer}{M+buffer-x}\right)$

Details

There are three ways to correct for the initial disequilibrium between the activity of ${}^{238}$U, ${}^{234}$U, ${}^{230}$Th, and ${}^{226}$Ra; or between ${}^{235}$U and ${}^{231}$Pa:

1. Specify the assumed initial activity ratios and calculate how much excess ${}^{206}$Pb and ${}^{207}$Pb these would have produced.

2. Measure the current activity ratios to infer the initial ratios. This approach only works for young samples.

3. The initial ${}^{230}$Th/${}^{238}$U activity ratio can also be estimated by providing the Th/U-ratio of the magma.

Value

a list with the following items:

U48, ThU, RaU, PaU

the same as the corresponding input arguments

equilibrium

a boolean flag indicating whether option=TRUE and/or x=1 for all activity ratios

Q

the eigenvectors of the disequilibrium matrix exponential

Qinv

the inverse of Q

L

a named vector of all the relevant decay constants

See Also

mclean, concordia, ludwig

Examples

d <- diseq(U48=list(x=0,option=1),ThU=list(x=2,option=1),
           RaU=list(x=2,option=1),PaU=list(x=2,option=1))
fn <- system.file("diseq.csv",package="IsoplotR")
UPb <- read.data(fn,method='U-Pb',format=2,d=d)
concordia(UPb,type=2,xlim=c(0,700),ylim=c(0.05,0.5))

---

Dissimilarity between detrital age distributions

Description

Calculates the pairwise dissimilarity between detrital age distributions, using either the Wasserstein-2 or Kolmogorov-Smirnov distance.

Usage

diss(x, ...)

## Default S3 method:
diss(x, y, method = "KS", ...)

## S3 method for class 'detritals'
diss(x, method = "W2", ...)

Arguments

x	an object of class detrital OR a vector of numbers
...	extra arguments (not used)
y	a vector of numbers
method	either 'KS' (for Kolmogorov-Smirnov distance), or 'W2' (for Wasserstein-2 distance).

Details

The Kolmogorov-Smirnov statistic is the maximum vertical difference between two empirical cumulative distribution functions. The Wasserstein distance is a function of the area between them. Both dissimilarity measures are useful for multidimensional scaling.

Value

an object of class dist.

Author(s)

Written by Pieter Vermeesch, using modified code from Mathieu Vrac's CDFt package (KolmogorovSmirnov function), and Dominic Schuhmacher's transport package (transport1d function).

See Also

mds

Examples

d <- diss(examples$DZ,method='KS')
mds(d)

---

Get error ellipse coordinates for plotting

Description

Constructs an error ellipse at a given confidence level from its centre and covariance matrix

Usage

ellipse(x, y, covmat, alpha = 0.05, n = 50)

Arguments

x	x-coordinate (scalar) for the centre of the ellipse
y	y-coordinate (scalar) for the centre of the ellipse
covmat	the [2x2] covariance matrix of the x-y coordinates
alpha	the probability cutoff for the error ellipses
n	the resolution (number of segments) of the error ellipses

Value

an [nx2] matrix of plot coordinates

Examples

x = 99; y = 101;
covmat <- matrix(c(1,0.9,0.9,1),nrow=2)
ell <- ellipse(x,y,covmat)
plot(c(90,110),c(90,110),type='l')
polygon(ell,col=rgb(0,1,0,0.5))
points(x,y,pch=21,bg='black')

---

Th-U evolution diagram

Description

Plots Th-U data on a $^{234}$U/$^{238}$U-$^{230}$Th/$^{238}$U evolution diagram, a $^{234}$U/$^{238}$U-age diagram, or (if $^{234}$U/$^{238}$U is assumed to be in secular equilibrium), a $^{230}$Th/$^{232}$Th-$^{238}$U/$^{232}$Th diagram; calculates isochron ages.

Usage

evolution(
  x,
  xlim = NULL,
  ylim = NULL,
  tticks = NULL,
  aticks = NULL,
  oerr = 3,
  transform = FALSE,
  Th0i = 0,
  show.numbers = FALSE,
  levels = NA,
  clabel = "",
  ellipse.fill = c("#00FF0080", "#FF000080"),
  ellipse.stroke = "black",
  line.col = "darksalmon",
  isochron = FALSE,
  model = 1,
  exterr = FALSE,
  sigdig = 2,
  hide = NULL,
  omit = NULL,
  omit.fill = NA,
  omit.stroke = "grey",
  ...
)

Arguments

x	an object of class ThU
xlim	x-axis limits
ylim	y-axis limits
tticks	time intervals of the evolution grid
aticks	initial activity ratio ticks of the evolution grid
oerr	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').
transform	if TRUE, plots $^{234}$U/$^{238}$U vs. Th-U age.
Th0i	initial $^{230}$Th correction. 0: no correction 1: if x$format is 1 or 2, project the data along an isochron fit. If x$format is 3 or 4, infer the initial $^{230}$Th/$^{238}$U activity ratio from the isochron. 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, anchor the isochrons to the equiline, based on 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).
show.numbers	label the error ellipses with the grain numbers?
levels	a vector with additional values to be displayed as different background colours within the error ellipses.
clabel	label of the colour legend.
ellipse.fill	fill colour for the error ellipses. This can either be a single colour or multiple colours to form a colour ramp. Examples: 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 ellipses, set ellipse.fill=NA
ellipse.stroke	the stroke colour for the error ellipses. Follows the same formatting guidelines as ellipse.fill
line.col	colour of the age grid
isochron	fit an isochron to the data?
model	if isochron=TRUE, choose one of three regression models: 1: maximum likelihood regression, using either the modified error weighted least squares algorithm of York et al. (2004) for 2-dimensional data, or the Maximum Likelihood formulation of Ludwig and Titterington (1994) for 3-dimensional data. These algorithms take into account the analytical uncertainties and error correlations, under the assumption that the scatter between the data points is solely caused by the analytical uncertainty. If this assumption is correct, then the MSWD value should be approximately equal to one. There are three strategies to deal with the case where MSWD>1. The first of these is to assume that the analytical uncertainties have been underestipmated by a factor $\sqrt{MSWD}$. 2: total least squares regression: a second way to deal with over- or underdispersed datasets is to simply ignore the analytical uncertainties. 3: maximum likelihood regression with overdispersion: instead of attributing any overdispersion (MSWD > 1) to underestimated analytical uncertainties (model 1), one can also attribute it to the presence of geological uncertainty, which manifests itself as an added (co)variance term.
exterr	propagate the decay constant uncertainty in the isochron age?
sigdig	number of significant digits for the isochron age
hide	vector with indices of aliquots that should be removed from the plot.
omit	vector with indices of aliquots that should be plotted but omitted from the isochron age calculation.
omit.fill	fill colour that should be used for the omitted aliquots.
omit.stroke	stroke colour that should be used for the omitted aliquots.
...	optional arguments to the generic plot function

Details

Similar to the concordia diagram (for U-Pb data) and the helioplot diagram (for U-Th-He data), the evolution diagram simultaneously displays the isotopic composition and age of U-series data. For carbonate data (Th-U formats 1 and 2), the Th-U evolution diagram consists of a scatter plot that sets out the $^{234}$U/$^{238}$U-activity ratios against the $^{230}$Th/$^{238}$U-activity ratios as error ellipses, and displays the initial $^{234}$U/$^{238}$U-activity ratios and ages as a set of intersecting lines. Alternatively, the $^{234}$U/$^{238}$U-ratios can also be set out against the $^{230}$Th-$^{234}$U-$^{238}$U-ages. In both types of evolution diagrams, IsoplotR provides the option to project the raw measurements along the best fitting isochron line and thereby remove the detrital $^{230}$Th-component. This procedure allows a visual assessment of the degree of homogeneity within a dataset, as is quantified by the MSWD.

Neither the U-series evolution diagram, nor the $^{234}$U/$^{238}$U vs. age plot is applicable to igneous datasets (Th-U formats 3 and 4), in which $^{234}$U and $^{238}$U are in secular equilibrium. For such datasets, IsoplotR produces an Osmond-style regression plot that is decorated with a fanning set of isochron lines.

References

Ludwig, K.R. and Titterington, D.M., 1994. Calculation of $^{230}$Th/U isochrons, ages, and errors. Geochimica et Cosmochimica Acta, 58(22), pp.5031-5042.

Ludwig, K.R., 2003. Mathematical-statistical treatment of data and errors for $^{230}$Th/U geochronology. Reviews in Mineralogy and Geochemistry, 52(1), pp.631-656.

See Also

isochron

Examples

attach(examples)
evolution(ThU)

dev.new()
evolution(ThU,transform=TRUE,isochron=TRUE,model=1)

---