Menurut Lawson (2014):
Response surface methods generally refer to a complete package of statistical design and analysis tools that are used for the following three steps. 1. Design and collection of data to fit an equation to approximate the rela- tionship between the factors and response. 2. Regression analysis to fit a model to describe the data. 3. Examination of the fitted relationship through graphical and numerical techniques.
Berikut ini adalah ilustrasi yang dirujuk dari Lawson (2014).
Central Composite Design (CCD)
#install.packages("rsm")
library(rsm)
ccd.pick(k=3) #find a good central-composite design## n.c n0.c blks.c n.s n0.s bbr.c wbr.s bbr.s N alpha.rot alpha.orth
## 1 8 9 1 6 6 1 1 1 29 1.681793 1.680336
## 2 8 2 1 6 1 1 1 1 17 1.681793 1.673320
## 3 8 6 1 6 4 1 1 1 24 1.681793 1.690309
## 4 8 5 1 6 3 1 1 1 22 1.681793 1.664101
## 5 8 10 1 6 7 1 1 1 31 1.681793 1.699673
## 6 8 8 1 6 5 1 1 1 27 1.681793 1.658312
## 7 8 3 1 6 2 1 1 1 19 1.681793 1.705606
## 8 8 7 1 6 5 1 1 1 26 1.681793 1.712698
## 9 8 4 1 6 2 1 1 1 20 1.681793 1.632993
## 10 8 4 1 6 3 1 1 1 21 1.681793 1.732051Untuk menyusun salah satu dari rancangan yang dihasilkan oleh fungsi ccd.pick() di atas, kita dapat menggunakan fungsi ccd() sebagai berikut.
library(rsm)
ccd.up<-ccd(y~x1+x2+x3, n0=c(4,2), alpha="rotatable", randomize=FALSE)
head(ccd.up)## run.order std.order x1.as.is x2.as.is x3.as.is y Block
## 1 1 1 -1 -1 -1 NA 1
## 2 2 2 1 -1 -1 NA 1
## 3 3 3 -1 1 -1 NA 1
## 4 4 4 1 1 -1 NA 1
## 5 5 5 -1 -1 1 NA 1
## 6 6 6 1 -1 1 NA 1
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ x1.as.is
## x2 ~ x2.as.is
## x3 ~ x3.as.isKita dapat pula menambahkan informasi taraf faktor yang digunakan seperti di bawah ini.
ccd.up<-ccd(y~x1+x2+x3,n0=c(4,2),alpha="rotatable", coding=list(x1~(Temp-150)/10,x2~(Press-50)/5, x3~(Rate-4/1)), randomize=FALSE)
head(ccd.up)## run.order std.order Temp Press Rate y Block
## 1 1 1 140 45 3 NA 1
## 2 2 2 160 45 3 NA 1
## 3 3 3 140 55 3 NA 1
## 4 4 4 160 55 3 NA 1
## 5 5 5 140 45 5 NA 1
## 6 6 6 160 45 5 NA 1
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (Temp - 150)/10
## x2 ~ (Press - 50)/5
## x3 ~ (Rate - 4/1)Box-Behnken Design
Treb<-bbd(y~x1+x2+x3, randomize=FALSE, n0=3,
coding=list(x1~(A-6)/2,
x2~(B-15)/5,
x3~(C-2.5)/.5))
head(Treb)## run.order std.order A B C y
## 1 1 1 4 10 2.5 NA
## 2 2 2 8 10 2.5 NA
## 3 3 3 4 20 2.5 NA
## 4 4 4 8 20 2.5 NA
## 5 5 5 4 15 2.0 NA
## 6 6 6 8 15 2.0 NA
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (A - 6)/2
## x2 ~ (B - 15)/5
## x3 ~ (C - 2.5)/0.5Rancangan di atas dapat ditampilkan pula dengan cara berikut.
#install.packages("Vdgraph")
library(Vdgraph)
data(D310)
D310## x1 x2 x3
## 1 0.0000 0.000 1.2906
## 2 0.0000 0.000 -0.1360
## 3 -1.0000 -1.000 0.6386
## 4 1.0000 -1.000 0.6386
## 5 -1.0000 1.000 0.6386
## 6 1.0000 1.000 0.6386
## 7 1.7636 0.000 -0.9273
## 8 -1.7636 0.000 -0.9273
## 9 0.0000 1.736 -0.9273
## 10 0.0000 -1.736 -0.9273des <- transform(D310, Temp=10*x1+150,
Press=5*x2+50,
Rate=x3+4)
des[sample(1:10) ,4:6]## Temp Press Rate
## 8 132.364 50.00 3.0727
## 4 160.000 45.00 4.6386
## 6 160.000 55.00 4.6386
## 10 150.000 41.32 3.0727
## 5 140.000 55.00 4.6386
## 1 150.000 50.00 5.2906
## 7 167.636 50.00 3.0727
## 2 150.000 50.00 3.8640
## 9 150.000 58.68 3.0727
## 3 140.000 45.00 4.6386Non-Standard Response Surface Designs
#install.packages("daewr")
library(daewr)## Registered S3 method overwritten by 'DoE.base':
## method from
## factorize.factor conf.design##
## Attaching package: 'daewr'## The following object is masked _by_ '.GlobalEnv':
##
## Trebdata(qsar)
#install.packages("AlgDesign")
library(AlgDesign)
desgn1<-optFederov(~quad(.),data=qsar,nTrials=15,center=TRUE, criterion="D",nRepeats=40)
desgn2<-optFederov(~quad(.),data=qsar,nTrials=15,center=TRUE,
criterion="I",nRepeats=40)
Compare2FDS(desgn1$design,desgn2$design,"D-optimal", "I-optimal",mod=2)
desgn2$design## Compound HE DMz S0K
## 1 1 -12.221 -0.162 64.138
## 4 4 -14.893 1.035 96.053
## 9 9 -11.813 1.219 77.020
## 12 12 -14.460 2.266 109.535
## 13 13 -8.519 -0.560 71.949
## 14 14 -10.287 -0.675 96.600
## 16 16 -11.167 0.418 104.047
## 19 19 -14.491 -0.561 88.547
## 22 22 -13.121 -1.692 101.978
## 28 28 -12.637 -2.762 112.492
## 29 29 -12.118 -2.994 81.106
## 32 32 -14.804 -3.780 113.856
## 33 33 -9.209 -0.423 74.871
## 34 34 -10.970 -0.302 99.603
## 36 36 -11.868 -1.322 107.010Fitting the Response Surface Model
Linear Model
library(daewr)
data(cement)
head(cement)## Block WatCem BlackL SNF y
## C1.1 1 0.33 0.12 0.08 109.5
## C1.2 1 0.35 0.12 0.08 117.0
## C1.3 1 0.33 0.18 0.08 110.5
## C1.4 1 0.35 0.18 0.08 121.0
## C1.5 1 0.33 0.12 0.12 120.0
## C1.6 1 0.35 0.12 0.12 130.0
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (WatCem - 0.34)/0.01
## x2 ~ (BlackL - 0.15)/0.03
## x3 ~ (SNF - 0.1)/0.02library(rsm)
grout.lin <- rsm(y ~ SO(x1, x2, x3),data = cement, subset = (Block == 1))
anova(grout.lin)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3) 3 465.13 155.042 80.3094 0.002307 **
## TWI(x1, x2, x3) 3 0.25 0.083 0.0432 0.985889
## PQ(x1, x2, x3) 1 37.88 37.879 19.6207 0.021377 *
## Residuals 3 5.79 1.931
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1General Quadratic Model
data(Treb)
head(Treb)## A B C y
## 1 4 10 2.5 33
## 2 8 10 2.5 85
## 3 4 20 2.5 86
## 4 8 20 2.5 113
## 5 4 15 2.0 75
## 6 8 15 2.0 105
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (A - 6)/2
## x2 ~ (B - 15)/5
## x3 ~ (C - 2.5)/0.5treb.quad <- rsm(y ~ SO(x1, x2, x3), data = Treb)
summary(treb.quad)##
## Call:
## rsm(formula = y ~ SO(x1, x2, x3), data = Treb)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 90.00000 1.16905 76.9859 7.006e-09 ***
## x1 19.75000 0.71589 27.5880 1.171e-06 ***
## x2 19.75000 0.71589 27.5880 1.171e-06 ***
## x3 -11.50000 0.71589 -16.0639 1.703e-05 ***
## x1:x2 -6.25000 1.01242 -6.1733 0.0016247 **
## x1:x3 4.75000 1.01242 4.6917 0.0053768 **
## x2:x3 6.75000 1.01242 6.6672 0.0011461 **
## x1^2 -9.37500 1.05376 -8.8967 0.0002986 ***
## x2^2 -1.37500 1.05376 -1.3048 0.2487686
## x3^2 -3.37500 1.05376 -3.2028 0.0239200 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9975, Adjusted R-squared: 0.9929
## F-statistic: 218.9 on 9 and 5 DF, p-value: 5.964e-06
##
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3) 3 7299.0 2433.00 593.4146 8.448e-07
## TWI(x1, x2, x3) 3 428.8 142.92 34.8577 0.0008912
## PQ(x1, x2, x3) 3 351.5 117.16 28.5759 0.0014236
## Residuals 5 20.5 4.10
## Lack of fit 3 14.5 4.83 1.6111 0.4051312
## Pure error 2 6.0 3.00
##
## Stationary point of response surface:
## x1 x2 x3
## 0.9236846 -1.7161183 -2.7698217
##
## Stationary point in original units:
## A B C
## 7.847369 6.419409 1.115089
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 1.280298 -3.551452 -11.853845
##
## $vectors
## [,1] [,2] [,3]
## x1 -0.1236692 0.5238084 0.8428112
## x2 0.8323200 -0.4077092 0.3755217
## x3 0.5403233 0.7479291 -0.3855551Berikut ini adalah contoh apabila menambahkan komponen blok.
data(cement)
grout.quad<-rsm(y ~ Block + SO(x1,x2,x3), data=cement)
summary(grout.quad)##
## Call:
## rsm(formula = y ~ Block + SO(x1, x2, x3), data = cement)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.1628e+02 1.0691e+00 108.7658 2.383e-15 ***
## Block2 4.4393e-01 1.0203e+00 0.4351 0.67375
## x1 5.4068e+00 6.1057e-01 8.8553 9.746e-06 ***
## x2 9.2860e-01 6.1057e-01 1.5209 0.16262
## x3 4.9925e+00 6.1057e-01 8.1767 1.858e-05 ***
## x1:x2 1.2500e-01 7.9775e-01 0.1567 0.87895
## x1:x3 -1.3443e-14 7.9775e-01 0.0000 1.00000
## x2:x3 1.2500e-01 7.9775e-01 0.1567 0.87895
## x1^2 1.4135e+00 5.9582e-01 2.3723 0.04175 *
## x2^2 1.3251e+00 5.9582e-01 2.2240 0.05322 .
## x3^2 1.5019e+00 5.9582e-01 2.5207 0.03273 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9473, Adjusted R-squared: 0.8887
## F-statistic: 16.17 on 10 and 9 DF, p-value: 0.0001414
##
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## Block 1 0.00 0.003 0.0006 0.98068
## FO(x1, x2, x3) 3 751.41 250.471 49.1962 6.607e-06
## TWI(x1, x2, x3) 3 0.25 0.083 0.0164 0.99693
## PQ(x1, x2, x3) 3 71.45 23.817 4.6779 0.03106
## Residuals 9 45.82 5.091
## Lack of fit 5 42.49 8.498 10.1972 0.02149
## Pure error 4 3.33 0.833
##
## Stationary point of response surface:
## x1 x2 x3
## -1.9045158 -0.1825251 -1.6544845
##
## Stationary point in original units:
## WatCem BlackL SNF
## 0.32095484 0.14452425 0.06691031
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 1.525478 1.436349 1.278634
##
## $vectors
## [,1] [,2] [,3]
## x1 0.1934409 0.8924556 0.4075580
## x2 0.3466186 0.3264506 -0.8793666
## x3 0.9178432 -0.3113726 0.2461928Nonlinear Mechanistic Model
data(Tet)
mod.nln1 <- nls( Conc ~ gamma0 * (exp( -k1 * (Time - t0)) - exp( -k2 * (Time - t0))), data = Tet, start = list( gamma0=10, k1 = .12, k2 = .5, t0 = .5))
summary(mod.nln1)##
## Formula: Conc ~ gamma0 * (exp(-k1 * (Time - t0)) - exp(-k2 * (Time - t0)))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## gamma0 2.64964 0.36446 7.270 0.000770 ***
## k1 0.14880 0.01441 10.327 0.000146 ***
## k2 0.71575 0.12605 5.678 0.002359 **
## t0 0.41224 0.09495 4.342 0.007416 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04482 on 5 degrees of freedom
##
## Number of iterations to convergence: 11
## Achieved convergence tolerance: 3.632e-06Optimum Operating Conditions
Contour Plots
data(Treb)
library(rsm)
treb.quad <- rsm(y ~ SO(x1, x2, x3), data = Treb)
par (mfrow=c(2,2))
contour(treb.quad, ~ x1+x2+x3 )
par (mfrow=c(2,2))
persp(treb.quad, ~ x1+x2+x3, zlab="Distance", contours=list(z="bottom") )
par (mfrow=c(1,2))
contour(treb.quad, x1~x3, at=list(x2=1))
persp(treb.quad, x1~x3, at=list(x2=1),zlab="Distance",
contours=list(z="bottom"))
Reference
Lawson, J. (2014). Design and Analysis of Experiments with R (Vol. 115). CRC press.