INTRODUCTION
The CRONUS-Earth and CRONUS-EU projects (Phillips et al. Reference Phillips, Argento, Balco, Caffee, Clem, Dunai, Finkel, Goehring, Gosse, Hudson, Jull, Kelly, Kurz, Lal, Lifton, Marrero, Nishiizumi, Reedy, Schaefer, Stone, Swanson and Zreda2016) constituted a coordinated effort among many scientists to address basic issues, such as accuracy and reliability, which arise in the use of cosmogenic nuclides in earth science research. One facet of this project involved an intercomparison study of various geologic reference materials. For each reference material, aliquots were distributed to several laboratories; at each lab, the concentration of a particular cosmogenic nuclide was measured. Besides serving as a check on how well the different labs agreed, the exercise also provided the opportunity to combine the results from the different labs to calculate a consensus value for the cosmogenic nuclide concentration of a given reference material. In this paper we discuss issues that arise when computing a consensus value, especially when overdisperison is present.
The results for many of the CRONUS samples suffered from overdispersion; this means that the measurements from the various labs differed by more than we would expect, given the quoted uncertainties. See Phillips et al. (Reference Phillips, Argento, Balco, Caffee, Clem, Dunai, Finkel, Goehring, Gosse, Hudson, Jull, Kelly, Kurz, Lal, Lifton, Marrero, Nishiizumi, Reedy, Schaefer, Stone, Swanson and Zreda2016) for a full discussion of the overdispersion issues in the CRONUS samples.
When faced with an overdispersion problem, there are several, quite different ways to proceed. One family of approaches involves increasing the variance of each result (Galbraith and Laslett Reference Galbraith and Laslett1993). One version of this adds an extra amount to the variance of each result; the size of the added variance is adjusted to fit the spread in the data. Vermeesch et al. (Reference Vermeesch, Balco, Blard, Dunai, Kober, Niedermann, Shuster, Strasky, Stuart, Wieler and Zimmermann2015), Blard et al. (Reference Blard, Balco, Burnard, Farley, Fenton, Friedrich, Jull, Niedermann, Pik, Schaefer, Scott, Shuster, Stuart, Stute, Tibari, Winckler and Zimmermann2015), and Phillips et al. (Reference Phillips, Argento, Balco, Caffee, Clem, Dunai, Finkel, Goehring, Gosse, Hudson, Jull, Kelly, Kurz, Lal, Lifton, Marrero, Nishiizumi, Reedy, Schaefer, Stone, Swanson and Zreda2016) used this method, which we will refer to as the added variance method. A different version of this approach multiplies the variance of each result by a common factor, in order to fit the spread in the data. Examples of this variance multiplier approach are discussed in Christen and Perez (Reference Christen and Perez2009) and Dose (Reference Dose2003).
An alternate approach uses the data to assess the reliability of each lab’s results, and then to give less weight to the results that differ too much from the others. There are many versions of this approach, of varying degrees of computational complexity, but many have a Bayesian character (e.g. Christen Reference Christen1994; Dose Reference Dose2003; Hanson Reference Hanson2005; Blaauw et al. Reference Blaauw, Bakker, Christen, Hall and van der Plicht2007; Bronk Ramsey Reference Bronk Ramsey2009). For general background on the Bayesian approach, see D’Agostini (Reference D’Agostini2003) and Dose (Reference Dose2003). Our goal in this paper is not to give a lengthy review of these approaches. Rather, we focus on one representative such approach, and compare its output to that of the more traditional, added variance approach. In this paper we review, and then apply, an approach to this problem presented in the paper Press (Reference Press1997).
The Press method, which is directly based on Bayesian thinking, has many attractive features: it automatically gives less weight to outliers, and it produces a realistic probability distribution for the quantity of interest. For us, one virtue of the Press approach is its computational simplicity relative to other versions of the Bayesian approach. The calculations presented in this paper required numerical evaluations of double integrals, which were quite straightforward using, for example Maple or MATLAB. Other Bayesian approaches, some of which are designed to handle more complicated situations, can require more extensive calculations. After explaining the Press method, we apply it to the CRONUS cosmogenic nuclide data sets. The Press approach does not seem to be widely known in the geochronology community, so we hope that this paper can serve as a useful introduction.
The question naturally arises as to which approach for combining the results is better, the Press method or the added variance method? We must realize that any model makes assumptions, either explicit or implicit, about what has gone wrong. Thus, the answer to our question depends on the source of the discrepancies in the results. If the overdispersion is due to all the labs underestimating their uncertainties by more or less the same amount, then the added variance approach makes sense. However, if the overdispersion is caused by one or two labs with inadequate results, while the other labs have accurate results, then the Press method is a better representation of the situation. Thus neither method is better in any absolute sense. In the final analysis, the scientist must use judgement to decide which method is suitable.
It is also worth noting that we can often gain insight into the data by applying both methods and seeing how the results differ. An example of this is given later in the section “CRONUS-P 3He Results.”
One way to view the Press methodology, as well other versions of this approach, is as follows. In the examples with overdispersion, we have a set of N results, and there is no a priori reason to exclude any particular data point. On the other hand, there is clearly a significant level of inconsistency among the various results. Thus, simply using a weighted average to determine the best estimate is not the correct approach. The Press method allows for the fact that each data point might be incorrect, with a probability that is determined self-consistently from the data. The method:
-
∙ produces a probability distribution for the value of the physical quantity that reflects the dispersion in the data, and so does not give an artificially small standard deviation to our estimate;
-
∙ automatically downgrades the influence of an outlier that quotes too small an uncertainty; and
-
∙ allows us to calculate the probability that a given data point is unsatisfactory.
The plan of this paper is as follows. In the next section, we review the Press approach. In a later section we then apply it to the CRONUS samples, some of which exhibit overdispersion. To show how the Press method works, we first give detailed discussions of two samples, the 3He results on CRONUS-P and the 14C measurements in CRONUS-A. Our results for the other CRONUS samples are then given in a more succinct fashion. The paper concludes with a discussion and summary.
THE PRESS APPROACH
Derivation of the Basic Equations
In this section we review the Press method (Press Reference Press1997) as applied to a set of cosmogenic nuclide measurements. We assume that we have results from N laboratories. For each lab, the output is a probability density function (pdf) for the concentration (x) of the nuclide being measured. Using this pdf we may compute, for example, x i and σ i , the average and standard deviation for lab i.
Following Press, we then assume that, for all the labs, there is a probability β that its reported result is correct, and a probability (1−β) that it is incorrect. The same β is used for all the labs. At this point we do not know the value of β. We use the data to generate a joint pdf for X and β, where X is the (unknown) true value of the nuclide concentration.
The key equation we now use is a version of Bayes’ theorem
$$f(X,\beta \!\mid\!results)f(results){\equals}f(results\!\mid\!X,\beta )f_{0} (X,\beta )$$
Here, f (X,β|results) is what we want, the joint pdf for X and β, given the results from the N labs. For our purposes, we take the results to be the set of values x i ,σ i . The pdf f(results) is the probability of obtaining the results, and may be omitted if we properly normalize our answer
$$f(X,\beta \!\mid\!results){\equals}cf(results\!\mid\!X,\beta )f_{0} (X,\beta )$$
Here c is chosen to enforce normalization
$${\int}_{\!\!\!0}^\infty dX{\int}_{\!\!\!0}^1 d\beta f(X,\beta \!\mid\!results){\equals}1$$
The crucial input comes in modeling the pdf f (results|X, β), which is the probability density for obtaining the results, given X and β. This pdf is known as the likelihood. Again following Press, we write it as follows:
$$f(results\!\mid\!X,\beta ){\equals}\prod\limits_{i{\equals}1}^N (\beta G_{{\sigma _{i} }} (x_{i} {\minus}X){\plus}(1{\minus}\beta )G_{s} (x_{i} {\minus}X))$$
We are using the following notation for the Gaussians:
$$G_{\sigma } (x){\equals}(2\pi \sigma ^{2} )^{{{\minus}1\!/2}} \exp ({\minus}{{x^{2} } \over {2\sigma ^{2} }})$$
To interpret Equation (4), we reason as follows. If the result from lab i is satisfactory (which has probability β), then the probabilty of measuring x i is a Gaussian centered on the true value X, with the quoted uncertainty σ i . If the result from lab i is unsatisfactory, then the pdf for x i is taken to be a very broad, uninformative Gaussian centered on X, with a large standard deviation given by s. Our choice for s will be discussed in the next section. Thus, the likelihood builds in the possibility that any given lab may have a result which is much poorer than expected.
Finally the pdf ƒ0 (X, β) is called the prior for this problem. It embodies the probabilities we assign to X and to β before we acquired the data. For example, if we have great confidence in the measurements done at all the labs, we might choose a prior in which most of the statistical weight was near β=1. Conversely, if we are not so confident, we may choose a prior which distributes the weight more evenly from β=0 to β=1. For our applications in this paper, we will take the prior to be a constant; so, we build in no prior knowledge as to the value of X, and make no strong judgement about the accuracy of the data.
Once we have used Equations (2) and (4) to get ƒ (X, β|results), we obtain the pdf for X by integrating out β
$$f(X\!\mid\!results){\equals}{\int}_{\!\!\!0}^1 d\beta f(X,\beta \!\mid\!results)$$
So the results for X include an average over possible values of β. As already mentioned, we may compute the probability that the result of any particular lab is correct; this is an important way to gain insight into the results. Note that this probability is not given simply by <β>, and in general will be different for different labs.
To calculate the probability that a given lab’s result is satisfactory, we proceed as follows (Press, Reference Press1997). Let p i denote the probability that lab i is satisfactory. Its useful to define the following quantities:
$$H_{i} {\equals}f_{0} (X,\beta )\prod\limits_{j\,\ne\,i} (\beta G_{{\sigma _{j} }} (x_{j} {\minus}X){\plus}(1{\minus}\beta )G_{s} (x_{j} {\minus}X))$$
$$T_{i} {\equals}{\int}_{\!\!\!0}^1 d\beta {\int}_{\!\!\!0}^\infty dXH_{i} \beta G_{{\sigma _{i} }} (x_{i} {\minus}X)$$
$$B_{i} {\equals}{\int}_{\!\!\!0}^1 d\beta {\int}_{\!\!\!0}^\infty dXH_{i} (1{\minus}\beta )G_{s} (x_{i} {\minus}X)$$
Then the result for p i is given by
$${{p_{i} } \over {1{\minus}p_{i} }}{\equals}{{T_{i} } \over {B_{i} }}$$
We should emphasize that we cannot claim, in any absolute sense, that the result of a given lab is correct or incorrect, since we do not know the true value X. Rather, the probability that a given result is correct reflects the extent to which that result is consistent with the measurements of the other labs.
One technical point should be mentioned. The Gaussians appearing in the likelihood (4) should each be normalized on x i , for fixed X; since x i is a nonnegative quantity, this means that the integral over x i , from zero to infinity, should be equal to one. If the width of the Gaussian (i.e. either s or σ i ) is much less than X, the factor (2πσ 2)–1/2 in Equation (5) will be accurate; however, if s is chosen to be large, there may be cases where the more careful normalization should be used. This is given by
$$G_{s} {\equals}{1 \over \gamma }(2\pi s^{2} )^{{{\minus}1\!/2}} \exp ({\minus}{{(x_{i} {\minus}X)^{2} } \over {2s^{2} }})$$
$$\gamma {\equals}1{\minus}{1 \over 2}er\,fc({X \over {s\sqrt 2 }})$$
APPLICATION TO CRONUS SAMPLES
CRONUS-P 3He Results
For our first example we will analyze the 3He concentration of the CRONUS-P pyroxene sample. This sample was measured at six different laboratories, and the results were reported and analyzed by Blard et al. (Reference Blard, Balco, Burnard, Farley, Fenton, Friedrich, Jull, Niedermann, Pik, Schaefer, Scott, Shuster, Stuart, Stute, Tibari, Winckler and Zimmermann2015). The results from the six labs differ by more than we would expect, given the uncertainties quoted by the labs. Blard et al. refer to this as overdispersion, and give a good discussion of the important issues. The data for the 3He concentration of CRONUS-P are given in Table 3 of Blard et al. (Reference Blard, Balco, Burnard, Farley, Fenton, Friedrich, Jull, Niedermann, Pik, Schaefer, Scott, Shuster, Stuart, Stute, Tibari, Winckler and Zimmermann2015). We see that the concentrations reported by the six laboratories range from 4.86 to 5.25 in units of 109 atoms per gram. (We will use these units for X, s, and σ X in our CRONUS-P discussions.) Since the quoted uncertainties range from 0.02 to 0.07, there is clearly too big a spread in the results.
When we apply the Press method, we need to specify the value of s, the standard deviation of the broad Gaussian in Equation (4). For the method to give sensible, convincing results, we need to think carefully about this value. First, the broad Gaussian should be much broader (and so less informative) than the good Gaussians. Thus we expect that a suitable value for s should be perhaps five to ten times larger than the typical uncertainties in the good Gaussians. Next, it would be unsatisfactory if the results depended too strongly on the particular value of s that we use. Thus, as a check, we can vary s around the value we choose, and see if the results change too drastically. Finally, the value of <β> that emerges from the analysis is important. If it is too small, that implies that the broad Gaussians are making the major contribution to the results, and this would be difficult to justify.
For the CRONUS-P sample, we find that s=0.35 works very well. As already mentioned, we take the prior ƒ0(X,β), to be a constant; thus we do not build in any strong presumptions as to the likely values of X or β. Figure 1 shows the pdf for X that we obtain at s=0.35. The Press method yields the following results for CRONUS-P:
$$\,\lt\!X\!\gt{\equals}4.94\quad \sigma _{X} {\equals}0.07\quad \,\lt\!\beta \!\gt{\equals}0.49$$
Note that the pdf for X has a two-peaked structure, with a small amount of statistical weight assigned to the peak near X=5.2 Thus, this pdf is distinctly un-Gaussian, and is not completely characterized by its average and variance.
We repeated the calculations for a range of values of s (Figure 2). We see that once s is greater than about 0.30, the value of <X> changes very little; the value of σ X changes a bit more, but is between 0.07 and 0.06 over a wide range. The value of β is the most sensitive to s, but is only of secondary interest to us.
This illustrates how our CRONUS-P 3He results depend on the parameter s. We plot, σ X , and <β> as functions of s. Note that <X> and σ X are in units of 109 atoms per gram.
It is interesting to compare the conclusions of the Press method with those obtained by Blard et al. (Reference Blard, Balco, Burnard, Farley, Fenton, Friedrich, Jull, Niedermann, Pik, Schaefer, Scott, Shuster, Stuart, Stute, Tibari, Winckler and Zimmermann2015). In order to deal with the overdispersion, they used the added variance method. For each lab, they replaced
$$\sigma _{i}^{2} $$
by
$$\sigma _{i}^{2} {\plus}\xi ^{2} $$
, where ξ is the same for all labs and is adjusted to match the actual dispersion of the results. They obtain X=5.02±0.06. Our interpretation of the difference between the results of the two methods is that the Press method gives less weight to the two labs that produced the high concentrations (5.20 and 5.25), and so gives a lower estimate for <X>. For example, using the equations presented in the “Derivation of the Basic Equations” section, we can compute the probability that any given measurement is correct. Figure 3 shows these results, again for s=0.35. Thus, the probability that the highest result (x=5.25±0.02) is correct is only 0.018, while the probability that the x=4.91±0.05 result is correct is 0.74.
CRONUS-A 14C Results
In this section we discuss the 14C measurements of the CRONUS-A quartz sample. Five labs measured the 14C concentration of this sample; the results are reported in Jull et al. (Reference Jull, Scott and Bierman2015) and Lifton et al. (Reference Lifton, Goehring, Wilson, Kubley and Caffee2015). To apply the Press analysis, we first combined the data from each lab into one final result for that lab. To combine the data from a given lab, we used the usual error-weighted mean. To evaluate the uncertainty, we took a conservative approach; we computed it using three different approaches, and used the largest of the three. These three methods are given by equations in Bevington and Robinson (Reference Bevington and Robinson2003:4.19, 4.22, 4.13–4.14). We then applied the Press method to combine the final results from all the labs.
In Table 1 we list the results for each lab. We see that the labs measure values ranging from 6.69 to 7.04, in units of 105 atoms per gram (we will use units of 105 atoms per gram for X, s, and σ X for the CRONUS-A results). A look at the uncertainties shows that overdispersion is present.
14C CRONUS-A data.
When we apply the Press method, we find that taking s=0.55 works well. Figure 4 shows the resulting pdf for X that we obtain; it has tails with a non-Gaussian character. The pdf yields
$$\,\lt\!X\!\gt\,{\equals}6.89\quad \sigma _{X} {\equals}0.06\quad \,\lt\!\beta \!\gt\,{\equals}0.65$$
This shows the (unnormalized) pdf for the concentration of 14C in the CRONUS-A sample. The x-axis is in units of 105 atoms per gram. We used a value of s=0.55, as discussed in the section “CRONUS-A 14C Results.”
Figure 5 shows the s dependence of the results. These are very similar to the findings for CRONUS-P: once s is greater than about 0.3, <X> changes very little, σ X is slightly more sensitive, and β is most sensitive.
This illustrates how our CRONUS-A 14C results depend on the parameter s. We plot <X>, σ X , and <β> as functions of s. Note that <X> and σ X are in units of 105 atoms per gram.
The most extreme outlier is the sample with x=6.69±0.06; its probability of being a correct result is p=0.27. An interesting contrast is the sample with x=7.03±0.14. The probability that this result is good is p=0.78. Figure 6 shows these results for all the labs.
As we did for CRONUS-P in the last section, we can compare the results of the Press method for CRONUS-A with the results of the added variance method. We used RadialPlotter (Vermeesch Reference Vermeesch2009) to do the calculation, and obtained <X>=6.86 and σ X =0.046. These are very close to the Press method results, although the pdf’s for X will differ in their detailed form.
Other CRONUS Samples
In the previous two subsections (“CRONUS-P 3He Results” and “CRONUS-A 14C Results”) we presented the detailed workings of the Press method in its application to two CRONUS samples. Here, we report the results of applying the method to other CRONUS samples. These include CRONUS-A (10Be and 26Al) and CRONUS-N (10Be and 26Al); we also analyzed CREU-1 (21Ne). For CREU-1, we used the 21Ne data presented in Figure 2 of Vermeesch et al. (Reference Vermeesch, Balco, Blard, Dunai, Kober, Niedermann, Shuster, Strasky, Stuart, Wieler and Zimmermann2015). For the other samples, we used a two-step procedure, similar to that described in the previous section on CHRONUS-A 14C results. That is, we first combined the data from a given lab to obtain a final, summary result for that lab; to do this, we used the data shown in Tables 2a and 2b of Jull et al. (Reference Jull, Scott and Bierman2015). We then used the Press method to combine the results from these labs. For comparison, we also combined the results from the various labs using the added variance method. Table 2 gives the results of our calculations. The overdispersion ratio is defined as
$$\xi /\!\mu $$
, where μ is the average value. Appendix A of Phillips et al. (Reference Phillips, Argento, Balco, Caffee, Clem, Dunai, Finkel, Goehring, Gosse, Hudson, Jull, Kelly, Kurz, Lal, Lifton, Marrero, Nishiizumi, Reedy, Schaefer, Stone, Swanson and Zreda2016) reports the results of a similar two step approach, in which the added variance method is used to combine the results from the different labs. Those authors used data only from labs that measured a given sample at least three times, and only considered a sample that was measured by at least three labs. For our calculations we used all the data.
Results for all CRONUS intercomparison materials.
We comment on a few significant aspects of our results. For the CRONUS-A 26Al samples, the differences among the various laboratories are quite small, and the calculated overdispersion is in fact zero: ξ=0. This is also reflected in the Press analysis; for these samples we have <β>=0.805, which the highest listed in Table 2. The CREU 21Ne is at the opposite extreme; this sample has the highest overdispersion ratio, equal to 3.5%. It also has <β>=0.30, which is the smallest in the table.
In Figure 7 we plot, for all the samples analyzed, the overdispersion ratio versus the value of <β> . Overall, there is a strong correlation; high values of <β> correspond to a low value of the dispersion ratio. The two 26Al samples, from CRONUS-A and CRONUS-N, have the smallest overdispersion ratios and the largest values of <β>. This reflects the fact that the uncertainties assigned to the results from the various labs were most consistent with the spread in the results from the labs.
This plots, for each sample material, the value of the overdispersion ratio versus the value of <β>. The overdispersion ratio is defined to be
$$\xi /\!\mu $$
, and one σ uncertainties are shown for β.
DISCUSSION AND CONCLUSION
We provided a detailed review of the Press method for combining data, and applied it to the CRONUS cosmogenic nuclide data sets. We stress that this method is broadly representative of a whole class of similar approaches, some of which are designed to cover more complex situations (Christen Reference Christen1994; Dose Reference Dose2003; Hanson Reference Hanson2005; Bronk Ramsey Reference Bronk Ramsey2009; Blaauw et al. Reference Blaauw, Bakker, Christen, Hall and van der Plicht2007). When we are combining results from several labs, and the results suffer from overdispersion, the Press method provides an alternative to the more usual procedure of adding an additional amount to the variance of each result. The Press approach is computationally simple and is based on an appealing logic. As discussed in the Introduction, to decide which procedure is optimal for a given data set, we must distinguish among several different causes for the overdispersion.
One point to emphasize is that there is no universal, model-free way to handle overdispersion. We must formulate a definite model for what has happened. The added variance procedure assumes that each lab has underestimated its dispersion by roughly the same amount, and so adds the same ξ
2 to each
$$\sigma _{i}^{2} $$
. This model could be realistic: for example, each lab could be using the same technique, and the technique could be subject to an extra, unsuspected uncertainty. The Press method formulates a quite different scenario to account for the overdispersion. Press assumes that each lab, with probability β, has made a correct measurement, in the sense that the good Gaussian
$$G_{{\sigma _{i} }} $$
in Equation (4) is suitable. However, Press assumes that for each lab, there is a probability (1−β) that its result is incorrect, and that its result is best included as the broad Gaussian G
s
in Equation (4). As Press stresses, the Bayesian nature of his method is reflected in the way the probability β is handled: it is averaged over, and its pdf is affected by the data.
We placed emphasis on the value of crucial parameter s; that is, it is best if the results are not very sensitive to the choice of s over a range of values. Its important to understand that in choosing a value of s to use in the broad, uninformative Gaussian of Equation (4), we are not trying to replace the quoted uncertainty σ
i
with a more accurate value. Rather, we are replacing the correct Gaussian with a much less informative one. Thus, both the role and the magnitude of s are very different from that of the parameter
$$\xi $$
used in the added variance method.
It is instructive to consider a situation in which the two methods would differ in a significant way (Press Reference Press1997). Suppose we are combining a set of results, and that all the results agree with each other very well except for one outlier which also quotes a very small uncertainty. We also assume that we have no prior reason to distrust the outlier. The Press approach will automatically reduce the influence of such an outlier on the final result. The added variance method will not do this; in fact, if the outlier claims a smaller uncertainty than the other samples, it will be given the most weight in determining the final answer. As we discussed in the Introduction, there is no universal answer as to which of these outcomes is best.
The samples we analyzed exhibited a wide range of overdispersion; the CRONUS-A 26Al had <β>=0.805 (and an overdispersion of zero), while the CREU 21Ne had <β>=0.302 (and an overdispersion ratio of 3.5%). For two of the CRONUS samples, we gave a detailed discussion of how to use the Press method. We found that for the CRONUS-P 3He data set, in which the results of six labs were combined to determine the 3He concentration of a pyroxene sample, the Press method yielded a somewhat smaller concentration than the added variance employed by Blard et al. For the CRONUS-A 14C data set, in which the results of five labs were combined, the Press approach gave an answer very close to the result given by the traditional method. This agreement between two quite different methods is a good indication of the robustness of the CRONUS-A 14C consensus result.
Acknowledgments
This work was supported by National Science Foundation grant EAR 0112480.
