# Can Excel® be used for gauge R&R study based on the analysis of variance (Anova)?

Alban Petitjean, Olivier Musset

► To cite this version:

Alban Petitjean, Olivier Musset. Can Excel® be used for gauge R&R study based on the analysis of variance (Anova)?. Measurement: Sensors, 2025, 42, pp.101978. 10.1016/j.measen.2025.101978. hal-05385738

HAL Id: hal-05385738

https://hal.science/hal-05385738v1

Submitted on 27 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Can Excel® be used for gauge R&R study based on the analysis of variance (Anova)?

Alban Petitjean* ID, Olivier Musset ID

Laboratoire Interdisciplinaire CARNOT, CNRS UMR, Université de Bourgogne Europe, Dijon, 6303, France

## ARTICLE INFO

### Keywords:

Excel®

Anova

Gauge R&R

Repeatability

Reproducibility

Range method

Bootstrap

## ABSTRACT

The Analysis of Variance (Anova) method, alongside the range method, is widely used to conduct Gauge Repeatability and Reproducibility (Gauge R&R) studies for measurement systems. While Anova provides more precise and robust statistical insights compared to the range-based approach, it is often perceived as more complex due to the necessity of calculating sum of squares, mean squares, and performing F-tests. Furthermore, this method can lead to misleading conclusions if the results—especially residuals—are not thoroughly and correctly analyzed. Despite its analytical complexity, the Anova method can be effectively implemented using spreadsheet tools such as Excel®. Excel® offers a flexible and accessible environment to perform rapid recalculations when input data are modified. In this paper, we demonstrate how to conduct a Gauge R&R study using the Anova approach within Excel®. The dual objective of this work is (i) to provide a user-friendly implementation within Excel®, and (ii) to enable rigorous and accurate interpretation of the statistical results produced.

## 1. Introduction

Measurement systems are employed across a wide range of industrial sectors to ensure the quality and reliability of manufacturing processes. However, the readings provided by these systems are not always dependable, due to inherent variations [1] associated with calibration, stability, repeatability, reproducibility, and linearity of the measuring instrument. In order to ensure the reliability of the measurements produced, it is therefore essential to characterize and account for these sources of variation. Such variations can be systematically analyzed and assessed through the application of statistical methodologies, such as the Gauge Repeatability and Reproducibility (Gauge R&R) study.

Gauge R&R is a statistical tool specifically designed to quantify both repeatability and reproducibility [1]. By definition [1], repeatability refers to the consistency of measurement outcomes under a uniform set of conditions within a controlled environment. In statistical terms, repeatability can be regarded as a manifestation of random error [1]. Within the context of a Gauge R&R analysis, this type of variation—termed equipment variation, or “intra-system variation” [2]—is observed when a single operator performs multiple measurements on the same part using the same instrument [3]. Reproducibility refers to the consistency of measurements obtained under a set of measurement

conditions across different environments [1]. Reproducibility is thus defined as the variability between operators [2]. Consequently, reproducibility is considered as the average variation across measurement conditions or systems [2]. In other words, within a Gauge R&R framework, a measurement depends on several parameters governed by a directional random effects model, which implies that each parameter exerts an influence on the measurement outcome [4]. The aim of this statistical study is to conduct measurements which, when grouped according to the conditions defined by repeatability and reproducibility, enable the determination of the uncertainty associated with repeatability and reproducibility.

The determination of these uncertainties is of paramount importance when this approach is part of a continuous improvement methodology [5], where the objective is to enhance the accuracy of the measurement system, as is the case in a Six Sigma initiative [6]. Six Sigma is a widely adopted methodology aimed at reducing defects and, consequently, improving the quality of a manufacturing or assembly process by minimizing process variation [6,7], particularly addressing the variables within the subprocesses [8]. The reduction of process variation relies on understanding the relative contributions of various input variables to the key performance indicators of the processes [9]. Therefore, it is essential to distinguish between process variations and those of the

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: alban.petitjean@u-bourgogne.fr (A. Petitjean).</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Block diagram for creating an Anova method.</figcaption></figure>

measurement system [10], by calculating the uncertainty associated with repeatability, reproducibility, and the measurement process variation.

Uncertainty is generally calculated using one of two methods [11]: either by evaluating the variance of the measurements through the Anova method (an acronym for Analysis of Variance), or by assessing the range of the measurements using the range method [5]. The results of these two methods can differ significantly [12]. In the literature, the Anova method is generally regarded as more accurate and robust than the range method, as it accounts for the interaction between the operator and the measured object [3,13–15]. However, residual analysis, which serves as an important diagnostic indicator in the Anova method [16], is often overlooked; when performed and the residuals are not normally distributed [12], the performance of the Anova method is reduced, and the range method may then provide more reliable results. In other words, the choice between these two methods for conducting an honest Gauge R&R analysis often depends on the evaluator's experience and training.

In cases where the evaluator is inexperienced, the range method is easily accessible [17], typically presented as a predefined formula in

tabular format [18,19], or even as a simple application [20]. Frequently, when the Anova method is used by a novice user [21], the calculations are either double-checked or compared [22] with the output of a statistical software program. The objective of this study is to demonstrate how to use Microsoft Excel® to conduct a Gauge R&R analysis. The study presents a pedagogical approach through a concrete example drawn from geochemistry. This methodological framework incorporates a practical case study to expose students to challenges that mirror real-world operational contexts. The intended audience includes undergraduate students in engineering and applied sciences, as well as professionals undergoing training in quality control and measurement systems analysis. These individuals are expected to develop a thorough understanding of the assumptions underlying the Gauge R&R methodology and to be capable of applying them across a range of practical scenarios. The advantage of using Excel is that when the measurement data are modified, the statistical calculations are automatically updated and the associated charts are redrawn. This allows the user to observe and analyze measurement system variability in real time.
<caption>**Table 1**</caption>

Measured resistance values, grouped according to the measurement series, operator are represented in the Excel® sheet for inputting measurements.

<table><thead><tr><th rowspan="2">Measurement series</th><th colspan="3">Operator 1</th><th colspan="3">Operator 2</th><th colspan="3">Operator 3</th></tr><tr><th>Measure 1</th><th>Measure 2</th><th>Measure 3</th><th>Measure 1</th><th>Measure 2</th><th>Measure 3</th><th>Measure 1</th><th>Measure 2</th><th>Measure 3</th></tr></thead><tbody><tr><td>1</td><td>55.23</td><td>55.20</td><td>55.27</td><td>55.28</td><td>55.30</td><td>55.34</td><td>55.31</td><td>55.33</td><td>55.29</td></tr><tr><td>2</td><td>55.81</td><td>55.79</td><td>55.80</td><td>55.83</td><td>55.88</td><td>55.82</td><td>55.77</td><td>55.86</td><td>55.79</td></tr><tr><td>3</td><td>55.78</td><td>55.70</td><td>55.74</td><td>55.67</td><td>55.72</td><td>55.71</td><td>55.72</td><td>55.78</td><td>55.72</td></tr><tr><td>4</td><td>55.55</td><td>55.49</td><td>55.42</td><td>55.47</td><td>55.50</td><td>55.55</td><td>55.43</td><td>55.42</td><td>55.53</td></tr><tr><td>5</td><td>56.23</td><td>56.17</td><td>56.23</td><td>56.18</td><td>56.25</td><td>56.27</td><td>56.19</td><td>56.17</td><td>56.26</td></tr><tr><td>6</td><td>56.48</td><td>56.35</td><td>56.39</td><td>56.35</td><td>56.34</td><td>56.28</td><td>56.26</td><td>56.35</td><td>56.35</td></tr><tr><td>7</td><td>56.06</td><td>56.07</td><td>55.94</td><td>56.07</td><td>56.14</td><td>56.07</td><td>56.07</td><td>56.05</td><td>56.01</td></tr><tr><td>8</td><td>56.19</td><td>56.16</td><td>56.11</td><td>56.07</td><td>56.15</td><td>56.16</td><td>56.20</td><td>56.14</td><td>56.12</td></tr><tr><td>9</td><td>55.91</td><td>55.87</td><td>55.91</td><td>55.90</td><td>55.89</td><td>55.84</td><td>55.82</td><td>55.84</td><td>55.75</td></tr><tr><td>10</td><td>55.62</td><td>55.65</td><td>55.70</td><td>55.68</td><td>55.74</td><td>55.68</td><td>55.73</td><td>55.68</td><td>55.73</td></tr></tbody></table>

## 2. Material and method

### 2.1. Mathematical model

The mathematical model of the Anova method [23–25] for a Gauge R&R study assumes that a measurement is a random variable that depends on three parameters: (i) the total number of measurements for a given series, denoted as $nb_m$, (ii) multiple series of measurements, with the total number referred to as $nb_s$, and (iii) several operators, with $nb_o$ representing the total number of distinct operators. In this case, the total number of measurements is $nb_m \cdot nb_o \cdot nb_s$.

Accordingly, a measurand $y$ [1] is denoted as $(y_{ijk})$, representing the k-th measurement taken by operator $j$ for the i-th measurement series. All indices fall within the following intervals (Equation (1)):

$$ \begin{cases} i \in [1, nb_s] \\ j \in [1, nb_o] \\ k \in [1, nb_m] \end{cases} \quad 1 $$

The probabilistic model associated with the measurand $(y_{ijk})$ [26] assumes a directional random effects model [4] meaning that it is necessary to consider the measurand $(y_{ijk})$ as the sum of six contributions: the true value, denoted $x$, the measurement instrument bias, referred to as $\mu$, the dispersion due to series $i$, denoted $\alpha$, the effect of operator $j$ on the measurement (reproducibility), denoted $\beta_j$, and the error due to the measurement instrument (repeatability), denoted $\epsilon_{ijk}$. In mathematical terms, the following Equation (2) is obtained:

$$ y_{ijk} = x + \mu + \alpha_i + \beta_j + \alpha\beta_{ij} + \epsilon_{ijk} \quad 2 $$

Equation (3) corresponds to the calculation of the variance, denoted Var, from Equation (2). Furthermore, the following assumptions are made.

1. $x, \mu$ are constants.

2. $(\alpha_i)$ are normal random variables with a mean of zero and variance $\sigma_a^2$ (where $\sigma$ is the standard deviation).

3. $(\beta_j)$ are normal random variables with a mean of zero and variance $\sigma_\beta^2$.

4. $(\alpha\beta_{ij})$ are normal random variables with a mean of zero and variance $\sigma_{\alpha\beta}^2$.

5. $(\epsilon_{ijk})$ are normal random variables with a mean of zero and variance $\sigma_\epsilon^2$.

6. $(\alpha_i), (\beta_j), (\alpha\beta_{ij})$ and $(\epsilon_{ijk})$ are mutually independent random variables.

Equation (3) corresponds to the calculation of the variance, denoted as Var, from Equation (2).

$$ \mathrm{Var}(y_{ijk}) = \sigma_{\alpha}^{2} + \sigma_{\beta}^{2} + \sigma_{\alpha\beta}^{2} + \sigma_{\epsilon}^{2} \quad 3 $$

Referring back to the definition [1,2] of repeatability, with the standard deviation denoted as $\sigma_r$, reproducibility, with the standard deviation denoted as $\sigma_R$, and Gauge R&R, with the standard deviation denoted as $\sigma_{R&R}$, Equation (3) becomes Equation (4):

$$ \mathrm{Var}(y_{ijk}) = \sigma_{\alpha}^{2} + \sigma_{R\&R}^{2} = \sigma_{\alpha}^{2} + \sigma_{R}^{2} + \sigma_{\epsilon}^{2} \quad \begin{cases} \sigma_{\epsilon} = \sigma_{\alpha} \\ \sigma_{R} = \sqrt{\sigma_{\beta}^{2} + \sigma_{\alpha\beta}^{2}} \\ \sigma_{R\&R} = \sqrt{\sigma_{R}^{2} + \sigma_{\epsilon}^{2}} \end{cases} \quad 4 $$

Fig. 1 outlines the general procedure for performing an Anova method and determining the uncertainty associated with repeatability and reproducibility. Each step of the procedure described in Fig. 1 is explained in the following sections.

### 2.2. Measurement

To determine the standard deviation associated with repeatability and reproducibility, while adhering to the assumptions of the mathematical model (see the Mathematical Model section), it is essential to carefully select the number of measurements per series ($nb_m$), the number of operators ($nb_o$), and the number of series ($nb_s$). If these three numbers are sufficiently large, typically around ten, the central limit theorem [27,28] and the law of large numbers [29] ensure the normality of the estimates for the linear regression parameters. Furthermore, many measurement instruments are computer-controlled and typically average multiple readings to again ensure the normality of the measurements.

The study [30] recommends, in order to ensure the normality of the measurements, using 3 operators, 10 measurement series, and 10 measurements per series. Depending on the study, these numbers may vary. For instance, the study [2] chooses 3 measurements per series. In Ref. [31], the number of operators, which is 3, does not allow a conclusive judgment on whether the series is erroneous or correct; it suggests doubling the number of operators ($nb_o = 6$).

#### 2.2.1. Resistance measurement

A first example employs a protocol consisting of 10 measurement series performed by 3 operators, with 3 measurements taken per series. The measurement device used in this example is a homemade ohmmeter. The objects measured are ten resistors, numbered from 1 to 10, each with a reference value provided by the manufacturer of 56 ± 2 Ω. Using the notation from the Mathematical Model section, the number of series $nb_s$ is 10, the number of operators $nb_o$ is 3, and the number of measurements per series $nb_m$ is 3. The total number of measurements, which is the product of the number of series, measurements, and operators $nb_m \cdot nb_o \cdot nb_s$, equals 90. All the results from these measurements are provided in Table 1.
<figure><img src="image_4.png" /><figcaption>Fig. 2. Excel® sheet for inputting measurements.</figcaption></figure>

### 2.2.2. Isotopic composition

Often, in a Gauge R&R study, other technical constraints require reducing the number of measurements. For example, in Ref. [24], the number of measurements must be determined based on the financial cost of the measurement. As the preparation time for the measurement is a significant constraint in this study [24], it conducted a small number of measurement series with few operators.

This reduction in the number of measurements can cause problems when the assumptions of the mathematical model are no longer met. This non-normality may also arise from outlier measurements [32] related to the measurement process, the operator, but most notably from the heterogeneity of the measured object. However, there is no consensus on how to filter out outlier values [32], particularly in the case of destructive measurements.

Destructive measurements alter the measured object. Therefore, the measured object is no longer identical over time. As such, the assumption of a time-invariant measured object is false [33]. In other words, during a destructive measurement, variations in the measurements are confounded with variations in the object being measured, and there is no systematic statistical solution to address this issue [33]. For example, if the internal variability of the measured object is significant, this could lead to an overestimation of repeatability [34]. Conversely, if the internal variability of the measured object is underestimated and does not reflect true variations, repeatability would be underestimated [6].

To ensure a reliable Gauge R&R, the choice of the measured object is crucial, as the effects of disruptive sources of variation related to the sample must be negligible [33]. One of the conditions of this protocol is that each operator performs the measurement twice on the same object [35]. In this case, the number of measured objects should be almost identical to the number of operators [36] in order to minimize the effects of non-random sampling. This reduction in the number of measurements presents a major challenge: the measurements must be representative of both the measured object and the measurement system [37,38]. These issues are encountered in the field of geochemistry, particularly when the reference material exhibits inhomogeneity [20]. Geochemistry focuses on the chemical composition of samples, incorporating both major and trace elements. With regard to trace elements, concentrations are generally expressed in parts per million (ppm). In addition to the measurement of elemental concentrations, it is frequently imperative to analyze their isotopic compositions, whether for major or trace elements. Isotopic measurements are utilized to determine the relative

abundance of heavy versus light isotopes in a given sample, typically expressed in per mil (‰). The present study focuses on carbon as the element of interest. Specifically, we consider two stable isotopes: ¹³C (a minor isotope, ~1.11‰ abundance) and ¹²C (the predominant stable isotope, ~98.89‰) [39]. The isotopic composition is determined through the measurement of the carbon isotopic ratio, denoted as δ¹³C (‰), which represents the standardized ratio of ¹³C-¹²C in the sample.

The experimental determination of δ¹³C involves two primary steps. Initially, CO₂ gas is produced from a carbonate sample, generally containing calcium. Secondly, the isotopic composition of the CO₂ is analyzed using spectroscopic techniques. Since the mid-20th century [40], the standard method for CO₂ production has relied on the utilization of orthophosphoric acid for the purpose of acid digestion, followed by purification of the resulting gas. This CO₂ is known to reflect the isotopic composition of the original solid sample. The gas is then analyzed either by mass spectrometry using a magnetic sector or by optical absorption spectroscopy in the infrared range. As an alternative to acid digestion, CO₂ can also be produced via laser-induced calcination [37,38,41,42].

Three operators measure the carbon isotopic composition of three samples of calcium carbonate derived from Comblanchien limestone, a homogeneous sedimentary rock originating from the Bourgogne region in France (latitude 47.10431, longitude 4.90692). This limestone is frequently used in our analytical procedures as a reference due to its minimal isotopic variability. Each operator analyses the δ¹³C values for three sub-samples, performing three replicate measurements per sample. The procedures are as follows: (i) operator 1 uses manual acid digestion for CO₂ production, followed by mass spectrometric analysis (dual-inlet IMRS, MultiCarb-GV Isoprime instrument from Elementar), (ii) operator 2 also employes acid digestion, but analyzes the CO₂ using optical spectroscopy (Delta Ray, ThermoFisher) and (iii) operator 3 utilizes laser-induced calcination to generate CO₂, which is then analyzed via optical spectroscopy (Delta Ray, ThermoFisher). These measurements are represented in Fig. 2 as an Excel® spreadsheet.

## 3. Result

### 3.1. Case 1: resistance measurement

Graphical tools help improve the effectiveness of a Gauge R&R analysis [43] by providing supplementary information related to the
<figure><img src="image_5.png" /><figcaption>Fig. 3. Box and whisker plots for the distribution of electrical resistance measurements. In the legend, the acronyms Op and Obj mean respectively Operator and series of the Object measured.</figcaption></figure>

assumptions of the mathematical model (see Mathematical Model section). To be effective, these graphical tools should be applied at the beginning of the Gauge R&R analysis in order to identify potential weaknesses in the measurement system and assess their causes.

### 3.1.1. Box and whisker plot

The boxplot makes it possible to (i) identify outlier measurements and (ii) characterize the symmetry of the measurement distribution [43]. The boxplot should be applied to homogeneous measurements, meaning that both the operator and the measured object must be the same: the measurements ($y_{ijk}$) are grouped according to operator $j$ and measurement series $i$. Due to these groupings, issues may arise related to

an insufficient number of measurements.

To create a box and whisker plot in Excel®, simply.

1. Select the grouped measurements

2. Select the **Insert** tab, and then select **Insert Statistical Chart**, and under **Box and Whisker**, pick a chart.

Fig. 3 presents the Box and Whisker plots for each measured object—initially in their original sequence. This visualization provides insights comparable to those offered by the scatter plot recommended in Ref. [2]. Consequently, both types of graphs can be displayed and interpreted together [43].

<figure><img src="image_6.png" /><figcaption>Fig. 4. Modified GRR diagrams for individual operators.</figcaption></figure>
The analysis of the Box and Whisker plots (Fig. 3) reveals that all measurements are contained within the range of [55.1; 56.4] Ω, which falls within the manufacturer's specified uncertainty of 56 ± 2 Ω. However, some measurement distributions, such as that of object 6, appear to be asymmetric. This graph thus highlights the potential for measurement errors that may not be apparent in this example. To ascertain whether these errors are systematic or random, further analyses are required to complement the current findings.

### 3.1.2. Modified GRR diagram

Instrumental drift [1,2], defined as the change in bias over time, must be evaluated as a priority. Indeed, before conducting a Gauge R&R, it is necessary to verify that the measurement system's properties do not change over time (i.e., the measurement system is stable and consistent). If this is not the case, the cause of the issue must be identified, as the measuring instrument is considered unreliable [1].

The Modified GRR diagram allows for easy verification of this property and also enables the analysis of the operator's effect on the measured object [44]. To create this plot, the following steps must be taken: (i) calculate the average of all measurements per series, and (ii) calculate the error, which is defined as the difference between each individual measurement and the previously calculated mean. In other words, the error, which is the endogenous variable $e_{ijk}$ as defined by Equation (5), is represented as a function of the exogenous variable $\bar{y}_i$ (Equation (5)).

$$ \left\{ \begin{array}{l} e_{ijk} = y_{ijk} - \bar{y}_i \\ \bar{y}_i = \frac{1}{nb_o \cdot nb_m} \sum_{k=1}^{nb_m} \sum_{j=1}^{nb_m} y_{ijk} \end{array} \right. \qquad (5) $$

To create this representation in Excel®, follow these steps.

1. Select the average $\bar{y}_i$ using the **Average** function on the grouped measurements by all measurements and operators.

2. Calculate the difference between the measurement $y_{ijk}$ and the average measurement $\bar{y}_i$.

3. Select the **Insert** tab, then choose **Insert Scatter (X, Y) or Bubble Chart**, and under **Scatter**, pick a chart.

4. Perform a linear regression: select the points on the graph, right-click, and press **Add Trendline**.

5. Add the average errors $e_{ijk}$ for each operator to the graph.

Fig. 4 shows this Modified GRR diagram. This representation can also be done in the form of a histogram [2] or an error chart [2]. While the analysis of these charts is similar, some studies prefer to visualize all of these charts together [43].

The x-axis of Fig. 4 provides information about the uniformity of the measurement system. Half of the measurements are distributed within

<table><caption>**Table 2**

Coefficient of direction a and y-intercept b obtained for the different operators.</caption><thead><tr><td>Operator</td><td>Slope a</td><td>y-intercept b</td></tr></thead><tbody><tr><td>1</td><td>0.89</td><td>-1.45</td></tr><tr><td>2</td><td>2.46</td><td>-1.45</td></tr><tr><td>3</td><td>0.23</td><td>-0.57</td></tr></tbody></table>

<table><caption>**Table 3**

Coefficient of direction a and y-intercept b obtained for the different operators. $t_1$ and $t_2$ values calculated from equation (7) and the data in Table 1.</caption><thead><tr><td>Operator</td><td>t<sub>1</sub></td><td>t<sub>2</sub></td></tr></thead><tbody><tr><td>1</td><td>5.50·10<sup>-13</sup></td><td>1.53·10<sup>-25</sup></td></tr><tr><td>2</td><td>7.81·10<sup>-15</sup></td><td>6.65·10<sup>-28</sup></td></tr><tr><td>3</td><td>1.22·10<sup>-14</sup></td><td>1.80·10<sup>-27</sup></td></tr></tbody></table>

the interval [55.60; 55.95] Ω. As a result, the majority of the resistances are distributed within 30% of the total measurement range. Fig. 4 also allows for the evaluation of the variability of individual measurements from a given operator and their comparison to the overall mean. By using the line representing the means of the product samples, it is possible to assess whether some operators exhibit larger deviations from these means. These deviations represent measurement biases [1]. For all three operators, these curves lie within an error of less than 1% of the resistance value, i.e., within ±0.05 Ω. Operator 1 shows the maximum variation, as this operator measured the two extreme values; the linear regression curve obtained for the entire set of values is increasing. The linear regression curves for Operators 2 and 3 are slightly decreasing and close to the average value. In other words, if there is any measurement bias, it appears to be constant.

To quantify these analyses, the manual [2] suggests using statistical tests on the linear regression curves, specifically on the slope coefficient a and the y-intercept b of equation (6). The number of measured objects and the number of measurements per series for conducting this study should be greater than or equal to 5 and 10, respectively [2].

$$ \exists j \in [1, nb_o], e_{ijk} = a \cdot \bar{y}_i + b \qquad 6 $$

These coefficients are given for the different operators in Table 2.

The first test [2] has the hypothesis that the slope coefficient a is zero. To verify this hypothesis, the quantity $t_1$ calculated according to equation (7) must be less than or equal to the Student's t-test value.

The second test [2] aims to verify if the y-intercept b or the instrument bias is zero. In this case, the quantity $t_2$ defined by equation (7) must be less than or equal to the Student's t-test value.

To use these tests, the first test must be performed first, i.e., verify that the bias is constant, and then the second test must be performed, which means checking if the bias is zero.

$$ \left\{ \begin{array}{l} t_1 = \frac{\frac{|a|}{s}}{\sqrt{\sum_{i=1}^{nb_s} \left( \bar{y}_i - \frac{1}{nb_s} \sum_{i=1}^{nb_s} \bar{y}_i \right)^2}} \\ t_2 = \frac{|b|}{s \sqrt{\frac{1}{nb_s \cdot nb_m} + \frac{\left( \frac{1}{nb_s} \sum_{i=1}^{nb_s} \bar{y}_i \right)^2}{\sum_{i=1}^{nb_s} \left( \bar{y}_i - \frac{1}{nb_s} \sum_{i=1}^{nb_s} \bar{y}_i \right)^2}}} \quad \text{where } s = \sqrt{\frac{\sum_{i=0}^{nb_s} e_{ijk}^2 - b \cdot \sum_{i=0}^{nb_s} e_{ijk} - a \cdot \sum_{i=0}^{nb_s} e_{ijk} \bar{y}_i}{nb_s \cdot nb_m - 2}} \end{array} \right. $$
<figure><img src="image_7.png" /><figcaption>Fig. 5. Range control chart.</figcaption></figure>

The formulas of equation (7) are easily implementable in Excel® using the formulas: SUM, SQRT, and AVERAGE. Table 3 provides the values of the quantities t₁ and t₂ for the three operators.

The value of the Student's t-test is obtained using the function T.INV,

where the first argument is defined as $1 - \alpha/2$ and the second is the degrees of freedom, which equals $nb_s \cdot nb_m - 2$. The parameter $\alpha$ represents the decision threshold, typically set at 5%. If the decision threshold is 5%, the value of the Student's t-test for a total number of measurements $nb_m$ and total number of series $nb_s$ equal to 3 is 2.05.

According to Table 3, the theoretical value of the Student's t-test is,

<figure><img src="image_8.png" /><figcaption>Fig. 6. Average control chart.</figcaption></figure>
for the different operators, higher than the value of the quantity $t_1$ or the quantity $t_2$. In other words, there is no bias dependent on the measured objects or, respectively, there is no operator that makes biased measurements. In the critical case where an operator seems to be a possible cause of unexpected measurement results, it is preferable to perform a new analysis of the repeatability and reproducibility by excluding the measurements taken by this operator [45].

### 3.1.3. Control chart

#### 3.1.3.1. Range control chart.
The modified Range control chart [43] allows for checking the variability of repeated measurements. This chart represents the average range $\bar{R}_{ij}$ associated with each operator as a function of the average of the measured objects (Equation (8)).

$$ \bar{R}_{ij} = \max_{k \in [1, nb_m]} (y_{ijk}) - \min_{k \in [1, nb_m]} (y_{ijk}) \qquad 8 $$

To implement Equation (8) in Excel®, it is necessary to select the relevant range of measurements and use the MAX and MIN functions.

The average range $\bar{R}_{ij}$ must be compared to the central line (CL) average, calculated using the AVERAGE function in Excel®, as well as the lower control limit (LCL) and upper control limit (UCL), as defined in Equation (9).

$$ \left\{ \begin{array}{l} CL = \bar{R}_i = \frac{1}{nb_o} \sum_{j=1}^{nb_o} \bar{R}_{ij} \\ LCL = D_3 \cdot CL \\ UCL = D_4 \cdot CL \end{array} \right. \qquad 9 $$

Where $D_3$ and $D_4$ are tabulated functions calculated based on the number $nb_m$ [46]. For this study, these constants are respectively 0 and 2.574. Fig. 5, generated by executing the following instructions: Select the Insert tab, then select Insert Scatter (X, Y) or Bubble Chart, and under Scatter, pick a chart, represents the range control chart.

The analysis of Fig. 5 shows that the system under study is under control concerning the variability of repeated measurements because all the ranges are within the LCL (Lower Control Limit) and UCL (Upper Control Limit).

#### 3.1.3.2. Average control chart.
The average control chart allows for the evaluation of the variability of the measured objects by calculating the average $\bar{y}_{ij}$ of the objects measured by each operator (Equation (10)).

$$ \bar{y}_{ij} = \frac{1}{nb_m} \sum_{k=1}^{nb_m} y_{ijk} \qquad 10 $$

Equation (10) can be implemented in Excel® by selecting the relevant range of measurements and using the AVERAGE function. The central line CL average of the operators (Equation (11)) is also displayed on the average control chart.

The averages $\bar{y}_{ij}$ obtained must be compared to the lower control limit (LCL) and upper control limit (UCL), as defined in Equation (11).

$$ \left\{ \begin{array}{l} CL = \frac{1}{nb_o} \sum_{j=1}^{nb_o} \bar{y}_{ij} \\ LCL = \bar{y}_{ij} - A_2 \cdot \bar{R}_i \\ UCL = \bar{y}_{ij} + A_2 \cdot \bar{R}_i \end{array} \right. \qquad 11 $$

Where $A_2$ is a tabulated function calculated based on the number $nb_m$ [46] and $\bar{R}_i$ is calculated using Equation (10). To display the average control chart (Fig. 6), the following steps should be followed: Select $\bar{y}_i$ and $\bar{y}_{ij}$, go to the Insert tab, then select Insert Scatter (X, Y) or Bubble Chart, and under Scatter, choose a chart.

As with the range control chart, all the means $\bar{y}_{ij}$ fall within the control limits. However, only 50 % of the measurements fall within the control limits. Provided that the measured object covers the entire variability of the process, this measurement system is considered acceptable for evaluating the variability between the measured samples if at least half of the plotted means lie outside the control limits [2]. Therefore, the measurements seem to follow the assumptions of the theoretical model of the Anova method (Section: Mathematical Model); thus, it remains to apply this method.

### 3.1.4. Computation for anova

The graphical tools presented in the Graphical Tools section allow for verification of the assumptions outlined in the Mathematical Model section. At this stage, the Anova method must be implemented, which—as the name suggests—is based on an extension of the variance calculation. Accordingly, several sums of squares must be computed for the three parameters involved: the measurement, the operator, and the measurement series. The sum of squares associated with each parameter is denoted as follows: SCS for the series, SCO for the operators, SCOS for the interaction between operators and series, SCE for the error, and TCS for the total (Equations (12)–(16)). An equivalent Excel® formula is provided for each of these sums of squares—SCS, SCO, SCOS, SCE and TCS—as long as a range of values is defined, noted $P_{V_u}$, where V represents the quantity and u denotes the corresponding index (i.e., operator, number of measurements, or measurement series).

$$ \left\{ \begin{array}{l} SCS = \frac{1}{nb_o \cdot nb_m} \sum_{i=1}^{nb_o} \left( \sum_{k=1}^{nb_m} \sum_{j=1}^{nb_o} y_{ijk} \right)^2 - \frac{1}{nb_m \cdot nb_o \cdot nb_s} \left( \sum_{i=1}^{nb_s} \sum_{j=1}^{nb_o} \sum_{k=1}^{nb_m} y_{ijk} \right)^2 \\ \quad \forall i \in [1, nb_s], S_{SCS_i} = \text{POWER}(\text{SUM}(P_{y_{sbi}}), 2) \\ SCS_{Excel} = \text{SUM}(P_{S_{SCM_i}}) / (P_{nb_o} \cdot P_{nb_m}) - \text{SUM}(P_{y_{sbi}}) \hat{2} / (P_{nb_o} \cdot P_{nb_m} \cdot P_{nb_s}) \end{array} \right. \qquad 12 $$

$$ \left\{ \begin{array}{l} SCO = \frac{1}{nb_m \cdot nb_s} \sum_{j=1}^{nb_o} \left( \sum_{k=1}^{nb_m} \sum_{i=1}^{nb_s} y_{ijk} \right)^2 - \frac{1}{nb_m \cdot nb_o \cdot nb_s} \left( \sum_{i=1}^{nb_s} \sum_{j=1}^{nb_o} \sum_{k=1}^{nb_m} y_{ijk} \right)^2 \\ \quad \forall j \in [1, nb_o], S_{SCO_j} = \text{POWER}(\text{SUM}(P_{y_{sbi}}), 2) \\ SCO_{Excel} = \text{SUM}(P_{SSCO_j}) / (P_{nb_m} \cdot P_{nb_s}) - \text{SUM}(P_{y_{sbi}}) \hat{2} / (P_{nb_o} \cdot P_{nb_m} \cdot P_{nb_s}) \end{array} \right. \qquad 13 $$

$$ \left\{ \begin{array}{l} SCOS = \frac{1}{nb_m} \sum_{i=1}^{nb_s} \sum_{j=1}^{nb_o} \left( \sum_{k=1}^{nb_m} y_{ijk} \right)^2 - \frac{1}{nb_m \cdot nb_o \cdot nb_s} \left( \sum_{i=1}^{nb_s} \sum_{j=1}^{nb_o} \sum_{k=1}^{nb_m} y_{ijk} \right)^2 - SSM - SCO \\ \quad \forall (i, j) \in [1, nb_s] \times [1, nb_o], S_{SCOS_{ij}} = \text{POWER}(\text{SUM}(P_{y_{sbi}}), 2) \\ SCOS_{Excel} = \text{SUM}(P_{SSCOS_{ij}}) / (P_{nb_m}) - \text{SUM}(P_{y_{sbi}}) \hat{2} / (P_{nb_o} \cdot P_{nb_m} \cdot P_{nb_s}) - P_{SSM} - P_{SCO} \end{array} \right. \qquad 14 $$
<caption>**Table 4**</caption>

New values obtained for the degrees of freedom, mean squares, the experimental Fisher test value, and the Fisher-Snedecor critical value when $F_{\text{expe}, SCOM}$ is smaller than $F_{F-S, SCOM}$.

<table><thead><tr><th></th><th>deg'</th><th>MS'</th><th>F<sub>F-S</sub>'</th><th>F<sub>expe'</sub>'</th></tr></thead><tbody><tr><td>SCS</td><td>deg<sub>SCS</sub></td><td>MSCS</td><td>F<sub>Fisher-Snedecor</sub>(α, deg<sub>SCS</sub>, deg'<sub>SCE</sub>)</td><td>MSCS<br>MSCE</td></tr><tr><td>SCO</td><td>deg<sub>SCO</sub></td><td>MSCO</td><td>F<sub>Fisher-Snedecor</sub>(α, deg<sub>SCO</sub>, deg'<sub>SCE</sub>)</td><td>MSCO<br>MSCE</td></tr><tr><td>SCOS</td><td>deg<sub>SCOS</sub> +<br>deg<sub>SCE</sub></td><td>0</td><td>F<sub>Fisher-Snedecor</sub>(α, deg<sub>SCOS</sub>,<br>deg<sub>SCE</sub>)</td><td>0</td></tr><tr><td>SCE</td><td>deg<sub>SCOS</sub> +<br>deg<sub>SCE</sub></td><td>SCOS + SCE<br>deg<sub>SCE</sub></td><td></td><td></td></tr></tbody></table>

$$ \left\{ \begin{aligned} \text{TCS} &= \sum_{i=1}^{nb_m} \sum_{j=1}^{nb_a} \sum_{k=1}^{nb_s} (y_{ijk})^2 - \frac{1}{nb_m \cdot nb_a \cdot nb_s} \left( \sum_{i=1}^{nb_m} \sum_{j=1}^{nb_a} \sum_{k=1}^{nb_s} y_{ijk} \right)^2 \\ \text{TCS}_{\text{Excel}} &= \text{SUMSQ}(P_{y_{ijk}}) - \text{SUM}(P_{y_{ijk}}) \frac{2}{P_{nb_a} \cdot P_{nb_m} \cdot P_{nb_s}} \end{aligned} \right. \quad 15 $$

$$ \text{SCE} = \text{TCS} - \text{SCS} - \text{SCO} - \text{SCOS} \quad 16 $$

From the sum of squares, the next step in calculating the variance—i. e., the mean square (using the same notation as the sum of squares, but with a preceding *M*)—is to determine the degrees of freedom, denoted *deg* (Equations (17)–(21)). The mean square is then obtained by dividing each sum of squares by its corresponding degrees of freedom (Equations (17)–(21)).

$$ \text{deg}_{\text{SCS}} = \text{nb}_s - 1 \Rightarrow \text{MSCS} = \frac{\text{SCS}}{\text{deg}_{\text{SCS}}} \quad 17 $$

$$ \text{deg}_{\text{SCO}} = \text{nb}_o - 1 \Rightarrow \text{MSCO} = \frac{\text{SCO}}{\text{deg}_{\text{SCO}}} \quad 18 $$

$$ \text{deg}_{\text{SCOS}} = \text{deg}_{\text{SCS}} \cdot \text{deg}_{\text{SCO}} \Rightarrow \text{MSCOS} = \frac{\text{SCOS}}{\text{deg}_{\text{SCOS}}} \quad 19 $$

$$ \text{deg}_{\text{SCE}} = \text{nb}_o \cdot \text{nb}_s \cdot (\text{nb}_m - 1) \Rightarrow \text{MSCE} = \frac{\text{SCE}}{\text{deg}_{\text{SCE}}} \quad 20 $$

$$ \text{deg}_{\text{TCS}} = \text{nb}_o \cdot \text{nb}_m \cdot \text{nb}_s - 1 \Rightarrow \text{MTCS} = \frac{\text{TCS}}{\text{deg}_{\text{TCS}}} \quad 21 $$

Formulas 17-21 are easily implementable in Excel®.

To compare the variances of two groups assumed to follow normal distributions (as per the model assumptions), the appropriate test is the Fisher test. This test is based on a comparison between the experimental value, denoted $F_{\text{expe}}$, and the critical value given by the Fisher-Snedecor distribution table at a 5% significance level, denoted $F_{F-S}$ (Equations (22)–(24)). The $F_{\text{Fisher-Snedecor}}$ value used to obtain $F_{F-S}$ can be computed in Excel® using the FINV function, where the first argument is the risk level $\alpha$, the second is the numerator degrees of freedom of the calculated quantity (SCS, SCO or SCOS), and the third is the denominator degrees of freedom of SCE. The threshold $\alpha$ is generally set at 5% [2].

<caption>**Table 5**</caption>

Results of the Anova method performed on Table 1.

<table><thead><tr><th>Error source</th><th>Sum of squares</th><th>Degrees of freedom</th><th>Average square</th><th>Experimental Fisher test</th><th>Fisher test</th><th>Risk</th></tr></thead><tbody><tr><td>Series</td><td>8.3412</td><td>9</td><td>0.9268</td><td>920.66</td><td>2.040</td><td>0.000</td></tr><tr><td>Operator</td><td>0.0002</td><td>2</td><td>0.0001</td><td>0.100</td><td>3.150</td><td>0.905</td></tr><tr><td>Interaction (operator; series)</td><td>0.0212</td><td>18</td><td>0.0012</td><td>1.172</td><td>1.778</td><td>0.301</td></tr><tr><td>Instrument (error)</td><td>0.0604</td><td>60</td><td>0.0010</td><td></td><td></td><td></td></tr><tr><td>Total</td><td>8.4231</td><td>89</td><td>0.0946</td><td></td><td></td><td></td></tr></tbody></table>

$$ \left\{ \begin{aligned} F_{\text{expe,SCM}} &= \frac{\text{MSCM}}{\text{MSCE}} \\ F_{F-S,SCS} &= F_{\text{Fisher-Snedecor}}(\alpha, \text{deg}_{\text{SCS}}, \text{deg}_{\text{SCE}}) \end{aligned} \right. \quad 22 $$

$$ \left\{ \begin{aligned} F_{\text{expe,SCO}} &= \frac{\text{MSCO}}{\text{MSCE}} \\ F_{F-S,SCO} &= F_{\text{Fisher-Snedecor}}(\alpha, \text{deg}_{\text{SCO}}, \text{deg}_{\text{SCE}}) \end{aligned} \right. \quad 23 $$

$$ \left\{ \begin{aligned} F_{\text{expe,SCOM}} &= \frac{\text{MSCOM}}{\text{MSCE}} \\ F_{F-S,SCOS} &= F_{\text{Fisher-Snedecor}}(\alpha, \text{deg}_{\text{SCOS}}, \text{deg}_{\text{SCE}}) \end{aligned} \right. \quad 24 $$

The critical case in comparing $F_{\text{expe}}$ and $F_{F-S}$ arises when $F_{\text{expe}}$ is smaller than $F_{F-S}$. In this case, the variance should be forced to zero, and this adjustment must be reflected in the error term. The test should then be repeated with the new values. To illustrate these two steps, let's assume that $F_{\text{expe}}$ for SCOS is smaller than $F_{F-S}$ for SCOS, and for all other comparisons, $F_{\text{expe}}$ exceeds $F_{F-S}$. The new values obtained in this case are shown in the following table (Table 4):

Repeating the comparison between $F_{F-S}$' and $F_{\text{expe}}'$ does not alter the result (due to the nature of variations in the Fisher-Snedecor distribution), but it allows for the calculation of the associated risk (denoted as $R_{\text{isk}}$). By definition, according to the Fisher distribution (denoted as F), this corresponds to the one-tailed probability to the right of the value $F_{\text{expe}}$ (Equations (25)–(27)). This Fisher distribution is implemented in Excel® as: FDIST(F'_{\text{expe}}; deg'_{\text{num}}; deg'_{\text{deno}}; FALSE).

$$ R_{\text{isk,SCM}} = F(F'_{\text{expe,SCM}}, \text{deg}'_{\text{SCM}}, \text{deg}'_{\text{SCE}}) \quad 25 $$

$$ R_{\text{isk,SCO}} = F(F'_{\text{expe,SCO}}, \text{deg}'_{\text{SCO}}, \text{deg}'_{\text{SCE}}) \quad 26 $$

$$ R_{\text{isk,SCOM}} = F(F_{\text{expe,SCOM}}, \text{deg}_{\text{SCOM}}, \text{deg}_{\text{SCE}}) \quad 27 $$

Table 5 presents the results of various formulas (Equations (12)–(27)) applied to the resistance measurements (Table 1).

For all sources of error, the experimental Fisher test is always lower than the theoretical Fisher test, which implies that the procedure explained earlier (Table 4) should be carried out. The result of this new procedure are shown in Table 6.

Finally, it is necessary to calculate, assuming existence, the square root of the mean squares to determine the uncertainty associated with repeatability $\sigma_r$, reproducibility $\sigma_R$, variations between operators and measurements $\sigma_{\text{OM}}$, variability between measurement series $\sigma_s$, and

<caption>**Table 6**</caption>

New results of the Anova method performed on Fig. 2.

<table><thead><tr><th>Error source</th><th>Sum of squares</th><th>Degrees of freedom</th><th>Average square</th><th>Experimental Fisher test</th></tr></thead><tbody><tr><td>Series</td><td>8.3412</td><td>9</td><td>0.9268</td><td>0</td></tr><tr><td>Operator</td><td>0.0000</td><td>2</td><td>0.0001</td><td>1</td></tr><tr><td>Interaction (operator; series)</td><td>0.0000</td><td>18</td><td>0.0012</td><td>1</td></tr><tr><td>Instrument (error)</td><td>0.0818</td><td>80</td><td>0.0010</td><td></td></tr><tr><td>Total</td><td>8.4231</td><td>89</td><td>0.0946</td><td></td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>Fig. 7. Histogram of the isotopic residuals.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Normal probability plot of the isotopic residuals.</figcaption></figure>

total variability σtot (Equations (28)–(33)). If the quantity calculated from Equations (28)–(33) is undefined, then the uncertainty value is equal to 0. In Excel®, the SQRT function should be used to compute the square root of Equations (28)–(33).

$$ \sigma_r = \sqrt{\text{MSCE}'} \qquad 28 $$

$$ \sigma_R = \sqrt{\frac{\text{MSCO}' - \text{MSCOS}'}{\text{nb}_m \cdot \text{nb}_s}} \qquad 29 $$

$$ \sigma_{\text{OS}} = \sqrt{\frac{\text{MSCOS}' - \text{MSCE}'}{\text{nb}_m}} \qquad 30 $$

$$ \sigma_{\text{R&R}} = \sqrt{\sigma_r^2 + \sigma_R^2 + \sigma_{\text{OS}}^2} \qquad 31 $$

$$ \sigma_s = \sqrt{\frac{\text{MSCS}' - \text{MSCOS}'}{\text{nb}_m \cdot \text{nb}_o}} \qquad 32 $$

$$ \sigma_{\text{tot}} = \sqrt{\sigma_{\text{R&R}}^2 + \sigma_s^2} \qquad 33 $$
<figure><img src="image_11.png" /><figcaption>Fig. 9. Box and whisker plots for An unarranged δ¹³C and for B ranged δ¹³C. In the legend, the acronyms Op and Obj mean respectively Operator and series of the Object measured.</figcaption></figure>

These various uncertainties for the resistance measurements are as follows: 0.032 for σ<sub>r</sub>, 0.000‰ for σ<sub>R</sub>, 0.000‰ for σ<sub>OS</sub>, 0.321‰ for σ<sub>s</sub>, and 0.323‰ for σ<sub>tot</sub>. It is recommended to present the confidence intervals by multiplying the standard deviation by a factor of 5.15 [1] and 6 for [31].

To interpret these uncertainties, it is necessary to calculate the
relative contributions in percentage (denoted as C) [25,30] (Equations
(34)–(38)).

$$
C_r = \frac{\sigma_r}{\sigma_{tot}} \qquad 34
$$

$$
C_R = \frac{\sigma_R}{\sigma_{tot}} \qquad 35
$$

$$
C_{OS} = \frac{\sigma_{OS}}{\sigma_{tot}} \qquad 36
$$

$$
C_{R\&R} = \frac{\sigma_{R\&R}}{\sigma_{tot}} \qquad 37
$$

$$
C_s = \frac{\sigma_s}{\sigma_{tot}} \qquad 38
$$

The contributions for the isotopic measurements are as follows: 9.92 % for C<sub>r</sub>, 0.00 % for C<sub>R</sub>, 9.92 % for C<sub>R&R</sub>, 0.0 % for C<sub>OS</sub>, and 99.51 % for C<sub>s</sub>. Different thresholds have been defined by evaluating these capabilities, particularly the gauge R&R capability C<sub>R&R</sub>. If the C<sub>R&R</sub> capability is: (i) less than 10 %, the process is satisfactory, (ii) between 10 % and 30 %, the process is improvable, (iii) greater than 30 %, the process is degraded [2]. As the value of capability C<sub>R&R</sub> is less than 10 % for the resistance measurement, the process is satisfactory. In other words, the primary source of variability originates predominantly from the measured object itself, rather than, for instance, from the operator.

3.1.5. Anova analysis

To validate an Anova method, the residuals must be calculated, and
the normality of these residuals must be verified [12,47]. The residuals
r<sub>jk</sub> are, by definition, the difference between the measurement y<sub>ijk</sub> and
the average of the measurements according to the number of measure-
ments (Equation (39)). In Excel®, Equation (39) can be implemented
using the AVERAGE function.

$$
r_{jk} = y_{ijk} - \frac{1}{nb_m} \sum_{k=1}^{nb_m} y_{ijk} \qquad 39
$$

<figure><img src="image_12.png" /><figcaption>Fig. 10. Modified GRR diagrams for individual operators.</figcaption></figure>
<table><caption>**Table 7**
Coefficient of direction *a* and y-intercept *b* obtained for the different operators. *t*₁ and *t*₂ values calculated from equation (7) and the data in Table 2.</caption><thead><tr><td>Operator</td><td>t<sub>1</sub></td><td>t<sub>2</sub></td></tr></thead><tbody><tr><td>1</td><td>0.70</td><td>1.79</td></tr><tr><td>2</td><td>1.30</td><td>3.43</td></tr><tr><td>3</td><td>0.21</td><td>0.78</td></tr></tbody></table>

The histogram (Fig. 7) allows for the visualization of these residuals $r_{jk}$, which are obtained in Excel® by following these steps: Select the **Insert** tab, then select **Insert Statistical Chart**, and under **Histogram**, choose a histogram.

The normal probability plot is a complementary graph to the histogram.

It is obtained (Fig. 8) by displaying the probability percentage as a function of the residuals, following these steps.

1. Sort the residuals in ascending order: Copy only the residual values into another column. Select the **Home** tab, then select **Sort & Filter**, and click on **Sort Smallest to Largest**.

2. Assign a rank to each residual, i.e., 1 for the smallest value, 2 for the second smallest value, and so on.

3. Calculate the percentile $p_u$ corresponding to the rank $r_u$ for a measurement *u* (Equation (40)) using the MAX and MIN functions.

4. Generate the plot: Select the **Insert** tab, then select **Insert Scatter (X, Y)** or **Bubble Chart**, and under **Scatter**, pick a chart.

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Range control chart.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Average control chart.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 13. Screenshot of Excel® spreadsheet with formulas associated with each box.</figcaption></figure>

$$p_u = \frac{r_u - \min_u(r_u) + 1}{\max_u(r_u) - \min_u(r_u) + 2}$$

Despite the presence of outlier coordinate values, specifically (0.04, 97.80) and (0.05, 98.90), Figs. 7 and 8 show that the residuals are normally distributed. In particular, the Pearson correlation coefficient of the linear regression performed on the residuals (Fig. 8) is close to one (0.9865). This analysis validates the Anova method by verifying the assumptions of the model.

## 3.2. Case 2: isotopic measurement

### 3.2.1. Graphical tools

This approach is once again applied to the isotopic measurements. Therefore, it is essential to first verify the assumptions of the mathematical model (see the Mathematical Model section). Consequently, the statistical distribution of the measurements must be analyzed through control charts, particularly using Boxplots, which facilitate the detection of outliers and the symmetry of the measurement distributions.

Fig. 9 presents the Box and Whisker plots for each measured object—initially in their original sequence in part A, and subsequently

<figure><img src="image_16.png" /><figcaption>Fig. 14. Histogram of the isotopic residuals.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 15. Normal probability plot of the isotopic residuals.</figcaption></figure>

sorted by standard deviation in part B. This visualization provides in-
sights comparable to those offered by the scatter plot recommended in
Ref. [2]. Consequently, both types of graphs can be displayed and
interpreted together [43].

The analysis of the Box and Whisker plots (Fig. 9) reveals an outlier
in the first series of measurements by Operator 1, with a value of 0.69 ‰
when examining the data from Fig. 2. Outliers can also lead to asym-
metric distributions. For example, this is observed in the first and second
measurement series of Operator 2, as well as in the third series of
Operator 3.

It is also worth noting that Operator 2 produces measurements with less dispersion compared to the other two operators, but with a higher average value. The Modified GRR diagram provides a means to further explore this analysis. Fig. 10 shows this Modified GRR diagram.

The x-axis of Fig. 10 provides information about the uniformity of the
measurement system. The overall means of the measured objects are
grouped into two clusters: one around 1.45 ‰ and the other around
1.56 ‰. This indicates significant uniformity originating from the object
under study. Indeed, the measurements taken during this study are
geological in nature and made on different samples. These measure-
ments are therefore likely to give variations in the order of 0.1 ‰ [48].

Fig. 10 also allows for the evaluation of the variability of individual
measurements from a given operator and their comparison to the overall
mean. By using the line representing the means of the product samples,
it is possible to assess whether some operators exhibit larger deviations
from these means. These deviations represent measurement biases [1].
Operator 1 shows the least variability compared to the mean, while
Operator 2 shows the most. Operator 2 tends to measure values higher
than the mean, as discussed in the Box and Whisker plot section. Addi-
tionally, for these two operators, it seems that this shift is constant, as
the average curve and the curve obtained through linear regression
almost overlap. In other words, if there is a measurement bias, it appears
to be constant.

To quantify these analyses, the manual [2] suggests using statistical
tests t₁ and t₂ (equation (7)) on the linear regression curves, specifically
on the slope coefficient *a* and the y-intercept *b*. Table 7 provides the
values of the quantities t₁ and t₂ for the three operators.

The value of the Student's t-test obtained using the function T.INV is
2.78. According to Table 7, the theoretical value of the Student's t-test is,
for the different operators, higher than the value of the quantity t₁. In
other words, there is no bias dependent on the measured objects.

However, the value of the Student's t-test is lower for the value of the
quantity t₂ for operator 2. This means that this operator makes biased
measurements. In the critical case where an operator seems to be a
possible cause of unexpected measurement results, it is preferable to
perform a new analysis of the repeatability and reproducibility by
excluding the measurements taken by this operator [45].

To pursue the analysis of the persistent bias associated with operator
2, it is necessary to assess the variability of repeated measurements using
the modified range control chart. Fig. 11 represents the range control
chart.

The analysis of Fig. 11 shows that the system under study is under control concerning the variability of repeated measurements because all the ranges are within the LCL (Lower Control Limit) and UCL (Upper Control Limit). Operator 1 exhibits the highest variability due to the outlier value. Operator 2 has an average range three times smaller than that of Operator 1 due to a low dispersion of the measurements. For the range analysis to be meaningful, it must be conducted in conjunction with the mean value [49], as illustrated Fig. 12.

As with the range control chart, all the means $\bar{y}_{ij}$ fall within the control limits. This case is exceptional because few measurements are made. In general, not all values fall within the control limits. Provided that the measured object covers the entire variability of the process, a measurement system is considered acceptable for evaluating the variability between the measured samples if at least half of the plotted means lie outside the control limits [2].

Ultimately, the various analyses carried out using control charts have
revealed that the measurements are biased, thereby invalidating the
assumptions of the Anova model (Mathematical Model section). At this
stage, it is therefore necessary to eliminate the outlier measurements
associated with Operator 2 and to initiate new measurement series,
which should be conducted in greater numbers.

3.2.2. Anova method

While the control charts have highlighted significant problems in the
measurement data, the Anova method is still implemented, with the sole
purpose of serving as a pedagogical exercise. Fig. 13 is a screenshot
summarizing the Anova method implemented in an Excel® spreadsheet.

These various uncertainties for the isotopic measurements are as
follows: 0.25 ‰ for σ<sub>r</sub>, 0.00 ‰ for σ<sub>R</sub>, 0.00 ‰ for σ<sub>OS</sub>, 0.00 ‰ for σ<sub>s</sub> , and
0.35 ‰ for σ<sub>tot</sub>. So, the contributions associated for the isotopic mea-
surements are as follows: 70.7 % for C<sub>r</sub>, 0.0 % for C<sub>R</sub>, 70.7 % for C<sub>R&R</sub>,
<figure><img src="image_18.png" /><figcaption>Fig. 16. Application user interface for the range method and bootstrap method.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 17. Application user interface illustrating imported data.</figcaption></figure>

0.0 % for $C_{OS}$, and 0.0 % for $C_s$. As the $C_{R&R}$ value is greater than 30 %, the process of isotopic measurements is degraded, which can be explained by the analysis of the graphical tools performed earlier. To validate an Anova method, the residuals must be calculated, and the normality of these residuals must be analyzed (Fig. 14) and verified (Fig. 15).

Figs. 14 and 15 show that the residuals are not normally distributed.

This analysis does not validate the Anova method by verifying the assumptions of the model and identifying excessively influential measurements. The outliers, for instance at the coordinates (– 0.59, 3.57), detected in the Graphical tools section, are highly influential.
<figure><img src="image_20.png" /><figcaption>**Fig. 18.** Results obtained using the bootstrap method on the measurements shown in Fig. 2.</figcaption></figure>

## 4. Discussion

The Excel® spreadsheet makes it easy to perform a Gauge R&R study. In particular, once the user has implemented a Gauge R&R study,

another study can be quickly and easily conducted by changing the values in the spreadsheet cells. However, the conditions for the Gauge R&R must remain consistent (same nb<sub>m</sub>, nb<sub>o</sub>, and nb<sub>s</sub>). The challenge arises in the implementation of the Gauge R&R, particularly when the number of measurements (nb<sub>m</sub>·nb<sub>o</sub>·nb<sub>s</sub>) is large. Indeed, as this number increases, the calculation steps quickly become tedious and time-consuming. In such cases, it is necessary to use statistical software products such as JMP or Minitab.

Another limitation of using the Excel® spreadsheet is performing a Gauge R&R on measurements that are not normally distributed, as highlighted by the issue of isotopic measurements. Indeed, if the measurement values are not normally distributed, while the Excel® spreadsheet can identify outliers, it does not allow for the execution of a non-parametric Anova method. Osma's study [12] suggests using the Range method as an alternative. Another possibility is to perform permutations on the measurements [50] using the Bootstrap method [51, 52]. The statistical test can also be modified to correspond to a non-parametric test, such as the Jonckheere Terpstra test [53]. Additionally, transformations of the measurements can be applied, such as aligned rank transformation [54] or transformations based on functions like the logarithm or the hyperbolic sine [55]. In cases where non-parametric methods are restrictive, a non-metric Anova can be performed to test the disparity of measurement groups [56]. Therefore, while several methods exist for performing Anova, they are very difficult to implement in Excel®. To overcome this limitation, a software developed in Python, available at the following link: https://github.com/ICBAlban/Can-Excel-be-used-for-Gauge-R-R-study-based-on-th-e-Analysis-of-Variance-Anova-.git, has been created to implement the Range and Bootstrap methods from an Excel® spreadsheet (Fig. 16). The Range method implemented is based on the code provided in the study

<figure><img src="image_21.png" /><figcaption>**Fig. 19.** Results obtained by performing the range method on the measurements in Fig. 2.</figcaption></figure>
[20]. The user interface of this application consists of three frames. The first frame, named 'Parameter', corresponds to the parameters nb<sub>m</sub>, nb<sub>o</sub>, and nb<sub>s</sub>, which must be correctly filled in to import the measurements from an Excel® spreadsheet. Measurement import is carried out by pressing the 'Open' button. The format of the imported Excel® file must match that presented in Table 1 or Fig. 2. If the format does not match, the measurement values cannot be imported, and an error message will appear.

Once the Excel® spreadsheet is opened (see Fig. 17), the measure-
ments are displayed in the second frame labeled 'Data'. The displayed
measurements can be directly modified within this application. The user
must then choose between the Bootstrap method and the Range method.
For the Range method, the calculation of the d<sub>2</sub> and d<sub>3</sub> functions is based
on a random sampling approach [20], the accuracy of which depends on
the number of iterations. The criterion used in this software [20] is that
the power of the number of iterations divided by two corresponds to the
precision of the calculated uncertainties. In this example, the number of
significant figures in the measurements is 2. Therefore, a minimum of
10<sup>4</sup> iterations should be performed to achieve an uncertainty of 1/10<sup>4</sup>/2 =
0.01, or 10<sup>6</sup> iterations for an uncertainty finer than the measurements, i.
e., below 0.01. The number of iterations is crucial to obtaining results
with appropriate precision [20]. Additionally, the values for the control
limits (USL and LSL) must be adjusted in order to perform an honest
gauge R&R [57]. Fig. 18 presents the results obtained using the Boot-
strap method, while the results of the Range method are presented in
Fig. 19.

According to Figs. 18 and 19, the uncertainty associated with
reproducibility σ<sub>R</sub> and the measurement series σ<sub>s</sub>, or σ<sub>p</sub> as noted in
Fig. 19, are non-zero. These non-parametric methods, therefore, allow
for the determination of uncertainty arising from factors other than
repeatability. Given that the uncertainty associated with repeatability σ<sub>r</sub>
is dominant, the isotopic measurement process is impaired. These
methods yield results consistent with the Anova approach. Conse-
quently, the entire set of measurements is fundamentally flawed. To
improve the measurement process, it is necessary to (i) repeat the
measurements by increasing their number to meet the recommendations
outlined in the Measurement section, for instance, by using 3 operators,
10 measurement series, and 3 measurements per series, or (ii) imple-
ment a continuous improvement approach as discussed in the intro-
duction to identify error sources and enhance measurement precision.

**5. Conclusion**

Excel® is an effective tool for conducting a Gauge R&R analysis using the Anova method. In particular, Excel® can be used by an operator, technician, or engineer, as this software does not require advanced programming expertise. The cells in this spreadsheet can also be frozen to prevent errors and, therefore, the loss of information. This method relies on a mathematical model based on several assumptions, as discussed in the Mathematical Model section. Therefore, before performing the Anova method, it is essential to verify these assumptions. Excel® facilitates this verification through various graphical representations. If these assumptions are not met, the root cause must be identified, which may, as mentioned in the Graphical Tools section, be related to the operator, the measurement system, or the measured object. Once the cause is identified, new measurements should be taken, or outlier data points should be filtered out.

Once the assumptions underlying the Gauge R&R are validated, the Anova method can be implemented. With the use of various built-in functions such as SUM, AVERAGE, SQRT, and POWER, along with statistical tests like FDIST, the Anova calculations can be performed within the software.

To interpret the results obtained through the Anova method, re-
siduals can be quickly computed and visualized in Excel® in the form of
histograms or Normal Probability Plots. If the residuals do not follow a

normal distribution, the reliability of the results is questionable. In such
cases, it is advisable to filter out the outliers identified using graphical
tools and to rerun the Anova method. If, even after filtering, the re-
siduals still do not follow a normal distribution, it is recommended to
redo the Gauge R&R analysis using non-parametric Anova method such
the range method or the bootstrap method. A software that implements
these methods in an Excel® spreadsheet is available at the following
address: https://github.com/ICBAlban/Can-Excel-be-used-for-Gauge-R-
R-study-based-on-the-Analysis-of-Variance-Anova-.git.

Declaration of competing interest

The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.

**Appendix A. Supplementary data**

Supplementary data to this article can be found online at https://doi.
org/10.1016/j.measen.2025.101978.

Data availability

Data will be made available on request.

**References**

[1] BIPM, Vocabulaire international de métrologie - Concepts fondamentaux et généraux et termes associés (VIM), 3er ed., 2008.

[2] Automotive Industry Action Group, Measurement Systems Analysis, 4rd ed., Southfield, 2010.

[3] J. Antony, G. Knowles, et al.P. Roberts, Gauge capability analysis: classical versus ANOVA, Qual. Assur. 6 (n° %13) (1999) 173-181.

[4] S.B. Vardeman, et al.J.M. Jobe, Basic Engineering Data Collection and Analysis, Iowa: Iowa State University, 2001.

[5] M.-H.C. Li, et al.A. Al-Refaie, Improving wooden parts' quality by adopting DMAIC procedure, Qual. Reliab. Eng. Int. 24 (2008) 251-360, n° %13.

[6] M. Sharma, S.P. Sahni, and S. Sharma, Validating a destructive measurement system using Gauge R&R—a case study, Eng. Manage. Prod. Services. 11 (2019) 34–42, n° %14.

[7] F.W. Breyfogle, et al.B. Meadows, The six sigma implementation process, SAE Tech. Paper. (2001).

[8] J.M. Juran, et al.A.B. Godfrey, The Quality Control Process, McGraw-Hill, New York: NY, USA, 1999.

[9] A. Sarkar, A.R. Mukhopadhyay, et al.S.K. Ghosh, Measurement system analysis for implementing design for Six Sigma, Int. J. Prod. Qual. Manag. 14 (2014) 373–386, n° %13.

[10] M. Persijn, et al.Y.V. Nuland, Relation between measurement system capability and process capability, Qual. Eng. 9 (1996) 95-98, n° %11.

[11] D.W. Hoffa, et al.C. Laux, Gauge R&R: an effective methodology for determining the adequacy of a new measurement system for micron-level metrology, J. Industrian Technol. 23 (2007) 1-9, n° %14.

[12] A. Osma, An assessment of the robustness of gauge repeatability and reproducibility analysis in automotive components, Proc. Inst. Mech. Eng. - Part D J. Automob. Eng. 225 (2011) 895-912, n° %17.

[13] D.C. Montgomery, et al.G.C. Rugner, Gauge capability and designed experiments. Part I: basic methods, Qual. Eng. 6 (1993) 115-135, n° %11.

[14] D.C. Montgomery, et al.G.C. Runger, Gauge capability analysis and designed experiments. Part II: experimental design models and variance component estimation, Qual. Eng. 6 (1993) 289-305, n° %12.

[15] R.K. Burdick, C.M. Borror, et al.D.C. Montgomery, A review of methods for measurement systems capability analysis, J. Qual. Technol. 35 (2003) 342–354, n° %14.

[16] G.P. Quinn, et al.M.J. Keough, Experimental Design and Data Analysis for Biologists, Cambridge university press, Cambridge, 2002.

[17] J.W. Joubert, et al.S. Meintjes, Repeatability & reproducibility: implications of using GPS data for freight activity chains, Transp. Res. Part B Methodol. 76 (2015) 81–92.

[18] C.S. Ackermann, «Evaluating destructive measurements using gage R&R, in: Proceedings. IEEE/SEMI Advanced Semiconductor Manufacturing Conference and Workshop, 1993, pp. 101–105.

[19] W. Zhan, Gauge R&R and Troubleshooting, ASEE Annual Conference & Exposition, 2016, 2016.

[20] A. Petitjean, et al.O. Musset, Application for determining repeatability and reproducibility in an honest gauge R&R study, Qual. Reliab. Eng. Int. (2025).

[21] M. Koretsky, I. Ntroducing ChE Sophomores to Measurement System Analysis and Analysis of Variance through Experiential Learning, chez 2005 Annual Conference, Portland, 2005.
[22] A. Waseem, A. Zulqarnain, R. Khalid, et al.S. Saleem, Gauge repeatability and reproducibility: an assessment of method and medium for a laboratory, Int. J. Six Sigma Compet. Advant. 9 (2015) 126-149, n° %12-4.

[23] J. Mandel, «Repeatability and reproducibility 4 (1972) 174-185.

[24] J.-N. Pan, «Determination of the optimal allocation of parameters for gauge repeatability and reproducibility study 21 (6) (2004) n° %1672-n° %1682.

[25] J.-N. Pan, Evaluating the Gauge Repeatability and Reproducibility for Different Industries, 40, Springer, 2006, pp. 499-518.

[26] S.B. Vardeman, et al.J.M. Jobe, Statistical Methods for Quality Assurance, 2e éd., Springer-Verlag, New York, 2016.

[27] P.-S. d. Laplace, Mémoire sur divers points d'analyse, Journal de l'Ecole Polytechnique, Cahier 15 (1809) 229-265, n° %18.

[28] P.-S. Laplace, Supplément au mémoire sur les approximations des formules qui sont fonction de très grands nombres, Mémoires de l'Académie Royale des Sciences de Paris (1810) 559-565.

[29] J. Bernoulli, Ars coniectandi, Impensis Thurnisiorum (2025).

[30] D.W. Hoffa, et al.C. Laux, Gauge R&R: an effective methodology for determining the adequacy of a new measurement system for micron-level metrology, J. Industrian Technol. 23 (4) (2007).

[31] R.K. Burdick, Y.-J. Park, D.C. Montgomery, et al.C.M. Borror, Confidence intervals for misclassification rates in a gauge R&R study, J. Qual. Technol. 37 (4) (2018) 294-303.

[32] A.A. Saupi, et al.H. Midi, «Outlier detection method in crossed gage repeatability and reproducibility (R&R) random effect model, Malaysian. J. Mathe. Sci. 15 (2021) 333-345, n° %13.

[33] J. De Mast, et al., Trip, «Gauge R&R studies for destructive measurements», J. Qual. Technol. 37 (2005) 40-49, n° %11.

[34] H. Ya-Juan, et H. Zhen, An applied study of destructive measurement system analysis. 2007 2nd IEEE Conference on Industrial Electronics and Applications, IEEE, 2007, pp. 30-34.

[35] G. Box, W.H. Hunter, et al.S. Hunter, Statistics for Experimenters, John Wiley and sons, New York, 1978.

[36] M.S. Hamada, «Alternative designs and analyses for destructive and non-replicable gauge R&R studies, Qual. Eng. 33 (2001) 687-694, n° %14.

[37] A. Petitjean, C. Thomazo, O. Musset, I. Jovovic, P. Sansjofre, et al.K. Kirsimäe, A laser-laser method for carbonate C and O isotope measurement, metrology assessment, and stratigraphic applications, Frontiers. Geochem. I (2024) 1334490.

[38] A. Petitjean, C. Thomazo, O. Musset, I. Jovovic, P. Sansjofre, et al.K. Kirsimäe, Corrigendum: a laser-laser method for carbonate C and O isotope measurement, metrology assessment, and stratigraphic applications, Frontiers. Geochem. 2 (2024) 1447889.

[39] B. Fry, Stable Isotope Ecology, 521, Springer, New York, 2006.

[40] J.M. McCrea, On the isotopic chemistry of carbonates and a paleotemperature scale, J. Chem. Phys. 18 (1950) 849-857, n° %16.

[41] C. Thomazo, P. Sansjofre, Musset Olivier, T. Cocquerez, et al.S. Lalonde, «In situ carbon and oxygen isotopes measurements in carbonates by fiber coupled laser diode-induced calcination: a step towards field isotopic characterization, Chem. Geol. 578 (2021) 120323.

[42] A. Petitjean, O. Musset, P. Sansjofre, C. Thomazo, A laser-laser method for In Situ C and O isotope measurements in carbonate materials, Geostand. Geoanal. Res. (2025) ggr.70023, https://doi.org/10.1111/ggr.70023. (Accessed 20 October 2025).

[43] J. Plura, D. Vykydal, F. Tøsenovský, et al.P. Klaput, Graphical tools for increasing the effectiveness of gage repeatability and reproducibility analysis, Processes. 11 (2022) 1, n° %11.

[44] R.D. Snee, Interpreting operator-part interaction in gage repeatability and reproducibility studies, Qual. Eng. 33 (2021) 395-403, n° %13.

[45] P. Mikulová, et al.J. Plura, Detailed analysis of GRR study results and their visualization, Quality Prod. Improvement-QPI 1 (2019).

[46] ISO 7870-2, Control Charts—Part 2: Shewhart Control Charts, International Organization for Standardization ISO, Geneva, 2013.

[47] A. Osma, A statistical approach to comparing wheel hub forging processes, Proc. Inst. Mech. Eng. - Part D J. Automob. Eng. 223 (2009) 1559-1576, n° %112.

[48] C. Malfilatre, P. Boulvais, M.-P. Dabard, S. Bourquin, E. Hallot, D. Pallix, et al. D. Gapais, Petrographical and geochemical characterization of Comblanchien limestone (Bourgogne, France): a fingerprint of the building stone provenance, C. R. Geosci. 344 (2012) 14-24, n° %11.

[49] D.J. Wheeler, A Better Way to Do R&R Studies, Quality Digest, 2013.

[50] M.J. Anderson, Permutation tests for univariate or multivariate analysis of variance and regression, Can. J. Fish. Aquat. Sci. 58 (2001) 626-639, n° %13.

[51] B. Zhou, et al.W.H. Wong, A bootstrap-based non-parametric ANOVA method with applications to factorial microarray data, Stat. Sin. (2011) 495-514.

[52] C. Giovanni, A bootstrap method for the measurement error estimation in Gauge R & R R\&R Studies, Qual. Reliab. Eng. Int. 38 (2022) 3404-3421, n° %17.

[53] A. Ali, A. Rasheed, A.A. Siddiqui, M. Naseer, S. Wasim, et al.W. Akhtar, Non-parametric test for ordered medians: the Jonckheere Terpstra test, Int. J. Stat. Med. Res. 4 (2015) 203, n° %12.

[54] J.O. Wobbrock, L. Findlater, D. Gergle, et al.J.J. Higgins, «The aligned rank transform for nonparametric factorial analyses using only anova procedures, Proceedings of the SIGCHI conference on human factors in computing systems (2011) 143-146.

[55] G. Genta, et al.M. Galetto, «Study of measurement process capability with non-normal data distributions, Proced. CIRP. 75 (2018) 385-390.

[56] A. Malyutina, J. Tang, et al.A. Amirousefi, «Resolving network clusters disparity based on dissimilarity measurements with nonmetric analysis of variance, iScience. 26 (2023) n° %111.

[57] D.J. Wheeler, An Honest Gauge R&R Study, 2009. Manuscript, n° %1189.