An Improvement of the Hotelling 𝑇^𝟐 Statistic in Monitoring Multivariate Quality Characteristics

Abstract

The Hotelling 𝑇² statistic is the most popular statistic used in multivariate control charts to monitor multiple qualities. However, this statistic is easily affected by the existence of more than one outlier in the data set. To rectify this problem, robust control charts, which are based on the minimum volume ellipsoid and the minimum covariance determinant, have been proposed. Most researchers assess the performance of multivariate control charts based on the number of signals without paying much attention to whether those signals are really outliers. With due respect, we propose to evaluate control charts not only based on the number of detected outliers but also with respect to their correct positions. In this paper, an Upper Control Limit based on the median and the median absolute deviation is also proposed. The results of this study signify that the proposed Upper Control Limit improves the detection of correct outliers but that it suffers from a swamping effect when the positions of outliers are not taken into consideration. Finally, a robust control chart based on the diagnostic robust generalised potential procedure is introduced to remedy this drawback.

1. Introduction

In statistical quality control, a process changes into an out-of-control situation when outliers appear in two different ways, namely, outliers that are randomly distributed within a data set and outliers that sequentially occur after a specific observation during a specific period of time in the data set. The former and the latter situations are referred to as scatter outliers and a sustained step shift, respectively, [1–4].

The detection of correct outliers in phase I of the monitoring scheme is crucial. If outliers are not correctly detected, the result leads to model misspecification and to incorrect results during phase II [5]. The Hotelling𝑇²statistic, which was first introduced by Hotelling in [6], is the most popular statistic used in multivariate control charts to monitor multiple quality characteristics [7–11].

Vargas [12] demonstrated that a𝑇²statistic based on the usual classical estimators fails to detect multiple scatter outliers for individual observations (𝑛=1), although this statistic is effective in the presence of a small number of outliers. It is now evident that the Hotelling𝑇²statistic, which is based on the usual classical sample mean vector and variance-covariance matrix, is easily affected by the existence of more than one outlier in the HDS (Historical Data Set). In addition, the 𝑇²statistic suffers from a masking or swamping effect [13, 14]. Sullivan and Woodall [10] showed that the 𝑇² statistic based on the usual sample variance-covariance matrix for individual observations is not only less effective in detecting scatter outliers in the HDS but also poor in sustained step shifts in the mean vector.

Robust methods for multivariate data, based on the MVE and MCD, have been widely used in regression contexts for diagnosing influential observations and high leverage points and outliers [14, 15] but have only recently been applied to multivariate quality control process applications. Vargas [12] employed the minimum volume ellipsoid (MVE) and the minimum covariance determinant (MCD) as two robust estimates of location and dispersion [16], instead of the usual classical sample mean vector and covariance matrix in the Hotelling𝑇²statistic. The application of robust control charts for individual observations based on the MVE and MCD has also been discussed extensively by Jensen et al. [17]. It is worth mentioning that there is no guarantee that the mathematical distribution of the 𝑇² statistic is preserved by replacing the location and scale estimators with robust versions. To remedy this problem, Vargas [12] and Jensen et al. [17] used an empirical distribution of the robust 𝑇²statistic for calculating the emprical upper control limits (UCLs) of their proposed robust control charts.

We have seen the application of MVE and MCD methods in the development of robust control charts. Other methods such as outlier identification in high dimensions [18] and some proposed multivariate outlier detection techniques [19] may also be considered in the control process applications. Our main aim in this paper is to propose a robust multivariate control chart based on the diagnostic robust generalised potential (DRGP), which was initiated by Habshah et al. [20]. This is the first attempt to introduce another robust control chart as an alternative to the two existing robust MVE-based and MCD-based control charts. Hence, our focus in this paper is only limited to the DRGP-based control chart and compare its performance with the two preceding charts. We do not wish to compare these charts with other charts based on other methods mentioned in Wilcox [19] and Filzmoser et al. [18].

Vargas [12] and Jensen et al. [17] evaluated the performance of each control chart based on the number of detected outliers, regardless of whether they came from a correct outlier positions. In other words, their work is only devoted to the detection of outliers regardless of whether the outliers are true (correct outlier position) or false outliers. Their work has motivated us to consider the identification of correct outliers when evaluating different control charts. In this regard, we introduce another empirical method for calculating the UCLs of the robust 𝑇²statistic based on the median and the MAD of the estimators.

2. Diagnostics for the Identification of High Leverage Points

High leverage points are a type of influential observation that is substantially different for one or more predictor variables [21, 22]. It is now evident that high leverage points are responsible for leading to model misspecifications and misleading results [23–25].

There are some methods in the literature for identifying high leverage points in linear regression models [14, 22, 25–28]. The 𝑖th diagonal elements of hat matrix 𝐻 are referred to as leverage points and are denoted byℎ_𝑖𝑖=𝑣^𝑇_𝑖𝑉^𝑇𝑉𝑣_𝑖,𝑖=1,2,…,𝑚, (2.1) where 𝑉 is an 𝑚×𝑘 matrix ofpredictor variables of regression model. Hoaglin and Welsch [29] considered observations to be high leverage points when ℎ_𝑖𝑖 are greater than 2𝑘/𝑚. Hadi [30] introduced another measure to diagnose high leverage points, which is known as “potential measures.” According to Hadi [30], the 𝑖th potential is defined as follows:𝑝_𝑖𝑖=𝑣^𝑇_𝑖𝑉^𝑇_(𝑖)𝑉_(𝑖)⁻¹𝑣_𝑖, (2.2) where 𝑉_(𝑖) is the data matrix with its 𝑖th row deleted. By using simple matrix algebra, it is easy to obtain a relationship between the potentials and the diagonal elements of 𝐻, as follows:𝑝_𝑖𝑖=ℎ_𝑖𝑖1−ℎ_𝑖𝑖. (2.3)

Hadi [30] suggested a confidence bound cutoff point for 𝑝_𝑖𝑖 as follows:Median𝑝_𝑖𝑖+𝑐MAD𝑝_𝑖𝑖, (2.4) whereMAD(𝑝_𝑖𝑖)=Median{|𝑝_𝑖𝑖−Median(𝑝_𝑖𝑖)|}/0.6745 and 𝑐 is a constant that is chosen between 2 or 3, as appropriate. Robust version of the Mahalanobis distance is also being used to identify high leverage points [14, 18, 20]. Habshah et al. [20] pointed out that although these robust diagnostic techniques can rectify the masking problem, they are affected by the swamping effect, which is not desirable either.

To remedy this problem, Habshah et al. [20] proposed a unified approach, which is called the diagnostic robust generalised potential (DRGP) method which accommodates both the diagnostic and robust approaches together. The robust approach is utilised to detect the suspected high leverage points and then diagnostic approach is utilised to confirm our suspicion. The DRGP partitions the data into two sets. The first set consists of suspicious cases, which are deleted from the original observations (denoted by 𝐷), and the second set contains the remaining data (denoted by 𝑅). It is clear that if 𝑑 is the number of cases that includes 𝐷, then the 𝑅 set contains 𝑚−𝑑 observations and𝑑<𝑚−𝑘. Without loss of generality, we assume that 𝑑 cases are placed in the last 𝑑 rows of 𝑌 and 𝑉 so that the hat matrix is partitioned and that the 𝑖th deletion leverage is defined as follows:ℎ^(−𝐷)_𝑖𝑖=𝑣^𝑇_𝑖𝑉^𝑇_𝑅𝑉_𝑅⁻¹𝑣_𝑖, (2.5) where 𝑉_𝑅 indicates the remaining observation matrix [31]. By considering (2.2), (2.3) and (2.5), the generalised potential is defined as follows:𝑝^∗_𝑖𝑖=ℎ^(−𝐷)_𝑖𝑖1−ℎ^(−𝐷)_𝑖𝑖forℎ𝑖∈𝑅_{(−𝐷)𝑖𝑖}for𝑖∈𝐷. (2.6)

Similar to the potential values and with regard to (2.4), the cutoff point for 𝑝^∗_𝑖𝑖 isMedian𝑝^∗_𝑖𝑖+𝑐MAD𝑝^∗_𝑖𝑖. (2.7)

It is worth mentioning that the DRGP employs the robust Mahalanobis distance in the first step to detect high leverage points as preliminary suspicious observations; these points are placed in the 𝐷 set. Next, in the second step, only the cases that are greater than (2.7) are reported as the final detections. In the next section, the DRGP is employed to effectively detect outliers in multivariate quality control charts for individual observations.

3. Multivariate Robust 𝑇²Control Charts

Suppose that there is an HDS in the phase I monitoring scheme that consists of 𝑚 time-ordered observation vectors of dimension𝑝, which are observed independently, where 𝑝 is the number of quality characteristics that are measured (𝑝<𝑚). It is assumed that each vector comes from a 𝑝-variate normal distribution. Thus, if 𝑋_𝑖∈ℝ^𝑝 is a vector in the HDS for the 𝑖th time period, 𝑋_𝑖∼𝑁_𝑝(𝜇,Σ), where 𝜇 and Σ are the population mean vector and the variance-covariance matrix, respectively. As mentioned earlier, the Hotelling𝑇²statistic is used to detect outliers in multivariate control charts. The general form of this statistic is𝑇²_𝑖=𝑋_𝑖−𝜇^𝑇−1𝑋_𝑖−𝜇,𝑖=1,2,…,𝑚. (3.1)

Because the parameters in (3.1) are usually unknown, the usual sample mean and variance-covariance matrix are used as the classical estimations of 𝜇 and Σ. In practice, these variables are expressed by1𝑋=𝑚^𝑚_𝑖=1𝑋_𝑖,1𝑆=𝑚−1^𝑚_𝑖=1𝑋_𝑖−𝑋𝑋_𝑖−𝑋^𝑇. (3.2)

In phase I, the parameters are retrospectively estimated based on the current HDS; as a result, the vector 𝑋_𝑖 is not independent of the estimators 𝑋 and 𝑆. In this situation, the statistical distribution of (3.1) is given as𝑇²∼(𝑚−1)²𝑚𝐵𝑝2,𝑚−𝑝−12, (3.3) where 𝐵(𝑝/2,(𝑚−𝑝−1)/2) represents a beta distribution with parameters 𝑝/2 and (𝑚−𝑝−1)/2 [8, 11]. From (3.3), the upper control limit (UCL) of 𝑇²is((𝑚−1)²/𝑚)𝐵(𝛼,𝑝/2,(𝑚−𝑝−1)/2), where𝛼 is the probability of a false alarm for each point plotted on the control chart and 𝐵(𝛼,𝑝/2,(𝑚−𝑝−1)/2) is the𝛼th upper quantile of the beta distribution with parameters 𝑝/2 and (𝑚−𝑝−1)/2.

The lower control limit (LCL) is often set to zero [8, 32]. It should be noted that the aforementioned UCLs are exact when applied to a single point in phase I, whereas phase I is a retrospective analysis of all observations. Therefore, the values of 𝛼 cannot be applied to a set of points. In this situation, if all of the statistics were distributed independently, then the overall probability of a false alarm would be𝛼^′=1−(1−𝛼)^𝑚, (3.4) where 𝛼 indicates the probability of a false alarm, which is assigned for each observation plotted on the control chart in a subgroup of size 𝑚. In practice, it is reasonable to determine the UCL by simulation to give a specified overall false alarm [5, 33, 34]. Hereafter, in this paper, we still refer to the overall false, 𝛼^′ as 𝛼 for simplicity.

The use of (3.1) is not effective in the presence of multiple outliers, so robust alternatives are proposed. In this regard, two of the recently proposed robust alternative approaches to 𝑇² are based on the MVE and the MCD estimators, which will be denoted by 𝑇²_mve and 𝑇²_mcd, respectively, and are defined as follows:𝑇²_mve,𝑖=𝑋_𝑖−𝑋_mve^𝑇𝑆⁻¹_mve𝑋_𝑖−𝑋_mve²,𝑖=1,2,…,𝑚,𝑇²_mcd,𝑖=𝑋_𝑖−𝑋_mcd^𝑇𝑆⁻¹_mcd𝑋_𝑖−𝑋_mcd²,𝑖=1,2,…,𝑚, (3.5)(3.6) where 𝑋_mve and 𝑋_mcd are the robust estimations of the sample mean and 𝑆⁻¹_mve and 𝑆⁻¹_mcd are the corresponding estimators of the sample variance-covariance matrix. As previously mentioned in Section 1, due to the unknown distribution of 𝑇²_mve and 𝑇²_mcd, empirical methods are used to determine the UCLs. The empirical simulated UCLs for (3.5) or (3.6) are usually determined by finding the 𝛼th upper quantile of the empirical distribution of the corresponding statistic. For this situation, Jensen et al. [17] and Vargas [12] defined 𝛼 as the overall false alarm. Following the idea of Habshah et al. [20], another empirical UCL is proposed as follows:Median𝑇²_mve/mcd,𝑖+𝑐MAD𝑇²_mve/mcd,𝑖, (3.7) where MAD is the median absolute deviation of either 𝑇²_mve,𝑖 or the 𝑇²_mcd,𝑖, as defined in (2.4). For simplicity, the first empirical UCL is referred to as Empr, and the proposed UCL will be denoted by Med-Mad. As will be shown later, (3.7) tends to declare too many observations as outliers, by detecting outliers regardless of whether the detected outliers are true (in a correct outlier position) or false, even though (3.7) has a better performance in detecting real outliers at their correct positions in the HDS than Empr UCL. In this regard, we propose to apply the DRGP based on the MVE and MCD with the Med-Mad UCLs to reduce the number of undue signals caused by the detection of outliers irrespective of their correct positions and, at the same time, to effectively detect correct outliers.

Vargas [12] and Jensen et al. [17] employed the probability of signals to evaluate and compare control charts based on 𝑇²_mve and 𝑇²_mcd. Their work is only based on the number of detections, without considering whether those detections are correct outliers. It is worth mentioning that the probability of signals cannot be properly judged if there is a swamping effect in the monitoring scheme. In the following section, we will present our proposed control scheme and explain how it can detect real outliers at their correct positions.

4. Simulation Study

In this section, a Monte Carlo simulation study is carried out to assess the performance of the control schemes discussed previously. The simulation is designed based on three subsamples, each of size 𝑚=30, 50, and 100, with a number of characteristics 𝑝=2, 3, 5, and 10. Let us assume that the in-control process is a 𝑝-variate normal distribution with mean vector 𝜇₀and covariance matrix Σ.

The simulated empirical UCLs are obtained by generating 5000 in-control data sets for each combination of 𝑚 and 𝑝. Due to the affine equivariant property of the 𝑇²statistics, these limits are applicable to any values of 𝜇 and Σ. Many researchers, such as Jensen et al. [17] and Vargas [12], have determined the Empr UCLs by calculating all of the𝑇²statistics for each observation in generated data sets of each subgroup of size 𝑚 and recording the maximum value of the𝑇²statistics. Subsequently, the upper 𝛼th percentile of the 5000 recorded maximum values of the𝑇²statistics is declared as the Empr UCL. In this manner, they defined 𝛼 as the overall false alarm. In this paper, we consider the probability of the overall false alarm, which is equal to𝛼=0.05.

To make the overall false alarm of the control charts based on the Med-Mad UCLs equivalent to Empr UCLs, which is equal to 𝛼=0.05, the following steps are considered. The Med-Mad UCLs are attained based on all 𝑚×5000 values of the simulated𝑇² statistics. It is worth mentioning here that, to keep the overall false alarm at 𝛼=0.05, the probability of a false alarm for each observation in each subgroup of sizem is calculated based on (3.4). Hence, to have the same overall false alarm, the values for 𝑐 in (3.7) and for 𝑚=30, 50, and 100 must be chosen as Φ⁻¹_(0.002)≃2.88, Φ⁻¹_(0.001)≃3, and Φ⁻¹_(0.001)≃3, respectively. The Empr and Med-Mad UCLs of 𝑇², 𝑇²_mve, and 𝑇²_mcd for each combination of 𝑚 and 𝑝 are presented in Table 1.

Table 1 

The simulated UCLs for all of the 𝑇² statistics.

Then, a contaminated data set of size 𝑚 in the 𝑝 dimension is generated for different values of the noncentrality parameter (ncp). The out-of-control process is a p-variate normal distribution with the same covariance matrix but with a shifted mean vector of 𝜇₁. Thus, the variation here remains stable. The magnitude of the shift is measured by a scalar defined as follows:𝜇₁−𝜇₀^𝑇Σ⁻¹𝜇₁−𝜇₀. (4.1)

This measure is called the noncentrality parameter and is hereafter referred to as ncp. Four outlier percentage levels are considered, which are denoted by𝜖=5%, 10%, 15%, and 20%. It is clear from (4.1) that the severity of the shift only depends on the values of 𝜇₁. Hence, without loss of generality, it can be assumed that 𝜇₀ is a zero vector and Σ is a 𝑝×𝑝 identity matrix.

The control charts are assessed based on the proposed criterion, which are based on the number of true detected outliers with regard to the correct position of the generated out-of-control observations in the data set. The number of outliers detected without regard to their positions is also presented for comparison.

Repeating this process 5000 times, the number of detections and the number of correctly detected outliers with correct positions are recorded for each replication. The numbers of detected outliers are determined by comparing each of the 𝑇² values with the respective UCLs given in Table 1. Each detected outlier is checked by its position in the data set to determine whether it can truly be generated from the intentional simulated contaminated points in the data set. The true or correctly detected outliers refer to the outliers detected at the correct position. The number of detected outliers simply indicates the outliers that have been detected irrespective of their correct position. The average number of detections over 5000 iterations is presented for 𝑝=2, 𝑚=30, 50, and 100 in Tables 2, 3, and 4. It is important to note that the presented values were rounded up to two digits. The values in parentheses represent the number of outliers detected regardless of their correct positions.

Table 2 

The number of correctly detected outliers and the number of detected outliers for 𝑝=2 and 𝑚=30.

Table 3 

The number of correctly detected outliers and the number of detected outliers for 𝑝=2 and 𝑚=50.

Table 4 

The number of correctly detected outliers and the number of detected outliers for 𝑝=2 and 𝑚=100.

Let us first focus on the results of 𝑇², 𝑇²_mve, and 𝑇²_mcd obtained by using the Empr and Med-Mad UCLs. Following these results, we will see later why the DRGP approach is proposed.

As can be seen from these tables, the UCLs based on the Med-Mad approach for all three control charts have better performance in detecting the real outliers at their correct positions, compared to Empr UCL. The classical 𝑇² chart, based on Empr UCL, performs very poorly. It can be seen that the Med-Mad UCLs are more reliable in detecting the correctly detected outliers for the robust control charts compared to Empr UCLs, particularly when the percentage of outliers increases. However, the results signify that both 𝑇²_mve and 𝑇²_mcd based on the Med-Mad UCLs suffer from a swamping effect due to the detection of more outliers without regard to their correct positions.

For example, at 10%, regarding outliers in the HDS with 𝑚=100, 𝑝=2, and ncp=25 and both 𝑇²_mve and 𝑇²_mcd based on the Med-Mad UCLs, the methods detect exactly 10 outliers at the correct positions, but they also detect 15 outliers irrespective of their correct positions. It is interesting to note that the performance of the robust control charts based on the Empr UCLs is reasonably close to that of the robust control charts based on the Med-Mad UCLs for very large values of ncp, such as ncp≥45. The 𝑇²_mve and 𝑇²_mcd parameters, based on the Med-Mad UCLs, are equally good at detecting correct outliers. Nonetheless, with increasing subgroup size, 𝑇²_mcd based on the Med-Mad UCL is slightly better than the 𝑇²_mve based on Med-Mad, particularly for a large percentage of outliers.

We can see that although the classical 𝑇²Med-Mad-based method is better than the 𝑇²Empr-based method, it detects a smaller number of exact outliers as the percentage of outliers increases. These results are consistent with other values of 𝑝 but are not reported here due to space constraints. The findings of Tables 2 to 4 seem to suggest that the Med-Mad UCLs are more reliable than Empr UCLs in detecting correct outliers. However, the Med-Mad UCLs detected more outliers irrespective of their correct positions due to swamping effects. In other words, we have shown that the robust control charts based on the Med-Mad UCLs effectively detect the number of correct outliers, but they overdetected outliers irrespective of their correct positions. As such, we need to employ control charts that can reduce such undue detections. In this regard, we suggest applying the DRGP procedure discussed in Section 2. The same simulation procedure was then carried out, and the DRGP was applied to the data sets.

The results obtained by using the DRGP approach are exhibited in Table 5. Due to space limitations, the results are presented only for 𝑚=30, 100 and 𝜖=5%, 20%. As can be seen from Table 5, there is a steady decrease in the number of undue observations when the DRGP approach is applied. For example, as shown in Table 4, the total number of detections by the Med-Mad UCL with ncp=55 is 24 for the MVE and MCD, while it decreases to 21 in Table 5. To simplify the presentation of results, the proportion of correctly detected outliers to detected outliers is calculated and referred to as the correct detection rate. The values of the correct detection rates for 𝑚=50 and 𝑝=2 are shown in Table 6. The results indicate that the DRGP approach provides higher correct detection rates compared to the other methods.

Table 5 

The number of correctly detected outliers and the number of outliers detected using the DRGP when 𝑝=2.

Table 6 

The correct detection rate for 𝑚=50 and 𝑝=2.

However, these results were not very encouraging for small shifts (ncp=5). It can be seen from Table 6 that there is a gradual and steady rise in the correct detection rate with increasing outlier percentage. The results are similar for other values of 𝑚 and 𝑝, which are not tabulated due to space limitations. For more clarification, the correct detection rates for various values of 𝑚 and 𝑝 are plotted in Figures 1 and 2. These figures confirm that the DRGP approach gives a higher correct detection rate.

(a)

(b)

Figure 1 

The correct detection rate for different outlier percentages for 𝑝=2: 𝑚=30 (a) and 𝑚=100 (b).

(a)

(b)

Figure 2 

The correct detection rate for different outlier percentages for 𝑝=3∶ 𝑚=30 (a) and 𝑚=100 (b).

5. Numerical Example

In this section, a numerical example is introduced to assess the performance of our method. This is a bivariate data set which is taken from Shewhart [35]. It presents the measurements of the depth of sapwood and the depth of penetration of creosote in telephone poles. The subgroup size of the original dataset is 10 and we only focus on the first column of the data set based on 20 subgroups. The𝑇²,𝑇²_mve,𝑇²_mcd, and the DRGP statistics were then applied to the data set. The Empr UCLs and the Med-Mad UCLs for 𝑇², 𝑇²_mve, and 𝑇²_mcd statistics are 9.010, 31.723.738, and 113.140 and 7.010, 5.910, and 9.461, respectively. These UCLs are calculated using the simulation, as discussed in Section 4, for 𝛼=0.05.

Figure 3 shows the different 𝑇²control charts for Empr UCLs. As can be seen from this graph, none of the control charts based on the Empr UCL is able to detect any outlier in the data set. On the other hand, robust 𝑇² statistics based on the Med-Mad UCLs can identify severed outliers (Figure 4).

Figure 3 

Control charts for the creosoting Telephone Poles Data, based on the Empr UCLs.

Figure 4 

Control charts for the creosoting Telephone Poles Data, based on the Med-Mad UCLs.

The points 6 and 16 are detected by 𝑇²_mve as outliers, and 𝑇²_mcd identified the cases of 1, 6, 8, 16, and 21. The robust control charts based on the DRGP approach detected one observation as outlier which is observation number 6. The results are not included here due to space limitations.

6. Conclusions

Most research studies evaluate the performance of robust multivariate 𝑇²_mve and 𝑇²_mcd control charts based on the number of outliers detected without regard to the correct position of outliers in the data set. The detection of real outliers (outliers at the correct position in the data set) is crucial to avoid making wrong inferences. This study has shown that although the 𝑇²_mve and the 𝑇²_mcd based on the Med-Mad UCLs are effective in the detection of correct outliers, they have the tendency to declare undue observations as outliers (irrespective of whether they are true or false outliers) due to a swamping effect. In this respect, in the evaluation of robust control charts, we suggest not only a consideration of the number of outlier detections but also a consideration of the correct position of the detected outliers.

Our findings also suggest that the proposed Med-Mad UCLs have better performance than the commonly used Empr UCLs in detecting outliers with regard to the correct position of outliers, especially for higher proportions of outliers. The practical finding of this paper is that the robust control chart based on the DRGP with the proposed Med-Mad UCLs gives credible performance.

The limitations of our study are that inferences or conclusions are only confined to the detection of multiple outliers for individual observations, scatter outlier situations, and moderate dimensional data sets (𝑝≤10).