Phase II monitoring of auto-correlated linear profiles using linear mixed model

In many circumstances, the quality of a process or product is best characterized by a given mathematical function between a response variable and one or more explanatory variables that is typically referred to as profile. There are some investigations to monitor auto-correlated linear and nonlinear profiles in recent years. In the present paper, we use the linear mixed models to account autocorrelation within observations which is gathered on phase II of the monitoring process. We undertake that the structure of correlated linear profiles simultaneously has both random and fixed effects. The work enhanced a Hotelling’s T statistic, a multivariate exponential weighted moving average (MEWMA), and a multivariate cumulative sum (MCUSUM) control charts to monitor process. We also compared their performances, in terms of average run length criterion, and designated that the proposed control charts schemes could effectively act in detecting shifts in process parameters. Finally, the results are applied on a real case study in an agricultural field.


Introduction
Control charts are used to detect anomalies in the processes. They are most often used to monitor production-related processes. In many business-related processes, the quality of a process or product can be characterized by a relationship between a response variable and one or more explanatory variables which is referred to as profile. The purpose of the analyzing of profile in phase I is to determine the stability of the process and estimate parameters; however, in phase II, analyzers are interested in rapidly detecting the significant shifts in the

A r c h i v e o f S I D www.SID.ir
process parameters. Phase I analysis of simple linear profiles has been investigated by a number of authors such as Stover and Brill (1998), Kang and Albin (2000), Kim et al. (2003) and Woodall 2004, 2007). Many authors including Kang and Albin (2000), Kim et al. (2003), Gupta et al. (2006), Zou et al. (2006), Saghaei et al. (2009), and Mahmoud et al. (2009) have investigated phase II monitoring of simple linear profiles. Noorossana et al. (2010a, b) investigated monitoring of multivariate simple linear profiles on phase II. Zou et al. (2007) and Kazemzadeh et al. (2009a, b) considered cases when the profiles can be characterized by multiple and polynomial regression models respectively. Mahmoud (2008) considered phase I monitoring of multiple linear profiles, and Kazemzadeh et al. (2008) proposed three methods for monitoring the kth-order polynomial profile in phase I. Ding et al. (2006), Moguerza et al. (2007), Williams et al. (2007), and Vaghefi et al. (2009) investigated nonlinear profiles. In these studies, it is implicitly assumed that the error terms within or between profiles is independently and identically normally distributed; however in some cases, these assumptions can be violated. Noorossana et al. (2010a, b) analyzed the effects of non-normality on the monitoring of simple linear profiles. Noorossana et al. (2008) and Kazemzadeh et al. (2009a, b) investigated autocorrelation between successive simple linear and polynomial profiles respectively. Soleimani et al. (2009) proposed a transformation to eliminate the autocorrelation between observations within a simple linear profile in phase II. Jensen et al. (2008) proposed two T 2 control charts based on linear mixed model (LMM) to account for the autocorrelation within linear profiles in phase I. They concluded that the linear mixed model is superior to the least square approach for unbalanced or missing data, especially when the number of observation within a profile is small and the correlation is weak. Jensen and Birch (2009) used nonlinear mixed model to account correlation within nonlinear profiles. Qie et al. (2010) investigated nonparametric profile monitoring with arbitrary design using mixed models. They proposed a control chart that combines the exponentially weighted moving average control chart based on local linear kernel smoothing and a nonparametric regression test under the assumption that observations within and between individual profiles are independent of each other.
The present study acts as an extension of the work of Jensen et al. (2008) in applying a linear mixed model on the presence of autocorrelation within linear profiles on phase I control chart applications; conversely, our focus is on phase II of profile monitoring in which one could use the proposed control charts to detect any departures from the given profile parameters.
The remainder of the paper is organized as follows. In 'Linear mixed model' section, the LMM is mathematically presented. In the 'Proposed methods' section, our methods including three modified multivariate control charts namely Hotelling T 2 (Hotelling 1947), multivariate exponential weighted moving average (MEWMA) and a multivariate cumulative sum control charts (MCUSUM) are illustrated. In 'Simulation studies'section, the results of simulation study to evaluate the performance of the methods are presented. In addition, a case study from an agriculture field is investigated on the section 'Case study.' The final section closes with concluding remarks.

Linear mixed model
Linear mixed models (Laird and Ware 1982) are popular for analysis of longitudinal data. A linear mixed model contains fixed and random effects and is linear in these effects. This model allows us to account autocorrelation within profiles. In matrix notation, a mixed model can be represented as , and X and Z are matrices of regressors relating the observations y to β and b.
In the 1950s, Charles Roy Henderson provided the best linear unbiased estimate (BLUE) of fixed effects and best linear unbiased predictions (BLUP) of random effects. Subsequently, mixed modeling has become a major area of statistical research, including work on the computation of maximum likelihood estimates, nonlinear mixed effect models, missing data in mixed effects models, and Bayesian estimation of mixed effects models (West et al. 2007).
Henderson's 'mixed model equations' (MME) are (Robinson 1991) as follows: The solutions to the MME, β , and b are BLUEs and BLUPs for β and b, respectively.
Mixed models require somewhat sophisticated computing algorithms to fit. Solutions to the MME are obtained by methods similar to those used for linear least squares. For complicated models and large datasets, iterative methods may be needed.
In profile monitoring, one could suppose that the jth response follows a LMM; therefore, where X i is a (n j × p) matrix of regressors, and Z j is a (n j × q) matrix associated with random effects. β is a (p × 1) vector of fixed effects, and y j is the (n j × 1) response vector for the jth profile. The coefficient vector of the random effect terms is b j ~ MN (0,D), and D is assumed to be a diagonal matrix; thus, the random effects are assumed not to be correlated with each other. In addition, it is assumed that and ε j is (n j × 1) vector of errors where If the errors are assumed to be independent, 2 j σ = R I , but correlated, the functional structure for the error terms may be used.
As noted before, it is considered that β is an estimator of β, and j b is a predictor of b j , then j j = y X β is the population average, and j j j j + = y X β Z b is the profile specific prediction; so if D and R j are known, then it can be shown as follows: (Schabenberger and Pierce 2002).

Proposed methods
In this paper, we propose a linear mixed model approach for accounting the correlation within linear profiles in phase II. It is assumed that profiles are correlated based on first-order autoregressive (AR(1)) structure.
If the errors follow an autocorrelated structure such as an AR (1) It is assumed that for the jth sample collected over time, our observations are (X i ,y ij ), i = 1,2,…,n and j = 1,2,…,m. We considered the case that all the fixed effects have a corresponding random effect, ( ) j j = X Z . If the process is in control, the problem can be formulated as follows: where ε ij are the correlated error terms and a ij are white noises as a ij ~ N(0,σ 2 ). The β 0 , β 1 ,…,β p − 1 are fixed effects that are the same for all profiles. The b 0j ,b 1j ,…,b p − 1j are random effects for the jth profile and they are normal random variables with zero mean and variance of 2 2 , respectively, which are not to be correlated with each other and also not to be correlated with the errors. The x values are fixed and constant from profile to profile. In this article we especially focused on phase II of the monitoring process, so all profile's parameters, process variance, and correlation coefficient are known in phase I. Accordingly, we utilized the modified Jensen et al. (2008) approach to monitor autocorrelation on phase II.

The Hotelling's T 2 statistic control chart
As a first proposed control chart, we use T 2 statistic to monitor the fixed effects for each sample. This statistic is given by where ( ) ( ) X V X X V y and β 0 denote the in-control value of β.
In Equation 9 the variance covariance matrix of fixed effects is ( ) The upper control limit, UCL, is chosen to achieve a specified in control average run length (ARL).

The MEWMA control chart
Our second proposed control chart is based on MEWMA for monitoring the vector of j β . Here the MEWMA statistics is as follows: where z 0 = 0 and θ(0 < θ < 1) is the smoothing parameter. Therefore, the chart statistic denotes by MEWMA j is given by This control chart gives a signal when EWMA j > UCL, where (UCL > 0) is chosen to achieve a specified in control ARL.

The MCUSUM control chart
The third suggested method is based on the MCUSUM control chart. In this method, the statistic is given by

A r c h i v e o f S I D www.SID.ir
∑ c s β β s β β and s 0 = 0 and k is a selected constant.
The estimator of variance covariance matrix is The chart gives a signal if ( )

Simulation studies
To show the performance of the proposed methods, we considered the underlying linear profile as Equation 14: and a ij ~ N(0,1),b 0j ~ N(0,.1),b 1j ~ N(0,.1). In our simulation investigation, we considered three significant different autocorrelation coefficients: a ρ = 0.1 to designate a weak type autocorrelated process, intermediate autocorrelation by ρ = 0.5, and strong autocorrelation by ρ = 0.9. The in-control ARL is roughly set equal to 200 and the ARL values were evaluated through 10,000 simulation replications under different shifts in intercept, slope, and errors (standard deviation). For MEWMA control chart, the smoothing parameter θ is chosen to be 0.2. As a general rule, to design MCUSUM control chart with the k approach, one chooses k to be half of the delta shift which is the amount of shift in the process that we wish to detect, expressed as a multiple of the standard deviation of the data points. Accordingly, we set k equal to 0.5. UCLs of control charts are designed to achieve a specified in control ARL of 200. The simulated UCLs for each proposed control chart are shown in Table 1. The three proposed control charts are compared on different scenarios of the example in terms of ARL, and the calculated amounts for the different changes in the intercept is shown in Table 2. According to the Table 2, under σ shift in the intercept, when autocorrelation is weak (ρ = 0.1), the MEWMA method performs relatively similar to MCUSUM control chart, and they also have better performance for detecting the small, moderate, and large shifts than the T 2 control chart. In the intermediate and strong autocorrelation circumstance (ρ = 0.5) and (ρ = 0.9), MCUSUM performs uniformly better than the other two methods. Moreover, MEWMA uniformly performs better than T 2 control chart. Figure 1 presented the derivative ARL under different shifts in intercept when autocorrelation is different in three levels. Table 3 shows the simulation results under different shifts in slope.  From Table 3, under βσ shift in slope, while the autocorrelation is weak (ρ = 0.1), the proposed MCUSUM method uniformly performs better than MEWMA method. Also, MEWMA performs consistently better than T 2 method. In addition, similar results are obtained when the autocorrelation is intermediate (ρ = 0.5). Once the amount of autocorrelation coefficient is high, MCUSUM and MEWMA methods perform uniformly better than the T 2 method and also, MCUSUM method performs relatively similar to MEWMA method. Figure 2 illustrates ARL under different shifts in slope once autocorrelation be changed in the aforementioned levels. Next comparisons of the proposed three control charts in terms of ARL under δσ shift in the standard deviation followed. Table 4 shows that the proposed T 2 chart performs significantly better than MEWMA and MCUSUM charts in different amount of correlation coefficients. In addition for strong and intermediate autocorrelation condition, MEWMA and MCUSUM

A r c h i v e o f S I D www.SID.ir
have similar manners and when the autocorrelation is weak, MEWMA relatively achieves better performance. Derivative ARL under different shifts of standard deviation is presented in Figure 3 when autocorrelation is different.  Based on the simulation results, it is evident that the proposed MEWMA and MCUSUM methods act relatively better than the T 2 chart in detecting shift in the parameters of profile; conversely, the proposed T 2 chart performs better than the MEWMA and MCUSUM in detecting shift in the variation.

Case study
Consider the case study carried out by Schabenberger and Pierce (2002). It was a real data set from ten apple trees which 25 apples are randomly chosen on each tree. Their focus was on the analysis of the apples in the largest size, with initial diameters exceeded 2.75 in. Totally there were 80 apples in aspiration size. Diameters of the apples were recorded in every 2 weeks during 12 weeks. Figure 4 shows 16 diameters out of 80 apples in the time domain. In their investigation, functional profile between time and diameter considered as quality characteristic that needs to be monitored over time. Schabenberger and Pierce (2002) and also later Soleimani et al. (2009) modeled such correlation between observations by a first-order autoregressive model of AR(1). Based on the preceding analysis, the following statements hold a linear mixed model equation for the declared case study: Consequences of simulation run in the previous section leads us to use MEWMA and MCUSUM which have relatively similar performance on detecting shift in the profile parameters rather than the T 2 method. Hence, the proposed MEWMA control chart was applied in monitoring the linear profile. The smoothing constant (θ) is set equal to 0.2. In order to achieve an in control ARL of 200, the upper control limit is set equal to 7 based on 10,000 simulation runs. In order to examine performance of the control chart, six random samples from the in control simple linear profile are initially generated. Formerly, three random samples are generated to show an out-of-control condition under the intercept shift coefficient of 0.6. Figure 5 illustrates sensitivity of the MEWMA control chart based on our proposed method which temperately depicts quick signal.

Concluding remarks
We have studied the sensitivity of three multivariate control charts to detect one-step permanent shift in any parameters of a mixed model linear profile. Our specially designed MCUSUM, MEWMA, and T 2 control charts were also studied as competitors of each other to depict shifts in intercept and slope parameters and also process variation while first-order autoregressive model describes correlations within observations.
The performances of the methods were compared in terms of average run length criteria. Table 5 shows the summarized results. The following summary recommendations are made: 1 The proposed approach has good performance across the range of possible shifts and it can be used in phase II of linear profile monitoring on the presence of autocorrelation within observations. 2 The anticipated MEWMA and MCUSUM methods almost uniformly perform better efficiency than the T 2 Hotelling control chart under different step shifts in the intercept and slope parameters of linear profile.