A Numerical Stability Analysis in the Inclusion of an Inverse Term in the Design of Experiments for Mixtures

Abstract

A mixture experiment is one where the response depends only on the relative proportions of the ingredients present in the mixture. Different regression models are used to analyze mixture experiments, such as the Scheffé model, the Slack Variable model, and models with inverse terms. Models with inverse terms are worthy of consideration in certain applications. These models have been analyzed considering their fit quality, but not their numerical stability. This article analyzes the numerical stability of the model with inverse terms and the use of pseudo components. Likewise, a criterion is defined for the selection of the regression model with inverse terms, based on the quality of fit and numerical stability.

1. Introduction

Many products are formed by mixing two or more ingredients. One or more properties of each product are generally of interest to the manufacturer, who is responsible for mixing the ingredients. In every case, the measured property of the final product depends on the percentages or proportions of the individual ingredients that are present in the formulation [1]. In the general mixture problem, the measured response is assumed to depend only on the proportions of the ingredients present in the mixture and not on the amount of the mixture. These proportions are connected by a linear constraint of the following type:

$$x_1+x_2+\dots +x_q=1$$

(1)

In some cases, the proportions of the mixture components may be subject to additional restrictions such as

$$a_i\le x_i\le b_i$$

(2)

for one or more components. Experiments with mixtures are generally modeled using Scheffé polynomial models [2]. The linear Scheffé model has the following form:

$$E\left(Y\right)=\sum _{i=1}^q{\beta }_i{x}_i$$

(3)

Likewise, Scheffé’s quadratic model can be described as follows:

$$E\left(Y\right)=\sum _{i=1}^q{\beta }_i{x}_i+\sum _{i=1}^{q-1}\sum _{j=i+1}^q{\beta }_{ij}{x}_i{x}_j$$

(4)

where ${\beta }_i$ and ${\beta }_{ij}$ are unknown parameters that must be estimated, generally using least squares. Due to Equation (1), the form of the quadratic Scheffé model only contains linear terms and cross-product terms.

A general mixture model, in matrix terms, can be presented as $Y=X\beta +\in$ or $E(Y)=X\beta$. To estimate the parameters in the β matrix via least squares, the following expression is commonly used $\widehat{\beta }={\left(X^{\prime }X\right)}^{-1}{X}^{\prime }y$, where the covariance matrix is $V\left(\widehat{\beta }\right)={\left(X^{\prime }X\right)}^{-1}{\sigma }^2$. The vector of fitted values is given by $\widehat{y}=X\beta$ and the residual vector is $\in =y-\widehat{y}=y-X\widehat{\beta }$. The vector $\in$ is assumed to follow a normal distribution, that is, $\in ~N(0,{\sigma }^2)$. In regression analysis, the least squares method is widely employed for parameter estimation due to its simplicity and efficiency under certain conditions. However, several alternative methods offer advantages in specific scenarios. One such method is Maximum Likelihood Estimation (MLE), which estimates parameters by maximizing the likelihood function, provided the data follow a specific probability distribution. MLE is highly applicable and efficient, particularly in complex models or non-normal distributions, though it demands more computational resources compared to the least squares method. Robust regression techniques, including Huber regression and M-estimators, are designed to handle outliers and violations of assumptions such as homoscedasticity. These methods diminish the influence of extreme values on parameter estimates, making them more reliable when outliers are present. Ridge regression is another alternative, particularly useful when multicollinearity exists among the predictor variables. By introducing a regularization term (λ) to the least squares objective function, ridge regression shrinks the regression coefficients, thereby reducing overfitting, a feature especially valuable in high-dimensional datasets.

If an exact linear dependence between the columns of $X$ exist, that is, if there is a set of all non-zero ${c_j}^{\prime }s$ such that $\sum _{j=1}^p{c}_j{x}_j=0$, then the matrix $X$ has a rank inferior to $p$ (predictor variable), and the inverse of $X^{\prime }X$ does not exist. In this case, many software packages would give an error message and would not calculate the inverse. However, if the linear dependence is only approximate, it is $\sum _{j=1}^p{c}_j{x}_j\approx 0$, and then we have the condition usually identified as collinearity or ill-conditioning (some authors use the term multicollinearity) [3]. In this case, many software packages proceed to calculate ${\left(X^{\prime }X\right)}^{-1}$ without any signal foreseeing to the potential problem [4]. When there is the presence of ill-conditioning, computer routines used to calculate ${\left(X^{\prime }X\right)}^{-1}$ can give erroneous results; this means that the least squares solution may be incorrect. Even if ${\left(X^{\prime }X\right)}^{-1}$ is correct, the variance of ${\widehat{\beta }}^{\prime }s$, given in the diagonal terms in $V\left(\widehat{\beta }\right)={\left(X’X\right)}^{-1}{\sigma }^2$, can be inflated by ill-conditioning [5].

Prescott et al. (2002) [3] have shown that the quadratic Kronecker model is a quadratic model specification that is less susceptible to ill-conditioning compared to the Scheffé model. They concluded that the quadratic Kronecker model always reduces the maximum eigenvalue of the information matrix ($X^{\prime }X$), reducing ill-conditioning. When the presence of collinearity among the terms in the Scheffé model is a possibility, and the appearance of the complete model form is of concern, using the model known as the Slack Variable model makes sense. Reference to the Slack Variable model form appears in [4,6,7,8,9,10,11,12,13,14,15,16]. The justification for using the Slack Variable model is that it offers a lower degree of collinearity in terms of the adjusted model, which represents less numerical instability in the estimation of the model coefficients [17,18]; likewise, it reduces the variances of individual estimated regression coefficients, reduces the correlation between the estimators, and makes the model less dependent on the precise location of the design points [4] have presented a literature review on Slack Variable models for mixture experiments.

Emphasizing the practical aspects of selecting appropriate mixture designs, numerous techniques have been developed for both the design and analysis of mixture experiments. These advancements focus on constructing and interpreting mixture models as well as analyzing mixture data. While substantial progress has been made in designing experiments for Scheffé models [19], the optimal design for more general blending models, particularly those extending the Becker class, has not yet been systematically addressed. These models, which account for nonlinear blending effects, are crucial in practical applications where the relationships between response variables, parameters, and mixture factors are more complex. To address this gap, ref. [20] proposed methods for determining D- and A-optimal designs for general blending models. By reformulating the optimal experimental design problem into an optimization problem—using either Nonlinear Programming or Mixed Integer Nonlinear Programming—their approach is applicable to quadratic and special cubic blending models, including those in the H2 class introduced by Becker. These methods are exemplified through applications in combustion science and fuel property characterization. Furthermore, a new class of models has been introduced to unify and extend existing statistical methodologies for analyzing mixture experiments [21]. This class integrates the Scheffé and Becker models, thereby broadening the range of mixture component effects that can be captured. These models maintain simplicity when appropriate but are flexible enough to accommodate more complex phenomena as required. The proposed methodology significantly extends traditional approaches, with supplementary material available for further exploration.

Draper and John (1977) [22] have presented another model alternative that proposes the inclusion of an inverse term. Models that include an inverse term are augmentations of the Scheffé polynomial with the additional terms of the form $x_i^{-1}$, included to account for the possible extreme change in the response as $x_i$ approaches zero [22]. Draper and John (1977) [22] have argued that the use of a model with inverse terms offers a better quality of fit than the Scheffé polynomial model, when an extreme change in the response behavior exists.

When any component of the mixture is equal to zero ($x_i=0$), to be able to include an inverse term such as $x_i^{-1}$ in the model, it is necessary to add a small positive quantity to each component of the mixture ($x_i$), say $c_i$ [22]. Draper and John (1977) [22] suggested that a general rule for defining the value of $c_i$ is simply to use a value between 0.02 and 0.05. In addition, they suggest that $c_1=c_2=\cdots =c_q=c$. Draper and John (1977) [22] have mentioned that values of $c_i$ in that range are suitable for most problems. It is important to highlight that the determination of the value of $c_i$ is not defined formally, or based on some mathematical criterion. The value of $c_i$ is completely at the discretion of the researcher, without having any reference or criteria. Therefore, the definition of $c_i$, through some formal criteria, is a research opportunity.

Although the model with inverse terms can offer a better fit, in this article, we show that this model has certain disadvantages regarding the conditioning of the information matrix used for estimating the parameters by the least squares method. Derived from the above, a strategy is proposed to minimize the problems derived from the poor conditioning of the information matrix. This strategy consists of determining the value of $c_i$ for which the model presents better numerical stability.

As a result of this research, it was observed that the conditioning of the information matrix, and consequently the numerical stability of the model, varies with different values of $c_i$. Therefore, it is proposed that the value of $c_i$ should be determined based on the model’s numerical stability. Specifically, the value of $c_i$ should be selected to optimize the conditioning of the information matrix, ensuring the most accurate least squares estimates.

2. Methods

In this section, we discuss various mathematical methods used in modeling mixture experiments, specifically focusing on the limitations and enhancements of the Scheffé polynomial model when dealing with extreme responses as certain components of the mixture approach boundary values (often zero). We also aim to introduce techniques such as pseudo components and how they can improve the numerical stability of the model, as well as provide a method to determine the adequacy of polynomial models in mixture designs.

In Section 2.1, we discuss how the Scheffé polynomial model can be modified using inverse terms to address extreme changes in the response when one or more components of a mixture approach zero. The equations provided extend the Scheffé model to account for these conditions, if no component value exactly reaches zero but can approach very small positive values. The need to handle these situations arises because the basic Scheffé model struggles when the response shows a drastic change near boundaries. Section 2.2 introduces the concept of pseudo components, which are used to simplify mixture models when components have lower limits (i.e., boundary constraints). The transformation of components into pseudo components helps to handle extreme changes in responses. The equations provided show how to transform the original component values into pseudo components to account for boundary effects, and this section emphasizes the need for adding a small value (denoted as $c_i$) when a component’s value is zero to stabilize the model. Next, Section 2.3 explains how to assess the adequacy of polynomial models used for mixtures in terms of the numerical stability and collinearity. It describes how to calculate the condition number (CN) of the information matrix to measure the numerical stability, with smaller CN values indicating better model conditioning. Additionally, it proposes the use of the CN to determine the appropriate value of $c_i$, which is essential for improving the stability of the pseudo component transformations. Finally, in Section 2.4, a discussion is presented on alternative problem formulations and regularization methods.

2.1. Mixture Model with Inverse Terms

As mentioned before, the Scheffé model can be re-parameterized; in fact, there are several ways to write a polynomial model of any order using the mixture restriction of Equation (1) [3].

When there exists an extreme change in the response as the value of a certain component ($x_i)$ tends to a boundary (usually zero), the Scheffé model cannot cope well with this situation [6]. For modeling extreme changes in the response as the value of some components tends to zero, [22] proposed the following models:

$$E\left(Y\right)=\sum _{i=1}^q{\beta }_i{x}_i+\sum _{i=1}^q{\beta }_i{x}_i^{-1},$$

(5)

and

$$E\left(Y\right)=\sum _{i=1}^q{\beta }_i{x}_i+\sum \sum _{i<j}^q{\beta }_{ij}{x}_i{x}_j+\sum _{i=1}^q{\beta }_i{x}_i^{-1}.$$

(6)

Equations (5) and (6) are a modification of the Scheffé polynomial model with the additional term of the form $x_i^{-1}$. From the above, it is assumed that the value of $x_i$ never reaches zero; however, the value can be very close to zero, that is, $x_i\to {\epsilon }_i>0$, where ${\epsilon }_i$ is some small quantity, which is defined for each application of this model. Likewise, if only a couple of components are as likely to produce extreme changes in the response as $x_i\to 0$, then only these terms are included in the model [1].

2.2. Use of Pseudo Components

When lower limit restrictions like those shown below are considered, $0\le L_i\le x_i$ for all $i$ and $\sum _{i=1}^q{L}_i<1$ lower pseudo components (L-pseudo components) are suggested for use in place of components in original units ($x_i$) [1]:

$${x^{\prime }}_i=(x_i-L_i)/(1-\sum _{i=1}^q{L}_i)$$

(7)

If an extreme change in the response occurs as $x_i$ approaches the bound $L_i$ (smallest component value $x_i$), we define ${L^{\prime }}_i=L_i-{\epsilon }_i$, where $L_i>{\epsilon }_i>0$, so instead of the previous lower pseudo component definition shown in Equation (7), the transformation is redefined as follows [1]:

$${x^{\prime }}_i={\displaystyle \frac{x_i-{L^{\prime }}_i}{1-\sum _{i=1}^q{L^{\prime }}_i}}$$

(8)

Draper and John (1977) [22] mentioned that when the value of any component of the mixture is equal to zero ($x_i=0$), then in order to include a term such as $x_i^{-1}$ in the model, a small positive amount must be added, say $c_i$, to each value of $x_i$. This is equivalent to working again with the pseudo components ${x_i}^{\prime }$, where this time, they are defined as follows [1]:

$${x^{\prime }}_i={\displaystyle \frac{{(x}_i+c_i)-{L^{\prime }}_i}{1-\sum _{i=1}^q({L^{\prime }}_i-c_i)}}$$

(9)

As mentioned in the Introduction, the value of $c_i$ has not been defined by a formal criterion, leaving it to the discretion of the researcher. In the Results and Discussion, we will describe how to define the value of $c_i$ so that an improvement in the numerical stability of the model is obtained.

2.3. Adequacy of Polynomial Models for Mixtures

Next, we describe the procedure for determining the adequacy of the polynomial model used in terms of the conditioning of the information matrix. In this case, the objective is to determine the presence of collinearity (see the works of Cornell and Gorman (2003) [23] and Prescott et al. (2002) [3]).

Let ${\lambda }_{max}>{\lambda }_2>\dots {\lambda }_{p-1}>{\lambda }_{min}$ be the $p$ eigenvalue of the information matrix ($X^{\prime }X$), that is, the $p$ solutions to the determinant equation $\left|X^{\prime }X-\lambda I\right|=0$, which is polynomial with $p$ roots.

There are several definitions of the conditional number (CN) of a matrix. The definition generally used in statistical applications is the square root of the ratio of the maximum to the minimum eigenvalue of $X^{\prime }X$, denoted by the following [3]:

$$CN\left(X^{\prime }X\right)=\sqrt{{\displaystyle \frac{{\lambda }_{max}}{{\lambda }_{min}}}}$$

(10)

(Some authors omit the square root and just take the ratio of the maximum to the minimum eigenvalue). If $X^{\prime }X$ is close to singular, then ${\lambda }_{min}$ is close to 0 and $CN\left(X^{\prime }X\right)$ is extremely large. Small values of ${\lambda }_{min}$ and large values of ${\lambda }_{max}$ indicate the presence of collinearity. A smaller CN indicates more stability (better conditioning) in the least squares estimates.

This article uses the CN to determine the level of numerical stability of the model. Likewise, it is proposed to define the value of $c_i$ based on the CN of the information matrix. That is, the value of $c_i$ which results in the lowest value of CN is the value that will be used for the transformation to pseudo components of the form (9).

When there is no constant column in the $X$ matrix, the adjusted total sum of the squares in the denominator of the multiple correlation coefficient (${\mathrm{R}}_{\mathrm{j}}^2)$ must be replaced [3]. We define $x_i$ as the j-th column of $X$ and ${\mathrm{X}}_{\mathrm{j}}$ as the matrix that results when column ${\mathrm{x}}_{\mathrm{j}}$ is deleted from $\mathrm{X}$. Then, ${\mathrm{R}}_{\mathrm{j}}^2$ is the multiple correlation coefficient obtained by regressing ${\mathrm{x}}_{\mathrm{j}}$ on $X_{\mathrm{j}}$. When the first column of $X$ is a constant column, ${\mathrm{R}}_{\mathrm{j}}^2$ is usually calculated, for $j=2,\dots ,p,$ as

$$R^2={\displaystyle \frac{x_j^{\prime }{X}_j{\left(X_j^{\prime }{X}_j\right)}^{-1}{X}_j^{\prime }{x}_j-x_j^{\prime }{11}^{\prime }{x}_j/n}{x_j^{\prime }{x}_j-x_j^{\prime }{11}^{\prime }{x}_j/n}},$$

When there is no constant column in the $X$ matrix, the unadjusted multiple correlation coefficient can be obtained by

$$R_j^2={\displaystyle \frac{x_j^{\prime }{X}_j{\left(X_j^{\prime }{X}_j\right)}^{-1}{X}_j^{\prime }{x}_j}{x_j^{\prime }{x}_j}},forj=1,\dots ,p$$

2.4. Alternative Problem Formulations and Regularization Methods

While model reformulation, such as the inclusion of inverse terms or the Slack Variable model, can help mitigate collinearity and improve stability, there are additional approaches that can further enhance the conditioning of the Fisher information matrix and ensure more reliable results. One of the primary concerns with mixture models, especially in the presence of ill-conditioning, is the direct inversion of the matrix $X^{\prime }X$, which may fail in the case of near-linear dependencies between the predictor variables. In such cases, reformulating the problem can be highly beneficial. For instance, instead of solving the linear equations directly, regularization methods can be introduced to modify the objective function and reduce the influence of poorly conditioned parameters [24].

As a first step, ridge regression can be applied to the least squares problem. This method involves adding a penalty term proportional to the square of the regression coefficients to the least squares objective function. The regularization parameter (λ) effectively shrinks the coefficients, which helps stabilize the inversion of the $X^{\prime }X$ matrix, improving the conditioning of the Fisher information matrix [25]. By incorporating ridge regression into the mixture model estimation, the overall numerical stability of the model is enhanced, particularly in the presence of collinearity among the mixture components. Alternatively, Lasso regression, which uses L1 regularization, can also help in reducing collinearity issues. This method introduces a sparsity constraint, where some coefficients are driven to zero, thereby effectively reducing the number of predictors and mitigating the risk of ill-conditioning. In the context of mixture models, where the components may be highly correlated, Lasso can help identify the most influential predictors, ensuring better stability and interpretability [26].

In addition to ridge and Lasso regression, robust regression techniques, such as Huber regression and M-estimators, can also be employed to enhance the numerical stability of mixture models. These methods are designed to handle outliers and deviations from the assumptions of homoscedasticity, which can improve the overall conditioning of the Fisher matrix by reducing the influence of extreme values. For mixture experiments where the response is sensitive to outliers or noise, robust regression techniques are highly recommended [27]. Moreover, quantile regression could be considered in cases where the central tendency of the data is not fully representative of the underlying process. Unlike the ordinary least squares method, which focuses on the mean of the distribution, quantile regression estimates conditional quantiles, offering a more flexible and robust approach that can improve the robustness and conditioning of the estimated parameters [28].

The conditioning of the Fisher information matrix can also be influenced by the choice of optimality criteria used to fit the model. For example, instead of minimizing the sum of squared residuals, more robust or alternative loss functions can be considered. Huber loss, for example, is less sensitive to outliers and provides a more stable estimation when the response is prone to extreme values or non-normal errors. By minimizing the Huber loss, the impact of outliers is reduced, which in turn helps in improving the conditioning of the Fisher matrix. Additionally, alternative design criteria, such as D-optimality or A-optimality, can be applied when selecting mixture designs. These criteria aim to maximize the information provided by the experimental design, improving the conditioning of the Fisher information matrix by ensuring that the mixture components are well spread and sufficiently informative. By optimizing the experimental design based on these criteria, it is possible to achieve better numerical stability and more reliable parameter estimates in the presence of collinearity [29].

3. Results and Discussion

In this section, we illustrate how to define the value of $c_i$ based on the CN. For this, we use an example described in [1] regarding the octane rating of different gasoline blends.

Example of Octane Number of Gasoline

In a study, the octane number of nine mixtures consisting of olefin ($x_1$), an aromatic ($x_2$), and a saturator ($x_3$) at 1.5 mL of lead per gallon was used as an example. The octane number and composition of the mixtures in their original units and coded in pseudo components, are shown in

Table 1. Octane numbers in original units and in pseudo components.

	Original Units			Coding in Pseudo Components			Octanes
	${\mathit{x}}_1$	${\mathit{x}}_2$	${\mathit{x}}_3$	${{\mathit{x}}^{\prime }}_1$	${{\mathit{x}}^{\prime }}_2$	${{\mathit{x}}^{\prime }}_3$	y
1	0.010	0.870	0.120	0.023	0.875	0.102	111.5
2	0.541	0.000	0.459	0.545	0.02	0.435	101.3
3	0.427	0.061	0.512	0.433	0.08	0.487	80.6
4	0.022	0.464	0.514	0.034	0.476	0.49	91.0
5	0.007	0.957	0.036	0.02	0.961	0.019	107.0
6	0.414	0.278	0.308	0.42	0.293	0.287	97.0
7	0.648	0.030	0.322	0.65	0.049	0.301	98.6
8	0.162	0.514	0.324	0.172	0.525	0.303	92.2
9	0.008	0.068	0.924	0.021	0.087	0.892	77.8
10	0.010	0.870	0.120	0.023	0.875	0.102	111.5

.

As can be noted in Table 1, this experiment presents lower limits. In these cases, the literature recommends using pseudo components of the form shown in Equation (7). However, point 2 in Table 1 shows that the lower limit of$x_2$ is equal to zero. The above would make it impossible to include an inverse term of the form $x_i^{-1}$ for a transformation on pseudo components shown in Equation (9). The solution proposed by [22] consisted of selecting a lower limit for $x_i^{\prime }$, slightly greater than zero, for example, 0.02. This value is the one corresponding to $c_i$.

Draper and John (1977) [22] recommend $c_i=0.02$. The transformation to pseudo components is carried out in the following way: ${x^{\prime }}_i={\displaystyle \frac{{(x}_i+c_i)-{L^{\prime }}_i}{1-\sum _{i=1}^q({L^{\prime }}_i-c_i)}}$.

By substituting $c_i=0.02$ and ${L^{\prime }}_i$, we have ${x^{\prime }}_i={\displaystyle \frac{{(x}_i+0.02)-L_i}{1.017}}$.

The transformation to pseudo components, for all values of $x_i^{\prime }$, is shown in Table 1.

The next six models were obtained following a sequential fit of the Scheffé polynomial models proposed in [22] with units in L-pseudo components:

Model 1: $E\left(y\right)={108.52x}_1+{109.46x}_2+{71.51x}_3$

Model 2: $E\left(y\right)=140.62x_1+112.26x_2+74.07x_3-87.79x_1{x}_2-71.46x_1{x}_3-1.04x_2{x}_3$

Model 3: $E\left(y\right)=117.25x_1+99.14x_2+61.62x_3+0.25x_1^{-1}$

Model 4: $E\left(y\right)=95.90x_1+111.28x_2+67.42x_3+0.34x_2^{-1}$

Model 5: $E\left(y\right)=95.4x_1+115.6x_2+65.7x_3+0.4x_2^{-1}-0.1x_3^{-1}$

Model 6: $E\left(y\right)=109.9x_1+107.0x_2+52.2x_3+0.3x_1^{-1}+0.3x_2^{-1}-0.3x_3^{-1}$

Table 2. Measurements of the quality of fit and numerical stability of the models.

Model	${\mathit{R}}_{\mathit{A}}^2$	${\mathit{S}}^2$	${\mathit{\lambda }}_{\mathit{m}\mathit{i}\mathit{n}}$	${\mathit{\lambda }}_{\mathit{m}\mathit{a}\mathit{x}}$	CN	p-Value
$1(c_i=$0.02)	0.65	42.5	0.50	3	2	0.000
$2(c_i=$0.02)	0.52	59.2	0.004	3	28	0.000
$3(c_i=$0.02)	0.69	38.4	0.20	7783	212	0.030
$4(c_i=$0.02)	0.79	25.4	0.40	3314	95	0.011
$5(c_i=$0.02)	0.78	27.1	0.40	2731	78	0.032
$6(c_i=$0.02)	0.88	14.8	0.10	9458	298	0.031

presents the values of the adjusted correlation coefficient $R_A^2$ and the mean square of the error $S^2$. Taking as reference the values of $R_A^2$ and $S^2$, the best fit is obtained with model 6, showing a value of $R_A^2=0.88$ and $S^2=14.8$. However, it can be seen in the last column of Table 2 that the value of the CN of model 4 offers a greater degree of numerical stability, since it has a value of CN equal to 95, while model 6 presents a CN equal to 298, more than double that of model 4. The single-inverse-term model, model 4, also appears to be a better fitted model than the Scheffé second-degree polynomial, model 2. Thus, with this set of data, we have an example where a second-degree equation does not improve the fit of the surface obtained with the first-degree model but the addition of inverse terms to the first-degree model does improve the fit.

Model 6 is selected, given that it presents a better fit based on the values of $R_A^2$ and $S^2$. Next, different values of $c_i$ are assigned in Equation (9). The result of this exercise, as well as the quality of fit and numerical stability measures of the model, are presented in

Table 3. Measures of the quality of fit and numerical stability of model 6.

$\mathbf{Model}6.\mathit{y}\left(\mathit{x}\right)=109.9{\mathit{x}}_1+107.0{\mathit{x}}_2+52.2{\mathit{x}}_3+0.3{\mathit{x}}_1^{-1}+0.3{\mathit{x}}_2^{-1}-0.3{\mathit{x}}_3^{-1}$
Model	${\mathit{R}}_{\mathit{A}}^2$	${\mathit{S}}^2$	${\mathit{\lambda }}_{\mathit{m}\mathit{i}\mathit{n}}$	${\mathit{\lambda }}_{\mathit{m}\mathit{a}\mathit{x}}$	CN	p-Value
6 ($c_i=$0.03)	0.87	16.6	0.10	1615	248	0.036
6 ($c_i=$0.04)	0.85	18.9	0.10	2797	224	0.044
6 ($c_i=$0.05)	0.83	21	0.10	1911	212	0.051
6 ($c_i=$0.06)	0.81	23	0.10	1407	204	0.058
6 ($c_i=$0.07)	0.80	25	0.00	1092	201	0.065
6 ($c_i=$0.08)	0.79	26	0.04	879	199	0.072

.

Table 2 shows that model 6 presents a CN value of 298, which suggests a certain degree of collinearity and, consequently, the poor numerical stability of the model. The above could result in an inflation of the variance in the model predictions. Likewise, Table 3 shows that as the value of $c_i$ increases, the CN decreases. In the case of the value of $c_i=0.07$, there is a satisfactory level of quality of fit with $R_A^2=0.80$ and $S^2=25$, achieving a decrease in the CN from 298 to 201. This represents an improvement in the numerical stability of the model, as well as a more formal criterion to determine the value of $c_i$.

Figure 1 shows a graphical representation of how both the quality of fit and the numerical stability of the model change with different values of $c_i$. As can be seen in Figure 1, as the value of $c_i$ increases, the numerical stability of the model (with respect to the value of CN) increases as well until reaching an inflection point, where the numerical stability of the model begins to worsen. In addition to the above, regarding the quality of fit, with higher values of $c_i$, worse values of $R_A^2$ and $S^2$ are obtained.

Now suppose that model 4 is selected since it is the model that presents greater numerical stability, based on the value of CN. Assigning different values of $c_i$ in Equation (9), the quality of fit and numerical stability measures of the model are obtained, as shown in

Table 4. Measures of the quality of fit and numerical stability of model 4.

$\mathbf{Model}4.\mathit{y}\left(\mathit{x}\right)=95.90{\mathit{x}}_1+111.28{\mathit{x}}_2+67.42{\mathit{x}}_3+0.34{\mathit{x}}_2^{-1}$
Model	${\mathit{R}}_{\mathit{A}}^2$	${\mathit{S}}^2$	${\mathit{\lambda }}_{\mathit{m}\mathit{i}\mathit{n}}$	${\mathit{\lambda }}_{\mathit{m}\mathit{a}\mathit{x}}$	CN	p-Value
4 ($c_i=$0.03)	0.79	25.7	0.30	1792	76	0.012
4 ($c_i=$ 0.04)	0.78	26.6	0.30	1199	68	0.013
4 ($c_i=$ 0.05)	0.77	27.7	0.20	893	63	0.014
4 ($c_i=$0.06)	0.76	28.9	0.19	710	60	0.015
4 ($c_i=$0.07)	0.76	29.9	0.16	589	59	0.017
4 ($c_i=$0.08)	0.75	30.9	0.14	504	58	0.018

.

Model 4 presents a value of CN equal to 95 (see Table 2), which, compared to the CN of model 6 (CN = 298), is much smaller. In Table 4, once again, it can be seen that the value of CN changes with different values of $c_i$. For $c_i=0.03$, there is an adequate quality of fit and a decrease in the value of CN (from 95 to 76). The above represents an improvement in the numerical stability of the model, and a more formal criterion for determining the value of $c_i$ as well.

Figure 2 shows a graphical representation of how both the quality of fit and the numerical stability of the model change with different values of $c_i$. In Figure 2, it is shown again that the numerical stability of the model changes with different values of $c_i$. The above can be used to decrease the value of the CN, and consequently, improve the numerical stability of the model.

4. Conclusions

In this article, it was found that the conditioning of the information matrix, which represents the numerical stability of the model adjusted by least squares, changes with different values of $c_i$ in the transformation into lower limit pseudo components. An information matrix with poor conditioning can result in highly correlated estimators with a high standard error. Likewise, the estimated model is highly dependent on the precise location of the experimental points, due to the inflation of the variance of the estimators. The above reflects the importance of taking into consideration the numerical stability of the model when determining the model to be used. For the data used in this article, the first-order inverse model is a useful mode to be considered. Likewise, it contains fewer terms.

We can conclude that the numerical stability of the model can be improved by defining the value of $c_i$, based on the CN. It is important to highlight that this determination offers a mathematical criterion for professionals who wish to include an inverse term in the regression model. Based on the results that we have derived in this article, we recommend its use in the mixture experiment setting both in design search algorithms and in model fitting.

There are several potential extensions of this work that could lead to further advancements in statistical analysis, particularly in the context of mixture experiments and regression models. Building on the findings related to the correlation of estimators and variance inflation, extensions of this work could integrate regularization techniques, such as ridge regression or Lasso, to further control for multicollinearity and stabilize parameter estimates. This could be particularly useful in complex mixture systems with many variables or when data are scarce. Investigating the robustness of the inverse models and their numerical stability in the presence of noisy or incomplete data could be another useful extension. Since experimental data are often imperfect, evaluating how well these models perform under different levels of noise or outliers could refine their applicability in real-world scenarios.