the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Errorchrons and anchored isochrons in IsoplotR
Isochrons are usually fitted by “York regression”, which uses a weighted leastsquares approach that accounts for correlated uncertainties in both variables. Despite its tremendous popularity in modern geochronology, the York algorithm has two important limitations that reduce its utility in several applications. First, it does not provide a satisfactory mechanism to deal with socalled “errorchrons”, i.e. datasets that are overdispersed with respect to the analytical uncertainties. Second, York regression is not readily amenable to anchoring, in which either the slope or the intercept of the isochron is fixed based on some external information. Anchored isochrons can be very useful in cases where the data are insufficiently spread out to constrain both the radiogenic and nonradiogenic isotopic composition.
This paper addresses both of these issues by extending a maximum likelihood algorithm that was first proposed by Titterington and Halliday (1979). The new algorithm offers the ability to attribute any excess dispersion to either the inherited component (“model 3a”) or diachronous closure of the isotopic system (“model 3b”). It provides an opportunity to anchor isochrons to either a fixed nonradiogenic composition or a fixed age. Last but not least, it allows the user to attach meaningful analytical uncertainty to the anchor. The new method has been implemented in IsoplotR for immediate use in $\mathrm{Ar}/\mathrm{Ar}$, $\mathrm{Pb}/\mathrm{Pb}$, $\mathrm{U}/\mathrm{Pb}$, $\mathrm{Th}/\mathrm{Pb}$, $\mathrm{Rb}/\mathrm{Sr}$, $\mathrm{Sm}/\mathrm{Nd}$, $\mathrm{Lu}/\mathrm{Hf}$, $\mathrm{Re}/\mathrm{Os}$, $\mathrm{K}/\mathrm{Ca}$, and U–Th–He geochronology.
 Article
(1778 KB)  Fulltext XML
 BibTeX
 EndNote
Isochrons are mixing lines between radiogenic and inherited isotopic endmembers. They are an essential component of radiometric geochronology and exist in several forms. Sections 1–7 of this paper will deal with the simple case of “$P/D$ isochrons”, which apply to geochronometers that are based on the decay of a single radioactive parent nuclide (P; e.g. ^{87}Rb, ^{40}K, ^{147}Sm) to a particular daughter nuclide (D; e.g. ^{87}Sr, ^{40}Ar, ^{143}Nd). Sections 8 and 9 will discuss $\mathrm{Pb}/\mathrm{Pb}$ and $\mathrm{U}/\mathrm{Pb}$ isochrons, respectively. These are based on the paired decay of two parents (^{238}U and ^{235}U) to two daughters (^{206}Pb and ^{207}Pb, respectively) and require linear regression in two or three dimensions.
$P/D$ isochrons come as two types. “Conventional” isochrons are straightline relationships of the following kind:
where i is the aliquot number ($\mathrm{1}\le i\le n$), d is a nonradiogenic sister isotope of D, $[D/d{]}_{\mathrm{0}}$ is the inherited component, λ is the decay constant, and t is the time elapsed since isotopic closure. A full list of definitions is provided in Appendix A. “Inverse” isochrons are obtained by permuting P, D, and d:
The choice between conventional and inverse isochrons depends on the relative precision of the mass spectrometer measurements. Inverse isochrons are preferred if d≪D. Otherwise, conventional isochrons are fine (Li and Vermeesch, 2021).
Equations (1) and (2) can be recast into a generic linear form:
Table A1 maps the parameters of this generic equation to the parameters of Eqs. (1) and (2) as well as subsequent variants thereof. Equation (3) is usually solved by York et al. (2004) regression (hereafter simply referred to as “York regression”). York regression uses a leastsquares algorithm to estimate the intercept (a) and slope (b) of the isochron line from an n×5 table of paired isotopic ratio measurements, along with their standard errors and their error correlations.
Although this paper will use the York parameters a and b to fit the isochron, all results will be presented in terms of the geologically more meaningful inherited endmember ratio $[D/d{]}_{\mathrm{0}}$ and radiogenic endmember ratio $[D/P{]}^{\ast}\equiv {e}^{\mathit{\lambda}t}\mathrm{1}$, where $[D/d{]}_{\mathrm{0}}=a$ and $[D/P{]}^{\ast}=b$ for conventional isochrons and $[D/d{]}_{\mathrm{0}}=\mathrm{1}/a$ and $[D/P{]}^{\ast}=b/a$ for inverse isochrons (Li and Vermeesch, 2021).
The most accurate and precise results are obtained from samples that are evenly spread along the isochron line, spanning the entire range of values from the inherited to the radiogenic endmember. Unfortunately, this condition is not always fulfilled. It is not uncommon for most or all aliquots in a sample to cluster together at one point along the isochron, making it difficult to accurately estimate the endmember compositions. For example, no precise isochron ages can be obtained from samples whose radiogenic daughter component is dwarfed by the inherited daughter component. Conversely, the composition of the inherited component cannot be precisely estimated in extremely radiogenic samples. Finally, when the data cluster together in the middle, then neither endmember component can be reliably determined (Fig. 1).
Sometimes, such poorly constrained isochrons can be fixed using external information. For example, if the composition of the inherited component is known through some independent means (e.g. by analysing a cogenetic mineral that is naturally poor in P and rich in D and d), then this information can be used to anchor the isochron. Anchoring reduces the number of unknown parameters by 1, benefitting numerical stability and precision. Conversely, if the age of the sample is known, then the inherited component can be estimated by anchoring to the radiogenic endmember.
There is currently no formally documented way to anchor York regression. A commonly used “hack” is to add an extra data point with infinite precision representing either the inherited or radiogenic endmember component. However, this hack does not provide a satisfactory mechanism to assign uncertainty to the anchor. This paper solves that problem. Section 2 introduces a maximum likelihood formulation of York regression, which is amenable to anchoring with uncertainty. Section 3 shows how this formulation can be used to model “errorchrons” that are overdispersed with respect to the analytical uncertainties.
Section 4 provides further details about the implementation of maximum likelihood regression, which can attribute overdispersion to either the intercept (“model 3a”) or slope (“model 3b”) of the linear fit. It turns out that model 3b is computationally more challenging than model 3a. Section 5 shows how this issue can be avoided by inverting and flipping the isochron axes.
Section 6 shows how anchored isochron regression represents a trivial special case of the maximum likelihood algorithm. Section 7 attaches two different geological interpretations to the uncertainty of isochron anchors. These two interpretations lead to two different fit models (“model 1” and “model 3”, which are so named for historical reasons). Sections 8 and 9 apply the same logic to $\mathrm{Pb}/\mathrm{Pb}$ and $\mathrm{U}/\mathrm{Pb}$ isochrons, respectively. Finally, Sect. 10 shows how the new algorithms can be used within the IsoplotR software package.
York et al. (2004) use the method of least squares to fit the general problem of weighted regression with correlated uncertainties in both variables. Titterington and Halliday (1979) obtained identical results using the method of maximum likelihood. The latter approach offers greater flexibility than least squares and will be used in the remainder of this paper.
The method of maximum likelihood is a standard statistical technique to estimate the parameters of a probability distribution from a set of measurements. In the case of twodimensional isochron regression, the parameters are the intercept a, the slope b, and the true (but unknown) x_{i} values. The data are the X_{i} and Y_{i} measurements along with their assumed uncertainties and error correlations (see Appendix A for a full list of definitions). For a dataset of n aliquots, there are 2n measurements (n X_{i} values and n Y_{i} values) and n+2 parameters (a, b, and n x_{i} values). That leaves $\mathrm{2}nn\mathrm{2}=n\mathrm{2}$ “degrees of freedom” to estimate the parameters.
This overconstrained problem can be solved by minimizing the paired differences (“residuals”) between the true x_{i} and y_{i} values (where ${y}_{i}=a+b{x}_{i}$) and between the measured X_{i} and Y_{i} values.
To make the problem tractable, the vector of residuals Δ_{i} is assumed to be drawn from a bivariate normal distribution with zero mean and covariance matrix Σ_{i}:
and
The product of the normal probabilities for all the aliquots (from $\mathrm{1}\le i\le n$) is called the “likelihood” of the parameters given the data. Maximizing this number produces the most accurate possible parameter estimates $\widehat{a}$, $\widehat{b}$, and ${\widehat{x}}_{i}$. Alternatively, and equivalently, the same estimates can be obtained by maximizing the logarithm of the likelihood. This is the preferred approach as it improves numerical stability and facilitates the estimation of the parameter uncertainties. The loglikelihood function can be formally defined as
This formulation allows for correlated uncertainties between the variables ($s[{X}_{i},{Y}_{i}]\ne \mathrm{0}$ in Eq. 5) but not between aliquots ($s[{X}_{i},{Y}_{j}]=\mathrm{0}$ if i≠j). Daëron and Vermeesch (2023) discuss the generalized case of “omnivariant” regression, where this requirement has been relaxed and Eq. (6) is replaced with a single matrix expression.
The degree to which the residuals Δ_{i} are consistent with the assumed analytical uncertainties Σ_{i} can be assessed using a chisquare test and the MSWD parameter. Readers who are not familiar with these concepts are referred to Appendix C for details. Isochrons that exhibit significant overdispersion with respect to the analytical uncertainties are colloquially referred to as errorchrons. Whether it is right to use such a pejorative term for this common phenomenon is debatable (Schaen et al., 2021). Four approaches may be used to deal with errorchrons.
 Model 1.

Inflate the analytical uncertainties by a factor $\sqrt{\text{MSWD}}$.
 Model 2.

Ignore the analytical uncertainties and replace York regression with orthogonal regression or a similar technique. This approach will not be discussed further in this paper.
 Model 3.

Quantify the dispersion as an additional free parameter. There are two options for doing so.

Model 3a. Attribute the excess dispersion to variability of the inherited component and, hence, the isochron intercept. For conventional isochrons, this means that the true isotopic ratios of the cogenetic aliquots do not belong to a single isochron line but to a family of parallel isochron lines whose intercepts are normally distributed (Titterington and Halliday, 1979). The standard deviation of this normal distribution (σ_{a}) can be used to quantify the dispersion of the data around the “central” isochron line.

Model 3b. Attribute the excess dispersion to diachronous closure of the isotopic system. This mechanism produces families of conventional isochrons that share a common intercept but differ in their slopes. Again, the distribution of the slopes can be assumed to follow a normal distribution, with mean b and standard deviation σ_{b}. The latter value can be used to quantify the dispersion of the data, which in this case captures the degree of diachroneity.

Model 3a isochron regression can be formalized by modifying Eq. (4) as follows:
The dispersion parameter σ_{a} can be estimated (as ${\widehat{\mathit{\sigma}}}_{a}$) by incorporating it into the covariance matrices:
The equivalent expressions for model 3b regression are
and
respectively.
Three specific cases of Eq. (6) are of interest in the present discussion.
Model 1.
Model 3a.
Model 3b.
Note that the (logged) determinant of the covariance matrices is not included in the constant for model 3 fits because Σ_{a,i} and Σ_{b,i} are both functions of x, whereas Σ_{i} is not. ℒℒ can be numerically maximized using an iterative twostep procedure. For model 1 regression, the first step maximizes $\mathcal{L}{\mathcal{L}}_{\mathrm{1}}\left(\mathbf{x}\righta,b,\mathbf{X},\mathbf{Y},\mathbf{\Sigma})$ to find x for any pair of a and b values. The second step repeats the first step for different values of a and b until the overall optimum is reached.
This procedure can easily be adapted to model 3 regression by adding σ_{a} or σ_{b} to the unknown parameters. It works well for model 3a, where the first step has a direct solution. However, it is slow for model 3b regression, in which finding x requires an additional level of iteration. See Appendix D for further details. Fortunately, in practice model 3b is rarely or never needed for geochronology, as explained next.
For inverse isochrons, model 3b regression does not actually provide any useful information. Unlike conventional isochrons, whose chronological information is contained in the slope, the chronological information of inverse isochrons is contained in their horizontal intercept. This information can be unlocked by flipping the axes of the isochron diagram around, inverting the isochron, and treating the $D/P$ ratio as the dependent variable. Applying model 3a regression to the flipped data yields a dispersion estimate σ_{a} that can be used to quantify the age dispersion.
The same trick can be used to estimate the slope uncertainty of a conventional isochron. A pragmatic way to avoid the slow convergence rate of model 3b regression is to invert the isochron, flip the dependent and independent variables around, invert the isochron a second time, and carry out a model 3a regression on the transformed data. The resulting σ_{a} estimate can be converted to a σ_{b} estimate:
Figure 2 applies model 3 regression to two synthetic Ar–Ar datasets with a nonradiogenic ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ ratio of 400 and an age of 100 Ma. The dataset of Fig. 2a exhibits overdispersion of the y intercept ([^{40}Ar $/$ ^{36}Ar${]}_{\mathrm{0}}=\mathrm{400},{\mathit{\sigma}}_{a}=\mathrm{40}$), whereas the dataset of Fig. 2b is overdispersed in the x intercept (t=100 Ma, σ_{t}=10 Ma). In both cases, the maximum likelihood algorithm has retrieved the correct solution from the noisy data.
Anchored isochron regression requires just a trivial modification of the maximum likelihood algorithm. It suffices to treat the anchored parameter as data. For example, the slope and intercept of a model 1 isochron can be anchored by maximizing $\mathcal{L}{\mathcal{L}}_{\mathrm{1}}(a,\mathbf{x}b,\mathbf{X},\mathbf{Y},\mathbf{\Sigma})$ and $\mathcal{L}{\mathcal{L}}_{\mathrm{1}}(b,\mathbf{x}a,\mathbf{X},\mathbf{Y},\mathbf{\Sigma})$, respectively. Note that these calculations assume that the anchor is known exactly. This is reflected in the zero uncertainties of the nonradiogenic component in Fig. 1b and e and of the radiogenic components in Fig. 1c and f. Such absolute certainty is unrealistic in the noisy world of geology. Section 7 introduces two ways to formally account for uncertainty in anchored isochron regression.
When assigning uncertainty to statistical parameters, it is important to clearly define the meaning of this uncertainty. In the case of anchored isochron regression, the uncertainty of the intercept or slope can carry two meanings. For example, when the intercept of an isochron is anchored at a value $\overline{a}\pm {\mathit{\sigma}}_{\overline{a}}$, this can mean either of the following.
 Model 1.

The data are underlain by a single isochron whose intercept is only approximately known, with a most likely value of $\overline{a}$ and a precision (standard error) of ${\mathit{\sigma}}_{\overline{a}}$.
 Model 3.

The data were drawn from a family of isochron lines whose intercepts follow a normal distribution with mean $\overline{a}$ and standard deviation ${\mathit{\sigma}}_{\overline{a}}$.
These two different approaches can be implemented by replacing the loglikelihood functions of Eq. (13) with the following.
The equivalent expressions for anchored slopes with uncertainty are as follows.
Note that model 1 treats the slope and intercept as unknowns, so their maximum likelihood estimates generally do not equal the anchored values $\overline{a}$ and $\overline{b}$. However, the smaller ${\mathit{\sigma}}_{\overline{a}}$ and ${\mathit{\sigma}}_{\overline{b}}$ are relative to $\overline{a}$ and $\overline{b}$, the closer that $\widehat{a}$ and $\widehat{b}$ are to $\overline{a}$ and $\overline{b}$. In contrast, anchored model 3 regression fixes the dispersion of the anchored parameters so that ${\widehat{\mathit{\sigma}}}_{a}={\mathit{\sigma}}_{\overline{a}}$ for anchored intercepts, whereas ${\widehat{\mathit{\sigma}}}_{b}={\mathit{\sigma}}_{\overline{b}}$ for anchored slopes.
Figure 3 applies the two approaches to the synthetic dataset of Fig. 1 using anchors of $[D/d{]}_{\mathrm{0}}=\mathrm{1}\pm \mathrm{0.05}$ and $[D/P{]}^{\ast}=\mathrm{1}\pm \mathrm{0.05}$ (1σ). As expected, models 1 and 3 produce slightly different outcomes, with model 1 resulting in $[D/d{]}_{\mathrm{0}}$ and $[D/P{]}^{\ast}$ estimates that are slightly different than the anchored values, whereas model 3 reproduces the anchors exactly. The model 3 fits produce more precise estimates for the unanchored quantities. This reflects the fact that some of the uncertainty in the isochron is captured by the dispersion parameter, which is absent from the model 1 fit.
Model 3 partitions the analytical and geological uncertainty between the standard errors of $\widehat{a}$ and $\widehat{b}$ on the one hand and the dispersion parameter on the other hand. It is not easy to visualize this partitioned uncertainty. Figure 3 uses a grey confidence envelope to show the analytical uncertainty and a black error bar to visualize the geological dispersion. Together these two items represent the total uncertainty budget.
It is, of course, also possible to treat the dispersion parameter as an unknown by maximizing $\mathcal{L}{\mathcal{L}}_{\mathrm{3}a}(b,{\mathit{\sigma}}_{a},\mathbf{x}\overline{a},\mathbf{X},\mathbf{Y},\mathbf{\Sigma})$ or $\mathcal{L}{\mathcal{L}}_{\mathrm{3}b}(a,{\mathit{\sigma}}_{b},\mathbf{x}\overline{b},\mathbf{X},\mathbf{Y},\mathbf{\Sigma})$. Section 9 will give an example of this in the context of U–Pb geochronology.
The ${}^{\mathrm{207}}\mathrm{Pb}{/}^{\mathrm{206}}\mathrm{Pb}$ method is based on the paired decay of two radioactive parents P_{1} and P_{2} (^{238}U and ^{235}U) to two daughters D_{1} and D_{2} (^{206}Pb and ^{207}Pb, respectively) in the presence of a nonradiogenic sister isotope d (^{204}Pb). This gives rise to the following implicit age equation:
where λ_{1} and λ_{2} are the decay constants of P_{1} and P_{2}, respectively, and $[{P}_{\mathrm{2}}/{P}_{\mathrm{1}}]$ is constant. Equation (17) can be recast into the generic form of Eq. (3) to form a “$D/D$ isochron” using the mapping of Table A1. This gives rise to conventional (${}^{\mathrm{207}}\mathrm{Pb}{/}^{\mathrm{204}}\mathrm{Pb}$ vs. ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{204}}\mathrm{Pb}$) and inverse (${}^{\mathrm{207}}\mathrm{Pb}{/}^{\mathrm{206}}\mathrm{Pb}$ vs. ${}^{\mathrm{204}}\mathrm{Pb}{/}^{\mathrm{206}}\mathrm{Pb}$) isochrons.
Unlike inverse $P/D$ isochrons, which are characterized by negative slopes, inverse $D/D$ isochrons have positive slopes. There is perfect symmetry between conventional and inverse $D/D$ isochrons in the sense that the intercept of a conventional isochron equals the slope of an inverse isochron and vice versa. This slightly changes and, in fact, simplifies the procedure for anchored isochron regression.
There is no need for flipped isochron regression of $\mathrm{Pb}/\mathrm{Pb}$ data. Instead, all possible scenarios can be handled by straightforward model 1 and model 3a regression. Isochrons can be anchored to the nonradiogenic component in conventional isochron space and to the radiogenic component in inverse isochron space.
The $\mathrm{U}/\mathrm{Pb}$ method, like the $\mathrm{Pb}/\mathrm{Pb}$ method, is based on the paired decay of ^{238}U to ^{206}Pb and of ^{235}U to ^{207}Pb. The two chronometers can be treated separately and plotted as $P/D$ isochrons. These twodimensional isochrons can be anchored to the radiogenic and nonradiogenic composition using the methods described in Sects. 2–7.
Alternatively, the two ingrowth equations can also be coupled, forming a threedimensional “total $\mathrm{Pb}/\mathrm{U}$ isochron”. This problem can be solved using the method of maximum likelihood (Ludwig, 1998). There is just one important difference between this solution and the Yorklike algorithm of Sect. 2. York regression is formulated in terms of a generic intercept (a) and slope (b). In contrast, the “Ludwig” algorithm is formulated in terms of the actual nonradiogenic ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{204}}\mathrm{Pb}$ (=α) and ${}^{\mathrm{207}}\mathrm{Pb}{/}^{\mathrm{204}}\mathrm{Pb}$ (=β) ratio, as well as the age (t) of the system.
It is relatively straightforward to generalize the concept of model 3 regression to total $\mathrm{Pb}/\mathrm{U}$ isochrons when the dispersion is attributed to the isochron age. However, it is not immediately clear how to do the same for the nonradiogenic composition because this would require partitioning the dispersion between the α and β parameters. This paper will not attempt to answer this question. Instead, let us conclude the technical discussion by shifting to a far more common type of $\mathrm{U}/\mathrm{Pb}$ data.
Most published $\mathrm{U}/\mathrm{Pb}$ studies do not report ^{204}Pb because this nuclide is rare and difficult to measure. Instead, “semitotal $\mathrm{Pb}/\mathrm{U}$” regression is done in Wetherill or Tera–Wasserburg concordia space, which are akin to conventional and inverse isochron space, respectively. In the absence of a nonradiogenic sister isotope of Pb, the isochron is then redefined as a mixing line between an inherited ${}^{\mathrm{207}}\mathrm{Pb}{/}^{\mathrm{206}}\mathrm{Pb}$ ratio and the radiogenic ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{238}}\mathrm{U}$ and ${}^{\mathrm{207}}\mathrm{Pb}{/}^{\mathrm{235}}\mathrm{U}$ ratios.
The semitotal $\mathrm{Pb}/\mathrm{U}$ problem can be cast into a maximum likelihood format by replacing Eq. (4) with
where X_{i} and Y_{i} are the measured ${}^{\mathrm{207}}\mathrm{Pb}{/}^{\mathrm{235}}\mathrm{U}$ and ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{238}}\mathrm{U}$ ratios of the ith aliquot, respectively, and x_{i} is the true ${}^{\mathrm{207}}{\mathrm{Pb}}_{\mathrm{0}}{/}^{\mathrm{235}}\mathrm{U}$ ratio of the ith aliquot, whilst $a={e}^{{\mathit{\lambda}}_{\mathrm{238}}t}\mathrm{1}$ and $b={[}^{\mathrm{206}}\text{Pb}{/}^{\mathrm{207}}\text{Pb}{]}_{\mathrm{0}}{[}^{\mathrm{235}}\text{U}{/}^{\mathrm{238}}\text{U}]$. Equation (19) can be solved using the same recipe as the York reformulation of Sects. 2–7, with the caveat that the isochron terminates on the concordia line at the $\mathrm{U}/\mathrm{Pb}$ composition corresponding to t.
Figure 4 shows a model 3a anchored semitotal $\mathrm{Pb}/\mathrm{U}$ isochron in Tera–Wasserburg space. The true age for this synthetic dataset is 1500 Ma. The unconstrained model 1 isochron produces an imprecise isochron age of 1491±100 Ma, which is improved to 1504±38 Ma by anchoring to ${[}^{\mathrm{207}}\text{Pb}{/}^{\mathrm{206}}\text{Pb}{]}_{\mathrm{0}}=\mathrm{1.1}$. Just like the $P/D$ isochrons of Fig. 3, the uncertainty budget for the $\mathrm{U}/\mathrm{Pb}$ isochrons of Fig. 4 is partitioned between the analytical and geological dispersion, where the latter is treated as a free parameter with a maximum likelihood estimate of ${\widehat{\mathit{\sigma}}}_{a}=\mathrm{0.053}\pm \mathrm{0.026}$.
All the algorithms presented in this paper have been implemented in the free and open geochronological toolbox IsoplotR (Vermeesch, 2018). The ability to fit isochrons with models 1, 2, and 3 dates back to IsoplotR's first public release. Versions 2.1 and 5.2 added model 1 anchors to $\mathrm{U}/\mathrm{Pb}$ isochrons without and with uncertainty, respectively. Anchored York regression and model 3 anchors with uncertainty were added to IsoplotR 6.0. These functions can be accessed via the GUI or from the command line.
In the GUI, anchored regression is available from the “isochron” menu (so not from the “concordia” menu for $\mathrm{U}/\mathrm{Pb}$ data). The relevant settings can be selected from the “options” menu and should be selfexplanatory. Moving on to the R console, use a builtin ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ dataset to illustrate the commandline functionality and start with an unanchored inverse isochron fit using York regression.
library(IsoplotR) attach(examples) isochron(ArAr)
Carry out model 3a isochron regression anchored to a variable initial ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ ratio of 300±5 (1σ).
settings('iratio','Ar40Ar36',300,5) isochron(ArAr,anchor=1,model=3,taxis=TRUE)
Here, settings(…)
is used to change the nonradiogenic argon composition, anchor=1
fixes the isochron to the inherited component, model=3
tells IsoplotR to assign the dispersion to the reported uncertainty of the atmospheric ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ ratio, and taxis=TRUE
changes the x axis to a timescale (as in Fig. 2).
Finally, carry out model 3b isochron regression anchored to a radiogenic composition that was set at 62±1 Ma (1σ).
isochron(ArAr,anchor=c(2,62,1), model=3,taxis=TRUE)
The dispersion can be changed to a free parameter by setting the third element of the anchor
argument to a nonpositive value. Documentation for additional options and settings can be accessed from the console using the help(isochron)
or ?isochron
commands.
This paper builds on previous work by McIntyre et al. (1966), Titterington and Halliday (1979), and Ludwig (2012). McIntyre et al. (1966) introduced the concept of model 3 isochron regression, attributing the excess dispersion of errorchrons to variability of the inherited component. This approach is implemented in Isoplot (Ludwig, 2012) and is referred to as model 3a regression in this paper. Note that the McIntyre et al. (1966) definition of model 2 regression differs from that of Ludwig (2012). It is more akin (but not identical) to the $\sqrt{\text{MSWD}}$ approach of Sect. 3.
Titterington and Halliday (1979) recast the model 3 algorithm of McIntyre et al. (1966) in a maximum likelihood framework. This paper extends the maximum likelihood approach to a second type of geological scenario, in which excess dispersion of the data around the isochron is not attributed to the inherited component but to the radiogenic component. This approach is referred to as model 3b regression.
The maximum likelihood formulation of isochron regression can also be used to anchor isochrons to either the inherited or the radiogenic component and to assign geologically meaningful uncertainty to the anchor. In reality it is possible that some samples are affected by both mechanisms so that different aliquots of the same sample differ in their initial ratios as well as the timing of their isotopic closure. Unfortunately, it is not possible to simultaneously capture both types of dispersion using the algorithms of this paper.
The difference between model 1 and model 3 regression represents two contrasting views of geological reality. The model 1 approach assumes that the isotopic composition of minerals represents discrete components recorded at distinct events. In contrast, model 3 isochrons represent a “fuzzier” reality, in which the initial composition or timing of isotopic closure is allowed to vary within a rock. Under the latter model, the concept of an errorchron no longer makes sense.
P  a radioactive parent (e.g. ^{87}Rb) 
D  the radiogenic daughter of P (e.g. ^{87}Sr) 
d  a nonradiogenic sister isotope of D (e.g. ^{86}Sr) 
a, b, x_{i}, y_{i}  the true (but unknown) values of the intercept, slope, and isotopic ratio of the ith aliquot (out of n aliquots) from a sample 
$\widehat{a}$, $\widehat{b}$, ${\widehat{x}}_{i}$, ${\widehat{y}}_{i}$  the estimated values of a, b, x_{i}, and y_{i}, respectively 
X_{i}, Y_{i}  the measured values of x_{i} and y_{i}, respectively 
s[X_{i}], s[Y_{i}], s[X_{i},Y_{i}]  the standard errors and covariance of X_{i} and Y_{i} 
${\mathit{\u03f5}}_{x,i}={X}_{i}{x}_{i}$, ${\mathit{\u03f5}}_{y,i}={Y}_{i}{y}_{i}$  the residuals of the data around the bestfit line 
Σ_{i}  the 2×2 covariance matrix of X_{i} and Y_{i} 
x, y, X, Y, Σ  arrays and sets containing all n x_{i}, y_{i}, X_{i}, Y_{i}, and Σ_{i} values 
All uncertainties in this paper are reported as 95 % confidence intervals except if noted otherwise.
The synthetic data in Figs. 1 and 3 were generated as follows.

Generate x by drawing 10 random numbers from a logitnormal distribution where the mean and standard deviation of the logits are 0 and 1$/$5, respectively.

Generate y by plugging x into Eq. (3) with $a=b=\mathrm{1}$.

Turn the resulting 10 pairs of {x_{i},y_{i}} values into 10 trios of $\mathit{\{}{P}_{i},{D}_{i},{d}_{i}\mathit{\}}$ values where d_{i} is a random number between 100 and 500.

Use each of the $\mathit{\{}{P}_{i},{D}_{i},{d}_{i}\mathit{\}}$ values obtained in the previous step as parameters for a Poisson experiment.

Use the Poisson values obtained in the previous step to form 10 pairs of X and Y measurements and their (correlated) uncertainties.

Add some excess dispersion by shrinking the uncertainties by 20 %.
Note that this procedure fits the definition of a model 1 isochron, even though panels (b) and (d) of Fig. 3 use the resulting data to illustrate model 3 regression. The synthetic data in Figs. 2 and 4 are generated using a slightly different algorithm that simulates model 3a and 3b overdispersion. This modified procedure uses steps 1 and 3–4 of the procedure for Fig. 1 whilst replacing step 2 with a random effects model, in which a and b are drawn from random normal distributions with standard deviations corresponding to the dispersion parameters ${\mathit{\sigma}}_{\overline{a}}$ (10 % for Fig. 2a and 5 % for Fig. 4) and ${\mathit{\sigma}}_{\overline{b}}$ (10 % for Fig. 2b).
Whether a dataset is overdispersed with respect to the analytical uncertainties can be formally assessed using the chisquare statistic and test. To this end, calculate the χ^{2} statistic using the sum of squares term of Eq. (6).
If analytical uncertainty is the only source of dispersion in the dataset, then χ^{2} is expected to follow a chisquare distribution with n−2 degrees of freedom ($\text{DF}=n\mathrm{2}$, or, for anchored regression, $\text{DF}=n\mathrm{1}$). This is called the “null distribution”. The probability of observing a value greater than χ^{2} under this distribution is called the p value. A p value cutoff of 0.05 is generally used to distinguish isochrons (p≥0.05) from errorchrons (p<0.05). Dividing Eq. (C1) by the number of degrees of freedom produces a “reduced chisquare statistic”, which is known to geologists as an MSWD (McIntyre et al., 1966).
In the absence of overdispersion and for sufficiently large samples (n>20), the expected null distribution of the MSWD is approximately normal with a mean of 1 and a standard deviation $\sqrt{\mathrm{2}/\text{DF}}$ (Wendt and Carl, 1991). This information can be used as an alternative way to test the null hypothesis.
Section 4 describes how the loglikelihood functions of Eq. (13) can be maximized in a twostep process: (1) find the x_{i} values that maximize ℒℒ for any pair of a and b values, and (2) search the space of a and b values until the overall maximum is found. For model 1 regression, the first step has a direct solution. Some shorthand notation for the inverse of the covariance matrix Σ_{i} is introduced:
where ${\mathrm{\Omega}}_{\mathrm{1},\mathrm{2}}={\mathrm{\Omega}}_{\mathrm{2},\mathrm{1}}$. The fitted points x_{i} can be found by solving
which yields
The same approach can be used for model 3a regression (with $\partial \mathcal{L}{\mathcal{L}}_{\mathrm{3}a}/\partial {x}_{i}=\mathrm{0}$). However, it does not work for model 3b isochrons because ${\partial}^{\mathrm{2}}\mathcal{L}{\mathcal{L}}_{\mathrm{3}b}/\partial {x}_{i}\partial {x}_{j}\ne \mathrm{0}$ when i≠j. Therefore, the x_{i} values are interdependent and cannot be solved separately. Although it is possible to jointly optimize the entire x vector, this takes considerably longer than model 1 and model 3a regression. Fortunately, as explained in Sect. 5, model 3b regression can be reformulated in terms of the model 3a algorithm.
IsoplotR is available from CRAN (https://CRAN.Rproject.org/package=IsoplotR, last access: 15 July 2024) and Zenodo (https://doi.org/10.5281/zenodo.11114310; Vermeesch and Lloyd, 2024).
This paper only uses synthetic data, which were generated using the procedures outlined in Appendix B.
The author is a member of the editorial board of Geochronology. The peerreview process was guided by an independent editor, and the author also has no other competing interests to declare.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
The author would like to thank reviewers Donald Davis and John Rudge for suggesting to expand the discussion of model 3 regression and for checking the equations. This work was completed during a sabbatical stay at the Institute of Geology and Geophysics, Chinese Academy of Sciences, which was supported by the CAS President's International Fellowship Initiative (PIFI). The author would like to thank XianHua Li of IGGCAS and Yang Li of Peking University for hosting him at IGGCAS. He would also like to thank IsoplotR users Guilhem Hoareau (CNRS, France), Shitou Wu (IGGCAS, China), and Stijn Glorie (Adelaide, Australia) for discussing the utility of anchored isochrons and moving this feature up the priority list for IsoplotR.
This research has been supported by the Natural Environment Research Council (grant no. NE/T001518/1).
This paper was edited by Noah M McLean and reviewed by Donald Davis and John Rudge.
Daëron, M. and Vermeesch, P.: Omnivariant generalized least squares regression: Theory, geochronological applications, and making the case for reconciled Δ_{47} calibrations, Chem. Geol., 647, 121881, https://doi.org/10.1016/j.chemgeo.2023.121881, 2023. a
Li, Y. and Vermeesch, P.: Short communication: Inverse isochron regression for Re–Os, K–Ca and other chronometers, Geochronology, 3, 415–420, https://doi.org/10.5194/gchron34152021, 2021. a, b
Ludwig, K. R.: On the treatment of concordant uraniumlead ages, Geochim. Cosmochim. Ac., 62, 665–676, https://doi.org/10.1016/S00167037(98)000593, 1998. a
Ludwig, K. R.: User's manual for Isoplot 3.75: a geochronological toolkit for Microsoft Excel, Berkeley Geochronology Center, 5, https://doi.org/10.5281/zenodo.12744084, 2012. a, b, c
McIntyre, G. A., Brooks, C., Compston, W., and Turek, A.: The Statistical Assessment of RbSr Isochrons, J. Geophys. Res., 71, 5459–5468, 1966. a, b, c, d, e
Schaen, A. J., Jicha, B. R., Hodges, K. V., Vermeesch, P., Stelten, M. E., Mercer, C. M., Phillips, D., Rivera, T. A., Jourdan, F., Matchan, E. L., Hemming, S. R., Morgan, L. E., Kelley, S. P., Cassata, W. S., Heizler, M. T., Vasconcelos, P. M., Benowitz, J. A., Koppers, A. A. P., Mark, D. F., Niespolo, E. M., Sprain, C. J., Hames, W. E., Kuiper, K. F., Turrin, B. D., Renne, P. R., Ross, J., Nomade, S., Guillou, H., Webb, L. E., Cohen, B. A., Calvert, A. T., Joyce, N., Ganerød, M., Wijbrans, J., Ishizuka, O., He, H., Ramirez, A., Pfänder, J. A., LopezMartínez, M., Qiu, H., and Singer, B. S.: Interpreting and reporting ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ geochronologic data, Bulletin, 133, 461–487, https://doi.org/10.1130/B35560.1, 2021. a
Titterington, D. M. and Halliday, A. N.: On the fitting of parallel isochrons and the method of maximum likelihood, Chem. Geol., 26, 183–195, 1979. a, b, c, d
Vermeesch, P.: IsoplotR
: a free and open toolbox for geochronology, Geosci. Front., 9, 1479–1493, 2018. a
Vermeesch, P. and Lloyd, J. C.: pvermees/IsoplotR: 6.2 (6.2), Zenodo [code], https://doi.org/10.5281/zenodo.11114310, 2024. a
Wendt, I. and Carl, C.: The statistical distribution of the mean squared weighted deviation, Chem. Geol.: Isotope Geoscience Section, 86, 275–285, 1991. a
York, D., Evensen, N. M., Martínez, M. L., and De Basabe Delgado, J.: Unified equations for the slope, intercept, and standard errors of the best straight line, Am. J. Phys., 72, 367–375, 2004. a, b
 Abstract
 Introduction
 Maximum likelihood formulation of York regression
 Dealing with overdispersion
 Specific cases
 Flipped isochron regression
 Anchored isochron regression
 Accounting for anchor uncertainty
 $\mathrm{Pb}/\mathrm{Pb}$ isochrons
 $\mathrm{U}/\mathrm{Pb}$ isochrons
 Implementation in IsoplotR
 Conclusions
 Appendix A: Definitions
 Appendix B: Data
 Appendix C: Testing overdispersion
 Appendix D: Further details about the maximum likelihood solution
 Code availability
 Data availability
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Maximum likelihood formulation of York regression
 Dealing with overdispersion
 Specific cases
 Flipped isochron regression
 Anchored isochron regression
 Accounting for anchor uncertainty
 $\mathrm{Pb}/\mathrm{Pb}$ isochrons
 $\mathrm{U}/\mathrm{Pb}$ isochrons
 Implementation in IsoplotR
 Conclusions
 Appendix A: Definitions
 Appendix B: Data
 Appendix C: Testing overdispersion
 Appendix D: Further details about the maximum likelihood solution
 Code availability
 Data availability
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References