Comparative statistical analysis of selected control charts for highly capable processes

68 2019, XXII, 2 Business Administration and Management DOI: 10.15240/tul/001/2019-2-005 Introduction The statistical process control (SPC) is widely used in industry. Its aim is to achieve process stability and improve capability through the reduction of variability and it is included for example in the Six Sigma or Lean Six Sigma methodologies. The control of attribute data represents a considerable part of it. Until recently, the same approach was used to monitor variables or attributes data. In this approach, subgroups of items are taken from a process and sample characteristics are plotted in a control chart to see whether their variation is only random or whether it is affected by an assignable cause. The presence of such cause is indicated by exceeding the control limits that are based on the sample characteristic distribution. When items are classifi ed as conforming or nonconforming, the proportion or the count of nonconforming units in a process (fraction nonconforming) is traditionally monitored by a p chart or np chart. Due to new manufacturing technologies and concepts, many processes are of such a high quality that the fraction nonconforming or the probability of observing a nonconforming unit is very small. This probability is a parameter of the background binomial distribution and the size of the subgroups would have to be enormously large to enable the normal approximation, on which the computation of control limits is based. The impossibility to accomplish such conditions has the following consequences: real properties of a chart such as the average run length (ARL) or the risk of false signal  differ from those assuming a normal distribution and the lower control limit of a p chart is located at zero and does not enable the recognition of the process improvement. Moreover, the appearance of a p chart containing most zero points is inappropriate. Consequently, the concept based on the sample characteristics is no longer useful for high-quality processes. Several alternatives based on the Bernoulli process rather than on the normal approximation have been recently presented. These methods assume continuous inspection. But not necessarily 100% items must be inspected. For example Reynolds and Stoumbos (1999) admit a situation when “the production rate is higher than the inspection rate”, Szarka and Woodall (2011) also mention an interval sampling when the items are inspected at scheduled periods. Bourke (2001) emphasizes that “the pattern of the sampling inspection can be quite haphazard without causing any diffi culty” for the performance of the control chart he had studied. In the Bernoulli process, random values xi = 0 or xi = 1, i = 1, 2, …, express whether an inspected item is conforming or nonconforming, respectively. The in-control state of a process is defi ned by a constant probability p0 of an occurrence of a nonconforming item. Usually a sustained shift to the out-of-control state is considered; the process repetitively produces units at level p0 until it suddenly shifts to an unacceptable level p and remains at this level until a remedial measure is taken. Especially the situation p > p0 must be detected as soon as possible and sometimes a minimum unacceptable value for p is given. Sometimes the inverse inequality p < p0 may be of interest, too, because it indicates a process improvement. In the new approach, the number of units Yi being inspected until the r-th nonconforming one has been observed is followed, where r ≥ 1. Either the geometric (for r = 1) or the more general negative binomial distribution of Yi are assumed. Apart from the Shewhart-type charts based on these two distributions (often called CCC or CCC-r charts) the CUSUM or EWMA charts can be constructed. COMPARATIVE STATISTICAL ANALYSIS OF SELECTED CONTROL CHARTS FOR HIGHLY CAPABLE PROCESSES Eva Jarošová, Darja Noskievičová EM_2_2019.indd 68EM_2_2019.indd 68 19.6.2019 15:11:0419.6.2019 15:11:04 69 2, XXII, 2019 Business Administration and Management The control limits of the CCC and CCC-r charts are the probability limits; they are constructed directly as percentiles for a given risk of false signal . Cumulative sum (CUSUM) charts use information from all the prior observations (in this case the run lengths Yi defi ned above) and they are considered to be an effi cient alternative to the Shewhart chart when small process shifts are of interest. Usually r = 1 is considered and the corresponding CUSUM chart is called a geometric (or CCC) CUSUM. Exponentially weighted moving average (EWMA) charts are also based on all the prior observations and have a similar effi ciency as the CUSUM charts. Numerous authors studied effi ciency of these control charts (see further) but only few case studies have been published, see (Chang & Gan, 2001) or (Di Bucchianico et al., 2005) In this paper the CCC-r charts and the geometric CUSUM chart are analyzed. The aim of the paper is to compare the CCC-r and geometric CUSUM charts by their speed of detecting an upward shift in the fraction nonconforming when they are applied to a process with the target fraction nonconforming of the order of 100 ppm. Both types of charts are then applied to a real process of the electronic assembly that has become a basis for the study. To our best knowledge no study dealing with such a small fraction nonconforming together with the low risk of false alarm has been published yet. 1. Control Charts Description For our comparative analysis we have chosen three types of the control charts based on the Bernoulli process: the CCC, CCC-r and CCC-CUSUM charts. Their theoretical basis will be described in this section. 1.1 CCC Chart The CCC chart was designed by Calvin (1983) and further studied and expanded by Goh (1987), Bourke (1991) or Xie et al. (1995). The name CCC stands for the cumulative count of conforming units, but more frequently, the variable Y that is monitored includes also the immediately following nonconforming unit. Other names such as the conforming run length CRL in (Wu, 2000) or run-length RL1 in (Bourke, 1991) can be met in literature. We consider the latter case but we will use the inaccurate but more known name CCC in our paper. Variable Y follows the geometric distribution G(p) with the probability function 1 () (1 )   y fy p p , y = 1, 2, ..., (1) where p is the probability of observing a nonconforming unit in any inspection. The centreline of the CCC chart corresponds to the median of the geometric distribution , (2) and the two-sided probability control limits (the upper control limit UCL and the lower control limit LCL) are based on the chosen risk of false alarm α (Xie et al., 2002): ln 2 ln(1 )   UCL p , ln 1 2 ln(1 )       LCL p . (3) As these limits are highly asymmetric, a logscale is sometimes used (see for instance Xie et al., 2002). As soon as a nonconforming unit is observed, the value Yi is plotted on the chart and counting starts again from zero. The interpretation of this control chart is quite different from the interpretation of the conventional Shewhart p chart: the point above UCL indicates the probable process improvement while the point below LCL shows the probable process deterioration. When only an upward shift in the process fraction nonconforming is of interest, the onesided lower control limit is given by  ln 1 ln(1 )    LCL p . (4) For the conventional Shewhart or CUSUM charts, when both the sample sizes and the sampling intervals are equal, the average run length ARL can be used to assess the effi ciency of a chart. It represents the expected number of samples before the fi rst signal occurs. This measure is inappropriate with charts based on the Bernoulli process. The observations Yi EM_2_2019.indd 69EM_2_2019.indd 69 19.6.2019 15:11:0419.6.2019 15:11:04 70 2019, XXII, 2 Business Administration and Management correspond to different numbers of inspected units and this fact must be taken into consideration. That is why another measure has been introduced. Reynolds and Stoumbos (1999) use the average number of observations to signal ANOS that represents the expected number of inspected items to the signal. The average time to signal ATS, the average number inspected ANI, the average number of items sampled ANIS are different names for the same measure used in Wu et al. (2000), Bourke (1991) and Chang and Gan (2001), respectively. In our paper we use the acronym ANOS, or ANOS0 for the in-control process with fraction nonconforming p0. 1.2 CCC-r Chart To improve the sensitivity of detecting small upward shifts in the fraction nonconforming, the chart based on r successive run lengths has been developed. Bourke (1991) fi rst considered RL2 chart for r = 2 with the moving sums RL2 = Yi-1 + Yi for i = 2, 3, ..., but later the separated sums Y1 + Y2 , Y3 + Y4 etc. were used, see e.g. (Wu et al., 2001). More general case 2r is considered in Xie et al. (1999), Kaminsky et al. (1992), Xie and Goh (1997). Often r values of 2 or 3 are recommended, see Wu et al. (2000) or Schwertman (2005). Ohta et al. (2001) constructed an economic model to fi nd the best value of r. All these charts are based on a negative binomial distribution of the monitored variable. Common name for these charts is the CCC-r chart, although other names may appear, e.g. SCRL chart in Wu et al. (2000). The cumulative count of units until the r-th observed nonconforming unit follows a negative binomial distribution with the probability function 1 () (1 ) 1        ryr y fy p p r, y = r, r + 1, r + 2,... (5) Assuming that the probability of a nonconforming unit is p, the two-sided control limits and the centerline are determined by solving the following equations (see Chen & Cheng, 2010): 1 () (1)1 12           UCL ryr yr y FUCL p p r, (6) 1 () (1) 12           LCL ryr yr y FLCL p p r, (7) . (8) When only an upward shift is followed, the one-sided lower control limit must satisfy the equation 1 () (1) 1           LCL ryr yr y FLCL p p r (9) for the specifi ed risk of false signal α. The interpretation of the CCC-r chart is similar to the CCC chart interpretation. With an increasing value of parameter r the CCC-r chart is more sensitive to small upward shifts in p. But on the other hand more and more Bernoulli observations are needed to obtain a point in the chart and therefore the inspection cost increase (see e.g. Ohta & Rahim, 2001). For that reason it is necessary to set the optimal parameter r. 1.3 Geometric CUSUM Chart The geometric CUSUM chart was proposed by Bourke (1991). It is based on the schemes introduced in (Page, 1954) that can be expressed by 1 max(0, )    iii SSYK , (10) 1 min(0, )    iii SSYK , i = 1, 2, ... These schemes can be used to detect downward or upward shifts in a distribution of Yi. If variables Yi follow a geometric distribution and the process deteriorates, i.e. if the fraction nonconforming p increases, values of Yi less than Kpredominate and  i S becomes more and more negative. As soon as  i S drops under the specifi ed limit Hfor some i, an out-of-control signal is given. The value of Hdetermines the chart performance, i.e. ANOS0 and ANOS. The more negative it is, the longer the time to signal. To avoid the negative values of  i S, the lower Page scheme is often presented in the form 1 max(0, )    iii SSYK . Here the original expression given by (10) is used, with an omitted minus sign in the upper index. EM_2_2019.indd 70EM_2_2019.indd 70 19.6.2019 15:11:0419.6.2019 15:11:04 71 2, XXII, 2019 Business Administration and Management Values of the CUSUM statistic are then zero or negative and the lower limit H is negative. When a shift from in-control p0 to a higher value p is to be detected, the minimum value p1 that is considered to be inappropriate must be chosen. The level of the fraction nonconforming is acceptable for p ≤ p0 , and unacceptable for p ≥ p1. As a matter of fact, hypotheses H0: p = p0 versus H1: p = p1 are tested. Using the Wald sequential probability test and the geometric distribution given by (1), the reference value K is 10 01 0 1 (1 ) ln (1 ) 1 ln 1     pp pp Kp p . (11) Frequently S0 = 0 is used but a head-start can be chosen to detect an initial out-of-control state quickly (see Bourke, 1991). The determination of the lower limit H is somewhat complicated. When the Markov chain approach is used, the number of H transient states is quite large and computational problems arise. Bourke (1991) suggests approximating the geometric distribution by an exponential distribution, grouping states into a substantially less number of categories and then using an approach for a continuous variable based on transitions among the categories by Brook and Evans (1972). Reynolds and Stoumbos (1999) introduced the corrected diffusion approximation and presented formulas for the Bernoulli CUSUM chart, which is closely related to the geometric CUSUM chart: ** 22 0 20 1 exp( ) 1 ()   BB Hr Hr ANOS p rp r , (12) where 1 1 0 1 ln 1       p rp 10 2 01 (1 ) ln (1 )      pp rpp , * 00 0 () (1 )    BB HH pp p , * B H can be obtained from (12) by means of the Excel function Goal Seek. Once the parameter HB of the Bernoulli CUSUM chart has been found, the parameter HG of the geometric CUSUM chart can be expressed as (Szarka & Woodall, 2012): HG = – m(HB–1). (13) In this paper, the lower limit H was set at fi rst approximately using the approach by Reynolds and Stoumbos (1999), and then adjusted by means of simulations: 1. For a specifi ed value of ANOS0, an approximate value of the limit HB of the Bernoulli CUSUM that is implicitly included in the formula (12) was calculated and the corresponding value of the limit HG of the geometric CUSUM was determined from (13). 2. The value of HG was refi ned iteratively through simulations so that the resulting ANOS0 is as close to the specifi ed value as possible. It means that run lengths Yi for i = 1, 2, ..., from G(p0) were generated and the CUSUM characteristic Si was calculated as well as the cumulative sum of run lengths (or the number of Bernoulli observations) NOSi. As soon as Si decreased under the lower limit H for some i = I, the simulation was interrupted and NOSI was stored. This cycle was repeated 10,000 times, i.e. ANOS0 was estimated as an average from 10,000 NOS values. Then the value of HG was modifi ed and another 10,000 cycles were performed. This procedure was repeated until the nearest possible value to the theoretical ANOS0 was attained. The fi nal value H of the parameter obtained  34 7 00 0 0 0 0 00 0 0 0 0.410 0.0842ln 0.0391(ln ) 0.00376(ln ) 0.000008(ln ) if 0.01 0.5 () (1 ) / / (1 ) / 3 if 0 0.01 pp p p p ppp p p p           EM_2_2019.indd 71EM_2_2019.indd 71 19.6.2019 15:11:0419.6.2019 15:11:04 72 2019, XXII, 2 Business Administration and Management by this iterative way was used in the subsequent simulations. 2. Aspects of the Control Chart Performance Assessment When the performance of different types of control charts is to be compared, a suitable criterion must be chosen. It is common to compare the effi ciency of charts that have roughly equal ANOS0 (see e.g. Reynolds & Stoumbos, 1999; Wu et al., 2000). In the Shewhart-type charts, the probability of a signal at the moment when a nonconforming unit has been observed depends only on the length of the current run, however, the value of the CUSUM statistic is affected by the past history of the process. Therefore the zero-state and steady-state analysis is distinguished at the CUSUM chart if ANOS is examined, see e.g. (Szarka & Woodall, 2011; 2012). In the zero-state analysis a shift is supposed to be present when the process monitoring begins. The steady-state analysis is based on the assumption that a large number of conforming units were inspected before a shift occurred. In addition, the moment of the process shift (the change point) affects ANOS in charts based on the conforming counts. The steadystate fi xed-shift model assumes that the shift occurs immediately after a nonconforming item is found. Wu and Spedding (1999) were the fi rst to point out the problem in the CCC chart and considered the alternative steady-state randomshift model, when a shift may occur at any time during monitoring. The time to the change point is assumed to follow an exponential distribution (Wu et al., 2000) or a uniform distribution (Stoumbos & Reynolds, 1997). The random shift model is further considered in Wu et al. (2001) or Bourke (2006). A number of comparative studies were performed. In Bourke (1991) the RL2 chart and the geometric CUSUM were compared and the zero-state ANOS was determined. The in-control fraction nonconforming p0 was 0.01 and the ANOS curves showing dependence on the fraction nonconforming p up to 15p0 were displayed for a chosen design of the geometric CUSUM chart. The CCC-3 chart and the geometric CUSUM chart (and np chart) with the random-shift model were considered in Wu et al. (2000). Small in-control p0 between 0.0001 and 0.00065 were assumed, the specifi ed p1 for the upward shift detection ranged from p0 to 10p0 and only ANOS(p1) was tabulated. Chang and Gan (2001) compared the CCC chart and the geometric, Bernoulli and binomial CUSUM charts. The zero-state and the case p0 = 0.0001, p1 = 0.0003 were considered. The ANOS values were tabulated for p up to 10,000 p0. The RL2 (or CCC-2) chart and the geometric CUSUM (apart from the Poisson CUSUM chart and the optimal np chart with curtailment) were compared in (Wu, 2006). The steadystate analysis with the random-shift model was considered for p0 ranging from 0.005 to 0.02, p1 equal to 3p0 or 5p0, and ANOS0 equal to 10/p0 and 100/p0. More extensive comparisons of the moving sums RL2 and the geometric CUSUM (among others) can be found in Bourke (2008). Here both the initial-state and steady-state random-shift analysis were performed. The p0 values were between 0.002 and 0.01. For some specifi ed values of ANOS0 expressed as multiples of 1/p0 and ranging between 1,000 and 50,000 the ANOS(p1) was given, p1 values were between 3p0 and 20p0. In most works the geometric CUSUM chart was considered to be the best method that is capable to detect an increase in p more quickly than other charts, except for very large shifts, when the ANOS in the CCC-r charts was smaller. The conclusions in Wu et al. (2000) were rather different; the higher quickness of the SCRL (CCC-3) chart was sometimes manifested even for smaller shifts. The difference may be caused by the fact that both the downward and upward shifts were considered in the study. As the surprising fi nding concerned the value of p0 similar to ours, the comparison with our results may be of interest. 2.1 Estimation of ANOS for CCC-CUSUM Chart Three scenarios mentioned above were considered: zero state, fi xed-shift steady state, and random-shift steady state. The shift of the fraction nonconforming was represented by the change from p0 to p. The moment of the change depended on a scenario. a) Zero state The process shift occurred immediately. All observations came from the process with the fraction nonconforming p, i.e. all generated observations followed the geometric distribution G(p). EM_2_2019.indd 72EM_2_2019.indd 72 19.6.2019 15:11:0519.6.2019 15:11:05 73 2, XXII, 2019 Business Administration and Management b) Steady state, fi xed shift Sometimes the necessity of a large number of in-control observations is emphasized when the steady-state model is to be analyzed. For example, Szarka and Woodall (2012) use 5,000 observations from G(p0) before a shift occurs. This requirement seems to be worthless with respect to the random in-control variation of the CUSUM statistic. Here another method was used. The number of observations from G(p0) was much smaller and not fi xed, but it was restricted by the limit of 50,000 Bernoulli observations. Geometric observations Yi from G(p0) were generated and at the same time, their cumulative sum was recorded. As soon as the cumulative sum exceeded the predetermined value of 50,000 for some i = I, the YI value was cancelled and the change point was set immediately after I  1 observations from G(p0). All subsequent observations were from G(p). c) Steady state, random shift Now the change point was set after a fi xed number of 50,000 inspected items, i.e. randomly between two subsequent observed nonconforming items. Again, observations from G(p0) were generated. As soon as their cumulative sum exceeded the value of 50,000 for some i = I, only a part ZI of this YI obtained as a difference between 50,000 and the prior cumulative sum was taken into account. Then instead of YI the sum ZI + Yi+1 was considered, where Yi+1 came from G(p). All subsequent observations were from G(p). In both cases b) and c) false signals before the process shift were excluded; when the CUSUM statistic exceeded the lower limit H, the current observation was cancelled and another one was generated. 2.2 Calculating ANOS for CCC-r Chart The values of ANOS for the CCC-r charts were determined according to the formula ()  r ANOS pFLCL (14) where F(LCL) can be set using equation (9). The theoretical value ANOS0 for p = p0 is 0 0   r ANOS p (15) In the case of the steady-state randomshift model, ANOS for CCC-r was estimated using simulation as follows: the run lengths Yi from G(p0) and G(p) were generated as in c) of the previous section 2.1 but now sums of three (or two) subsequent observations were calculated and after each new triple (pair) the sum was compared with the corresponding lower control limit LCL. The cumulative sum of run lengths NOSi was calculated as before and as soon as the current sum of three (two) run lengths decreased under LCL for some i = I, the simulation was interrupted and NOSi was stored. To arrange the same sequence of observations as in 2c), the part of the former procedure for the geometric CUSUM with cancelling observations that led to false alarms was used. 3. Comparative Statistical Analysis Based on the experience with the manufacturing process which is described in detail in Section 4, the in-control p0 = 0.0002 and the out-of-control p1 = 5p0 were chosen. Three values of the risk α were considered: 0.0027, 0.005, and 0.01. The semi-economic model by Brodecká (2013) resulted in the CCC-3 chart for α = 0.0027 and the CCC-2 chart for α equal to 0.005 or 0.01. The lower limits for both types are displayed in Tab. 1. The fraction nonconforming p ranged from 0.0004 to 0.003. The values of ANOS for all charts are shown in Tab. 2, 3 and 4 and in Fig. 1, 2 and 3. CCC-r Geometric CUSUM αr LCL HGH 0.0027 3 1,354 -7,171 -7,179 0.0050 2 518 -5,901 -5,904 0.0100 2 744 -5,044 -5,061 Source: own Tab. 1: Lower limits of CCC-r chart and geometric CUSUM chart EM_2_2019.indd 73EM_2_2019.indd 73 19.6.2019 15:11:0519.6.2019 15:11:05 74 2019, XXII, 2 Business Administration and Management CUSUM CCC-3 pzero fi xed random traditional random 0.0002 5,557,857 5,557,857 5,557,857 5,566,358 5,557,857 0.0004 116,645 115,474 117,576 424,203 245,093 0.0006 24,954 24,727 25,981 101,870 64,965 0.0008 11,867 11,675 12,637 39,066 26,577 0.0010 7,624 7,482 8,204 19,332 14,328 0.0012 5,561 5,446 6,076 11,226 9,088 0.0014 4,390 4,286 4,837 7,271 6,379 0.0016 3,613 3,519 3,994 5,097 4,885 0.0018 3,075 2,996 3,419 3,791 3,904 0.0020 2,674 2,606 2,979 2,952 3,201 0.0022 2,366 2,308 2,643 2,385 2,726 0.0024 2,123 2,073 2,365 1,984 2,370 0.0026 1,923 1,873 2,146 1,690 2,124 0.0028 1,757 1,714 1,972 1,468 1,912 0.0030 1,618 1,580 1,820 1,297 1,743 Source: own CUSUM CCC-2 p zero fi xed random traditional random 0.0002 1,999,842 1,999 842 1,999,842 2,006,896 1,999,842 0.0004 77,281 77,227 79,117 268,458 273,083 0.0006 19,729 19,481 20,728 85,071 86,770 0.0008 9,850 9,707 10,639 38,361 40,558 0.0010 6,372 6,197 6,996 20,980 22,525 0.0012 4,694 4,554 5,197 12,961 14,153 0.0014 3,702 3,594 4,133 8,708 9,781 0.0016 3,053 2,972 3,442 6,220 7,189 0.0018 2,595 2,538 2,940 4,655 5,480 0.0020 2,260 2,208 2,572 3,613 4,308 0.0022 2,003 1,955 2,282 2,889 3,527 0.0024 1,798 1,756 2,059 2,366 2,971 0.0026 1,639 1,596 1,881 1,978 2,552 0.0028 1,503 1,460 1,725 1,682 2,216 0.0030 1,388 1,348 1,596 1,451 1,933 Source: own Tab. 2: ANOS for the geometric CUSUM (three models of the shift occurrence) and the CCC-3 chart (formula (14) and simulation), α = 0.0027 Tab. 3: ANOS for the geometric CUSUM (three models of the shift occurrence) and the CCC-2 chart (formula (14) and simulation), α = 0.005 EM_2_2019.indd 74EM_2_2019.indd 74 19.6.2019 15:11:0519.6.2019 15:11:05 75 2, XXII, 2019 Business Administration and Management CUSUM CCC-2 p zero fi xed random traditional random 0.0002 1,000,268 1,000,268 1,000,268 1,000,511 1,000,268 0.0004 57,341 57,770 59,423 137,798 140,892 0.0006 16,510 16,264 17,619 44,930 47,776 0.0008 8,482 8,280 9,270 20,832 22,720 0.0010 5,558 5,425 6,190 11,707 13,244 0.0012 4,092 3,987 4,608 7,426 8,661 0.0014 3,223 3,148 3,676 5,120 6,101 0.0016 2,662 2,599 3,050 3,750 4,632 0.0018 2,249 2,202 2,625 2,875 3,699 0.0020 1,958 1,919 2,297 2,285 3,028 0.0022 1,728 1,694 2,039 1,869 2,527 0.0024 1,545 1,508 1,830 1,565 2,157 0.0026 1,394 1,363 1,665 1,336 1,894 0.0028 1,269 1,242 1,524 1,160 1,676 0.0030 1,168 1,142 1,406 1,021 1,504 Source: own Tab. 4: ANOS for the geometric CUSUM (three models of the shift occurrence) and the CCC-2 chart (formula (14) and simulation), α = 0.01 Fig. 1: ANOS curves for the CCC-3 and the geometric CUSUM chart, α = 0.0027 Source: own EM_2_2019.indd 75EM_2_2019.indd 75 19.6.2019 15:11:0519.6.2019 15:11:05 76 2019, XXII, 2 Business Administration and Management Fig. 2: ANOS curves for the CCC-2 and the geometric CUSUM chart, α = 0.005 Source: own Fig. 3: ANOS curves for the CCC-2 and the geometric CUSUM chart, α = 0.01 Source: own EM_2_2019.indd 76EM_2_2019.indd 76 19.6.2019 15:11:0619.6.2019 15:11:06