Delineamento de Experimentos: Projetos Fatoriais

Aula 1 da disciplina métodos de otimização e simulação - MOS do Programa de Pós-Graduação em Engenharia unesp-FEG Guaratinguetá-SP: DoE planejamento fatorial completo com 2 níveis”

library(car) # Pacotes para análise estatística
library(rsm) # Pacotes para Metodologia da superfície
#de resposta - RSM
library(FrF2) #pacote para DOE fatorial
library(nortest) #pacote para normalidade
library(olsrr) #pacote para homocedasticidade
library(lmtest)
library(rio)
library(nortest)
library(Rsolnp) # Otimização não linear

DOE<- expand.grid(x1=c(-1,1),
                  x2=c(-1,1))
DOE

DOE$y<-c(500,1100,900,1900)
DOE

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

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

Residuals:
ALL 4 residuals are 0: no residual degrees of freedom!

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)     1100        NaN     NaN      NaN
x1               400        NaN     NaN      NaN
x2               300        NaN     NaN      NaN
x1:x2            100        NaN     NaN      NaN

Residual standard error: NaN on 0 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:    NaN 
F-statistic:   NaN on 3 and 0 DF,  p-value: NA

# DOE com Replicação
DOE<- expand.grid(x1=c(-1,1),
                  x2=c(-1,1))
DOE

DOE<- rbind(DOE,DOE) # Replicação
DOE

DOE$y<-c(500,1100,900,1900,
         510,1000,940,1920)
funcao<- lm(y~ x1*x2,data = DOE)
summary(funcao)

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

Residuals:
  1   2   3   4   5   6   7   8 
 -5  50 -20 -10   5 -50  20  10 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1096.25      13.75  79.727 1.48e-07 ***
x1            383.75      13.75  27.909 9.81e-06 ***
x2            318.75      13.75  23.182 2.05e-05 ***
x1:x2         111.25      13.75   8.091  0.00127 ** 
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 38.89 on 4 degrees of freedom
Multiple R-squared:  0.9971,    Adjusted R-squared:  0.9949 
F-statistic: 460.6 on 3 and 4 DF,  p-value: 1.561e-05

summary(aov(funcao))

            Df  Sum Sq Mean Sq F value   Pr(>F)    
x1           1 1178113 1178113  778.92 9.81e-06 ***
x2           1  812812  812812  537.40 2.05e-05 ***
x1:x2        1   99012   99012   65.46  0.00127 ** 
Residuals    4    6050    1513                     
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

hist(funcao$residuals,freq = F,
     col="blue", lwd=2, 
     main="Histograma dos Resídos",
     xlab = "Observações")
curve(dnorm(x,mean(funcao$residuals),
            sd(funcao$residuals)),add =T,
      col="magenta", lwd=4)

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

shapiro.test(funcao$residuals)

    Shapiro-Wilk normality test

data:  funcao$residuals
W = 0.98574, p-value = 0.9856

library(car)
library(lmtest)
library(olsrr)
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 x1:x2 

        Test Summary (Unadjusted p values)       
 ----------------------------------------------
  Variable           chi2      df        p      
 ----------------------------------------------
  x1               2.067892     1    0.15042936 
  x2               1.792501     1    0.18062180 
  x1:x2            3.366164     1    0.06654856 
 ----------------------------------------------
  simultaneous     7.226556     3    0.06501674 
 ----------------------------------------------

bptest(funcao)

    studentized Breusch-Pagan test

data:  funcao
BP = 8, df = 3, p-value = 0.04601

dwt(funcao)

 lag Autocorrelation D-W Statistic p-value
   1      -0.3553719      2.690083   0.198
 Alternative hypothesis: rho != 0

dwtest(funcao)

    Durbin-Watson test

data:  funcao
DW = 2.6901, p-value = 0.8778
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.036, df = 1, p-value = 0.3088

library(FrF2)
funcao2<- lm(y~ x1+x2,data = DOE)
MEPlot(funcao2)

funcao<- lm(y~ x1*x2,data = DOE)
IAPlot(funcao)

predict(funcao, newdata = data.frame(
  x1=1, x2=1), se.fit = T,
  interval = "conf")

$fit
   fit      lwr      upr
1 1910 1833.648 1986.352

$se.fit
[1] 27.5

$df
[1] 4

$residual.scale
[1] 38.89087

predict(funcao, newdata = data.frame(
  x1=1, x2=1), se.fit = T,
  interval = "pred")

$fit
   fit      lwr      upr
1 1910 1777.754 2042.246

$se.fit
[1] 27.5

$df
[1] 4

$residual.scale
[1] 38.89087

#========== Otimização=========================================================
obj=function(x){a<- predict(funcao,
newdata = data.frame(x1=x[1],
                     x2=x[2]))
return(-a)}
x0<-c(0,0)
optim(par = x0,fn=obj,lower = -1,upper = 1,
      method = "L-BFGS-B")

$par
[1] 1 1

$value
[1] -1910

$counts
function gradient 
       2        2 

$convergence
[1] 0

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

solnp(pars = x0,fun = obj,LB = rep(-1,2),
      UB = rep(1,2))

Iter: 1 fn: -1910.0000   Pars:  1.00000 1.00000
Iter: 2 fn: -1910.0000   Pars:  1.00000 1.00000
solnp--> Completed in 2 iterations
$pars
[1] 1 1

$convergence
[1] 0

$values
[1] -1096.25 -1910.00 -1910.00

$lagrange
[1] 0

$hessian
     [,1] [,2]
[1,]    1    0
[2,]    0    1

$ineqx0
NULL

$nfuneval
[1] 100

$outer.iter
[1] 2

$elapsed
Time difference of 0.04138803 secs

$vscale
[1] 1 1 1

#========== Gráficos de Superfície e Contorno===================================
library(rsm)
persp(funcao, ~x1+x2,zlab = "y", col=rainbow(50),
      contours = "colors",lwd=1)

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

#================ Gerando Experimentos com Pacotes =============================
library(FrF2)
DOE<-FrF2(nruns = 4,
          nfactors = 2,
          factor.names = c("x1","x2"),
          randomize = F,
          replications = 2)

creating full factorial with 4 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(500,1100,900,1900,
         510,1000,940,1920)
DOE

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

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

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

Residuals:
  1   2   3   4   5   6   7   8 
 -5  50 -20 -10   5 -50  20  10 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1096.25      13.75  79.727 1.48e-07 ***
x11           383.75      13.75  27.909 9.81e-06 ***
x21           318.75      13.75  23.182 2.05e-05 ***
x11:x21       111.25      13.75   8.091  0.00127 ** 
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 38.89 on 4 degrees of freedom
Multiple R-squared:  0.9971,    Adjusted R-squared:  0.9949 
F-statistic: 460.6 on 3 and 4 DF,  p-value: 1.561e-05

summary(aov(funcao))

            Df  Sum Sq Mean Sq F value   Pr(>F)    
x1           1 1178113 1178113  778.92 9.81e-06 ***
x2           1  812812  812812  537.40 2.05e-05 ***
x1:x2        1   99012   99012   65.46  0.00127 ** 
Residuals    4    6050    1513                     
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Esse comando é para inserir as variáveis naturais
DOE<- FrF2(nruns = 4,
     factor.names = list(x1=c(1,2),
                         x2=c(1,4)),
     randomize = F,
     replications = 2)

creating full factorial with 4 runs ...

DOE$y<-c(500,1100,900,1900,
         510,1000,940,1920)

DOE

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

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

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

Residuals:
  1   2   3   4   5   6   7   8 
 -5  50 -20 -10   5 -50  20  10 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1096.25      13.75  79.727 1.48e-07 ***
x11           383.75      13.75  27.909 9.81e-06 ***
x21           318.75      13.75  23.182 2.05e-05 ***
x11:x21       111.25      13.75   8.091  0.00127 ** 
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 38.89 on 4 degrees of freedom
Multiple R-squared:  0.9971,    Adjusted R-squared:  0.9949 
F-statistic: 460.6 on 3 and 4 DF,  p-value: 1.561e-05

summary(aov(funcao))

            Df  Sum Sq Mean Sq F value   Pr(>F)    
x1           1 1178113 1178113  778.92 9.81e-06 ***
x2           1  812812  812812  537.40 2.05e-05 ***
x1:x2        1   99012   99012   65.46  0.00127 ** 
Residuals    4    6050    1513                     
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

library(lmtest)
library(car)
residualPlots(funcao)

residualPlot(funcao)

curve(dnorm(x,mean(funcao$residuals),sd(funcao$residuals)),
      from = min(funcao$residuals),to=max(funcao$residuals),
      col="purple",lwd=3,ylab="Densidade",
      xlab = "Observações")
lines(density(funcao$residuals),col="blue",lwd=2)

legend("topleft",legend = c("Real","Estimada"),
       col=c("purple","blue"),lwd=3)