Experimento RSM design CCD - Rotacional: Livro Montgomery (DOE): Exemplo 11.2

library(rsm)
library(lmtest)
library(car)
library(nortest)
library(FrF2)
library(rio)
DOE<- ccd(basis = 2,
          alpha = "rotatable",
          n0=c(5,0),
          oneblock = T,
          randomize = F,
          coding = list(x1~(a-85)/5,
                        x2~(b-175)/5))
DOE
   run.order std.order        a        b
1          1         1 80.00000 170.0000
2          2         2 90.00000 170.0000
3          3         3 80.00000 180.0000
4          4         4 90.00000 180.0000
5          5         5 85.00000 175.0000
6          6         6 85.00000 175.0000
7          7         7 85.00000 175.0000
8          8         8 85.00000 175.0000
9          9         9 85.00000 175.0000
10         1         1 77.92893 175.0000
11         2         2 92.07107 175.0000
12         3         3 85.00000 167.9289
13         4         4 85.00000 182.0711

Data are stored in coded form using these coding formulas ...
x1 ~ (a - 85)/5
x2 ~ (b - 175)/5
DOE$y1<-c(76.5,77,78,79.5,79.9,80.3,80,79.7,79.8,
          75.6,78.4,77,78.5)
DOE$y2<-c(62,60,66,59,72,69,68,70,71,71,68,57,58)
DOE$y3<-c(2940,3470,3680,3890,3480,3200,3410,3290,3500,
          3020,3360,3150,3630)
DOE
   run.order std.order        a        b   y1 y2   y3
1          1         1 80.00000 170.0000 76.5 62 2940
2          2         2 90.00000 170.0000 77.0 60 3470
3          3         3 80.00000 180.0000 78.0 66 3680
4          4         4 90.00000 180.0000 79.5 59 3890
5          5         5 85.00000 175.0000 79.9 72 3480
6          6         6 85.00000 175.0000 80.3 69 3200
7          7         7 85.00000 175.0000 80.0 68 3410
8          8         8 85.00000 175.0000 79.7 70 3290
9          9         9 85.00000 175.0000 79.8 71 3500
10         1         1 77.92893 175.0000 75.6 71 3020
11         2         2 92.07107 175.0000 78.4 68 3360
12         3         3 85.00000 167.9289 77.0 57 3150
13         4         4 85.00000 182.0711 78.5 58 3630

Data are stored in coded form using these coding formulas ...
x1 ~ (a - 85)/5
x2 ~ (b - 175)/5
funcao<- rsm(y1~SO(x1,x2),data = DOE)
summary(funcao)

Call:
rsm(formula = y1 ~ SO(x1, x2), data = DOE)

            Estimate Std. Error  t value  Pr(>|t|)    
(Intercept) 79.94000    0.20104 397.6418 < 2.2e-16 ***
x1           0.74497    0.15893   4.6874  0.002241 ** 
x2           0.76517    0.15893   4.8144  0.001934 ** 
x1:x2        0.25000    0.22476   1.1123  0.302754    
x1^2        -1.37625    0.17044  -8.0749 8.588e-05 ***
x2^2        -1.00125    0.17044  -5.8746  0.000615 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple R-squared:  0.9508,    Adjusted R-squared:  0.9156 
F-statistic: 27.05 on 5 and 7 DF,  p-value: 0.0001935

Analysis of Variance Table

Response: y1
            Df  Sum Sq Mean Sq F value    Pr(>F)
FO(x1, x2)   2  9.1237  4.5619 22.5750 0.0008860
TWI(x1, x2)  1  0.2500  0.2500  1.2372 0.3027535
PQ(x1, x2)   2 17.9548  8.9774 44.4260 0.0001053
Residuals    7  1.4145  0.2021                  
Lack of fit  3  1.2025  0.4008  7.5631 0.0399491
Pure error   4  0.2120  0.0530                  

Stationary point of response surface:
       x1        x2 
0.3088613 0.4206644 

Stationary point in original units:
        a         b 
 86.54431 177.10332 

Eigenanalysis:
eigen() decomposition
$values
[1] -0.963403 -1.414097

$vectors
         [,1]       [,2]
x1 -0.2897841 -0.9570920
x2 -0.9570920  0.2897841
funcao<-rsm(y1~FO(x1,x2)+PQ(x1,x2),data = DOE)
summary(funcao)

Call:
rsm(formula = y1 ~ FO(x1, x2) + PQ(x1, x2), data = DOE)

            Estimate Std. Error  t value  Pr(>|t|)    
(Intercept) 79.94000    0.20399 391.8756 < 2.2e-16 ***
x1           0.74497    0.16127   4.6194 0.0017115 ** 
x2           0.76517    0.16127   4.7446 0.0014553 ** 
x1^2        -1.37625    0.17294  -7.9578 4.536e-05 ***
x2^2        -1.00125    0.17294  -5.7895 0.0004101 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple R-squared:  0.9421,    Adjusted R-squared:  0.9131 
F-statistic: 32.54 on 4 and 8 DF,  p-value: 5.363e-05

Analysis of Variance Table

Response: y1
            Df  Sum Sq Mean Sq F value    Pr(>F)
FO(x1, x2)   2  9.1237  4.5619 21.9250 0.0005667
PQ(x1, x2)   2 17.9548  8.9774 43.1469 5.181e-05
Residuals    8  1.6645  0.2081                  
Lack of fit  4  1.4525  0.3631  6.8516 0.0445322
Pure error   4  0.2120  0.0530                  

Stationary point of response surface:
       x1        x2 
0.2706539 0.3821049 

Stationary point in original units:
        a         b 
 86.35327 176.91052 

Eigenanalysis:
eigen() decomposition
$values
[1] -1.00125 -1.37625

$vectors
   [,1] [,2]
x1    0   -1
x2   -1    0
par(mfrow=c(1,2))
curve(dnorm(x,mean(funcao$residuals),sd(funcao$residuals)),
from = min(funcao$residuals),to=max(funcao$residuals),
col="blue",lwd=2,ylab = "Densidade",main="Curva Normal")
qqnorm(funcao$residuals,lwd=2,col="turquoise")
qqline(funcao$residuals,col="red",lwd=2)
shapiro.test(funcao$residuals)

    Shapiro-Wilk normality test

data:  funcao$residuals
W = 0.89607, p-value = 0.1181
ad.test(funcao$residuals)

    Anderson-Darling normality test

data:  funcao$residuals
A = 0.44205, p-value = 0.2423
par(mfrow=c(1,1))

residualPlot(funcao)

bptest(funcao)

    studentized Breusch-Pagan test

data:  funcao
BP = 6.3089, df = 4, p-value = 0.1772
dwt(funcao)
 lag Autocorrelation D-W Statistic p-value
   1      -0.2064809      2.130376   0.944
 Alternative hypothesis: rho != 0
dwtest(funcao,alternative = c("two.sided"))

    Durbin-Watson test

data:  funcao
DW = 2.1304, p-value = 0.9718
alternative hypothesis: true autocorrelation is not 0
bgtest(funcao,order = 1)

    Breusch-Godfrey test for serial correlation of order up
    to 1

data:  funcao
LM test = 0.9643, df = 1, p-value = 0.3261
persp(funcao,~x1+x2,zlab = "y",col=rainbow(510),
      contours = "color")

contour(funcao,~x1+x2,image = T)

contour(funcao,~x1+x2,col=topo.colors(12))

funcao2<-rsm(y2~SO(x1,x2),data = DOE)
summary(funcao2)

Call:
rsm(formula = y2 ~ SO(x1, x2), data = DOE)

            Estimate Std. Error  t value  Pr(>|t|)    
(Intercept) 70.00000    0.74197  94.3429 3.960e-12 ***
x1          -1.65533    0.58658  -2.8220    0.0257 *  
x2           0.55178    0.58658   0.9407    0.3782    
x1:x2       -1.25000    0.82955  -1.5068    0.1756    
x1^2        -0.68750    0.62904  -1.0929    0.3106    
x2^2        -6.68750    0.62904 -10.6313 1.427e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple R-squared:  0.9467,    Adjusted R-squared:  0.9086 
F-statistic: 24.85 on 5 and 7 DF,  p-value: 0.0002553

Analysis of Variance Table

Response: y2
            Df  Sum Sq Mean Sq F value    Pr(>F)
FO(x1, x2)   2  24.357  12.178  4.4242   0.05726
TWI(x1, x2)  1   6.250   6.250  2.2706   0.17558
PQ(x1, x2)   2 311.356 155.678 56.5561 4.778e-05
Residuals    7  19.268   2.753                  
Lack of fit  3   9.268   3.089  1.2358   0.40668
Pure error   4  10.000   2.500                  

Stationary point of response surface:
        x1         x2 
-1.3566432  0.1680434 

Stationary point in original units:
        a         b 
 78.21678 175.84022 

Eigenanalysis:
eigen() decomposition
$values
[1] -0.6230873 -6.7519127

$vectors
         [,1]      [,2]
x1 -0.9947312 0.1025173
x2  0.1025173 0.9947312
funcao2<-rsm(y2~FO(x1,x2)+PQ(x2),data = DOE)
summary(funcao2)
Near-stationary-ridge situation detected -- stationary point altered
 Change 'threshold' if this is not what you intend

Call:
rsm(formula = y2 ~ FO(x1, x2) + PQ(x2), data = DOE)

            Estimate Std. Error  t value  Pr(>|t|)    
(Intercept) 69.52174    0.64613 107.5971 2.626e-15 ***
x1          -1.65533    0.63253  -2.6170   0.02795 *  
x2           0.55178    0.63253   0.8723   0.40569    
x2^2        -6.59783    0.67251  -9.8107 4.194e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple R-squared:  0.9203,    Adjusted R-squared:  0.8937 
F-statistic: 34.62 on 3 and 9 DF,  p-value: 2.856e-05

Analysis of Variance Table

Response: y2
            Df  Sum Sq Mean Sq F value    Pr(>F)
FO(x1, x2)   2  24.357  12.178  3.8049   0.06345
PQ(x2)       1 308.068 308.068 96.2496 4.194e-06
Residuals    9  28.806   3.201                  
Lack of fit  5  18.806   3.761  1.5045   0.35674
Pure error   4  10.000   2.500                  

Stationary point of response surface:
        x1         x2 
0.00000000 0.04181504 

Stationary point in original units:
       a        b 
 85.0000 175.2091 

Eigenanalysis:
eigen() decomposition
$values
[1]  0.000000 -6.597826

$vectors
   [,1] [,2]
x1   -1    0
x2    0    1
curve(dnorm(x,mean(funcao2$residuals),sd(funcao2$residuals)),
      from = min(funcao2$residuals),to=max(funcao2$residuals),
      col="blue",lwd=2,ylab = "Densidade",main="Curva Normal")

qqnorm(funcao2$residuals,lwd=2,col="turquoise")
qqline(funcao2$residuals,col="red",lwd=2)

shapiro.test(funcao2$residuals)

    Shapiro-Wilk normality test

data:  funcao2$residuals
W = 0.96382, p-value = 0.8113
ad.test(funcao2$residuals)

    Anderson-Darling normality test

data:  funcao2$residuals
A = 0.26901, p-value = 0.6184
residualPlot(funcao2)

bptest(funcao2)

    studentized Breusch-Pagan test

data:  funcao2
BP = 0.44484, df = 3, p-value = 0.9308
dwt(funcao2)
Near-stationary-ridge situation detected -- stationary point altered
 Change 'threshold' if this is not what you intend
 lag Autocorrelation D-W Statistic p-value
   1      -0.2992774      2.428138    0.46
 Alternative hypothesis: rho != 0
dwtest(funcao2,alternative = c("two.sided"))

    Durbin-Watson test

data:  funcao2
DW = 2.4281, p-value = 0.453
alternative hypothesis: true autocorrelation is not 0
bgtest(funcao2,order = 1)

    Breusch-Godfrey test for serial correlation of order up
    to 1

data:  funcao2
LM test = 1.4251, df = 1, p-value = 0.2326
persp(funcao2,~x1+x2,zlab = "y",col=rainbow(510),
      contours = "color")

contour(funcao2,~x1+x2,image = T)

contour(funcao2,~x1+x2,col=topo.colors(12))

funcao3<-rsm(y3~FO(x1,x2),data = DOE)
summary(funcao3)

Call:
rsm(formula = y3 ~ FO(x1, x2), data = DOE)

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept) 3386.154     44.163 76.6745 3.477e-15 ***
x1           152.604     56.297  2.7107  0.021907 *  
x2           229.853     56.297  4.0829  0.002204 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple R-squared:  0.706, Adjusted R-squared:  0.6472 
F-statistic: 12.01 on 2 and 10 DF,  p-value: 0.002195

Analysis of Variance Table

Response: y3
            Df Sum Sq Mean Sq F value   Pr(>F)
FO(x1, x2)   2 608963  304481 12.0090 0.002195
Residuals   10 253545   25355                 
Lack of fit  6 187825   31304  1.9053 0.277343
Pure error   4  65720   16430                 

Direction of steepest ascent (at radius 1):
       x1        x2 
0.5531155 0.8331046 

Corresponding increment in original units:
       a        b 
2.765578 4.165523 
curve(dnorm(x,mean(funcao3$residuals),sd(funcao3$residuals)),
      from = min(funcao3$residuals),to=max(funcao3$residuals),
      col="blue",lwd=2,ylab = "Densidade",main="Curva Normal")

qqnorm(funcao3$residuals,lwd=2,col="turquoise")
qqline(funcao3$residuals,col="red",lwd=2)

shapiro.test(funcao3$residuals)

    Shapiro-Wilk normality test

data:  funcao3$residuals
W = 0.94724, p-value = 0.5571
ad.test(funcao3$residuals)

    Anderson-Darling normality test

data:  funcao3$residuals
A = 0.331, p-value = 0.4618
residualPlot(funcao3)

bptest(funcao3)

    studentized Breusch-Pagan test

data:  funcao3
BP = 1.1001, df = 2, p-value = 0.5769
dwt(funcao3)
 lag Autocorrelation D-W Statistic p-value
   1      0.06979422      1.818395   0.884
 Alternative hypothesis: rho != 0
dwtest(funcao3,alternative = c("two.sided"))

    Durbin-Watson test

data:  funcao3
DW = 1.8184, p-value = 0.8827
alternative hypothesis: true autocorrelation is not 0
bgtest(funcao3,order = 1)

    Breusch-Godfrey test for serial correlation of order up
    to 1

data:  funcao3
LM test = 0.12859, df = 1, p-value = 0.7199
persp(funcao3,~x1+x2,zlab = "y",col=rainbow(510),
      contours = "color")

contour(funcao3,~x1+x2,image = T)

contour(funcao3,~x1+x2,col=topo.colors(12))

LS0tDQp0aXRsZTogIlJTTV9CT1hfQkVITktFTiINCm91dHB1dDogaHRtbF9ub3RlYm9vaw0KLS0tDQpFeHBlcmltZW50byBSU00gZGVzaWduIENDRCAtIFJvdGFjaW9uYWw6IExpdnJvIE1vbnRnb21lcnkgKERPRSk6IEV4ZW1wbG8gMTEuMg0KLS0tDQogDQoNCmBgYHtyfQ0KbGlicmFyeShyc20pDQpsaWJyYXJ5KGxtdGVzdCkNCmxpYnJhcnkoY2FyKQ0KbGlicmFyeShub3J0ZXN0KQ0KbGlicmFyeShGckYyKQ0KbGlicmFyeShyaW8pDQpET0U8LSBjY2QoYmFzaXMgPSAyLA0KICAgICAgICAgIGFscGhhID0gInJvdGF0YWJsZSIsDQogICAgICAgICAgbjA9Yyg1LDApLA0KICAgICAgICAgIG9uZWJsb2NrID0gVCwNCiAgICAgICAgICByYW5kb21pemUgPSBGLA0KICAgICAgICAgIGNvZGluZyA9IGxpc3QoeDF+KGEtODUpLzUsDQogICAgICAgICAgICAgICAgICAgICAgICB4Mn4oYi0xNzUpLzUpKQ0KRE9FDQoNCkRPRSR5MTwtYyg3Ni41LDc3LDc4LDc5LjUsNzkuOSw4MC4zLDgwLDc5LjcsNzkuOCwNCiAgICAgICAgICA3NS42LDc4LjQsNzcsNzguNSkNCkRPRSR5MjwtYyg2Miw2MCw2Niw1OSw3Miw2OSw2OCw3MCw3MSw3MSw2OCw1Nyw1OCkNCkRPRSR5MzwtYygyOTQwLDM0NzAsMzY4MCwzODkwLDM0ODAsMzIwMCwzNDEwLDMyOTAsMzUwMCwNCiAgICAgICAgICAzMDIwLDMzNjAsMzE1MCwzNjMwKQ0KRE9FDQpmdW5jYW88LSByc20oeTF+U08oeDEseDIpLGRhdGEgPSBET0UpDQpzdW1tYXJ5KGZ1bmNhbykNCmZ1bmNhbzwtcnNtKHkxfkZPKHgxLHgyKStQUSh4MSx4MiksZGF0YSA9IERPRSkNCnN1bW1hcnkoZnVuY2FvKQ0KcGFyKG1mcm93PWMoMSwyKSkNCmN1cnZlKGRub3JtKHgsbWVhbihmdW5jYW8kcmVzaWR1YWxzKSxzZChmdW5jYW8kcmVzaWR1YWxzKSksDQpmcm9tID0gbWluKGZ1bmNhbyRyZXNpZHVhbHMpLHRvPW1heChmdW5jYW8kcmVzaWR1YWxzKSwNCmNvbD0iYmx1ZSIsbHdkPTIseWxhYiA9ICJEZW5zaWRhZGUiLG1haW49IkN1cnZhIE5vcm1hbCIpDQpxcW5vcm0oZnVuY2FvJHJlc2lkdWFscyxsd2Q9Mixjb2w9InR1cnF1b2lzZSIpDQpxcWxpbmUoZnVuY2FvJHJlc2lkdWFscyxjb2w9InJlZCIsbHdkPTIpDQpzaGFwaXJvLnRlc3QoZnVuY2FvJHJlc2lkdWFscykNCmFkLnRlc3QoZnVuY2FvJHJlc2lkdWFscykNCnBhcihtZnJvdz1jKDEsMSkpDQpyZXNpZHVhbFBsb3QoZnVuY2FvKQ0KYnB0ZXN0KGZ1bmNhbykNCmR3dChmdW5jYW8pDQpkd3Rlc3QoZnVuY2FvLGFsdGVybmF0aXZlID0gYygidHdvLnNpZGVkIikpDQpiZ3Rlc3QoZnVuY2FvLG9yZGVyID0gMSkNCnBlcnNwKGZ1bmNhbyx+eDEreDIsemxhYiA9ICJ5Iixjb2w9cmFpbmJvdyg1MTApLA0KICAgICAgY29udG91cnMgPSAiY29sb3IiKQ0KY29udG91cihmdW5jYW8sfngxK3gyLGltYWdlID0gVCkNCmNvbnRvdXIoZnVuY2FvLH54MSt4Mixjb2w9dG9wby5jb2xvcnMoMTIpKQ0KYGBgDQoNCmBgYHtyfQ0KZnVuY2FvMjwtcnNtKHkyflNPKHgxLHgyKSxkYXRhID0gRE9FKQ0Kc3VtbWFyeShmdW5jYW8yKQ0KZnVuY2FvMjwtcnNtKHkyfkZPKHgxLHgyKStQUSh4MiksZGF0YSA9IERPRSkNCnN1bW1hcnkoZnVuY2FvMikNCmN1cnZlKGRub3JtKHgsbWVhbihmdW5jYW8yJHJlc2lkdWFscyksc2QoZnVuY2FvMiRyZXNpZHVhbHMpKSwNCiAgICAgIGZyb20gPSBtaW4oZnVuY2FvMiRyZXNpZHVhbHMpLHRvPW1heChmdW5jYW8yJHJlc2lkdWFscyksDQogICAgICBjb2w9ImJsdWUiLGx3ZD0yLHlsYWIgPSAiRGVuc2lkYWRlIixtYWluPSJDdXJ2YSBOb3JtYWwiKQ0KcXFub3JtKGZ1bmNhbzIkcmVzaWR1YWxzLGx3ZD0yLGNvbD0idHVycXVvaXNlIikNCnFxbGluZShmdW5jYW8yJHJlc2lkdWFscyxjb2w9InJlZCIsbHdkPTIpDQpzaGFwaXJvLnRlc3QoZnVuY2FvMiRyZXNpZHVhbHMpDQphZC50ZXN0KGZ1bmNhbzIkcmVzaWR1YWxzKQ0KcmVzaWR1YWxQbG90KGZ1bmNhbzIpDQpicHRlc3QoZnVuY2FvMikNCmR3dChmdW5jYW8yKQ0KZHd0ZXN0KGZ1bmNhbzIsYWx0ZXJuYXRpdmUgPSBjKCJ0d28uc2lkZWQiKSkNCmJndGVzdChmdW5jYW8yLG9yZGVyID0gMSkNCnBlcnNwKGZ1bmNhbzIsfngxK3gyLHpsYWIgPSAieSIsY29sPXJhaW5ib3coNTEwKSwNCiAgICAgIGNvbnRvdXJzID0gImNvbG9yIikNCmNvbnRvdXIoZnVuY2FvMix+eDEreDIsaW1hZ2UgPSBUKQ0KY29udG91cihmdW5jYW8yLH54MSt4Mixjb2w9dG9wby5jb2xvcnMoMTIpKQ0KYGBgDQpgYGB7cn0NCmZ1bmNhbzM8LXJzbSh5M35GTyh4MSx4MiksZGF0YSA9IERPRSkNCnN1bW1hcnkoZnVuY2FvMykNCmN1cnZlKGRub3JtKHgsbWVhbihmdW5jYW8zJHJlc2lkdWFscyksc2QoZnVuY2FvMyRyZXNpZHVhbHMpKSwNCiAgICAgIGZyb20gPSBtaW4oZnVuY2FvMyRyZXNpZHVhbHMpLHRvPW1heChmdW5jYW8zJHJlc2lkdWFscyksDQogICAgICBjb2w9ImJsdWUiLGx3ZD0yLHlsYWIgPSAiRGVuc2lkYWRlIixtYWluPSJDdXJ2YSBOb3JtYWwiKQ0KcXFub3JtKGZ1bmNhbzMkcmVzaWR1YWxzLGx3ZD0yLGNvbD0idHVycXVvaXNlIikNCnFxbGluZShmdW5jYW8zJHJlc2lkdWFscyxjb2w9InJlZCIsbHdkPTIpDQpzaGFwaXJvLnRlc3QoZnVuY2FvMyRyZXNpZHVhbHMpDQphZC50ZXN0KGZ1bmNhbzMkcmVzaWR1YWxzKQ0KcmVzaWR1YWxQbG90KGZ1bmNhbzMpDQpicHRlc3QoZnVuY2FvMykNCmR3dChmdW5jYW8zKQ0KZHd0ZXN0KGZ1bmNhbzMsYWx0ZXJuYXRpdmUgPSBjKCJ0d28uc2lkZWQiKSkNCmJndGVzdChmdW5jYW8zLG9yZGVyID0gMSkNCnBlcnNwKGZ1bmNhbzMsfngxK3gyLHpsYWIgPSAieSIsY29sPXJhaW5ib3coNTEwKSwNCiAgICAgIGNvbnRvdXJzID0gImNvbG9yIikNCmNvbnRvdXIoZnVuY2FvMyx+eDEreDIsaW1hZ2UgPSBUKQ0KY29udG91cihmdW5jYW8zLH54MSt4Mixjb2w9dG9wby5jb2xvcnMoMTIpKQ0KYGBgDQoNCg==