Spectrochimica Acta Part B 63 (2008) 142 – 153 www.elsevier.com/locate/sab
Limits of detection and decision. Part 3 ☆ E. Voigtman Department of Chemistry, University of Massachusetts-Amherst, 710 N. Pleasant St., Amherst, MA 01003-9336, United States Received 28 August 2007; accepted 12 November 2007 Available online 21 November 2007
Abstract It has been shown that the MARLAP (Multi-Agency Radiological Laboratory Analytical Protocols) for estimating the Currie detection limit, which is based on ‘critical values of the non-centrality parameter of the non-central t distribution’, is intrinsically biased, even if no calibration curve or regression is used. This completed the refutation of the method, begun in Part 2. With the field cleared of obstructions, the true theory underlying Currie's limits of decision, detection and quantification, as they apply in a simple linear chemical measurement system (CMS) having heteroscedastic, Gaussian measurement noise and using weighted least squares (WLS) processing, was then derived. Extensive Monte Carlo simulations were performed, on 900 million independent calibration curves, for linear, “hockey stick” and quadratic noise precision models (NPMs). With errorless NPM parameters, all the simulation results were found to be in excellent agreement with the derived theoretical expressions. Even with as much as 30% noise on all of the relevant NPM parameters, the worst absolute errors in rates of false positives and false negatives, was only 0.3%. © 2007 Elsevier B.V. All rights reserved. Keywords: Limit of detection; Detection limit; Heteroscedastic; Bias
1. Introduction
2. Background
Parts 1 and 2 dealt with the detailed statistical exploration of a simple univariate chemical measurement system (CMS) with homoscedastic, Gaussian measurement noise and ordinary least squares (OLS) processing of the calibration curve data. In practice, it is often the case that real chemical measurement systems have measurement noises that are dependent on analyte content. Accordingly, in this paper, the model system is extended to account for heteroscedastic, Gaussian measurement errors and weighted least squares (WLS) processing of calibration curve data. A major goal is elucidation of the experimental decision, detection and quantification limits that arise from Currie's classic hypothesis testing formulation of detection theory [1].
In Part 1, it was shown that Hubaux and Vos's instantiation [2] of Currie's decision and detection limit schema resulted in negatively biased detection limits as a consequence of the variable width prediction intervals about the calibration curve regression line. This result was obtained with the simplest of measurement noises (homoscedastic, Gaussian white noise), the simplest calibration curve data processing (OLS) and the standard measurement protocol: perform replicate measurements on a specimen under test, average the results, subtract the calibration curve intercept and finally divide by the calibration curve slope to obtain the point estimate of the content of the specimen under test. With heteroscedastic, Gaussian white noise, the only modification necessary in Hubaux and Vos's treatment is replacement of OLS with weighted least squares (WLS). Unfortunately, the prediction intervals still vary in width as a function of the content variable, so the resulting detection limits are still biased. As a consequence, the method must be rejected. In Part 2, it was shown that there is never any need for critical values of the non-central t distribution when the noise, calibration
☆
This article is published in a special honor issue dedicated to Jim Winefordner on the occasion of his retirement, in recognition of his outstanding accomplishments in analytical atomic and molecular spectroscopy. E-mail address:
[email protected]. 0584-8547/$ - see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.sab.2007.11.012
E. Voigtman / Spectrochimica Acta Part B 63 (2008) 142–153
curve data processing and measurement protocol are as in the preceding paragraph. Indeed, Tables 3 and 4 and Fig. 8 in Part 2 definitively show that use of Currie detection limit expressions based upon critical values of the non-central t distribution are certain to result in bias, even using the MARLAP sanctioned correction factor c4 [[3], eqn 20.65], for all degrees of freedom from 1 to 10. As a further demonstration of the futility of using this method, consider Fig. 1, which shows what happens when the CMS is simplified to the maximum extent, discarding even the calibration curve itself. In this Monte Carlo simulation, N = 2 independent replicate measurements are made from a true blank and their sample mean, y¯0, and sample standard deviation, s0, are then computed. This procedure is repeated 10 million times, in 10 simulations of 1 million steps each. For each simulation step, the χ dis-
143
tributed variate yC is recomputed via yC = tpη1/2 s0, where tp = 6.313751514 for p = 0.05 and ν = 1, and η1/2 is simply 1=2 g1=2 u M01 þ N 1 ¼ ½ð1=1Þ þ ð1=2Þ1=2 c1:224744871
ð1Þ
with M0 = 1 being the number of replicates averaged in constructing the 4 net test responses shown in Fig. 1. Testing for false positives simply involves repeatedly making single measurements of a true blank, subtracting the associated y¯0 value to generate the net test blank and then comparing with the associated yC decision limit. Perfect performance would be 50 000 ± 0 false positives. The actual values were 50 036.9 ±210.6, where the ± term is the standard deviation for the 10 results. Clearly, the numerical discrepancy is statistically insignificant. Testing for false negatives involves generating additional net test blanks and then shifting them by yD (for NR1),
Fig. 1. Monte Carlo simulation model demonstrating failure of the MARLAP detection limit expression [[3], eqn 20.66] for one degree of freedom.
144
E. Voigtman / Spectrochimica Acta Part B 63 (2008) 142–153
δ0.05,0.05,1η1/2σ0 (for NR2) and δ0.05,0.05,1η1/2s0/c4 (for NR3). The critical delta value, found by numerical integration of the relevant non-central t distribution, is approximately 12.527 for p = 0.05 = q and ν = 1. Since p = q and the noise is homoscedastic and Gaussian, yD = 2yC. Thus yD is χ distributed and the NR1 test values, judged against yC, yield 50 004.9 ± 347.8 false negatives. Similarly good performance is obtained with the NR2 test values, which absolutely depend upon σ0 being known. In this case, YD (=δ0.05,0.05,1η1/2σ0) is approximately 1.5342, since σ0 = 0.1 in the simulation. But this result is much inferior to the equally valid YD expression: YD ¼ zp g1=2 r0 þ zq g1=2 r0 c0:40291
ð2Þ
given that zp = zq = 1.644853627. Finally, the results for NR3, computed as per MARLAP [[3], eqn 20–66], are unacceptably poor: 33 809.5 ± 281.0 false negatives. Note, too, that δ0.05,0.05,1η1/2s0/c4 is χ distributed and therefore qualitatively the same as yD: it simply is not usable because of its intrinsic bias. In 1997, Currie attempted to devise viable detection limit expressions, containing critical values of the non-central t distribution, for heteroscedastic Gaussian noise in the CMS [[4], pp. 158–9]. The reported results were poor, with a rate of false negatives (0.023 ± 0.005) that was more than a factor of two lower than the nominally specified q = 0.05. In his footnote 9 he stated “the conjectural solution overestimated the detection limit in this particular test by about 25%. Work is underway both to explore the problem theoretically, and to create a correction function c(ν,m), where m is the heteroscedasticity parameter”. As has been shown above, however, the fundamental problem was in the use of detection limit expressions containing critical values of the non-central t distribution. The MARLAP documents also attempt the heteroscedastic extension of the critical values of the non-central t distribution scheme, even going so far as to detail an iterative procedure for its implementation [[3], pp. 20-54–20-58]. However, in view of the serious failure of the scheme for homoscedastic noise, regardless of whether calibration curves are used (see Part 2) or not, and in view of the unqualified success of the “tp + tq” method, for all degrees of freedom, no further use is made of the ill-fated “critical delta values” scheme and it is rejected. In what follows, it will be shown that the “tp + tq” method holds the key to the valid instantiation of Currie's detection limit scheme when the measurement noise is additive, zero mean, heteroscedastic, Gaussian white noise.
Assume a univariate chemical measurement system (CMS) characterized by the response domain linear calibration model f ð X ÞuF ð X Þ þ eð bX Þua þ bX þ eð bX Þua þ Y ð X Þ þ eð bX Þ
where X is the errorless independent variable in the content domain, α is the true intercept, β is the true slope, F(X) is the theoretical dependent variable, ɛ(βX) is additive, zero mean, white, Gaussian (normal) noise, f(X) is the experimental dependent variable and Y(X) is the theoretical net response. As before, X will be a concentration, quantity or amount, in appropriate units, and X ≥ 0. It is assumed that there are no systematic errors and, without loss of generality, β N 0. The units of Y and σ0 are the same, but all units are ignored below, since they are irrelevant. Further details concerning heteroscedasticity assumptions for ɛ(βX) are given below. It is assumed that weighted least squares (WLS) is performed on N data pairs of the form Xi, f¯i, where each f¯i is the sample mean of Mi independent replicate measurements per Xi standard. The number of degrees of freedom, ν, is N − 2 throughout and the Xi standards are assumed to have zero error. Once the N measurements necessary to construct a calibration curve have been properly made and WLS has been performed, it is possible to compute the WLS slope, b, the WLS intercept, a, their sample standard errors, sb and sa, respectively, the sample standard error about regression, swr, and so on, with the WLS calibration curve expression being simply f ð xÞ ¼ a þ bx ¼ a þ yð xÞ:
P
In Parts 1 and 2, the nomenclature was entirely conventional, as befitted the assumptions concerning the CMS, the measurement noise, the OLS processing of the calibration curve data, the standard measurement protocol and so on. However, one of the problems that arise when the extension to heteroscedastic noise is made is the easy confusion of the CMS system's response and its net response. Accordingly, a more carefully developed nomenclature will be employed below, to protect against this error.
ð4Þ
Comparing Eqs. (3) and (4), it is evident that f(x), a, b and y (x) are experimental domain point estimates of F(X), α, β and Y (X), respectively. Both a and b are Gaussian distributed: a ∼ N: α,σa and b ∼ N:β,σb, with σa and σb being their population standard errors, respectively. Unlike X, which is either a theoretical value or one of the errorless standards, Xi, used to construct the calibration curve, x is a random variate point estimate of X. Let Xu be the unknown true value of the content of a specimen under test. Then from Eq. (3), the response of the CMS to Xu is f`(Xu). Hence, if ɛ(βX) ∼ N:0,σ(βX), then f(Xu) ∼ N:α + βXu,σ(βXu). If M independent replicate measurements are made at X = Xu, then their sample mean, denoted f¯(Xu), is distributed as N:α + βXu,σ(βXu)/ √M. Thus, since ɛ(βX) is Gaussian distributed, so is f¯(Xu). The next step is seemingly trivial, but vitally important: f ðxu Þu f ðXu Þ
3. Theoretical
ð3Þ
ð5Þ
since now WLS calibration Eq. (4) may be used, yielding P
yðxu Þ ¼ f ðxu Þ a ¼ f ðXu Þ a
ð6Þ
with distribution yðxu ÞfN : bXu ; rd ðXu Þ
ð7Þ
where r2d ðbXu Þur2a þ r2 ðbXu Þ=M :
ð8Þ
E. Voigtman / Spectrochimica Acta Part B 63 (2008) 142–153
The Gaussian distributed random variate y(xu) is the experimental domain net response variate. Division by b then yields the experimental domain content variate xu: xu ¼ yðxu Þ=b
ð9Þ
whereas the corresponding theoretical domain content value is Xu ¼ Y ðXu Þ=b:
ð10Þ
Since Yu = βXu, Eq. (8) can equivalently be written as r2d ðYu Þ
¼
r2a
þ r ðYu Þ=M
ð11Þ
2
where σ (Yu) is the variance of ɛ(Yu). As a consequence of Eq. (11), the variance of the net response never equals that of the mean of M responses unless M = 1 and α is substituted for a in Eq. (6). This will rarely, if ever, be true since α is almost never known. 2
3.1. Homoscedastic noise with WLS or OLS processing Before proceeding to the detailed consideration of various noise precision models (NPMs) for heteroscedastic, Gaussian noise, it is useful to consider the “bridging” case where the noise is homoscedastic and Gaussian, but the number of replicate measurements per standard, denoted by Mi, is non-constant. In this case, WLS is more efficient than OLS, providing better estimates of σa and σb. This situation is optimal for WLS because the raw weights, denoted by wiraw, are simply proportional to the Mi values and therefore are errorless. The weight normalization factor, Ns, is defined as Ns uN =
N X
wraw i uN =
i¼1
N X
N X P Mi =r20 ¼ N r20 = Mi ur20 = M i
i¼1
ð12Þ
i¼1
where σ 0 is the population standard deviation of the homoscedastic noise and M¯i is the average number of replicates per standard. Thus the normalized weights are simply wi uNs wraw i
ð13Þ
and the normalization ensures that the weights conveniently sum to N, the number of standards. Then the weighted mean and weighted mean square of the standards are given by P
X w uN 1
N X
wi Xi
ð14Þ
i¼1
and Swxx uN 1
N X
P wi Xi X w
ð15Þ
i¼1
and the resulting “η1/2” factor is g1=2 w ¼
h P P 2 i1=2 M i =M0 þ ð1=N Þ þ X w =Swxx
ð16Þ
145
where M0 is the number of independent replicates of the blank used to define the Currie decision limit. The “w” subscript on η signifies that it is computed from weighted quantities. However, since the noise is homoscedastic and all the weights are true weights, Eq. (16) is errorless and ηw1/2 simply a constant determined by the specific values of the calibration standards and the measurement protocol. This will not be true when heteroscedastic noise is considered. As in Parts 1 and 2, it immediately follows that the quadrant 1 results are YC ¼ zp rd ¼ zp rwr g1=2 w
ð17Þ
and YD ¼ zp rd þ zq rd ¼ zp þ zq rd ¼ zp þ zq rwr g1=2 w
ð18Þ
where σwr is the population standard error about regression. It further happens that qffiffiffiffiffiffi P ð19Þ rwr ¼ Ns1=2 ¼ r0 = M i so YC and YD could immediately be computed if σ0 were known. In quadrant 2, the corresponding expressions are yC ¼ tp sd ¼ tp swr g1=2 w
ð20Þ
and yD ¼ tp sd þ tq sd ¼ tp þ tq sd ¼ tp þ tq swr g1=2 w
ð21Þ
where swr is the sample standard error about regression produced by the WLS processing and is related to the sample standard deviation of the blank by qffiffiffiffiffiffi P ð22Þ swr ¼ s0 = M i Since swr is χ distributed and tp, tq and ηw1/2 are constants, both yC and yD are χ distributed. Hence, both quadrant 4 variates, xC and xD, are distributed as modified non-central t variates. Note that it is always true that XC = YC/β and XD = YD/β, in quadrant 3, and likewise xC = yC/b and xD = yD/b, in quadrant 4. If the number of replicates per standard is a constant, denoted MX, then WLS with normalized weights reduces to OLS since all the normalized weights are unity. Then Eqs. (12)–(22) simplify by replacing all factors of wi by unity, dropping all “w” subscripts, and replacing M¯i by MX. Further, if MX is unity, then the situation is as in Parts 1 and 2 and sr = s0. 3.2. Heteroscedastic noise with WLS processing Now suppose the noise is heteroscedastic and Gaussian distributed. By far the most commonly encountered situation in experimental work is when the noise increases with X, having its minimum value for the blank, i.e., at X = 0, so this will be assumed below. There must exist a true underlying noise precision relationship for ɛ(βX) and this is based on the fundamental physico-chemical principles upon which the CMS is based. The best case scenario is when the true noise precision
146
E. Voigtman / Spectrochimica Acta Part B 63 (2008) 142–153
relationship (model) is known, as are all of its controlling parameters, while the worst case scenario is that little is known of the noise precision model (NPM) or its parameters. In the discussion that follows, it is assumed initially that the NPM is known, through whatever means, and that the true parameters of the NPM are also known. In reality, this is virtually certain to never happen, so the assumption that the parameters are known will be relaxed subsequently, but it will still be assumed that the true NPM is known. If this assumption proves false, then the situation is uncomfortably close to an exercise in curve fitting. The three simplest physically plausible NPMs are the linear, quadratic and “hockey stick” NPMs, the latter so-called because of its evocative ɛ(βX) vs X curve. This NPM is referred to by Gibbons and Coleman as the “approximate Rocke and Lorenzato” NPM [[5], p. 68, eqn 7.22]. For each of these NPMs, the Currie decision and detection limits, in both quadrants 1 and 2, are separately derived below. Repeating a vitally important statement made above, the corresponding quadrant 3 and 4 results are immediately obtained from XC = YC/ β and XD = YD/β, in quadrant 3, and xC = yC/b and xD = yD/b, in quadrant 4. 3.2.1. Linear noise precision model In this model, the NPM is simply eð bX ÞfN : 0; r0 þ AbX ¼ N : 0; r0 þ AY
ð23Þ
where μ is the heteroscedasticity population parameter and σ0 is the population standard deviation of the blank. The raw 2 weights are thus wraw i u1=ðr0 þ AbXi Þ . At small values of X, the noise has almost constant population standard deviation, while μ is the asymptotic population RSD: A ¼ lim eð bX Þ=bX :
ð24Þ
X Yl
with the positive root being the correct one. To prevent B = 0, thereby averting “infinite” YD, it is necessary that 1=2
AbM0 =zp
so, as Currie noted [4], the detection limit fails to exist for sufficiently high values of μ. Note that for M0 = 1 and 95% confidence, zp = 1.644853627, so μ b 0.607957 is required. Eq. (27) may be expanded to yield r2d ðYD Þ ¼ r2d þ M01 2r0 AYD þ A2 YD2
ð25Þ
ð31Þ
Clearly, even though the assumed NPM (Eq. (23)) is remarkably simple, Eq. (31) is necessarily more complicated due to the addition of variances when the net response is computed. Unfortunately, it is a common mistake [[4], p. 160, eq. 16] to assume that σd(YD), once Eq. (23) is assumed, can then be re-modeled by an equation of the form rd ðYD Þ ¼ rd þ AYD
ð32Þ
which, as Eq. (31) proves, is completely wrong. Indeed, once the NPM is specified, σd(YD) is implicitly defined and cannot then be arbitrarily re-modeled in any way. In quadrant 2, σ0 and μ are estimated by s0 and m, respectively, and Eqs. (25)–(31) become yD ¼ yC þ t q s d ð yD Þ ¼ t p s d þ t q s d ð yD Þ
ð33Þ
s2d ¼ s2a þ s20 =M0 ¼ s2d ðyu u0Þ
ð34Þ
s2d ðyD Þ ¼ s2a þ ðs0 þ myD Þ2 =M0 ¼ s2d ðyu uyD Þ
ð35Þ
In quadrant 1, YC will still be given by Eq. (17), but Eq. (18) is replaced by the following general equation: YD ¼ YC þ zq rd ðYD Þ ¼ zp rd þ zq rd ðYD Þ
ð30Þ
1
(
yD ¼ B ̂
) 2 1=2 2 2 2 yC þ A ̂ þ yC þ A ̂ yC B ̂ 1 tq =tp
ð36Þ
where r2d ur2a þ r20 =M0 ¼ r2d ðYu u0Þ
ð26Þ
A ̂ us0 tq2 m=M0 & B ̂ u1 m2 tq2 =M0
ð37Þ
and 1=2
r2d ðYD Þ
¼
r2a
2
þ ðr0 þ AYD Þ =M0 ¼
r2d ðYu uYD Þ
ð27Þ
mbM0 =tp
ð38Þ
and Eqs. (26) and (27) follow directly from Eq. (11), with M = M0 and the indicated values of Yu. Combining Eqs. (25)–(27) yields a quadratic equation in YD, having the following solution: h i1=2 YD ¼ B1 Yc þ A þ ðYc þ AÞ2 YC2 B 1 z2q =z2p ð28Þ where Aur0 z2q A=M0 & Bu1 A2 z2q =M0
ð29Þ
1=2 : sd ðyD Þ ¼ s2d þ M01 2s0 myD þ m2 y2D
ð39Þ
Note that Eqs. (33) and (34) are always true, regardless of NPM. A sample of linear NPM noise is shown in Fig. 2 (upper left), where σ0 = 0.1, μ = 0.3 and β = 2. For each of the six standard values (X = 1, 2,…, 6) and the blank (at X = 0), 400 samples of the noise are plotted. Clearly, the noise is strongly heteroscedastic and use of OLS would be ill-advised.
E. Voigtman / Spectrochimica Acta Part B 63 (2008) 142–153
147
Fig. 2. Four traces showing 2800 independent replicates each (400 for the blanks and for each of the 6 standards) for the three noise precision models (NPMs) under study.
3.2.2. “Hockey stick” noise precision model In this model, the NPM is simply
eð bX ÞfN : 0; r20 þ A2 b2 X
2 1=2
¼ N : 0; r20 þ A2 Y
and 2 1=2
ð40Þ
1 h i1=2 yD ¼ yC 1 tq2 m2 =M0 1 þ 1 1 tq2 m2 =M0 1 tq2 =tp2 :
ð44Þ
Note that the upper limits on μ and m are the same as for the linear NPM: Eqs. (30) and (38), respectively. A sample of “hockey stick” NPM noise is shown in Fig. 2 (upper right), where σ0 = 0.1, μ = 0.3 and β = 2.
where, again, μ is the heteroscedasticity population parameter and σ0 is the population standard of 2deviation the blank. The 2 2 2 raw weights are thus wraw u1= r þ A b X i 0 i . Eqs. (25) and (26) must always be true, but now the NPM in Eq. (40) causes Eq. (27) to become slightly simpler: r2d ðYD Þ ¼ r2a þ r20 þ A2 YD2 =M0 ¼ r2d þ A2 YD2 =M0 ð41Þ
3.2.3. Quadratic noise precision model In this model, the NPM is
Combining Eqs. (25), (26) and (41) yields the following solution for YD:
eð bX ÞeN : 0; r0 þ AbX þ qð bX Þ2 ¼ N : 0; r0 þ AY þ qY 2
1 h i1=2 : YD ¼ YC 1 z2q μ2 =M0 1 þ 1 1 z2q μ2 =M0 1 z2q =z2p
where μ and ρ are heteroscedasticity population parameters and σ0 is the population standard deviation of the blank and the raw weights should be obvious. Again Eqs. (25) and (26) are true, but now
ð42Þ
As for the linear NPM, in quadrant 2, σ0 and μ are estimated by s0 and m, respectively, and Eqs. (41) and (42) become ð43Þ s2d ðyD Þ ¼ s2a þ s20 þ m2 y2D =M0 ¼ s2d þ m2 y2D =M0
2 r2d ðYD Þ ¼ r2a þ r0 þ AYD þ qYD2 =M0 ¼ r2d ðYu uYD Þ:
ð45Þ
ð46Þ
148
E. Voigtman / Spectrochimica Acta Part B 63 (2008) 142–153
Fig. 3. Monte Carlo simulation model showing how false positives and false negatives are computed for the NPMs under study.
Expansion of Eq. (46) yields the following quartic: 1
r2d ðYD Þ ¼ r2d þ M0
2
2r0 AYD þ 2r0 q þ A YD2 þ 2AqYD3 þ q2 YD4 :
ð47Þ
There are a number of ways to find YD. One method is to assume YD = 2YC, solve Eq. (47) for σd(YD) and then use Eq. (25) to iteratively refine the estimate of YD. This was conveniently implemented in a Microsoft Excel spreadsheet, using the built-in Solver function. In quadrant 2, σ0, μ and ρ are estimated by s0, m and r, respectively, and Eq. (47) becomes s2d ðyD Þ ¼ s2d þ M01 2s0 myD þ 2s0 r þ m2 y2D þ 2mry3D þ r2 y4D :
ð48Þ
In the Monte Carlo computations described below, Newton– Raphson root finding was used, since the polynomial nature of Eq. (48) facilitated efficient generation of the necessary derivative and convergence is quadratic. Note that YD and yD will fail to exist under sufficiently heteroscedastic conditions, but upper limits are not readily obtained given that the root finding also is not absolutely guaranteed to converge. A sample of quadratic NPM noise is shown in Fig. 2 (lower left), where σ0 = 0.1, μ = 0.2, ρ = 0.002 and β = 2. 3.2.4. Limits of quantification Although not strictly necessary to the theory of decision and detection limits, limits of quantification were an integral part of
Currie's original formulation and have remained so. Since the theory follows immediately from the above sections, it is simple to exhibit the expressions explicitly. Let Q be the desired RSD at the limit of quantification. This can be arbitrary and Currie simply used Q = 0.1 as a default [4]. The population RSD is then defined as follows: RSDpop at Y uYQ urd YQ =YQ
ð49Þ
where YQ is the quadrant 1 limit of quantification. The sample RSD is similarly defined as RSDsample at yuyQ usd yQ =yQ
ð50Þ
where yQ is the quadrant 2 limit of quantification. These expressions are set equal to Q and then solved for the limits of quantification. For homoscedastic, Gaussian noise, σd (YQ) = σd and sd (yQ) = sd, so YQ = σd/Q and yQ = sd/Q, as is well known. With Q = 0.1, these are simply YQ = 10σd and yQ = 10sd. For the linear NPM, Eqs. (49) and (50) are set equal to Q, each is squared and then Eqs. (31) and (39) are substituted, resulting in a pair of quadratic equations that have the following solutions: h 1 i1=2 YQ ¼ Q2 A2 =M0 ðr0 A=M0 Þ þ ðr0 A=M0 Þ2 þr2d Q2 A2 =M0
ð51Þ
E. Voigtman / Spectrochimica Acta Part B 63 (2008) 142–153
149
Fig. 4. Dialog box for the “Linear calibration curves 1” block in Fig. 3.
and h 1 i1=2 yQ ¼ Q2 m2 =M0 ðs0 m=M0 Þ þ ðs0 m=M0 Þ2 þs2d Q2 m2 =M0
ð52Þ
where the positive roots were the correct ones. Clearly, for a given value of Q, the limits of quantification will fail to exist unless μ b Q√M0 and m b Q√M0, respectively. For the “hockey stick” NPM, Eqs. (41) and (43) are used in place of Eqs. (31) and (39) and the results are 1=2 YQ ¼ rd Q2 A2 =M0 ð53Þ and 1=2 yQ ¼ sd Q2 m2 =M0
ð54Þ
with the same limits on μ and m. For the quadratic NPM, a pair of quartic equations results: 0 ¼ YQ4 þ 2Aq1 YQ3 þ 2r0 q1 þ A2 Q2 M0 q2 YQ2 ð55Þ þ 2r0 Aq2 YQ þ r2d M0 q2 and 0 ¼ y4Q þ 2mr1 y3Q þ 2s0 r1 þ m2 Q2 M0 r2 y2Q þ 2s0 mr2 yQ þ s2d M0 r2 :
ð56Þ
As ρ → 0, Eq. (55) reduces to Eq. (51) and as r → 0, Eq. (56) reduces to Eq. (52). As for the equations for YD and yD in the quadratic NPM, these equations can be solved by a variety of means. For “one-off” evaluation, Microsoft Excel's Solver function works well, but real solutions to Eqs. (55) and (56) will fail to exist if the noise is sufficiently heteroscedastic. 3.2.5. Specifying the noise precision model Given that a CMS must have a true underlying NPM, failure to use the correct NPM automatically calls into question the validity of any decision, detection and quantification limits computed for the CMS. Indeed, if an incorrect NPM is used, it could be argued that such limits are mere formulaic prescriptions for certain error. Hence, whenever the noise is heteroscedastic, it is vitally important to know as much about the noise structure as possible and be as certain as possible that the NPM has been correctly identified. Typically, this involves the usual mix of theoretical considerations, prior literature concerning the CMS, experience with the CMS or related systems, and so on. It cannot be emphasized too strongly that if a mistake occurs in the determination of the NPM, the consequences can be seriously detrimental with regard to all subsequent results and conclusions. Assuming that the true NPM has been identified, the next step is estimation of the population parameters that control it. In the
150
E. Voigtman / Spectrochimica Acta Part B 63 (2008) 142–153
preceding sections, the population parameters for the linear and “hockey stick” NPMs were σ0, the population standard deviation of the analytical blank, and μ, the asymptotic population RSD. For the quadratic NPM, an additional heteroscedasticity parameter, ρ, is also needed. These parameters are estimated by s0, the sample standard deviation of the analytical blank, m, the asymptotic sample RSD, and r, respectively. In practice, values of s0, m and r (if the quadratic NPM is used) are found by making many replicate measurements at a number of appropriately chosen standard X values and then unweighted regression is performed on the Xi,si data pairs, where si is the sample standard deviation of the replicates at Xi. Since s0, m and r are random variates, so, too, are all quantities computed from them, including wiraw, Ns, wi, X¯w, Swxx, ηw1/2, swr, and sd. Note that the random variate ηw1/2 is now given by g1=2 w ¼
h
P i1=2 2 s20 =M0 Ns þ ð1=N Þ þ X w =Swxx
ð57Þ
rather than by Eq. (16). Any given set of Xi,si data pairs will result in a fixed value of ηw1/2, but it is, nevertheless, a sample from a distribution of such values. Consequently, swr is no longer a χ variate, though it may be approximately so, and, since sd = ηw1/2swr, the same applies to sd. Thus the Currie decision limit yC, given by tpsd, is not χ distributed and xC is not distributed as a modified non-central t variate. For the Currie detection limits yD and xD, the PDFs are even more complicated than those of yC and xC, respectively, and prospects are dim for finding analytic PDFs for yC, xC, yD or xD when the
noise is heteroscedastic. Fortunately, Monte Carlo simulations are so fast and accurate now that there is no compelling need for the analytic PDFs, though they would never be spurned. 4. Experimental To test the above theoretical developments, extensive Monte Carlo calculations were performed. All of the software and hardware details are unchanged from Part 1, quo vide. Given the number of variables involved, all of the simulation results presented below were based on simulations having the following fixed values: N = 6 standards (X = 1, 2,…, 6), each measured only once, σ0 = 0.1, α = 1, and β = 2. The present studies used p = 0.05 = q, p = 0.05 = 2q and 2p = 0.05 = q. Various values of μ and ρ were used, as reported below. It is important to note that one complete calibration curve, and all associated test statistics and computations, were performed for each simulation step, for one million steps per simulation and 10 simulations in a simulation set. This requires several minutes of computation time, but provides the advantage of diagnostic quality histograms of yC, yD and other variates of interest. Fig. 3 shows a typical Monte Carlo simulation model that was used to test the theory presented above. As for the simulations in Parts 1 and 2, this is both the computer model itself and also its own self-documenting block diagram. For each simulation step, M0 = 3 independent true blank replicates were averaged and then blank-corrected by subtracting the WLS intercept from their sample mean. To test for false positives, the resulting net response variates, labeled “NR1”, were then compared with the yC variates computed using Eq. (20).
Fig. 5. Experimental, net response domain 10 million event histograms of yC and yD for an extreme case (μ = 0.45) of the linear NPM.
E. Voigtman / Spectrochimica Acta Part B 63 (2008) 142–153
Testing for false negatives is only slightly more complicated. First, M0 = 3 independent true blank replicates were averaged and then blank-corrected by subtracting the WLS intercept from their sample mean. Next, the resulting raw net response variates were then variance adjusted so that they were samples of a wider Gaussian distribution having sample standard deviation sd(yD) centered at the random variate yD. In the model, this was accomplished by multiplying the raw net response variates by sd(yD)/sd and then adding yD. The resulting net response variates, labeled “NR2”, were then compared with the yC variates to test for false negatives. Regardless of NPM, sd is computed via Eq. (34). For the linear NPM, yD was computed via Eq. (36) and sd(yD) via Eq. (39). For the “hockey stick” NPM, yD was computed via Eq. (44) and sd(yD) via Eq. (43). Lastly, for the quadratic NPM, five iterations of Newton– Raphson root finding were employed, using Eq. (33) for yD, Eq. (48) for sd(yD) and an initial guess of yD = 2yC. As shown in Fig. 3, the heart of the model is the “Linear calibration curves 1” block at left center. The dialog box of this block is shown in Fig. 4. This block is where the desired NPM is specified along with the requisite NPM parameters and other essential parameters and critical values. Note that under each of the σ0, μ, and ρ parameter entry boxes are entry boxes allowing specification of up to 30% normally distributed variation in the nominal parameter values. More specifically, if zeroes are entered in Table 1 Monte Carlo error rates for p = 0.05 and q = 0.05 μ
ρ or % False positives NPM RSD
% Errors False negatives
% Errors
0.3 0.1 0.03 0 0.3 0.1 0.03 0.01 0.003 0.3 0.1 0.03 0.3 0.1 0.03 0 0.3 0.1 0.03 0.01 0.003 0.3 0.1 0.03 0.2 0.2 0.2 0.05 0.05 0.05
LIN LIN LIN LIN LIN LIN LIN LIN LIN LIN LIN LIN HS HS HS HS HS HS HS HS HS HS HS HS 0.002 0.002 0.002 0.003 0.003 0.003
− 0.16 − 0.04 − 0.10 0.04 − 0.01 − 0.06 0.28 0.55 0.00 − 1.37 − 0.93 2.21 − 0.05 − 0.06 0.05 0.03 − 0.08 − 0.69 0.79 0.61 0.02 − 1.27 − 3.71 5.91 − 0.17 0.05 − 0.66 − 0.10 0.03 1.57
−0.13 0.05 0.13 0.02 −0.06 −0.25 0.37 0.27 0.22 −1.49 −0.97 2.23 −0.08 −0.28 0.27 −0.09 −0.17 −0.38 0.99 0.85 0.15 −1.15 −3.83 5.70 −0.07 −0.02 −0.80 −0.17 0.32 1.74
0 0 0 0 10 10 10 10 10 30 30 30 0 0 0 0 10 10 10 10 10 30 30 30 0 10 30 0 10 30
49 918.3 ± 173.2 49 981.3 ± 141.1 49 952.2 ± 164.2 50 019.1 ± 189.4 49 995.6 ± 259.2 49 970.1 ± 231.3 50 138.2 ± 172.9 50 277.0 ± 255.7 50 002.0 ± 150.8 49 313.6 ± 258.1 49 536.9 ± 190.1 51 102.9 ± 155.6 49 976.7 ± 184.4 49 968.1 ± 230.9 50 025.7 ± 223.4 50 015.0 ± 266.9 49 958.4 ± 149.1 49 655.1 ± 198.7 50 393.7 ± 132.0 50 306.9 ± 178.5 50 010.2 ± 175.8 49 364.9 ± 336.5 48 144.5 ± 211.6 52 956.6 ± 164.7 49 916.3 ± 177.2 50 025.6 ± 244.1 49 667.9 ± 191.5 49 951.3 ± 210.1 50 012.9 ± 177.1 50 785.2 ± 142.0
49 934.7 ± 220.1 50 024.9 ± 170.6 50 065.1 ± 266.7 50 008.7 ± 215.1 49 968.9 ± 141.3 49 877.5 ± 216.3 50 184.6 ± 203.0 50 135.4 ± 250.3 50 108.7 ± 236.5 49 254.8 ± 238.2 49 517.1 ± 303.7 51 115.5 ± 206.8 49 958.1 ± 183.1 49 859.5 ± 223.1 50 136.3 ± 231.4 49 953.8 ± 215.7 49 913.2 ± 364.5 49 808.6 ± 133.1 50 494.0 ± 156.1 50 423.6 ± 200.8 50 075.1 ± 243.7 49 425.7 ± 323.4 48 084.5 ± 146.0 52 848.9 ± 217.7 49 962.7 ± 214.3 49 987.9 ± 185.1 49 601.6 ± 117.6 49 913.8 ± 240.3 50 161.3 ± 248.8 50 871.1 ± 207.2
151
these boxes, then σ0, μ, and ρ are used as is. However, if a nonzero value (no greater than 30) is entered in one of these boxes, the corresponding parameter simply becomes the location parameter of a Gaussian distribution for the associated sample test statistic. Thus s0 fN : r0 ; 0:01 kRSD r0
ð58Þ
mfN : A; 0:01 kRSD A
ð59Þ
and rfN : q; 0:01 kRSD q:
ð60Þ
The variability was specified as being Gaussian distributed because s0, m and r (if used) are initially obtained via regression on the noise precision Xi,si data pairs. In Fig. 4, the NPM is shown as being quadratic and the three NPM parameters have 10%RSDs. This models the situation that is probably the best that can be achieved in the real world: the true NPM is known, but its parameters can only be estimated to perhaps one significant figure, if that. A real CMS is highly unlikely to exhibit heteroscedasticities nearly as large as those shown in Fig. 2, especially that shown in the lower right. This case is for a linear NPM with μ = 0.45, which is not far below the theoretical limit, for M0 = 1 and 95% confidence, of 0.607957. 5. Results and discussion Before presenting the main group of results, it is useful to examine one model in some detail and the chosen one is the extreme linear NPM model with σ0 =0.1, μ= 0.45, M0 =1 and p =0.05 =q. At the bottom of the dialog box is a checkbox that allows the “Linear calibration curves 1” block to compute and repeatedly output the theoretical errorless values of YC, YD, XC, XD, σd, and σd (YD) instead of the random variates yC, yD, xC, xD, sd, and sd (yD), respectively. In quadrant 1, the theoretical results were YC = 2.3000, YD = 10.713, XC = 1.1500, XD = 5.3563, σd = 1.3983, σd (YD) =5.1145 and YQ is undefined for any Q ≤ 0.45. Note that YD is only 4.66 times higher than YC, but the noise is so great that XD is only a little below the highest standard used, i.e., X =6. In this case, the NR1 and NR2 variates are Gaussian distributed: NR1fN : 0; rd and NR2fN : YD ; rd ðYD Þ
ð61Þ
and the NR1 and NR2 histogram results were both Gaussian distributed and in excellent agreement with theory: NR1hist ∼ N:0.0002743, 1.3980 and NR2hist ∼N:10.714, 5.1134. In quadrant 2, the 10 million event histograms for yC and yD are shown in Fig. 5, with yC plotted on the bottom and left scales and yD plotted on the top and right scales. Since there was no noise on σ0 or μ, the weights are true and both yC and yD are χ distributed variates. The average values of yC and yD are y¯C = 2.802 and y¯D = 75.46, so y¯D is almost 27 times higher than y¯C and x¯D is far above the highest standard used. This model
152
E. Voigtman / Spectrochimica Acta Part B 63 (2008) 142–153
was deliberately made extreme in almost every regard, yet its performance in terms of false positives (49 980.7 ± 222.5) and false negatives (50 054.2 ± 199.5), was excellent. Table 1 gives the numbers of false positives and false negatives for 12 linear (“LIN”), 12 “hockey stick” (“HS”) and 6 quadratic NPMs, for p = 0.05 = q, σ0 = 0.1, and either 0, 10 or 30%RSDs on each of the operative NPM parameters. Table 2 differs from Table 1 only in that it uses p = 0.05 and q = 0.025 while Table 3 uses p = 0.025 and q = 0.05. In each table, the ± terms are the sample standard deviations for 10 simulations of 1 million calibration curves each, so the three tables summarize the error rate performance of 900 million independent calibration curves. Note that the top NPMs in each NPM category in Table 1 have the typical noises shown in Fig. 2. Consider first the 30 cases where the %RSDs on the relevant NPM parameters are zero. Then, regardless of NPM and regardless of p and q values, the results were in excellent agreement with the theory presented above, with no relative error greater than 0.29% and no statistically significant deviation from the expected error rates. For the linear and “hockey stick” NPMs with 10%RSD on σ0 and μ, regardless of p and q, it was generally found that the error rates were slightly negatively biased for μ = 0.1 and 0.3, switched to positive bias for μ ≤ 0.03 and the positive bias declined as μ decreased. With 30%RSD on σ0 and μ, again regardless of p and q, the error rates were more negatively biased Table 2 Monte Carlo error rates for p = 0.05 and q = 0.025 μ
ρ or % False positives NPM RSD
% Errors False negatives
% Errors
0.3 0.1 0.03 0 0.3 0.1 0.03 0.01 0.003 0.3 0.1 0.03 0.3 0.1 0.03 0 0.3 0.1 0.03 0.01 0.003 0.3 0.1 0.03 0.2 0.2 0.2 0.05 0.05 0.05
LIN LIN LIN LIN LIN LIN LIN LIN LIN LIN LIN LIN HS HS HS HS HS HS HS HS HS HS HS HS 0.002 0.002 0.002 0.003 0.003 0.003
0.22 0.03 0.12 0.29 − 0.18 − 0.23 0.38 0.23 0.24 − 1.49 − 1.12 2.29 0.00 − 0.16 0.11 0.09 − 0.45 − 0.44 0.77 0.56 0.12 − 1.44 − 3.91 6.05 0.10 − 0.06 − 0.83 − 0.04 0.26 1.36
− 0.04 − 0.15 − 0.05 0.06 − 0.05 0.02 0.32 0.59 0.35 − 1.31 − 0.51 3.56 0.17 0.01 0.17 − 0.18 − 0.09 − 0.73 1.30 1.23 − 0.06 − 1.53 − 4.18 9.46 − 0.02 0.13 − 0.58 0.01 0.18 2.72
0 0 0 0 10 10 10 10 10 30 30 30 0 0 0 0 10 10 10 10 10 30 30 30 0 10 30 0 10 30
50 112.0 ± 210.8 50 013.7 ± 191.6 50 058.6 ± 149.3 50 146.1 ± 224.5 49 909.6 ± 210.5 49 885.3 ± 239.9 50 187.8 ± 135.5 50 113.9 ± 165.8 50 117.5 ± 236.6 49 253.8 ± 183.2 49 440.2 ± 286.4 51 146.4 ± 200.5 50 001.0 ± 192.9 49 919.4 ± 274.7 50 054.8 ± 160.6 50 045.0 ± 152.3 49 773.6 ± 192.9 49 780.0 ± 242.5 50 383.3 ± 192.4 50 279.6 ± 242.9 50 058.8 ± 281.3 49 281.8 ± 234.5 48 043.8 ± 148.2 53 027.1 ± 211.1 50 052.1 ± 284.4 49 969.9 ± 312.2 49 585.9 ± 132.2 49 981.6 ± 215.2 50 132.3 ± 256.8 50 678.8 ± 171.8
24 990.1 ± 172.9 24 963.5 ± 148.9 24 988.6 ± 153.9 25 016.2 ± 136.3 24 988.7 ± 158.3 25 005.7 ± 156.3 25 079.7 ± 188.8 25 148.7 ± 126.9 25 088.3 ± 105.9 24 672.4 ± 135.6 24 871.7 ± 134.6 25 889.5 ± 210.3 25 042.4 ± 147.9 25 002.5 ± 104.4 25 043.0 ± 148.1 24 954.9 ± 155.4 24 976.9 ± 157.1 24 817.5 ± 179.9 25 324.9 ± 184.2 25 307.4 ± 198.8 24 986.1 ± 150.1 24 616.3 ± 162.4 23 955.8 ± 78.4 27 364.9 ± 140.3 24 996.2 ± 115.1 25 033.6 ± 133.0 24 854.8 ± 170.2 25 001.7 ± 125.9 25 046.0 ± 107.1 25 681.0 ± 205.8
Table 3 Monte Carlo error rates for p = 0.025 and q = 0.05 μ
ρ or % False positives NPM RSD
% Errors False negatives
% Errors
0.3 0.1 0.03 0 0.3 0.1 0.03 0.01 0.003 0.3 0.1 0.03 0.3 0.1 0.03 0 0.3 0.1 0.03 0.01 0.003 0.3 0.1 0.03 0.2 0.2 0.2 0.05 0.05 0.05
LIN LIN LIN LIN LIN LIN LIN LIN LIN LIN LIN LIN HS HS HS HS HS HS HS HS HS HS HS HS 0.002 0.002 0.002 0.003 0.003 0.003
− 0.12 0.15 − 0.01 − 0.15 − 0.25 0.10 0.08 0.66 0.24 − 1.57 − 0.82 3.25 0.01 − 0.14 − 0.20 0.01 − 0.45 − 0.44 0.94 0.98 − 0.10 − 1.68 − 4.01 8.94 0.16 − 0.18 − 0.86 − 0.19 0.25 2.06
0.03 − 0.04 − 0.07 − 0.01 − 0.04 − 0.05 0.29 0.23 0.23 − 1.20 − 1.17 2.30 0.12 − 0.08 − 0.07 0.03 − 0.02 − 0.72 0.67 0.74 0.20 − 1.20 − 3.89 6.01 0.03 0.17 − 1.10 0.26 0.45 1.54
0 0 0 0 10 10 10 10 10 30 30 30 0 0 0 0 10 10 10 10 10 30 30 30 0 10 30 0 10 30
24 970.8 ± 175.0 25 036.5 ± 125.0 24 998.1 ± 192.3 24 963.6 ± 182.6 24 937.7 ± 127.4 25 024.4 ± 199.9 25 019.9 ± 176.8 25 164.4 ± 153.8 25 058.8 ± 145.7 24 606.4 ± 138.6 24 793.9 ± 186.5 25 813.1 ± 197.3 25 003.6 ± 184.3 24 964.7 ± 220.6 24 950.9 ± 130.2 25 001.5 ± 197.9 24 887.1 ± 159.0 24 889.3 ± 143.5 25 234.8 ± 152.5 25 244.1 ± 105.5 24 974.1 ± 133.5 24 580.9 ± 124.4 23 997.3 ± 170.6 27 234.0 ± 96.9 25 040.1 ± 151.5 24 956.2 ± 209.9 24 786.1 ± 175.6 24 951.4 ± 298.8 25 061.6 ± 142.0 25 515.3 ± 156.7
50 016.2 ± 283.3 49 981.5 ± 190.5 49 963.8 ± 207.0 49 996.4 ± 226.7 49 980.8 ± 216.1 49 973.3 ± 203.2 50 145.5 ± 302.5 50 113.3 ± 197.4 50 115.3 ± 214.9 49 399.5 ± 226.7 49 413.3 ± 252.7 51 150.7 ± 309.5 50 058.4 ± 285.6 49 958.4 ± 138.0 49 964.9 ± 217.1 50 012.8 ± 250.4 49 991.0 ± 138.3 49 639.8 ± 307.5 50 336.6 ± 253.3 50 371.8 ± 182.9 50 101.9 ± 165.4 49 400.8 ± 157.9 48 056.1 ± 221.0 53 007.1 ± 124.5 50 013.2 ± 205.5 50 082.6 ± 148.8 49 449.8 ± 148.0 50 129.5 ± 205.1 50 227.0 ± 313.9 50 769.6 ± 268.5
for μ = 0.1 and 0.3 and switched to more positively biased for μ = 0.03. The cause of this bias switching is unknown at present, but it would not be a surprising minor consequence of the fundamental non-linearities in the detection limit expressions. For the linear NPM with 10%RSD, the maximum relative errors were both at μ = 0.01, with 0.66% relative error in false positives (p = 0.025, q = 0.05) and 0.59% relative error in false negatives (p = 0.05, q = 0.025). To put this in perspective, the absolute errors were only 0.016% and 0.015%, respectively, i.e., the expected error rates of exactly 2.5% were found to be 2.516% and 2.515%, respectively. For the “hockey stick” NPM with 10%RSD, the maximum relative errors were 0.98% relative error in false positives (μ = 0.01, p = 0.025, q = 0.05) and 1.3% relative error in false negatives (μ = 0.03, p = 0.05, q = 0.025). These correspond to absolute errors of 0.024% and 0.032%, respectively. For the quadratic NPMs, even having 10%RSD on σ0, μ, and ρ, performance was as predicted by theory, with no relative error greater than 0.45%. With 30%RSDs on the NPM parameters, performance degrades, as would be expected, but not by very much in absolute terms. For the linear NPMs, the maximum relative errors were at μ values of 0.03, with 3.3% relative error in false positives (p = 0.025, q = 0.05) and 3.6% relative error in false negatives (p = 0.05, q = 0.025). These correspond to absolute errors of 0.081% and 0.089%, respectively. For the “hockey
E. Voigtman / Spectrochimica Acta Part B 63 (2008) 142–153
stick” NPMs, the maximum relative errors were at μ values of 0.03, with 8.9% relative error in false positives (p = 0.025, q = 0.05) and 9.5% relative error in false negatives (p = 0.05, q = 0.025). These correspond to absolute errors of 0.22% and 0.24%, respectively. For the quadratic NPMs, the maximum relative errors were 2.1% (p = 0.025, q = 0.05) and 2.7% (p = 0.05, q = 0.025), for μ = 0.05 and ρ = 0.003, corresponding to absolute errors of 0.052% and 0.068%, respectively. Finally, note that the largest absolute errors were obtained with the “hockey stick” NPMs with μ = 0.03: regardless of p and q, whenever the relevant true error rate was exactly 5%, the obtained value was about 5.3%. Further investigation of the low μ region is planned due to its importance in typical calibration situations. However, these error rates are remarkably good considering the magnitude of the NPM parameter variations. 6. Conclusions The ‘critical values of the non-central t distribution’ methodology has been shown to fail for homoscedastic, Gaussian measurement noise, even when no calibration curve is used. This means there can be no viable extension to heteroscedastic, Gaussian measurement noise since the homoscedastic regime is simply the unavoidable limiting case of the heteroscedastic regime. Despite the fact that the MARLAP agencies, and a number of other respected international standards organizations [4,6,7], have adopted versions of the ‘critical values of the non-central t distribution’ methodology, this intrinsically flawed method should be discarded. A second important conclusion is that it is not difficult to validly instantiate Currie's decision, detection and quantification limits schema for heteroscedastic, Gaussian measurement noise. This was clearly demonstrated for the three most important physically plausible NPMs and all of them performed well, even when the NPM parameters were randomly varied with 30%RSDs. This is of practical importance because the it will almost certainly never happen that the true NPM is known and its true parameter values are known as well. In this latter case, in fact, it would be possible to compute the true Currie limits immediately, without need for any calibration curve, since all the necessary information is already in the NPM and its parameters. A third conclusion is that, because true weights are almost certainly not known, the “sum of squares” minimized in the WLS process cannot be a true gamma variate, so the sample standard error about regression, swr, cannot truly be χ distributed. Furthermore, since ηw1/2 is a random variate, rather than a constant, both sd and yC can only be approximated as χ variates, with the accuracy of the approximation depending on the accuracy of the weights actually used in the WLS calibration curve processing and on how narrow the distribution of ηw1/2 is relative to that of swr. Hence, in the experimental content domain (quadrant 4), the modified non-central t distribution becomes only a useful approximation as well. Finally, if the NPM is unknown for whatever reason, it is debatable whether it makes any sense to attempt to formulaically compute Currie schema limits. It is clearly dangerous to
153
curve fit noise precision data to arbitrary model equations and then assume the true NPM is whichever model equation gives the best curve fit performance. This can lead to situations where there is no compelling reason to chose one NPM over another and thereby result in multiple sets of very different decision, detection and quantification limits [[5], pp. 42–47, 69–70]. The temptation would then be to pick the lowest limits as the “real” ones and this can even be an iterative slippery slope. This completely subverts Currie's original conception of his schema. Related to this dilemma is that of ignoring the NPM entirely and just setting the raw weights equal to the squared reciprocals of the sample standard errors of the noise precision data. This is a standard method, but it has several drawbacks. First, error rate performance, evaluated via Monte Carlo simulations not shown, tends to be rather poor unless the number of replicate measurements per standard, in the noise precision experiment, is unrealistically large. Typically, this means high double digits, or even triple digits, of replicates per standard. Second, theory is of little help because the “sum of squares” becomes an elegant infinite sum of gamma variates [[8], eqn 2.9] that is tantamount to useless in finding the PDF of swr and subsequent variates such as yC. In Part 4, Boumans' RSDB–EC methodology [9], and various simple LODs that neglect type 2 errors, are investigated and comprehensive recommendations are made with regard to both future harmonization of LOD methodology and standardized nomenclature. In the meantime, it is recommended that no further harmonization takes place, unless it involves discarding methodologies involving either prediction intervals of any kind or ‘critical values of the non-centrality parameter of the noncentral t distribution’ or both. The author thanks S. Nadarajah and S. Kotz for inspiration for the present work [10]. References [1] L.A. Currie, Limits for qualitative and quantitative determination — application to radiochemistry, Anal. Chem. 40 (1968) 586–593. [2] A. Hubaux, G. Vos, Decision and detection limits for linear calibration curves, Anal. Chem. 42 (1970) 849–855. [3] Multi-Agency Radiological Laboratory Analytical Protocols (MARLAP) Manual Volume III, Section 20A.3.1 20-54–20-58. United States agencies in MARLAP are EPA, DoD, DoE, DHS, NRC, FDA, USGS and NIST. [4] L.A. Currie, Detection: international update, and some emerging di-lemmas involving calibration, the blank, and multiple detection decisions, Chemom. Intell. Lab. Syst. 37 (1997) 151–181. [5] R.D. Gibbons, D.E. Coleman, Statistical Methods for Detection and Quantification of Environmental Contamination, first ed., John Wiley & Sons, New York, 2001. [6] ISO standard 11843, part 2. [7] Implementing council directive 96/23/EC concerning the performance of analytical methods and the interpretation of results, 2002/657/EC Commission Decision of August 12, 2002. [8] P.G. Moschopoulos, The distribution of the sum of independent gamma random variates, Ann. Inst. Stat. Math. 37 (1985) Part A 541–544. [9] P.W.J.M. Boumans, Detection limits and spectral interferences in atomic emission spectrometry, Anal. Chem. 66 (1994) 459A–467A. [10] S. Nadarajah, S. Kotz, Computation of signal-to-noise ratios, MATCH Commun. Math. Comput. Chem. 57 (2007) 105–110.