STAT 345/834: Design and Analysis of Experiments
Unit 8: Introduction Response Surface Design
Longhai Li
March 2018
1 Draw Some Response Surfaces
1.1 An Additive Model
add.model = function(x1,x2)
{
x1+x2
}
# response surface plot
x2<- x1 <- seq(0,1,by=0.01)
z <- outer(x1,x2, add.model)
persp (x1, x2, z, theta = 30)
# contour plot
contour(x1,x2,z)
1.2 An Interaction Model
inter.model = function(x1,x2)
{
x1+x2+5*x1*x2
}
# response surface plot
x2<- x1 <- seq(0,1,by=0.01)
z <- outer(x1,x2, inter.model)
persp(x1, x2, z, theta=45)
# contour plot
contour(x1,x2,z)
1.3 A Circular Quadratic Model
quad.model = function(x1,x2)
{
x1+x2-10*(x1-0.5)^2-10*(x2-0.5)^2
}
# response surface plot
x2<- x1 <- seq(0,1,by=0.01)
z <- outer(x1,x2, quad.model)
persp(x1, x2, z, theta=45)
# contour plot
contour(x1,x2,z)
1.4 An Elliptical Quadratic Model
quad.model = function(x1,x2)
{
-x1-x2-10*(x1-0.5)^2-40*(x2-0.5)^2
}
# response surface plot
x2<- x1 <- seq(0,1,by=0.01)
z <- outer(x1,x2, quad.model)
persp(x1, x2, z, theta=45)
# contour plot
contour(x1,x2,z)
1.5 A Complex Model
quad.model = function(x1,x2)
{
x1+x2+5*x1*x2-10*(x1-.5)^2-20*(x2-.5)^2
}
# response surface plot
x2<- x1 <- seq(0,1,by=0.01)
z <- outer(x1,x2, quad.model)
persp(x1, x2, z, theta = 45)
# contour plot
contour(x1,x2,z)
2 Design and Analysis of Second-order response surface around optimum
2.1 Read in the data
yield<- read.table("data/yield.txt", header=T)
yield## y
## 1 76.5
## 2 77.0
## 3 78.0
## 4 79.5
## 5 79.9
## 6 80.3
## 7 80.0
## 8 79.7
## 9 79.8
## 10 78.4
## 11 75.6
## 12 78.5
## 13 77.0attach(yield)2.2 Without Central Composite Design
2.2.1 Fit a 2nd-order model
# define x1 and x2
x1<-c(-1,-1,1,1,rep(0,5))
x2<-c(-1,1,-1,1,rep(0,5))
data.frame (x1,x2, y=y[1:9]) # using only the first 9 observations## x1 x2 y
## 1 -1 -1 76.5
## 2 -1 1 77.0
## 3 1 -1 78.0
## 4 1 1 79.5
## 5 0 0 79.9
## 6 0 0 80.3
## 7 0 0 80.0
## 8 0 0 79.7
## 9 0 0 79.8g<-lm(y[1:9]~x1+x2+I(x1^2)+x1*x2,data=yield)
summary(g)##
## Call:
## lm(formula = y[1:9] ~ x1 + x2 + I(x1^2) + x1 * x2, data = yield)
##
## Residuals:
## 1 2 3 4 5 6
## 0.000e+00 -6.939e-18 -7.254e-18 4.889e-18 -4.000e-02 3.600e-01
## 7 8 9
## 6.000e-02 -2.400e-01 -1.400e-01
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 79.9400 0.1030 776.446 1.65e-11 ***
## x1 1.0000 0.1151 8.687 0.000966 ***
## x2 0.5000 0.1151 4.344 0.012217 *
## I(x1^2) -2.1900 0.1544 -14.181 0.000144 ***
## x1:x2 0.2500 0.1151 2.172 0.095611 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2302 on 4 degrees of freedom
## Multiple R-squared: 0.9868, Adjusted R-squared: 0.9737
## F-statistic: 75.04 on 4 and 4 DF, p-value: 0.00051432.2.2 ANOVA
anova (g)## Analysis of Variance Table
##
## Response: y[1:9]
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 4.000 4.000 75.472 0.0009664 ***
## x2 1 1.000 1.000 18.868 0.0122172 *
## I(x1^2) 1 10.658 10.658 201.094 0.0001436 ***
## x1:x2 1 0.250 0.250 4.717 0.0956108 .
## Residuals 4 0.212 0.053
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.3 Analysis of Central Composite Design
In CCD,we take 4 more observations so that we can estimate pure quadratic effect of \(x_2\).
2.3.1 Fit a 2nd order model with pure quadratic terms
x1<-c(-1,-1,1,1,rep(0,5),1.414,-1.414,0,0)
x2<-c(-1,1,-1,1,rep(0,5),0,0,1.414,-1.414)
data.frame (x1,x2, y)## x1 x2 y
## 1 -1.000 -1.000 76.5
## 2 -1.000 1.000 77.0
## 3 1.000 -1.000 78.0
## 4 1.000 1.000 79.5
## 5 0.000 0.000 79.9
## 6 0.000 0.000 80.3
## 7 0.000 0.000 80.0
## 8 0.000 0.000 79.7
## 9 0.000 0.000 79.8
## 10 1.414 0.000 78.4
## 11 -1.414 0.000 75.6
## 12 0.000 1.414 78.5
## 13 0.000 -1.414 77.0h<-lm(y~x1+x2+I(x1^2)+I(x2^2)+x1*x2,data=yield)
summary(h)##
## Call:
## lm(formula = y ~ x1 + x2 + I(x1^2) + I(x2^2) + x1 * x2, data = yield)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.23995 -0.18089 -0.03995 0.17758 0.36005
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 79.93995 0.11909 671.264 < 2e-16 ***
## x1 0.99505 0.09415 10.568 1.48e-05 ***
## x2 0.51520 0.09415 5.472 0.000934 ***
## I(x1^2) -1.37645 0.10098 -13.630 2.69e-06 ***
## I(x2^2) -1.00134 0.10098 -9.916 2.26e-05 ***
## x1:x2 0.25000 0.13315 1.878 0.102519
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2663 on 7 degrees of freedom
## Multiple R-squared: 0.9827, Adjusted R-squared: 0.9704
## F-statistic: 79.67 on 5 and 7 DF, p-value: 5.147e-062.3.2 ANOVA
anova (h)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 7.9198 7.9198 111.6873 1.484e-05 ***
## x2 1 2.1232 2.1232 29.9413 0.000934 ***
## I(x1^2) 1 10.9816 10.9816 154.8663 4.979e-06 ***
## I(x2^2) 1 6.9721 6.9721 98.3225 2.262e-05 ***
## x1:x2 1 0.2500 0.2500 3.5256 0.102519
## Residuals 7 0.4964 0.0709
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.3.3 Analysis after removing interaction term (non-significant)
h2<-lm(y~x1+x2+I(x1^2)+I(x2^2),data=yield)
summary(h2)##
## Call:
## lm(formula = y ~ x1 + x2 + I(x1^2) + I(x2^2), data = yield)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.23995 -0.18089 -0.08232 0.06005 0.44808
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 79.9400 0.1366 585.215 < 2e-16 ***
## x1 0.9951 0.1080 9.213 1.56e-05 ***
## x2 0.5152 0.1080 4.770 0.00141 **
## I(x1^2) -1.3764 0.1158 -11.883 2.31e-06 ***
## I(x2^2) -1.0013 0.1158 -8.645 2.49e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3054 on 8 degrees of freedom
## Multiple R-squared: 0.974, Adjusted R-squared: 0.961
## F-statistic: 75.02 on 4 and 8 DF, p-value: 2.226e-06anova(h2)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 7.9198 7.9198 84.888 1.559e-05 ***
## x2 1 2.1232 2.1232 22.757 0.001408 **
## I(x1^2) 1 10.9816 10.9816 117.707 4.604e-06 ***
## I(x2^2) 1 6.9721 6.9721 74.730 2.489e-05 ***
## Residuals 8 0.7464 0.0933
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1detach ("yield")