The analysis of residuals variation and outliers to obtain robust response surface

In this paper, the main idea is to compute the robust regression model, derived by experimentation, in order to achieve a model with minimum effects of outliers and fixed variation among different experimental runs. Both outliers and nonequalit y of residual variation can affect the response surface parameter estimation. The common way to estimate the regression model coefficients is the ordinary least squares m ethod. The weakness of this method is its sensitivity to outliers and specific residual behavi or, so we pursue the modified robust method to solve this problem. Many papers have proposed different robus t methods to decrease the effect of outliers, but trends in residual behaviors pose another important issue that should be taken into account. The trends in residuals can cause faul ty estimations and thus faulty future decisions and outcomes, so in this paper, an iterat ive weighting method is used to modify both the outliers and the residuals that follow abnormal t ends in variation, like descending or ascending trends, so they will have less effect on the coefficient estimation. Finally, a numerical example illustrates the proposed approach.


Introduction
In many cases, especially in experimental results, some of the data are wrong and should be treated as outliers. These points, which may occur because of operator reading faults and the like, may have a confusing effect on the total interpretation of the results. A common method of explaining and analyzing the results of experiments is by response surface design. This term is used for a regression equation that shows the whole behavior of the control variables, the nuisance factors, and the response or responses. We can use the estimated function to predict the response to changes in the values of specific controllable factors. After determining an experimental design and performing experiments, the next steps are generally statistical analysis and then the selection of values for the input variables so as to optimize the output. This can be done by fitting a regression model between the controllable factors and the response variables. Future interpretations are based on this regression model, so the exact model is very important and may affect the optimization stage. This model is generally

A r c h i v e o f S I D www.SID.ir
constructed by the ordinary least squares (OLS) method. But basic OLS is very sensitive to outliers, and they may have an inordinate effect on the ultimate conclusion. So a robust method or a modified OLS should be used for decreasing the outliers' sensitivity.
Our goal in this study is to decrease the destructive effect of outliers. In order to do so, at the first stage, the robust regression or modified regression model should be computed. Then the trend of residuals in response surface design is another aspect which should be considered. The trend behaviors among residuals, both descending and ascending trends, can cause faulty interpretations. However, there are some assumptions in estimating the regression coefficients that should not be violated. It seems that by decreasing the effects of some residuals in the coefficient estimation stage, the initial assumptions can be satisfied, and moreover, because this decreases the overall variability, the robustness of the model will be increased. The main purpose of this paper is to decrease the effects of outliers that violate the variance equality test of residuals. An example of such trends and abnormal behavior is shown in Figure 1. As mentioned before, the OLS method is very sensitive to outliers. To diminish the effect of these points, some alternative methods of model fitting, such as least absolute deviations and other robust approaches that simplify the task of outlier identification by weighting the large residuals, are used instead of OLS. Response surfaces have been studied by many researchers, and many approaches have been proposed either to obtain efficient response surfaces or optimize the response surface by different models. Hejazi et al. (2010) proposed a novel approach based on goal programming to find the best combination of factors to optimize multi-response-multi-covariate surfaces by considering location and dispersion effects. Kazemzadeh et al. (2008) proposed a method to optimize multi-response surfaces based on a goal programming method. Robust regression approaches have also been surveyed by many researchers. Huber (Bertsimas and Shioda 2007) proposed M-estimator methods to obtain robust regression. Morgenthaler and Schumacher (1999) discussed robust response surfaces in chemistry based on design of experiment. Because of the weakness of the previous methods in compensating for outliers, the redescending M-estimators (also named GM-estimators) were proposed by Andrews et al. (1972), which are able to reject extreme outliers entirely. Hund et al. (2002) presented various methods of outlier detection and evaluated robustness tests with different experimental designs. Bickela and Frühwirthb (2006) compared different robust estimators with their applications. The M-and GMestimators work by an iterative procedure. As a consequence, several authors (e.g., Cummins and Andrews 1995) have called these estimators as iteratively reweighted least squares, or IRLS methods. Ortiz et al. (2006) discussed some of the robust methods used for robust regression in analytical chemistry. To obtain a more efficient and yet robust method, Siegel (1982) proposed the repeated median estimator. Also another useful robust method is the least median squared (LMS) method proposed by Rousseeuw (1984). Massart et al. (1986) showed the advantages of its use in chemical analysis. The other useful method is the least trimmed squared (LTS) that was proposed by Rousseeuw and Leroy (1987). Nguyena and Welsch (2010) studied outlier detection and proposed a new least trimmed squares approximation. Both the LMS and the LTS are defined by minimizing a robust measure of the scatter of the residuals. Generalizing this, Rousseeuw and Yohai (1984) introduced Sestimators which are significantly more efficient than the previous estimators. A more recent suggestion is the constrained M-estimates, or CM, proposed by Mendes and Tyler (1995), which combines the good local properties of the M-estimates and the good global robustness properties of the S-estimates. A 'partial' version of the M-estimator based on the 'fair' ψ function and an appropriate weighting scheme was recently proposed by Serneels et al. (2005). The authors believed that the partial robust M-regression outperforms existing methods for robust partial least square regression. Bertsimas and Shioda (2007) presented mixed integer programming or MIP models for the classification and robust regression problems. Zioutas and Avramidis (2005) examined the effect of deleting outliers in the regression model obtained by mixed integer programming and compared the performance of this model with that of least squares, or LS, and LMS. Another new method in robust regression is the mixed linear model surveyed by Dornheim and Brazauskas (2011). (Pop and Sârbu 1996) proposed a new fuzzy regression algorithm to obtain robust models. Maronna et al. (2006) proposed many M-estimators using robust regression methods in both single response and multiple responses. Shahriari et al. (2011) proposed a novel two-step robust estimation of the process mean method based on M-estimator and their method is less sensitive to the presence of outliers. For better illustration of proposed method, the literature review has been classified in Table 1.  Huber (1981) ✓ ✓ Siegel (1982) ✓ ✓ Rousseeuw (1984) ✓ ✓ Massart et al. (1986) ✓ ✓ Rousseeuw and Leroy (1987) ✓ ✓ Cummins and Andrews (1995) ✓ ✓ Pop and Sârbu (1996)  ✓ ✓ Bickela and Frühwirthb (2006) ✓ ✓ Ortiz et al. (2006) ✓ ✓ ✓ Bertsimas and Shioda (2007) ✓ Nguyena and Welsch (2010) ✓ ✓ Dornheim and Brazauskas (2011) ✓ Proposed REVIM

✓ ✓ ✓
In this paper, a novel robust approach considering both outlier data and trends in residuals variations which do not violate the normality assumption is discussed.
This paper is organized as follows. The section 'Robust estimation of the coefficients by iterative weighting methods' presents the robust modification of the response surface by an iterative weighting procedure. The proposed method is defined in section 'Robust estimation of coefficients by testing equality of variations in specified intervals'. To illustrate the proposed method, a numerical example is presented in section 'Numerical example'. Finally, the last section is the 'Conclusion' of this paper.

Robust estimation of the coefficients by iterative weighting methods
To compensate for the effects of the outlier values, we can either remove the outlier data or modifying them by weighting the residuals. The first approach is not rational, so we choose to modify them in order to decrease the effect of outliers in the coefficient estimation stage. The proposed idea is as follows: In this equation, µ i is a function defined by unknown coefficients (β i ). For example, if µ 1 = β 1 + β 2 x 1 and x i are constants, the response y i can be obtained from the experimental results, and the regression model describes the relation between the variables and the expected values of the y i .
If all the measurements are good, then the OLS method provides a reasonable model and the coefficients are estimated by minimizing the following equation: The weights are not pre-assigned values because the quality of each y i is not known in advance. The reasonable values for the weights are based on the residuals defined by the following equation: To make the estimator invariant with respect to the scale of the residuals, the r i is divided by 's,' which is a robust estimation of the scale. The value of 's' is often taken to be equal to 1.4826 MAD, where MAD is the median of the absolute deviations of the residuals from their median and 1.4826 is a bias adjustment for the standard deviation under the normal distribution.
The weights should be inversely proportional to the value of the residuals, where c is a constant. The procedure is as follows: compute the first coefficients of the regression model, compute the residuals and weights, and then compute the new coefficients by the equation. This procedure can be repeated until a good solution is obtained, because the values of the coefficients and the values of the residuals and weights are different. This procedure is known as iterative weighting OLS. The procedure terminates when the change in the estimation from one iteration to the next is sufficiently small. This iterative method is good for modifying outliers, but the trends in residual behavior are not considered. As illustrated before, another approach, in addition to taking outliers into account, is the equality of variation between residuals which is the main idea of the rest of this paper.

Robust estimation of coefficients by testing equality of variations in specified intervals
First, normality assumption is checked. If the normality assumption is violated, this robust approach based on Huber function cannot be applied. The method proposed in this part begins by dividing the experiment runs into a intervals to examine the hypothesis of equality in variations in these intervals. The equality test used in this paper is Bartlett's test (Anderson and McLean 1974). The number of points in each interval can be chosen in the analysis stage. This stage satisfies one of the OLS hypotheses. If this parameter is small, the variation between points might be large and if the number of the points is large, the equality test of the variances may not be reliable, so this value should be determined rationally.

Bartlett's test
Although residual plots are frequently used to diagnose inequality of variance, statistical tests have also been proposed. One widely used procedure is Bartlett's test. The procedure involves computing a statistic whose sampling distribution is closely approximated by the chi-square distribution with a-1 degrees of freedom, where the a random samples are from independent normal populations. This statistic is defined as The hypothesis of equality of variances is rejected if , 1 x α α − is the upper α percentiles of the chi-square distribution with a-1 degrees of freedom.

Proposed robust approach (robust equal variances iterative method)
The following steps are proposed as an iterative method to decrease the effect of trends in the residuals and improve the robust estimation of coefficients. The proposed model is based on OLS method which should be modified.
First of all, our goal is that the residuals that violate the hypothesis of equality in variances should have less effect on the estimates of the coefficients of the regression model, so we should consider modifying these points to have equal variances. Therefore, the residuals derived by experiments are divided into a intervals, and then the variances of each interval are calculated and denoted by 2 i s . The next step is to test the equality of variances with Bartlett's test; if the result shows that the variances in a intervals do not have significant differences, this part of the procedure is stopped, whereas if the result shows that the variances in the a intervals differ significantly, the iterative weighting procedure is used to modify. As the next step, the critical q statistic in Bartlett's test, for which the hypothesis will be rejected, is computed. The critical q for Bartlett's test is denoted by q c , and in our case, c q q > . The ( ) 2 1 s is the ith rank-sorted variances of the intervals, in descending order. The ( ) ( ) 2 j i s is the variances of the points in ith ranked interval in the jth. In the proposed approach, in higher and higher iterations approach, the residual variances approach equality. The maximum feasible variance of each intervals for which Bartlett's test is fulfilled is denoted as 2 max i s .
Because the q value is greater than q c after its computation, it should be decreased iteratively until the both values are equal. To do this, we select the largest variance of all intervals. This value is decreased in each iteration. If the value of ( ) s in this decreasing procedure, all three variances decrease in the next iteration, and this procedure continues until q is equal to q c . All maximum feasible variances of intervals that satisfy the hypothesis of Bartlett's test are computed with this iterative procedure. Next, we consider two parallel lines l 1 :y = +ω and l 2 :y = −ω, with slope zero and parallel to the x-axis. We decrease the ω values iteratively and compute the variances of points between these two

A r c h i v e o f S I D www.SID.ir
lines. Until the variances of the points are equal to the maximum variance of intervals derived from the last step, this decreasing procedure continues. After that, the points outside these lines are weighted by function in (10).
The pseudo code of the proposed method illustrates this approach: 1. Divide residuals into a intervals.
2. Compute the variances of outliers in each interval, and denote these variances by 2 i s .
3. Sort the variances of each interval in descending order and denote them by ( ) 2 i s .

4.
Calculate the critical value of Bartlett's test in terms of confidence level and the number of intervals. 5. Compute the q value by formula in (7). 6. Compute the statistic of Bartlett's test by formula in (6). 7. Compute the q c by considering the critical value of Bartlett's test based on formula (6).

8.
Do while q ≥ q c , ( Determine s i 2 max as the maximum feasible variances of each intervals to satisfy the hypothesis of equal variances.

13.
Determine the weight of the outside residuals by function below (formula 10) 2 1 2 1 1 , After this step, by the robust iterative weighting method, the outliers are modified. After each iteration, the test of equality in variances is checked, and if the hypothesis is violated, the above mentioned method is applied. After the equality of variances is satisfied, the iterative weighting method continues, and this procedure continues as long as a minor change in the estimation occurs from each iteration to the next one. This procedure is called the robust equal variances iterative method procedure (REVIM).
The flowchart in Figure 2 illustrates this method.

Numerical example
This is a hypothetical numerical experiment. Suppose that we have an experiment containing one response variable and four explanatory control variables, each of which has three levels, and the objective of the study is to optimize the yield of a product. The data to be used are shown in Table 2. We want to explore the yield response surface by using a second-order regression model. A Box-Behnken design with 27 treatments is used for this experiment. The blocking is used to diminish the effect of nuisance factors, and the blocks are assigned, for example, to 3 days. The primary fitted response regression model is as follows: The normality assumptions are checked in this example and the results are given in Figure 3. The P value obtained by normality test is 0.414. This value shows that the residuals follow normal distribution.

A r c h i v e o f S I D www.SID.ir
The analysis of variance (ANOVA) results are shown in Table 3.  Figure 4 shows the residuals of the model in the order of runs. As shown in the figure, the residuals have a rough trend, and if we divide the runs into 3 equal intervals, each containing 9 runs, the second interval has larger variance. This can be proved by Bartlett's test. The result of Bartlett's test is illustrated in Figure 5 and Table 4.  In this case, we have three intervals, so a-1 = 2. If we consider the significance level 0.95, the critical statistic x 0.05,2 2 is equal to 5.99 and the test statistic = 9.09 is greater than 5.99, then the hypothesis is rejected. Therefore, we want to compute the maximum standard deviation of each interval that satisfies the hypothesis of Bartlett's test. By the proposed method, the maximum values of these standard variations are ordered as follows: s 1 = 2.07, s 2 = 4.95, s 3 = 4.95. Based on these values, we can compute the limits, l 1 = 8.7 and l 1 = −8.75. Two points 5 and 16 are outside these limits (circled in the figure), so the weighting procedure is applied and new coefficients are calculated by the robust weighting method. Figure 6 shows the method graphically. The residuals are computed by this procedure, and the hypothesis of equal variance in three intervals is satisfied. Bartlett's test is applied to the residuals obtained from the proposed method, and the results are given in Figure 7.

Figure 7 Proposed method residuals.
The value of Bartlett's test is 4.32, and the hypothesis of equal variances is not rejected. Figure 8 illustrates the result better.

Figure 8 Bartlett's test results after applying REVIM.
Therefore, the iterative weighting method based on these coefficients to modify the effect of outliers by Huber function with c = 2 is applied. This process is applied in each iteration, and if the equality test is not satisfied, the modification is applied. The final robust coefficients and residuals are presented in Tables 5 and 6.  So the residuals in final iteration show that residuals which hinted to be an outlier in first iteration are really outliers, and the model estimation is more accurate and close to the model with no outlier data that we call it actual model. These results can be compared with the results obtained by the model with no outlier data. The comparison shows that the proposed model is more precise and accurate than the OLS method in estimation of the regression coefficients by considering unequal variation between residuals. The results are given in Table 7.

Conclusions
A robust estimation of the response surface is the primary goal of this paper. To this end, the proposed method is defined, instead of the common ordinary least squares method of estimating coefficients of the response surface, to decrease the effects of two main causes of the imprecise estimation of coefficients, outliers, and trends in residuals. As the effect of trends in residuals should be taken into account, the proposed method simultaneously modifies the effects of trends and outliers. For each iteration, an equality test of residual variances is performed, and after this hypothesis is satisfied, the outliers are modified. A goal for future research may be to examine the weighting; instead of computing the distance of the residuals from base line after plotting the responses in a normal probability plot (NPP), the weighting function may be proportional to the distance of the response from the NPP regression line.     i v e o f S I  Figure 5 i v e o f S www.SID.ir Figure 6 i v e o f S www.SID.ir Figure 7 i v e o f S I  Figure 8