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 
and R charts. The
chart plots the subsample averages as a function of time and the R chart plots
the subsample ranges.
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,
and R charts are used to detect mean shifts, changes in variability and outliers.
These phenomena are known as disturbances of i.i.d. N (μ,σ) core process.
In order to set up and
R charts, one must either have priori knowledge of μ and σ or estimate
these parameters by examining the process over a period of time. Such estimates
are usually made on the basis of group of subsamples gathered during an initialization
period. With the similar conceptual framework, we focus our attention on the
question of detecting a change in variance where the S and R chart are known
to be two most widely used control charts for detecting changes in process variability
with S chart preferable to R chart for moderate-to-large sample sizes. Further
details about the S chart for non-Gaussian variables were discussed by Sim
(2000). The general considerations for the performance of shewhart control
chart can be referred to Jamali et al. (2006,
2007).
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:
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 
.
Following an unknown point in time τ (known as the process change point),
the process variance is assumed to have changed from 
to 
where δ is the unknown magnitude of the change. It is supposed that once
this step change in the process variance had occurred, the new level of 
will remain until the special cause has been observed and eliminated.
By letting Si be the standard deviation calculated for the ith subgroup, where
Xij is the jth observation in subgroup i and 
is the ith subgroup average. Analogously, Ri = max{Xij}-min{Xij}
is defined as the range of the ith subgroup. In the S chart framework we will
suppose that ST is the first subgroup standard deviation to beyond
the control limit and that this signal is not false alarm. However, if R chart
is used, then we will assume that RT is the first subgroup range
to exceed the control limit. Accordingly, S1, S2,...,
Sτ (or R1, R2,..., Rτ)
are the subgroup standard deviations (ranges) from in-control process while
Sτ+1, Sτ+2,..., ST (or Rτ+1,
Rτ+2,..., RT) are from the changed process.
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
where
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
where, 
is subsample median, k = 1, 2, 3,..., g and g is total number of subsamples.
When n is an odd number, subtracting out the median will produce one zero value
which is dropped, giving a centered subsample of size n–1 which is described
by Rocke (1983) as
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
where, kn,g,c is a factor chosen in such a way that 
We will call this biweight A estimator S'7 as B7 with tuning constant
c = 7 in the computation of the residuals. The value of c = 7 is chosen in order
to make an estimator sacrifices some efficiency in the absence of disturbance
and increase its efficiency when disturbances do occur. By means of Monte Carlo
studies, this choice of c performs the best and acceptably when viewed on the
whole.
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
where
and
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 
as modified biweight A estimator having the same appearance as in Eq.
3 with the altered form of uk, j given by Eq. 9.
To simplify, B*7 will represent 
to mean the Modified biweight A estimator with tuning constant c = 7 applied
to the residuals. Hence, the process standard deviation σ is evaluated
by
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
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
and
In general, the average of the g subsamples scale estimates is defined as
Consequently, the process standard deviation σ is then estimated by
| Table 1: |
Simulated constant values c4 for different scale
estimators |
 |
|
Let us now denote the proposed 
as robust Modified biweight A (D7) estimator thereafter with tuning constant
c = 7 is applied in the computation of the residuals. Again dn, c
is a factor chosen such a way that 
Tabled values of dn, c which is in analogy to A = c4 for
n = 5 and n = 10 are shown in Table 1 determined from 100,000
simulations for each setting (Montgomery, 2005).
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
Generally, the control limits for a scale chart is defined as
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 
estimate. For instance, the
for the scale estimator of MAD chart is
Similarly, the
for the scale estimator of the proposed robust Modified biweight A, denoted
as D7 previously henceforth is given by
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
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
.
Subsequent to an unknown point in time τ (known as the process change point),
we presume that the process variance changes from
to 
where δ is the unknown magnitude of the change. Supposedly, once this step
change in the process variance occurs, it will remain at the new level of 
until the special cause has been observed and eliminated. In our study we will
focus mainly on δ = 2.0 in virtue of the significance of the results and
ease of comparison.
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 
and expected length from the 10,000 simulation runs will be computed |
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
and expected lengths based on the three different charts under investigation
are given in Table 2 to 5.
Type I: Symmetric variance disturbance: The data under symmetric variance
disturbances were considered in which average change point estimates
are tabulated for various sizes of magnitude a ranging from 2.0 to 7.0 in Table
2. As the actual change point for the simulations is at time 100, the average
estimated time of the process change,
should be close to 100. For the process step change of standardized magnitude
δ = 2.0, the averages estimated time of the process change for three charts
considered are in most cases distributed around 100 which is fairly close to
the actual change point of 100. However, in terms of expected length, the proposed
robust Modified biweight A chart seems to perform consistently better compared
to the other two charts.
| 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 
thus, for observations x11, x12,..., x1n,...,
xT1,..., xTn, the MLE of τ is the value of τ
that maximizes the likelihood function or, equivalently, its logarithm. Then,
the likelihood function is given by
and the logarithm of the likelihood function is
In which case, for any given τ = t, the likelihood function is maximized by
By substituting 
in the log-likelihood function, the result is
Hence, the maximum likelihood of τ is the value of t that maximizes the log-likelihood function defined as follows:
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.
