Effect of Measurement Error on the Multivariate CUSUM Control Chart for Compositional Data

1  Introduction

A widely used method is statistical process control (SPC) to monitor a process. SPC is a decision-making tool used in almost all industrial manufacturing processes to achieve process stability and sustainable quality product. W. A. Shewhart introduced the CCs in the 1920s. The goal of a CC is to find assignable causes of shifts and to minimize the production of a large number of defective items. In practice, determining which CC is appropriate for specific data can be difficult. It is possible to figure out by looking at the distribution of the underlying data process. CCs with compound or simple rules can be computed using the run-length distribution using a Markov chain approach [1]. Based on concomitant order statistics and long-run dependence, a new bivariate semiparametric CC is proposed by Koutras et al. [2]. Integrating techniques for preventive maintenance and quality control policies with the time-between-event chart for reliability-oriented design has been investigated by He et al. [3]. Recently under different classical estimators, control charts for shape parameters of power function distributions were studied by Zaka et al. [4]. Hybrid EWMA CCs are designed to monitor process dispersion efficiently for the automotive industry [5]. For trimethylchlorosilane purification, simulations and multivariate CCs are utilized to increase the effectiveness [6].

Multivariate EWMA CCs are proposed using multiple variables based on an AIC calculation for multivariate SPC [7]. The frequency of unblinded and blinded adverse events accumulating from single and multiple trials using CCs was monitored in [8]. The ARL comparison of the confidence interval chart using MEWMA CC using the inertial impact of estimating the actual value of the mean vector (MV) has been made in [9]. Further, the economic design of the MEWMA CC using the Markov chain method, a hybrid Taguchi loss method and a genetic algorithm approach is also used to obtain the optimal design [10]. Phase-II EWMA CC has been proposed to monitor exponentially distributed process based on Type-II censored data to detect the shift in location and scale parameters [11]. For big data, Ahmad et al. [12] studied the distribution of T2 statistics to analyze high dimensional data for SPC.

Compositional data (CoDa) vectors are positive components presented as percentages, ratios, proportions or parts of a whole. Aitchison’s principles and statistical tools from [13] have been used to solve many practical problems. For more information about CoDa, readers are referred to [14] and [15]. Because the CoDa variable aggregates are limited to constant values, they can not be used the same way as standard multivariate data [16]. CoDa has applications in various fields, including chemical research surveys, geology, engineering sciences and econometric data analysis. Some of the most recent articles dealing with statistical methods and CoDa processing are discussed here. An automobile market application in which the authors model brand share of the market as a function of media investments while taking account of the brand’s cost using CoDa regression has been examined [17]. The sample space’s algebraic-geometric structure was introduced to adhere to the principles of sample space and the structure of CoDa [18]. Recently, Zaidi et al. [19] investigated the impact of M.E’s on Hotelling T2 CoDa CC using Average run length (ARL) as a performance measure. The error model with covariate M.E is supposed to be linear. The MEWMA CC for CoDa has been studied by Tran et al. [20]. After that, Zaidi et al. [21] examined the impact of M.E’s on the MEWMA-CoDa CC’s based on the ARL performance.

The CUSUM charts are well-known for detecting smaller and more frequent changes than Shewhart charts and are among the most common and widely used charts in practice [22]. Risk monitoring for CUSUM CCs with model error suggested by Knoth et al. [23]. Further, Shu et al. [24] analyzed the CUSUM accounting process for monotone population changes. It is essential to optimize several multiple objective functions simultaneously when performing multi-objective optimization and the multivariate CC method is unaffected by small and moderate variable shifts [25]. Such multivariate charts are commonly used to detect changes in process parameters around their corresponding axes, which must be illustrated in a specific direction. The MCUSUM charts for detecting shifts in MV and VCM have been discussed in [26]. A double progressive mean CC for monitoring Poisson process observations has been discussed in [27]. Furthermore, Alevizakos et al. [28] proposed a triple moving average CC to improve the chart’s detection ability and found them to be more effective than the moving average and double moving average CCs in detecting small to moderate shifts.

CCs are often created with the hypothesis that the factor of the quality feature is assessed accurately. However, M.E’s can occur due to various circumstances. Many researchers examined the impact of M.E’s on CCs. Imprecise measurement instruments or human errors are two main causes contributing to M.E’s. The presence of M.E’s must be taken into account while implementing CCs because it has a huge impact on their performance [29]. A cumulative sum (CUSUM) CCs has been proposed for the drilling process to monitor the variations of discharge pulses for sensing the actual condition of the whole drilling process taking M.E’s into account [30]. For a finite horizon production process, Nguyen et al. [31] suggested one-sided Shewhart-type and EWMA T2 charts handle the ratio of two normal random variables (R.V). For the X and S2 charts, Linna et al. [32] presented LM with covariate M.E X=A+BY+e, where e is a random error (R.E) due to measurement imprecision. The model X=Y+e to investigate the effects of M.E’s on the Shewhart X chart has been proposed by Bennett [33], where Y represents the actual value of the quality characteristic, and X represents the analyzed value provided by the measurement equipment.

Typically, practitioners monitor processes without considering the possibility of M.E’s. Despite sophisticated measurement techniques and methodologies, M.E’s is commonly caused by operational and environmental factors. As Montgomery et al. [34] outlined, Gage capacity tests commonly quantify M.E’s variability. The impacts of M.E’s on the functionality of a multivariate CC for the combined monitoring of normal processes MV and VCM have been studied by Maleki et al. [35]. A thorough examination of the impact of M.E’s in SPC has been conducted by Maravelakis et al. [36]. The performance analysis of univariate adaptive CCs in the existence of M.E’s has been examined by Hu et al. [37]. Based on an additive covariate model, Maleki et al. [38] investigated the adverse effects of M.E’s on the detection performance of triple sampling X¯ control charts. Attribute control chart performance in the presence of M.E. to monitor the parameter of a modified power function distribution is studied by Zaka et al. [39]. Two-parameter exponential distribution with M.E.s can be monitored using a general weighted moving average CC [40]. Multivariate homogeneously weighted moving average charts were used to monitor the process mean in the presence of M.E’s [41].

Although many studies on the assessment of CCs in the presence of M.E’s have been conducted, few have focused on multivariate data, and very few researchers have contributed to CoDa. But as far as the authors are aware, none of them has been dedicated to MCUSUM CC for CoDa in the presence of M.E. Hence this paper aims to fill this gap and investigate the impact of M.E’s on the MCUSUM-CoDa CC. The paper is organized as follows: In Section 2, the modelling and transformation of CoDa are briefly explained. Then, Section 3 introduces the LM with covariate M.E for CoDa. Section 4 details the MCUSUM-CoDa CC in the presence of M.E’s, and Section 5 deals with the performance evaluation of the proposed CC. Finally, two very detailed illustrative examples are provided in Section 6, and conclusions and future research directions are presented in Section 7.

2  Compositional Data (CoDa)

A row vector y=(y1,…,yp) consists of p−part composition defined on the simplex space 𝒮p, known as CoDa, where,

𝒮p={y=(y1,y2,…,yp)|yi>0,i=1,2,…,p&∑i=1pyi=κ}.(1)

where κ is a positive constant. There are many possible values for κ, e.g., if the data of interest is in percentages, we will use κ=100, while if the data is in the form of probabilities or proportions, κ=1 will be used. All the CoDa vectors are supposed to be row vectors throughout this paper. One of the major properties of CoDa is relative property; two CoDa vectors are supposed to be equivalent if they carry the same relative information regardless of the actual values. The closure function of CoDa is defined as

𝒞(y)=(κy1∑i=1pyi,κy2∑i=1pyi,…,κyp∑i=1pyi).(2)

As CoDa deals with the constraint of constant sum, Euclidean geometry is not supposed to be useful for CoDa. Let us take an example to make it more clear, Suppose there are two CoDa vectors, y=(0.2,0.5,0.3)∈𝒮p and z=(0.1,0.25,0.15)∈𝒮p, when we add them using Euclidean geometry, we get y+z=(0.3,0.75,0.45)∈𝒮p. In another case, by multiplying CoDa with a scalar, we get 5×y=(1,2.5,1.5)∈𝒮p. So, it is concluded that using Euclidean geometry operators for CoDa is not a suitable option. To deal with CoDa, Aitchison [14] introduced a new geometry known as Aitchison geometry. These new operators defined by Aitchison are

• the perturbation operator ⊕ of y∈𝒮p by z∈𝒮p (for adding CoDa vectors) defined as

y⊕z=𝒞(y1z1,y2z2,…,ypzp),(3)

• the powering operator ⊙ of y∈𝒮p by a constant c∈R (for multiplying CoDa vectors with a scalar) defined as

c⊙y=𝒞(y1c,y2c,…,ypc).(4)

The constant sum constraint makes CoDa difficult to deal with in its original form, hence pre-defined log-ratio transformations (LRT) can be used to transform CoDa from 𝒮p to R. One of the common LRT for CoDa vector y∈𝒮p is the centered log-ratio transformation (CLRT),

clr(y)=(ln⁡y1y¯G,ln⁡y2y¯G,…,ln⁡ypy¯G),(5)

where y¯G is the component-wise geometric mean of y, i.e.,

y¯G=(∏i=1p⁡yi)1p=exp⁡(1p∑i=1pln⁡yi).(6)

Another LRT for a CoDa vector y∈𝒮p is the isometric log-ratio transformation (ILRT),

ilr(y)=y∗=clr(y)B⊺,(7)

where for B, we have many options (see [15]). B is a matrix of size (p−1,p). In this paper, B is taken as

Bi,j={1(p−i)(p−i+1)j≤p−i−p−ip−i+1j=p−i+10j>p−i+1.

Also, there is inverse ILRT, which can transform the ilr-coordinates y∗ into the composition coordinates y. The inverse ILRT is defined as

ilr−1(y∗)=y=𝒞(exp⁡(y∗B)).(8)

In this paper, the CoDa vectors are denoted without a “*” sign, while the transformed vectors are denoted with a “*” in the superscript.

3  Additive Measurement Errors Model for CoDa

Assume that there are p-part compositions yi=(yi,1,…,yi,p)∈𝒮p at time i=1,2,…, where y1,y2,… are independent multivariate normal (MNOR) random compositions with MV μ∗ and VCM Σ∗, i.e., yi∼MNOR𝒮p(μ∗,Σ∗). Suppose yi cannot be directly observable; instead, we use quality characteristics xi,1,…,xi,m to measure yi with xi,j, j=1,…,m, following the L.M with covariate M.E defined below:

xi,j=a⊕(b⊙yi)⊕εi,j,(9)

where a∈𝒮p is a perturbation and b∈R is a powering constant. With εi,j random error term follow a multivariate normal (MNOR) distribution MNOR𝒮p(0,ΣM∗), where ΣM∗ is the known M.E’s VCM. This M.E’s model is defined in [32]. Following [19], the researchers advocated averaging the m measurements to minimize the impact of M.E’s and keep the M.E’s component’s variance to a minimum. In this paper, following ([15]), the authors define the composition MV x¯i at the time i=1,2,…:

x¯i=1m⊙(xi,1⊕⋯⊕xi,m)=a⊕(b⊙yi)⊕(1m⊙(εi,1⊕⋯⊕εi,m)).

where a∗=ilr(a) and x¯i∼MNOR𝒮p(μx¯∗,Σx¯∗) with

μx¯∗=a∗+bμ∗,(10)

Σx¯∗=b2Σ∗+1mΣM∗.(11)

One can refer to Appendix A for the estimation of a and b.

4  Multivariate CUSUM CC for Compositional Data

Assume we have yt∼MNOR𝒮p(μ0∗,Σ∗) and yt∼MNOR𝒮p(μ1∗,Σ∗), where the IC and OOC process mean is defined by μ0∗ and μ1∗, respectively. The MCUSUM-CoDa CC corresponding to the LM with covariate M.E monitors the following statistic:

Ct=[(st⊺(Σx¯m)−1st)]1/2,t=1,2,…(12)

where st is the vector in Rp−1 defined as

st=0ifQt≤kst=(st−1+x¯t∗−a∗−bμ0∗)(1−k/Qt)ifQt>k,

with s0=0 and k>0. The value of k is chosen relative to the size of the shift involved, such that k=12δ, where δ is the size of the shift in standard deviation units. This approach comes very close to minimizing the ARL1 value for detecting a shift of size δ for fixed ARL0. It is necessary to note that a widely used value of k in practice is k=12. Let

Qt=[(st−1+x¯t∗−a∗−bμ0∗)⊺(Σx¯m)−1(st−1+x¯t∗−a∗−bμ0∗)]1/2.

The MCUSUM-CoDa CC displays a signal when Ct>h, where h is the upper control limit (UCL). A pre-defined ICARL0 is used to select the suitable value of h.

To assess the MCUSUM-CoDa CC’s run-length efficiency, the authors implement a Markov chain approximation suggested by Crosier [42]. The Markov Chain model needed to measure the MCUSUM-CoDa chart’s ARL is given in [43].

Lowry et al. [44] proposed ARL of MCUSUM charts depends on the parameter δ for non-centrality. Where δ can be written as

δ=(μ1−μ0)(Σ)−1(μ1−μ0)⊺.

When M.E’s is included in the process, the non-centrality parameter is denoted by δM and is equal to

δM=b2(μ1∗−μ0∗)(b2Σ∗+1mΣM∗)−1(μ1∗−μ0∗)⊺.(13)

Linna et al. [32] examined the effectiveness of the shift detection of multivariate charts and found them to be more effective in detecting a shift in one direction and ineffective in detecting a shift in the other direction. For this reason, Linna et al. [32] proposed using the minimum and maximum values of δM, i.e., δB (minimum value of δM) and δW (maximum value of δM). Using (7), we get

δB=δλ1,δW=δλp−1,

where λp−1 and λ1 are the largest and the smallest eigenvalues of the (p−1,p−1) matrix b2Σ∗(b2Σ∗+1mΣM∗)−1, respectively.

The ARL of the MCUSUM-CoDa chart can be obtained using the Markov chain technique suggested by Lowry et al. [44].

The authors assumed the unit sample size n=1; in that case, we can only detect the shift in the MV of the MCUSUM-CoDa chart. But if we use n>p, it is also possible to monitor variability using the statistic defined in [45] (see, for instance, [21]).

5  Performance Evaluation of the MCUSUM-CoDa CC in the Presence of M.E’s

This section aims to evaluate the performance of the MCUSUM-CoDa CC in the presence of M.E’s. The M.E’s VCM equals ΣM∗=σM2I, where σM is the standard-deviation M.E’s, and I is the (2,2) identity matrix. The optimal couple (k,H) that minimizes the OOCARL for a fixed value of the shift δM subject to a constrained value for the ICARL is chosen first in the design of the MCUSUM-CoDa CC. The procedure for determining the optimal couple consists of three steps:

•   Select the most appropriate value for the ICARL0. In this case, we choose ARL0=370.

•   Select a particular value of k. In this case, following [46], we choose k=0.5.

•   Obtain the value of H such that ARL=ARL0 when no shift occurs, i.e., when δ=δM=0.

•   Select the optimal value of H* from all the H values such that the OOCARL value is minimum for a given value of the shift δ for μ∗ (to attain the best statistical performance).

The following six cases for the CoDa VCM Σ∗ are considered:

Case #I uncorrelated with equal variances

Σ∗=(3003),

Case #II negatively correlated with equal variances

Σ∗=(3−1/2−1/23),

Case #III uncorrelated with unequal variances

Σ∗=(1.5003),

Case #IV positively correlated with unequal variances

Σ∗=(1.51/21/23),

Case #V negatively correlated with unequal variances

Σ∗=(1.5−1/2−1/23),

Case #VI positively correlated with equal variances

Σ∗=(31/21/23).

The goal is to investigate the effect of involved parameters σM, b, p and m to evaluate the performance of the MCUSUM-CoDa CC in the presence of M.E’s for the above-defined cases.

5.1 Case I

This subsection investigates the effect of the involved parameters σM, b and m corresponding to Case #1 (i.e., uncorrelated case with equal variances). For shifts δ∈[0,2], different ARL values are given in Table 1 for the MCUSUM-CoDa CC without and in the presence of M.E’s. The different ARL curves are presented in Fig. 1. From Table 1, we can conclude:

•   The M.E VCM σM has a negative impact on the ARL performance of the proposed CC. When we increase the values of σM, the OOCARL also indicates a slight increase. For example, when we select σM=0.1 and δ=1.5 the ARL values are (ARLB=6.95,ARLW=6.97). But, when σM=0.6 (i.e., σM increased), the ARL values are (ARLB=7.45,ARLW=7.57).

•   The subgroup size m has a positive impact on the ARL performance of the proposed CC. An increase in the value of m imposes a slight decrease in the ARL values. For example, for the fixed values of σM=0.3, b=1 and p=3, when m=1 and δ=1.5, (ARLB=7.76,ARLW=7.94). But, when we increase m=6, (ARLB=6.99,ARLW=7.02).

•   The powering constant b has a positive impact on the ARL performance of the proposed CC. An increase in the value of b imposes a significant decrease in the ARL values. For example, for the fixed values of σM=0.3, m=3 and p=3, when b=0.25 and δ=1.5, (ARLB=12.32,ARLW=13.57). But, when we increase b=8, (ARLB=ARLW=6.85).

Figure 1: ARL curves corresponding to the without measurement error case (solid line), Best ARL (dotted lines) and Worst ARL (dashed lines), for Case #1

5.2 Case II

This subsection investigates the effect of the involved parameters σM, b and m corresponding to Case #2 (i.e., negatively correlated case with equal variances). For shifts δ∈[0,2], different ARL values are given in Table 2 for the MCUSUM-CoDa CC without and in the presence of M.E’s. The different ARL curves are presented in Fig. 2. From Table 2, we can conclude:

•   The M.E VCM σM has a negative impact on the ARL performance of the proposed CC. When we increase the values of σM, the OOCARL also indicates an increase. For example, when we select σM=0.1 and δ=1.5 the ARL values are (ARLB=6.94,ARLW=6.98). But, when σM=0.6 (i.e., σM increased), the ARL values are (ARLB=7.39,ARLW=7.67).

•   The subgroup size m has a positive impact on the ARL performance of the proposed CC. An increase in the value of m imposes a slight decrease in the ARL values. For example, for the fixed values of σM=0.3, b=1 and p=3, when m=1 and δ=1.5, (ARLB=7.67,ARLW=8.11). But, when we increase m=6, (ARLB=6.98,ARLW=7.05).

•   The powering constant b has a positive impact on the ARL performance of the proposed CC. An increase in the value of b imposes a significant decrease in the ARL values. For example, for the fixed values of σM=0.3, m=3 and p=3, when b=0.25 and δ=1.5, (ARLB=11.75,ARLW=14.66). But, when we increase b=8, (ARLB=ARLW=6.85).

Figure 2: ARL curves corresponding to the without measurement error case (solid line), Best ARL (dotted lines) and Worst ARL (dashed lines), for Case #2

5.3 Case III

This subsection investigates the effect of the involved parameters σM, b and m corresponding to Case #1 (i.e., uncorrelated case with unequal variances). For shifts δ∈[0,2], different ARL values are given in Table 3 for the MCUSUM-CoDa CC without and in the presence of M.E’s. The different ARL curves are presented in Fig. 3. From Table 3, we can conclude:

•   The M.E VCM σM has a negative impact on the ARL performance of the proposed CC. When we increase the values of σM, the OOCARL also indicates an increase. For example, when we select σM=0.1 and δ=1.5 the ARL values are (ARLB=6.94,ARLW=7.04). But, when σM=0.6 (i.e., σM increased), the ARL values are (ARLB=7.44,ARLW=8.07).

•   The subgroup size m has a positive impact on the ARL performance of the proposed CC. An increase in the value of m imposes a slight decrease in the ARL values. For example, for the fixed values of σM=0.3, b=1 and p=3, when m=1 and δ=1.5, (ARLB=8.72,ARLW=7.75). But, when we increase m=6, (ARLB=7.14,ARLW=6.99).

•   The powering constant b has a positive impact on the ARL performance of the proposed CC. An increase in the value of b imposes a significant decrease in the ARL values. For example, for the fixed values of σM=0.3, m=3 and p=3, when b=0.25 and δ=1.5, (ARLB=18.93,ARLW=12.30). But, when we increase b=8, (ARLB=ARLW=6.85).

Figure 3: ARL curves corresponding to the without measurement error case (solid line), Best ARL (dotted lines) and Worst ARL (dashed lines), for Case #3

5.4 Case IV

This subsection investigates the effect of the involved parameters σM, b and m corresponding to Case #1 (i.e., positively correlated case with unequal variances). For shifts δ∈[0,2], different ARL values are given in Table 4 for the MCUSUM-CoDa CC without and in the presence of M.E’s. The different ARL curves are presented in Fig. 4. From Table 4, we can conclude:

•   The M.E VCM σM has a negative impact on the ARL performance of the proposed CC. When we increase the values of σM, the OOCARL also indicates an increase. For example, when we select σM=0.1 and δ=1.5 the ARL values are (ARLB=7.06,ARLW=6.94). But, when σM=0.6 (i.e., σM increased), the ARL values are (ARLB=8.21,ARLW=7.41).

•   The subgroup size m has a positive impact on the ARL performance of the proposed CC. An increase in the value of m imposes a slight decrease in the ARL values. For example, for the fixed values of σM=0.3, b=1 and p=3, when m=1 and δ=1.5, (ARLB=8.94,ARLW=7.71). But, when we increase m=6, (ARLB=7.17,ARLW=6.85).

•   The powering constant b has a positive impact on the ARL performance of the proposed CC. An increase in the value of b imposes a significant decrease in the ARL values. For example, for the fixed values of σM=0.3, m=3 and p=3, when b=0.25 and δ=1.5, (ARLB=20.53,ARLW=12.01). But, when we increase b=8, (ARLB=ARLW=6.85).

Figure 4: ARL curves corresponding to the without measurement error case (solid line), Best ARL (dotted lines) and Worst ARL (dashed lines), for Case #4

5.5 Case V

This subsection investigates the effect of the involved parameters σM, b, p and m corresponding to Case #5 (i.e., negatively correlated with unequal variances). For shifts δ∈[0,2], different ARL values are given in Table 5 for the MCUSUM-CoDa CC without and in the presence of M.E’s. The different ARL curves are presented in Fig. 5. From Table 5, we can conclude:

•   The ARL in case #5 are greater than the above four cases for all the parameters σM, m and b.

•   The M.E VCM σM has a negative impact on the ARL performance of the proposed CC. When we increase the values of σM, the OOCARL also indicates an increase. For example, when we select σM=0.1 and δ=1.5 the ARL values are (ARLB=7.09,ARLW=7.23). But, when σM=0.6 (i.e., σM increased), the ARL values are (ARLB=7.60,ARLW=8.46).

•   The subgroup size m has a positive impact on the ARL performance of the proposed CC. An increase in the value of m imposes a slight decrease in the ARL values. For example, for the fixed values of σM=0.3, b=1 and p=3, when m=1 and δ=1.5, (ARLB=7.92,ARLW=9.25). But, when we increase m=6, (ARLB=7.14,ARLW=7.35).

•   The powering constant b has a positive impact on the ARL performance of the proposed CC. An increase in the value of b imposes a significant decrease in the ARL values. For example, for the fixed values of σM=0.3, m=3 and p=3, when b=0.25 and δ=1.5, (ARLB=12.64,ARLW=22.57). But, when we increase b=8, (ARLB=7.00,ARLW=7.01).

Figure 5: ARL curves corresponding to the without measurement error case (solid line), Best ARL (dotted lines) and Worst ARL (dashed lines), for Case #5

5.6 Case VI

This subsection investigates the effect of the involved parameters σM, b, p and m corresponding to Case #6 (i.e., positively correlated with equal variances). For shifts δ∈[0,2], different ARL values are given in Table 6 for the MCUSUM-CoDa CC without and in the presence of M.E’s. The different ARL curves are presented in Fig. 6. From Table 6, we can conclude:

•   Similar to Case #5, the ARL in case #6 are also greater than the first four cases for all the parameters σM, m and b.

•   The M.E VCM σM has a negative impact on the ARL performance of the proposed CC. When we increase the values of σM, the OOCARL also indicates an increase. For example, when we select σM=0.1 and δ=1.5 the ARL values are (ARLB=7.08,ARLW=7.12). But, when σM=0.6 (i.e., σM increased), the ARL values are (ARLB=7.54,ARLW=7.76).

•   The subgroup size m has a positive impact on the ARL performance of the proposed CC. An increase in the value of m imposes a slight decrease in the ARL values. For example, for the fixed values of σM=0.3, b=1 and p=3, when m=1 and δ=1.5, (ARLB=7.82,ARLW=8.17). But, when we increase m=6, (ARLB=7.13,ARLW=7.18).

•   The powering constant b has a positive impact on the ARL performance of the proposed CC. An increase in the value of b imposes a significant decrease in the ARL values. For example, for the fixed values of σM=0.3, m=3 and p=3, when b=0.25 and δ=1.5, (ARLB=12.00,ARLW=14.38). But, when we increase b=8, (ARLB=ARLW=7.00).

Figure 6: ARL curves corresponding to the without measurement error case (solid line), Best ARL (dotted lines) and Worst ARL (dashed lines), for Case #6

5.7 Effect of Number of Variables p

This subsection investigates the effect of the number of variables involved p for all 6 cases mentioned above. The values of Σ∗ corresponding to the Cases #1 to #6 are given in Table 7, when p=3, p=5 and p=10 are considered.

For shifts δ∈[0,2], different ARL values are given in Table 8 for the MCUSUM-CoDa CC in the presence of M.E's. The ARL curves for different values of shift for all the six cases are also presented in Figs. 7 and 8. From Table 8, we can conclude:

•   For Case #1, when we increase the values of p, the OOCARL also indicates an increase. For example, when we select p=3 and δ=1.5 the ARL values are (ARLB=7.14,ARLW=7.2). But, when p=10 (i.e., p increased), the ARL values are (ARLB=ARLW=10.40).

•   For Case #2, when we increase the values of p, the OOCARL also indicates an increase. For example, when we select p=3 and δ=1.5 the ARL values are (ARLB=7.12,ARLW=7.25). But, when p=5 (i.e., p increased), the ARL values are (ARLB=8.19,ARLW=10.32).

•   For Case #3, when we increase the values of p, the OOCARL also indicates an increase. For example, when we select p=3 and δ=1.5 the ARL values are (ARLB=7.44,ARLW=7.14). But, when p=5 (i.e., p increased), the ARL values are (ARLB=11.23,ARLW=11.6).

•   For Case #4, when we increase the values of p, the OOCARL also indicates an increase. For example, when we select p=3 and δ=1.5 the ARL values are (ARLB=7.51,ARLW=7.12). But, when p=5 (i.e., p increased), the ARL values are (ARLB=11.97,ARLW=13.04).

•   For Case #5, when we increase the values of p, the OOCARL also indicates an increase. For example, when we select p=3 and δ=1.5 the ARL values are (ARLB=7.29,ARLW=7.71). But, when p=5 (i.e., p increased), the ARL values are (ARLB=9.03,ARLW=10.93).

•   For Case #6, when we increase the values of p, the OOCARL also indicates an increase. For example, when we select p=3 and δ=1.5 the ARL values are (ARLB=7.26,ARLW=7.37). But, when p=5 (i.e., p increased), the ARL values are (ARLB=10.07,ARLW=10.52).

Figure 7: ARL curves corresponding to the without measurement error case (solid line), Best ARL (dotted lines) and Worst ARL (dashed lines), for Cases #1 to #3

Figure 8: ARL curves corresponding to the without measurement error case (solid line), Best ARL (dotted lines) and Worst ARL (dashed lines), for Cases #4 to #6

6  Practical Implementation of the Proposal in the Manufacturing Scenario

For the proposal implementation, two illustrated examples have been incorporated in this study: one involving manufacturing uncoated aspirin tablets with the composition of aspirin, micro-crystalline cellulose, talcum powder and magnesium stearate and the other is the monitoring of machines that are responsible for the percentages of three components (whole-grain cereal, dried fruits, and nuts) in muesli production.

6.1 Proposal Implementation Related to Aspirin Manufacturing Process

Pharmaceutical manufacturing is the process of synthesizing pharmaceutical medications on a large scale as part of the pharmaceutical business. The medication production process may be divided into several unit activities, such as milling, granulation, coating, tablet pressing, etc. Here the authors simulated an example of a company producing uncoated aspirin tablets. The base parameters have been taken from [47]. According to Tiwari et al. [47], the composition of the tablet consists of 150 mg of aspirin, 60 mg of micro crystalline cellulose, 20 mg of talcum powder and magnesium stearate. The company wants to check the measuring unit in charge of measuring raw material (i.e., estimate a, b and ΣM). For this purpose, the firm cautiously prepared k=4 samples of aspirin and measured them n=10 times each. The data on manufactured aspirin and the ARL values are given in Table 9 and Fig. 9. By using the data in Table 9, the estimated values of a^, b^ and Σ^M have been obtained. The results of estimation are as follows:

•   a^=(0.339,0.327,0.334),

•   b^=1.074,

•   Σ^M={0.0010.0010.0010.009}.

Figure 9: Data used to check machine: yi and xi,j in 𝒮p (on left), yi∗ and xi,j∗ in R2 (on right)

The estimated values of μ0^ and Σ0^ of control parameters, 20 sub batches of aspirin were sampled and are measured thrice each. The simulation results are in Table 10 and Fig. 10. Using the data in Table 10, one can easily obtain:

•   μx¯^=(1.376,0.719)

•   Σ^x¯=(0.00090.00240.00240.0219).

Figure 10: Phase-I data for aspirin example: x¯i (on left) and x¯i∗ (on right)

Using the above values, we can obtain μ0^ and Σ^.

•   μ0^=(1.284,0.646)

•   Σ^={0.00050.00180.00180.0164}.

The results of the MCUSUM-CoDa charts are also presented in Table 10. Fig. 11 also shows the values of the MCUSUM-CoDa chart along with UCL. In Phase-II of the production process, 20 batches of aspirin have been produced and sampled three times each. These values are shown in Table 11, and Fig. 12 also shows the MCUSUM-CoDa chart and UCL values. It can be seen from Fig. 13 that the process was IC till sample #3, but sample #4 went OOC because of an error. Hence the process stops at sample #4 to study the reason for the OOC point and restarts again after fixing the error, and it is shown in Fig. 12 that after fixing the error, the process remains IC throughout all the samples.

Figure 11: MCUSUM-CoDa chart for aspirin Phase-I data

Figure 12: Phase-II data for aspirin example: x¯i (on left) and x¯i∗ (on right)

Figure 13: MCUSUM-CoDa chart for aspirin Phase-II data

6.2 Proposal Implementation Related to the Muesli Manufacturing Process

For the second example, the authors are using the same example used in [19] and [21] of a company that manufactures breakfast muesli with (A) 66 percent of cereals, (B) 24 percent dried fruits, and (C) 10 percent nuts in every 100 grams. The company has decided to begin by calibrating the measurement device used to determine the percentage of each component in the manufactured muesli (i.e., estimate a, b and ΣM∗). The company meticulously prepared k=4 samples with perfectly-known percentages for each component and measured them m=7 times. Table 12 shows the results obtained and the ilr transformed values to perform this calibration. yi∗=ilr(yi) and xi,j∗=ilr(xi,j).

The following estimates are obtained using Table 12 with d=2, vi=yi∗ and ui,j=xi,j∗ for i=1,…,4 and j=1,…,7, we get the following estimators a^, b^ and Σ^M∗ for a, b and ΣM∗:

•   a^=(0.3354,0.3357,0.3289) or, equivalently, a^∗=(0.0163,−0.0006),

•   b^=1.1070,

•   Σ^M∗={0.00140.00080.00080.0103}.

To estimate the IC composition parameters (Phase-I) μ^0∗ and Σ^∗, 20 batches of muesli have been measured m=3 times. Table 13 shows the results with the values of xi,j, x¯i and x¯i∗, i=1,…,20. The results for x¯i and x¯i∗ are plotted in Fig. 14. From the values x¯i∗, it is easy to obtain:

•   μ^x¯∗=(1.2766,0.7657),

•   Σ^x¯∗={0.01460.01060.01060.0511}.

Figure 14: Phase-I data for the muesli example: x¯i (on left) and x¯i∗ (on right)

Then, using (10) and (11), we have

•   μ^0∗=1b^(μ^x¯∗−a^∗)=(1.1385,0.6922),

•   Σ^∗=1b^2(Σ^x¯∗−1mΣ^M∗)=(0.01160.00840.00840.0389).

The authors have computed the MCUSUM-CoDa statistics Ci∗ for i=1,…,20 using Eq. (12). The UCL for δ=1.5 is H=9.715. The MCUSUM-CoDa statistics Ci∗ values are listed in Table 13 and plotted in Fig. 15 with the UCL=H=9.715. One can see that the process is IC, as all the values in Fig. 15 are smaller than the UCL.

Figure 15: MCUSUM-CoDa CC for muesli Phase-I data

Concerning Phase-II, i=1,…,20 batches of muesli have also been measured m=3 times. Table 14 shows the results with the values xi,j, x¯i and x¯i∗, i=1,…,20. The results x¯i and x¯i∗ are also plotted in Fig. 16. The MCUSUM-CoDa statistics Ci∗ values in the presence of M.E’s are listed in Table 14 and plotted in Fig. 17 with the UCL=H=9.715 obtained in Phase-I.

Figure 16: Phase-II data for the muesli example: x¯i (on left) and x¯i∗ (on right)

Figure 17: MCUSUM-CoDa CC for muesli Phase-II data

The process seems to be IC up to sample #14, but sample #15 and sample #16 are OOC. Then, in samples #15, #16 and #17, it is discovered that a hatch malfunction had caused a sudden drop in the quantity of whole-grain cereals, causing a shift from μ^0∗=(1.1385,0.6922) to μ^x¯∗=(1.2766,0.7657). Using Eq. (13), we found that δ=1.62. Here, a shift of size δ=1.5 in μ^x¯∗ has been interpreted. For m=3 and δ=1.5, the UCL of MCUSUM-CoDa chart is H=9.715. The process restarts without incident after the hatch has been repaired.

7  Conclusions

The CCs are an effective statistical process monitoring tool frequently used to assess the stability of manufacturing processes. They monitor the industrial process by giving practitioners an early signal about an OOC process to keep the process IC, but M.E negatively impacts the performance of CCs in quality control applications. In this paper, the authors investigate the impact of M.E on the MCUSUM-CoDa CC’s performance. The identity matrix assumption has been incorporated into the diagonal matrix in the LM with covariate M.E. This makes it easier to see how the M.E effects the proposed chart’s performance. Six cases have been considered for the VCM of CoDa. Different combinations of the involved parameters were investigated to inspect the effect of the parameters on the CC. The following are some key findings of the study from this research. For all six cases, (i). when we increase the values of σM, the OOCARL also indicates an increase (ii). an increase in the value of m imposes a slight decrease in the ARL values (iii). an increase in the value of b imposes a significant decrease in the ARL values (iv). an increase in the values of p imposes a significant increase in the ARL. All these findings are relevant to the existing studies in the literature. In the end, two illustrated examples are given to demonstrate the implementation of the proposed CCs in the presence of M.E. One is based on the manufacturing process of uncoated aspirin tablets, and the other is based on monitoring the machines used for manufacturing the muesli. For future works, the proposed chart can be studied with a variable sampling interval or variable sample size scheme. Also, the effect of estimated parameters can be seen on the previously defined CCs for CoDa.

Data and Code Availability Statement: All data, models, and codes supporting the results of this study are accessible on reasonable request from the corresponding author.

Ethics Approval: This manuscript has not been simultaneously submitted to more than one journal. So the proposed work is original and has not been published anywhere else in any form or language.

Consent to Participate and Publish: Consent was obtained from all authors for their participation in completing this manuscript and their consent to publish it.

Author’s Contribution: All the authors contributed to this paper.

Funding Statement: This work is supported by the National Natural Science Foundation of China (Grant No. 71802110); the Humanity and Social Science Foundation of the Ministry of Education of China (Grant No. 19YJA630061).

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

Appendix A

Let u1,…,uk be k row vectors in Rp. For each vector ui, i=1,…,k, we observe n random vectors vi,1,…,vi,n and we assume that vi,j=a+bui+εi,j, where a∈Rp is a constant row vector and b∈R is a constant scalar, while a random error term following multivariate normal distribution is defined as εi,j with MNORRd(0,Σ). We define the (p×k×n,1) column vector v and the (p×k×n,p+1) matrix U for the purpose of estimation of unknown parameters a and b as

v=(v1,1⊺⋮v1,n⊺v2,1⊺⋮v2,n⊺⋮vk,1⊺⋮vk,n⊺)andU=(Iu1⊺⋮⋮Iu1⊺Iu2⊺⋮⋮Iu2⊺⋮⋮Iuk⊺⋮⋮Iuk⊺),

where I is the (p, p) identity matrix. Ordinary least square method has been used to compute the (d+1,1) column vector c^ for minimizing the sum of square of residuals Σ that is,

c^=(V⊺V)−1V⊺u.

Then, an estimation of a will be the first p components of c^=(c1,c2,…,cp,cp+1)⊺, while estimation of b will be the last component (i.e., p+1). We get the estimation of Σ as

Σ^=1k×n∑i=1k∑j=1n(vi,j−a^−b^ui)⊺(vi,j−a^−b^ui).

Appendix B

Let D be a symmetric matrix of size (p, p) and E>0 be a positive definite matrix of same size. Suppose ymax is the row vector of size (1,p), satisfying

ymax=argminy(yGy⊺),

subject to the constraint

ymaxEymax⊺=a>0.

As presented in [48] (Theorem A.9.2 page no. 479), one can prove that if λ1≤λ2≤⋯≤λp and u1,u2,…,up will be eigenvalues in ordered manner and eigenvectors E−1D, respectively, then

ymax=u1cu1Eu1⊺,

and

maxy(yDy⊺)=ymaxGymax⊺=aλ1.

Similarly

maxx(xGx⊺)=xmaxGxmax⊺=cλp.

Following Mardia 1979, the author also used c = 1.

Appendix C

The Scilab Codes for ARL, when p=3, b=1, m=3, σM=0.3 for Case #1 (i.e., Uncorrelated with equal variances)

ARL0 = 370;

d1 = 0.25; d2 = 0.5; d3 = 0.75; d4 = 1; d5 = 1.25; d6 = 1.5; d7 = 1.75; d8 = 2;

p = 3;

b = 1;

mm = 3;

S = [3, 0;0, 3];

SM = 0.3 * eye(p−1, p−1);

b2 = b * b

AA = b2 * inv(b2 * S + SM/mm)

BB = inv(S)

invBB = S [eve, eva] = spec(invBB * AA) if eva(1, 1) < eva(2, 2) then lam1 = eva(1, 1) lam2 = eva(2, 2) v1 = eve(:, 1)

vp = eve(:, 2)

else

lam1 = eva(2, 2)

lam2 = eva(1, 1)

v1 = eve(:, 2)

vp =eve(:, 1)

end

m1 = 30

m2 = 30

dmi1 = d1 * lam1 dma1 = d1 * lam2

dmi2 = d2 * lam1 dma2 = d2 * lam2

dmi3 = d3 * lam1 dma3 = d3 * lam2

dmi4 = d4 * lam1 dma4 = d4 * lam2

dmi5 = d5 * lam1 dma5 = d5 * lam2

dmi6 = d6 * lam1 dma6 = d6 * lam2

dmi7 = d7 * lam1 dma7 = d7 * lam2

dmi8 = d8 * lam1 dma8 = d8 * lam2

[r21,H21,ARL21,SDRL21] = optMEWMAUTM(m1, m2, p−1,ARL0, dmi1)

Z21 = [dmi1;r21;H21;ARL21] [r31,H31,ARL31,SDRL31] = optMEWMAUTM(m1, m2,

p−1,ARL0, dma1)

Z31 = [dma1;r31;H31;ARL31]

disp(Z21) disp(Z31)

[r22,H22,ARL22,SDRL22] = optMEWMAUTM(m1, m2, p−1,ARL0, dmi2)

Z22 = [dmi2;r22;H22;ARL22] [r32,H32,ARL32,SDRL32] = optMEWMAUTM(m1, m2,

p−1,ARL0, dma2)

Z32 = [dma2;r32;H32;ARL32]

disp(Z22) disp(Z32) [r23,H23,ARL23,SDRL23] = optMEWMAUTM(m1, m2, p−1,ARL0, dmi3)

Z23 = [dmi3;r23;H23;ARL23] [r33,H33,ARL33,SDRL33] = optMEWMAUTM(m1, m2,p−1,ARL0, dma3)

Z33 = [dma3;r33;H33;ARL33]

disp(Z23) disp(Z33)

[r24,H24,ARL24,SDRL24] = optMEWMAUTM(m1, m2, p−1,ARL0, dmi4)

Z24 = [dmi4;r24;H24;ARL24] [r34,H34,ARL34,SDRL34] = optMEWMAUTM(m1, m2,

p−1,ARL0, dma4)

Z34 = [dma4;r34;H34;ARL34]

disp(Z24) disp(Z34)

[r25,H25,ARL25,SDRL25] = optMEWMAUTM(m1, m2, p−1,ARL0, dmi5)

Z25 = [dmi5;r25;H25;ARL25] [r35,H35,ARL35,SDRL35] = optMEWMAUTM(m1, m2,

p−1,ARL0, dma5)

Z35 = [dma5;r35;H35;ARL35]

disp(Z25) disp(Z35)

[r26,H26,ARL26,SDRL26] = optMEWMAUT(m1, m2, p−1,ARL0, dmi6)

Z26 = [dmi6;r26;H26;ARL26] [r36,H36,ARL36,SDRL36] = optMEWMAUT(m1, m2, p−1,ARL0, dma6)

Z36 = [dma6;r36;H36;ARL36]

disp(Z26) disp(Z36)

[r27,H27,ARL27,SDRL27] = optMEWMAUTM(m1, m2, p−1,ARL0, dmi7)

Z27 = [dmi7;r27;H27;ARL27] [r37,H37,ARL37,SDRL37] = optMEWMAUTM(m1, m2,

p−1,ARL0, dma7)

Z37 = [dma7;r37;H37;ARL37]

disp(Z27) disp(Z37)

[r28,H28,ARL28,SDRL28] = optMEWMAUTM(m1, m2, p−1,ARL0, dmi8)

Z28 = [dmi8;r28;H28;ARL28] [r38,H38,ARL38,SDRL38] = optMEWMAUTM(m1, m2,

p−1,ARL0, dma8)

Z38 = [dma8;r38;H38;ARL38]

disp(Z28) disp(Z38)

Where the base function “optCUSUMUTM” to find the ARL is given below,

//———————————— function [Q0, q0] = QCUSUM0M(H, k, p) //————————–

// in control [argout, argin] = argn() if (argin<3)|(argin>4) error(“incorrect number of arguments”) end if H<=0 error(“argument “H” must be > 0”)

end if k<=0 error(“argument “k” must be > 0”)

end if p<=0 error(“argument “p” must be > 0”)

end

UCL = H

m = 100

g = 2 * UCL/(2 * m + 1)

Q0 = ones(m + 1, m + 1)

t = [0:1:m]

c=((1−k)∗t∗g/k)2x0=(0.5∗g/k)2Q01=cdfchn(‘‘PQ″,x0∗ones(c),p∗ones(c),c)

tt=[1:1:m]′upi=((tt+0.5)∗g/k)2up=repmat(upi,1,m+1)downi=((tt−0.5)∗g/k)2

down = repmat(downi, 1, m + 1)

C = repmat(c, m, 1)

Q02 = cdfchn(“PQ”, up, p * ones(m, m + 1),C)-cdfchn(“PQ”, down, p * ones(m, m + 1),C)

Q0 = [Q01;Q02]’

q0 = zeros(m + 1, 1)

q0(1) = 1

endfunction

//————————— function [ARL,SDRL] = rlCUSUM0M(H, k, p) //—————————

// in-control ARL [argout, argin] = argn() if (argin<3)|(argin>4) error(”incorrect number of arguments”) end if H<=0 error(“argument “H” must be > 0”)

end if k<=0 error(“argument “k” must be > 0”)

end if p<=0 error(“argument “p” must be > 0”)

end [Q0, q0] = QCUSUM0(H, k, p)

W = inv(eye(Q0)-Q0)

Z = q0’ * W

ARL = sum(Z)

if argout == 2

Z = Z * W * Q0

SDRL=sqrt(2∗sum(Z)−ARL2+ARL) end endfunction

//—————————- function [Q, q] = QCUSUM1M(H, m1, m2, k, p, delta) //——————

//Out-of-control [argout, argin] = argn() if (argin<5)|(argin>6) error(“incorrect number of arguments”) end if H<=0 error(“argument “H” must be > 0”)

end if k<=0 error(“argument “k” must be > 0”)

end if p<=0 error(“argument “p” must be > 0”)

end

UCL = H;

g1 = 2 * UCL/(2 * m1 + 1);

Q1 = ones(2 * m1 + 1, 2 * m1 + 1);

t1 = [1:1:2 * m1 + 1];

c1 = -UCL + (t1–0.5) * g1;

C1 = repmat(c1, 2 * m1 + 1, 1);

tt1 = [1:1:2 * m1 + 1]’;

Tt1 = repmat(tt1, 1, 2 * m1 + 1);

up1 = (-UCL + Tt1 * g1−(1−k) * C1)/k−delta;

down1 = (−UCL + (Tt1–1) * g1−(1−k) * C1)/k-delta;

Q11 = cdfnor(“PQ”, up1, 0 * ones(C1), 1 * ones(C1))-cdfnor(“PQ”, down1, 0 * ones(C1), 1 * ones(C1)); Q1 = Q11’;

g2 = 2 * UCL/(2 * m2 + 1);

Q2 = ones(m2 + 1, m2 + 1);

t2 = [0:1:m2];

c2=((1−k)∗t2∗g2/k)2;x0=(0.5∗g2/k)2;Q21=cdfchn(“PQ”,x0∗ones(t2),(p−1)∗ones(t2),c2);tt2=[1:1:m2]′;

up2i=((tt2+0.5)∗g2/k)2;up2=repmat(up2i,1,m2+1);down2i=((tt2−0.5)∗g2/k)2;down2=repmat(down2i,1,m2+1);

C2=repmat(c2,m2,1);Q22=cdfchn(“PQ”,up2,(p−1)∗ones(m2,m2+1),C2)−cdfchn(“PQ”,down2,(p−1)∗ones(m2,m2+1),C2);

Q2=[Q21;Q22]′;Q=Q1.∗.Q2counter=1;foral=1:2∗m1+1forbe=0:m2a=(al−(m1+1))2∗g12+(be∗g2)2;ifa>=UCL2

Q(:, counter) = 0;

Q(counter,:) = 0;

end

counter = counter + 1;

end

end

q = zeros((2 * m1 + 1) * (m2 + 1), 1)

start = m1 * (m2 + 1) + 1

q(start) = 1

endfunction

//——————— function [ARL,SDRL] = rlCUSUM1M(H, m1, m2, k, p, delta) //—————

//Out-of-Control ARL [argout, argin] = argn() if (argin<5)|(argin>6) error(“incorrect number of arguments”) end if H<=0 error(“argument “H” must be > 0”)

end if k<=0 error(“argument “k” must be > 0”)

end if p<=0 error(“argument “p” must be > 0”)

end [Q, q] = QCUSUM1M(H, m1, m2, k, p, delta)

W = inv(eye(Q)-Q)

Z = q’ * W

ARL = sum(Z)

if argout == 2

Z = Z * W * Q

SDRL=sqrt(2∗sum(Z)−ARL2+ARL)

end endfunction

//——————————- function f = HCUSUMM(H, k, arl0, p) //——————————-

// function of H

if (H<=0)

f = else

ARL = rlCUSUM0M(H, k, p)

f = (arl0-ARL)

end endfunction

//————- function [H,ARL,SDRL] = optMEWMAUTM(m1, m2, p, arl0, delta) //————-

// Final ARL with optimal H [argout, argin] = argn()

if argin = 5 error(”incorrect number of arguments”)

end

lines(0)

k = 0.5

H0 = 7

H = simplexolve(H0,HCUSUMM, list(k, arl0, p), tol = 1e-4); [ARL,SDRL] = rlCUSUM1M(H, m1, m2, k, p, delta)

endfunction