A Robust Method for Calculating Precision for Interlaboratory Studies with a Staggered-Nested Design

1. Introduction

In the case of quantitative methods, the aim of interlaboratory validation studies is to characterize method performance in terms of trueness and precision. ISO 5725-2/-3/-4 ([1], [2], [3]) provide a number of different designs allowing the evaluation of these two performance characteristics.

A particularly powerful design is the staggered-nested design described in ISO 5725-3 [2]. The simplest such design is the two-factor staggered-nested design, where each laboratory obtains three test results. For laboratory $i$, test results $y_{i1}$ and $y_{i2}$ are obtained under repeatability conditions, and $y_{i3}$ under intermediate conditions, e.g. on a different day.

The standard calculation method for precision is analysis of variance (ANOVA), preceded by outlier testing. While ANOVA, under the normal distribution assumption, is very efficient in the case of balanced designs, in the case of unbalanced designs, such as the staggered-nested design, the usual ANOVA does not use all the information available in the data. This is related to the fact that the reproducibility standard deviation is not calculated directly, but rather from the laboratory and intermediate standard deviation values; while the latter, in turn, is calculated from the standard deviation under intermediate conditions and the repeatability standard deviation. Another disadvantage of ANOVA and other conventional methods is that they very much depend on the data following a normal distribution and are thus highly sensitive to outliers. This can only be partially offset via outlier tests, since even if conspicuous values are not yet statistically significant outliers, considerable deviations can result in the determined precision data.

The advantage of robust estimators is that no outlier testing is required. For the estimation of means and standard deviation values, different robust estimators exist. The performance of a given robust estimators can be characterized via breakdown point (proportion of data that can be outliers without the estimate being affected) and efficiency (ratio of the statistical uncertainty of the estimate to that of the classical estimator under the normal distribution assumption). An overview of different breakdown point and efficiency values is provided in ISO 13528 [4].

The Q/Hampel and Q_n methods have both the highest breakdown point and the best efficiency.

The Q/Hampel method uses the Q method for the calculation of the robust reproducibility standard deviation $s_R$ and repeatability standard deviation $s_r$ together with the Hampel estimator for the calculation of the location parameter $x^{*}$ as described in ISO 13528 [4]. The theoretical basis for the Q method, including asymptotic performance and finite sample breakdown, is described in Müller et al. [5] and Uhlig [6].

The Q method is not only robust against outlying results, but also against a situation where many test results are identical, e.g. due to quantitative data on a discontinuous scale or due to rounding distortions. In such a situation other Q-like methods (e.g. the Q_n method originally introduced by Rousseeuw et al. [8] for univariate data) can fail because many pairwise differences are zero.

The Q method was introduced for the one-way random effect model and later modified – for ISO 13528 [4] – to deal with the situation that that there are many equal data, e.g., for the case where the number of significant digits is too small.

The Q/Hampel method is typically used for conventional designs, but it can also be applied for the staggered-nested design with two factors according to ISO 5725-3 [2] – in particular for estimating the intermediate standard deviation in addition to the reproducibility and repeatability standard deviation.

2. Robust Statistical Analysis of Results by Means of the Q/Hampel Method in a Staggered-Nested Design with Two Factors

For each level, the data obtained in the experiment are denoted $y_{ikl}$ (with $i$ representing factor 0, i.e. laboratory, $i=1,\dots ,p$; $k$ representing factor 1, $k=\mathrm{1,2}$; and $l$ representing the replicate, with $l$ ranging from $1$ to $n\left(k\right)$, with $n\left(1\right)=2$ and $n\left(2\right)=1$), i.e. for laboratory $i$ there are three measurement results $y_{i11}$, $y_{i12}$, $y_{i21}.$In summary, test results, grouped by laboratory, are denoted as follows:

$$\underset{Lab 1}{\underbrace{y_{111},y_{112},y_{121}}}, \underset{Lab 2}{\underbrace{y_{211},y_{212},y_{221}}},\dots ,\underset{Lab p}{\underbrace{y_{p11},y_{p12},y_{p21}}}$$

(1)

2.1. Determination of the Robust Reproducibility Standard Deviation ${\mathit{s}}_{\mathit{R}}$

Using the Q Method

The calculation relies on the use of pairwise differences within the data set and, thus, does not depend on the estimate of the mean or median.

The algorithm can be described as follows.

Based on the measurement results as structured in equation (1), the cumulative distribution function of all absolute between-laboratory differences is calculated as follows:

$$H_1\left(x\right)={\displaystyle \frac{2}{9p(p-1)}}{\sum }_{1\le i<j\le p}\left(\mathbf{I}\left\{\left|y_{i11}-y_{j11}\right|\le x\right\}+\mathbf{I}\left\{\left|y_{i11}-y_{j12}\right|\le x\right\}+\mathbf{I}\left\{\left|y_{i11}-y_{j21}\right|\le x\right\}+\mathbf{I}\left\{\left|y_{i12}-y_{j11}\right|\le x\right\}+\mathbf{I}\left\{\left|y_{i12}-y_{j12}\right|\le x\right\}+\mathbf{I}\left\{\left|y_{i12}-y_{j21}\right|\le x\right\}+\mathbf{I}\left\{\left|y_{i21}-y_{j11}\right|\le x\right\}+\mathbf{I}\left\{\left|y_{i21}-y_{j12}\right|\le x\right\}+\mathbf{I}\left\{\left|y_{i21}-y_{j21}\right|\le x\right\}\right)$$

(2)

where $\mathbf{I}\left\{\left|y_{ikl}-y_{jkl}\right|\le x\right\}=\left\{\begin{array}{cc}1& \mathrm{if}\left|y_{ikl}-y_{jkl}\right|\le x\\ 0& \mathrm{otherwise}\end{array}\right\}$ denotes the indicator function.

Discontinuity points of $H_1\left(x\right)$ are denoted

$$x_1,\dots ,x_m, \mathrm{where} x_1<x_2<\dots <x_m$$

For each positive discontinuity points $x_1,\dots ,x_m$, define

$$G_1\left(x_i\right)=\left\{\begin{array}{ccc}\mathrm{0,5}\bullet \left(H_1\left(x_i\right)+H_1\left(x_{i-1}\right)\right)& \mathrm{if}& i\ge 2\\ \mathrm{0,5}\bullet H_1\left(x_1\right)& \mathrm{if}& i=1;x_1>0\end{array}\right.$$

(3)

and let

$$G_1\left(0\right)=0$$

For each $x$ within the interval $\left[0,x_m\right]$, $G_1\left(x\right)$ is obtained by linear interpolation between discontinuity points ${0<x}_1<x_2<\dots <x_m$.

Finally, the robust reproducibility standard deviation $s_R$ is obtained as

$$s_R={\displaystyle \frac{G_1^{-1}\left(\mathrm{0,25}+\mathrm{0,75}\bullet H_1\left(0\right)\right)}{\sqrt{2}{\mathsf{\Phi }}^{-1}\left(\mathrm{0,625}+\mathrm{0,375}\bullet H_1\left(0\right)\right)}}\bullet b_p$$

(4)

where $H_1\left(0\right)$ is calculated as in equation (2) and is set equal to zero unless there are identical values in the data set.

In Equation (4), ${\mathsf{\Phi }}^{-1}\left(q\right)$ denotes the $q^{th}$ quantile of the standard normal distribution and $b_p$ denotes the correction factor corresponding to the number of laboratories $p$.

The correction factors $b_p$ were obtained via a simulation study, which will now be briefly described. In each simulation step, normally ($N\left(\mathrm{0,1}\right)$) distributed data corresponding to $p$ laboratories were generated and the robust reproducibility standard deviation $s_R$ was calculated in accordance with the Q method from the formulas given above. Taking the mean value across ${10}^6$ simulation steps – separately for each $p$ – it was possible to calculate the expected value for the reproducibility standard deviation $s_R$. For a given $p$, the correction factor $b_p$ was then obtained by taking the reciprocate of the expected value. For each value $p$ between 4 and 100, Table 1 provides the expected value for the reproducibility standard deviation $s_R$ along with the corresponding relative standard error and correction factor $b_p$.

For $p>12$, it was possible to derive a functional relationship between $p$ and the correction factor $b_p$ via nonlinear optimization. This functional relationship is given in equation (5) and presented in Figure 1 for $12<p\le 100$.

$$b_p={\left(\mathrm{0,2680}{\displaystyle \frac{1}{p^{\mathrm{2,3363}}}}+\mathrm{0,5810}{\displaystyle \frac{1}{p}}+\mathrm{0,9998}\right)}^{-1}$$

(5)

2.2. Determination of the Robust Intermediate Standard Deviation ${\mathit{s}}_{\mathit{I}\left(1\right)}$Using the Q Method

Take the case that factor 1 is day. As explained above, in the two-factor staggered-nested design, there are two results for day 1, and one for day 2. For a given laboratory, there are thus two within-laboratory differences corresponding to factor 1. Accordingly, the cumulative distribution function corresponding to factor 1 is calculated as follows:

$$H_{2,I\left(1\right)}\left(x\right)={\displaystyle \frac{1}{2p}}{\sum }_{i=1}^p\left(\mathbf{I}\left\{\left|y_{i11}-y_{i21}\right|\le x\right\}+\mathbf{I}\left\{\left|y_{i12}-y_{i21}\right|\le x\right\}\right)$$

(6)

where $\mathbf{I}\left\{\left|y_{i1l}-y_{i21}\right|\le x\right\}=\left\{\begin{array}{cc}1& \mathrm{if}\left|y_{i1l}-y_{i21}\right|\le x\\ 0& \mathrm{otherwise}\end{array}\right\}(l=\mathrm{1,2})$ denotes the indicator function.

As above, the discontinuity points of $H_{2,I\left(1\right)}\left(x\right)$ are denoted

$$x_1,\dots ,x_m, \mathrm{where} x_1<x_2<\dots <x_m$$

For each positive discontinuity points $x_1,\dots ,x_m$, define

$$G_{2,I\left(1\right)}\left(x_i\right)=\left\{\begin{array}{ccc}\mathrm{0,5}\bullet \left(H_{2,I\left(1\right)}\left(x_i\right)+H_{2,I(1}\left(x_{i-1}\right)\right)& \mathrm{if}& i\ge 2\\ \mathrm{0,5} \bullet H_{2,I(1}\left(x_1\right)& \mathrm{if}& i=1;x_1>0 \end{array}\right.$$

(7)

and let

$$G_{2,I\left(1\right)}\left(0\right)=0$$

For each $x$ within the interval $\left[0,x_m \right]$, $G_{2,I\left(1\right)}\left(x\right)$ is obtained via linear interpolation between discontinuity points ${0<x}_1<x_2<\dots <x_m$.

Finally, the robust intermediate standard deviation $s_{I\left(1\right)}$ is obtained as

$$s_{I\left(1\right)}={\displaystyle \frac{G_{2,I\left(1\right)}^{-1}\left(\mathrm{0,5}+\mathrm{0,5}\bullet H_{2,I\left(1\right)}\left(0\right)\right)}{\sqrt{2}{\mathsf{\Phi }}^{-1}\left(\mathrm{0,75}+\mathrm{0,25}\bullet H_{2,I\left(1\right)}\left(0\right)\right)}}\bullet c_p$$

(8)

where $H_{2,I\left(1\right)}\left(0\right)$ is calculated as in equation (6) and is set equal to zero unless there are identical values in the data set.

As above, ${\mathsf{\Phi }}^{-1}\left(q\right)$ denotes the $q^{th}$ quantile of the standard normal distribution.

Similarly to $b_p$ in the previous section, for a given number of laboratories $p$, the correction factor $c_p$ was calculated via a simulation study.

For each value $p$ between 4 and 100, Table 2 provides the expected value for the intermediate standard deviation $s_{I\left(1\right)}$ along with the corresponding relative standard error and correction factor $c_p$.

For $p>12$, it was possible to derive a functional relationship between the number $p$ of laboratories and the correction factor $c_p$ via nonlinear optimization. The functional relationship for $12<p\le 100$ is given in equation (9) and shown in Figure 2.

$$c_p=\left\{\begin{array}{cc}{\left(\mathrm{2,1251}{\displaystyle \frac{1}{p^{\mathrm{11,3592}}}}+\mathrm{0,3051}{\displaystyle \frac{1}{p}}+\mathrm{0,9999}\right)}^{-1}& p \mathrm{odd}\\ {\left(\mathrm{2,9723}{\displaystyle \frac{1}{p^{\mathrm{4,6860}}}}+\mathrm{0,3199}{\displaystyle \frac{1}{p}}+\mathrm{0,9998}\right)}^{-1}& p \mathrm{even}\end{array}\right\}$$

(9)

If $s_{I\left(1\right)}$results in a value greater than $s_R$, then $s_{I\left(1\right)}$ is set to $s_R$.

2.3. Determination of the Robust Repeatability Standard Deviation ${\mathit{s}}_{\mathit{r}}$

Using the Q Method

As far as repeatability is concerned, there is only one difference which can be calculated per laboratory in the two-factor staggered-nested design: namely, that corresponding to the this first level of factor 1 (e.g. the two results obtained on day 1, if factor 1 is day). Accordingly, the cumulative distribution function corresponding to repeatability precision is calculated as follows::

$$H_{2,r}\left(x\right)={\displaystyle \frac{1}{p}}{\sum }_{i=1}^p\mathbf{I}\left\{\left|y_{i11}-y_{i12}\right|\le x\right\}$$

(10)

where $\mathbf{I}\left\{\left|y_{i11}-y_{i12}\right|\le x\right\}=\left\{\begin{array}{cc}1& \mathrm{if}\left|y_{i11}-y_{i12}\right|\le x\\ 0& \mathrm{otherwise}\end{array}\right\}$ denotes the indicator function.

As above, the discontinuity points, the discontinuity points of $H_{2,r}\left(x\right)$ are denoted

$$x_1,\dots ,x_m, \mathrm{where} x_1<x_2<\dots <x_m$$

For each positive discontinuity points $x_1,\dots ,x_m$, define

$$G_{2,r}\left(x_i\right)=\left\{\begin{array}{ccc}\mathrm{0,5}\bullet \left(H_{2,r}\left(x_i\right)+H_{2,r}\left(x_{i-1}\right)\right)& \mathrm{if}& i\ge 2\\ \mathrm{0,5} \bullet H_{2,r}\left(x_1\right)& \mathrm{if}& i=1;x_1>0 \end{array}\right.$$

(11)

and let

$$G_{2,r}\left(0\right)=0$$

For each $x$ within the interval $\left[0,x_m \right]$, the function $G_{2,r}\left(x\right)$ is calculated via linear interpolation between discontinuity points ${0<x}_1<x_2<\dots <x_m$.

Finally, the robust repeatability standard deviation $s_r$ is obtained as

$$s_r={\displaystyle \frac{G_{2,r}^{-1}\left(\mathrm{0,5}+\mathrm{0,5}\bullet H_{2,r}\left(0\right)\right)}{\sqrt{2}{\mathsf{\Phi }}^{-1}\left(\mathrm{0,75}+\mathrm{0,25}\bullet H_{2,r}\left(0\right)\right)}}\bullet c_p$$

(12)

where $H_{2,r}\left(0\right)$ is calculated as in equation (10) and is equal to zero unless there are identical values in the data set.

As above, ${\mathsf{\Phi }}^{-1}\left(q\right)$ denotes the $q^{th}$ quantile of the standard normal distribution.

The correction factor $c_p$ is the same as in the previous section, see Table 2 and equation (9).

If $s_r$ is greater than $s_{I\left(1\right)}$, then $s_r$ is set equal to $s_{I\left(1\right)}$.

2.4. Determination of the Robust Mean ${\mathit{x}}^{*}$

Using the Hampel Estimator

Calculate the weighted means for each laboratory, denoted $y_1, \dots ,y_p$, i.e.

$$y_i={\displaystyle \frac{1}{4}}\left(y_{i11}+y_{i12}+{2y}_{i21}\right)$$

(13)

Calculate the robust mean, $x^{*}$, by solving the equation

$$\sum _{i=1}^p\Psi \left({\displaystyle \frac{y_i-x^{*} }{s^{*}}}\right)=0$$

(14)

where

$$\Psi \left(q\right)=\left\{\begin{array}{c}0, q\le -\mathrm{4,5}\\ -\mathrm{4,5}-q, -\mathrm{4,5}<q\le -3\\ -\mathrm{1,5}, -3<q\le \mathrm{1,5}\\ q, -\mathrm{1,5}<q\le \mathrm{1,5}\\ \mathrm{1,5}, \mathrm{1,5}<q\le 3\\ \mathrm{4,5}-q, 3<q\le \mathrm{4,5}\\ 0, q>\mathrm{4,5}\end{array}\right.$$

(15)

and

$$s^{*}=\sqrt{s_R^2-{\displaystyle \frac{1}{2}}{s}_{I\left(1\right)}^2-{\displaystyle \frac{1}{8}}{s}_r^2}$$

(16)

where $s_R$, $s_{I\left(1\right)}$ and $s_r$ denote the robust reproducibility, intermediate and repeatability standard deviations obtained in accordance with the Q method (as described in 2.1, 2.2 and 2.3), respectively.

The exact solution may be obtained in a finite number of steps (not iteratively) using the property that $\mathsf{\Psi }$ is partially linear in $x^{*}$ and by means of the interpolation nodes of the left side of equation (14) (interpreted here as a function of $x^{*}$).

The interpolation nodes are obtained as follows

- for the first value $y_1$:

$d_1=y_1-\mathrm{4,5}\bullet s^{*}$, $d_2=y_1-3\bullet s^{*}$, $d_3=y_1-\mathrm{1,5}\bullet s^{*}$, $d_4=y_1+\mathrm{1,5}\bullet s^{*}$, $d_5=y_1+3\bullet s^{*}$,

$d_6=y_1+\mathrm{4,5}\bullet s^{*}$- for the first value $y_2$:

$d_1=y_2-\mathrm{4,5}\bullet s^{*}$, $d_2=y_2-3\bullet s^{*}$, $d_3=y_2-\mathrm{1,5}\bullet s^{*}$, $d_4=y_2+\mathrm{1,5}\bullet s^{*}$, $d_5=y_2+3\bullet s^{*}$,

$d_6=y_2+\mathrm{4,5}\bullet s^{*}$- and so on for all values $y_3, \dots ,y_p$.

The notes $d_1,d_2,d_3,\dots ,d_{6\bullet p}$ are sorted in ascending order: $d_{\left\{1\right\}},d_{\left\{2\right\}},d_{\left\{3\right\}},\dots ,d_{\left\{6\bullet p\right\}}$.

For each $m=1,\dots ,(6\bullet p-1)$, the following quantity is then calculated:

$$p_m={\sum }_{i=1}^p\Psi \left({\displaystyle \frac{y_i-d_{\left\{m\right\}} }{s^{*}}}\right)$$

(17)

It is then checked whether

$p_m=0$. If so, $d_{\left\{m\right\}}$ is a solution of equation (14).

$p_{m+1}=0$. If so, $d_{\left\{m+1\right\}}$ is a solution of equation (14).

$p_m\bullet p_{m+1}<0$. If so, $x_m=d_{\left\{m\right\}}-{\displaystyle \frac{p_m}{\frac{p_{m+1}-p_m}{d_{\left\{m+1\right\}}-d_{\left\{m\right\}}}}}$ is a solution of equation (14).

Let $S$ denote the set of all solutions of equation (14).

The solution $x^{*}\in S$ nearest to the median is taken as the location parameter $x^{*}$, i.e

$$\left|x^{*}-\mathrm{median}(y_1,y_2, \dots ,y_p)\right|=\min \left\{\left|x-\mathrm{median}(y_1,y_2, \dots ,y_p)\right|;x\in S\right\}$$

(18)

Several solutions may exist. If there are two solutions nearest the median, or if there is no solution at all, the median itself is taken as the location parameter $x^{*}$.

This implementation of Hampel’s estimator has approximately 96 % efficiency for normally distributed data.

If this estimation method is used, laboratory results differing from the mean by more than 4,5 times the reproducibility standard deviation no longer have any effect on the calculation result, i.e. they are treated as outliers.

3. Conclusion

The Q/Hampel procedure described in this paper extends the range of robust statistical methods to staggered-nested designs, providing for the first time an adequate approach for handling outliers in such complex experimental formats. In conventional ISO 5725 approaches, especially under unbalanced conditions, the absence of a straightforward outlier identification process at the intermediate stage often necessitates the exclusion of all results from a laboratory, resulting in significant information loss. By contrast, the robust approach presented here allows the retention of valuable data from all laboratories, maximizing the information available for the estimation of precision parameters and means. Thus, the introduction of robust estimators in staggered-nested designs ensures both the integrity of statistical analysis and the efficient use of all available data in method validation studies.