# Robust Control Charts for Monitoring Process Mean of Phase-I Multivariate Individual Observations

Authors: Asokan Mulayath Variyath; Jayasankar Vattathoor
Journal: Journal of Quality and Reliability Engineering
 (2013)

---

## Abstract

Hoteling'sT2 control charts are widely used in industries to monitor multivariate processes. The classical estimators, sample mean, and the sample covariance used in T2 control charts are highly sensitive to the outliers in the data. In Phase-I monitoring, control limits
are arrived at using historical data after identifying and removing the multivariate outliers. We propose Hoteling's T2 control charts with high-breakdown robust estimators based on
the reweighted minimum covariance determinant (RMCD) and the reweighted minimum
volume ellipsoid (RMVE) to monitor multivariate observations in Phase-I data. We assessed the performance of these robust control charts based on a large number of Monte
Carlo simulations by considering different data scenarios and found that the proposed control charts have better performance compared to existing methods.


---

## Body

## 1. Introduction

Control charts are widely used in industries to monitor/control processes. Generally, the construction of a control chart is carried out in two phases. The Phase-I data is analyzed to determine whether the data indicates a stable (or in-control) process and to estimate the process parameters and thereby the construction of control limits. The Phase-II data analysis consists of monitoring future observations based on control limits derived from the Phase-I estimates to determine whether the process continues to be in control or not. But trends, step changes, outliers, and other unusual data points in the Phase-I data can have an adverse effect on the estimation of parameters and the resulting control limits. That is, any deviation from the main assumption (in our case, identically and independently distributed from normal distribution) may lead to an out-of-control situation. Therefore, it becomes very important to identify and eliminate these data points prior to calculating the control limits. In this paper, all these unusual data points are referred to as “outliers.”Multivariate quality characteristics are often correlated, and to monitor the multivariate process mean Hoteling’sT2 control chart [1, 2] is widely used. To implement Hoteling's T2 control chart for individual observations in Phase-I, for each observation xj we calculate
(1)T2(xj)=(xj-x-)′S1-1(xj-x-),
where xj = (xj1,xj2,…,xjp)′ is the jth p-variate observation, (j=1,2,…,m) and the sample mean x-, sample covariance matrix S1 are based on m Phase-I observations. In Phase-I monitoring, the T2(xj) values are compared with the T2 control limit derived by assuming that the xj’s are multivariate normal so that the T2 control limits are based on the beta distribution with the parameters p/2 and (m-p-1)/2. However, the classical estimators, sample mean, and sample covariance are highly sensitive to the outliers, and hence robust estimation methods are preferred as they have the advantage of not being unduly influenced by the outliers. The use of robust estimation methods is well suited to detect multivariate outliers because of their high breakdown points which ensure that the control limits are reasonably accurate. Sullivan and Woodall [3] proposed a T2 chart with an estimate of the covariance matrix based on the successive differences of observations and showed that it is effective in detecting process shift. However, these charts are not effective in detecting multiple multivariate outliers because of their low breakdown point.Vargas [4] introduced two robust T2 control charts based on robust estimators of location and scatter, namely, the minimum covariance determinant (MCD) and minimum volume ellipsoid (MVE) for identifying the outliers in Phase-I multivariate individual observations. Jensen et al. [5] showed that TMCD2 and TMVE2  control charts have better performance when outliers are present in the Phase-I data. Chenouri et al. [6] used reweighted MCD estimators for monitoring the Phase-II data, without constructing Phase-I control charts. However, in many situations Phase-I control charts are necessary to assess the performance of the process and also to identify the outliers. We propose T2 control charts based on the reweighted minimum covariance determinant (RMCD)/reweighted minimum volume ellipsoid (RMVE) (TRMCD2/TRMVE2) for monitoring Phase-I multivariate individual observations. RMCD/RMVE estimators are statistically more efficient than MCD/MVE estimators and have a manageable asymptotic distribution. We empirically arrive at Phase-I control limits for the TRMCD2/TRMVE2 control chart for some specific sample sizes and fitted a nonlinear model to determine control limits for any sample size for dimensions 2 to 10. Our simulation studies show that TRMCD2/TRMVE2 control charts are performing well compared to TMCD2/TMVE2 control charts for monitoring the Phase-I data.The organization of the remaining part of the paper is as follows. In Section2, we discuss the properties of a good robust estimator and we briefly explain the MCD/MVE estimators and their reweighted versions. The proposed TRMCD2/TRMVE2 control charts are given in Section 3 along with the control limits arrived at based on Monte Carlo simulations. We assess the performance of the proposed control charts in Section 4, and the implementation of the proposed methods is illustrated in a case example in Section 5. Our conclusions are given in Section 6.


## 2. Robust Estimators

The affine equivariance property of the estimator is important because it makes the analysis independent of the measurement scale of the variables as well as the transformations or rotations of the data. The breakdown point concept introduced by Donoho and Huber [7] is often used to assess the robustness. The breakdown point is the smallest proportion of the observations which can render an estimator meaningless. A higher breakdown point implies a more robust estimator, and the highest attainable breakdown point is 1/2 in the case of median in the univariate case. For more details on affine equivariance and breakdown points one may refer to Chenouri et al. [6] or Jensen et al. [5].An estimator is said to be relatively efficient compared to any other estimator if the mean square error for the estimator is the least for at least some values of the parameter compared to others. A robust estimator is considered to be good if it carries the property of affine equivariance along with a higher breakdown point and greater efficiency. In addition to the above three properties of a good robust estimator, it should be possible to calculate the estimator in a reasonable amount of time to make it computationally efficient.It is difficult to get an affine equivariant and robust estimator as affine equivariance and high breakdown will not come simultaneously. Lopuhaä and Rousseeuw [8] and Donoho and Gasko [9] showed that the finite sample breakdown point of (m-p+1)/(2m-p+1) is difficult for an affine equivariant estimator. The largest attainable finite sample breakdown point of any affine equivariant estimator of the location and scatter matrix with a sample size m and dimension p is (m-p+1)/2m [10]. Therefore relaxing the affine equivariance condition of the estimators to invariance under the orthogonal transformation makes it easy to find an estimator with the highest breakdown point.The classical estimators, sample mean vector, and covariance matrix of location and scatter parameters are affine equivariant but their sample breakdown point is as low as1/m. The MCD and MVE estimators have the highest possible finite sample breakdown point (m-p+1)/2m. However, both of these estimators have very low asymptotic efficiency under normality. But the reweighted versions of MCD and MVE estimators have better efficiency without compromising on the breakdown point and rate of convergence compared to MCD and MVE. In the next two subsections, we discuss in detail about the MCD and MVE estimators and their reweighted versions.


### 2.1. MCD and RMCD Estimators

The MCD estimators of location and scatter parameters of the distribution are determined by a two-step procedure. In step 1, all possible subsets of observations of sizeh=(m*γ),  where0.5≤γ≤1 are obtained. In step 2, the subset whose covariance matrix has the smallest possible determinant is selected. The MCD location estimator x-MCD is defined as the average of this selected subset of h points, and the MCD scatter estimator is given by SMCD=aγ,p*bγ,pm*CMCD, where CMCD is the covariance matrix of the selected subset, the constant aγ,p is the multiplication factor for consistency [11], and bγ,pm is the finite sample correction factor [12]. Here (1-γ) represents the breakdown point of the MCD estimators. The MCD estimator has its highest possible finite sample breakdown point when h=(m+p+1)/2 and has an m-1/2 rate of convergence but has a very low asymptotic efficiency under normality. Computing the exact MCD estimators (x-MCD, SMCD) is computationally expensive or even impossible for large sample sizes in high dimensions [13], and hence various algorithms have been suggested for approximating the MCD. Hawkins and Olive [14] and Rousseeuw and van Driessen [15] independently proposed a fast algorithm for approximating MCD. The FAST-MCD algorithm of Rousseeuw and van Driessen finds the exact MCD for small datasets and gives a good approximation for larger datasets, which is available in the standard statistical software SPLUS, R, SAS, and Matlab.MCD estimators are highly robust, carry equivariance properties, and can be calculated in a reasonable time using the FAST-MCD algorithm; however, they are statistically not efficient. The reweighted procedure will help to carry both robustness and efficiency. That is, first a highly robust but perhaps an inefficient estimator is computed, which is used as a starting point to find a local solution for detecting outliers and computing the sample mean and covariance of the cleaned data set as in Rousseeuw and van Zomeren [16]. This consists of discarding those observations whose Mahalanobis distances exceed a certain fixed threshold value. MCD is the current best choice for the initial estimator of a two-step procedure as it contains the robustness, equivariance, and computational efficiency properties along with its m-1/2 rate of convergence. Hence RMCD estimators are the weighted mean vector(2)x-RMCD=(∑j=1mwjxj)(∑j=1mwj),
and the weighted covariance matrix
(3)SRMCD=cα,p*dγ,αm,p*∑j=1mwj(xj-x-RMCD)(xj-x-RMCD)′∑j=1mwj,
where cα,p is the multiplication factors for consistency [11], dγ,αm,p is the finite sample correction factor [12], and the weights wj are defined as
(4)wj={1ifRD(xj)≤qα,0otherwise,
where the robust distance RD(xj)=(xj-x-MCD)′SMCD-1(xj-x-MCD) and qα is (1-α)100% quantile of the chi-square distribution with p degrees of freedom.This reweighting technique improves the efficiency of the initial MCD estimator while retaining (most of) its robustness. Hence the RMCD estimator inherits the affine equivariance, robustness, and asymptotic normality properties of the MCD estimators with an improved efficiency.


### 2.2. MVE and RMVE Estimators

Determining the MVE estimators of location and scatter parameters of the distribution is almost in line with that of the MCD estimator. As in the case of MCD, all the possible subsets of data points with sizeh=(m*γ)  (where0.5≤γ≤1) is obtained first. Then the ellipsoid of minimum volume that covers the subsets are obtained to determine the MVE estimators. The MVE location estimator is the geometrical center of the ellipsoid, and the MVE scatter estimator is the matrix defining the ellipsoid itself, multiplied by an appropriate constant to ensure consistency [13, 16]. Thus MVE estimator does not correspond to the sample mean vector and the sample covariance matrix as in the case of the MCD estimator. Here (1-γ) represents the breakdown point of the MVE estimators, as in the case of MCD, and it has the highest possible finite sample breakdown point when h=(m+p+1)/2m [8, 17]. The MVE estimator has an m-1/3 rate of convergence and a nonnormal asymptotic distribution [17].As in the case for MCD estimators, MVE estimators are also not efficient. Hence, a reweighted version similar to that for MCD has been proposed by Rousseeuw and van Zomeren [16]. Note that it has been shown more recently that the RMVE estimators do not improve on the convergence rate (and thus the 0% asymptotic efficiency) of the initial MVE estimator [8, 12]. Therefore, as an alternative, a one-step M-estimator can be calculated with the MVE estimators as the initial solution [13, 18] which results in an estimator with the standard m-1/2 convergence rate to a normal asymptotic distribution. For more details on MCD/MVE estimators one may refer to Chenouri et al. [6] or Jensen et al. [5]. The algorithm to determine the MVE/RMVE estimators is available in the statistical software SPLUS, R, SAS, and Matlab.


## 2.1. MCD and RMCD Estimators

The MCD estimators of location and scatter parameters of the distribution are determined by a two-step procedure. In step 1, all possible subsets of observations of sizeh=(m*γ),  where0.5≤γ≤1 are obtained. In step 2, the subset whose covariance matrix has the smallest possible determinant is selected. The MCD location estimator x-MCD is defined as the average of this selected subset of h points, and the MCD scatter estimator is given by SMCD=aγ,p*bγ,pm*CMCD, where CMCD is the covariance matrix of the selected subset, the constant aγ,p is the multiplication factor for consistency [11], and bγ,pm is the finite sample correction factor [12]. Here (1-γ) represents the breakdown point of the MCD estimators. The MCD estimator has its highest possible finite sample breakdown point when h=(m+p+1)/2 and has an m-1/2 rate of convergence but has a very low asymptotic efficiency under normality. Computing the exact MCD estimators (x-MCD, SMCD) is computationally expensive or even impossible for large sample sizes in high dimensions [13], and hence various algorithms have been suggested for approximating the MCD. Hawkins and Olive [14] and Rousseeuw and van Driessen [15] independently proposed a fast algorithm for approximating MCD. The FAST-MCD algorithm of Rousseeuw and van Driessen finds the exact MCD for small datasets and gives a good approximation for larger datasets, which is available in the standard statistical software SPLUS, R, SAS, and Matlab.MCD estimators are highly robust, carry equivariance properties, and can be calculated in a reasonable time using the FAST-MCD algorithm; however, they are statistically not efficient. The reweighted procedure will help to carry both robustness and efficiency. That is, first a highly robust but perhaps an inefficient estimator is computed, which is used as a starting point to find a local solution for detecting outliers and computing the sample mean and covariance of the cleaned data set as in Rousseeuw and van Zomeren [16]. This consists of discarding those observations whose Mahalanobis distances exceed a certain fixed threshold value. MCD is the current best choice for the initial estimator of a two-step procedure as it contains the robustness, equivariance, and computational efficiency properties along with its m-1/2 rate of convergence. Hence RMCD estimators are the weighted mean vector(2)x-RMCD=(∑j=1mwjxj)(∑j=1mwj),
and the weighted covariance matrix
(3)SRMCD=cα,p*dγ,αm,p*∑j=1mwj(xj-x-RMCD)(xj-x-RMCD)′∑j=1mwj,
where cα,p is the multiplication factors for consistency [11], dγ,αm,p is the finite sample correction factor [12], and the weights wj are defined as
(4)wj={1ifRD(xj)≤qα,0otherwise,
where the robust distance RD(xj)=(xj-x-MCD)′SMCD-1(xj-x-MCD) and qα is (1-α)100% quantile of the chi-square distribution with p degrees of freedom.This reweighting technique improves the efficiency of the initial MCD estimator while retaining (most of) its robustness. Hence the RMCD estimator inherits the affine equivariance, robustness, and asymptotic normality properties of the MCD estimators with an improved efficiency.


## 2.2. MVE and RMVE Estimators

Determining the MVE estimators of location and scatter parameters of the distribution is almost in line with that of the MCD estimator. As in the case of MCD, all the possible subsets of data points with sizeh=(m*γ)  (where0.5≤γ≤1) is obtained first. Then the ellipsoid of minimum volume that covers the subsets are obtained to determine the MVE estimators. The MVE location estimator is the geometrical center of the ellipsoid, and the MVE scatter estimator is the matrix defining the ellipsoid itself, multiplied by an appropriate constant to ensure consistency [13, 16]. Thus MVE estimator does not correspond to the sample mean vector and the sample covariance matrix as in the case of the MCD estimator. Here (1-γ) represents the breakdown point of the MVE estimators, as in the case of MCD, and it has the highest possible finite sample breakdown point when h=(m+p+1)/2m [8, 17]. The MVE estimator has an m-1/3 rate of convergence and a nonnormal asymptotic distribution [17].As in the case for MCD estimators, MVE estimators are also not efficient. Hence, a reweighted version similar to that for MCD has been proposed by Rousseeuw and van Zomeren [16]. Note that it has been shown more recently that the RMVE estimators do not improve on the convergence rate (and thus the 0% asymptotic efficiency) of the initial MVE estimator [8, 12]. Therefore, as an alternative, a one-step M-estimator can be calculated with the MVE estimators as the initial solution [13, 18] which results in an estimator with the standard m-1/2 convergence rate to a normal asymptotic distribution. For more details on MCD/MVE estimators one may refer to Chenouri et al. [6] or Jensen et al. [5]. The algorithm to determine the MVE/RMVE estimators is available in the statistical software SPLUS, R, SAS, and Matlab.


## 3. Robust Control Charts

We propose to useT2 charts with robust estimators of location and dispersion parameters based on RMCD/RMVE for monitoring the process mean of Phase-I multivariate individual observations. RMCD/RMVE estimators inherit the nice properties of initial MCD estimators such as affine equivariance, robustness, and asymptotic normality while achieving a higher efficiency. We now define a robust T2 control chart with RMCD and RMVE estimators for ith multivariate observation as
(5)TRMCD2(xi)=(xi-x-RMCD)′SRMCD-1(xi-x-RMCD),TRMVE2(xi)=(xi-x-RMVE)′SRMVE-1(xi-x-RMVE),
where x-RMCD, x-RMVE are the mean vectors and SRMCD, SRMVE are the dispersion matrices under the RMCD/RMVE methods based on m multivariate observations.The exact distribution ofTRMCD2/TRMVE2 estimators not available, hence the control limits for Phase-I data are obtained empirically. In the next subsection we apply Monte Carlo simulation to estimate quantiles of the distribution of TRMCD2 and TRMVE2 for several combinations of sample sizes and dimensions. For each dimension, we further introduce a method to fit a smooth nonlinear model to arrive, the control limits for any given sample size.


### 3.1. Computation of Control Limits

We performed a large number of Monte Carlo simulations to obtain the control limits. We generatedn=200,000 samples of size m from a standard multivariate normal distribution MVN(0, Ip) with dimension p. Due to the invariance of the TRMCD2 and TRMVE2 statistics, these limits will be applicable for any values of μ and Σ. Using the reweighted MCD/MVE estimators x-RMCD, SRMCD, x-RMVE, and SRMVE with a breakdown value of γ=0.50,TRMCD2/TRMVE2 statistics for each observation in the data set were calculated using (5), and the maximum value attained for each data set of size m was recorded. The empirical distribution of maximum of TRMCD2 and TRMVE2 was inverted to determine the (1-α)100% quantiles. We used the R-function “CovMcd()” in the “rrcov” package written by Torodov [19] to ascertain the RMCD/RMVE estimators.We have constructed the empirical distribution ofTR2MCD/TR2MVE as above for m=[30(1)50,55(5)100,110(10)200],p=(2,3,…10) when γ=0.50 and arrived at the control limits for α = (0.05, 0.01, and 0.001). The scatter plots of the quantiles and sample sizes for different dimensions suggest a family of nonlinear models of the form
(6)fp,α,γ,m=a1,(p,α,γ)+a2,(p,α,γ)ma3,(p,α,γ),
where  a1(p,α,γ),a2(p,α,γ), and a3(p,α,γ) are the model parameters. For clarity, the scatter plot of the actual and the fitted values of the quantiles of TRMCD2 and TRMVE2 for p=2,6, and 10 are given in Figures 1, 2, and 3; other plots are omitted to save space.Scatter plot ofTRMCD2/TRMVE2 control limits and the fitted curve for p=2.
(a)
(b)
(c)
(d)Scatter plot ofTRMCD2/TRMVE2 control limits and the fitted curve for p=6.
(a)
(b)
(c)
(d)Scatter plot ofTRMCD2/TRMVE2 control limits and the fitted curve for p=10.
(a)
(b)
(c)
(d)From Figures1, 2, and 3, we can see that the nonlinear fit is very well supported by the high R2 values, which help us to determine the TRMCD2 and TRMVE2 control limits for any given sample size. The least square estimates of the parameters a1(p,α),a2(p,α), and a3(p,α) when γ=0.50 for dimensions p=(2,3,…,10) and α = (0.05, 0.01 and 0.001) for TRMCD2/TRMCD2 control charts are given in Table 1. Using these estimates, the control limits for TRMCD2 and TRMVE2 can be found using (6) for any sample size.Table 1
Estimates of the model parametersa1(p,α),a2(p,α),a3(p,α) for TRMCD2/TRMVE2 control charts.
p
α = 0.05
α = 0.01
α = 0.001
a
^
1
a
^
2
a
^
3
a
^
1
a
^
2
a
^
3
a
^
1
a
^
2
a
^
3
T
RMCD
2
2
17.223
41102
2.647
21.134
38170
2.329
27.051
192909
2.508
3
20.134
35844
2.209
24.287
128924
2.344
31.350
1144947
2.718
4
23.152
269357
2.548
28.181
1272773
2.773
35.575
5989325
2.973
5
24.685
467949
2.524
28.437
1417059
2.632
31.013
2666196
2.593
6
26.962
1762051
2.746
29.654
3061216
2.711
31.662
5414248
2.669
7
24.892
1099128
2.493
22.882
1585224
2.416
19.058
3465278
2.444
8
27.236
2908821
2.667
27.245
4922576
2.644
28.326
12134778
2.710
9
23.974
2447649
2.534
21.420
4726835
2.554
18.772
14096595
2.676
10
31.894
12572909
2.914
37.085
34375654
3.033
56.573
172176786
3.301
T
RMVE
2
2
17.442
29553
2.494
21.365
31571
2.244
27.594
148747
2.434
3
20.286
22497
2.066
24.387
59096
2.13
31.326
338665
2.402
4
23.095
108855
2.286
27.549
291064
2.372
35.109
1255429
2.576
5
24.796
238966
2.334
28.302
508097
2.367
32.008
1063783
2.377
6
27.585
1041090
2.606
31.126
1882888
2.601
37.136
4714353
2.671
7
28.151
1541634
2.598
30.936
3183762
2.635
39.357
12199414
2.827
8
34.917
14798692
3.127
45.767
75616029
3.419
70.875
840512379
3.904
9
39.191
59094377
3.415
50.271
275604839
3.679
72.768
1960966919
4.039
10
50.733
950607720
4.099
68.154
4696452032
4.379
110.587
56398461817
4.881For the implementation of a robust control chart, first collect a sample ofm multivariate individual observations with dimension p. Compute robust estimates of mean and covariance matrix using R or any other software with γ=0.50, and determine TRMCD2/TRMVE2. Outliers can be determined by comparing the TRMCD2/TRMVE2 values with control limits obtained using (6) for specific values of α,m,p, and the constants given in Table 1. The outlier free data can be used to construct the standard T2 control chart for monitoring the Phase-II observations.


## 3.1. Computation of Control Limits

We performed a large number of Monte Carlo simulations to obtain the control limits. We generatedn=200,000 samples of size m from a standard multivariate normal distribution MVN(0, Ip) with dimension p. Due to the invariance of the TRMCD2 and TRMVE2 statistics, these limits will be applicable for any values of μ and Σ. Using the reweighted MCD/MVE estimators x-RMCD, SRMCD, x-RMVE, and SRMVE with a breakdown value of γ=0.50,TRMCD2/TRMVE2 statistics for each observation in the data set were calculated using (5), and the maximum value attained for each data set of size m was recorded. The empirical distribution of maximum of TRMCD2 and TRMVE2 was inverted to determine the (1-α)100% quantiles. We used the R-function “CovMcd()” in the “rrcov” package written by Torodov [19] to ascertain the RMCD/RMVE estimators.We have constructed the empirical distribution ofTR2MCD/TR2MVE as above for m=[30(1)50,55(5)100,110(10)200],p=(2,3,…10) when γ=0.50 and arrived at the control limits for α = (0.05, 0.01, and 0.001). The scatter plots of the quantiles and sample sizes for different dimensions suggest a family of nonlinear models of the form
(6)fp,α,γ,m=a1,(p,α,γ)+a2,(p,α,γ)ma3,(p,α,γ),
where  a1(p,α,γ),a2(p,α,γ), and a3(p,α,γ) are the model parameters. For clarity, the scatter plot of the actual and the fitted values of the quantiles of TRMCD2 and TRMVE2 for p=2,6, and 10 are given in Figures 1, 2, and 3; other plots are omitted to save space.Scatter plot ofTRMCD2/TRMVE2 control limits and the fitted curve for p=2.
(a)
(b)
(c)
(d)Scatter plot ofTRMCD2/TRMVE2 control limits and the fitted curve for p=6.
(a)
(b)
(c)
(d)Scatter plot ofTRMCD2/TRMVE2 control limits and the fitted curve for p=10.
(a)
(b)
(c)
(d)From Figures1, 2, and 3, we can see that the nonlinear fit is very well supported by the high R2 values, which help us to determine the TRMCD2 and TRMVE2 control limits for any given sample size. The least square estimates of the parameters a1(p,α),a2(p,α), and a3(p,α) when γ=0.50 for dimensions p=(2,3,…,10) and α = (0.05, 0.01 and 0.001) for TRMCD2/TRMCD2 control charts are given in Table 1. Using these estimates, the control limits for TRMCD2 and TRMVE2 can be found using (6) for any sample size.Table 1
Estimates of the model parametersa1(p,α),a2(p,α),a3(p,α) for TRMCD2/TRMVE2 control charts.
p
α = 0.05
α = 0.01
α = 0.001
a
^
1
a
^
2
a
^
3
a
^
1
a
^
2
a
^
3
a
^
1
a
^
2
a
^
3
T
RMCD
2
2
17.223
41102
2.647
21.134
38170
2.329
27.051
192909
2.508
3
20.134
35844
2.209
24.287
128924
2.344
31.350
1144947
2.718
4
23.152
269357
2.548
28.181
1272773
2.773
35.575
5989325
2.973
5
24.685
467949
2.524
28.437
1417059
2.632
31.013
2666196
2.593
6
26.962
1762051
2.746
29.654
3061216
2.711
31.662
5414248
2.669
7
24.892
1099128
2.493
22.882
1585224
2.416
19.058
3465278
2.444
8
27.236
2908821
2.667
27.245
4922576
2.644
28.326
12134778
2.710
9
23.974
2447649
2.534
21.420
4726835
2.554
18.772
14096595
2.676
10
31.894
12572909
2.914
37.085
34375654
3.033
56.573
172176786
3.301
T
RMVE
2
2
17.442
29553
2.494
21.365
31571
2.244
27.594
148747
2.434
3
20.286
22497
2.066
24.387
59096
2.13
31.326
338665
2.402
4
23.095
108855
2.286
27.549
291064
2.372
35.109
1255429
2.576
5
24.796
238966
2.334
28.302
508097
2.367
32.008
1063783
2.377
6
27.585
1041090
2.606
31.126
1882888
2.601
37.136
4714353
2.671
7
28.151
1541634
2.598
30.936
3183762
2.635
39.357
12199414
2.827
8
34.917
14798692
3.127
45.767
75616029
3.419
70.875
840512379
3.904
9
39.191
59094377
3.415
50.271
275604839
3.679
72.768
1960966919
4.039
10
50.733
950607720
4.099
68.154
4696452032
4.379
110.587
56398461817
4.881For the implementation of a robust control chart, first collect a sample ofm multivariate individual observations with dimension p. Compute robust estimates of mean and covariance matrix using R or any other software with γ=0.50, and determine TRMCD2/TRMVE2. Outliers can be determined by comparing the TRMCD2/TRMVE2 values with control limits obtained using (6) for specific values of α,m,p, and the constants given in Table 1. The outlier free data can be used to construct the standard T2 control chart for monitoring the Phase-II observations.


## 4. Performance Analysis

We assess the performance of the proposed charts when outliers are present due to the shift in the process mean. In their study, Jensen et al. [5] concluded that the TMCD2/TMVE2 control charts had better performance in terms of probability of signal. Hence, we compare the performance of our proposed method with TMCD2/TMVE2 charts as well as the standard T2 charts based on classical estimators. Our study compares more combinations of dimension p, sample size  m,  and  π. For a particular combination of p,  m,  and  π,  a number of datasets are generated. Out of the  m observations generated, m*π of them are random data points generated from the out-of-control distribution, and the remaining  m*(1-π) observations are generated from the in-control distribution so that the sample of m data points may contain some outliers. We set π=0.10 and 0.20 to ensure that the sample contains few outliers. Without loss of generality, we consider the in-control distribution as N(0,Ip). The out-of-control distribution is a multivariate normal with a small shift in the mean vector with same covariance matrix. The amount of mean shift is defined through a noncentrality parameter (δ), which is given by
(7)δ=(μ1-μ)′Σ-1(μ1-μ),
where (μ1-μ) is the shift in the mean vector. The larger the value of δ is, the more extreme the outliers are. The proportion of datasets that had at least one TRMCD2 or  TRMVE2 statistic greater than the control limit was calculated, and this proportion becomes the estimated probability of signal. We compared the performance of these charts with standard T2 charts, TMCD2, and TMVE2 charts. The standard T2 chart was included in our performance study as a reference because of its common usage.The probability of a signal for different values ofδ=(0,5,10,15,20,25,30) and for some of the values of  m=(30,50,100,150),p=(2,6,10) and π=(10%,20%) was considered in our study. Fifty thousand datasets of size m were generated for each combination of p, π, and δ, and the probability of signal was estimated for α=0.05,0.01, and 0.001. We considered various combinations of μ1, μ2, and ρ which determine δ as per (7) and found that the probability of signal is the same irrespective of the combination of μ1, μ2 and ρ. Hence we have considered μ1=μ2 and ρ=0 for various values of δ. We have presented only a selected set of plots to save space. The plots of probability of signal for α=0.05 and 0.01, p=2 and 6, and m=50 and 100 are given in Figures 4, 5, 6, and 7 for easier understanding. For dimension p=10, we used m=100 and 150, and the plots of probability of signal are given in Figures 8 and 9.Probability of signal forT2 control chart with different estimation methods for p=2,m=50.
(a)
(b)
(c)
(d)Probability of signal forT2 control chart with different estimation methods for p=2,   m=100.
(a)
(b)
(c)
(d)Probability of signal forT2 control chart with different estimation methods for p=6,   m=50.
(a)
(b)
(c)
(d)Probability of signal for forT2 control chart with different estimation methods for p=6,   m=100.
(a)
(b)
(c)
(d)Probability of signal forT2 control chart with different estimation methods for p=10,   m=100.
(a)
(b)
(c)
(d)Probability of signal forT2 control chart with different estimation methods for p=10,   m=150.
(a)
(b)
(c)
(d)From Figures4–9, we can see that when the value of the noncentrality parameter is zero or close to zero, the probability of signal is close to α which is expected for an in-control process. As the value of the noncentrality parameter increases the probability of signals also increases. Using this criterion, we select the best method for identifying the outliers. If the probability of signal does not increase for increase in noncentrality parameter, then it is clear that the estimator has broken down and is not capable of detecting the outliers.A careful examination of these plots of probability of signals corresponding to various values ofp, m, and π indicates that for small values of p and m, TRMVE2 performs well. As  m  and  p  increase, TRMCD2 chart is superior. For example, from Figures 4 and 5 we see that TRMVE2 has slight advantage over TRMCD2. But compared to TMCD2/TMVE2 charts, TRMCD2/TRMVE2 charts are performing well which is evident from all the plots presented here. When p is large (see Figures 8 and 9), the TRMCD2 has clear advantage compared to TRMVE2. From these figures, we see that standard T2 control chart possesses little ability to detect the outliers and the TMVE2, and TMVE2 stands below the TRMCD2/TRMVE2 charts throughout all the values of δ.Asp increases for a fixed value of m, the breakdown points of RMCD and RMVE get smaller as the breakdown value is given by (m-p+1)/2m. This suggests that the larger p is, the larger m will need to be in order to maintain the breakdown point, which is very well demonstrated in Figures 8 and 9. In general, there was always one estimator, RMCD or RMVE, that was found to be superior across all the values of the noncentrality parameter as long as the proportion of outliers was not so big as to cause the estimators to break down. This greatly simplifies the conclusions that can be made about when the RMCD or RMVE estimators are preferred to the MCD and MVE estimators.Nevertheless,TRMCD2 and TRMCD2 charts are preferred for the various combinations of m, p, and π, and some broad recommendations can be made on the selection among these two charts. When m<100, the TRMVE2 will be the best for small dimension. When m≥100, the TRMCD2 is preferred. As p increases, then the percentage of outliers that can be detected by the TRMVE2 chart decreases. It is true for both the charts that when p is higher, the number of outliers that can be detected decreases for smaller sample sizes. Thus for Phase-I applications where the number of outliers is unknown, TRMVE2 should be used only for smaller sample sizes, and it is also computationally feasible. TRMCD2 should be used for larger sample sizes or when it is believed that there is a large number of outliers. When the dimension is large, larger sample sizes are needed to ensure that the estimator does not break down and lose its ability to detect outliers. Hence for larger dimension cases, TRMCD2 is preferred with large sample sizes. For very small samples (m<30), one may opt for higher values of γ, for which control limits need to be developed.


## 5. Case Example

To illustrate the applicability of the proposed control chart method, we discuss a real case example taken from an electronic industry. The data gives 105 measurements of 3 axial components of acceleration measured by accelerometer on a e-compass unit fixed on the objects. The mean vector and covariance matrix under the classical, RMCD, and RMVE methods of the sample data considered are given by(8)X-=(6.31435.73395.7527),S=(3.4022-1.1524-1.0746-1.15241.92491.1209-1.07461.12092.2004),X-RMCD=(5.71256.26436.0837),SRMCD=(2.85490.1901-0.29260.19011.11750.4433-0.29260.44332.3115),X-RMVE=(5.77906.13795.9894),SRMVE=(2.7837-0.2642-0.5745-0.26421.27780.6080-0.57450.60802.2909).
A simple comparison of these estimators indicates that there are outliers in the Phase-I data. The plots of T2,TRMCD2, and TRMVE2 values along with the respective control limits at 99% confidence level for the sample data are given in Figure 10.Figure 10
T
2
,
T
RMCD
2
,
and
T
RMVE
2 control charts for the sample data.The control limits forT2 are arrived at based on beta distribution, and TRMCD2/TRMVE2 are calculated using (6) for p=3 and m=105. From Figure 10, it is very clear that both TRMCD2 and TRMVE2 control chart alarms signal for 3 outliers whereas the standard T2 control chart alarm signals for none even though all the charts are having the same pattern. This indicates the effectiveness of the proposed robust control charts in identifying the outliers.


## 6. Conclusions

Use of robust control chart in Phase-I monitoring is very important to assess the performance of the process as well as detecting outliers. We proposeTRMCD2/TRMVE2 control charts for Phase-I monitoring of multivariate individual observations. The control limits for these charts are arrived empirically and a non-linear regression model is used for arriving control limits for any sample size. The performance of the proposed charts were compared under various data scenarios using large number of Monte Carlo simulations. Our simulation studies indicate that TRMVE2 control charts are performing well for smaller sample sizes and smaller dimension where as TRMCD2 control charts are performing well for larger sample sizes and larger dimensions. We illustrated our proposed robust control chart methodology using a case study from the electronic industry.



---

*Source: 542305-2013-06-09.xml*