Abstract
Until now, we have considered how a dependent variable, yield, or response depends on specific levels of independent variables or factors. The factors could be categorical or numerical; however, we did note that they often differ in how the sum of squares for the factor is more usefully partitioned into orthogonal components. For example, a numerical factor might be broken down into orthogonal polynomials (introduced in Chap. 12). For categorical factors, methods introduced in Chap. 5 are typically employed. In the past two chapters, we have considered linear relationships and fitting optimal straight lines to the data, usually for situations in which the data values are not derived from designed experiments. Now, we consider experimental design techniques that find the optimal combination of factor levels for situations in which the feasible levels of each factor are continuous. (Throughout the text, the dependent variable, Y, has been assumed to be continuous.) The techniques are called response-surface methods or response-surface methodology (RSM).
This is a preview of subscription content, log in via an institution.
Notes
- 1.
G. E. P. Box and K. B. Wilson (1951), “On the Experimental Attainment of Optimum Conditions.” Journal of the Royal Statistical Society, Series B, vol. 13, pp. 1–45.
- 2.
By “flat,” we mean from a curvature-of-the-earth perspective; we are not referring to the issue of hills and valleys.
- 3.
In general, the only time that stage two requires more than one experiment is the unfortunate, but not unlikely, situation of a “saddle point,” or some other undesirable happenstance. This issue is discussed later in the chapter.
- 4.
This example, and a broad outline of some of the features used for the illustration of the RSM process, were suggested by the discussion in C. R. Hicks, Fundamental Concepts in the Design of Experiments, 3rd edition, New York, Holt, Rinehart and Winston, 1982.
- 5.
In the real world, issues involving the packaging and other aspects of banana preservation are quite complex. There are many considerations and a large number of factors involved. However, as we noted, our example is greatly simplified, but captures the features we wish to illustrate. The example is loosely based on one author’s experience designing a complex experiment addressing some of these issues for a well-known harvester and shipper of bananas. Incidentally, in general, brown spots on bananas are not unhealthy or inedible – only unsightly!
- 6.
We obtain the transformation for U as follows: let U = a + bX 1. When X 1 = 1/4, U = −1; this gives −1 = a + b(1/4) = a + b/4; when X 1 = 1/2 , U = +1; this gives 1 = a + b(1/2) = a + b/2. We have two equations with two unknowns, a and b. We solve these and find a = −3 and b = 8. The second equation is found in a similar fashion.
- 7.
Latitude refers to the angle that measures north-south position from the equatorial plane, whereas longitude is the angle that measures west-east position from the Greenwich Meridian.
- 8.
We repeat that this will be the last experiment, unless the situation is an “unlucky” one, in which there is either a saddle point, or local maxima in addition to the global maximum.
- 9.
G. E. P. Box and J. S. Hunter (1957), “Multi-factor Experimental Designs for Exploring Response Surfaces.” Annals of Mathematical Statistics, vol. 28, n. 1, pp. 195–241.
- 10.
In theory, the point could be a minimum point, but, unless the surface is very irregular, this is very unlikely after being directionally guided by a series of steepest-ascent analyses. A saddle point represents a stationary situation where there is neither a maximum nor a minimum – this term is used due to the resemblance of the surface to a saddle. That is, if you are on a horse, sitting on a saddle, facing north-south, then from a north-south direction, you are sitting at a minimum point, while from an east-west direction, you are sitting at a maximum point.
- 11.
Note that the fractional-factorial design described is of resolution five.
- 12.
The authors don’t know where this world is; if we did, we’d buy property there!
- 13.
Remember: We are describing what the NASA team did, not necessarily what we would have done; then again, we weren’t there, and they made the judgment call. This thought applies to many of the subsequent steps the NASA researchers chose to follow.
- 14.
Actually, it was the taxpayers’ money.
- 15.
Of course, search techniques were available back in the late 1960s on mainframe computers. However, they were not readily accessible nor user-friendly. Today, one might question whether the calculus step shouldn’t always be replaced by Solver or its equivalent.
- 16.
G. Derringer and R. Suich (1980), “Simultaneous Optimization of Several Response Variables.” Journal of Quality Technology, vol. 12, n. 4, pp. 214–219.
Author information
Authors and Affiliations
1 Electronic Supplementary Material
Example 16.5
– NASA study (XLS 69 kb)
Example 16.6
– Extraction yield (XLS 59 kb)
Appendix
Appendix
Example 16.9 NASA Example in R
In these final demonstrations, we will use rsm and DoE.wrapper packages.
# Option 1: using rsm package
First, we create the design matrix and an object with the responses, which are later combined.
> NASA <- ccd(4, n0=2, randomize=FALSE, alpha="rotatable")
# n0 is the number of center points per block.
> y <- c(2.180699, 2.278868, 2.198657,…, 2.371160) > NASA <- cbind(NASA, y) > NASA
run.order | std.order | x1 | x2 | x3 | x4 | Block | y | |
1 | 1 | 1 | -1 | -1 | -1 | -1 | 1 | 2.180699 |
2 | 2 | 2 | 1 | -1 | -1 | -1 | 1 | 2.278868 |
3 | 3 | 3 | -1 | 1 | -1 | -1 | 1 | 2.198657 |
⋮ | ⋮ | |||||||
2 8 | 10 | 10 | 0 | 0 | 0 | 0 | 2 | 2.371160 |
# Option 2: using DoE.wrapper package
> NASA <- ccd.design(4, ncenter=2, randomize=FALSE, alpha= +"rotatable") full factorial design needed creating full factorial with 16 runs ... > y <- c(2.180699, 2.278868, 2.198657,..., 2.371160) > NASA <- add.response(NASA, y) > NASA
Block.ccd | X1 | X2 | X3 | X4 | y | |
C1.1 | 1 | -1 | -1 | -1 | -1 | 2.180699 |
C1.2 | 1 | 1 | -1 | -1 | -1 | 2.278868 |
C1.3 | 1 | -1 | 1 | -1 | -1 | 2.198657 |
⋮ | ⋮ | |||||
C2.10 | 2 | 0 | 0 | 0 | 0 | 2.371160 |
class = design, type = ccd
We can analyze either of the options using the rsm() function described below, where SO stands for “second-order,” FO for “first-order,” PQ for “pure quadratic,” and TWI for “two-way interaction ”:
> NASA_rsm <- rsm(y~SO(x1, x2, x3, x4), data=NASA) > summary(NASA_rsm) Call: rsm(formula = y ~ SO(x1, x2, x3, x4), data = NASA)
Estimate | Std. Error | t value | Pr(>|t|) | ||
(Intercept) | 2.354352 | 0.047993 | 49.0558 | 3.838e-16 | *** |
x1 | -0.011898 | 0.019593 | -0.6073 | 0.5541321 | |
x2 | 0.013357 | 0.019593 | 0.6817 | 0.5073746 | |
x3 | 0.029071 | 0.019593 | 1.4837 | 0.1617127 | |
x4 | -0.036931 | 0.019593 | -1.8849 | 0.0819968 | . |
x1:x2 | -0.003351 | 0.023997 | -0.1396 | 0.8910812 | |
x1:x3 | 0.024473 | 0.023997 | 1.0199 | 0.3263951 | |
x1:x4 | -0.053768 | 0.023997 | -2.2406 | 0.0431455 | * |
x2:x3 | -0.010251 | 0.023997 | -0.4272 | 0.6762391 | |
x2:x4 | -0.021033 | 0.023997 | -0.8765 | 0.3966624 | |
x3:x4 | 0.077319 | 0.023997 | 3.2221 | 0.0066779 | ** |
x1^2 | -0.084317 | 0.019593 | -4.3034 | 0.0008579 | *** |
x2^2 | -0.026245 | 0.019593 | -1.3395 | 0.2033652 | |
x3^2 | -0.025982 | 0.019593 | -1.3261 | 0.2076414 | |
x4^2 | -0.059971 | 0.019593 | -3.0608 | 0.0091087 | ** |
--- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Multiple R-squared: 0.7825, Adjusted R-squared: 0.5482 F-statistic: 3.34 on 14 and 13 DF, p-value: 0.0182 Analysis of Variance Table Response: y
Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
FO(x1, x2, x3, x4) | 4 | 0.060696 | 0.015174 | 1.6469 | 0.22181 |
TWI(x1, x2, x3, x4) | 6 | 0.160430 | 0.026738 | 2.9021 | 0.05070 |
PQ(x1, x2, x3, x4) | 4 | 0.209692 | 0.052423 | 5.6899 | 0.00715 |
Residuals | 13 | 0.119775 | 0.009213 | ||
Lack of fit | 10 | 0.104698 | 0.010470 | 2.0833 | 0.29661 |
Pure error | 3 | 0.015076 | 0.005025 |
Stationary point of response surface:
x1 | x2 | x3 | x4 |
0.006267162 | 0.176237359 | 0.473383841 | -0.036460362 |
Eigenanalysis: $values
[1] | 0.003014912 | -0.029150191 | -0.055951382 | -0.114428173 |
$vectors
[,1] | [,2] | [,3] | [,4] | |
x1 | 0.05527429 | 0.002767506 | 0.71328348 | 0.69868718 |
x2 | 0.32843898 | -0.938320173 | -0.08597682 | 0.06550631 |
x3 | -0.76662083 | -0.331545323 | 0.41456151 | -0.36126024 |
x4 | -0.54896730 | -0.098108573 | -0.55854581 | 0.61403272 |
Below, we show the command to obtain the coefficients of our model, and using the predict() function we can obtain the solution using the values found under “Stationary point of response surface” on the output.
> NASA_rsm Call: rsm(formula = y ~ SO(x1, x2, x3, x4), data = NASA) Coefficients:
(Intercept) | FO(x1, x2, x3, x4)x1 | FO(x1, x2, x3, x4)x2 |
2.354352 | -0.011898 | 0.013357 |
FO(x1, x2, x3, x4)x3 | FO(x1, x2, x3, x4)x4 | TWI(x1, x2, x3, x4)x1:x2 |
0.029071 | -0.036931 | -0.003351 |
TWI(x1, x2, x3, x4)x1:x3 | TWI(x1, x2, x3, x4)x1:x4 | TWI(x1, x2, x3, x4)x2:x3 |
0.024474 | -0.053768 | -0.010251 |
TWI(x1, x2, x3, x4)x2:x4 | TWI(x1, x2, x3, x4)x3:x4 | PQ(x1, x2, x3, x4)x1^2 |
-0.021033 | 0.077319 | -0.084317 |
PQ(x1, x2, x3, x4)x2^2 | PQ(x1, x2, x3, x4)x3^2 | PQ(x1, x2, x3, x4)x4^2 |
-0.026245 | -0.025982 | -0.059971 |
> predict(NASA_rsm, newdata=data.frame(X1=0.006267162, +X2=0.176237359, X3=0.473383841, X4=-0.036460362))
1 |
2.363046 |
Next, using the following commands, we can generate contour plots for the pairs of X’s, shown in Fig. 16.14:
> par(mfrow=c(2,3)) > contour(NASA_rsm, ~x1+x2+x3+x4, at=summary(NASA_rsm)+$canonical$xs)
Example 16.10 Extraction Yield of Bioactive Compounds in R
Here, we will also demonstrate how to set up and analyze a Box-Behnken design using rsm and DoE.wrapper packages; this is similar to what we did in the previous example.
# Option 1: using rsm package
First, we create the design matrix and an object with the responses, which are later combined.
> extract <- bbd(4, n0=3, block=FALSE, randomize=FALSE) > y <- c(7.3, 17.6, 13, …, 22.2) > extract <- cbind(extract, y) > extract
run.order | std.order | x1 | x2 | x3 | x4 | y | |
1 | 1 | 1 | -1 | -1 | 0 | 0 | 7.3 |
2 | 2 | 2 | 1 | -1 | 0 | 0 | 17.6 |
3 | 3 | 3 | -1 | 1 | 0 | 0 | 13.0 |
⋮ | ⋮ | ||||||
27 | 27 | 27 | 0 | 0 | 0 | 0 | 22.2 |
# Option 2: using DoE.wrapper package
> extract <- bbd.design(4, ncenter=2, randomize=FALSE) > y <- c(7.3, 17.6, 13, …, 22.2) > extract <- add.response(extract, y) > extract
A | B | C | D | y | |
1 | -1 | -1 | 0 | 0 | 7.3 |
2 | 1 | -1 | 0 | 0 | 17.6 |
3 | -1 | 1 | 0 | 0 | 13.0 |
⋮ | ⋮ | ||||
27 | 0 | 0 | 0 | 0 | 22.2 |
class = design, type= bbd
As with CCD, we can analyze either of the options using the rsm() function described below:
> extract_rsm <- rsm(y~SO(A, B, C, D), data=extract) > summary(extract_rsm) Call: rsm(formula = y ~ SO(A, B, C, D), data = extract)
Estimate | Std. Error | t value | Pr(>|t|) | ||
(Intercept) | 22.53333 | 0.36723 | 61.3595 | 2.323e-16 | *** |
A | 1.90833 | 0.18362 | 10.3930 | 2.358e-07 | *** |
B | -1.10833 | 0.18362 | -6.0361 | 5.882e-05 | *** |
C | 0.22500 | 0.18362 | 1.2254 | 0.24394 | |
D | 1.17500 | 0.18362 | 6.3992 | 3.407e-05 | *** |
A:B | -3.65000 | 0.31803 | -11.4767 | 7.934e-08 | *** |
A:C | 0.17500 | 0.31803 | 0.5503 | 0.59224 | |
A:D | 0.15000 | 0.31803 | 0.4716 | 0.64564 | |
B:C | -0.22500 | 0.31803 | -0.7075 | 0.49279 | |
B:D | -0.80000 | 0.31803 | -2.5155 | 0.02713 | * |
C:D | -0.17500 | 0.31803 | -0.5503 | 0.59224 | |
A^2 | -6.01667 | 0.27543 | -21.8450 | 4.967e-11 | *** |
B^2 | -4.74167 | 0.27543 | -17.2158 | 7.966e-10 | *** |
C^2 | -2.54167 | 0.27543 | -9.2281 | 8.470e-07 | *** |
D^2 | -3.69167 | 0.27543 | -13.4035 | 1.397e-08 | *** |
--- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Multiple R-squared: 0.9873, Adjusted R-squared: 0.9725 F-statistic: 66.6 on 14 and 12 DF, p-value: 3.642e-09 Analysis of Variance Table Response: y
Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
FO(A, B, C, D) | 4 | 75.617 | 18.904 | 46.7250 | 3.201e-07 |
TWI(A, B, C, D) | 6 | 56.387 | 9.398 | 23.2286 | 6.038e-06 |
PQ(A, B, C, D) | 4 | 245.222 | 61.305 | 151.5272 | 3.686e-10 |
Residuals | 12 | 4.855 | 0.405 | ||
Lack of fit | 10 | 4.668 | 0.467 | 5.0018 | 0.178 |
Pure error | 2 | 0.187 | 0.093 |
Stationary point of response surface:
A | B | C | D |
0.22908869 | -0.22209945 | 0.05555756 | 0.18654451 |
Eigenanalysis: $values
[1] | -2.518852 | -3.182962 | -3.969465 | -7.320388 |
$vectors
[,1] | [,2] | [,3] | [,4] | |
A | -0.08296588 | 0.43894755 | 0.3759016 | 0.8118741835 |
B | 0.11279773 | -0.66136843 | -0.4596205 | 0.5819084989 |
C | -0.98969169 | -0.09392243 | -0.1081149 | -0.0002993512 |
D | 0.03006143 | 0.60091217 | -0.7973444 | 0.0473573592 |
> extract_rsm Call: rsm(formula = y ~ SO(A, B, C, D), data = extract) Coefficients:
(Intercept) | FO(A, B, C, D)A | FO(A, B, C, D)B |
22.533 | 1.908 | -1.108 |
FO(A, B, C, D)C | FO(A, B, C, D)D | TWI(A, B, C, D)A:B |
0.225 | 1.175 | -3.650 |
TWI(A, B, C, D)A:C | TWI(A, B, C, D)A:D | TWI(A, B, C, D)B:C |
0.175 | 0.150 | -0.225 |
TWI(A, B, C, D)B:D | TWI(A, B, C, D)C:D | PQ(A, B, C, D)A^2 |
-0.800 | -0.175 | -6.017 |
PQ(A, B, C, D)B^2 | PQ(A, B, C, D)C^2 | PQ(A, B, C, D)D^2 |
-4.742 | -2.542 | -3.692 |
Using the predict() function, we can obtain the solution using the values found under “Stationary point of response surface” in the output.
> predict(extract_rsm, newdata=data.frame(A=0.22908869, +B=-0.22209945, C=0.05555756, D=0.18654451))
1 |
22.99085 |
Six contour plots are shown in Fig. 16.15.
> par(mfrow=c(2,3)) > contour(extract_rsm, ~A+B+C+D, at=summary(extract_rsm)+$canonical$xs)
Rights and permissions
Copyright information
© 2018 Springer International Publishing AG
About this chapter
Cite this chapter
Berger, P.D., Maurer, R.E., Celli, G.B. (2018). Introduction to Response-Surface Methodology. In: Experimental Design. Springer, Cham. https://doi.org/10.1007/978-3-319-64583-4_16
Download citation
DOI: https://doi.org/10.1007/978-3-319-64583-4_16
Published:
Publisher Name: Springer, Cham
Print ISBN: 978-3-319-64582-7
Online ISBN: 978-3-319-64583-4
eBook Packages: Mathematics and StatisticsMathematics and Statistics (R0)