Abstract: A control chart for detecting shifts in the variance of a process is developed for the case where the nominal value of variance is unknown. The Shewhart S control chart is one of the most extensively used statistical process control techniques developed for control process variability based on essential assumption that the underlying distribution of the quality characteristic is normal. However, this chart is very sensitive to the occurrence of occasional outliers. As an alternative, robust control charts are put forward when the underlying normality assumption is not met. In this study, a robust control chart for the process standard deviation σ by means of a robust scale estimator is proposed. The presented robust method seems to yield a better performance than the Shewhart method and had good properties for the contaminated and heavy-tailed distribution for moderate sample sizes. The proposed robust Modified biweight A chart acts as an alternative for practitioners who are interested in the detection of permanent shifts in the process variance whereby the presence of occasional outliers is usually associated with the occurrence of common causes.
INTRODUCTION
The most commonly used charting framework is to set up means and variances
charts in order to monitor the possible shift in location and scale (Hawkins
and Zamba, 2005). In this circumstance, a model of independent process data
collection from underlying distribution is assumed to be operating under an
adequate working model. An ordinary approach for controlling manufacturing process
is to take measurements on subsamples (subgroups) of items at regular intervals
and to analyze the results using
The crucial aspects of each chart are control limits which are horizontal lines used to constitute a simple decision criterion. In the event that the subsample statistic falls outside these limits, then the process control is deemed to be out-of-control, otherwise if the values fall within the control limits, it is deemed to be in-control (Montgomery, 2005).
Typically when the process is operating under optimal conditions, the measurements
taken on by the process are typically modeled as independent draws from a normal
distribution with mean μ and standard deviation σ. Subsequently,
However, it is now evident that both S and R charts are easily affected by outliers. An outlier is the observation far from the rest of the data defined by Staudte and Sheather (1990). While Williams et al. (2002) and Liu et al. (2004) defined the detected outliers as candidates for aberrant data that may otherwise adversely lead to model misspecification, biased parameter estimation and incorrect results.
The robust approach is incorporated in this study for constructing the control chart to provide methods that offer reliable parameter estimates. Further details regarding some inferences about robust estimation can be found in Fraiman et al. (2001), Hoaglin et al. (2000), Salibian-Barrera and Zamar (2004), Riazoshams et al. (2010), Midi and Hua (2009) and Willems and Van Aelst (2005). In the meaning of informal data-oriented characterization, the robust method should fit the bulk of the data well in such a way that when the data contain no outliers, it gives approximately similar results as the classical method. As opposed to that when a small amount of outliers occur they should perform more satisfactory than the classical method. Rocke (1992), Quesenberry (1993), Chakraborti (2000), Jensen et al. (2006) and Al-Nasser and Al-Rawwash (2007) studied the performance of the usual control charts from the effect of the estimation of the control limits in non-normality and with the presence of outliers. Further details regarding some inferences about the performance of robust multivariate statistical control charts can be found in Shabbak et al. (2011) and Mohammadi et al. (2011).
In this study, a robust process standard deviation σ control chart is proposed to remedy the problem of outliers in a process data. The proposed chart incorporated robust scale estimator to detect the permanent shift(s) in process variance.
SAMPLE MEDIAN AND MEDIAN ABSOLUTE DEVIATION
Sample median is probably the most resistant estimator that has widespread used in process control charts (Abu-Shawiesh, 2008). It is insusceptible to behavior in the tails of the distribution. The efficiency of the sample median drops off rapidly towards its asymptotic value of 0.64 as sample size increases under the normal distribution. It is defined as the middle value when a set of data values are organized in order of magnitude and denoted by MD. The sample median for a random sample of size n observations X1, X2,..., Xn can be defined as follows:
| (1) |
The advantages of using the sample median MD, are that, it is easy to compute, not easily affected by outliers and has a maximal of 50% breakdown point. In this situation, it appears to be an appropriate general purpose estimator and typically considered as an alternative choice to the sample mean especially when data are contaminated. While the median absolute deviation, denoted as MAD is an example of robust scale estimator which is more robust compared to the classical standard deviation, S (Hampel, 1974).
PROCESS STEP CHANGE MODEL
To illustrate the model, assuming that the process is initially in control
and the observations come from a normal distribution with a known mean of μ0
and a known variance of
By letting Si be the standard deviation calculated for the ith subgroup, where
| (2) |
Xij is the jth observation in subgroup i and
ROBUST MODIFIED BIWEIGHT A SCALE ESTIMATOR
It is a common practice to employ the S chart to monitor process variability because of tradition and ease of computation. Unfortunately such chart relies on the estimate of a standard deviation of a process which is highly biased in the present of outliers. This is also due to that the estimate has traditionally been computed from average standard deviation of 20-40 subsamples. To remedy these problems, Tatum (1997) proposed a procedure in which all the subsamples will initially be centered on its own median and then employed a biweight A estimator to the pooled residuals
| (3) |
where
| (4) |
| (5) |
and c is a tuning constant, k = 1, 2, 3,..., g, m’ = ng when n is even and m’ = (n–1)g when n is odd. It is commenced by removing the median value of each subsample from that corresponding subsample.
When n is an even number, this will yield a median-centered kth subsample given by
| (6) |
where,
| (7) |
The purpose of subtracting the subsample medians as a preface for measuring variability are such that the variability is measured within the subsamples and not between subsamples. By doing so, the centering step is intended to be protected against outliers. Indeed, this protection must be compensated by reduction in efficiency for the i.i.d. normal setting. Hence, this modification of the biweight A scale estimator in Eq. 3 will forgo some efficiency under the i.i.d. normal setting both by individually centering each subsample and by centering with the subsample median instead of the subsample mean. The process standard deviation σ is then estimated as
| (8) |
where, kn,g,c is a factor chosen in such a way that
B7 seems to perform very well in simulation studies when operating under diffuse disturbances but in the occurrence of localized disturbance, nonetheless, the relative performance of B7 suffers as it does not utilize the information on the location of the disturbed subsample(s). To rectify this problem, a term hk is incorporated in the definition of uk, j so as to downweight subsamples for which the IQR defined in Eq. 7 is unusually large
| (9) |
where
| (10) |
and
| (11) |
The size of IQR in the kth subsample relative to M is evaluated by Ek in which M is the MAD of subsample-centered observations. When Ek is uncommonly large, the hk will downweight all the observations in the kth subsample by growing the values of uk, j. The downweighting initiates when Ek>4.5 and continue until Ek>7.5, at which point the observations in the subsample will be assigned the weight 0. A simulation study is carried out to determine the cutoff values of the distribution Ek in the i.i.d. normal setting.
Now, let denote
| (12) |
where dn, g, c again is a factor chosen in such a way that
In this study, an alternative scale estimator which is based on the Modified biweight A estimator is proposed. Instead of considering the pooled residuals given in Eq. 3, we propose to compute subsample scale estimate on the basis of the corresponding subsample defined as
| (13) |
with n is subsample size, comparing to Eq. 3, in fact m' = n when n is even and m' = (n–1) when n is odd where
| (14) |
and
| (15) |
In general, the average of the g subsamples scale estimates is defined as
| (16) |
Consequently, the process standard deviation σ is then estimated by
| (17) |
| Table 1: | Simulated constant values c4 for different scale estimators |
Let us now denote the proposed
CONTROL LIMITS
In this section we have incorporated the general equations for constructing control limits (Kocherlakota and Kocherlakota, 1995). For illustrative purpose, in order to monitor the variability, for instance, the S chart using the standard deviation, the scale parameter is estimated by averaging sample standard deviation over the m subgroups
| (18) |
Generally, the control limits for a scale chart is defined as
| (19) |
To present an algorithm for setting up a control chart for the proposed robust
scale estimator we first need to determine the constant c4 appropriate
to the particular estimator. This constant was obtained using computer simulations.
For illustration, let a sample of size n taken from N (0,1) and then corresponding
scale estimator value is calculated. Here, for instance, if we use sample standard
deviation S as our scale estimator, the constant c4 is determined
by averaging over 100,000 repetitions of these estimated values. This will give
E(S/c4) = σ for any σ. Table 1 exhibits
the constant c4 for the estimators in present study such that E(scale
estimator/c4) = σ. It can be seen that if scale estimator is
the sample standard deviation (E1), the simulated values agree closely
with the standard tabled values (Montgomery, 2005). The
estimates E2 and E3 in Table 1 represent
the determined constant values for the MAD chart and the proposed robust Modified
biweight A chart. Hence, the corresponding control limits for different charts
can be constructed by applying Eq. 19 with the desired
Similarly, the
To summarize, there are three scale estimators under study:
| • | E1: Standard Deviation (S) |
| • | E2: Median Absolute Deviation (MAD) |
| • | E3: Proposed robust Modified biweight A (D7) |
CHANGE POINT ESTIMATOR
We adopt change point estimator which is the maximum likelihood estimator (MLE) of the time of the change in the normal process variance (Samuel et al., 1998). The derivation of the estimator can be found in the appendix. In the following we illustrate the procedure of estimating the time of step change when the process variance does occur. Once the corresponding control chart (S, MAD or the proposed robust Modified biweight A chart) signals, the change point estimator is then applied to the data to estimate the time of the change at which one must search for the value of i in the range 1≤i≤T which maximizes
| (20) |
With the reverse cumulative sum of squares,
is the sum of the T–i+1 most recent subgroups sums of squares. The subgroup corresponding to the largest Di value which will indicate the first subgroup from the changed process.
SIMULATION STUDY
We will now analyze the performance of the change point estimator adopting
the three different control charts through a simulation study. Assuming that
the process is initially in-control and the observations come from a normal
distribution with a known mean of μ0 and a known variance of
Four types of disturbances were employed in this simulation study:
| • | I 95-5 mixture of symmetric variance disturbances model in which each observation xk, i has a 95% probability of being drawn from N (0,1) distribution and a 5% probability of being drawn from N (0,a) distribution where a = 2.0, 3.0,...,7.0 |
| • | II Asymmetric variance disturbances model in which each observation is drawn from the N (0,1) distribution and has a 5% probability of having a multiple of a χ2(1) random variable added to it (the multiplier λ takes values of 0.01, 0.05 and 0.10) |
| • | III Localized model, symmetric variance disturbances in which one subsample are drawn from the N (0,a) distribution where a = 2.0, 3.0,..., 7.0 and the remaining subsamples contain observations drawn from the N (0,1) distribution |
| • | IV Mixture model of diffuse mean disturbances in which each observation has a 95% probability and a 5% probability of being drawn from N (b,1) distribution with b = 1.0, 2.0,…, 6.0 |
Summarize for the steps of simulation:
| • | For each type of contaminated models discussed above, data consisting of m = 30 subgroups of size n = 5 were generated. Then, the corresponding control limits for the three charts (S, MAD and the proposed robust Modified biweight A chart) under investigation were constructed |
| • | Once the corresponding control chart signals, the change point estimator was then applied to the data to estimate the time of the change at which one must look for the value of i in the range 1≤i≤T which maximizes Eq. 20. In what follows, the change point estimate will be evaluated. The procedure was repeated for a total of 10,000 times for each of the contaminated models under studied |
| • | Soon after the end of the simulations, the average of the change point
estimate |
The three charts are evaluated based on expected length. In fact, the expected
length E(T) is the expected time at which the control chart signals a change
in the process variance which is actually occurred following a chosen fixed
subgroup of 100. It is generally noticed that Shewhart’s S control chart
occasionally issues a signal of a change in a process variance monitoring a
considerable amount of time after the change in the process variance that really
occurred. Estimating the time of process change with the time that a control
chart issues a signal would lead to a badly biased and, therefore, probably
misleading estimate of the time of the process change. The bias could be stemmed
from potentially large delay in generating a signal from the control chart.
If the time of change can be ascertained correctly, process engineers would
have a smaller search window within which to detect the special cause. Hence,
it is served as another criterion of assessing the performance of a control
chart which measures how quick the chart would signal. In this sense we would
prefer to have an expected length which is as close as possible to the actual
change point of 100. To further illustrate, for instance, the expected length
of 107.83 implies that the corresponding chart required about 8 points before
issuing a signal. The performance of the estimators in terms of average change
point estimates
Type I: Symmetric variance disturbance: The data under symmetric variance
disturbances were considered in which average change point estimates
| Table 2: | Expected lengths and average of change point estimates for type of disturbance I with magnitude of shift in variance of δ = 2.0 |
| Table 3: | Expected lengths and average of change point estimates for type of disturbance II with magnitude of shift in variance of δ = 2.0 |
For instance with the magnitude of a = 5.0, the expected length using S chart was 107.83 and the MAD chart was 112.32 while the proposed robust Modified biweight A scale chart (Proposed Scale Chart, D7) had an expected length of 103.15, it just needs about three points to issue a signal. Process engineers monitoring the process should view this as an important aspect due to the fact that they would have a smaller search window within which to look for the special cause.
Type II: Asymmetric variance disturbance: Now we turn to the control charts constructed when the data were generated from asymmetric variance disturbance model. Table 3 gives summary of the performance of average change point estimates and the expected lengths from the simulation runs for different magnitudes of λ values which include 0.01, 0.05 and 0.10. On the whole, the estimator manages to estimate the process change point correctly for all charts under investigation which essentially should be close to 100. Even so, the performance among the charts in terms of expected length is quite apparent. As magnitude of λ increases, the S chart deteriorates quickly. This can be noticed from Table 3 that with λ = 0.05 and δ = 2.0, S chart demands about 267.66 points before issuing a signal. By contrast, both the MAD and the proposed robust Modified biweight A charts require approximately 9.95 points and 2.82 points, respectively. The overall performance is remarkably consistent for each level magnitude of shifts using the proposed robust Modified biweight A chart.
| Table 4: | Expected lengths and average of change point estimates for type of disturbance III with magnitude of shift in variance of δ = 2.0 |
| Table 5: | Expected lengths and average of change point estimates for type of disturbance IV with magnitude of shift in variance of δ = 2.0 |
Type III: Localized variance disturbance: Again, for the process step change of standardized magnitude δ = 2.0, a comparable results can be observed from Table 4 in terms of the average estimated time of the process change. The estimates seem to scatter around closely to the actual change point of 100. However, while considering the performance of expected length, the MAD chart appears to perform poorly compared to the other two charts. For instance, it can be seen that with magnitude a = 6.0, the expected length by means of the MAD chart is 126.00. In the meantime, S chart and the proposed robust Modified biweight A chart both required expected length of about 104.00 and 103.00, respectively. The relative performances become apparent with the increases of magnitude of shift a. On the whole, by comparison, the proposed robust Modified biweight A chart appears to exhibit slightly better in overall performance in series of magnitudes of a under investigation.
Type IV: Mean shift disturbance: Lastly we focus on the performance of the three charts on the occasion that mean shift disturbance does occur in the preliminary data for the control charting process. As can be expected, on average, each chart is able to detect the position of the process change considerably near to the actual time of the change. The proposed robust Modified biweight A chart gives the best result among all the charts examined associated with expected length while the MAD chart is the worst as magnitude of b increases. For a magnitude b of 5.0, for instance, the proposed robust Modified biweight A chart outperforms the other two charts. In general, it needs approximately 3.87 points to indicate a signal. On the other hand, both S and MAD charts on average attain expected length of 109.75 and 121.02, respectively (Table 5).
CONCLUSION
In this study, a robust Modified biweight A chart is proposed and the performance of the proposed chart was compared to some existing control charts, namely the Standard Deviation chart (S) and the MAD chart (MAD). The change point estimator was employed to assess their performances by identifying the change point of a step change following a signal from the respective dispersion control charts. Overall, the results show that the estimator provides accurate and precise estimates of the time of a permanent step change in the process variance. With respect to expected length, the proposed robust Modified biweight A chart performs consistently and gives the best overall performances for the Type II, III and IV disturbances. Moreover, it is also comparable for Type I disturbance with the magnitudes of shift 2 or 3. Through extensive numerical simulations and performance evaluations, the proposed robust Modified biweight A control chart which incorporates outlier-resistant framework is more efficient than the existing charts. It is expected that the proposed approach which is relatively simple to implement, can play an importance part in quality improvement effort.
APPENDIX
The derivation of the Maximum Likelihood Estimator (MLE) of τ, that is
the process dispersion change point estimator, is presented in this section
(Samuel et al., 1998). Let denote the MLE of
the change point τ be
| (21) |
and the logarithm of the likelihood function is
| (22) |
In which case, for any given τ = t, the likelihood function is maximized by
| (23) |
By substituting
| (24) |
Hence, the maximum likelihood of τ is the value of t that maximizes the log-likelihood function defined as follows:
| (25) |
ACKNOWLEDGMENTS
The authors would like to express his sincere appreciation to the referees who suggest valuable and constructive advices to improve the quality of this article.