DoE Fatorial Completo

Dinâmica da Catapulta

library(FrF2)
library(rsm)
library(lmtest)
library(car)
library(Rsolnp)

DOE<-expand.grid(x1=c(3,6),
                 x2=c(3,5),
                 x3=c(2,4))
DOE

DOE<- rbind(DOE,DOE) #Replicação
DOE$Bloco<- as.factor(c(rep(1,8),rep(2,8)))
DOE

DOE$y<-c(284,258,437,420,547,532,707,663,
         264,239,445,425,521,522,691,628)
DOE

funcao<-lm(y~ x1*x2*x3+Bloco,data = DOE)
summary(funcao)

Call:
lm.default(formula = y ~ x1 * x2 * x3 + Bloco, data = DOE)

Residuals:
    Min      1Q  Median      3Q     Max 
-11.062  -3.688   0.000   3.688  11.062 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  -97.687    107.113  -0.912   0.3921  
x1           -44.917     22.575  -1.990   0.0869 .
x2            54.250     25.971   2.089   0.0751 .
x3            82.125     33.862   2.425   0.0457 *
Bloco2       -14.125      5.194  -2.719   0.0298 *
x1:x2         10.083      5.475   1.842   0.1081  
x1:x3         16.458      7.139   2.305   0.0545 .
x2:x3         12.875      8.213   1.568   0.1609  
x1:x2:x3      -4.458      1.731  -2.575   0.0367 *
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.39 on 7 degrees of freedom
Multiple R-squared:  0.9979,    Adjusted R-squared:  0.9955 
F-statistic: 420.3 on 8 and 7 DF,  p-value: 1.09e-08

summary(aov(funcao))

            Df Sum Sq Mean Sq  F value   Pr(>F)    
x1           1   2730    2730   25.297  0.00151 ** 
x2           1  97500   97500  903.451 1.16e-08 ***
x3           1 259845  259845 2407.764 3.82e-10 ***
Bloco        1    798     798    7.395  0.02979 *  
x1:x2        1    390     390    3.614  0.09904 .  
x1:x3        1     68      68    0.631  0.45318    
x2:x3        1    827     827    7.659  0.02779 *  
x1:x2:x3     1    716     716    6.631  0.03674 *  
Residuals    7    755     108                      
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

funcao<-lm(y~ x1+x2+x3+Bloco+x2*x3,data = DOE)
summary(funcao)

Call:
lm.default(formula = y ~ x1 + x2 + x3 + Bloco + x2 * x3, data = DOE)

Residuals:
    Min      1Q  Median      3Q     Max 
-24.125  -4.156   0.250   8.531  14.625 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -260.625     46.586  -5.594 0.000229 ***
x1            -8.708      2.315  -3.762 0.003711 ** 
x2            99.625     10.980   9.073 3.85e-06 ***
x3           156.187     14.317  10.909 7.12e-07 ***
Bloco2       -14.125      6.945  -2.034 0.069346 .  
x2:x3         -7.187      3.472  -2.070 0.065290 .  
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 13.89 on 10 degrees of freedom
Multiple R-squared:  0.9947,    Adjusted R-squared:  0.992 
F-statistic:   375 on 5 and 10 DF,  p-value: 4.897e-11

summary(aov(funcao))

            Df Sum Sq Mean Sq  F value   Pr(>F)    
x1           1   2730    2730   14.152  0.00371 ** 
x2           1  97500   97500  505.411 6.82e-10 ***
x3           1 259845  259845 1346.958 5.37e-12 ***
Bloco        1    798     798    4.137  0.06935 .  
x2:x3        1    827     827    4.285  0.06529 .  
Residuals   10   1929     193                      
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#=================== Teste de Normalidade=======================================
hist(funcao$residuals,lwd=4,freq = F,col = "magenta",
     main="Resíduos")
curve(dnorm(x,mean(funcao$residuals),sd(funcao$residuals)),
      add = T,col="blue",lwd=6)

curve(dnorm(x,mean(funcao$residuals),sd(funcao$residuals)),
      from = min(funcao$residuals),to=max(funcao$residuals),
      col="purple", main="Curva normal dos resíduos")

  qqnorm(funcao$residuals,col="blue",lwd=4)
  qqline(funcao$residuals,col="red",lwd=2)

shapiro.test(funcao$residuals)

    Shapiro-Wilk normality test

data:  funcao$residuals
W = 0.93796, p-value = 0.3248

library(FrF2)
library(rsm)
library(lmtest)
library(olsrr)
library(car)
residualPlots(funcao)

           Test stat Pr(>|Test stat|)  
x1           -0.2148          0.83471  
x2           -2.2174          0.05379 .
x3           -0.0827          0.93593  
Bloco                                  
Tukey test    1.3075          0.19105  
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

residualPlot(funcao)

ols_test_breusch_pagan(funcao,rhs = T,multiple = T)

 Breusch Pagan Test for Heteroskedasticity
 -----------------------------------------
 Ho: the variance is constant            
 Ha: the variance is not constant        

              Data               
 --------------------------------
 Response : y 
 Variables: x1 x2 x3 Bloco2 x2:x3 

         Test Summary (Unadjusted p values)        
 ------------------------------------------------
  Variable            chi2       df        p      
 ------------------------------------------------
  x1               0.01622404     1    0.89864471 
  x2               2.06663816     1    0.15055307 
  x3               1.73542325     1    0.18771978 
  Bloco2           1.38914089     1    0.23855024 
  x2:x3            3.60532371     1    0.05759485 
 ------------------------------------------------
  simultaneous     5.21736635     5    0.38993270 
 ------------------------------------------------

bptest(funcao)

    studentized Breusch-Pagan test

data:  funcao
BP = 7.08, df = 5, p-value = 0.2148

dwt(funcao)

 lag Autocorrelation D-W Statistic p-value
   1      -0.1763429      2.047415   0.634
 Alternative hypothesis: rho != 0

dwtest(funcao)

    Durbin-Watson test

data:  funcao
DW = 2.0474, p-value = 0.3185
alternative hypothesis: true autocorrelation is greater than 0

bgtest(funcao)

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

data:  funcao
LM test = 1.3679, df = 1, p-value = 0.2422

#========= Gráficos de Superfície e de Contorno=================================
persp(funcao,~x1+x2,zlab = "Distância",col=rainbow(250),
      contours = "colors")

persp(funcao,~x1+x3,zlab = "Distância",col=rainbow(250),
      contours = "colors")

persp(funcao,~x2+x3,zlab = "Distância",col=rainbow(250),
      contours = "colors")

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

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

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

#========= Gráfico das interações ==============================================
IAPlot(funcao)

#===== Gráfico das médias dos efeitos ==========================================
funcao2<-lm(y~x1+x2+x3+Bloco,data = DOE)
MEPlot(funcao2)

#========== Predição ===========================================================
predict(funcao,newdata = data.frame(x1=3,
                    x2=5,
                    x3=4,
 Bloco="1"),se.fit = T,interval = "pred")

$fit
      fit      lwr      upr
1 692.375 656.0861 728.6639

$se.fit
[1] 8.505421

$df
[1] 10

$residual.scale
[1] 13.88929

predict(funcao,newdata = data.frame(x1=3,
                                    x2=5,
                                    x3=4,
    Bloco="2"),se.fit = T,interval = "pred")

$fit
     fit      lwr      upr
1 678.25 641.9611 714.5389

$se.fit
[1] 8.505421

$df
[1] 10

$residual.scale
[1] 13.88929

obj=function(x){a<-
  predict(funcao,newdata = data.frame(x1=x[1],
                                      x2=x[2],
                                      x3=x[3],
                                      Bloco="2"))
return(-a)}
x0<-c(0,0,0)
optim(par = x0,fn=obj,lower = c(3,3,2),
      upper = c(6,5,4),method = "L-BFGS-B")

$par
[1] 3 5 4

$value
[1] -678.25

$counts
function gradient 
       2        2 

$convergence
[1] 0

$message
[1] "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"

#======= Usando pacotes do R para DoE fatorial =================================

DOE<- FrF2(nruns = 8,
           nfactors = 3,
           factor.names = c("x1","x2","x3"),
           replications = 2,
           randomize = F,
           blocks = 1)

creating full factorial with 8 runs ...

  
  DOE$y<-c(284,258,437,420,547,532,707,663,
           264,239,445,425,521,522,691,628)
DOE

class=design, type= full factorial 
NOTE: columns run.no and run.no.std.rp 
 are annotation,  not part of the data frame

funcao3<-lm(y~x1*x2*x3+Blocks,data = DOE)
summary(funcao3)

Call:
lm.default(formula = y ~ x1 * x2 * x3 + Blocks, data = DOE)

Residuals:
    Min      1Q  Median      3Q     Max 
-11.062  -3.688   0.000   3.688  11.062 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  481.000      3.673 130.960 3.99e-13 ***
x11          -13.062      2.597  -5.030  0.00151 ** 
x21           78.062      2.597  30.057 1.16e-08 ***
x31          127.438      2.597  49.069 3.82e-10 ***
Blocks.2     -14.125      5.194  -2.719  0.02979 *  
x11:x21       -4.938      2.597  -1.901  0.09904 .  
x11:x31       -2.062      2.597  -0.794  0.45318    
x21:x31       -7.187      2.597  -2.767  0.02779 *  
x11:x21:x31   -6.687      2.597  -2.575  0.03674 *  
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.39 on 7 degrees of freedom
Multiple R-squared:  0.9979,    Adjusted R-squared:  0.9955 
F-statistic: 420.3 on 8 and 7 DF,  p-value: 1.09e-08

funcao3<-lm(y~x1+x2+x3+Blocks+x2*x3,data = DOE)
summary(funcao3)

Call:
lm.default(formula = y ~ x1 + x2 + x3 + Blocks + x2 * x3, data = DOE)

Residuals:
    Min      1Q  Median      3Q     Max 
-24.125  -4.156   0.250   8.531  14.625 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  481.000      4.911  97.951 3.01e-16 ***
x11          -13.062      3.472  -3.762  0.00371 ** 
x21           78.062      3.472  22.481 6.82e-10 ***
x31          127.438      3.472  36.701 5.37e-12 ***
Blocks.2     -14.125      6.945  -2.034  0.06935 .  
x21:x31       -7.187      3.472  -2.070  0.06529 .  
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 13.89 on 10 degrees of freedom
Multiple R-squared:  0.9947,    Adjusted R-squared:  0.992 
F-statistic:   375 on 5 and 10 DF,  p-value: 4.897e-11

IAPlot(funcao3)

funcao4<-lm(y~x1+x2+x3+Blocks,data = DOE) 
MEPlot(funcao4)

DOE<- FrF2(nruns = 8,
           nfactors = 3,
           factor.names = list(x1=c(3,6),
                                 x2=c(3,5),
                                 x3=c(2,4)),
           replications = 2,
           randomize = F,
           blocks = 1)

creating full factorial with 8 runs ...

DOE

class=design, type= full factorial 
NOTE: columns run.no and run.no.std.rp 
 are annotation,  not part of the data frame

  
  DOE$y<-c(284,258,437,420,547,532,707,663,
           264,239,445,425,521,522,691,628)
DOE

class=design, type= full factorial 
NOTE: columns run.no and run.no.std.rp 
 are annotation,  not part of the data frame

funcao3<-lm(y~x1*x2*x3+Blocks,data = DOE)
summary(funcao3)

Call:
lm.default(formula = y ~ x1 * x2 * x3 + Blocks, data = DOE)

Residuals:
    Min      1Q  Median      3Q     Max 
-11.062  -3.688   0.000   3.688  11.062 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  481.000      3.673 130.960 3.99e-13 ***
x11          -13.062      2.597  -5.030  0.00151 ** 
x21           78.062      2.597  30.057 1.16e-08 ***
x31          127.438      2.597  49.069 3.82e-10 ***
Blocks.2     -14.125      5.194  -2.719  0.02979 *  
x11:x21       -4.938      2.597  -1.901  0.09904 .  
x11:x31       -2.062      2.597  -0.794  0.45318    
x21:x31       -7.187      2.597  -2.767  0.02779 *  
x11:x21:x31   -6.687      2.597  -2.575  0.03674 *  
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.39 on 7 degrees of freedom
Multiple R-squared:  0.9979,    Adjusted R-squared:  0.9955 
F-statistic: 420.3 on 8 and 7 DF,  p-value: 1.09e-08

funcao3<-lm(y~x1+x2+x3+Blocks+x2*x3,data = DOE)
summary(funcao3)

Call:
lm.default(formula = y ~ x1 + x2 + x3 + Blocks + x2 * x3, data = DOE)

Residuals:
    Min      1Q  Median      3Q     Max 
-24.125  -4.156   0.250   8.531  14.625 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  481.000      4.911  97.951 3.01e-16 ***
x11          -13.062      3.472  -3.762  0.00371 ** 
x21           78.062      3.472  22.481 6.82e-10 ***
x31          127.438      3.472  36.701 5.37e-12 ***
Blocks.2     -14.125      6.945  -2.034  0.06935 .  
x21:x31       -7.187      3.472  -2.070  0.06529 .  
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 13.89 on 10 degrees of freedom
Multiple R-squared:  0.9947,    Adjusted R-squared:  0.992 
F-statistic:   375 on 5 and 10 DF,  p-value: 4.897e-11

IAPlot(funcao3)

funcao4<-lm(y~x1+x2+x3+Blocks,data = DOE) 
MEPlot(funcao4)