Isoplot R Manual

User Manual:

Open the PDF directly: View PDF .
Page Count: 63

age
agespectrum
cad
central
concordia
data2york
ellipse
evolution
examples
helioplot
isochron
IsoplotR
kde
ludwig
mds
peakfit
radialplot
read.data
set.zeta
settings
titterington
weightedmean
york
Index

Package ‘IsoplotR’

August 29, 2018

Title Statistical Toolbox for Radiometric Geochronology

Version 1.5

Date 2018-08-21

Description Plots U-Pb data on Wetherill and Tera-Wasserburg concordia diagrams. Calculates con-

cordia and discordia ages. Performs linear regression of measurements with correlated errors us-

ing 'York', 'Titterington' and 'Ludwig' approaches. Generates Kernel Density Esti-

mates (KDEs) and Cumulative Age Distributions (CADs). Produces Multidimensional Scal-

ing (MDS) conﬁgurations and Shepard plots of multi-sample detrital datasets using the Kol-

mogorov-Smirnov distance as a dissimilarity measure. Calcu-

lates 40Ar/39Ar ages, isochrons, and age spectra. Computes weighted means account-

ing for overdispersion. Calculates U-Th-He (single grain and central) ages, logra-

tio plots and ternary diagrams. Processes ﬁssion track data using the external detec-

tor method and LA-ICP-MS, calculates central ages and plots ﬁssion track and other data on ra-

dial (a.k.a. 'Galbraith') plots. Constructs total Pb-U, Pb-Pb, Re-Os, Sm-Nd, Lu-Hf, Rb-

Sr and 230Th-U isochrons as well as 230Th-U evolution plots.

Author Pieter Vermeesch [aut, cre]

Maintainer Pieter Vermeesch <p.vermeesch@ucl.ac.uk>

Depends R (>= 3.0.0)

Imports MASS, grDevices, graphics, stats, utils

License GPL-3

URL http://isoplotr.london-geochron.com

LazyData true

RoxygenNote 6.0.1

Rtopics documented:

age.............................................. 2

agespectrum......................................... 6

cad.............................................. 8

central............................................ 11

concordia .......................................... 13

data2york .......................................... 16

2age

ellipse............................................ 17

evolution .......................................... 18

examples .......................................... 21

helioplot........................................... 23

isochron........................................... 25

IsoplotR........................................... 31

kde.............................................. 31

ludwig............................................ 36

mds ............................................. 37

peakﬁt............................................ 39

radialplot .......................................... 43

read.data........................................... 47

set.zeta ........................................... 50

settings ........................................... 52

titterington.......................................... 55

weightedmean........................................ 56

york ............................................. 60

Index 63

age Calculate isotopic ages

Description

Calculates U-Pb, Pb-Pb, Ar-Ar, Re-Os, Sm-Nd, Rb-Sr, Lu-Hf, U-Th-He, Th-U and ﬁssion 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", exterr = TRUE, J = c(NA,

NA), zeta = c(NA, NA), rhoD = c(NA, NA), ...)

## S3 method for class 'UPb'

age(x, type = 1, wetherill = TRUE, exterr = TRUE, i = NA,

sigdig = NA, common.Pb = 0, ...)

## S3 method for class 'PbPb'

age(x, isochron = TRUE, common.Pb = 1, exterr = TRUE,

i = NA, sigdig = NA, ...)

## S3 method for class 'ArAr'

age(x, isochron = FALSE, i2i = TRUE, exterr = TRUE,

i = NA, sigdig = NA, ...)

age 3

## S3 method for class 'UThHe'

age(x, isochron = FALSE, central = FALSE, i = NA,

sigdig = NA, ...)

## S3 method for class 'fissiontracks'

age(x, central = FALSE, i = NA, sigdig = NA,

exterr = TRUE, ...)

## S3 method for class 'ThU'

age(x, isochron = FALSE, i2i = TRUE, exterr = TRUE,

i = NA, sigdig = NA, detritus = 0, Th02 = c(0, 0), Th02U48 = c(0, 0,

1e+06, 0, 0, 0, 0, 0, 0), ...)

## S3 method for class 'ReOs'

age(x, isochron = TRUE, i2i = TRUE, exterr = TRUE,

i = NA, sigdig = NA, ...)

## S3 method for class 'SmNd'

age(x, isochron = TRUE, i2i = TRUE, exterr = TRUE,

i = NA, sigdig = NA, ...)

## S3 method for class 'RbSr'

age(x, isochron = TRUE, i2i = TRUE, exterr = TRUE,

i = NA, sigdig = NA, ...)

## S3 method for class 'LuHf'

age(x, isochron = TRUE, i2i = TRUE, exterr = TRUE,

i = NA, sigdig = NA, ...)

Arguments

xcan 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 (if method='fissiontracks'),

• a four element vector containing Ar40Ar39,s[Ar40Ar39],J,s[J],

• 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 six element vector containing Rb,s[Rb],Sr,s[Sr],Sr87Sr86, and s[Sr87Sr86]

• a six element vector containing Re,s[Re],Os,s[Os],Os187Os188, and

s[Os187Os188]

• a six element vector containing Sm,s[Sm],Nd,s[Nd],Nd143Nd144, and

s[Nd144Nd143]

• a six element vector containing Lu,s[Lu],Hf,s[Hf],Hf176Hf177, and

s[Hf176Hf177]

4age

• a ﬁve element vector containing 0/8,s[0/8],4/8,s[4/8] and cov[0/8,4/8]

• an object of class UPb,PbPb,ArAr,ThU,RbSr,SmNd,ReOs,LuHf,UThHe or

fissiontracks.

... additional arguments

method one of either 'U238-Pb206','U235-Pb207','Pb207-Pb206','Ar-Ar','Th-U',

'Re-Os','Sm-Nd','Rb-Sr','Lu-Hf','U-Th-He'or 'fissiontracks'

exterr propagate the external (decay constant and calibration factor) uncertainties?

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

type scalar ﬂag 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 ﬁtted through all the U-Pb analyses using the max-

imum 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 ﬁtted ignoring the analytical uncertainties.

5: a discordia line should be ﬁtted using a modiﬁed maximum likelihood al-

gorithm that accounts for overdispersion by adding a geological (co)variance

term.

wetherill logical ﬂag to indicate whether the data should be evaluated in Wetherill (TRUE)

or Tera-Wasserburg (FALSE) space. This option is only used when type=2

i(optional) index of a particular aliquot

sigdig number of signiﬁcant digits for the uncertainty estimate (only used if type=1,

isochron=FALSE and central=FALSE).

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

isochron logical ﬂag indicating whether each Ar-Ar analysis should be considered sepa-

rately (isochron=FALSE) or an isochron age should be calculated from all Ar-Ar

analyses together (isochron=TRUE).

i2i ‘isochron to intercept’: calculates the initial (aka ‘inherited’, ‘excess’, or ‘com-

mon’) 40Ar/36Ar, 207Pb/204Pb, 87Sr/86Sr, 143Nd/144Nd, 187Os/188Os or 176Hf/177Hf

ratio from an isochron ﬁt. 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.

central logical ﬂag indicating whether each analysis should be considered separately

(central=FALSE) or a central age should be calculated from all analyses to-

gether (central=TRUE).

age 5

detritus detrital 230Th correction (only applicable when x$format == 1 or 2.

0: no correction

1: project the data along an isochron ﬁt

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

3: correct the data using the measured present day 230Th/238U, 232Th/238U and

234U/238U-ratios in the detritus.

Th02 2-element vector with the assumed initial 230Th/232Th-ratio of the detritus and

its standard error. Only used if isochron==FALSE and 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

Value

1. if xis a scalar or a vector, returns the age using the geochronometer given by method and its

standard error.

2. if xhas 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.conc,err[t.conc],err[p.conc], containing the

207Pb/235U-age and standard error, the 206Pb/238U-age and standard error, the 207Pb/206Pb-

age and standard error, the single grain concordia age and standard error, and the p-value for

concordance, respectively.

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

4. if xhas class PbPb,ArAr,RbSr,SmNd,ReOs,LuHf,ThU or UThHe and isochron=FALSE, returns

a table of Pb-Pb, Ar-Ar, Rb-Sr, Sm-Nd, Re-Os, Lu-Hf, Th-U or U-Th-He ages and their

standard errors.

5. if xhas class ThU and isochron=FALSE, returns a 5-column table with the Th-U ages, their

standard errors, the initial 234U/238U-ratios, their standard errors, and the correlation coefﬁ-

cient between the ages and the initial ratios.

6. if xhas class PbPb,ArAr,RbSr,SmNd,ReOs,LuHf,UThHe or ThU and isochron=TRUE, returns

the output of the isochron function.

7. if xhas class fissiontracks and central=FALSE, returns a table of ﬁssion track ages and

standard errors.

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

function.

See Also

concordia,isochron,central

Examples

data(examples)

tUPb <- age(examples$UPb,type=1)

tconc <- age(examples$UPb,type=2)

tdisc <- age(examples$UPb,type=3)

tArAr <- age(examples$ArAr)

tiso <- age(examples$ArAr,isochron=TRUE,i2i=TRUE)

tcentral <- age(examples$FT1,central=TRUE)

6agespectrum

agespectrum Plot a (40Ar/39Ar) release spectrum

Description

Produces a plot of boxes whose widths correspond to the cumulative amount of 39Ar (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, alpha = 0.05, plateau = TRUE,

random.effects = TRUE, plateau.col = rgb(0, 1, 0, 0.5),

non.plateau.col = rgb(0, 1, 1, 0.5), sigdig = 2, line.col = "red",

lwd = 2, title = TRUE, show.ci = TRUE, xlab = "cumulative fraction",

ylab = "age [Ma]", ...)

## S3 method for class 'ArAr'

agespectrum(x, alpha = 0.05, plateau = TRUE,

random.effects = TRUE, plateau.col = rgb(0, 1, 0, 0.5),

non.plateau.col = rgb(0, 1, 1, 0.5), sigdig = 2, exterr = TRUE,

line.col = "red", lwd = 2, i2i = FALSE, ...)

Arguments

xa three-column matrix whose ﬁrst column gives the amount of 39Ar in each

aliquot, and whose second and third columns give the age and its uncertainty.

an object of class ArAr

... optional parameters to the generic plot function

alpha the conﬁdence level of the error bars/boxes and conﬁdence intervals.

plateau logical ﬂag indicating whether a plateau age should be calculated. If plateau=TRUE,

the function will compute the weighted mean of the largest succession of steps

that pass the Chi-square test for age homogeneity. If TRUE, 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.

plateau.col the ﬁll colour of the rectangles used to mark the steps belonging to the age

plateau.

agespectrum 7

non.plateau.col

if plateau=TRUE, the steps that do NOT belong to the plateau are given a differ-

ent colour.

sigdig the number of signiﬁcant digits of the numerical values reported in the title of

the graphical output (only used if plateau=TRUE).

line.col colour of the average age line

lwd width of the average age line

title add a title to the plot?

show.ci show a 100(1-α)% conﬁdence interval for the plateau age as a grey band

xlab x-axis label

ylab y-axis label

exterr propagate the external (decay constant and calibration factor) uncertainties?

i2i ‘isochron to intercept’: calculates the initial (aka ‘inherited’, ‘excess’, or ‘com-

mon’) 40Ar/36Ar ratio from an isochron ﬁt. Setting i2i to FALSE uses the default

values stored in settings('iratio',...)

Details

IsoplotR deﬁnes the ‘plateau age’ as the weighted mean age of the longest sequence (in terms

of cumulative 39Ar content) of consecutive heating steps that pass the modiﬁed Chauvenet crite-

rion (see weightedmean). Note that this deﬁnition is different (and simpler) than the one used by

Isoplot (Ludwig, 2003). However, it is important to mention that all deﬁnitions of an age plateau

are heuristic by nature and should not be used for quantitative inference.

Value

If plateau=TRUE, returns a list with the following items:

mean a 3-element vector with:

x: the plateau mean

s[x]: the estimated standard deviation of x

ci[x]: the width of a 100(1-α)% conﬁdence interval of t

disp a 3-element vector with:

w: the overdispersion, i.e. the standard deviation of the Normal distribution that is assumed to

describe the true ages.

ll: the width of the lower half of a 100(1-α)% conﬁdence interval for the overdispersion

ul: the width of the upper half of a 100(1-α)% conﬁdence interval for the overdispersion

df the degrees of freedom for the weighted mean plateau ﬁt

mswd the mean square of the weighted deviates of the plateau

p.value the p-value of a Chi-square test with df =n−2degrees of freedom, where nis the number

of steps in the plateau and 2 degrees of freedom have been removed to estimate the mean and

the dispersion.

fract the fraction of 39Ar contained in the plateau

plotpar plot parameters for the weighted mean (see weightedmean), which are not used in the age

spectrum

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

8cad

See Also

weightedmean

Examples

data(examples)

agespectrum(examples$ArAr,ylim=c(0,80))

cad 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]",

colmap = "heat.colors", col = "black", ...)

## S3 method for class 'detritals'

cad(x, pch = NA, verticals = TRUE, xlab = "age [Ma]",

colmap = "heat.colors", ...)

## S3 method for class 'UPb'

cad(x, pch = NA, verticals = TRUE, xlab = "age [Ma]",

col = "black", type = 4, cutoff.76 = 1100, cutoff.disc = c(-15, 5),

common.Pb = 0, ...)

## S3 method for class 'PbPb'

cad(x, pch = NA, verticals = TRUE, xlab = "age [Ma]",

col = "black", common.Pb = 1, ...)

## S3 method for class 'ArAr'

cad(x, pch = NA, verticals = TRUE, xlab = "age [Ma]",

col = "black", i2i = FALSE, ...)

## S3 method for class 'ThU'

cad(x, pch = NA, verticals = TRUE, xlab = "age [ka]",

col = "black", i2i = FALSE, detritus = 0, Th02 = c(0, 0),

Th02U48 = c(0, 0, 1e+06, 0, 0, 0, 0, 0, 0), ...)

## S3 method for class 'ReOs'

cad 9

cad(x, pch = NA, verticals = TRUE, xlab = "age [Ma]",

col = "black", i2i = TRUE, ...)

## S3 method for class 'SmNd'

cad(x, pch = NA, verticals = TRUE, xlab = "age [Ma]",

col = "black", i2i = TRUE, ...)

## S3 method for class 'RbSr'

cad(x, pch = NA, verticals = TRUE, xlab = "age [Ma]",

col = "black", i2i = TRUE, ...)

## S3 method for class 'LuHf'

cad(x, pch = NA, verticals = TRUE, xlab = "age [Ma]",

col = "black", i2i = TRUE, ...)

## S3 method for class 'UThHe'

cad(x, pch = NA, verticals = TRUE, xlab = "age [Ma]",

col = "black", ...)

## S3 method for class 'fissiontracks'

cad(x, pch = NA, verticals = TRUE,

xlab = "age [Ma]", col = "black", ...)

Arguments

xa numerical vector OR an object of class UPb,PbPb,ArAr,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 ﬂag indicating if the horizontal lines of the CAD should be connected by

vertical lines

xlab x-axis label

colmap an optional string with the name of one of R’s built-in colour palettes (e.g.,

heat.colors,terrain.colors,topo.colors,cm.colors), which are to be

used for plotting data of class detritals.

col colour to give to single sample datasets (not applicable if xhas class detritals)

type scalar indicating whether to plot the 207Pb/235U age (type=1), the 206Pb/238U

age (type=2), the 207Pb/206Pb age (type=3), the 207Pb/206Pb-206Pb/238U age

(type=4), or the (Wetherill) concordia age (type=5)

cutoff.76 the age (in Ma) below which the 206Pb/238U-age and above which the 207Pb/206Pb-

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 207Pb/235U and 206Pb/238U age (if 206Pb/238U < cutoff.76)

or between the 206Pb/238U and 207Pb/206Pb age (if 206Pb/238U > cutoff.76). Set

cutoff.disc=NA if you do not want to use this ﬁlter.

10 cad

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 ‘com-

mon’) 40Ar/36Ar, 207Pb/204Pb, 87Sr/86Sr, 143Nd/144Nd, 187Os/188Os, 230Th/232Th

or 176Hf/177Hf ratio from an isochron ﬁt. Setting i2i to FALSE uses the default

values stored in settings('iratio',...) or zero (for the Pb-Pb method).

When applied to data of class ThU, setting i2i to TRUE applies a detrital Th-

correction.

detritus detrital 230Th correction (only applicable when x$format == 1 or 2.

0: no correction

1: project the data along an isochron ﬁt

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

3: correct the data using the measured present day 230Th/238U, 232Th/238U and

234U/238U-ratios in the detritus.

Th02 2-element vector with the assumed initial 230Th/232Th-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

Empirical cumulative distribution functions or cumulative age distributions CADs are the most

straightforward way to visualise the probability distribution of multiple dates. Suppose that we

have a set of ndates ti. The the CAD is a step function that sets out the rank order of the dates

against their numerical value:

CAD(t) = Pi1(t < ti)/n

where 1(∗)=1if∗is true and 1(∗)=0if∗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 de-

trital apatite ﬁssion track thermochronology. Journal of Geophysical Research: Earth Surface,

112(F3).

See Also

kde,radialplot

central 11

Examples

data(examples)

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

central Calculate U-Th-He and ﬁssion track central ages and compositions

Description

Computes the geometric mean composition of a continuous mixture of ﬁssion track or U-Th-He

data and returns the corresponding age and ﬁtting parameters.

Usage

central(x, ...)

## Default S3 method:

central(x, alpha = 0.05, ...)

## S3 method for class 'UThHe'

central(x, alpha = 0.05, model = 1, ...)

## S3 method for class 'fissiontracks'

central(x, mineral = NA, alpha = 0.05, ...)

Arguments

xan object of class UThHe or fissiontracks, OR a 2-column matrix with (strictly

positive) values and uncertainties

... optional arguments

alpha cutoff value for conﬁdence intervals

model choose one of the following statistical models:

1: weighted mean. This model assumes that the scatter between the data points is

solely caused by the analytical uncertainty. If the 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 ﬁrst of these is to assume that the

analytical uncertainties have been underestimated by a factor √M SW D.

2: unweighted mean. A second way to deal with over- or underdispersed datasets

is to simply ignore the analytical uncertainties.

3: weighted mean with overdispersion: instead of attributing any overdispersion

(MSWD > 1) to underestimated analytical uncertainties (model 1), one could

also attribute it to the presence of geological uncertainty, which manifests itself

as an added (co)variance term.

mineral setting this parameter to either apatite or zircon changes the default efﬁciency

factor, initial ﬁssion track length and density to preset values (only affects results

if x$format=2)

12 central

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 ﬁssion track data, the analytical uncertainty is assumed to obey Poisson 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.

Value

If xhas class UThHe, 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 ﬁtting model.

df the degrees of freedom (2n−2) of the ﬁt (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 three- or four-element vector with:

t: the central age.

s[t]: the standard error of t.

ci[t]: the width of a 100(1 −α)% conﬁdence interval for t.

disp[t]: the studentised 100(1 −α)% conﬁdence interval enhanced by a factor

of √mswd (only reported if model=1).

wthe geological overdispersion term. If model=3, this is a three-element vector

with the standard deviation of the (assumedly) Normal dispersion and the lower

and upper half-widths of its 100(1 −α)% conﬁdence interval. w=0 if code-

model<3.

OR, otherwise:

age a three-element vector with:

t: the central age.

s[t]: the standard error of t.

ci[t]: the width of a 100(1 −α)% conﬁdence interval for t.

disp a three-element vector with the overdispersion (standard deviation) of the excess scatter, and

the upper and lower half-widths of its 100(1 −α)% conﬁdence interval.

concordia 13

mswd the reduced Chi-square statistic of data concordance, i.e. mswd =X2/df, where X2is 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 ﬁssion 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

data(examples)

print(central(examples$UThHe)$age)

concordia Concordia diagram

Description

Plots U-Pb data on Wetherill and Tera-Wasserburg concordia diagrams, calculate concordia ages

and compositions, evaluates the equivalence of multiple (206Pb/238U-207Pb/235U or 207Pb/206Pb-

206Pb/238U) compositions, computes the weighted mean isotopic composition and the correspond-

ing concordia age using the method of maximum likelihood, computes the MSWD of equiva-

lence 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

207Pb/206Pb intercept (for Tera-Wasserburg), taking into account error correlations and decay con-

stant uncertainties.

Usage

concordia(x, tlim = NULL, alpha = 0.05, wetherill = TRUE,

show.numbers = FALSE, levels = NA, clabel = clabel,

ellipse.col = c("#00FF0080", "#FF000080"), concordia.col = "darksalmon",

exterr = FALSE, show.age = 0, sigdig = 2, common.Pb = 0,

ticks = NULL, ...)

14 concordia

Arguments

xan object of class UPb

tlim age limits of the concordia line

alpha probability cutoff for the error ellipses and conﬁdence intervals

wetherill logical ﬂag (FALSE for Tera-Wasserburg)

show.numbers logical ﬂag (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.col a vector of two background colours for the error ellipses. If levels=NA, then

only the ﬁrst colour is used. If levels is a vector of numbers, then ellipse.col

is used to construct a colour ramp.

concordia.col colour of the concordia line

exterr show decay constant uncertainty?

show.age one of either:

0: plot the data without calculating an age

1: ﬁt a concordia composition and age

2: ﬁt 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 (207Pb/206Pb)-ratio intercept (for Tera-Wasserburg). If mswd>0, then

the analytical uncertainties are augmented by a factor √mswd.

3: ﬁt a discordia line ignoring the analytical uncertainties

4: ﬁt a discordia line using a modiﬁed maximum likelihood algorithm that in-

cludes accounts for any overdispersion by adding a geological (co)variance term.

sigdig number of signiﬁcant digits for the concordia/discordia age

common.Pb apply a common lead correction using one of three methods:

1: use the Stacey-Kramer two-stage model to infer the initial Pb-composition

2: use the isochron intercept as the initial Pb-composition

3: use the Pb-composition stored in settings('iratio','Pb206Pb204')and

settings('iratio','Pb207Pb204')

ticks an optional vector of age ticks to be added to the concordia line to override

IsoplotR’s default spacing, which is based on R’s pretty function.

... optional arguments to the generic plot function

Details

The concordia diagram is a graphical means of assessing the internal consistency of U-Pb data. It

sets out the measured 206Pb/238U- and 207Pb/235U-ratios against each other (‘Wetherill’ diagram)

or, equivalently, the 207Pb/206Pb- and 206Pb/238U-ratios (‘Tera-Wasserburg’ diagram). The space

of concordant isotopic compositions is marked by a curve, the ‘concordia line’. Isotopic ratio

measurements are shown as 100(1-alpha)% conﬁdence ellipses. Concordant samples plot near

concordia 15

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:

xa 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 ﬁt, 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 ﬁt, respectively.

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

These values are useful when expanding the analytical uncertainties if mswd>1.

age a 4-element vector with:

t: the concordia age (in Ma)

s[t]: the estimated uncertainty of t

ci[t]: the studentised 100(1 −α)% conﬁdence interval of tfor the appropriate degrees of

freedom

disp[t]: the studentised 100(1 −α)% conﬁdence interval for taugmented by √mswd to

account for overdispersed datasets.

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

model the ﬁtting model (=show.age-1).

xa two-element vector with the upper and lower intercept ages (if wetherill=TRUE) or the lower

intercept age and 207Pb/206Pb intercept (if wetherill=FALSE).

cov the covariance matrix of the elements in x.

err a[2 x 2] or [3 x 2] matrix with the following rows:

s: the estimated standard deviation for x

ci: the studentised 100(1 −α)% conﬁdence interval of xfor the appropriate degrees of free-

dom

disp[t]: the studentised 100(1 −α)% conﬁdence interval for xaugmented by √mswd to

account for overdispersed datasets (only reported if show.age=2).

df the degrees of freedom of the concordia ﬁt (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-ﬁt measure. mswd > 1 indicates

overdispersion w.r.t the analytical uncertainties (not reported if show.age=3).

wthree-element vector with the standard deviation of the (assumedly) Normal overdispersion term

and the lower and upper half-widths of its 100(1 −α)% conﬁdence interval (only important if

show.age=4).

nthe number of aliquots in the dataset

16 data2york

References

Ludwig, K.R., 1998. On the treatment of concordant uranium-lead ages. Geochimica et Cos-

mochimica Acta, 62(4), pp.665-676.

Examples

data(examples)

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

dev.new()

concordia(examples$UPb,wetherill=FALSE,

xlim=c(24.9,25.4),ylim=c(0.0508,0.0518),

ticks=249:254,exterr=TRUE)#'data(examples)

data2york Prepare geochronological data for York regression

Description

Takes geochronology data as input and produces a ﬁve-column table as output, which can be used

for York regression.

Usage

data2york(x, ...)

## Default S3 method:

data2york(x, format = 1, ...)

## S3 method for class 'UPb'

data2york(x, wetherill = TRUE, ...)

## S3 method for class 'ArAr'

data2york(x, inverse = TRUE, ...)

## S3 method for class 'PbPb'

data2york(x, inverse = TRUE, ...)

## S3 method for class 'PD'

data2york(x, exterr = FALSE, common = FALSE, ...)

## S3 method for class 'UThHe'

data2york(x, ...)

## S3 method for class 'ThU'

data2york(x, type = 2, generic = TRUE, ...)

ellipse 17

Arguments

xa ﬁve or six column matrix OR an object of class UPb,PbPb,ArAr,ThU,PD or

UThHe, generated by the read.data(...) function

... optional arguments

format one of

1. X,s[X],Y,s[Y],rho

where rho is the error correlation between Xand Y, or

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

wetherill If TRUE, returns a table with X=7/5,sX=s[7/5],Y=6/8,sY=s[6/8],rXY

If FALSE, returns a table with X=8/6,sX=s[8/6],Y=7/6,sY=s[7/6],rho=rXY

inverse If xhas class ArAr and inverse=FALSE, returns a table with columns X=9/6,

sX=s[9/6],Y=0/6, codesY=s[0/6], rXY

If xhas class ArAr and inverse=TRUE, returns a table with columns X=9/0,

sX=s[9/0],Y=6/0, codesY=s[6/0], rXY

If xhas class PbPb and inverse=FALSE, returns a table with columns X=6/4,

sX=s[6/4],Y=7/4, codesY=s[7/4], rXY

If xhas class PbPb and inverse=TRUE, returns a table with columns X=4/6,

sX=s[4/6],Y=7/6,sY=s[7/6],rXY

exterr If TRUE, propagates the external uncertainties (e.g. decay constants) into the

output errors.

common Applies a common lead correction to the data.

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,

rho

or if FALSE and type=2, uses Th232U238,errTh232U238,Th230U238,errTh230U238,

rho.

Value

f <- system.ﬁle("RbSr.csv",package="IsoplotR") dat <- read.csv(f) yorkdat <- data2york(dat) isochron(yorkdat)

ellipse Get coordinates of error ellipse for plotting

Description

Constructs an error ellipse at a given conﬁdence level from its centre and covariance matrix

18 evolution

Usage

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

Arguments

xx-coordinate (scalar) for the centre of the ellipse

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

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

evolution Th-U evolution diagram

Description

Plots Th-U data on a 234U/238U-230Th/238U evolution diagram, a 234U/238U-age diagram, or (if

234U/238U is assumed to be in secular equilibrium), a 230Th/232Th-238U/232Th diagram, calculates

isochron ages.

Usage

evolution(x, xlim = NA, ylim = NA, alpha = 0.05, transform = FALSE,

detritus = 0, Th02 = c(0, 0), Th02U48 = c(0, 0, 1e+06, 0, 0, 0, 0, 0,

0), show.numbers = FALSE, levels = NA, clabel = "",

ellipse.col = c("#00FF0080", "#FF000080"), line.col = "darksalmon",

isochron = FALSE, model = 1, exterr = TRUE, sigdig = 2, ...)

evolution 19

Arguments

xan object of class ThU

xlim x-axis limits

ylim y-axis limits

alpha probability cutoff for the error ellipses and conﬁdence intervals

transform if TRUE, plots 234U/238U vs. Th-U age.

detritus detrital 230Th correction (only applicable when x$format is 2or 3.

0: no correction

1: project the data along an isochron ﬁt

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

3: correct the data using the measured present day 230Th/238U, 232Th/238U and

234U/238U-ratios in the detritus.

Th02 2-element vector with the assumed initial 230Th/232Th-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

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.col a vector of two background colours for the error ellipses. If levels=NA, then

only the ﬁrst colour will be used. If levels is a vector of numbers, then

ellipse.col is used to construct a colour ramp.

line.col colour of the age grid

isochron ﬁt a 3D isochron to the data?

model if isochron=TRUE, choose one of three regression models:

1: maximum likelihood regression, using either the modiﬁed 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 uncertain-

ties 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 ﬁrst of these is to

assume that the analytical uncertainties have been underestipmated by a factor

√MSW D.

2: ordinary least squares regression: a second way to deal with over- or under-

dispersed 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 man-

ifests itself as an added (co)variance term.

20 evolution

exterr propagate the decay constant uncertainty in the isochron age?

sigdig number of signiﬁcant digits for the isochron age

... 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 234U/238U-activity ratios against the 230Th/238U-activity ratios as error ellipses, and

displays the initial 234U/238U-activity ratios and ages as a set of intersecting lines. Alternatively, the

234U/238U-ratios can also be set out against the 230Th-234U-238U-ages. In both types of evolution

diagrams, IsoplotR provides the option to project the raw measurements along the best ﬁtting

isochron line and thereby remove the detrital 230Th-component. This procedure allows a visual

assessment of the degree of homogeneity within a dataset, as is quantiﬁed by the MSWD.

Neither the U-series evolution diagram, nor the 234U/238U vs. age plot is applicable to igneous

datasets (Th-U formats 3 and 4), in which 234U and 238U 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 230Th/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 230Th/U geochronol-

ogy. Reviews in Mineralogy and Geochemistry, 52(1), pp.631-656.

See Also

isochron

Examples

data(examples)

evolution(examples$ThU)

dev.new()

evolution(examples$ThU,transform=TRUE,

isochron=TRUE,model=1)

examples 21

examples Example datasets for testing IsoplotR

Description

U-Pb, Pb-Pb, Ar-Ar, Re-Os, Sm-Nd, Rb-Sr, Lu-Hf, U-Th-He, Th-U, ﬁssion track and detrital

datasets

Details

examples an 18-item list containing:

UPb: an object of class UPb containing a high precision U-Pb dataset of Kamo et al. (1996) packaged

with Ken Ludwig (2003)’s Isoplot program.

PbPb: an object of class PbPb containing a Pb-Pb dataset from Connelley et al. (2017).

DZ: an object of class detrital containing a detrital zircon U-Pb dataset from Namibia (Vermeesch

et al., 2015).

ArAr: an object of class ArAr containing a 40Ar/39Ar spectrum of Skye basalt produced by Sarah

Sherlock (Open University).

UThHe: an object of class UThHe containing a U-Th-Sm-He dataset of Fish Lake apatite produced

by Daniel Stockli (UT Austin).

FT1: an object of class fissiontracks containing a synthetic external detector dataset.

FT2: an object of class fissiontracks containing a synthetic LA-ICP-MS-based ﬁssion track

dataset using the zeta calibration method.

FT3: an object of class fissiontracks containing a synthetic LA-ICP-MS-based ﬁssion track

dataset using the absolute dating approach.

ReOs: an object of class ReOs containing a 187Os/187Re-dataset from Selby (2007).

SmNd: an object of class SmNd containing a 143Nd/147Sm-dataset from Lugmair et al. (1975).

RbSr: an object of class RbSr containing an 87Rb/86Sr-dataset from Compston et al. (1971).

LuHf: an object of class LuHf containing an 176Lu/177Hf-dataset from Barfod et al. (2002).

ThU: an object of class ThU containing a synthetic ‘Osmond-type’ dataset from Titterington and

Ludwig (1994).

LudwigMean: an object of class other containing a collection of 206Pb/238U-ages and errors of the

example dataset by Ludwig (2003).

LudwigKDE: an object of class 'other'containing the 206Pb/238U-ages (but not the errors) of the

example dataset by Ludwig (2003).

LudwigSpectrum: an object of class 'other'containing the 39Ar abundances, 40Ar/39Ar-ages and

errors of the example dataset by Ludwig (2003).

LudwigMixture: an object of class 'other'containing a dataset of dispersed zircon ﬁssion track

ages of the example dataset by Ludwig (2003).

22 examples

References

Barfod, G.H., Albarede, F., Knoll, A.H., Xiao, S., Telouk, P., Frei, R. and Baker, J., 2002. New

Lu-Hf and Pb-Pb age constraints on the earliest animal fossils. Earth and Planetary Science Letters,

201(1), pp.203-212.

Compston, W., Berry, H., Vernon, M.J., Chappell, B.W. and Kaye, M.J., 1971. Rubidium-strontium

chronology and chemistry of lunar material from the Ocean of Storms. In Lunar and Planetary

Science Conference Proceedings (Vol. 2, p. 1471).

Connelly, J.N., Bollard, J. and Bizzarro, M., 2017. Pb-Pb chronometry and the early Solar System.

Geochimica et Cosmochimica Acta, 201, pp.345-363.

Galbraith, R. F. and Green, P. F., 1990: Estimating the component ages in a ﬁnite mixture, Nuclear

Tracks and Radiation Measurements, 17, 197-206.

Kamo, S.L., Czamanske, G.K. and Krogh, T.E., 1996. A minimum U-Pb age for Siberian ﬂood-

basalt volcanism. Geochimica et Cosmochimica Acta, 60(18), 3505-3511.

Ludwig, K. R., and D. M. Titterington., 1994. "Calculation of 230Th/U isochrons, ages, and errors."

Geochimica et Cosmochimica Acta 58.22, 5031-5042.

Ludwig, K. R., 2003. User’s manual for Isoplot 3.00: a geochronological toolkit for Microsoft

Excel. No. 4.

Lugmair, G.W., Scheinin, N.B. and Marti, K., 1975. Sm-Nd age and history of Apollo 17 basalt

75075-Evidence for early differentiation of the lunar exterior. In Lunar and Planetary Science

Conference Proceedings (Vol. 6, pp. 1419-1429).

Selby, D., 2007. Direct Rhenium-Osmium age of the Oxfordian-Kimmeridgian boundary, Stafﬁn

bay, Isle of Skye, UK, and the Late Jurassic time scale. Norsk Geologisk Tidsskrift, 87(3), p.291.

Vermeesch, P. and Garzanti, E., 2015. Making geological sense of ‘Big Data’ in sedimentary prove-

nance analysis. Chemical Geology, 409, pp.20-27.

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

249(3),pp.339-347.

Examples

data(examples)

concordia(examples$UPb)

agespectrum(examples$ArAr)

isochron(examples$ReOs)

radialplot(examples$FT1)

helioplot(examples$UThHe)

evolution(examples$ThU)

kde(examples$DZ)

radialplot(examples$LudwigMixture)

helioplot 23

agespectrum(examples$LudwigSpectrum)

weightedmean(examples$LudwigMean)

helioplot Visualise U-Th-He data on a logratio plot or ternary diagram

Description

Plot U-Th(-Sm)-He data on a (log[He/Th] vs. log[U/He]) logratio plot or U-Th-He ternary diagram

Usage

helioplot(x, logratio = TRUE, model = 1, show.central.comp = TRUE,

show.numbers = FALSE, alpha = 0.05, contour.col = c("white", "red"),

levels = NA, clabel = "", ellipse.col = c("#00FF0080", "#0000FF80"),

sigdig = 2, xlim = NA, ylim = NA, fact = NA, ...)

Arguments

xan object of class UThHe

logratio Boolean ﬂag indicating whether the data should be shown on bivariate log[He/Th]

vs. log[U/He] diagram, or a U-Th-He ternary diagram.

model choose one of the following statistical models:

1: weighted mean. This model assumes that the scatter between the data points is

solely caused by the analytical uncertainty. If the 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 ﬁrst of these is to assume that the

analytical uncertainties have been underestimated by a factor √M SW D.

2: unweighted mean. A second way to deal with over- or underdispersed datasets

is to simply ignore the analytical uncertainties.

3: weighted mean with overdispersion: instead of attributing any overdispersion

(MSWD > 1) to underestimated analytical uncertainties (model 1), it can also be

attributed to the presence of geological uncertainty, which manifests itself as an

added (co)variance term.

show.central.comp

show the geometric mean composition as a white ellipse?

show.numbers show the grain numbers inside the error ellipses?

alpha probability cutoff for the error ellipses and conﬁdence intervals

contour.col two-element vector with the ﬁll colours to be assigned to the minimum and

maximum age contour

levels a vector with additional values to be displayed as different background colours

within the error ellipses.

24 helioplot

clabel label of the colour scale

ellipse.col a vector of two background colours for the error ellipses. If levels=NA, then

only the ﬁrst colour will be used. If levels is a vector of numbers, then

ellipse.col is used to construct a colour ramp.

sigdig number of signiﬁcant digits for the central age

xlim optional limits of the x-axis (log[U/He]) of the logratio plot. If xlim=NA, the

axis limits are determined automatically.

ylim optional limits of the y-axis (log[Th/He]) of the logratio plot. If ylim=NA, the

axis limits are determined automatically.

fact three-element vector with scaling factors of the ternary diagram if fact=NA,

these will be determined automatically

... optional arguments to the generic plot function

Details

U, Th, Sm and He are compositional data. This means that it is not so much the absolute concentra-

tions of these elements that bear the chronological information, but rather their relative proportions.

The space of all possible U-Th-He compositions ﬁts within the constraints of a ternary diagram

or ‘helioplot’ (Vermeesch, 2008, 2010). If Sm is included as well, then this expands to a three-

dimensional tetrahaedral space (Vermeesch, 2008). Data that ﬁt within these constrained spaces

must be subjected to a logratio transformation prior to statistical analysis (Aitchison, 1986). In the

case of the U-Th-He-(Sm)-He system, this is achieved by ﬁrst deﬁning two (or three) new variables:

u≡ln[U/He]v≡ln[T h/He] (, w ≡ln[Sm/He])

and then performing the desired statistical analysis (averaging, uncertainty propagation, ...) on the

transformed data. Upon completion of the mathematical operations, the results can then be mapped

back to U-Th-(Sm)-He space using an inverse logratio transformation:

[He] = 1/[eu+ev+ (ew) + 1],[U] = eu/[eu+ev+ (ew) + 1]

[T h] = ev/[eu+ev+ (ew) + 1],([Sm] = ew/[eu+ev+ (ew) + 1]).

where [He] + [U] + [T h](+[Sm]) = 1. In the context of U-Th-(Sm)-He dating, the central age

is deﬁned as the age that corresponds to the arithmetic mean composition in logratio space, which

is equivalent to the geometric mean in compositional dataspace (Vermeesch, 2008). IsoplotR’s

helioplot function performs this calculation using the same algorithm that is used to obtain the

weighted mean U-Pb composition for the concordia age calculation. Overdispersion is treated sim-

ilarly as in a regression context (see isochron). Thus, there are options to augment the uncertainties

with a factor √MSW D (model 1); to ignore the analytical uncertainties altogether (model 2); or

to add a constant overdispersion term to the analytical uncertainties (model 3). The helioplot

function visualises U-Th-(Sm)-He data on either a ternary diagram or a bivariate ln[T h/U]vs.

ln[U/He]contour plot. These diagrams provide a convenient way to simultaneously display the

isotopic composition of samples as well as their chronological meaning. In this respect, they fulﬁl

the same purpose as the U-Pb concordia diagram and the U-series evolution plot.

References

Aitchison, J., 1986, The statistical analysis of compositional data: London, Chapman and Hall, 416

isochron 25

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

249(3), pp.339-347.

Vermeesch, P., 2010. HelioPlot, and the treatment of overdispersed (U-Th-Sm)/He data. Chemical

Geology, 271(3), pp.108-111.

See Also

radialplot

Examples

data(examples)

helioplot(examples$UThHe)

dev.new()

helioplot(examples$UThHe,logratio=FALSE)

isochron Calculate and plot isochrons

Description

Plots cogenetic Ar-Ar, Pb-Pb, Rb-Sr, Sm-Nd, Re-Os, Lu-Hf, U-Th-He or Th-U data as X-Y scatter-

plots, ﬁts an isochron curve through them using the york function, and computes the corresponding

isochron age, including decay constant uncertainties.

Usage

isochron(x, ...)

## Default S3 method:

isochron(x, xlim = NA, ylim = NA, alpha = 0.05,

sigdig = 2, show.numbers = FALSE, levels = NA, clabel = "",

ellipse.col = c("#00FF0080", "#FF000080"), ci.col = "gray80",

line.col = "black", lwd = 1, title = TRUE, model = 1, xlab = "x",

ylab = "y", ...)

## S3 method for class 'ArAr'

isochron(x, xlim = NA, ylim = NA, alpha = 0.05,

sigdig = 2, show.numbers = FALSE, levels = NA, clabel = "",

ellipse.col = c("#00FF0080", "#FF000080"), inverse = TRUE,

ci.col = "gray80", line.col = "black", lwd = 1, plot = TRUE,

exterr = TRUE, model = 1, ...)

## S3 method for class 'PbPb'

isochron(x, xlim = NA, ylim = NA, alpha = 0.05,

sigdig = 2, show.numbers = FALSE, levels = NA, clabel = "",

ellipse.col = c("#00FF0080", "#FF000080"), inverse = TRUE,

26 isochron

ci.col = "gray80", line.col = "black", lwd = 1, plot = TRUE,

exterr = TRUE, model = 1, growth = FALSE, ...)

## S3 method for class 'RbSr'

isochron(x, xlim = NA, ylim = NA, alpha = 0.05,

sigdig = 2, show.numbers = FALSE, levels = NA, clabel = "",

ellipse.col = c("#00FF0080", "#FF000080"), ci.col = "gray80",

line.col = "black", lwd = 1, plot = TRUE, exterr = TRUE, model = 1,

...)

## S3 method for class 'ReOs'

isochron(x, xlim = NA, ylim = NA, alpha = 0.05,

sigdig = 2, show.numbers = FALSE, levels = NA, clabel = "",

ellipse.col = c("#00FF0080", "#FF000080"), ci.col = "gray80",

line.col = "black", lwd = 1, plot = TRUE, exterr = TRUE, model = 1,

...)

## S3 method for class 'SmNd'

isochron(x, xlim = NA, ylim = NA, alpha = 0.05,

sigdig = 2, show.numbers = FALSE, levels = NA, clabel = "",

ellipse.col = c("#00FF0080", "#FF000080"), ci.col = "gray80",

line.col = "black", lwd = 1, plot = TRUE, exterr = TRUE, model = 1,

...)

## S3 method for class 'LuHf'

isochron(x, xlim = NA, ylim = NA, alpha = 0.05,

sigdig = 2, show.numbers = FALSE, levels = NA, clabel = "",

ellipse.col = c("#00FF0080", "#FF000080"), ci.col = "gray80",

line.col = "black", lwd = 1, plot = TRUE, exterr = TRUE, model = 1,

...)

## S3 method for class 'ThU'

isochron(x, type = 2, xlim = NA, ylim = NA, alpha = 0.05,

sigdig = 2, show.numbers = FALSE, levels = NA, clabel = "",

ellipse.col = c("#00FF0080", "#FF000080"), ci.col = "gray80",

line.col = "black", lwd = 1, plot = TRUE, exterr = TRUE, model = 1,

...)

## S3 method for class 'UThHe'

isochron(x, xlim = NA, ylim = NA, alpha = 0.05,

sigdig = 2, show.numbers = FALSE, ci.col = "gray80",

line.col = "black", lwd = 1, plot = TRUE, model = 1, ...)

Arguments

xEITHER a matrix with the following ﬁve columns:

Xthe x-variable

sX the standard error of X

isochron 27

Ythe y-variable

sY the standard error of Y

rXY the correlation coefﬁcient of Xand Y

an object of class ArAr,PbPb,ReOs,RbSr,SmNd,LuHf,UThHe or ThU.

... optional arguments to be passed on to the generic plot function if model=2

xlim 2-element vector with the x-axis limits

ylim 2-element vector with the y-axis limits

alpha conﬁdence cutoff for the error ellipses and conﬁdence intervals

sigdig the number of signiﬁcant digits of the numerical values reported in the title of

the graphical output

show.numbers logical ﬂag (TRUE to show grain numbers)

levels a vector with additional values to be displayed as different background colours

within the error ellipses.

clabel label for the colour scale

ellipse.col a vector of two background colours for the error ellipses. If levels=NA, then

only the ﬁrst colour will be used. If levels is a vector of numbers, then

ellipse.col is used to construct a colour ramp.

ci.col the ﬁll colour for the conﬁdence interval of the intercept and slope.

line.col colour of the isochron line

lwd line width

title add a title to the plot?

model construct the isochron using either:

1. Error-weighted least squares regression

2. Ordinary least squares regression

3. Error-weighted least squares with overdispersion term

xlab text label for the horizontal plot axis

ylab text label for the vertical plot axis

inverse if FALSE and xhas class ArAr, plots 40Ar/36Ar vs. 39Ar/36Ar.

if FALSE and xhas class PbPb, plots 207Pb/204Pb vs. 206Pb/204Pb.

if TRUE and xhas class ArAr, plots 36Ar/40Ar vs. 39Ar/40Ar.

if TRUE and xhas class PbPb, plots 207Pb/206Pb vs. 204Pb/206Pb.

plot if FALSE, suppresses the graphical output

exterr propagate external sources of uncertainty (J, decay constant)?

growth add Stacey-Kramers Pb-evolution curve to the plot?

type following the classiﬁcation of Ludwig and Titterington (1994), one of either:

1. ‘Rosholt type-II’ isochron, setting out 230Th/232Th vs. 238U/232Th

2. ‘Osmond type-II’ isochron, setting out 230Th/238U vs. 232Th/238U

3. ‘Rosholt type-II’ isochron, setting out 234U/232Th vs. 238U/232Th

4. ‘Osmond type-II’ isochron, setting out 234U/238U vs. 232Th/238U

28 isochron

Details

Given several aliquots from a single sample, isochrons allow the non-radiogenic component of the

daughter nuclide to be quantiﬁed and separated from the radiogenic component. In its simplest

form, an isochron is obtained by setting out the amount of radiogenic daughter against the amount

of radioactive parent, both normalised to a non-radiogenic isotope of the daughter element, and

ﬁtting a straight line through these points by least squares regression (Nicolaysen, 1961). The slope

and intercept then yield the radiogenic daughter-parent ratio and the non-radiogenic daughter com-

position, respectively. There are several ways to ﬁt an isochron. The easiest of these is ordinary

least squares regression, which weighs all data points equally. In the presence of quantiﬁable an-

alytical uncertainty, it is equally straightforward to use the inverse of the y-errors as weights. It

is signiﬁcantly more difﬁcult to take into account uncertainties in both the x- and the y-variable

(York, 1966). IsoplotR does so for its U-Th-He isochron calculations. The York (1966) method

assumes that the analytical uncertainties of the x- and y-variables are independent from each other.

This assumption is rarely met in geochronology. York (1968) addresses this issue with a bivariate

error weighted linear least squares algorithm that accounts for covariant errors in both variables.

This algorithm was further improved by York et al. (2004) to ensure consistency with the maximum

likelihood approach of Titterington and Halliday (1979).

IsoplotR uses the York et al. (2004) algorithm for its 40Ar/39Ar, Pb-Pb, Rb-Sr, Sm-Nd, Re-

Os and Lu-Hf isochrons. The maximum likelihood algorithm of Titterington and Halliday (1979)

was generalised from two to three dimensions by Ludwig and Titterington (1994) for U-series

disequilibrium dating. Also this algorithm is implemented in IsoplotR. The extent to which the

observed scatter in the data can be explained by the analytical uncertainties can be assessed using

the Mean Square of the Weighted Deviates (MSWD, McIntyre et al., 1966), which is deﬁned as:

MSW D = ([X−ˆ

X]Σ−1

X[X−ˆ

X]T)/df

where Xare the data, ˆ

Xare the ﬁtted values, and ΣXis the covariance matrix of X, and df =

k(n−1) are the degrees of freedom, where kis the dimensionality of the linear ﬁt. MSWD values

that are far smaller or greater than 1 indicate under- or overdispersed measurements, respectively.

Underdispersion can be attributed to overestimated analytical uncertainties. IsoplotR provides

three alternative strategies to deal with overdispersed data:

1. Attribute the overdispersion to an underestimation of the analytical uncertainties. In this case,

the excess scatter can be accounted for by inﬂating those uncertainties by a factor √MSW D.

2. Ignore the analytical uncertainties and perform an ordinary least squares regression.

3. Attribute the overdispersion to the presence of ‘geological scatter’. In this case, the excess

scatter can be accounted for by adding an overdispersion term that lowers the MSWD to unity.

Value

If xhas class PbPb,ArAr,RbSr,SmNd,ReOs or LuHf, or UThHe, returns a list with the following

items:

athe intercept of the straight line ﬁt and its standard error.

bthe slope of the ﬁt and its standard error.

cov.ab the covariance of the slope and intercept

df the degrees of freedom of the linear ﬁt (df =n−2)

isochron 29

y0 a four-element list containing:

y: the atmospheric 40Ar/36Ar or initial 207Pb/204Pb, 187Os/188Os, 87Sr/87Rb, 143Nd/144Nd or

176Hf/177Hf ratio.

s[y]: the propagated uncertainty of y

ci[y]: the studentised 100(1 −α)% conﬁdence interval for y.

disp[y]: the studentised 100(1 −α)% conﬁdence interval for yenhanced by √mswd (only

applicable if model=1).

age a four-element list containing:

t: the 207Pb/206Pb, 40Ar/39Ar, 187Os/187Re, 87Sr/87Rb, 143Nd/144Nd or 176Hf/177Hf age.

s[t]: the propagated uncertainty of t

ci[t]: the studentised 100(1 −α)% conﬁdence interval for t.

disp[t]: the studentised 100(1 −α)% conﬁdence interval for tenhanced by √mswd (only

applicable if model=1).

mswd the mean square of the residuals (a.k.a ‘reduced Chi-square’) statistic (omitted if model=2).

p.value the p-value of a Chi-square test for linearity (omitted if model=2)

wthe overdispersion term, i.e. a three-element vector with the standard deviation of the (as-

sumedly) Normally distributed geological scatter that underlies the measurements, and the

lower and upper half-widths of its 100(1−α)% conﬁdence interval (only returned if model=3).

OR, if xhas class ThU:

par if x$type=1 or x$type=3: the best ﬁtting 230Th/232Th intercept, 230Th/238U slope, 234U/232Th

intercept and 234U/238U slope, OR, if x$type=2 or x$type=4: the best ﬁtting 234U/238U

intercept, 230Th/232Th slope, 234U/238U intercept and 234U/232Th slope.

cov the covariance matrix of par.

df the degrees of freedom for the linear ﬁt, i.e. (3n−3) if x$format=1 or x$format=2, and (2n−2)

if x$format=3 or x$format=4

aif type=1: the 230Th/232Th intercept; if type=2: the 230Th/238U intercept; if type=3: the

234Th/232Th intercept; if type=4: the 234Th/238U intercept and its propagated uncertainty.

bif type=1: the 230Th/238U slope; if type=2: the 230Th/232Th slope; if type=3: the 234U/238U

slope; if type=4: the 234U/232Th slope and its propagated uncertainty.

cov.ab the covariance between aand b.

mswd the mean square of the residuals (a.k.a ‘reduced Chi-square’) statistic.

p.value the p-value of a Chi-square test for linearity.

tfact the 100(1 −α/2)% percentile of a t-distribution with df degrees of freedom.

y0 a four-element vector containing:

y: the initial 234U/238U-ratio

s[y]: the propagated uncertainty of y

ci[y]: the studentised 100(1 −α)% conﬁdence interval for y.

disp[y]: the studentised 100(1 −α)% conﬁdence interval for yenhanced by √mswd.

30 isochron

age a three (or four) element vector containing:

t: the initial 234U/238U-ratio

s[t]: the propagated uncertainty of t

ci[t]: the studentised 100(1 −α)% conﬁdence interval for t

disp[t]: the studentised 100(1 −α)% conﬁdence interval for tenhanced by √mswd (only

reported if model=1).

wthe overdispersion term, i.e. a three-element vector with the standard deviation of the (as-

sumedly) Normally distributed geological scatter that underlies the measurements, and the

lower and upper half-width of its 100(1−α)% conﬁdence interval (only returned if model=3).

da matrix with the following columns: the X-variable for the isochron plot, the analytical un-

certainty of X, the Y-variable for the isochron plot, the analytical uncertainty of Y, and the

correlation coefﬁcient between X and Y.

xlab the x-label of the isochron plot

ylab the y-label of the isochron plot

References

Ludwig, K.R. and Titterington, D.M., 1994. Calculation of 230Th/U isochrons, ages, and errors.

Geochimica et Cosmochimica Acta, 58(22), pp.5031-5042.

Nicolaysen, L.O., 1961. Graphic interpretation of discordant age measurements on metamorphic

rocks. Annals of the New York Academy of Sciences, 91(1), pp.198-206.

Titterington, D.M. and Halliday, A.N., 1979. On the ﬁtting of parallel isochrons and the method of

maximum likelihood. Chemical Geology, 26(3), pp.183-195.

York, D., 1966. Least-squares ﬁtting of a straight line. Canadian Journal of Physics, 44(5), pp.1079-

1086.

York, D., 1968. Least squares ﬁtting of a straight line with correlated errors. Earth and Planetary

Science Letters, 5, pp.320-324.

York, D., Evensen, N.M., Martinez, M.L. and De Basebe Delgado, J., 2004. Uniﬁed equations

for the slope, intercept, and standard errors of the best straight line. American Journal of Physics,

72(3), pp.367-375.

See Also

york,titterington,ludwig

Examples

data(examples)

isochron(examples$ArAr)

fit <- isochron(examples$PbPb,inverse=FALSE,plot=FALSE)

dev.new()

isochron(examples$ThU,type=4)

IsoplotR 31

IsoplotR library(IsoplotR)

Description

A list of documented functions may be viewed by typing help(package='IsoplotR'). Detailed

instructions are provided at http://isoplotr.london-geochron.com. A manuscript with the

theoretical background is in review for publication in Geoscience Frontiers.

Author(s)

Maintainer: Pieter Vermeesch <p.vermeesch@ucl.ac.uk>

See Also

Useful links:

•http://isoplotr.london-geochron.com

kde Create (a) kernel density estimate(s)

Description

Creates one or more kernel density estimates using a combination of the Botev (2010) bandwidth

selector and the Abramson (1982) adaptive kernel bandwidth modiﬁer.

Usage

kde(x, ...)

## Default S3 method:

kde(x, from = NA, to = NA, bw = NA, adaptive = TRUE,

log = FALSE, n = 512, plot = TRUE, pch = NA, xlab = "age [Ma]",

ylab = "", kde.col = rgb(1, 0, 1, 0.6), hist.col = rgb(0, 1, 0, 0.2),

show.hist = TRUE, bty = "n", binwidth = NA, ...)

## S3 method for class 'UPb'

kde(x, from = NA, to = NA, bw = NA, adaptive = TRUE,

log = FALSE, n = 512, plot = TRUE, pch = NA, xlab = "age [Ma]",

ylab = "", kde.col = rgb(1, 0, 1, 0.6), hist.col = rgb(0, 1, 0, 0.2),

show.hist = TRUE, bty = "n", binwidth = NA, type = 4,

cutoff.76 = 1100, cutoff.disc = c(-15, 5), common.Pb = 0, ...)

## S3 method for class 'detritals'

kde(x, from = NA, to = NA, bw = NA, adaptive = TRUE,

32 kde

log = FALSE, n = 512, plot = TRUE, pch = NA, xlab = "age [Ma]",

ylab = "", kde.col = rgb(1, 0, 1, 0.6), hist.col = rgb(0, 1, 0, 0.2),

show.hist = TRUE, bty = "n", binwidth = NA, ncol = NA,

samebandwidth = TRUE, normalise = TRUE, ...)

## S3 method for class 'PbPb'

kde(x, from = NA, to = NA, bw = NA, adaptive = TRUE,

log = FALSE, n = 512, plot = TRUE, pch = NA, xlab = "age [Ma]",

ylab = "", kde.col = rgb(1, 0, 1, 0.6), hist.col = rgb(0, 1, 0, 0.2),

show.hist = TRUE, bty = "n", binwidth = NA, common.Pb = 1, ...)

## S3 method for class 'ArAr'

kde(x, from = NA, to = NA, bw = NA, adaptive = TRUE,

log = FALSE, n = 512, plot = TRUE, pch = NA, xlab = "age [Ma]",

ylab = "", kde.col = rgb(1, 0, 1, 0.6), hist.col = rgb(0, 1, 0, 0.2),

show.hist = TRUE, bty = "n", binwidth = NA, i2i = FALSE, ...)

## S3 method for class 'ThU'

kde(x, from = NA, to = NA, bw = NA, adaptive = TRUE,

log = FALSE, n = 512, plot = TRUE, pch = NA, xlab = "age [ka]",

ylab = "", kde.col = rgb(1, 0, 1, 0.6), hist.col = rgb(0, 1, 0, 0.2),

show.hist = TRUE, bty = "n", binwidth = NA, i2i = FALSE,

detritus = 0, Th02 = c(0, 0), Th02U48 = c(0, 0, 1e+06, 0, 0, 0, 0, 0,

0), ...)

## S3 method for class 'ReOs'

kde(x, from = NA, to = NA, bw = NA, adaptive = TRUE,

log = FALSE, n = 512, plot = TRUE, pch = NA, xlab = "age [Ma]",

ylab = "", kde.col = rgb(1, 0, 1, 0.6), hist.col = rgb(0, 1, 0, 0.2),

show.hist = TRUE, bty = "n", binwidth = NA, i2i = TRUE, ...)

## S3 method for class 'SmNd'

kde(x, from = NA, to = NA, bw = NA, adaptive = TRUE,

log = FALSE, n = 512, plot = TRUE, pch = NA, xlab = "age [Ma]",

ylab = "", kde.col = rgb(1, 0, 1, 0.6), hist.col = rgb(0, 1, 0, 0.2),

show.hist = TRUE, bty = "n", binwidth = NA, i2i = TRUE, ...)

## S3 method for class 'RbSr'

kde(x, from = NA, to = NA, bw = NA, adaptive = TRUE,

log = FALSE, n = 512, plot = TRUE, pch = NA, xlab = "age [Ma]",

ylab = "", kde.col = rgb(1, 0, 1, 0.6), hist.col = rgb(0, 1, 0, 0.2),

show.hist = TRUE, bty = "n", binwidth = NA, i2i = TRUE, ...)

## S3 method for class 'LuHf'

kde(x, from = NA, to = NA, bw = NA, adaptive = TRUE,

log = FALSE, n = 512, plot = TRUE, pch = NA, xlab = "age [Ma]",

ylab = "", kde.col = rgb(1, 0, 1, 0.6), hist.col = rgb(0, 1, 0, 0.2),

show.hist = TRUE, bty = "n", binwidth = NA, i2i = TRUE, ...)

kde 33

## S3 method for class 'UThHe'

kde(x, from = NA, to = NA, bw = NA, adaptive = TRUE,

log = FALSE, n = 512, plot = TRUE, pch = NA, xlab = "age [Ma]",

ylab = "", kde.col = rgb(1, 0, 1, 0.6), hist.col = rgb(0, 1, 0, 0.2),

show.hist = TRUE, bty = "n", binwidth = NA, ...)

## S3 method for class 'fissiontracks'

kde(x, from = NA, to = NA, bw = NA,

adaptive = TRUE, log = FALSE, n = 512, plot = TRUE, pch = NA,

xlab = "age [Ma]", ylab = "", kde.col = rgb(1, 0, 1, 0.6),

hist.col = rgb(0, 1, 0, 0.2), show.hist = TRUE, bty = "n",

binwidth = NA, ...)

Arguments

xa vector of numbers OR an object of class UPb,PbPb,ArAr,ReOs,SmNd,RbSr,

UThHe,fissiontracks,ThU or detrital

... optional arguments to be passed on to R’s density function.

from minimum age of the time axis. If NULL, this is set automatically

to maximum age of the time axis. If NULL, this is set automatically

bw the bandwidth of the KDE. If NULL,bw will be calculated automatically using

the algorithm by Botev et al. (2010).

adaptive logical ﬂag controlling if the adaptive KDE modiﬁer of Abramson (1982) is used

log transform the ages to a log scale if TRUE

nhorizontal resolution (i.e., the number of segments) of the density estimate.

plot show the KDE as a plot

pch the symbol used to show the samples. May be a vector. Set pch=NA to turn them

off.

xlab the x-axis label

ylab the y-axis label

kde.col the ﬁll colour of the KDE speciﬁed as a four element vector of r, g, b, alpha

values

hist.col the ﬁll colour of the histogram speciﬁed as a four element vector of r, g, b, alpha

values

show.hist logical ﬂag indicating whether a histogram should be added to the KDE

bty change to "o","l","7","c","u", or "]" if you want to draw a box around the

plot

binwidth scalar width of the histogram bins, in Myr if log = FALSE, or as a fractional

value if log = TRUE. Sturges’ Rule (log2[n]+1, where nis the number of data

points) is used if binwidth = NA

type scalar indicating whether to plot the 207Pb/235U age (type=1), the 206Pb/238U

age (type=2), the 207Pb/206Pb age (type=3), the 207Pb/206Pb-206Pb/238U age

(type=4), or the (Wetherill) concordia age (type=5)

34 kde

cutoff.76 the age (in Ma) below which the 206Pb/238U and above which the 207Pb/206Pb

age is used. This parameter is only used if type=4.

cutoff.disc two element vector with the minimum (negative) and maximum (positive) per-

centage discordance allowed between the 207Pb/235U and 206Pb/238U age (if

206Pb/238U < cutoff.76) or between the 206Pb/238U and 207Pb/206Pb age (if

206Pb/238U > cutoff.76). Set cutoff.disc=NA if you do not want to use this

ﬁlter.

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

ncol scalar value indicating the number of columns over which the KDEs should be

divided.

samebandwidth logical ﬂag indicating whether the same bandwidth should be used for all sam-

ples. If samebandwidth = TRUE and bw = NULL, then the function will use the

median bandwidth of all the samples.

normalise logical ﬂag indicating whether or not the KDEs should all integrate to the same

value.

i2i ‘isochron to intercept’: calculates the initial (aka ‘inherited’, ‘excess’, or ‘com-

mon’) 40Ar/36Ar, 87Sr/86Sr, 143Nd/144Nd, 187Os/188Os, 230Th/232Th or 176Hf/177Hf

ratio from an isochron ﬁt. Setting i2i to FALSE uses the default values stored in

settings('iratio',...).

detritus detrital 230Th correction (only applicable when x$format == 1 or 2.

0: no correction

1: project the data along an isochron ﬁt

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

3: correct the data using the measured present day 230Th/238U, 232Th/238U and

234U/238U-ratios in the detritus.

Th02 2-element vector with the assumed initial 230Th/232Th-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

Given a set of nage estimates {t1, t2, ..., tn}, histograms and KDEs are probability density esti-

mators that display age distributions by smoothing. Histograms do this by grouping the data into a

number of regularly spaced bins. Alternatively, kernel density estimates (KDEs; Vermeesch, 2012)

smooth data by applying a (Gaussian) kernel:

KDE(t) = Pn

i=1 N(t|µ=ti, σ =h[t])/n

where N(t|µ, σ)is the probability of observing a value tunder a Normal distribution with mean

µand standard deviation σ.h[t]is the smoothing parameter or ‘bandwidth’ of the kernel density

kde 35

estimate, which may or may not depend on the age t. If h[t]depends on t, then KDE(t)is known

as an ‘adaptive’ KDE. The default bandwidth used by IsoplotR is calculated using the algorithm

of Botev et al. (2010) and modulated by the adaptive smoothing approach of Abramson (1982).

The rationale behind adaptive kernel density estimation is to use a narrower bandwidth near the

peaks of the sampling distribution (where the ordered dates are closely spaced in time), and a

wider bandwidth in the distribution’s sparsely sampled troughs. Thus, the resolution of the density

estimate is optimised according to data availability.

Value

If xhas class UPb,PbPb,ArAr,ReOs,SmNd,RbSr,UThHe,fissiontracks or ThU, returns an object

of class KDE, i.e. a list containing the following items:

xhorizontal plot coordinates

yvertical plot coordinates

bw the base bandwidth of the density estimate

ages the data values from the input to the kde function

log copied from the input

or, if xhas class =detritals, an object of class KDEs, i.e. a list containing the following items:

kdes a named list with objects of class KDE

from the beginning of the common time scale

to the end of the common time scale

themax the maximum probability density of all the KDEs

xlabel the x-axis label to be used by plot.KDEs(...)

References

Abramson, I.S., 1982. On bandwidth variation in kernel estimates-a square root law. The annals of

Statistics, pp.1217-1223.

Botev, Z. I., J. F. Grotowski, and D. P. Kroese. "Kernel density estimation via diffusion." The

Annals of Statistics 38.5 (2010): 2916-2957.

Vermeesch, P., 2012. On the visualisation of detrital age distributions. Chemical Geology, 312,

pp.190-194.

See Also

radialplot,cad

Examples

kde(examples$UPb)

dev.new()

kde(examples$FT1,log=TRUE)

dev.new()

kde(examples$DZ,from=1,to=3000,kernel="epanechnikov")

36 ludwig

ludwig Linear regression of U-Pb data with correlated errors, taking into ac-

count decay constant uncertainties.

Description

Implements the maximum likelihood algorithm for Total-Pb/U isochron regression of Ludwig (1998)

Usage

ludwig(x, ...)

## Default S3 method:

ludwig(x, ...)

## S3 method for class 'UPb'

ludwig(x, exterr = FALSE, alpha = 0.05, model = 1, ...)

Arguments

xan object of class UPb

... optional arguments

exterr propagate external sources of uncertainty (e.g., decay constants)?

alpha cutoff value for conﬁdence intervals

model one of three regression models:

1: ﬁt 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 (207Pb/206Pb)◦-ratio intercept (for Tera-Wasserburg). If MSW D>0,

then the analytical uncertainties are augmented by a factor √MSW D.

2: ﬁt a discordia line ignoring the analytical uncertainties

3: ﬁt a discordia line using a modiﬁed maximum likelihood algorithm that in-

cludes accounts for any overdispersion by adding a geological (co)variance term.

Details

The 3-dimensional regression algorithm of Ludwig and Titterington (1994) was modiﬁed by Lud-

wig (1998) to ﬁt so-called ‘Total Pb-U isochrons’. These are constrained to a radiogenic endmember

composition that falls on the concordia line. In its most sophisticated form, this algorithm does not

only allow for correlated errors between variables, but also between aliquots. IsoplotR currently

uses this algorithm to propagate decay constant uncertainties in the total Pb-U isochron ages. Future

versions of the program will generalise this approach to other chronometers as well.

mds 37

Value

par a two-element vector with the lower concordia intercept and initial 207Pb/206Pb-ratio.

cov the covariance matrix of par

df the degrees of freedom of the model ﬁt (3n−3, where nis the number of aliquots).

mswd the mean square of weighted deviates (a.k.a. reduced Chi-square statistic) for the ﬁt.

p.value p-value of a Chi-square test for the linear ﬁt

wthe overdispersion, i.e., a three-element vector with the estimated standard deviation of the (as-

sumedly) Normal distribution that underlies the true isochron; and the lower and upper half-

widths of its 100(1 −α)% conﬁdence interval (only relevant if model = 3).

References

Ludwig, K.R., 1998. On the treatment of concordant uranium-lead ages. Geochimica et Cos-

mochimica Acta, 62(4), pp.665-676.

Ludwig, K.R. and Titterington, D.M., 1994. Calculation of 230Th/U isochrons, ages, and errors.

Geochimica et Cosmochimica Acta, 58(22), pp.5031-5042.

See Also

concordia,titterington,isochron

Examples

f <- system.file("UPb4.csv",package="IsoplotR")

d <- read.data(f,method="U-Pb",format=4)

fit <- ludwig(d)

mds Multidimensional Scaling

Description

Performs classical or nonmetric Multidimensional Scaling analysis

Usage

mds(x, ...)

## Default S3 method:

mds(x, classical = FALSE, plot = TRUE, shepard = FALSE,

nnlines = FALSE, pos = NULL, col = "black", bg = "white", xlab = "",

ylab = "", ...)

## S3 method for class 'detritals'

mds(x, classical = FALSE, plot = TRUE,

shepard = FALSE, nnlines = FALSE, pos = NULL, col = "black",

bg = "white", xlab = "", ylab = "", ...)

38 mds

Arguments

xa dissimilarity matrix OR an object of class detrital

... optional arguments to the generic plot function

classical logical ﬂag indicating whether classical (TRUE) or nonmetric (FALSE) MDS should

be used

plot show the MDS conﬁguration (if shepard=FALSE) or Shepard plot (if shepard=TRUE)

on a graphical device

shepard logical ﬂag indicating whether the graphical output should show the MDS con-

ﬁguration (shepard=FALSE) or a Shepard plot with the ’stress’ value. This ar-

gument is only used if plot=TRUE.

nnlines if TRUE, draws nearest neighbour lines

pos a position speciﬁer for the labels (if par('pch')!=NA). Values of 1, 2, 3 and

4 indicate positions below, to the left of, above and to the right of the MDS

coordinates, respectively.

col plot colour (may be a vector)

bg background colour (may be a vector)

xlab a string with the label of the x axis

ylab a string with the label of the y axis

Details

Multidimensional Scaling (MDS) is a dimension-reducting technique that takes a matrix of pair-

wise ‘dissimilarities’ between objects (e.g., age distributions) as input and produces a conﬁguration

of two (or higher-) dimensional coordinates as output, so that the Euclidean distances between

these coordinates approximate the dissimilarities of the input matrix. Thus, an MDS-conﬁguration

serves as a ‘map’ in which similar samples cluster closely together and dissimilar samples plot

far apart. In the context of detrital geochronology, the dissimilarity between samples is given by

the statistical distance between age distributions. There are many ways to deﬁne this statistical

distance. IsoplotR uses the Kolmogorov-Smirnov (KS) statistic due to its simplicity and the fact

that it behaves like a true distance in the mathematical sense of the word (Vermeesch, 2013). The

KS-distance is given by the maximum vertical distance between two cad step functions. Thus, the

KS-distance takes on values between zero (perfect match between two age distributions) and one

(no overlap between two distributions). Calculating the KS-distance between samples two at a time

populates a symmetric dissimilarity matrix with positive values and a zero diagonal. IsoplotR im-

plements two algorithms to convert this matrix into a conﬁguration. The ﬁrst (‘classical’) approach

uses a sequence of basic matrix manipulations developed by Young and Householder (1938) and

Torgerson (1952) to achieve a linear ﬁt between the KS-distances and the ﬁtted distances on the

MDS conﬁguration. The second, more sophisticated (‘nonmetric’) approach subjects the input dis-

tances to a transformation fprior to ﬁtting a conﬁguration:

δi,j =f(KSi,j )

where KSi,j is the KS-distance between samples iand j(for 1≤i6=j≤n) and δi,j is the

‘disparity’ (Kruskal, 1964). Fitting an MDS conﬁguration then involves ﬁnding the disparity trans-

formation that maximises the goodness of ﬁt (or minimises the ‘stress’) between the disparities and

peakﬁt 39

the ﬁtted distances. The latter two quantities can also be plotted against each other as a ‘Shepard

plot’.

Value

Returns an object of class MDS, i.e. a list containing the following items:

points a two-column vector of the ﬁtted conﬁguration

classical a logical ﬂag indicating whether the MDS conﬁguration was obtained by classical (TRUE)

or nonmetric (FALSE) MDS

diss the dissimilarity matrix used for the MDS analysis

stress (only if classical=TRUE) the ﬁnal stress achieved (in percent)

References

Kruskal, J., 1964. Multidimensional scaling by optimizing goodness of ﬁt to a nonmetric hypothe-

sis. Psychometrika 29 (1), 1-27.

Torgerson, W. S. Multidimensional scaling: I. Theory and method. Psychometrika, 17(4): 401-419,

1952.

Vermeesch, P., 2013. Multi-sample comparison of detrital age distributions. Chemical Geology,

341, pp.140-146.

Young, G. and Householder, A. S. Discussion of a set of points in terms of their mutual distances.

Psychometrika, 3(1):19-22, 1938.

See Also

cad,kde

Examples

data(examples)

mds(examples$DZ,nnlines=TRUE,pch=21,cex=5)

dev.new()

mds(examples$DZ,shepard=TRUE)

peakfit Finite mixture modelling of geochronological datasets

Description

Implements the discrete mixture modelling algorithms of Galbraith and Laslett (1993) and applies

them to ﬁssion track and other geochronological datasets.

40 peakﬁt

Usage

peakfit(x, ...)

## Default S3 method:

peakfit(x, k = "auto", sigdig = 2, log = TRUE,

alpha = 0.05, ...)

## S3 method for class 'fissiontracks'

peakfit(x, k = 1, exterr = TRUE, sigdig = 2,

log = TRUE, alpha = 0.05, ...)

## S3 method for class 'UPb'

peakfit(x, k = 1, type = 4, cutoff.76 = 1100,

cutoff.disc = c(-15, 5), exterr = TRUE, sigdig = 2, log = TRUE,

alpha = 0.05, ...)

## S3 method for class 'PbPb'

peakfit(x, k = 1, exterr = TRUE, sigdig = 2, log = TRUE,

i2i = TRUE, alpha = 0.05, ...)

## S3 method for class 'ArAr'

peakfit(x, k = 1, exterr = TRUE, sigdig = 2, log = TRUE,

i2i = FALSE, alpha = 0.05, ...)

## S3 method for class 'ReOs'

peakfit(x, k = 1, exterr = TRUE, sigdig = 2, log = TRUE,

i2i = TRUE, alpha = 0.05, ...)

## S3 method for class 'SmNd'

peakfit(x, k = 1, exterr = TRUE, sigdig = 2, log = TRUE,

i2i = TRUE, alpha = 0.05, ...)

## S3 method for class 'RbSr'

peakfit(x, k = 1, exterr = TRUE, sigdig = 2, log = TRUE,

i2i = TRUE, alpha = 0.05, ...)

## S3 method for class 'LuHf'

peakfit(x, k = 1, exterr = TRUE, sigdig = 2, log = TRUE,

i2i = TRUE, alpha = 0.05, ...)

## S3 method for class 'ThU'

peakfit(x, k = 1, exterr = FALSE, sigdig = 2, log = TRUE,

i2i = TRUE, alpha = 0.05, detritus = 0, Th02 = c(0, 0),

Th02U48 = c(0, 0, 1e+06, 0, 0, 0, 0, 0, 0), ...)

## S3 method for class 'UThHe'

peakfit(x, k = 1, sigdig = 2, log = TRUE, alpha = 0.05,

...)

peakﬁt 41

Arguments

xeither an [n x 2] matrix with measurements and their standard errors, or an

object of class fissiontracks,UPb,PbPb,ArAr,ReOs,SmNd,RbSr,LuHf,ThU

or UThHe

... optional arguments (not used)

kthe number of discrete age components to be sought. Setting this parameter

to 'auto'automatically selects the optimal number of components (up to a

maximum of 5) using the Bayes Information Criterion (BIC).

sigdig number of signiﬁcant digits to be used for any legend in which the peak ﬁtting

results are to be displayed.

log take the logs of the data before applying the mixture model?

alpha cutoff value for conﬁdence intervals

exterr propagate the external sources of uncertainty into the component age errors?

type scalar valueindicating whether to plot the 207Pb/235U age (type=1), the 206Pb/238U

age (type=2), the 207Pb/206Pb age (type=3), the 207Pb/206Pb-206Pb/238U age

(type=4), or the (Wetherill) concordia age (type=5)

cutoff.76 the age (in Ma) below which the 206Pb/238U and above which the 207Pb/206Pb

age is used. This parameter is only used if type=4.

cutoff.disc two element vector with the maximum and minimum percentage discordance al-

lowed between the 207Pb/235U and 206Pb/238U age (if 206Pb/238U < cutoff.76)

or between the 206Pb/238U and 207Pb/206Pb age (if 206Pb/238U > cutoff.76).

Set cutoff.disc=NA if you do not want to use this ﬁlter.

i2i ‘isochron to intercept’: calculates the initial (aka ‘inherited’, ‘excess’, or ‘com-

mon’) 40Ar/36Ar, 207Pb/204Pb, 87Sr/86Sr, 143Nd/144Nd, 187Os/188Os or 176Hf/177Hf

ratio from an isochron ﬁt. 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.

detritus detrital 230Th correction (only applicable when x$format == 1 or 2.

0: no correction

1: project the data along an isochron ﬁt

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

3: correct the data using the measured present day 230Th/238U, 232Th/238U and

234U/238U-ratios in the detritus.

Th02 2-element vector with the assumed initial 230Th/232Th-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

Consider a dataset of ndates {t1, t2, ..., tn}with analytical uncertainties {s[t1], s[t2], ..., s[tn]}.

Deﬁne zi= log(ti)and s[zi] = s[ti]/ti. Suppose that these nvalues are derived from a mixture

42 peakﬁt

of k > 2populations with means {µ1, ..., µk}. Such a discrete mixture may be mathematically

described by:

P(zi|µ, ω) = Pk

j=1 πjN(zi|µj, s[zj]2)

where πjis the proportion of the population that belongs to the jth component, and πk= 1 −

Pk−1

j=1 πj. This equation can be solved by the method of maximum likelihood (Galbraith and

Laslett, 1993). IsoplotR implements the Bayes Information Criterion (BIC) as a means of au-

tomatically choosing k. This option should be used with caution, as the number of peaks steadily

rises with sample size (n). If one is mainly interested in the youngest age component, then it is more

productive to use an alternative parameterisation, in which all grains are assumed to come from one

of two components, whereby the ﬁrst component is a single discrete age peak (exp(m), say) and

the second component is a continuous distribution (as descibed by the central age model), but

truncated at this discrete value (Van der Touw et al., 1997).

Value

Returns a list with the following items:

peaks a3xkmatrix with the following rows:

t: the ages of the kpeaks

s[t]: the estimated uncertainties of t

ci[t]: the widths of approximate 100(1 −α)% conﬁdence intervals for t

props a2xkmatrix with the following rows:

p: the proportions of the kpeaks

s[p]: the estimated uncertainties (standard errors) of p

Lthe log-likelihood of the ﬁt

legend a vector of text expressions to be used in a ﬁgure legend

References

Galbraith, R.F. and Laslett, G.M., 1993. Statistical models for mixed ﬁssion track ages. Nuclear

Tracks and Radiation Measurements, 21(4), pp.459-470.

van der Touw, J., Galbraith, R., and Laslett, G. A logistic truncated normal mixture model for

overdispersed binomial data. Journal of Statistical Computation and Simulation, 59(4):349-373,

1997.

See Also

radialplot,central

Examples

data(examples)

peakfit(examples$FT1,k=2)

peakfit(examples$LudwigMixture,k='min')

radialplot 43

radialplot Visualise heteroscedastic data on a radial plot

Description

Implementation of a graphical device developed by Rex Galbraith to display several estimates of

the same quantity that have different standard errors.

Usage

radialplot(x, ...)

## Default S3 method: