Introduction
Response surface methodology (RSM) refers to a set of statistical methods for the design of experiments, modeling of responses, and optimization of factor levels (Myers et al., 2009). RSM has extensively been used in experimental sciences, including the food science of animal resources. In RSM, a satisfactory modeling of responses is crucial for an accurate understanding of the response system and a precise optimization of explanatory factors. A model can be said to be satisfactory when: (1) the model is significant at the 5% level (the pvalue of the model≤0.05), (2) the lack of fit is nonsignificant at the 5% level (the pvalue of the lack of fit>0.05) or the model has no lack of fit, and (3) the adjusted rsquare is as large as at least 0.8 (≥0.8).
In RSM, after obtaining data from an experiment, we fit a standard secondorder model to the data and then check whether the above three criteria are met. If not all three criteria are met, the secondorder model is unsatisfactory. Subsequently, we fit a “balanced higherorder model” to the data and check whether the above three criteria are met. If still not all three criteria are met, the balanced higherorder model too is unsatisfactory, and we fit a “balanced highestorder model” to the data. This threestep sequential modeling can be summarized as follows:

Step 1: Fit a secondorder model to the data. If this model is satisfactory, modeling stops here. Otherwise, go to step 2.

Step 2: Fit a balanced higherorder model to the data. If this model is satisfactory, modeling stops here. Otherwise, go to step 3.

Step 3: Fit a balanced highestorder model to the data.
When each factor of the experimental design has five levels, a balanced higherorder model in step 2 is a thirdorder model, and a balanced highestorder model in step 3 is the fullest balanced model proposed by Rheem et al. (2017). These models used in Rheem et al. (2017) contain cubic terms in the form of X_{j}^{3}. When each factor of the design has three levels that are coded as –1, 0, and 1, these models cannot be used, as cubic terms cannot be entered into the model because they are the same as linear terms [(–1)^{3}=–1, 0^{3}=0, and 1^{3}=1, that is, X_{j}^{3}=X_{j})]. Thus, when each factor has three levels, the balanced higherorder and balanced highestorder models that have no cubic terms in them are required. The construction of such models will be explained later.
The present study proposes a threestep modeling strategy for experiments where each factor of the design has three levels, and subsequently optimizes factor levels for maximizing the response based on the balanced highestorder model. The study was conducted using the data from a threelevel, twofactor experiment on a coffeesupplemented milk beverage (Ahn et al., 2017; Rheem et al., 2019).
Materials and Methods
The construction of a balanced highestorder model will be explained through reanalysis of a dataset from the experimental data in the article titled “Improving the Quality of Response Surface Analysis of an Experiment for Coffeesupplemented Milk Beverage: II. Heterogeneous Thirdorder Models and Multiresponse Optimization” authored by Rheem et al. (2019). In that study, two threelevel factors were used in an experiment. Among the responses in that experiment, the Zetapotential, for which the objective of optimization is maximization, is the Y variable in this study. To illustrate the model for the data from a threelevel, twofactor design, the factors (the X_{1} and X_{2} variables) in this experiment and their coded and actual levels are presented in Table 1.
The dataset used for reanalysis was the data in Rheem et al. (2019) from which the first center run was deleted. The data in Rheem et al. (2019) was obtained through the screening (Rheem and Oh, 2019) of the original data in Ahn et al. (2017). The dataset that was analyzed in this study is presented in Table 2. In this dataset, the experimental design is the central composite design by Box and Wilson (1951) for two factors having –1, 0, and 1 as their coded levels. With this dataset, we performed a response surface analysis using balanced higherorder and balanced highestorder models.
Data were analyzed using the SAS software. SAS (2016b) procedures were used for the modeling of the response and the drawing of a contour plot. The optimization of factor levels was performed by SAS datastep programming. A threedimensional plot was generated using SAS (2016a).
Results and Discussion
First, the secondorder polynomial regression model containing two linear, two quadratic, and one interaction term was fitted to the data by using the RSREG procedure of SAS/STAT. Results of the analysis of variance for the secondorder model are shown in Table 3.
Test of lack of fit  

Source  Degrees of freedom  Sum of squares  Mean square  Fvalue  pvalue 
Lack of fit  3  37.67  12.56  140.69  0.0071 
Pure error  2  0.18  0.09 
Table 3 says that: (1) the model is nonsignificant at the 5% level (the pvalue of the model=0.3338>0.05), (2) the lack of fit is significant at the 5% level (the pvalue of the lack of fit=0.0071<0.05), and (3) the adjusted rsquare=0.1388<0.8; none of the three criteria were satisfied, and thus, this model was unsatisfactory. As the secondorder model had a poor fit, we proceeded to build a balanced higherorder model.
In a multifactor polynomial regression model, if each factor appears the same number of times in all the terms in the model, such a model is said to be “balanced.” In a balanced response surface model, each factor has an equal opportunity of being evaluated through the coefficients of the model terms that contain it.
When there are two factors, a secondorder model, which consists of X_{1}, X_{2}; X_{1}^{2}, X_{2}^{2}; and X_{1}X_{2}, is balanced, because for j=1 and 2, each X_{j} appears three times in all the terms (1 time in X_{1}, X_{2}; 2 times in X_{1}^{2}, X_{2}^{2}; and 1 time in X_{1}X_{2}) in the model. Thus, the balanced higherorder model consists of X_{1}, X_{2}; X_{1}^{2}, X_{2}^{2}; X_{1}X_{2}; X_{1}X_{2}^{2}, X_{1}^{2}X_{2}, where for j=1 and 2, each X_{j} appears 7 times (1 time in X_{1}, X_{2}, 2 times in X_{1}^{2}, X_{2}^{2}, and 1 time in X_{1}X_{2}, and 3 times in X_{1}X_{2}^{2}, X_{1}^{2}X_{2}). Table 4 shows the results of the analysis of variance for the balanced higherorder model.
Test of lack of fit  

Source  Degrees of freedom  Sum of squares  Mean square  Fvalue  pvalue 
Lack of fit  1  5.66  5.66  63.45  0.0154 
Pure error  2  0.18  0.09 
Table 4 says that: (1) the model is nonsignificant at the 5% level (the pvalue of the model=0.0841>0.05), (2) the lack of fit is significant at the 5% level (the pvalue of the lack of fit=0.0154<0.05), and (3) the adjusted rsquare=0.7785<0.8; none of the three criteria were satisfied, and thus, the balanced higherorder model was also unsatisfactory. Therefore, for the next step of model fitting, a balanced highestorder model was fitted to the data.
In Table 4, the lack of fit had one degree of freedom, implying that we can add one more term to the model. For the model to be balanced, this additional term needs to contain all of X_{1} and X_{2}. As the last terms in the model are X_{1}X_{2}^{2} and X_{1}^{2}X_{2}, the next term to be entered into the model should be X_{1}^{2}X_{2}^{2}. When we add this term to the model, the model becomes a balanced highestorder model, which has no lack of fit, because the number of coefficients in the model is the same as the number of design points in the experimental design. Design points are indicated in Table 2. Results of the analysis of variance for the balanced highestorder model are presented in Table 5.
Test of lack of fit  

Source  Degrees of freedom  Sum of squares  Mean square  Fvalue  pvalue 
Lack of fit  0  0  
Pure error  2  0.18  0.09 
Table 5 says that: (1) the model is significant at the 5% level (the pvalue of the model=0.0081<0.05), (2) the model has no lack of fit, and (3) the adjusted rsquare=0.9898>0.8; all three criteria were satisfied, and the rsquare was 0.9980, almost equal to 1. Therefore, the balanced highestorder model was satisfactory and accepted as our final model. With Ŷ denoting the predicted value of Y, this model can be specified as:
Table 6 shows the coefficient estimates and their SEs, tvalues, and pvalues of the balanced higherorder model.
The pvalues in Table 6 say that the X_{1}, X_{1}^{2}, X_{1}X_{2}, X_{1}X_{2}^{2}, X_{1}^{2}X_{2}, and X_{1}^{2}X_{2}^{2} terms, which all involve the X_{1} factor, are all significant at the 5% level, whereas the X_{2} and X_{2}^{2} terms, which involve only the X_{2} factor, are nonsignificant at the 5% level. From this, we see that the X_{1} factor plays a more important role than the X_{2} factor in the model.
According to Ahn et al. (2017), the objective of optimization for Y was maximization. Thus, using a search on a grid (Oh et al., 1995), we optimized the model for maximization in the form of the above equation. In this experiment, the bounds were –1≤X_{j}≤1 for j=1, 2. Thus, within these bounds, we conducted a search on a grid using the SAS data step programming. A search for the maximum on a grid was performed by calculating the Ŷ function over a grid for the values of X_{1} and X_{2}, with an increment of 0.01 within the bounds –1≤X_{j}≤1 for j=1, 2, and then sorting the calculated function values in descending order. The optimum point at which Ŷ was maximized was found this way and is presented in Table 7.
X_{1}  X_{2}  F_{1}=Speed of primary homogenization (rpm)  F_{2}=Concentration of emulsifier (%)  Predicted maximum of Y (=Zetapotential) 

−1  0.08  5,000  0.208  36.1515 
In Rheem et el. (2019), the predicted maximum Zetapotential based on their model was 35.2957, which is lower than 36.1515, our predicted maximum Zetapotential in Table 7. And, according to the data and the model in Rheem et el. (2019), for their predicted maximum (35.2957) of Zetapotential, their predicted optimum points of (X_{1}, X_{2}) in coded levels and (F_{1}, F_{2}) in actual levels are calculated as (X_{1}=–1, X_{2}=0.10) and (F_{1}=5.000, F_{1}=0.210), which are different from (X_{1}=–1, X_{2}=0.08) and (F_{1}=5.000, F_{1}=0.208), our predicted optimum points in coded and actual levels in Table 7. Thus, it can be said that our approach has predicted improved results of optimization, which are different from the predicted optimization results in the previous research. A validation experiment is needed after the optimum point is predicted.
For the X_{1} and X_{2} factors, the threedimensional (3D) response surface plot was drawn, with the vertical axis representing the predicted response and two horizontal axes indicating the coded levels of the X_{1} and X_{2} factors. Fig. 1 is such a plot.
The twodimensional contour plot of the response surface was also drawn, with two axes indicating two coded factors. Fig. 2 is such a plot.
These figures visualize the optimization results shown in Table 7; they show that the predicted response was maximized at the point of (X_{1}=–1, X_{2}=0.08).