INTRODUCTION
One of the most popular tools used to construct a multivariate control chart using individual observations is the Hotelling’s T^{2} control chart. The main aim of a multivariate control chart is detecting the cause of the process change. In the literature, construction of this control chart is carried out in two distinct phases. In phase I operation, Hotelling's T^{2} chart is often used to purge outliers whereas signal detection is done in the phase II operation. In the first phase, Historical Data Set (HDS) is used to establish the control limits in order to detect outliers. If one or more of these initial observations in HDS is out of control, the phase I control limits are recomputed based on the remaining observations. This step is repeated until the process comes to the state of control and the incontrol parameters of the process are estimated. Then clean, phase I data is used to set control limits for monitoring future (phase II) data. Generally, for a retrospective phase I analysis of a Historical Data Set (HDS) our objective is defined in two steps.
In the first step, mean shifts that might result to ruins the estimation of the incontrol mean vector and variancecovariance matrix are discovered. In the second step, after identification and elimination of multivariate outliers, the incontrol parameters estimates which are established in the control subset of the HDS, to be used.
In the computation of the Hotelling's T^{2} the classical mean and
covariance are utilized, hence the results can be heavily influenced by outliers.
In recent years, several multivariate statistical process control techniques
have been proposed to analyze and monitor multivariate data (Tracy
et al., 1992; Sullivan and Woodall, 1996;
Vargas, 2003). In this study, in order to decrease the
impact of outliers we propose the use of the reweighted robust estimator instead
of the classical estimator in the computation of Hotelling's T^{2}.
HOTELLING’S T^{2}
Hotelling et al. (1974) introduced a statistic
that combines information from the dispersion and mean of several variables.
The T^{2} statistic may be computed using a single observation from
p components, or it may be computed using the mean from a sample of size n.
In this study, a subgroup of size 1 (i.e., a single observation) will be assumed
for the T^{2} computations. The T^{2} statistic for a px1 multivariate
normal vector, X = (X_{1},..., X_{P}) is defined as:
Where:
and
are the common estimators of the mean vector and covariance matrix obtained from the historical data set. Usually, we assume that the X_{i}’s are independent multivariate normal, MVN_{p}(μ, Σ) with μ and covariance matrix Σ. In order to detect any possible signal, for each individual observation, we compare the T^{2} with control limits.
Hotelling’s T^{2} statistic can be described by different distribution based on different situations. Suppose the parameters of the MVN distribution, μ and Σ are known. The T^{2} statistic for an individual observation vector X follows the chisquare distribution:
Assuming the parameters of the MVN distribution are unknown, the estimators
and
S must be used. These values are calculated from HDS, consisting of n observations.
The distribution of the T^{2} statistic for an individual observation
vector X, independent of
and S is:
where, F is distributed as F distribution with p and (np) degrees of freedom.
In another case, let us say the observation vector X is not independent of the
estimators
and S but is included in their computation. Then the distribution of the T^{2}
statistic is given as follows (Tracy et al., 1992):
where, B is distributed as Beta distribution with parameters p/2 and (np1)/2.
In the computation of Hotelling's T^{2} the classical estimators of
location and dispersion have been used to estimate the population mean vector
and covariance matrix. However, this classical estimators are affected by outliers
in the phase I data. Therefore, we construct an alternative estimator based
on Reweighted Minimum Covariance Determinant (RMCD) and Reweighted Minimum
Volume Ellipsoid Estimators (RMVE) which have a higher efficiency and reduce
the influence of outlying observations. Vargas (2003)
and Jensen et al. (2007) recommended using the
Minimum Covariance Determinants (MCD) and Minimum Volume Ellipsoid (MVE) estimators
of mean vector and covariance matrix in the Hotelling's T^{2} charts
(Wisnowski et al., 2002; Chenouri
et al., 2009). The exact distribution of the T^{2} based
on the RMCD and RMVE are not tractable, so we used Monte Carlo method to obtain
the appropriate control limits.
ROBUST ESTIMATORS
Modified version of the MCD and MVE estimators are RMCD and RMVE. Reweighting step that is included in the computation of RMCD and RMVE, greatly increase the finitesample efficiency of these estimators.
RMCD estimators inherit the good properties of MCD estimators such as affine
equivariance, the maximal asymptotic breakdown point and asymptotic normality.
In addition, a fast and efficient approximate algorithm to compute Reweighted
Minimum Covariance Determinant (RMCD) is available. The above mentioned estimators,
have been studied by several authors (Rousseeuw and Van
Zomeren, 1990; Lopuhaa and Rousseeuw, 1991; Willems
et al., 2002). In contrast, the reweighted MVE estimators of mean
vector and covariance matrix with ordinary MVE, reweighted MVE estimators are
more efficient, akin to the reweighted MCDs. So, we expect that the new robust
control charts (reweighted MVE) outperforms those based on the ordinary MVEs.
RMCD ESTIMATOR
Rousseeuw (1985) introduced the Minimum Covariance
Determinant (MCD) estimator that has finite sample and asymptotic breakdown
points 1/2, which is based on the subset of h = ⌊nγ⌋ (where
0.5<γ<1) data points whose covariance matrix has the smallest possible
determinant, where 1γ is the asymptotic breakdown point. The MCD location
estimate
is defined as the average of these h data points and the MCD scatter estimate
is given by S_{MCD} = aC_{MCD} where, C_{MCD} is the
covariance matrix of the subset of h points and a is the multiple of the consistency
and the finite sample correction factors (Willems et
al., 2002; Chenouri et al., 2009). The
MCD estimators have the highest finite sample breakdown point when h = ⌊(n+p+1)/2⌋
(Rousseeuw and Leroy, 1987). In order to estimate the
RMCD estimator, for each individual observation we computed the weight based
on the robust distance:
The weight (w_{t}) equals 1 when the squared robust distance D(X_{i})^{2}
is smaller than the cutoff value χ^{2}_{ρ,η}
(we use the value η = 0.975 which was advocated and used by Rousseeuw
and Van Driessen (1999). Otherwise the weight equals zero. Afterward we
can compute the and
S_{RMCD}:
Here, c_{η,p} = η/p (χ^{2}_{(p+2)}≤q_{η}),
makes S_{RMCD} consistent under the multivariate normal distribution
and yields more reliable outlier identification. The factor
is a finite sample correction given by Pison et al.
(2002). The most commonly used algorithm for computational purposes is FASTMCD
algorithm (Rousseeuw and Van Driessen, 1999) which has
been used in this study.
RMVE ESTIMATOR
In addition to the MCD, Rousseeuw (1985) introduced
the minimum volume ellipsoid (MVE) as affine equivariant and high breakdown
estimators of location and dispersion. The MVE estimators are based on the smallest
volume ellipsoid that covers at least h = ⌊nγ⌋ (where 0.5<γ<1)
observations. The MVE location estimator t∈R_{p} and dispersion
estimator C minimize the determinant of positive definite symmetric matrix C
of size p, subject to the condition:
The constant c determines the magnitude of C and is usually chosen to make
C, a consistent estimator of the covariance matrix under multivariate normal
model, i.e., .
Davies (1992) has shown that the MVE location estimator
has a slow n^{1/3} rate of convergence and a nonnormal asymptotic
distribution. This low rate of convergence implies that the asymptotic efficiency
of the MVE estimators is 0%. Therefore, to increase the efficiency of MVE, we
propose to employ the reweighted MVE (RMVE) similarly with the reweighted
MCD. For more details on MVE and RMVE estimators refer to the Van
Aelst and Rousseeuw (2009).
DESIGNING THE ROBUST T^{2} CONTROL LIMITS
As mentioned before, in phase I analysis the usual control charts use standard mean and covariance matrix estimators, which are sensitive to outliers. It is known that the Mahalanobis distance based on the classical estimators is effective in detection of a single outlier, but is not appropriate in the case of multiple outliers.
This is due to the fact that, outliers mask each other and estimator fails to detect them, which is described as the masking effect. In the literature, a variety of robust estimators have been proposed to overcome this problem and minimize the impact of outliers. To this end, we constructed a new statistic by substituting the classical estimators in the Hotelling's T^{2} by the reweighted MCD and reweighted MVE estimators.
Consider that the phase I historical data set consists of n timeordered vectors that are independent of each other. Each vector is of dimension p and there are no subsamples in the observations, the reweighted MCD based Hotelling T^{2} statistic for phase II observation X_{f}∉{X_{1},..., X_{n}} is define as follows:
The finitesample distributions of the T^{2}_{RMCD} (f) is
unknown, therefore, we computed the control limits based on the empirical distributions
of respective robust T^{2} charts. Several studies have been investigated
the asymptotic properties of these estimators (Lopuhaa and
Rousseeuw, 1991; Butler et al., 1993; Croux
and Haesbroec, 1999). It should be noted that the
and S_{RMCD} are a good approximation to the parameters μ and Σ
and X_{f}~MVN_{p} (μ, Σ). In this case, if we use
the Slutsky theorem, as n→4 the asymptotic distribution of robust T^{2}
is χ^{2}_{p}. This asymptotic distribution is only applicable
when n is large. In the case of small sample size, we apply Monte Carlo simulations
to estimate the quantiles of the T^{2}_{RMCD}(f), for several
combinations of sample sizes and dimensions. The robust T^{2} statistic
based on the reweighted MVE estimators of location and dispersion is defined
as:
where,
and S_{RMVE} are the reweighted MVE estimates of mean vector and covariance
matrix of phase I process, respectively. Since, the distribution or even asymptotic
distribution of T^{2}_{RMCD}(f) is not known, so, we obtain
the quantiles of the T^{2}_{RMVE}(f) for different combination
of dimensions and sample sizes.
SIMULATION
We generated 5000 samples of size n from a pmultivariate standard normal distribution MVN_{p}(0, I_{p}). Our simulation process in phase I and phase II was constructed as follows:
Phase I: We computed the reweighted MCD and MVE mean vector and covariance
matrix estimates (
and S_{RMCD}, S_{RMVE}) for each data set of size n.
Phase II: In addition, as phase II observation, for each data set we
randomly generated a new observation X_{f} from MVN_{p}(0, I_{p})
and calculated the respective T^{2}_{RMCD} and T^{2}_{RMVE}
value as given by Eq. 9, 10. Scatter plots
of the empirical 99% quantiles of T^{2}_{RMCD} and T^{2}_{RMVE}
versus the sample size n for different dimensions are presented in Fig.
1 and 2.
Since, the creation of control limits for each sample size is weariful and time consuming, it would be preferable to have a formula for the calculating of control limits.
Then we sketched the graph of the simulated quantiles for the T^{2}_{RMCD}(f) and T^{2}_{RMVE}(f) estimators versus the size of the data set n. A regression curve was fitted to smoothly predict the control limits for any phase I sample size n, dimensions p = 2, 6, 10 and confidence level 1α = 0.99. With respect to the graph, we utilized the regression equation f(n) = β_{1}+β_{2}/nβ_{3} for the quantiles. As mentioned earlier, since the T^{2}_{RMCD}(f) is asymptotically distributed as χ^{2}_{p}, it is reasonable to use the χ^{2}_{(p,1α)} instead of the β_{1,p,1α,γ} in Eq. 11:
where, χ^{2}_{(p,1α)} is the 1α quantile of the χ^{2}_{p} distribution with p degrees of freedom and β_{2,p,1α,γ} and β_{3,p,1α,γ} are constants. Three parameter curves give better fit for control limits of RMVE estimators based on T^{2} charts.
Table 1: 
The least square estimates of the regression parameters β_{1,p,1α,γ}
and β_{2,p,1α,γ} and β_{3,p,1α,γ}
for p = 2, 6, 10 confidence levels 1" = 0.99 and breakdown points
( = 0.5, 0.75 for RMCD estimators 


Fig. 1: 
The 99% simulated quantiles of T^{2}_{RMCD}
and the fitted curves for p = 2; 6; 10, γ = 0.5 (upper panel) and γ
= 0.75 (lower panel) 

Fig. 2: 
The 99% simulated quantiles of T^{2}_{RMVE}
and the fitted curves for p = 2; 6; 10, γ = 0.5 (upper panel) and γ
= 0.75 (lower panel) 
Table 2: 
The least square estimates of the regression parameters β_{1,p,1α,γ}
and β_{2,p,1α,γ} and β_{3,p,1α,γ}
for p = 2, 6, 10 confidence levels 1α = 0.99 and breakdown points
γ = 0.5, 0.75 RMVE estimators 

Non linear least squares method was applied to estimate the parameters β_{1,p,1α,γ},
β_{2,p,1α,γ} and β_{3,p,1α,γ}
which is shown in Table 1 and 2. Finally,
in order to compute the robust control limits for p = 2, 6, 10 with varying
sample size n, Table 1, 2 and Eq.
9 and 10 were used.
PERFORMANCE COMPARISON
Here, we designed a simulation study to assess the performance of our proposed
methods. We considered different number of variables (p = 2, 6, 10) and the
number of observations (n = 50, 150), breakdown point (γ = 0.5, 0.75),
the proportion of outliers (π = 0.10, 0.20) and the amount of the shift
in the process mean (δ = 3, 5). Our simulation studies include a nooutlier
pattern and a pattern with multiple outliers. We can decide about the performance
of the control chart by computing the probability of changes detection based
on the Phase II data. The efficiency of control chart is determined by probability
of signal that is the proportion of the T^{2}_{RMCD }and T^{2}_{RMVE}
that is located over the control limit, using 1500 replications. We let δ^{2}
= (μμ_{A})'Σ^{1} (μμ_{A})
be the non centrality parameter that measures the severity of a shift of the
out of control mean vector μ_{A} compared to the in control mean
vector μ. We can without loss of generality, use a zero mean vector as
μ and the identity covariance matrix as Σ, due to affine equivariance
property. To simulate changes in the process mean, the following cases are considered:
• 
In phase I we generated clean data set (no outlier) from the
standard multivariate normal distribution MVN_{p}(0, I_{p}),
π of them are random data points generated from the outofcontrol
distribution (MVN_{p}(μ_{I}, I_{p})) and the
other 1π observations were generated from the incontrol distribution
(MVN_{p}(0, I_{p})) 
• 
Phase II data were generated from MVN_{p}(0, μ_{II})
where, δ^{2}_{II} = μ_{II}^{2}.
To evaluate the performance of various methods in each simulation, we generated
a sample of size n and computed the T^{2}_{RMCD}(f) and
T^{2}_{RMVE}(f) for breakdown points of γ = 0.5, 0.75
from Eq. 9 and 10, for each observation
in phase II data. Then, we compared the T^{2}_{RMCD}(f)
and T^{2}_{RMVE}(f) to the approximate control limits 
The purpose of this study is to compare the proposed method with other techniques
such as standard T^{2} chart, robust T^{2}_{MCD} and
T^{2}_{MVE} estimators discussed in Vargas
(2003) and Jensen et al. (2007).
Table 3: 
Probability of signal in phase II, when phase I data is outlier
free, for p =2, 6, 10 and sample size n = 50, 150 with different values
of shift in phase II process mean vector μ_{II} 

Table 4: 
Probability of signal in phase II, when there is a slight
shift in phase I process mean vector (μ_{I} = 5), for p = 2,
6, 10 and sample size n = 50 with different values of shift in phase II
process mean vector (μ_{II}) 

The probability of signal on the same data sets is computed. The MCD and MVE
techniques are applied in phase I to detect outliers. Then we made the standard
T^{2} chart based on the clean phase I data set and compared the corresponding
T^{2} with an appropriate quantile of F distribution to monitor phase
II data. The probability of signal for detecting outliers for different phase
II process shifts in the mean vector are depicted in Table 37
δ_{I} = 0 (there was no outliers in the phase I): In all
cases, from Table 37, it is visible, by
increasing the noncentrality parameter the probability of signal increases
as well.
In the case of small sample sizes, we observed that the best estimator was
the standard T^{2}, which is also supported by the previous studies
(Wisnowski et al., 2002; Vargas,
2003; Jensen et al., 2007; Chenouri
et al., 2009). Conversely, as the sample size increased to 150, the
performance of all robust control charts are similar to the standard T^{2}
chart and the performance was satisfactory.
δ_{I} = 5 (small proportion of outlier in phase I): As shown in Table 4 and 5, when there was a slight shift in the phase I process mean vector, for small sample sizes T^{2}_{RMVE} slightly outperformed the T^{2}_{RMCD} and they perform better than the classical T^{2} and ordinary MCD and MVE. On the other hand, as the sample size increased, the performance of T^{2}_{RMCD} is slightly more superior to T^{2}_{RMVE} and other methods.
δ_{I} = 30 (large proportion of outlier in phase I): Based
on Table 6 and 7, when the noncentrality
parameter in phase I was large (δ_{I} = 30), the Reweighted robust
control charts for the breakdown point of γ = 0.5 and 0.75 surpassed the
MCD, MVE and the classical T^{2} charts. In the case of small sample
size, the performance of the T^{2}_{RMVE} based chart was the
best. Conversely as the sample size increased T^{2}_{RMCD} slightly
outperformed the T^{2}_{RMVE}.
Table 5: 
Probability of signal in phase II, when there is a slight
shift in phase I process mean vector (μ_{I} = 5), for p = 2,
6, 10 and sample size n =150 with different values of shift in phase II
process mean vector (μ_{II}) 

Table 6: 
Probability of signal in phase II, when there is a large amount
of shift in phase I process mean vector (μ_{I} = 30), for p
= 2, 6, 10 and sample size n =50 with different values of shift in phase
II process mean vector (μ_{II}) 

Table 7: 
Probability of signal in phase II, when there is a large amount
of shift in phase I process mean vector (μ_{I} = 30), for p
= 2, 6, 10 and sample size n =150 with different values of shift in phase
II process mean vector μ_{II} 

The T^{2} did not work well, even for small sample size. In order to clarify the performance of the reweighted robust control chart for different breakdown points of γ = 0.5 and 0.75 we classify the result based on the sample sizes.
Small sample size: Reweighted robust control chart with γ = 0.75 worked better than γ = 0.5 in high dimensions.
Large sample size: For the small proportion of outlier, the Reweighted robust control chart for both amount of γ = 0.5 and 0.75 worked almost similarly. However, if the percentage of outlier in data set was increased the Reweighted robust control chart with γ = 0.5 outperformed the other methods. It is worthwhile noting that, when the phase I sample contained higher proportion of outliers, higher value of breakdown point was preferred. This is the main reason for the failure of robust control charts with γ = 0.75 when there are large numbers of outliers in phase I with large δ^{2}.
CONCLUSIONS
Standard control charts are widely used in industry to detect the special causes
of variation. The common outofcontrol status is the occurrence of several
outliers in the process. It is well known that the usual parameters estimations
are sensitive to the presence of outliers, so the T^{2} chart based
on these estimators performs poorly. Our advice in this study is replacing the
classical mean vector and covariance matrix of the data in the Hotelling's T^{2}
statistic by the Reweighted MCD and Reweighted MVE estimators. These estimators
have many advantages, like affine equivariance and better efficiency than the
ordinary MCD and MVE estimators used in Vargas (2003),
Hardin and Rock (2004) and Jensen
et al. (2007).
We also recommend generating the T^{2}_{RMCD} and T^{2}_{RMVE} control limit for combination of different sample size and dimension via Monte Carlo simulation. Our simulation studies showed that when the process was incontrol and the sample size was small, the best estimator was the standard T^{2}, as noted in the literatures. Nonetheless for large sample size the T^{2}_{RMCD} and the T^{2}_{RMVE} performed similar to the classical T^{2}, chart. On the other hand, when there was outlier in phase I, T^{2}_{RMCD} and T^{2}_{RMVE} were more effective than the standard T^{2} and ordinary MCD and MVE charts.
In summary, we suggest that when phase I sample has outliers, the RMVE chart
is suitable for small sample size and the RMCD chart is the best in the case
of large sample size. Moreover, it is better to use γ = 0.5 with a large
sample size (at least 10 times greater than p) to ensure better and consistent
performance.
ACKNOWLEDGMENTS
We are grateful to Dr. Chenouri at University of Waterloo, Canada for his constructive comments.