Aula 1: Dinâmica Catapulta: Turma 1 ano

Aula DOE em Linguagem R Aula 1. Experimento da Catapulta com projetos fatoriais completos 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)

DOE = expand.grid(x1=c(3,6),
                  x2=c(3,5),
                  x3=c(2,4))
DOE<- rbind(DOE,DOE)
DOE$y<- c(1.52,1.28,3.27,3.41,3.39,3.80,5.80,5.48,
      1.63,1.68,3.70,3.62,3.51,3.26,5.36,5.53)
funcao<- lm(y~ x1*x2*x3,data = DOE)
summary(funcao)

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

Residuals:
    Min      1Q  Median      3Q     Max 
-0.2700 -0.1288  0.0000  0.1288  0.2700 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.95750    2.44904  -0.799    0.447
x1          -0.29250    0.51630  -0.567    0.587
x2           0.64250    0.59398   1.082    0.311
x3           0.47500    0.77446   0.613    0.557
x1:x2        0.06750    0.12522   0.539    0.605
x1:x3        0.09917    0.16327   0.607    0.560
x2:x3        0.12500    0.18783   0.665    0.524
x1:x2:x3    -0.02333    0.03960  -0.589    0.572

Residual standard error: 0.2376 on 8 degrees of freedom
Multiple R-squared:  0.9862,    Adjusted R-squared:  0.9741 
F-statistic: 81.66 on 7 and 8 DF,  p-value: 8.28e-07

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

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

Residuals:
     Min       1Q   Median       3Q      Max 
-0.29303 -0.14882 -0.01256  0.12763  0.27027 

Coefficients:
         Estimate Std. Error t value Pr(>|t|)   
x1       -0.68400    0.15996  -4.276  0.00206 **
x2        0.18191    0.14114   1.289  0.22959   
x3       -0.11225    0.23994  -0.468  0.65103   
x1:x2     0.15962    0.04798   3.327  0.00884 **
x1:x3     0.21662    0.06973   3.107  0.01259 * 
x2:x3     0.26318    0.07197   3.657  0.00526 **
x1:x2:x3 -0.05097    0.01891  -2.695  0.02460 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2328 on 9 degrees of freedom
Multiple R-squared:  0.9979,    Adjusted R-squared:  0.9962 
F-statistic: 606.2 on 7 and 9 DF,  p-value: 2.663e-11

funcao2<- lm(y~ x1*x2*x3,data = DOE)

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

$fit
       fit      lwr      upr
1 1.667118 1.074934 2.259302

$se.fit
[1] 0.1197625

$df
[1] 9

$residual.scale
[1] 0.2327765

1.52- predict(funcao,newdata = data.frame(x1=3,
                                    x2=3,
                                    x3=2)) # Resíduos

         1 
-0.1471176 

curve(dnorm(x,mean(funcao$residuals),sd(funcao$residuals)),
      from = min(funcao$residuals),to=max(funcao$residuals),
      col="blue",lwd=3,main="Curva dos Resíduos",
      xlab = "Observações",ylab = "Densidade")
lines(density(funcao$residuals),col="red",lwd=3)

legend("topleft",  legend = c("Real","Estimada"),
       col = c("blue","red"),lwd=2)

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

shapiro.test(funcao$residuals)

    Shapiro-Wilk normality test

data:  funcao$residuals
W = 0.94628, p-value = 0.4332

residualPlot(funcao)

ols_test_breusch_pagan(funcao2,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 x1:x2 x1:x3 x2:x3 x1:x2:x3 

         Test Summary (Unadjusted p values)        
 ------------------------------------------------
  Variable            chi2        df        p     
 ------------------------------------------------
  x1               0.085183201     1    0.7703924 
  x2               0.027546972     1    0.8681784 
  x3               0.100037971     1    0.7517841 
  x1:x2            0.111293116     1    0.7386767 
  x1:x3            0.225922519     1    0.6345638 
  x2:x3            0.002306547     1    0.9616951 
  x1:x2:x3         0.014828932     1    0.9030779 
 ------------------------------------------------
  simultaneous     6.482851250     7    0.4846300 
 ------------------------------------------------

bptest(funcao)

    studentized Breusch-Pagan test

data:  funcao
BP = 12.4, df = 6, p-value = 0.05361

dwt(funcao)

 lag Autocorrelation D-W Statistic p-value
   1       0.4736006      1.005327   0.082
 Alternative hypothesis: rho != 0

dwtest(funcao)

    Durbin-Watson test

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

durbinWatsonTest(funcao)

 lag Autocorrelation D-W Statistic p-value
   1       0.4736006      1.005327   0.098
 Alternative hypothesis: rho != 0

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

MEPlot(funcao)

persp(funcao,~x1+x2,zlab="y",col = rainbow(100),
      contours ="colors")

persp(funcao,~x1+x3,zlab="y",col = topo.colors(10),
      contours ="colors")

persp(funcao,~x2+x3,zlab="y",col = rainbow(100),
      contours ="colors")

contour(funcao,~x1+x2,image = T)
points(3,5,col="red",lwd=5)

contour(funcao,~x1+x3,col=topo.colors(10))
points(3,3,col="red",lwd=5)

contour(funcao,~x2+x3,col=rainbow(10))
points(3,3,col="red",lwd=5)

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

$par
[1] 3 5 4

$value
[1] -5.607635

$counts
function gradient 
       5        5 

$convergence
[1] 0

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