Improving the Quality of Response Surface Analysis of an Experiment for Coffee-Supplemented Milk Beverage: I. Data Screening at the Center Point and Maximum Possible R-Square

How the elimination of outliers at center runs can improve the statistical model will be illustrated through re-analysis of a dataset described in the article entitled “Optimization of Manufacturing Conditions for Improving Storage Stability of Coffee-Supplemented Milk Beverage Using Response Surface Methodology” authored by Ahn et al. (2017). This article modeled two responses, using two factors. The two responses are the particle size and the zeta-potential of milk beverage. The two factors are the speed of primary homogenization (unit: rpm) and the concentration of emulsifier (unit: %), for which X₁ and X₂ are used as the coded factor. X₁=−1, 0, and 1 correspond to the speed of primary homogenization=5,000 rpm, 10,000 rpm, and 15,000 rpm, respectively, and X₂=−1, 0, and 1 correspond to the concentration of emulsifier=0.1%, 0.2%, and 0.3%, respectively. Among the two responses, the first response, which is the particle size, had an extreme outlier at the center point. Thus, this response is used as the Y variable in this research article.

The dataset to be re-analyzed is shown in Table 1. In this dataset, the experimental design is the CCD for two factors with an axial value of 1 and five center runs. Using this design, we can fit to the data statistical models including a full second-order model.

Data were analyzed using SAS software. SAS/STAT (2013) procedures were used for statistical modeling. Graphs were drawn using Minitab (2017).

First, let us look at the responses from center runs to find outliers. In Table 1A, which is in Ahn et al. (2017), displays observations according to the run number, and Table 1B displays observations according to the standard order. Which display makes it easier to find out the outliers at the center point? Obviously, the second display! Thus, we suggest presenting the data in the form of Table 1B before data modeling.

In Table 1B, runs with standard order numbers 9 through 13 are the center runs, and the responses from these runs are 260.5, 186.433, 181.933, 175.633, and 180.333. Among these values, obviously 260.5 is an outlier, and if the reason for this observation is not an error in measuring or typing, this value should be eliminated. 3D scatterplots in Fig. 1, which graphically compares the data containing an outlier and the data without an outlier, say that at the center point, 260.5 is an extreme outlier.

Fig. 1. 3D scatterplots of original data containing an outlier at the center point (A) and data after outlier elimination (B).

Download Original Figure

Now, let us see what happens if statistical models are fitted to the dataset that contains this extreme outlier.

First, let us fit to the data the second-order polynomial regression model containing 2 linear, 2 quadratic, and 1 interaction terms by using the RSREG procedure of SAS/STAT. Results of analysis of variance for this model are shown in Table 2A.

The model in Table 2A has a very poor fit; its model p-value=0.9460 is so large that it is close to 1, whereas the p-value of an acceptable model is no larger than 0.05, and its r²=0.1322 is so small!

Now, let us find the maximum r² that can be obtained through the statistical modeling of this original dataset. The r² of the one-way ANOVA (analysis of variance) model on 9 design points, which are designated in Table 1B, is such an r². This one-way ANOVA model is the fullest model among the models that can be fitted to this dataset. Table 2B displays the results from this modeling.

The model in Table 2B has also a poor fit; its model p-value=0.9622 is so high that it is near 1, and its r²=0.3173 is still too low!

Thus, there is no way for the model displayed in Table 2A to be improved. Now, our remedy is to remove the observation with run number 2 that has the outlier 260.5 at the center point. Let us try it. Table 1C displays the dataset from which this outlier has been deleted.

Now, let us fit to the data in Table 1C the second-order polynomial regression model containing 2 linear, 2 quadratic, and 1 interaction terms by using RSREG procedure of SAS/STAT. Results of analysis of variance for this model are shown in Table 1C.

The model in Table 2C is better than that in Table 2A, but still unsatisfactory; its model p-value=0.5962 is still large, and its r²=0.3952 is still small! This implies that the second-order model is insufficient to represent this dataset. But, it may be possible for the model to be improved enough to well explain data. To check this possibility, let us also find the maximum r² that can be obtained through the statistical modeling of this reduced dataset. The r² of the one-way ANOVA model on 9 design points, which are designated in Table 1C, is such an r². This one-way ANOVA model is the fullest model among the models that can be fitted to this dataset. Table 2D displays the results from this modeling.

Now, the model displayed in Table 2D is satisfactory in its ability to explain the data; its model p-value=0.0230 is lower than 0.05, meeting the criterion on the p-value for an acceptable model, and its r²=0.9762 is large enough!

This implies that the model in Table 2C can be improved by adding some proper terms to the model. Now, such an augmentation of the model to improve it in the case of a cubic central composite design, which is our case, is another topic for research. This research has been done, and a satisfactory model has been found! However, to focus on one topic at a time, the presentation of the result of this research will be given in the next article.