A Multivariate CUSUM Procedure and a Most Active Bivariate Chart

Quality Technology & Quantitative Management Vol. 3, No. 4, pp. 401-414, 2006

Q T Q M ?

ICAQM 2006

A Multivariate CUSUM Procedure and a Most Active Bivariate Chart

Ruojia Li 1 and Richard A. Johnson 2

1

Eli Lilly and Company, Lilly Corporate Center, Indianapolis, IN, USA 2Department of Statistics, University of Wisconsin-Madison, Madison, WI, USA

(Received March 2005, accepted October 2005) ______________________________________________________________________ Abstract: We propose a new multivariate cumulative sum (CUSUM) scheme whose components are the difference between the Page [10] univariate CUSUM statistic for detecting an increase and that for detecting a decrease. Our procedure also applies when monitoring principal components. We then suggest a new quality chart procedure where, at each new observation, the results for the two most active components are graphed. These two components are selected to have the largest value of the appropriate bivariate CUSUM statistic. This novel graphical approach, intended to augment the information in the chart based on all of the variables, is illustrated with an application to data from an automobile assembly process. We compare our new scheme with the classic charts and a number of other existing multivariate CUSUM schemes, including the scheme proposed in Crosier[2]. In our simulations, among existing procedures, Crosier[2] has the best overall average run length (ARL) performance. Our new multivariate scheme has ARL performance that is usually comparable with that of Crosier’s scheme. 2T Keywords: CUSUM statistic, graphical procedure, most active bivariate plot process capability plots, multivariate quality monitoring.

______________________________________________________________________

1. Introduction

n this paper, we address the general problem of detecting a change in the level of a sequence of independent random vectors. Nowadays, with automated data collection a common practice, several characteristics are continually monitored in order to detect a change in the quality of a manufactured product.

I Consider a sequence of random vectors where the components of each random vector correspond to the measured characteristics being monitored. The sequence is “in control” if all of the observed variation is due to chance and the process mean is at some acceptable value. The sequence is considered to be “out of control” if, from some time point on, the mean has shifted to some other value. We are interested in sequential quality monitoring schemes for detecting a small or moderate shift in the mean of the process.

Sequential control schemes are evaluated in terms of their average run length (ARL). An effective scheme should have a large ARL when the process mean is at its desired level and have a small ARL when the process mean deviates from this value.

We propose a new multivariate CUSUM scheme in Section 2. In Section 3, we present the proposed most active bivariate chart. At each new observation, a figure is created for the two components having the largest value of the corresponding bivariate CUSUM statistic. Our CUSUM scheme for all of the variables and the most active bivariate chart are both illustrated with an application.

402 Li and Johnson

In Section 4, we take a simulation approach to compare the exact properties of our procedure and several competing monitoring schemes. The simulation studies show that, among the existing multivariate procedures, the CUSUM scheme proposed by Crosier[2] has the best overall performance. The ARL performance of our new scheme is usually comparable with Crosier’s scheme and better in a few circumstances.

2. Our New CUSUM Scheme

All of the multivariate procedures we treat are based on the sequence of independent 1p × random vectors 1X , ..., n X , ... where each has the target mean value a and the same covariance matrix . If is unknown, it will be estimated from a long series of observations taken when the process is stable. As is common with most of the literature, we describe our multivariate monitoring scheme in terms of a known covariance matrix .

ΣΣΣInitially, we define our new procedure in a univariate setting where , , ... is a sequence of independent random variables and a is the target mean value. We first define two cumulative sums with drifts:

1X 2X (1)1

()n

n i i S X k ?==?∑,

(2)1()n n i i S X k +==?∑,

As in Page’s CUSUM scheme, the two constants satisfy and . When the sequence of random variables , , ... are normally distributed, typically we can set k +>a a k ?<1X 2X k a k σ+=+ and k a k σ?=?, where σ is the standard deviation of and is one-half of the specified mean shift (expressed in standard deviations) that should be detected by the scheme.

X k Following Page [10], we then define

(1)(1)max ,n m n m S S ?≤=n S ?n ? (1)

(2)(2)max .n n m n m S S S +≤=?

(2)

With an eye towards multivariate extensions, we consider the difference . In particular, our new univariate CUSUM scheme signals when n S S +?n n S S h σ+??≥ (3)

for some specified positive constant .

h In the multivariate setting, each p -dimensional random vector , has covariance matrix =i X ,1,2,[,,...,]'i i i p X X X Σ and target mean vector 12=a [,,...,]'p a a a . To generalize our univariate CUSUM procedure to the multivariate setting, we begin by defining the multivariate CUSUM statistics and .

(1)n S (2)n S (1)(1)(1),1,[,...,]'n n n p S S S =1

(),?==?∑n i i X k (4) (2)(2)(2),1,[,...,]'n n n p 1

(S S S =);+==?∑n i i X k (5)

A Multivariate CUSUM Procedure 403

where ?=k 1[,...,]',p k k ??+=k 1[,...,]'p k k ++ and j j k a +> j j k a ?< for

j =1, 2, …, p . When each of the random vectors has multivariate normal distribution, typically, we take

12,,...X X ?=k ?a k σ and ?=k +a k σ (6)

so that (4) and (5) can be written as (1)=n S 1=( ?+∑n i )k σ and 1=i (2)=n S (??∑n

i X a X a i )k σ, Here =σ1[,...,]'p σσ, where j σ is the standard deviation of the j -th component of the i X ’s, and the constant k can be chosen as one-half of the specified shift in mean on the j -th component (expressed in standard deviations).

For each component j = 1, ..., p , we define the two sequences of statistics

(1)(1),,max n j m n m j n j S S ?≤=?,S ,,

(2)(2),,max n j n j m n m j S S S +≤=?.

and then the two multivariate statistics n S + and n S ? where

+=n S ,1,[,...,]'n n p S S ++ (7)

?=n S ,1,[,...,]'n n p S S ?? (8)

We now define our CUSUM statistic in terms of the differences . ,D n S n n S S +??,D n S =()()112n n n n S S S S ['+??+??Σ?]

(9) Our multivariate CUSUM scheme signals when ,D n S =()()112[']n n n n S S S S +??+??Σ?≥ h (10)

for some specified control limit h . For any given Σ and , simulation can be used to determine the control limit h to attain a specified ARL when the process is on target.

k When the dimension p is large, it is common practice to monitor principal components. Johnson and Wichern[6] describe an widely applied scheme due to Kourti and MacGregor

[9] where the first two principal components are monitored visually and the remaining principal components are combined into a single statistic for a second chart. Alternatively, we can apply our statistic to all or to just a few of the principal components. The principal components are based on the eigenvalue and eigenvector pairs (k λ, ) determined from

k e Σλ=j e j e (11)

and ordered so that 12...p λλ≥≥≥λ. The value of the first principal component, for random vector , is

i X ,j i Y =()1'?i e X μ

404 Li and Johnson and those of the other principal components, for random vector i X , are

,j i Y =()'j i e X μ? for j =1, …, p

As long as the smallest eigenvalue is not too small, it is convenient to standardize the

values of the principal components so that they all have variance 1. That is, we use

.,j i j i Z ==()'μ?j i e X for j

=1, …, p (12) The value for the statistic calculated from the ,D n S ,j i Z is the same as its value calculated from the ,j i Y but the value for only needs to be determined once for a given dimension since the resulting covariance matrix is always the identity.

h 3. A Most Active Bivariate Chart

We will introduce most active bivariate charts in the context of an application where we also give an application of our multivariate CUSUM procedure based on the statistic . We refer to the automobile sheet metal assembly data used in Johnson and Li [5], Section 4. More detail on the body assembly process is given in Ceglarek and Shi[1].

,D n S The data, collected by Dr. Darek Ceglarek, consist of 50 observations on 6 variables of which four were measured when the car body was complete and two were measured on the underbody at an earlier stage of assembly.

All measurements were taken by sensors that recorded the deviation from the nominal value in millimeters.

1x =deviation at mid right hand side – body complete

2x =deviation at mid left hand side – body complete

3x =deviation at back right hand side – body complete

4x =deviation at back left hand side – body complete

5x =deviation at mid right hand side of underbody

6x =deviation at mid left hand side of underbody

From the individual variable plots and chart, given in the appendix, the first 30 observations appear to be stable. Although this sample size is rather small, for illustrative purposes, we obtain the sample covariance matrix and the mean from these 30 observations.

2T S = 0.0626 0.0616 0.0474 0.0083 0.0197 0.00310.0616 0.0924 0.0268 0.0008 0.0228 0.01550.0474 0.0268 0.1446 0.0078 0.0211 0.004??90.0083 0.0008 0.0078 0.1086 0.0221 0.00660.0197 0.0228 0.0211 0.0221 0.3428 0.01460.0031 0.0155 0.0049 0.0066 0.0146 0.03??66????????????????????

[0.5063 0.2070 0.0620 0.0317 0.6980 0.0650]'=??????x

A Multivariate CUSUM Procedure 405

Using the estimated covariance matrix S in place of Σ, we chose = 0.5 and then used simulation to obtain the value = 7.54 so the scheme has an in-control ARL of about 200. In this simulation, sequences of independent normal random variables having mean vector and covariance matrix k h =0a Σ were generated. For this example, the value for h was estimated using 1000 sequences.

The values for our statistic, are shown in Figure 1. A shift is indicated at the 42nd observation and from 44th observation on.

Figure 1. Chart of based on original data. ,D n S Since six dimensions cannot be simultaneously graphed, we suggest that at each new observation, the results for the two components having the largest value of the bivariate CUSUM statistic be graphed. We call this the most active bivariate dimension and the most active bivariate chart shows the ellipse, determined from squaring both sides of (10), for all of the values of the appropriate 2 components of the differences +??m S S m , for m = 1, 2, ..., n . The bivariate value for the last point, the current position +??n n S S +??S S , is shown as a solid dot to distinguish it from the earlier points. When component j is active, the axes is simply marked Sj but the plotted values employ the j ?th components of the differences . m

m The most active bivariate chart, for the automotive data, is given in Figure 2. The first chart is for observation number 31, the first after those used to estimate the covariance matrix. It presents components 2 and 3. The second chart, for observation number 32, presents component 1 and 3 and so on. We immediately see which component or two are contributing heavily to the multivariate CUSUM statistic.

If the process is stable for long periods of time, for visual clarity, it may be best to retain just the most recent 10-30 observations. It may even be desirable to shade the observations from black for the most recent to white for the oldest to reflect time order.

We can also apply our procedure to the principal components. First choosing k = 0.5, we used simulation to produce the value h = 6.94 for an in control ARL of about 200.

406 Li and Johnson

Figure 2. Most active bivariate chart for using the original data. ,D n S

The eigenvectors and eigenvalues of the sample covariance matrix S , which define the sample principal components and their variances, are given in Table 1.

It is clear that the first three principal components explain a large proportion of the variance but we will use all six standardized principal components to determine our statistic The chart is shown in Figure 3. This chart signals at the 38th and 39th observations as well as from the 46th on.

,.D n S

A Multivariate CUSUM Procedure 407

Table 1. Eigenvalues and eigenvectors of S . Variable

1e 2e 3e 4e 5e 6e

1x 0.119 0.469 0.075 0.291 -0.267 -0.777

2x

0.13 0.458 0.251 0.624 -0.037 0.566 3x

0.143 0.717 -0.116 -0.632 0.175 0.147 4x

0.096 0.053 -0.956 0.257 -0.062 0.072 5x

0.968 -0.231 0.070 -0.064 -0.032 0.003 6x

0.052 0.013 -0.008 0.238 0.944 -0.22 i λ

0.3544 0.1864 0.1076 0.0972 0.0333 0.0088

Figure 3. Chart of based on the standardized principal components.

,D n S The most active bivariate chart for the principal components, shown in Figure 4, begins with the 31st observation and that case is based on the fourth and fifth principal components. The most active bivariate chart in Figure 4 uses the standardized values of the principal components (12) so that the ellipses reduce to circles. This most active chart signals a shift at the 39th observation and the fourth and fifth principal components are the two involved.

408 Li and Johnson

Figure 4. Most active bivariate chart based on principal components. The most active bivariate chart is not tied exclusively to our CUSUM procedure. For instance, the classical chart could be used in place of the chart in Figure 1 or Figure 3. We can then construct the most active bivariate chart using the usual ellipse format chart (see Johnson and Wichern [6]) We create the ellipse based on the constant for the p = 6 variable chart. This chart shows directly which one or two components had shifted in level.

2T Figure 5 illustrates the most active bivariate chart where the bivariate data themselves are plotted on the two axes. The first ellipse, for observation 31, shows the components and . Observation number 39, see the ninth ellipse displayed, is out of control mostly because of the sixth variable.

2T 1x 3x

A Multivariate CUSUM Procedure 409

Figure 5. Most active bivariate chart for using the original data.

2T 4. Comparisons of the ARL Performance of Multivariate Schemes

We begin by reviewing four popular multivariate procedures (see also Fuchs and Kenett[3], Johnson and Li[5], Lowry and Montgomery [7], and Yang and Trewn [13]). All of the competing multivariate procedures are also based on the sequence of independent 1p × random vectors 1X , ..., n X , ... where each has the target mean value and the same covariance matrix . Most of the multivariate schemes discussed below involve a constant k and a positive control limit .

a Σh The traditional multivariate chart [4] reduces each multivariate observations to a scalar by defining

2T

410 Li and Johnson

()()21'?=?n n n T X a X a Σ?k ? (13)

The multivariate scheme signals a shift in mean when first exceed a specified level . That is, the scheme signals when .

2T 2n T h 2n T h ≥The CUSUM of (COT) scheme forms a CUSUM of the scalar statistics defined in (13).

n T n T Let and be specified constants. Then, the COT statistic is iteratively defined as

00S ≥0k >1max(0,),n n n S S T ?=+ for 1,2,...n = (14)

where k is a specified positive constant.

The COT scheme signals when for a specified level .

n S h ≥h The Crosier [2] multivariate CUSUM scheme starts at 00=S . For a specified constant , the statistic is iteratively defined as k n 10 if ()(1) otherwise

n n n n C k S S X a k C ?≤?=?+??? (15) where 1111[()'()].???=+?Σ+?n n n n n C S X a S X a

For another specified constant , the Crosier multivariate scheme signals a shift in mean when

h =n Y '112≥h []n n S S ?Σ

Lowry et al . [8] proposed a multivariate exponentially weighted moving average (EWMA) that begins from . The multivariate EWMA statistic is iteratively defined as

00=Z n Z 1()()?=?+?n n Z R X a I R Z n for n =1, 2, (16)

where the weight matrix 1(,...,)p diag r r =R , 01j r ≤≤, j = 1, ..., p . This reduces to the situation where a univariate EWMA is applied to each individual component. When there is no priori reason to weight the p quality characteristics differently, it is suggested to use a common value , where 1...p r r ===r 01r ≤≤ is a constant.

The multivariate EWMA scheme signals, for a specified constant h , when

11[]n n 'h ?≥Z Z Σ.

The choice of weight matrix R has a considerable influence on the resulting ARL behavior.

To compare different schemes, we set the in-control ARL’s to be nearly equal and compare the ARL’s when there is a shift. Several authors have also compared various existing multivariate monitoring procedures including Pignatiello[11] and Lowry[7].

In our simulation, we chose the reference value and decision value so that the resulting schemes, which are designed to detect a shift of one standard deviation from the k h

A Multivariate CUSUM Procedure 411

target, have an in-control ARL of about 200.

We first generate bivariate normal observations with either =ΣΙ (uncorrelated)

or =Σ 1.0.6.6 1.0????????

(correlated) to compare the performances of different schemes. Table 2 and Table 3 show the results of our simulation for our CUSUM procedure based on , the multivariate chart, the CUSUM of T (COT) scheme, Crosier’s multivariate scheme, and the multivariate EWMA schemes with weight parameter ,D n S 2T r =0.1, 0.4, and 0.8.

Table 2. ARL Comparison under Bivariate Normal Observations(uncorrelated).

The estimated standard error is given below each estimated ARL. Shift h 2T 10.63 COT 4.04 Cros 5.55 .D n S 5.1

EW(.1) 8.73 EW(.4) 10.37 EW(.8) 10.6 (0, 0) 198.8 197.3 204.5 200.8 200.6 203.9 202.9

2.8 2.7 2.9 2.7 2.8 2.8 2.8 (1, 0) 42.0 20.4 9.5 9.4 9.8 1

3.1 28.1

0.6 0.2 0.1 0.1 0.1 0.1 0.4 (0, -1) 41.0 20.3 9.6 9.5 9.9 13.5 27.8

0.6 0.2 0.1 0.1 0.1 0.2 0.4 (.7, .7) 42.0 20.0 9.4 10.3 9.7 13.1 28.5

0.6 0.2 0.1 0.1 0.1 0.2 0.4 (.7, -.7) 41.8 19.9 9.3 10.1 9.6 12.9 28.9

0.6 0.2 0.1 0.1 0.1 0.2 0.4 (.5, 0) 117.2 81.0 28.8 32.2 27.5 54.0 94.9

1.6 1.1 0.3 0.4 0.3 0.7 1.3 (.4, -.4) 116.3 80.2 28.6 35.7 27.1 5

2.6 96.4

1.6 1.1 0.3 0.5 0.3 0.7 1.3 (2, 0) 7.0 4.7 4.0 3.6 4.2 3.4 4.9 0.1 0.0 0.0 0.0 0.0 0.0 0.1 (1.4, -1.4) 7.1 4.7 4.0 3.8 4.2 3.4 4.7

0.1 0.0 0.0 0.0 0.0 0.0 0.1 For each choice of and shift in mean, a series of observations was generated and the multivariate statistics were calculated until a shift is signaled. This procedure was repeated 5000 times and we calculated the ARL and its estimated standard error for each scheme. The values of estimated ARL and standard error are presented in Table 2 and Table 3. The estimated standard errors are given below the corresponding ARL.

ΣFrom our simulation and existing literature, we conclude that, due to the fact that the value of the statistic only depends on the most current observation, it is not sensitive to small and moderate shifts in the mean of a process even if the shift is persistent. By taking the CUSUM of , the COT procedure has ARL performance that is significantly improved over that of the chart. The ARL of Crosier’s scheme is considerably better than that of the COT scheme. The performance of multivariate EWMA schemes, depend heavily on the value of the weight parameters. If the weight is appropriately selected, the 2T T 2T

412 Li and Johnson multivariate EWMA will have very good ARL performance, which is comparable with Crosier’s scheme. Our proposed new procedure provides good overall ARL performance that is, in most cases, comparable with that of Crosier’s scheme and could be better in a few situations. Additional simulation results on ARL, for the existing tests, are given in Johnson and Li [5].

Table 3. ARL Comparison with Bivariate Normal data(correlated).

The estimated standard error is given below each estimated ARL

Shift h

2

T

10.63

COT

4.04

Cros

5.6

.D n

S

5.4

EW(.1)

8.71

EW(.4)

10.36

EW(.8)

10.6

(0,

0) 198.8 197.3 212.8 206.4 198.1 202.9 202.9

2.8 2.7

3.0 2.9 2.8 2.8 2.8

(.8,

0) 41.9 20.4 9.6 11.6 9.7 13.2 28.7

0.6 0.2 0.1 0.1 0.1 0.2 0.4

(0,

-.8) 41.0 20.4 9.6 11.7 9.9 13.4 27.7

0.6 0.2 0.1 0.1 0.1 0.2 0.4

(.4,

.4) 42 20.0 9.5 15.1 9.7 13.1 28.5

0.6 0.2 0.1 0.2 0.1 0.2 0.4

(.4,

-.4) 118.0 82.3 29.0 30.5 27.0 54.3 95.2

1.6 1.1 0.3 0.4 0.3 0.7 1.3

(.4,

0) 116.8 82.2 29.1 42.2 27.3 52.9 97.0

1.6 1.1 0.3 0.5 0.3 0.7 1.4

(.2, -.2) 173.6 154.7 85.5 89.9 74.0 126.8 157.1

2.4 2.1 1.1 1.2 1.0 1.7 2.2

(1.6,

0) 6.9 4.7 4.0 4.1 4.2 3.4 4.8

0.1 0.0 0.0 0.0 0.0 0.0 0.1

(.9, -.9) 41.6 20.1 9.6 9.5 9.8 13.2 28.3

0.6 0.2 0.1 0.1 0.1 0.2 0.4 References

1. Ceglarek, D. and Shi, J. (1995). “Dimensional Variation Reduction for Automotive

Body Assembly”, Manufacturing Review, Vol. 8, No. 2, 139-154.

2. C rosier, Ronald B. (1988). “Multivariate Generalizations of Cumulative Sum Quality

Control Schemes”, Technometrics, Vol.30, No. 3, 291-303.

3. Fuchs, C. and Kenett, R. S. (1998). Multivariate Quality Control: Theory and Applications.

New York: Marcel Dekker.

4. Jackson, J. E. (1959). “Quality Control Methods for Several Related Variables”,

Technometrics, Vol. 1, 359-377.

5. Johnson, Richard A. and Li, Roujia (2005). Multivariate Statistical Process Control

Schemes for Controlling a Mean. Editors H. Pham and C. F. Wu, Handbook of Engineering Statistics, Springer.

6. Johnson, Richard A. and Wichern, Dean W. (2002). Applied Multivariate Statistical

Analysis. New Jersey: Prentice Hall.

7. Lowry, C.A. and Montgomery, D.C. (1995). “A review of multivariate control charts”,

IIE Transactions, Vol.27, 800-810.

A Multivariate CUSUM Procedure 413

8. Lowry, C. A., Woodall, W . H., Champ, C. W . and Rigdon, S. E. (1992). “A multivariate exponentially weighted moving average control chart”, Technometrics , Vol. 34, No. 1, 46-53.

9. Kourti, Theodora and MacGregor John F . (1996). “Multivariate SPC Methods for Process and Product Monitoring”, Journal of Quality Technology , Vol. 28, No. 4, 409-428.

10. Page, E.S. (1954). “Continuous Inspection schemes”, Biometrics , Vol. 41.

11. Pignatiello, J. J. and Runger, G. C. (1990). “Comparisons of Multivariate CUSUM Charts”, Journal of Quality Technology , 22(3), 173-186.

12. Tracy, N. D., Young, J. C., and Mason, R. L. (1992) “Multivariate Quality Control Charts for Individual Observations”, J. Quality Technol . Vol. 24, 88-95.

13. Yang, K. and Trewn, J. (2004). Multivariate Statistical Methods in Quality Management. New York: McGraw-Hill.

Appendix: Stability of the Car body Data

As in Johnson and Li [5], we conclude that the first 30 observations are stable because the individual variable time plots and the graph of the statistic do not exhibit noticeable lack of stability for the first 30 observations. These plots are 2T

Figure 6 . Individual variable plots of the car body observations

414 Li and Johnson

Figure 7. chart based on the first 30 car body observations.

2T

Authors’ Biographies:

Richard A. Johnson is a Professor of Statistics at the University of Wisconsin. His research and consulting interests include reliability and life length analysis, applied multivariate analysis, and applications to engineering. He is a Fellow of the American Statistical Association, a Fellow of the Institute of Mathematical Statistics, and an elected member of the International Statistical Institute. He has been editor of Statistics and Probability Letters, since in began 25 years ago, and co-author of six books and several book chapters, and over one-hundred papers in the statistical and engineering literature. Ruojia Li is a research scientist at Eli Lilly & Company. She recently completed her Ph. D. degree in Statistics from the Department of Statistics at the University of Wisconsin. This research was conducted while she was a student, under the supervision of Professor Johnson. Her research interests include multivariate quality monitoring schemes and statistical applications in phamecutical research.