Statistical bias in isotope ratios

Abstract

This paper presents the mathematics of the systematic bias in the expected value of the ratio of two noise-corrected Poisson-distributed variables, such as ion counting measurements. Such bias can lead to the reporting of incorrect ratios and, in some cases, systematic correlations with other measurements which can impact the scientific interpretation. We describe a novel method of treating such measurements which results in a negligible, exponentially small bias. We also re-examine the conventional approach deriving an exact expression for the bias including the noise correction explicitly.

---

1 Introduction

Analysts measuring a time series of ion beam intensities are faced with choices of how to process the data in order to estimate intensity ratios. Is it preferable, for example, to calculate the mean of the ratios, the ratio of the means or is there a better approach entirely? At the heart of the problem is the fact that the mean value a ratio estimator returns, over the long term, is biased relative to the true ratio being estimated. A simple example serves to illustrate the phenomenon of bias as follows. Suppose we write a 1 and a 3 on the faces of coin A and a 3 and 5 on the faces of coin B. Clearly, the mean value of a coin flip is 2 and 4 for coins A and B respectively. We now flip the two coins and record the ratio, B/A, as indicated by the upturned faces. In such an experiment, there are four equally probable outcomes: 3/1, 3/3, 5/1, and 5/3. Hence, the expectation value of the ratio, E(B/A) = (3 + 1 + 5 + 5/3)/4 or 8/3. The expectation value of B/A is the mean value of B/A over the long term. If we had hoped to estimate the ‘true’ ratio equal to the ratio of the means for each coin, 4/2 = 2, then clearly the estimator B/A is biased. We define the bias as the relative difference between E(B/A) and the true ratio. In this example, therefore, the bias is (8/3–2)/2 = 1/3.

The situation with ratios of ion-counting signals in isotope ratio measurements is directly analogous. The number of ions counted is accurately modelled by Poisson statistics which states that the probability of detecting X ions is given by

Pois(X;μ_x) = e^−μ_xμ^X_x/X!, (1)

where μ_x is the population mean or expectation value of X. Note that this single parameter, μ_x, characterises the Poisson distribution. An analyst may wish to calculate an estimate of μ_x/μ_y, the ratio of the expectation values of two such signals, X and Y. More generally, it is often the case that analyte signals are superposed on ‘noise’, whose origin may be in the detector (dark noise), spectral or contamination during sample preparation. Regardless of origin the mean noise is often assumed to be constant, is determined empirically and the measured signals corrected for its contribution. Let the mean noise contribution to X and Y be μ_x0 and μ_y0 respectively. The ratio to be estimated, R, is, therefore, given by

In principal, any function, f(X,Y), may have it's bias relative to R determined by deriving the expectation valve, E(f), given by the double summation over all X and Y of f(X,Y)·Pois(X;μ_x)·Pois(Y;μ_y). For example, the commonly used ratio estimator,

is considered in Section 2 but, even without detailed analysis, r₀ is clearly problematic for any integer value of μ_y0, including zero, as E(r₀) does not exist due to the non-zero probability of events where the denominator, Y − μ_y0, is zero. In Section 3 a novel, quasi-unbiased ratio estimator is proposed, which is well behaved for all Y and μ_y0.

1.1. Previous work

Coakley et al.¹ and Ogliore et al.² have shown that X/Y is a biased estimate of R. Similar approaches by both^1,2 yield

E(X/Y) ≈ μ_x/μ_y(1 + 1/μ_y + 2/μ²_y) (4)

in the limit of large μ_y. As we shall show in Section 2, this approximation holds providing Y = 0 events are excluded from the distribution. Noise correction is not considered explicitly but, since the method used by these authors is rather general, it is straightforward to include it as a modification of the distribution. Note that eqn (4) is not a solution to the bias problem, it is only a means of estimating its magnitude. If it can be assumed that μ_y is constant then μ_y, and hence the bias, may be estimated and the ratio corrected but such a procedure lacks generality.

Ogliore et al.² compare ratio estimators r₁ = x̄/ȳ (the ‘ratio of the means’) and (the ‘mean of the ratios’) where x̄ = n⁻¹∑ⁿ_i=1x_i, and ȳ = n⁻¹∑ⁿ_i=1y_i are the means of samples x_i and y_i, and by deriving expressions for E(r₁) and E(r₂). Note, however, that these reduce to a single problem as follows. Let x_i and y_i be instances of independent random variables X and Y respectively where X has a Poisson distribution with mean μ_x, which we shall write X ∼ Pois(μ_x), and similarly Y ∼ Pois(μ_y). Since the sum of Poisson variables is itself a Poisson variable, we have nx̄∼ Pois(nμ_x) and similarly for nȳ. With this nomenclature, E(r₂) = E(X/Y) and is given by approximation 4 above and E(r₁) has the same form with μ_x and μ_y replaced by nμ_x and nμ_y respectively (see Ogliore et al.²eqn (19) and (22)). The biases in r₂ and r₁ are, therefore, O(μ_y⁻¹) and O{(nμ_y)⁻¹} respectively. Note that both r₁ and r₂ have, to first order, biases proportional to the reciprocal of the mean number of counts in the denominator when the ratio is taken.

Ogliore et al.² also consider Beale's ratio estimator, r₃, given by

where cov and var return the sample covariance and variance respectively. Beale's estimator reduces the bias to O{(nμ_y)⁻²} but, in common with r₁, all n data are reduced to a single value so any true within-analysis variations in R, which may be of interest, are obscured.

2 The ratio of noise-corrected Poisson variables

We will assume that the distribution of the ion counts, X and Y, obeys Poisson statistics. Furthermore, we will consider only cases where X and Y are independent. This latter restriction is not so severe as it might seem since by far the most important correlated variations in ion counting signals are common mode ‘intensity’ fluctuations, that is, proportional changes in both μ_x and μ_y which largely cancel by taking the ratio. Of course, common mode variations will not be completely rejected on account of any intensity dependence of the bias.

Let R be defined as given by eqn (2) and let Z be distributed like Y but truncated at y₀, that is, the probability distribution function, Pr(Z), is given by

where N is a normalisation constant and y₀ is an integer. To satisfy ∑ Pr(Z) = 1 the normalisation constant required is
where γ is the incomplete gamma function and P the normalised incomplete gamma function.³ We shall choose y₀ to be sufficiently large to ensure Z − μ_y0 cannot be zero or negative, i.e.
. Note that, in cases of large signal to noise ratio, the probability of Y ≤ μ_y0 may be so small as to make excluding these events notional in practice. We shall now derive an expression for the expectation value ofwhich is the conventional expression for the noise-corrected ratio but with rejections to avoid zeroes or negative values in the denominator.

For independent random variables we can separate them thus

reducing the problem to one in a single variable, Z, so we can drop the subscript y, i.e. μ_y0 → μ₀ and μ_y → μ.

The Taylor series expansion of 1/(Z − μ₀) about μ₀ = 0 is,

Note that the truncation of the distribution ensures ∣μ₀/Z∣ < 1, guaranteeing convergence. To take the expectation value of the right-hand side (RHS) of eqn (8) requires an expression for the expectation value of 1/Z^k. This problem is addressed in the appendix and given by eqn (25). Substituting into eqn (8) gives,

where coefficients a_k(j) are given by eqn (26) (or, more conveniently for computational purposes, using a₁(j) = (j − 1)!, a_j(j) = 1 and the recursion 11). Changing the order of the summation,

The first few values of a_k(j) are,

j =	1	2	3	4	5
a1(j)	1	1	2	6	24
a2(j)		1	3	11	50
a3(j)			1	6	35
a4(j)				1	10

From the exact form eqn (9), we can derive an asymptotic form in the limit as μ → ∞ and bounded μ₀ by replacing the infinite sum with a finite sum, and using P(y₀ + j + 1,μ) → 1, which is true so long as y₀ and j are also bounded. Hence,

Define the bias, B(μ,μ₀), as the relative difference between E(1/(Z − μ₀)) and 1/(μ − μ₀), i.e.

B(μ,μ₀) = (μ − μ₀)E(1/(Z − μ₀)) − 1.

Substituting from eqn (10) and after some cancellation we have

From eqn (27) it follows that

a_k+1(j + 1) − a_k(j) = ja_k+1(j), (11)

hence
For example, putting n = 3 we have

B = 1/μ + (2 + 2μ₀)/μ² + O(1/μ³) (12)

in agreement with Coakley et al.¹ and Ogliore et al.² for the case μ₀ = 0.

3 Quasi-unbiased ratios

Here we present an alternative ratio estimator which reduces the bias to a factor which is exponentially small or quasi-unbiased. Furthermore, the method does not require truncation of the distribution so all the measured data can be used.

Let Y be distributed as before. We define a new random variable,

Y′(Y,μ₀) = (Y + 1)/M(1,Y + 2,μ₀) (13)

where M(a,b,z) is the Kummer confluent hypergeometric function. The series expansion of M is⁴where (a)_k is the rising factorial or Pochhammer's symbol,
Explicitly,

E(1/Y′) = (μ − μ₀)⁻¹(1 − e^{−μ + μ₀}).

Proof. From eqn (13) the expectation value is given by
A summation of this form is given by Prudnikov et al.⁵ which is
Substituting a → 1 and b → 2 gives,
where we have used Gradshteyn and Ryzhik⁶ equation 9.212 and M(1,1,z) = e^z.

Let B′ be the bias of E(1/Y′) relative to (μ − μ₀)⁻¹,

Therefore, noise-corrected quasi-unbiased ratios, r′, can be computed from measurements of Poisson events X and Y with,

Using eqn (15) a ratio can be calculated at each cycle of data collection and, if required, statistics on these such as the mean and standard deviation. This contrasts using either the ratio of the means, r₁, or Beale's estimator, r₃, both of which return a single ratio from a set of n measurements.

Fig. 1 compares conventional and quasi-unbiased ratios from 128 cycles of measurements on boron isotopes ¹⁰B and ¹¹B. The ion counts have been summed in ‘blocks’ of p cycles, the ratio taken for each block and the mean over all blocks plotted. The value of p, therefore, controls integration period and the total number of blocks, q = 128/p. Let the total counts in each block be x_i and y_i for the two isotopes, where i = 1,2,…q. We use p = 1,2,4,8 and 16 and plot the means of r′(x_i,y_i)(t_y/t_x) (red) and (x_i/y_i)(t_y/t_x) (blue) against mean of 1/y_i, where t_x and t_y are the cycle integration times for the two isotopes. The ‘ratio of the means’, equivalent to putting p = 128, is also shown and is indistinguishable from Beale's estimator. The noise is negligible for these data and has been set to zero. Both ¹¹B/¹⁰B and its reciprocal are plotted to demonstrate that the effect is independent of which isotope is chosen as the denominator. The plots show that, to first order and, as predicted by theory, the conventionally processed data, shown in blue, plot as a straight line with equal slope and intercept. The y-intercept, which corresponds to counts → ∞, should equal the unbiased ratio. Our novel quasi-unbiased ratio estimator (eqn (15)), shown by the red data have, as expected, no bias regardless of the number of counts.

Fig. 1 Boron isotope ratios (a) ¹¹B/¹⁰B and (b) ¹⁰B/¹¹B as a function of the mean number of counts per integration period showing the conventionally computed ratios in blue and the novel quasi-unbiased ratios (eqn (15)) in red. Ratios are calculated for each integration period in the analysis and the mean value plotted. Note that the quasi-unbiased ratios (red) show no trend and give the desired result regardless of the integration period. All plotted points are computed from the same 128 cycles of data by summing the counts in blocks of p adjacent cycles for p = 1,2,4,8, and 16, dividing the analysis into 128/p integration periods. The greatest bias in the blue data corresponds to p = 1 which plots at the far right. The ‘ratio of the means’ is also shown (open blue symbol), which is equivalent to p = 128, close to Beale's estimator (not shown) which lies below the open symbol in both (a) and (b) by 0.0056% and 0.0035% respectively. See main text for further details. Data are raw secondary-ion mass spectrometry (SIMS) data from an analysis of a foraminifera using a CAMECA IMS 1270.

Fig. 2 shows the results of a Monte-Carlo simulation with μ_x = μ_y = 20 and noise from 0 to 15. The simulation shows the noise-corrected ratio calculated in four ways: (i) conventionally without data rejection (eqn (3)), (ii) conventionally with rejection when the denominator is zero or negative (eqn (7)), (iii) Beale's estimator, and (iv) quasi-unbiased (eqn (15)) all with μ_x0 = μ_y0 = b. The mean of 10⁶ ratio estimations is plotted for each value of the noise, which is incremented in steps of 0.01. The simulated quasi-unbiased ratios are closer to unity (no bias) than either of the conventional ratios or Beale's estimator over the entire range of b, although significant scatter does occur for large b (see caption). Beale's estimator (eqn (5)) on the noise-corrected data requires, for each of the 10⁶ simulated ratio estimations, the variance, covariance and mean intensities to be calculated. We calculate these statistics from a dataset of n simulated data pairs, {x_i,y_i; i = 1…n}, and to make the comparison with the other estimators fair we draw these from a Poisson distribution with mean value 20/n, from which noise of b/n is subtracted. Beale's estimator is now

where

x′_i = x_i − b/n (17)

and similarly for y′_i. A value of n = 10 is used for the simulation of Fig. 2. A simplified variation on Beale's estimator for Poisson data, r′₃, which is independent of n, may be obtained by substituting n → 1, cov(x_i,y_i) → 0, var(y_i) → Y, x′_i → X − b and y′_i → Y − b giving

Fig. 2 Monte-Carlo simulation of (X − b)/(Y − b) (blue), (X − b)/(Z − b) (green), Beale's estimator (grey, see main text), simplified Beale's estimator r′₃ (cyan, eqn (18)) and the novel quasi-unbiased ratio, r′ (red, eqn (15)) with μ_x0 = μ_y0 = b, where X and Y are independent Poisson variables with means, μ_x and μ_y, equal to 20 and Z = Y for Y > b and rejected otherwise. For each data point the simulation computes the mean over 10⁶ samples. Rare single events where Y is small can significantly shift the mean r′ for large values of b giving rise to the observed scatter in the red data for b ≳ 10. Black lines show the theoretical behaviour in the cases of the green and red data. Much of the blue data are obscured behind the green. The noise, b, is incremented in steps of 0.01.

These substitutions are arrived at by noting that, for independent Poisson distributions, E{cov(X,Y)} = 0 and E{var(Y)} = E(Y). Simulations of both r₃ and r′₃ are included in Fig. 2. Note that for zero noise, r′₃ = X/(Y + 1), i.e. identical to r′ (eqn (15)), and for b > 0 the (absolute) bias of r′₃ is marginally greater (more negative) than that of r₃.

Fig. 3 compares the variance of r′ (eqn (15)) with that of r (eqn (7)) as a function of mean signal, μ_x = μ_y = μ, for three different noise to signal ratios. Over the plotted parameter range the variance of r′ is always the smaller (var(r′)/var(r) < 1) and, therefore, more efficient ratio estimator.

Fig. 3 Variance of quasi-unbiased ratio, r′ (eqn (15)), divided by the variance of r = (X − b)/(Z − b) as a function of mean signal, μ_x = μ_y = μ, for relative noise, b/μ = 0, 5% and 10% showing that the variance in r′ is always smaller than that of the conventional ratio, r, over this range of parameters.

4 Discussion

Quasi-unbiased ratios offer advantages over conventional methods of calculating noise-corrected ratios of ion-counting measurements, namely, (i) an exponentially small statistical bias, (ii) no need to sacrifice within-analysis information by summing counts over entire analysis before taking the ratio, (iii) insensitivity to common-mode changes in signal intensity, (iv) no mathematical singularities, and (v) good stability even with low signal to noise ratios (≳2).

Ratio bias has increasing importance in isotope ratio measurements since the scientific demands lead researchers to strive for ever higher precisions on small quantities of sample. This is particularly the case in secondary-ion mass spectrometry (SIMS) where, because of low blanks, very low count rates are acceptable. In studies on short-lived radionuclides (SLRs) (e.g. see ref. 7–10) the bias is particularly insidious as it can give rise to an apparent linear relationship between measurements of the daughter nuclide and a proxy for the parent nuclide on a so-called ‘isochron’ plot. To illustrate, consider the SLR ⁶⁰Fe which decays to ⁶⁰Ni. A suite of measurements are made on phases with a range of Fe/Ni ratios and an isochron plot is made of the daughter, ⁶⁰Ni, against the parent element, Fe. Both are plotted as ratios using some stable isotope of Ni as the denominator. A straight line with positive slope (the isochron) indicates that ⁶⁰Fe was present at the time the Fe and Ni were fractionated between the phases. It is sometimes the case that most or all of the variation in Fe/Ni is controlled by the Ni content which, since it appears in the denominator on both axes, will give rise to a linear relationship in the data as a consequence of the bias, adding to the positive slope due to in-growth of ⁶⁰Ni. Huss et al.¹¹ have recently reported this problem with some of their own published data on the Fe–Ni system concluding that in some, but not all, samples their published estimates of initial ⁶⁰Fe can no longer be distinguished from zero. They also discuss the likely size of any corrections to other published work on¹⁰Be, ²⁶Al and ⁵³Mn concluding that any changes to the published conclusions are small except for one older study on the Mn–Cr system in pallasites.

Note that, where internally normalised ratios are calculated, i.e. where a third isotope is used to correct for mass-bias, the magnitude of the statistical bias can increase or decrease, depending on the relative mass differences, by propagation of the bias of the normalising ratio. For example and assuming a linear mass bias law, in the case of ⁶⁰Ni/⁶¹Ni normalised to ⁶²Ni/⁶¹Ni the statistical bias on ⁶⁰Ni/⁶¹Ni increases by a factor of two, whereas using ⁶²Ni as the denominator the statistical bias on ⁶⁰Ni/⁶²Ni stays the same magnitude but changes sign as it does in the case of ²⁶Mg/²⁴Mg normalised to ²⁵Mg/²⁴Mg. For ⁵³Cr/⁵²Cr normalised to ⁵⁰Cr/⁵²Cr the statistical bias increases by a factor of 1.5.

Studies where mass-bias correction is made by sample – standard bracketing are potentially susceptible to statistical bias in cases where there are differences in analyte concentration (ion count rate) between sample and standard. Standards are usually chosen to have analyte concentrations high enough that good precision can be achieved in a short time under the same analytical conditions employed on the sample. Where ion counts are higher on standards than samples and ratios are calculated conventionally, the statistical bias will result in systematically high ratios reported on samples corrected by sample-standard bracketing. In short, sample-standard bracketing does not necessarily eliminate the bias.

Huss et al.¹¹ rightly point out that ratio bias is a problem that the community will have to be aware of to avoid this source of systematic error in future work. However, we disagree that the best solution is necessarily to sum the counts over the entire analysis before taking the ratio (with or without using Beale's ratio estimator), or to correct for the bias based upon eqn (4) or (12), for the reasons (i)–(v) given at the beginning of this discussion, but rather to use eqn (15) to compute the ratio r′ at each measurement cycle.

It may seem laborious to have to evaluate the Kummer confluent hypergeometric function for every measurement cycle but this should not be particularly so if (i) a good library of special functions is available to the software developers or (ii) the signal to noise ratio is sufficiently high to be able to truncate the infinite series of eqn (14) to yield an approximate value for M(1,Y + 2,μ₀). The truncation error, ε_n, using an upper summation limit of n − 1 in eqn (14) is given by,

Let α = μ₀/(Y + 2) be subject to the constraint 0 ≤ α < 1. This constraint is not severe: it is sufficient only that the signal is at least as large as the mean noise (and both are positive). Since (Y + 2)_k ≥ (Y + 2)^k we may write,

Making n the subject of the inequality,

Therefore, if we wish to calculate M(1,Y + 2,μ₀) with a truncation error no larger than ν we can use

ande.g. v = 10⁻⁴, Y = 10 and μ₀ = 0.5 gives n = 3 and truncation error, ε_n = 6 × 10⁻⁵, smaller than ν as required.

Appendix

A.1 The expectation value of 1/Z^l

For Z distributed as a truncated Poisson distribution (eqn (6)) we havewhere,

P(a,z) = γ(a,z)/(a − 1)! a = 1,2,…

is the normalised incomplete gamma function and γ is the incomplete gamma function.³ The coefficients a_l(j) are given by l − 1 nested summationsor, alternatively, by a_l(j) = j!b_l(j) and the recursionThe proof of eqn (25) can be subdivided into proofs ofandwhere we have used Pochhammer's symbol, (Z + 1)_j = (Z + 1)(Z + 2)⋯(Z + j). The proof of eqn (29) yields easily as follows. From the probability distribution function defined in eqn (6) it follows that
where we have used Arfken¹² equation 10.70 in the final step completing the proof.

Eqn (28) yields as follows. Denote

{p}_j = 1/(p)_j.

Given integers s ≥ 0 and t ≥ s we have

Proof: for s = t, obviously RHS = LHS = 0. For t > s, by induction on t − s using

Definethat is, l nested summations where each lower limit is the next outer summation variable. We have explicitly

h_l(j) = 1/n^l{j + 1}_n (32)

Proof: by induction on l. From the definition (eqn (31)) follows the recursion

The result is trivial for l = 0; for l > 0, substituting the RHS value for h_l−1(k) (eqn (32)) into the recursion (eqn (33)),

using eqn (30) with t → ∞, s → j, and j → n + 1.Q.E.D.

Changing the order of summation and combining eqn (31) and (32) gives

with j → 0 and substituting using the identity ,
With n → Z this completes the proof of eqn (28) and, hence, also of eqn (25).

Acknowledgements

The SIMS data in Fig. 1 were collected at the Edinburgh Ion Microprobe NERC Facility on a sample provided by Dr Laura Foster of the University of Bristol. Dr Coath is supported in this work by STFC grant ref. ST/F002734/1.

References

1. K. Coakley, D. Simons and A. Leifer, Int. J. Mass Spectrom., 2005, 240, 107–120.
2. R. C. Ogliore, G. R. Huss and K. Nagashima, Nucl. Instrum. Methods Phys. Res., Sect. B, 2011, 269, 1910–1918.
3. Digital Library of Mathematical Functions, National Institute of Standards and Technology, 2010, ch. 8.
4. Digital Library of Mathematical Functions, National Institute of Standards and Technology, 2010, ch. 13.
5. A. Prudnikov, Y. Brychkov and O. Marichev, in More Special Functions, Gordon and Breach, 1986, vol. 3, p. 409.
6. I. Gradshteyn and I. Ryzhik, in Tables of Integrals, Series and Products, Academic Press, 7th edn, 2007, p. 1023.
7. K. McKeegan, M. Chaussidon and F. Robert, Science, 2000, 289, 1334–1337.
8. E. Kurahashi, N. T. Kita, H. Nagahara and Y. Morishita, Geochim. Cosmochim. Acta, 2008, 72, 3865–3882.
9. P. Hoppe, D. Macdougall and G. W. Lugmair, Meteorit. Planet. Sci., 2007, 42, 1309–1320.
10. S. Tachibana and G. Huss, Astrophys. J., 2003, 588, L41–L44.
11. G. Huss, R. Ogliore, K. Nagashima, M. Telus and C. Jilly, 42nd Lunar and Planetary Science Conference, #2608, 2011.
12. G. Arfken, in Mathematical Methods for Physicists, Academic Press, Inc., 3rd edn, 1985, p. 566.

---

Footnote

† Present address: University of California, Los Angeles, Department of Earth and Space Sciences, 595 Charles Young Dr E., Los Angeles, CA 90095-1567, USA.

---