Robust Hotelling 𝑇² Control Chart with Consistent Minimum Vector Variance

Abstract

Recently, an alternative robust control chart based on a new robust estimator known as minimum vector variance (MVV) estimator, 𝑇²_MVV, was introduced in Phase II. 𝑇²_MVV was able to detect out-of-control signal and simultaneously control false alarm rate even as the dimension increased. However, the estimated UCLs of 𝑇²_MVV are large as compared to the traditional chart. In this study, we improved the MVV estimators in terms of consistency and bias. The result showed great improvement in the control limit values while maintaining its good performance in terms of false alarm and probability of detection.

1. Introduction

Hotelling 𝑇² statistic was the first statistic known to be used in multivariate control chart. The control chart is referred to as Hotelling 𝑇² control chart. This statistic is used to measure the significance of the shifted distance from the out-of-control mean vector, 𝜇_𝑠, to the nominal mean vector, 𝜇₀, with the assumption that the covariance matrix remains constant at Σ₀. The purpose of the control chart is to monitor the stability of a multivariate process in Phase I and II. Analysis in Phase I seeks to identify a stable historical data set (HDS). From this dataset, the in-control mean vector and the in-control variance-covariance matrix are estimated and later will be used in the Phase II analysis. A successful process monitoring in Phase II totally depends on the estimates of the parameters obtained from a stable HDS. However, the estimators are easily affected by unstable process, that is, multivariate outliers. The existence of outliers can violate the normality assumption. This violation may lead to the inflation of control limits and reduction of the probability of detection in Phase I, which consequently will cause the level of false alarm to be distorted, and the power to detect changes will be reduced in Phase II process [1]. False alarm rate is the probability of out-of-control signal when a process is in control. The value becomes large if the process is unstable due to the increase in variability. Inflated false alarm rate can lead to unnecessary process adjustments and loss of confidence in the control chart as a monitoring tool [2]. Hence, a method which can control the false alarm rate to the desired (nominal) level is necessary.

However, the traditional Hotelling 𝑇² control chart is only effective in eliminating extreme outliers in small sample sizes, but it fails to detect moderate outliers particularly when the number of variables is large [3–5]. To overcome the problem, alternative estimation methods have been proposed in the literature. One of the approaches is to calculate the 𝑇² statistic based on successive differences variance-covariance matrix estimator [6–9]. Though this approach is effective in detecting shifts in the mean vector, it fails to detect other outliers as shown in Vargas [3]. Another approach is to use robust estimators in place of the classical estimators (𝑥 and 𝑆). Robust estimators are known to be more effective in detecting the deviation of data, or outliers as compared to the classical estimators [10]. A wide range of robust estimators of multivariate location and scatter is available; see [11, 12] for a review. However, MCD estimator is more attractive than others because it has good theoretical properties with affine equivariance, high breakdown value, bounded influence function, and better convergence rate [13, 14]. The study on the significant role of MCD estimators in the construction of robust Hotelling 𝑇² chart can be easily found in the literature. Vargas [3] and Jensen et al. [4] introduced robust control chart based on MCD estimator for multivariate individual observations. They identified and removed the outliers in Phase I analysis by using robust estimator and then calculated the classical estimator using the remaining good data points for Phase II analysis. They noticed some drawbacks when MCD was used in Phase I. Hotelling’s 𝑇² issued from MCD needed a larger sample size if large number of outliers was suspected to ensure that MCD estimator did not breakdown and lose its ability especially when monitoring with more quality characteristics (𝑝). To abate the problems, Chenouri et al. [5] proposed robust Hotelling 𝑇² chart based on reweighted MCD estimator. Besides possessing the nice properties of MCD estimator, the estimator was not unduly influenced by outliers and was more efficient than MCD. However, their approach was different from Vargas [3] and Jensen et al. [4] whereby they used RMCD estimator in place of classical estimators in constructing Hotelling 𝑇² chart for Phase II data directly. Using the same approach as Chenouri et al. [5], Alfaro and Ortega [15] conducted a comparison study on the performance of Hotelling 𝑇² control chart in Phase II process based on MCD, MVE, reweighted MCD, and trimmed estimator. Their finding showed a conflict between the percentage of outliers detection and the ability of the robust control charts in controlling the overall false alarm rate under certain conditions. To alleviate this conflict, Yahaya et al. [16] introduced the MVV estimator in Hotelling 𝑇² chart (𝑇²_MVV) in Phase II. In general, the result showed that 𝑇²_MVV chart was able to increase the detection of out-of-control signals and simultaneously control false alarm rates even with large number of quality characteristics. In contrast, the MCD charts performed well in detecting out-of-control signals but failed in controlling false alarm rates. The traditional chart, however, was able to control false alarm rates but not effective in detecting out-of-control signals. Despite the good performance of 𝑇²_MVV, the estimated UCLs for Hotelling 𝑇² chart issued from MVV estimators were large as compared to the traditional and MCD charts. Thus, this study attempts to improve the MVV estimators in achieving the desired UCLs by making the estimators consistent at normal model. Since in practice we always deal with finite samples, therefore the issue of bias in a finite sample is also considered in this study. The advantage of having unbiased estimator for a finite sample is that this estimator remains unbiased even though the sample size becomes larger [17]. With respect to the latter issue, this paper will also seek to improve the performance of MVV by making it unbiased for finite samples.

The organization of the remaining part of this paper is as follows. The formal definition of MVV estimator and the adjustment done on the MVV scatter estimator to ensure that it is consistent and unbiased will be discussed next, followed by the investigation on the improved MVV estimator through simulation study. The discussion continues with the computation of control limits for the traditional, MCD, and the improved Hotelling 𝑇²_MVV charts, and the improvement of the proposed chart is revealed in this section. A real data analysis from aircraft industry is presented to illustrate the applicability of the proposed charts before arriving to the conclusion in the last section.

2. Minimum Vector Variance (MVV) Estimator

Herwindiati et al. [18] had proved that MVV estimators possess three major properties of a good robust estimator, that is, high breakdown point, affine equivariance, and computational efficiency. The main method used in the estimation of MVV is the Mahalanobis squared distances (MSDs). Let 𝑋={𝑥₁,𝑥₂,…,𝑥_𝑛} be a data set of 𝑝-variate observations. Denote the MVV estimators for the location parameter and scatter by 𝑀_MVV, and 𝑆_MVV respectively. Now let 𝐻⊆𝑋, the 𝑚_MVV and 𝑆_MVV are determined based on the set 𝐻 consisting of ℎ=⌊(𝑛+𝑝+1)/2⌋ data such that 𝑆_MVV has minimum trace of 𝑆²_MVV, denoted as Tr(𝑆²_MVV), among all possible sets of ℎ data. To compute the estimates of MVV, we used the MVV algorithm proposed in Yahaya et al. [16]. The location and scatter estimators are defined as 𝑚_MVV=1ℎ^ℎ_𝑖=1𝑥_𝑖,𝑆_MVV=1ℎ^ℎ_𝑖=1𝑥_𝑖−𝑀_MVV𝑥_𝑖−𝑀_MVV^𝑡. (1)(2)

2.1. Consistency Factor

The aim of Hotelling 𝑇² chart in Phase I is to estimate the in-control parameters of location, 𝜇 and scatter, Σ. The usual estimators for these parameters are the normal maximum likelihood estimators (MLE). The estimation of parameters is based on the data set 𝑥={𝑥₁,𝑥₂,…,𝑥_𝑛} from multivariate normal distribution with density 𝑓_𝜇,Σ(𝑥)=1(2𝜋)^𝑝/2|Σ|^1/2𝑒^{(−(𝑥−𝜇)𝑡Σ−1(𝑥−𝜇)/2)}, (3) with 𝜇∈ℝ^𝑝 and Σ∈ℤ⁺. However, the distribution of (3) is only an approximation because a portion of the data may be contaminated by outliers [19]. With the existence of outliers, MLE which are known to be sensitive to outliers will not be able to precisely estimate the parameters. To address this problem, we propose MVV estimators, that is, robust estimators with highest breakdown point (50%) proposed by Herwindianti [20] to replace the MLE. We compute the MVV estimators in Phase I data sets, with location and scatter estimators as defined in (1) and (2), respectively. The MVV estimator has a fixed integer ℎ such that 𝑛+𝑝+12≤ℎ<𝑛. (4)

The preferred choice of ℎ for outlier detection is its lower bound, which yields the breakdown value, BP=(𝑛−2(𝑝−1))/2𝑛. Let 𝑚_MVV and 𝑆_MVV be the mean and the scatter matrix calculated from the ℎ observations out of 𝑥_𝑖, whose classical scatter matrix has the lowest vector variance resulting from ℎ smallest MSD. The 𝑆_MVV is a scatter 𝑝×𝑝 matrix which is positive definite, symmetric (PDS), and affine equivariant [20]. However, this estimator is not consistent under normal model. Robust scatter estimator is typically calibrated to be consistent for normal model. Known as Fisher consistency, this is a standard concept in robust statistics which denotes that the functionals evaluated at the model distribution return the true parameter value, Σ [19]. In order to achieve consistency under the normal model, 𝑆_MVV (in (2)) is multiplied by a consistency factor, 𝑐(ℎ), as follows: 𝑐(ℎ)𝑆_MVV=𝑐(ℎ)ℎ^ℎ_𝑖=1𝑥_𝑖−𝑚_MVV𝑥_𝑖−𝑚_MVV^𝑡. (5) The approximation of consistency factor can be obtained from elliptical truncation in the multivariate normal distribution based on squared distance. If 𝑥_𝑖∼𝑁(𝜇,Σ),  𝑐(ℎ) is defined as 𝑐(ℎ)=ℎ/𝑛𝑃𝜒²_𝑝+2<𝜒²_{𝑝,ℎ/𝑛}, (6) where 𝜒²_{𝑝,ℎ/𝑛} is the ℎ/𝑛-quantile of 𝜒²_𝑝 distribution. This formula is derived by Butler et al. [13] and further discussed in Croux and Haesbroeck [14] based on the functional form of the MCD estimator. Since MVV have the same functional form with the MCD estimator, we used (6) as the consistency factor for 𝑆_MVV. Albeit guaranteed consistency under normality distribution, Pison et al. [17] cautioned that MCD estimators were biased for small sample sizes. Thus, the consistency factor in (6) only might not be sufficient to make MVV estimator unbiased for small sample sizes. For that reason, we also include the computation of correction factor at any sample size 𝑛 and dimension 𝑝.

2.2. Correction Factor

A simulation study on the effect of correction factor on the MVV estimator is carried out for several sample sizes 𝑛 and dimension 𝑝 = 2, 5, 10, 15, and 20. We generated data sets 𝑋^(𝑗)∈ℝ^𝑛×𝑝 from standard multivariate normal distribution. For each data set 𝑋^(𝑗), 𝑗=1,…,𝑚 we then determine the 𝑐(ℎ)𝑆^(𝑗)_MVV in (5). If the estimator is unbiased, 𝐸[𝑐(ℎ)𝑆_MVV]=𝐼_𝑝, therefore the 𝑝th root of the determinant of 𝑐(ℎ)𝑆_MVV equals 1. Given the mean of the 𝑝th root of the determinant as mean(|𝑐(ℎ)𝑆_MVV|)=(1/𝑚)∑^𝑚_𝑗=1(|𝑐(ℎ)𝑆^(𝑗)_MVV|)^1/𝑝. To determine the correction factor, we performed 𝑚=1000 simulations for different sample sizes 𝑛 and dimensions 𝑝, with 𝛼=0.05 such that 𝜗^𝛼_𝑝,𝑛=1mean𝑐(ℎ)𝑆_MVV. (7) The computed values are displayed in Table 1. Then, using 𝜗^𝛼_𝑝,𝑛 in (7) as the correction factor for 𝑐(ℎ)𝑆_MVV, we obtain 𝜗^𝛼_𝑝,𝑛𝑐(ℎ)𝑆_MVV=𝜗^𝛼_𝑝,𝑛𝑐(ℎ)ℎ^ℎ_𝑖=1𝑥_𝑖−𝑚_MVV𝑥_𝑖−𝑚_MVV^𝑡. (8) Since 𝜗^𝛼_𝑝,𝑛𝑐(ℎ)𝑆_MVV can be considered consistent and unbiased, the determinant of 𝜗^𝛼_𝑝,𝑛𝑐(ℎ)𝑆_MVV should approach 1.

Table 1 

Values of the 𝜗^𝛼_𝑝,𝑛 for 𝑐(ℎ)𝑆_MVV.

3. Investigation through Simulation Experiment

Gather and Becker [21] have emphasized that robust estimators to be used in the method of outliers detection should have sufficient rate of convergence to some true underlying model parameter for consistency and unbiased. A sequence of asymptotically unbiased estimators for parameter 𝜃 is called consistent if lim_𝑛→∞𝑃(|𝜃_𝑛−𝜃|≥𝜀)=0. To illustrate the analysis on the consistency of MVV estimator at multivariate normal, data are randomly generated from 𝑁(0,𝐼_𝑝). An experiment is carried out for several values of sample sizes 𝑛 until convergent for a fixed moderate dimension such that 𝑝=10. Figure 1 shows the determinants of 𝜗^𝛼_𝑝,𝑛𝑐(ℎ)𝑆_MVV corresponding to the sample size, 𝑛. As the value of 𝑛 increases, we can observe that the determinant approaches 1 which implies that the 𝜗^𝛼_𝑝,𝑛𝑐(ℎ)𝑆_MVV is consistent. Next, the investigation using simulation experiment continues to show that 𝑚_MVV and 𝜗^𝛼_𝑝,𝑛𝑐(ℎ)𝑆_MVV which replaced the MLE, 𝜇 and Σ, in Hotelling 𝑇² are consistent and unbiased. The squared distances using any affine-equivariant robust location and scatter estimators which are consistent and unbiased under normal model are asymptotically 𝜒² distributed [21]. Therefore, if 𝑚_MVV and 𝜗^𝛼_𝑝,𝑛𝑐(ℎ)𝑆_MVV are consistent and unbiased estimators for 𝜇 andΣ, then with observations 𝑥_𝑖 i.i.d in ℝ^𝑝∼𝑁_𝑝(𝜇,Σ), it follows that 𝑑²_𝑖=(𝑥_𝑖−𝑚_MVV)𝜗^𝛼_𝑝,𝑛𝑐(ℎ)𝑆⁻¹_MVV(𝑥_𝑖−𝑚_MVV)^𝑡 is asymptotically 𝜒²_𝑝 distributed. Since 𝑑²_𝑖 is similar to Hotelling 𝑇², the asymptotic distribution of the improved Hotelling 𝑇²_MVV when 𝑛→∞ should follow 𝜒²_𝑝 distribution if the estimators are consistent. If we consider a sample of 𝑝 quality characteristics such that 𝑥_𝑖={𝑥_𝑖1,𝑥_𝑖2,…,𝑥_𝑖𝑝} where 𝑖=1,2,…,𝑛 as a phase I data set, then the improved 𝑇²_MVV statistic for 𝑥_𝑖 can be constructed in the following manner: 𝑇²_MVV(𝐼)(𝑖)=𝑥_𝑖−𝑚_MVV𝜗^𝛼_𝑝,𝑛𝑐(ℎ)𝑆_MVV⁻¹𝑥_𝑖−𝑚_MVV^𝑡. (9) To check on the distributions of the improved 𝑇²_MVV, we employed the QQ plots and evaluated the goodness of fit on those plots based on the slope and the 𝑅-square of the straight line as shown in Table 2. The hypothetical distribution represents the 𝜒²_𝑝 without error if all points are in a straight line with slope equals 1 and 𝑅-square also equals 1 [22]. Random data were generated from multivariate standard normal distribution MVN(0,𝐼_𝑝). This study is carried out for the sample size of 𝑛 = 10,000 with dimensions of 𝑝 = 2, 5, 10, 15, and 20. From this table we observe that the 𝑅-square values for all 𝑝’s are 0.999. With regard to the slopes, we can see a considerable difference in the values between the Hotelling’s 𝑇² with original MVV(𝑇²_MVV(𝑜)) and Hotelling’s 𝑇² with improved MVV(𝑇²_MVV(𝐼)) especially when 𝑝 = 2. The slopes for 𝑇²_MVV(𝐼) are consistent and approximately equal to 1 regardless of the dimensions (𝑝). In contrast, the slopes for 𝑇²_MVV(𝑜) are quite a distance away from 1 even though the pattern shows a declining in values towards 1 as 𝑝 increases. We observe that the values for the two measurements (𝑅² and slopes) are very close to the ideal value, which signify that the 𝜒²_𝑝 distribution fits well with the simulated 𝑇²_MVV(𝐼) values. The result implies that the constant 𝜗^𝛼_𝑝,𝑛𝑐(ℎ) fulfills the condition of the multiplicative factors to make the 𝑆_MVV estimators consistent and unbiased for Σ.

Table 2 

The slope and 𝑅-square for 𝑇²_MVV(𝑜) and 𝑇²_MVV(𝐼).

Figure 1 

Determinant of 𝜗𝛼𝑝,𝑛𝑐(ℎ)𝑆MVV when 𝑛→∞ and 𝑝 = 10.

4. MVV Hotelling 𝑇² Control Chart

Let 𝑥_𝑖={𝑥_𝑖1,𝑥_𝑖2,…,𝑥_𝑖𝑝} be the 𝑝-variate random sample of 𝑛 observations of preliminary data set in Phase I. Calculate the 𝑀_MVV and 𝜗^𝛼_𝑝,𝑛𝑐(ℎ)𝑆_MVV estimators. Since the estimators are known to be free from outliers due to their estimation process, they could be readily used as in-control estimators in Phase II. By using these estimates, 𝑇²_MVV(𝐼) statistic in (9) is computed for Phase II observation, 𝑥_𝑔={𝑥_𝑛+1,𝑥_𝑛+2,…} where 𝑥_𝑔∉𝑥_𝑖.

4.1. Estimation of Control Limits

In this section, we present the control limit of the improved 𝑇²_MVV(𝐼) control chart by using simulated data with different combinations of sample sizes, 𝑛, and number of dimensions, 𝑝. The control limit of 𝑇²_MVV(𝐼) chart is then compared with the control limit of 𝑇²_MVV(𝑜) chart, robust Hotelling 𝑇² chart using MCD (𝑇²_MCD) and the traditional Hotelling 𝑇² charts. The application of robust estimators in place of the mean and covariance structure in traditional Phase II Hotelling 𝑇² statistic will cause the distributional properties of the traditional chart to change [9]. To demonstrate the performance of 𝑇²_MVV(𝑜) and 𝑇²_MVV(𝐼) control charts, we need to identify the distribution of each method in order to obtain appropriate control limits, that is, UCL. Since the exact distribution of 𝑇²_MVV(𝑜) is unknown, we apply Monte Carlo method to estimate the quantiles of the 𝑇²_MVV(𝑜) and 𝑇²_MVV(𝐼), for several combinations of sample sizes and dimensions. In order to estimate the 95% quantile of 𝑇²_MVV(𝐼) for a given Phase I of sample size 𝑛 and dimension 𝑝, we generate 𝐾 = 5000 samples of size 𝑛 from a standard multivariate normal distribution, MVN_𝑝(0,𝐼_𝑝). For each data set of size 𝑛, we compute the MVV mean vector and the modified covariance matrix estimates, 𝑀_MVV(𝑘) and 𝜗^𝛼_𝑝,𝑛𝑐(ℎ)𝑆_MVV(𝑘), respectively, from 𝑘=1,…,𝐾. In addition, for each data set, we randomly generate a new observation 𝑥_𝑔,𝑘 treated as a Phase II observation from MVN_𝑝(0,𝐼_𝑝) and calculate the corresponding 𝑇²_MVV(𝐼)(𝑔,𝑘) values. The empirical distribution function of 𝑇²_MVV(𝐼)(𝑔) is based on the simulated values 𝑇²_MVV(𝐼)𝑔,1,𝑇²_MVV𝑔,2,…,𝑇²_MVV𝑔,𝐾. (10) We sort 𝑇²_MVV(𝐼)(𝑔,𝑘) values in ascending order, and the UCL is the 95% quantile of the 5000 statistics. The results of the investigation are presented in Table 3. We observe that the estimated UCLs for 𝑇²_MVV(𝑜) are large as compared to the traditional control charts 𝑇²_𝑆, 𝑇²_MCD, and 𝑇²_RMCD. However, after making the MVV scatter estimator consistent and unbiased as shown in (8), the results improved immensely. As we can see here, the UCLs are closer to the traditional UCLs.

Table 3 

Control limits of the investigated control charts for various combinations of sample sizes and dimensions.

4.2. Real Data Analysis

The application of the improved method 𝑇²_MVV(𝐼) on real data is illustrated using data furnished by Asian Composites Manufacturing Sdn. Bhd. (ACM) which involves in the production of advanced composite panels for the aircraft industry. ACM produces flat and contoured primary (Aileron Skins, Spoilers and Spars) and secondary (Flat Panels, Leading Edges and MISC: Components) structure composite bond assemblies and sub-assemblies for aerospace industries. The company provided us with data on spoilers as shown in Table 4. For the purpose of this study, a sample of 47 spoilers (𝑛=47) which consists of several features, namely, trim edge (𝑋₁), trim edge spar (𝑋₂), and drill hole (𝑋₃) were furnished to us by the company. Out of the total, 21 spoilers were collected from 2009, while the rest were from 2010. Hence, we decided to use the 2009 spoilers as Phase I historical data and considered the spoilers from 2010 as future data in this study. Estimates for the location vector (𝑋) and scatter matrix (𝑆) are presented in Table 5. In the last column of Table 5, we could clearly observe that the upper control limit (UCL) for 𝑇²_MVV(𝐼) is the closest to the traditional Hotelling 𝑇² with values of 11.5513 and 11.035, respectively, whereas the other control charts produce large UCL values especially the original 𝑇²_MVV(𝑇²_MVV(𝑜)). When we compare the improved with the original 𝑇²_MVV, we observe a large disparity between the two values such that 𝑇²_MVV(𝑜) = 41.298, and 𝑇²_MVV(𝐼) = 11.5513. The result indicates great improvement in the UCL values from 𝑇²_MVV(𝑜) to 𝑇²_MVV(𝐼). Table 6 identifies the out-of-control data (bold font) using the different 𝑇² statistics. Among the four statistics, 𝑇²_MVV(𝑜), 𝑇²_MVV(𝐼) and 𝑇²_MCD consider signal observations 20, 22, and 25 as out-of-control but 𝑇² fails to signal observation 22 and only consider, observations 20 and 25 as out-of-control. This is expected due to the low probability of detection in the traditional control chart [16]. For a clearer visualisation on the performance of the control charts in detecting out or control observations, graphical presentation of the corresponding control charts are put on view in Figure 2.

Table 4 

List of Phase I and Phase II real data.

Table 5 

Estimates of location vector, covariance matrix, and UCL.

Table 6 

The Hotelling 𝑇² values for the future (Phase II) data.

(a)

(b)

(c)

(d)

Figure 2 

Hotelling 𝑇2 Control Charts.

5. Conclusion

The UCL value for the Hotelling 𝑇² control chart using consistent and unbiased MVV estimators seemed to improve significantly from the Hotelling 𝑇² control chart based on the original MVV estimators. The improved control chart (𝑇²_MVV(𝐼)) was put to test on real data. Even though the performance of the improved 𝑇²_MVV control chart was on par with the original 𝑇²_MVV chart, the improved estimators have successfully reduced the inflated UCL of the original 𝑇²_MVV close to the UCL of the traditional Hotelling 𝑇² (𝑇²) control chart. However, when the improved control chart was compared with the traditional chart based on their almost equal UCLs, the finding showed that the improved control chart performed better in detecting out-of-control observations. With the good properties and performance, this improved MVV estimators that should be considered as alternative estimators to replace the usual mean and variance vector in the construction of the robust Hotelling 𝑇² control chart as well as other multivariate statistical procedures.