Spatial Risk Assessment: A Case of Multivariate Linear Regression

1. Introduction and Basic Information

In today’s manufacturing environment, where every product must comply with established standards before reaching consumers, the assessment of product quality plays a critical role. Evaluating product conformity involves measuring one or more characteristic properties relevant to the product in question [1]. Under a simple acceptance decision rule, the measured value is compared against the permissible limits defined by the tolerance interval ($T_{\mathrm{I}}$). It is checked whether the measured value of the quantity of interest is within or outside the limits of the tolerance interval [2]. Decision-making becomes particularly challenging when the measured value lies near the boundaries of the tolerance interval, as this can lead to significant economic repercussions for both producers and consumers [3,4,5]. The risk and associated costs arise from mistakenly rejecting a compliant product whose measured value falls within the tolerance limits. Conversely, a product may be accepted as conforming to specifications even though it is not. These potential mistakes are quantified as specific risks for each characteristic property [6,7,8]. When multiple characteristic properties are considered, an overall total specific risk is calculated [9,10].

The measurement result and measurement uncertainty are inextricably linked. Measurement uncertainty refers to the range of values that can be assigned to a measurand [11]. The measured value, along with its associated uncertainty, is used to establish a guard band around the specification limits, with acceptance and rejection zones [12,13,14]. Incorporating a guard band into the conformity assessment process, the risk of making wrong decisions can be effectively controlled, assessed and reduced [15,16,17,18]. In this context, in addition to the tolerance interval, risk assessment also includes an acceptance interval ($A_{\mathrm{I}}$), which specifies the permissible upper and lower bounds for measured values [19,20]. In this study, closed tolerance and acceptance intervals were observed. Within the producer’s risk minimization framework, the tolerance interval is enclosed by the acceptance interval. Conversely, under the consumer’s risk minimization approach, the acceptance interval is contained within the tolerance interval [14,21,22]. In both models, the tolerance and the acceptance interval are separated by a guard band of width w. In the so-called shared risk model, the limits of the tolerance interval and the acceptance interval coincide.

By incorporating the measurement uncertainty of a future measurement process, a probabilistic model to assess risk is obtained based on Bayes' theorem, with the application of decision rules with a guard band [23,24]. The core concept of "point-based" risk assessment was adopted by the Joint Committee for Guides in Metrology (JCGM) and formalized in guide 106:2012 [19]. This procedure yields global producers' and consumers' risks for a single value within the tolerance interval and acceptance interval. Global producer’s and consumer's risks are obtained by combining information about possible, assumed values that the quantity being measured can take with data obtained from future measurements. The true value of the product characteristic under measurement is treated as an unknown parameter associated with a prior probability distribution. Selection of this prior distribution follows the principle of maximum entropy (PME) and is guided by the nature of the measurand and available data sources [25,26,27,28,29,30]. Prior knowledge can be informed by historical measurement data, relevant standards, technical manuals, or expert judgement based on the measurer’s experience [31,32,33]. Essential inputs to define the prior include the best estimate of the measured quantity, denoted $y_0$, and the corresponding measurement uncertainty, $u_0$. Parameters of the prior distribution are estimated using techniques such as the method of moments or maximum likelihood estimation, leveraging information from $y_0$ and $u_0$ [19,34,35]. Depending on the number of available data (both $y_0$ and $u_0$, only $y_0$, or none), the prior can be represented as a two-parameter, one-parameter, or non-parametric distribution [25,35,36]. Measurement values obtained in the future are modeled via a likelihood function, which requires as input the measurement uncertainty, $u_m$ associated with the anticipated measurement process.

The point-based approach for risk estimation, as delineated in Guide 106:2012, yields a single, definitive value for both the global producer's risk and the global consumer's risk [19]. This assumes determining the risk for exactly one value of the best estimate $y_0$. It is used for risk quantification and quality assessment of various products, such as food, water, pharmaceuticals, medical treatments, natural gas, fuels, and other industrial goods [37,38,39,40,41,42,43]. Quality assessment can also be performed for multiple characteristics of an item of interest [44,45,46]. In the case of multi-component quality control, the total risk of false decisions is calculated. This requires determining a guard band for each controlled component [47,48].

By discretizing a naturally selected domain and transforming the argument of the prior distribution in Bayes' theorem, this fundamental point-based method for risk assessment can be two-dimensionally and three-dimensionally extended. Assume the best estimate $y_0$ is fixed and placed within both the tolerance interval and the acceptance interval. It is not required that the position of the best estimate $y_0$ be centered within either the tolerance or acceptance intervals [49]. The ${\mathrm{A}}_{{\mathrm{L}}_{-}}$, and ${\mathrm{A}}_{{\mathrm{U}}_{+}}$denote the lower and upper limits of the acceptance interval of the global producer's risk minimization model, respectively. The notations ${\mathrm{A}}_{{\mathrm{L}}_{+}}$ and ${\mathrm{A}}_{{\mathrm{U}}_{-}}$ refer to the lower and upper limits of the acceptance interval of the global consumer's risk minimization model. Holds that $\left[{\mathrm{A}}_{{\mathrm{L}}_{+}},{\mathrm{A}}_{{\mathrm{U}}_{-}}\right]\mathrm{⸦}\left[{\mathrm{T}}_{\mathrm{L}},{\mathrm{T}}_{\mathrm{U}}\right]\mathrm{⸦}\left[{\mathrm{A}}_{{\mathrm{L}}_{-}},{\mathrm{A}}_{{\mathrm{U}}_{+}}\right]$, Figure 1.

Figure 1 shows an interval representation of the risk assessment domain along the guard band axis. This interpretation allows the construction of risk curves on the interval $\left[-w,w\right]$. Holds that $w={\mathrm{A}}_{{\mathrm{L}}_{+}}-T_L=T_L-{\mathrm{A}}_{{\mathrm{L}}_{-}}={\mathrm{A}}_{{\mathrm{U}}_{+}}-T_U=T_U-{\mathrm{A}}_{{\mathrm{U}}_{-}}$. The two-dimensional extension of the point-based method for risk assessment given in JCGM 106:2012 is obtained by discretizing the interval $\left[-w,w\right]$ [35,43,50]. Values of global consumer's and producer's risk are calculated in each point of subdivision. The limits of the acceptance interval change with the change in the width of the guard band $w$ so that it simultaneously holds that ${\mathrm{A}}_{{\mathrm{L}}_{-}}\to$ ${\mathrm{A}}_{{\mathrm{L}}_{+}}$ and ${\mathrm{A}}_{{\mathrm{U}}_{+}}\to {\mathrm{A}}_{{\mathrm{U}}_{-}}$. In this process, the values ${\mathrm{T}}_{\mathrm{L}}$ and ${\mathrm{T}}_{\mathrm{U}}$ for the lower and upper limit of the tolerance interval are fixed. It is shown that the curves of the global producer's risk along the guard band axis on the interval $\left[-w,w\right]$ increase, while the curves of the global consumer's risk at the same time decrease [43,50,51,52].

In addition to the previously discussed extension, the point-based method for risk assessment can be two-dimensionally extended in a manner that the arguments of prior distribution are functionally related data. This approach is applicable for regression and calibration, allowing the construction of risk curves along the moderate scale [52,53,54]. Risk assessment in regression is carried out naturally, considering the line $y=x$. When the regression line intersects $y=x$ line, the resulting risk curves on the moderated scale take the form of upward-opening parabolas, with their minima located at the intersection points [52]. If the regression line coincides with or runs parallel to $y=x$, the risk curves become linear and run parallel to the axis of the moderated scale [55]. Calibration risk assessment models are derived from these regression-based risk frameworks [53]. These models require defining tolerance and acceptance intervals for the explanatory variable corresponding to the given values of the response variable [53,54]. Within the domain enclosed by the interval on the guard band axis and with an interval on the moderate scale, it is possible to observe the surfaces for global consumers' and global producers' risk in both regression and calibration settings [52,53,54].

In this study, a step further has been taken. This study presents an advancement in risk assessment methodology by introducing a three-dimensional extension of the point-based approach. The method enables spatial risk assessment and is suitable for application in the case of multiple linear regression, resulting in the construction of a 3D risk space. The extension is demonstrated through a practical example involving the risk assessment of a plane that was determined in estimating the flatness error of a coordinate-measuring machine (CMM) worktable. This control is important because deviations in the flatness of the CMM worktable from the limits specified by the manufacturer can significantly affect the measurement results in procedures such as calibration of standards, measuring equipment, etc. [56]. The equation of the plane through the sampled points was determined using the minimum zone (MZ) method. Several variants of the MZ method are prevalent in practice [57]. In this study the One Point Plane Bundle Method (OPPBM) was used, which identifies an optimal plane through the sampled points, minimizing the distance between points with the most extreme values measured [58]. The method is recommended by the ISO 1101:2017 standard for verifying form errors within a specified zone [59]. The risk space was defined by the coordinates of sampled points of the CMM worktable, whose flatness is being assessed. The third axis of the risk space is represented by the values of the global risk of producers and consumers computed at subdivision points of the interval $\left[-w,w\right]$, as illustrated in Figure 1.

The spatial risk assessment model was evaluated using confusion matrix-based metrics. There are four possible outcomes in conformity assessment of a product with the required specifications, Figure 2.

The true value (T_V) of a characteristic property of an item of interest may fall within the tolerance interval, while the measured value (M_V) may fall within the acceptance interval. This is related to the valid acceptance of measurements, i.e., the probability of true positive values (TP), in the conformity assessment process. If the M_V falls outside the acceptance interval, assuming that the T_V is outside the tolerance interval, a valid rejection of the measurement occurs, resulting in a true negative probability (TN). If the M_V of a property of an item of interest is within the acceptance interval, but the T_V of that property is assumed to fall outside the tolerance interval, this results in a global consumer's risk (Rc) or false acceptance probability (FA). A scenario where M_V is outside the acceptance interval while the T_V is within the tolerance interval represents the global producer's risk (Rp) or the probability of false rejection (FR) [19,35,51].

2. Materials and Methods

2.1. Measurement Description

The granite surface of a CMM worktable can be mechanically damaged during exploitation. This can significantly affect the accuracy of measurements taken with the given device [56,60,61]. The risk of measurements on the CMM worktable concerning the flatness of its granite surface was assessed in this study. Measurements were conducted at the Department of Production, Faculty of Technical Sciences, University of Novi Sad, by using a CMM device whose worktable surface is the object of interest. A Carl Zeiss Contura G2 CMM with an RDS rotating head and a VAST XXT passive measuring sensor was used. With this CMM, measurements can be performed in the range of 1000×1200×600 mm. According to the manufacturer's specifications, the maximum permissible error is ${\mathrm{M}\mathrm{P}\mathrm{E}}_{\mathrm{E}}=1.9+L/330$ µm. The value for L is given in mm. The measurement was performed with a silicon nitride (SiN) probe with a single tip, 5 mm in diameter and 75 mm in body width [56]. Sample collection was carried out at $n=40$ measurement sites, and measurements were collected point by point. Three measurements were performed at each of the sample points. The coordinates of the sampled points $\left(x_i,y_i\right),i=1,\dots ,n$ of the granite table of the CMM and the mean value calculated for three measurements for the flatness of the surface along the z-axis are given in Table 1 in [56]. The coordinate values of the sampled points and the values for the surface flatness in that table are given in millimeters.

2.2. Best Estimate and Measurement Uncertainty

The equation of the $\alpha$ plane of the granite table of the CMM was determined based on the measurements given in Table 1 in [56]. The software solution for determining the plane using the MZ method was developed at the Faculty of Technical Sciences in Novi Sad [58]. The equation is of the form $z=Ax+By+D=0$ and is given by the expression:

$$\alpha \dots z=1.71·{10}^{-6}x+1.715·{10}^{-6}-506.29043.$$

(1)

The parameters A, B, and D in equation (1) are given in millimeters. One of the input parameters of the model for estimating the global producer's risk $R_P$ and the global consumer's risk $R_C$ is the best estimate $y_0$ of an item of interest. In this case, the best estimate is represented as an n-tuple $y_0=\left(z\left(x_1,y_1\right),z\left(x_2,y_2\right),\dots ,z\left(x_n,y_n\right)\right),n=40$ whose components are obtained by inserting the coordinates of the sampled points $\left(x_i,y_i\right),i=1,\dots ,n$ into equation (1). For simplicity, the notations $z_i=z\left(x_i,y_i\right),i=1,\dots ,n$ for the components of the best estimate $y_0$ were introduced. The global producer's risk and global consumer's risk are computed for each value $z_i,i=\mathrm{1,2},\dots ,n.$ For the calculation of global producer's risk and global consumer's risk, the point-based method for risk assessment is still used, but the functional relationship of the data given by equation (1) is considered. The measurement uncertainty $u_0$ associated with the best estimates $z_i,i=\mathrm{1,2},\dots ,n$ is the same for all sampled points $\left(x_i,y_i\right),i=1,\dots ,n$ and is ${\mathrm{u}}_0=1.9\mathsf{\mu }\mathrm{m}$. This measurement uncertainty $u_0$ is determined according to research [62,63].

Each observed point $\left(x_i,y_i\right),i=1,\dots ,n$, near the regression plane $\alpha$ is assumed to be a realization from a normal distribution centered at the predicted value of the regression plane. Therefore, the prior has the form $g_0\left(\eta \right)=\phi \left(\eta ;z_i,{u_0}^2\right),i=\mathrm{1,2},\dots ,n$.

The measurement uncertainty $u_m$ required for the definition of the likelihood function is determined in the future measurement process, independently of the current procedure used to derive the plane from equation (1). The likelihood function $h\left({\eta }_m|\eta \right)$ of a measured quantity value ${\eta }_m$, with a given value $\eta$ is also modelled as a normal distribution and holds that $h\left({\eta }_m|\eta \right)=\phi \left({\eta }_m;\eta ,{u_m}^2\right)$ [19]. It is presumed that the measurement uncertainty $u_m$ of a future measurement process controlling the quality of the CMM granite worktable, after exposure to exploitation, has a value $u_m=u_0/2$. In addition to the stated values, the study also provides a comparison of risks for the scenarios where $u_m=u_0$ and $u_m=2u_0.$

2.3. Selection of the Reference Plane

The assessment of global producer's and consumer's risk for the linear regression model is conducted in relation to the $y=x$ line [52]. The estimate provides the risk of the regression line deviating from the $y=x$ line. The line $y=x$ naturally imposes itself as a reference line for risk assessment in linear regression models because it represents an ideal measurement. In this study, for spatial risk assessment purposes, it is assumed that data are functionally connected by a multiple linear regression model. However, in this instance, there is no inherent reference plane to serve as a benchmark for evaluating risk. Therefore, that plane must be defined. Since the risk is determined in the process of assessing the flatness of a surface, the reference plane is assumed to be parallel to the xy-plane. Considering that the MZ method was employed to derive the equation for the flatness testing of a granite table, the plane taken as the reference plane is as follows:

$$\beta \dots z_r={\displaystyle \frac{z_{min}+z_{max}}{2}}=-506.29085\mathrm{m}\mathrm{m},$$

(2)

where is $z_{min}=-506.2926$ mm, $z_{max}=-506.2891$ mm, and $z_r$ is the value of a reference plane. The choice of plane $\beta$ from equation (2) is based on the process of determining the equation of the plane using the MZ method. The reference plane $\beta$ is established through the center point of the parallelogram defined by the maximum and minimum value of measured surface flatness reported by a CMM [58].

2.4. Tolerance Plane and Acceptance Plane

In spatial risk assessment, the role of tolerance intervals and acceptance intervals is taken over by tolerance planes and acceptance planes. The tolerance space is bounded by the lower tolerance plane ${\tau }_L$ and the upper tolerance plane ${\tau }_U$. These two planes are parallel to the plane $\beta$ as well as to one another. Since it is a class 0 granite surface plate, the width of the tolerance space is equal to $∆T=7.5$ μm [64]. For each value $z_i,i=\mathrm{1,2},\dots ,n$ the upper tolerance plane has the equation ${{\tau }_U\dots T}_U=-506.2871$ mm, and the lower tolerance plane has the equation ${{\tau }_L\dots T}_L=-506.2946$ mm. The width of the acceptance space in the model of minimization of global producer's risk is equal to $∆A=1.2∆T$. In the model of minimization of the global consumer's risk, the width of the acceptance space is $∆A=0.8∆T$. The discretization of the domain $I_w=\left[-w,w\right]$ along the guard band axis is carried out by introducing the multiplicative factor $r_k\in \left[-\mathrm{1,1}\right]$ of the form $r_k=-1+0.1\left(k-1\right),k=\mathrm{1,2},\dots ,m.$ An interval subdivision method was performed with a step of 0.1. This leads to $m=21$ subdivision nodes. A finer step would increase the number of nodes. However, it is shown that a subdivision step of 0.1 is sufficient for determining the global producer's and consumer's risk, and the obtained risk surfaces are sufficiently smooth. The subdivision nodes of the interval $I_w$ are of the following form:

$$w_k={\displaystyle \frac{r_k{w}_U}{2}},k=\mathrm{1,2},\dots ,m,$$

(3)

where $w_U=0.2∆T=1.5$ µm is the total width of the guard band. The acceptance planes for fixed values of $T_U$ and $T_L$ are calculated from the following equations:

$$A_{U_{i,k}}=T_{U_i}-w_k,i=\mathrm{1,2},\dots ,n,k=\mathrm{1,2},\dots ,m,$$

(4)

$$A_{L_{i,k}}=T_{L_i}+w_k,i=\mathrm{1,2},\dots ,n,k=\mathrm{1,2},\dots ,m.$$

(5)

The index $i$ in equations (4-5) for the values of the lower and upper limits of the acceptance interval indicates that the interval refers to the points $z_i,i=\mathrm{1,2},\dots ,n$, and it holds that $T_{L_i}=T_L$ and $T_{U_i}=T_U$ for each $i=\mathrm{1,2},..,n$.

For $k=1$ it is true that $r_1=-1$, i.e., $A_{L_{i,1}}={\mathrm{A}}_{{\mathrm{L}}_{-}}$ and $A_{U_{i,1}}={\mathrm{A}}_{U_{+}}$. These are the conditions of the model of minimization of the global producer's risk. For $k=11$ it is true that $r_{11}=0$, i.e., $A_{L_{i,11}}={\mathrm{T}}_L$ and $A_{U_{i,11}}={\mathrm{T}}_U$. This therefore represents a shared risk model. If $k=21$ then $r_{21}=1$. This represents a model of minimization of global consumer's risk, so that $A_{L_{i,21}}={\mathrm{A}}_{{\mathrm{L}}_{+}}$ and $A_{U_{i,21}}={\mathrm{A}}_{U_{-}}$. To maintain simplicity, only the value of the multiplicative factor $r$ will be provided in the subsequent text without referencing the index $k$. The plane $\alpha$ and the reference plane $\beta$ together with the tolerance planes ${\tau }_L$ and ${\tau }_U$ and the acceptance planes ${\mathcal{A}}_L\dots A_{L_{i,21}}=-506.29385$mm and ${{\mathcal{A}}_U\dots A}_{U_{i,21}}=-506.28785$ mm, $i=\mathrm{1,2},..,n$, are shown in Figure 3.

Figure 3 illustrates the model of minimizing global consumer risk with linearized geometry. A model of minimization of global producer's risk can be presented similarly. The line highlighted in red represents the intersection of planes $\alpha$ and $\beta$. The equation of that line is given by:

$$p\dots y=-0.9970845x-244.898.$$

(6)

The points of the line $p$ are the only points of the plane $\alpha$ that are exactly in the middle of the tolerance field and the acceptance field [49]. It should be noticed that the plane α on the domain defined by the plane in which the sampled points $\left(x_i,y_i\right),i=1,\dots ,n$ lie does not intersect the tolerance planes or acceptance planes anywhere. This indicates that the model will not exhibit anomalies that may arise when such a mutual relationship between the tolerance and acceptance planes and the regression plane is present [52].

2.5. Risk Calculation

Considering the previously defined model parameters, assuming that the prior function and the likelihood function are normally distributed, the global consumer's risk $R_C$ in linear multivariate regression case, at each measurement point $\left(x_i,y_i\right),i=1,\dots ,n$, for each associated best estimate value $z_i,i=\mathrm{1,2},\dots ,n$, can be calculated from the equation:

$${\mathrm{R}}_{{\mathrm{C}}_{i,k}}=\underset{-\mathrm{\infty }}{\overset{{\displaystyle \frac{T_{L_i}-z_i}{{\mathrm{u}}_0}}}{\int }}{\mathsf{\phi }}_0\left(\mathrm{t}\right){\mathrm{F}}_{i,k}\left(\mathrm{t}\right)\mathrm{d}\mathrm{t}+\underset{{\displaystyle \frac{{\mathrm{T}}_{U_i}-z_i}{{\mathrm{u}}_0}}}{\overset{\mathrm{\infty }}{\int }}{\mathsf{\phi }}_0\left(\mathrm{t}\right){\mathrm{F}}_{i,k}\left(\mathrm{t}\right)\mathrm{d}\mathrm{t},i=\mathrm{1,2},...,n,k=\mathrm{1,2},...,m.$$

(7)

The global producer's risk $R_P$ is calculated from the equation:

$${\mathrm{R}}_{{\mathrm{P}}_{i,k}}=\underset{{\displaystyle \frac{{\mathrm{T}}_{L_i}-z_i}{{\mathrm{u}}_0}}}{\overset{{\displaystyle \frac{{\mathrm{T}}_{U_i}-z_i}{{\mathrm{u}}_0}}}{\int }}{\mathsf{\phi }}_0\left(\mathrm{t}\right)\left(1-{\mathrm{F}}_{i,k}\left(\mathrm{t}\right)\right)\mathrm{d}\mathrm{t},i=\mathrm{1,2},\dots ,n,k=\mathrm{1,2},\dots ,m.$$

(8)

In equations (7–8) ${\mathsf{\phi }}_0$ denotes the standard normal probability density function (PDF). The expression ${\mathrm{F}}_{i,k}\left(\mathrm{t}\right)$ is given by the equation:

$${\mathrm{F}}_{i,k}\left(\mathrm{t}\right)=\mathsf{\varphi }\left({\displaystyle \frac{{\mathrm{A}}_{U_{i,k}}-z_i-{\mathrm{u}}_0\mathrm{t}}{{\mathrm{u}}_{\mathrm{m}}}}\right)-\mathsf{\varphi }\left({\displaystyle \frac{{\mathrm{A}}_{L_{i,k}}-z_i-{\mathrm{u}}_0\mathrm{t}}{{\mathrm{u}}_{\mathrm{m}}}}\right),i=\mathrm{1,2},\dots ,n,k=\mathrm{1,2},\dots ,m.$$

(9)

The symbol $\mathsf{\varphi }$ in equation (9) stands for the cumulative distribution function for the normal distribution (CDF). Equation (9) and the integrals in equations (7–8) were obtained by transforming the argument of the prior in the point-based risk assessment method and by discretizing the constructed domain [19]. In transforming the prior argument, the functional connection of the data, given by equation (1), is considered. Since the function from equation (1) is continuous at the points $\left(x_i,y_i\right),i=1,\dots ,n$ and $\mathsf{\varphi }$ is continuous in the arguments from equation (9), the data obtained from equations (7–8) must also be functionally connected by a continuous function. The integrals from equations (7–8) were solved numerically. This was accomplished using packages from the R programming language [65,66]. The graphs were created using either R or Octave software [67].

3. Results

3.1. Risk Curves

Although it may happen that the measurement values for the flatness of the surface at points $\left(x_i,y_i\right),i=1,\dots ,n$ are equal, the resulting risks do not have to be equal. The values of global consumer's risk shown in Figure 4a correspond to the measured surface flatness values listed in Table 1 of [56] for $r=1$. The flatness of the surface of $z=-506.2903$ mm was measured at three different locations on the worktop of the CMM table, but the associated risks for these measurements, highlighted in red, vary. The representation in Figure 4a is given for the measured values of surface flatness before smoothing the data with the regression plane $\alpha$. When the risk is presented considering the functional relationship of the data described by the equation of the plane $\alpha$, it becomes evident that the risk values highly depend on the position of the sample points $\left(x_i,y_i\right),i=1,\dots ,n$. The dependence is expressed relative to the line $p$ from equation (6), and relative to the plane $\beta$ in which the line $p$ lies. Figures 4c and 4d illustrate the curves of global consumers' and producers' risk along the axis representing the values of the regression plane α calculated at the points $\left(x_i,y_i\right),i=1,\dots ,n$. These curves are upward-opening parabolas which have a minimum at a point that matches the value of the reference plane β from equation (2). Points that are equidistant from the axis of symmetry of the parabola, i.e., from the value specified in equation (2), have equal risks.

The curve that represents the global producer's risk increases along the guard band axis, while the curve of the global consumer's risk decreases along the guard band axis, Figure 4b. The curves shown are calculated for the value of the plane $\beta$. For that value, the risks are the lowest, so the curves in Figure 4b are referred to as curves of minimum. The curves of minimum intersect each other for $w_s=-0.44\mathsf{\mu }\mathrm{m}$. The risk values are then $R_C=R_P\approx 2.03\%$. On these curves, the points corresponding to the risk values when $r=1$, $r=0$ and $r=-1$ are particularly highlighted. The highlighted points have the same value as the points of minimum of the parabolas shown in Figures 4c, d. If $r\in \left[-1,\right.0⟩$, the guard band width values $w$ of the guard band width shown in Figure 4b are negative, corresponding to the conditions for a model that minimizes the global producer's risk. The values $r=0$ and $w=0$ indicate that it is a shared risk model. For $r\in ⟨0,\left.1\right]$, it holds that $w>0$. This represents the conditions for the model of minimization of global consumer's risk.

3.2. Risk Surfaces

The risk surfaces displayed in Figure 5a are the outcome of the three-dimensional extension of the standard, point-based method for risk assessment.

The risk surfaces presented in Figure 5a are displayed within the domain defined by the interval $I_w=\left[-0.75,0.75\right]\mathsf{\mu }\mathrm{m}$ on the guard band axis, and the range of measured values for surface flatness $I_z=\left[-506.292,-506.289\right]\mathrm{m}\mathrm{m}$. The maximum values of both the global producer's risk and global consumer's risk occur at the boundaries of interval $I_z$. The maximum values of the global producer's risk along the guard band axis are achieved for $w=0.75\mathsf{\mu }\mathrm{m}$. Simultaneously, the values of the global consumer's risk are minimal, Figure 5b. The risk surfaces in Figure 5a include the curves of minimum depicted in Figure 4b. All results analyzed so far obtained through the spatial risk assessment method are consistent with the results of previous research [43,51,52].

In Figure 4b, among others, the points for which $r=1$ are particularly highlighted. The global producer's risk along the guard band axis is then the highest and amounts to $R_{P1}=11.40\%$. Concurrently, the global consumer's risk along the guard band axis is at its minimum, with a value of $R_{C1}=0.45\%$. Figures 6a and 6b illustrate the risk surfaces for $r=1$. In contrast to the risk surfaces shown in Figure 5a, these surfaces are displayed within the $I_x\times I_y$ domain. The interval $I_x=\left[11.70031,984.7228\right]\mathrm{m}\mathrm{m}$ lies on the x-axis, while $I_y=\left[-1086.1684,-395.6763\right]\mathrm{m}\mathrm{m}$ is the interval on the y-axis. The lines $p_1$ and $p_2$ highlighted in Figure 6a and 6b pass through the points $\left(x,y,R_{C1}\right)$ and $\left(x,y,R_{p1}\right)$ respectively. The value of the variable $y$ for a given $x\in I_x$ is obtained from equation (6). The orthogonal projections of the lines $p_1$ and $p_2$ onto the xy-plane have the same direction and position as the projection of the line $p$ onto the xy-plane. This orthogonal projection of the line $p$ onto the xy-plane is hereinafter denoted by $p^{\text{'}}.$The risk of flatness deviation of the CMM worktable surface from the specified standards is smallest in the direction of line $p\text{'}$, Figure 6c. On the $I_x\times I_y$ domain, the risks are greatest at the point of the CMM workbench with coordinates $A\left(984.6031,-395.6776\right)$. At that point the risk values are $R_C\approx 0.89\%$ and $R_p\approx 14.30\%$.

The risk surfaces displayed in Figure 5a and Figure 6a,b have different domains. The depiction of the risk surface on the domain $I_w\times I_z$ enables the identification of the minimum curves along the guard band axis and the assessment of the surface flatness value that yields the lowest risks. The representation of the risk surfaces on the $I_x\times I_y$ domain has a practical application. This representation illustrates how to position the measurement object on the workbench surface so that the flatness of the surface and potential damage to the CMM device's workbench surface have the least impact on the measurement results.

3.3. Risk Spaces

Risk surfaces on the domain $I_x\times I_y$ can be obtained for each $k=\mathrm{1,2},\dots ,m$. The global producer's risk space and the global consumer's risk space are formed by stacking these layers on top of one another, Figure 7a,b.

The stacking of the consumer's global risk layers concerning the multiplicative factor $r$ is descending. For $r_l<r_s$,$l,s\in \left\{\mathrm{1,2},\dots ,m\right\},l\ne s,l<s$ holds $R_{Cl}>R_{Cs}$. The construction of the global producer's risk space is ascending. For $r_l<r_s$, $l,s\in \left\{\mathrm{1,2},\dots ,m\right\},l\ne s,l<s$ holds $R_{Pl}<R_{Ps}$. The minimum of each layer is along the direction of line $p$. The global consumer's risk space reaches its maximum at point A, with $r=-1$, and it has a value of $R_C=5.35\%$. As the behaviors of global producer's and consumer's risks are reversed, the maximum of the space of global risk of producers is reached at the same point, but for $r=1$ and amounts, as already stated, to $R_p=14.30\%$. It should be noticed that the spatial risk assessment method allows risk assessments for any guard band width and at locations where surface flatness measurements have not been made.

In this multivariate linear regression model, the output variable $z$ depends on two input variables $x$ and $y$. A representation of risk assessment in a model with more than two input variables would be given through combinations of pairs of variables.

3.4. Dependence of the Results on the Choice of the Reference Plane

This study presents the spatial risk assessment method for the first time. Therefore, there is currently no established professional consensus regarding the selection of the reference plane $\beta$. The choice of reference plane in this context is specific, as it is derived from the MZ method. The reference plane was determined based on the extreme values obtained from the surface flatness sample. In some problems, such values would be discarded as outliers and excluded. Thus, when selecting a reference plane for multiple linear regression models, it is advisable to use the sample mean or another relevant value, depending on the specific practical problem being addressed. In this case, there is no significant difference between the mean value of the data $\overline{z}=-506.290835$ mm and the value taken by plane $\beta .$ Therefore, the difference in risks is minimal, Figure 8d,e. When this disparity is substantial, there are also discernible differences between the estimated values for producers' and consumers' global risk.

Figure 8 illustrates the variations in values of global consumers' and producers' risk for $r=1$ and different choices of reference planes. The values of the reference planes are selected to match the characteristic points of the surface flatness sample: (a) minimum value, $z_{min}=-506.2926$ mm, (b) mode1, $z=-506.2915$ mm, (c) mode2, $z=-506.291$ mm (d) median, $z=-506.29095$ mm, (e) MZ method, the value taken by the plane $\beta ,$ $z=-506.29085$mm, (f) mean, $z=-506.290835$mm, (g) maximum value, $z_{max}=-506.2891$ mm. As the reference plane values shift from minimum to maximum, the region indicating the highest global consumer or producer risk transitions from point A on the workbench surface to the lower-left corner. It is essential to highlight the centrally symmetrical relationship between the graphs when the minimum and maximum surface flatness values are used as reference planes, Figure 8a and 8g. The least-risk direction and position correspond to case (e), where the reference plane is derived from the MZ method. This direction and position also correspond to the direction and position of the line $p^{\text{'}}.$ Considering a line that would be the intersection of plane α and any chosen reference plane, its projection onto the xy-plane would be parallel to line $p^{\text{'}}$.

3.5. Dependence of Results on Measurement Uncertainty $u_m$

Risk calculation requires information about the prior distribution and the measurement uncertainty, $u_m$, of a future measurement process. The spatial risk assessment method allows for the prediction of the risk value for the desired measurement uncertainty $u_m$ even without performing a measurement. All results so far are calculated for $u_m=u_0/2$. Figure 9 shows a comparison of the risk space from Figure 7 with the spaces obtained for $u_m=u_0$ and $u_m=2u_0.$ To facilitate comparison, the risk space is projected onto a plane to which the line $p$ is perpendicular. The observed lattice of space is limited by the layers obtained for $r=-1$ and $r=1$. The vertical lines on the graphs correspond to the lines connecting the layer vertices.

The space of global consumer’s risk contracts as measurement uncertainty $u_m$ increases. The values of the global consumer's risk for the layer of space where $r=1$ increase. Simultaneously, the values of global consumer's risk for the risk layer where $r=-1$ are declining. As a result of the change in the value of measurement uncertainty $u_m$, the increase in the risk values of the lower layer is significantly more pronounced than the decline in the risk value of the upper layer of the space of global consumer’s risk, Figure 9.

The contraction of the space of global consumer’s risk cannot be observed independently of the behavior of the space of global producer’s risk. With the increase in measurement uncertainty $u_m$, there is a translation of the space of global producer’s risk along the risk axis toward higher values, Figure 9. The space of the global producer’s risk expands until the point of alignment of the risk layers, that is, until the point at which the risk layers are almost parallel to the xy-plane. The critical value at which a change in the shape of the layers occurs is $u_k\approx 1.2842u0=2.44$ μm. That is approximately $32.53\%∆T$ in terms of the tolerance field range. After the critical value, the layers of the space of global producer’s risk change shape in projection, from an upward-opening parabola to a downward-opening parabola. The space of global consumer’s risk is contracting, while risk values are rising. Although for $u_m\ge u_k$ the values of the global consumer's risk are below $R_C=5.4\%$, the values of the global producer's risk are extremely high and go above 30%. Due to such unfavorable results for the producer, models that have a measurement uncertainty greater than the critical value should not be considered. If the method were used for quality control in the production of CMM granite worktables, under such conditions, the producer would have numerous non-compliant worktables and extremely high production costs.

4. Model Evaluation

4.1. Conformance Probability

To evaluate the model, the conformance probability $p_C$ is used. This is the probability that the item of interest, in this case, the flatness of the surface, is within the specification limits. For a measured value ${ɳ}_m$, outside the acceptance interval, conformance probability represents specific producer’s risk [19]. In machine learning terms, conformance probability represents the prevalence [68].

Evaluation of the model is required because the conformance probability may be small despite the risk values being small. These anomalies arise when a regression line or plane intersects a tolerance line or plane [52]. According to the two-dimensional interpretation of the risk assessment procedure depicted in Figure 2, the expression for calculating the conformance probability presented in [19] (p. 27) should be slightly redefined. Instead of the expression:

$$p_C={\int }_{T_{L_i}}^{T_{U_i}}{g}_0\left(\eta \right)d\eta ,i=\mathrm{1,2},\dots ,n,$$

(10)

that can lead to misunderstanding, the next expression should be observed:

$$p_C={\int }_{T_{L_i}}^{T_{U_i}}{\int }_{A_{L_{i,k}}}^{A_{U_{i,k}}}{g}_0\left(\eta \right)h\left({\eta }_m|\eta \right)d\eta d{\eta }_m,i=\mathrm{1,2},\dots ,n,k=\mathrm{1,2},\dots ,m.$$

(11)

Applying transformations that are equivalent to those that are implemented in the calculation of the global producer's and consumer's risk [19,51], equation (11) takes the form:

$$p_C=\underset{{\displaystyle \frac{T_{L_i}-y_0}{u_0}}}{\overset{{\displaystyle \frac{T_{U_i}-y_0}{u_0}}}{\int }}{\mathsf{\phi }}_0\left(\mathrm{t}\right){\mathrm{F}}_{i,k}\left(\mathrm{t}\right)\mathrm{d}\mathrm{t},\mathrm{i}=\mathrm{1,2},\dots ,\mathrm{n},\mathrm{k}=\mathrm{1,2},\dots ,\mathrm{m}.$$

(12)

Assuming that $A_{\mathrm{I}}=\mathbb{R}$, for $z_i,u_0,u_m<\infty$ and by applying the following properties related to CDF:

$$\underset{A_{L_{i,k}}\to -\infty }{\lim }\mathsf{\varphi }\left({\displaystyle \frac{A_{L_{i,k}}-z_i-{\mathrm{u}}_0\mathrm{t}}{{\mathrm{u}}_{\mathrm{m}}}}\right)=0,$$

(13)

$$\underset{A_{U_{i,k}}\to \infty }{\lim }\mathsf{\varphi }\left({\displaystyle \frac{A_{U_{i,k}}-z_i-{\mathrm{u}}_0\mathrm{t}}{{\mathrm{u}}_{\mathrm{m}}}}\right)=1,$$

(14)

it is obtained that equation (10) is valid for each $i=\mathrm{1,2},\dots ,n$ and for each $k=\mathrm{1,2},\dots ,m$. From equation (10), it is obvious that the value for the conformance probability depends on the best estimate $y_0$ and its associated measurement uncertainty $u_0$. Each value $z_i,i=\mathrm{1,2},\dots ,n$ has the same width of tolerance space. Therefore, the values of conformance probability over the domain $I_w\times I_z$ are the same for each $r_k,k=\mathrm{1,2},\dots ,m$, i.e., for each width of the guard band $w_k$. The same applies to the values of conformance probability over the domain $I_x\times I_y$, Figure 10.

The line of maximum values on the domain $I_w\times I_z$ passes through the points $\left(w_k,z_r,p_{Cmax}\right),k=\mathrm{1,2},\dots ,m$, where $p_{Cmax}\approx 95.15\%$. The line of maximum values on $I_x\times I_y$ domain follows the CMM workbench's least-risk direction of line $p^{\text{'}}$and has the value $p_{Cmax}$. The line of $p_{Cmin}\approx 88.62\%$ displays the lowest value in the $I_w\times I_z$ domain. On $I_x\times I_y$ domain, the minimum $p_{CA}=p_{Cmin}$ is achieved at point A, where the risks are at their highest.

4.2. Confusion Matrix

The lines $T_U$, $T_L$, $A_L$, and $A_U$ partition the area, representing possible values of T_V and M_V, into nine distinct fields, Figure 11a [69].

From Figure 11a, the definition areas for the quantities shown in Figure 2 can be reconstructed. Given that, as well as the definitions of $R_C$, $R_P$, TP and TN probabilities, it is simple to demonstrate that $p_C=TP+R_P$ and $1-p_c=R_c+TN$. It also held that $TP+TN+R_C+R_P=1$. The probabilities $R_C$, $R_P$, TP and TN represent the classes of the confusion matrix, Figure 11b [35,51]. By using the mentioned connections between the global producer's and consumer's risk, and conformance probability, in addition to the spatial representation of $R_C$ and $R_P$ classes, as shown in Figure 7, it is possible to construct spaces of classes TP and TN (Figure S1, S2).

4.3. Probability of Frequent and Rare Events

It is clear from the confusion matrix that another marginal probability can be constructed in addition to the conformance probability. This is the probability of frequent events $p_0=TP+R_C$. It represents the probability that the measured value is within the acceptance interval. The probability of rare events is defined as $1-p_0=TN+R_P$ [69]. Assuming $T_{\mathrm{I}}=\mathbb{R}$ and applying the same procedure that shows that the expressions from equations (10) and (11) are equivalent, the following equation is obtained:

$$p_{0_{i,k}}=\underset{-\infty }{\overset{\infty }{\int }}{\phi }_0\left(t\right)F_{i,k}\left(t\right)dt,i=\mathrm{1,2},\dots ,n,k=\mathrm{1,2},\dots ,m.$$

(15)

For each value of the multiplicative factor $r_k,k=\mathrm{1,2},\dots ,m$, a probability surface for frequent events is generated. By stacking these layers one on top of the other, the overall probability space for frequent events across the domain $I_x\times I_y$ is obtained, Figure 12a. The stacking of the layers of space, concerning the multiplicative factor $r$, is descending. For $r_l<r_s$, $l,s\in \left\{\mathrm{1,2},\dots ,m\right\},l\ne s,l<s$ holds $p_{0l}>p_{0s}$. The probability of frequent events can also be represented on the $I_w\times I_z$ domain, Figure 12b. The maximum probability of frequent events in both domains is $p_{0max}=96.58\%$. This value is achieved by the layer of the multiplicative factor $r=-1$, Figure 12a, i.e., under the conditions of the model of minimization of global producer's risk, Figure 12b.

A narrower spacing between the layers depicted in Figure 12a would be achieved by implementing a finer subdivision step of the interval $\left[-\mathrm{1,1}\right]$. At point A, where the risks are highest, the probability of frequent events is lowest and amounts to $p_{0A}=75.22\%$ on both domains. It should also be noted that, in machine learning terms, the probability of frequent events is also referred to as bias [68].

Given the stated range of values that have conformance probability and probability of frequent events, it can be inferred that the model performs satisfactorily on the provided data and that there are no anomalies in the model.

4.4. Curves, Surfaces and Spaces of Metrics Associated with the Confusion Matrix

The integration of all elements used in metrology for conformity assessment enables the evaluation of risk assessment models in multiple linear regression using metrics associated with confusion matrices [35,51].

Among the many tested metrics are the metrics that illustrate the relationship between the global producer's and global consumer's risk in relation to conformance and non-conformance probabilities, as well as the probabilities of frequent and rare events. In relation to the multiplicative factor $r$, the stacking order of layers for the False Positive Rate (FPR) and False Discovery Rate (FDR) metrics is a descending (Figure S3, S5). These metrics compare global consumer's risk with non-compliant measurement and probability of frequent events, respectively. Given that the curves of global consumer's risk along the guard band axis fall, the descending order of layers for these two metrics was expected. This confirms the previously mentioned results. In contrast, the False Negative Rate (FNR) and the False Omission Rate (FOR) are metrics that compare the global producer's risk with conformance probability, i.e., probability of rare events; thus, the ascending order of layers was obtained (Figure S4, S6). There are no intersections of layers for any of the metrics listed. Therefore, their behavior is considered stable for each $r\in \left[-\mathrm{1,1}\right].$The metrics for model evaluation that are most often used in practice are Accuracy, Precision, Recall and Specificity [70]. In the field of metrology, imbalanced data are the outcome of well-conducted measurement. Consequently, the probability of the TP class in the confusion matrix is the dominant one. For assessing imbalanced datasets, the F1-score metric serves as an effective evaluation tool. As the harmonic mean of Precision and Recall, the F1 score provides a balanced measure of these performance indicators [71,72]. Among the metrics listed, Precision, Recall and Specificity have stable behavior, while Accuracy and F1 score have unstable behavior (Figures S7–S11). Unstable behavior refers to the intersection of metrics layers.

Additionally, spaces for several less common metrics have also been constructed. Negative predictive value (NPV) and Markedness (MK) show stable behavior with respect to the intersection of metrics layers (Figures S12, S13). Balanced Accuracy (BA), G-mean, Bookmarked informedness (BM) and Cohen's kappa coefficient are metrics with unstable behavior (Figures S14–S17). The ranges for values of the multiplicative factor, $r_S$, for which stability in the behavior of tested metrics is achieved, are provided in the supplementary materials. These ranges are determined by observing all values of $r\in \left[-\mathrm{1,1}\right]$ and identifying a domain without layer intersection.

The most pronounced unstable behavior was observed with the MCC and DOR metrics. The surfaces of these metrics on the domain $I_w\times I_z$ are hyperbolic paraboloids, Figure 13a,b. A hyperbolic paraboloid is a translation surface that is created by moving one parabola over another. The surface of DOR is formed by translating a parabola with an opening upwards along a parabola with an opening downwards and have a classic saddle shape. The MCC surface is formed by translating a downward-opening parabola along an upward-opening one, with its orientation opposite to the DOR. In Figure 13a, b, translation parabolas are highlighted red. On the MCC surface, this parabola is a curve of the minimum, with a maximum value of ${MCC}_{max}=0.56$. This value is achieved for $r_{MCC}=-0.44$ and guard band width $w_{MCC}=-3.3·{10}^{-4}\mathrm{m}\mathrm{m}$. The values of the TNR and NPV metrics can be very small (Figures S11, S12), which results in relatively low values of the MCC metric on the domain $I_w\times I_z$, Figure 13c [68]. For the DOR metric, the parabola along which the translation is performed is the curve of the maximum with a minimum value ${DOR}_{min}=57.9$, achieved for $r_{DOR}=0.05$ and the guard band width of $w_{DOR}=3.75·{10}^{-5}\mathrm{m}\mathrm{m}$.

The surfaces of MCC and DOR metrics are generated in a variety of ways, and consequently, the creation of the space of these metrics on the domain $I_x\times I_y$ differs. For $r\in \left[-1,-0.44\right]$ the stacking of the MCC layers is ascending, i.e. for $r_l<r_s$, $l,s\in \left\{\mathrm{1,2},\dots ,m\right\},l\ne s,l<s$ holds ${\mathrm{m}\mathrm{i}\mathrm{n}\left(MCC\right)}_l<{\mathrm{m}\mathrm{i}\mathrm{n}\left(MCC\right)}_s$. For $r\in ⟨-0.44,\left.1\right]$ the values taken by the layers of the MCC metric are descending with respect to the multiplicative factor $r$ and it holds that ${\mathrm{m}\mathrm{i}\mathrm{n}\left(MCC\right)}_l>\mathrm{m}\mathrm{i}\mathrm{n}({MCC)}_s$. Here, the ascending and descending orders of the MCC minimums are discussed because the layers of MCC intersect at the edges of the domain $I_x\times I_y$ , Figure 14a, 14c.

The MCC layers' minima lie on the curve of minimum highlighted red, Figure 13a. The red vertical lines in Figures 16a and 16c connect the stability area, i.e. the area where there is no intersection of layers. The stability area is achieved for $r_S\in \left[0.45,1\right]$.

For the DOR metric, the situation is reversed. If it holds $r\in \left[-1,0.05\right]$, a descending trend in the values of the maximums of the DOR layers is clear. For $r\in ⟨0.05,\left.1\right]$, the maximums of the DOR layers are ascending. Therefore, the layers of the DOR metric also intersect each other, Figure 14b, d. The maximum values of the DOR layers lie on the curve of maximum highlighted red, Figure 13b. Unlike the MCC metric, if all $r\in \left[-\mathrm{1,1}\right]$ are considered, the DOR metric does not have a stable area on the domain $I_x\times I_y$, Figure 14d.

The stable behavior of the tested metrics on the $I_x\times I_y$ domain can be ensured by choosing the range of the multiplicative factor $r$, for which the injectivity property is fulfilled in the sense that $\forall r_1,r_2\in \left[-\mathrm{1,1}\right],\left(r_1\ne r_2\right)\Rightarrow M\left(r_1\right)\ne M\left(r_2\right).$ Labels $M\left(r_1\right)$ and $M\left(r_2\right)$ stands for values of the corresponding tested metrics calculated for given values of the multiplicative factor $r$. The values $r_{\mathrm{I}}$, of the multiplicative factor, which ensures the injectivity of the tested metrics, are displayed in Table 1.

The guard band widths from Table 1 are calculated from the following equations: $w_{-}=T_L-A_{L_{-}}=A_{U_{+}}-T_U$ and $w_{+}=A_{L_{+}}-T_L=T_U-A_{U_{-}}$. For the metrics Cohen’s kappa coefficient, MCC, and DOR, it is possible to determine two distinct intervals without intersection of the layers of the metrics. These intervals are separated and marked with L and R, Table 1. Both boundaries of the injectivity intervals, 6L, 7L and 8L, fall within the domain of the model of minimization for global producer's risk and represent the value of $A_{L_{-}}$. For that reason, their guard band widths $w_{-}$ are within the range given in Table 1, while the width of the guard band $w_{+}$ does not exist. Similarly, both limits of the interval 8R are within the domain of the model of minimization of consumer's global risk and represent the value of $A_{L_{+}}$. Consequently, the guard band $w_{-}$ does not exist, and the guard band $w_{+}$ is defined by the range provided in Table 1.

In addition to the injectivity intervals for individual metrics, it is also possible to determine the common injectivity interval, $r_C,$ of all tested metrics, where the layers of metrics do not intersect. Three possible outcomes are outlined in Table 2.

Intervals $I_1$ and $I_2$ are situated within the domain of the model for minimizing the global producer’s risk, while interval $I_3$ is in domain of the model for minimizing the global consumer’s risk. Since the curves of minimum for global producer's and consumer's risks intersect for $w_S=-0.44$ µm, Figure 4b, 5a, on the interval $I_2$ holds $R_C<R_P$. The same applies to interval $I_3$, while for interval $I_1$ it holds $R_C>R_P$. Choosing interval $I_1$ achieves the smallest values of global producer risk and choosing interval $I_3$ achieves the smallest values of global consumer risk. In this manner, optimal guard band widths were obtained.

5. Conclusions

This study introduces a spatial risk assessment method based on a multiple linear regression model and demonstrates its application in the case of CMM workbench flatness. This method enables calculation of global producer's and consumer's risks even for locations where no direct measurements are taken and for any value of the guard band width $w$. Through the case study, the least-risk area of the workbench surface was identified, and a critical measurement-uncertainty threshold was established, marking the point at which the global producer's risk landscape changes significantly. In addition, it is shown that by evaluating the model using metrics associated with confusion matrices, it is possible to determine the optimal width of the guard band.

With an appropriate selection of the reference plane, the method can be extended to risk assessment scenarios involving any input data type that can be represented through multiple linear regression. Consequently, this spatial risk assessment framework is applicable beyond metrology, offering potential utility across diverse scientific disciplines.