DOE COM BLOCAGEM: CATAPULTA_2024

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

DOE<-rbind(DOE,DOE)
DOE

DOE$Bloco<- as.factor(c(rep(1,8),rep(2,8)))

DOE

DOE$y<-c(1.71,  1.83,3.96, 3.54,4.22, 4.42, 5.80, 5.60,
         1.72,  1.87,3.83, 3.57,4.18, 4.40, 5.90, 5.90 )
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 
-0.13187 -0.02906  0.00000  0.02906  0.13187 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -6.38312    1.19684  -5.333 0.001084 ** 
x1           0.51000    0.29023   1.757 0.122293    
x2           1.83500    0.29023   6.323 0.000395 ***
x3           1.76937    0.37841   4.676 0.002272 ** 
Bloco2       0.03625    0.04452   0.814 0.442309    
x1:x2       -0.16000    0.07039  -2.273 0.057220 .  
x1:x3       -0.04312    0.09178  -0.470 0.652724    
x2:x3       -0.19437    0.09178  -2.118 0.071946 .  
x1:x2:x3     0.02062    0.02226   0.927 0.384980    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.08904 on 7 degrees of freedom
Multiple R-squared:  0.9983,    Adjusted R-squared:  0.9964 
F-statistic: 525.8 on 8 and 7 DF,  p-value: 4.992e-09

summary(aov(funcao))

            Df Sum Sq Mean Sq  F value   Pr(>F)    
x1           1  0.002   0.002    0.285  0.61021    
x2           1 11.816  11.816 1490.525 2.04e-09 ***
x3           1 21.137  21.137 2666.229 2.68e-10 ***
Bloco        1  0.005   0.005    0.663  0.44231    
x1:x2        1  0.154   0.154   19.433  0.00313 ** 
x1:x3        1  0.025   0.025    3.129  0.12023    
x2:x3        1  0.200   0.200   25.260  0.00152 ** 
x1:x2:x3     1  0.007   0.007    0.859  0.38498    
Residuals    7  0.055   0.008                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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.22524 -0.09830 -0.03607  0.11664  0.19714 

Coefficients:
         Estimate Std. Error t value Pr(>|t|)    
x1       -0.98765    0.14069  -7.020 6.19e-05 ***
x2        0.33735    0.14069   2.398 0.040043 *  
x3       -0.14013    0.23918  -0.586 0.572368    
x1:x2     0.19239    0.04754   4.047 0.002900 ** 
x1:x3     0.40617    0.07174   5.662 0.000309 ***
x2:x3     0.25492    0.07174   3.553 0.006183 ** 
x1:x2:x3 -0.08509    0.02003  -4.247 0.002151 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.178 on 9 degrees of freedom
Multiple R-squared:  0.999, Adjusted R-squared:  0.9982 
F-statistic:  1249 on 7 and 9 DF,  p-value: 1.037e-12

summary(aov(funcao))

          Df Sum Sq Mean Sq  F value   Pr(>F)    
x1         1 229.06  229.06 7232.437 2.18e-14 ***
x2         1  26.38   26.38  832.862 3.51e-10 ***
x3         1  17.98   17.98  567.797 1.93e-09 ***
x1:x2      1   2.35    2.35   74.073 1.23e-05 ***
x1:x3      1   0.53    0.53   16.603  0.00278 ** 
x2:x3      1   0.00    0.00    0.002  0.96373    
x1:x2:x3   1   0.57    0.57   18.039  0.00215 ** 
Residuals  9   0.29    0.03                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Resíduos
# Análise da normalidade dos resíduos
curve(dnorm(x,mean(funcao$residuals),sd(funcao$residuals)),
      from = min(funcao$residuals),to=max(funcao$residuals),
      lwd=3,col="violet",main="Curva Normal dos Resíduos",
      xlab="Observações",ylab = "Densidade")

shapiro.test(funcao$residuals)

    Shapiro-Wilk normality test

data:  funcao$residuals
W = 0.9352, p-value = 0.2943

# Testes de homocedasticidade e autocorrelação residual
library(lmtest)
library(car)
residualPlot(funcao)

bptest(funcao)

    studentized Breusch-Pagan test

data:  funcao
BP = 10.47, df = 6, p-value = 0.1062

durbinWatsonTest(funcao)

 lag Autocorrelation D-W Statistic p-value
   1      -0.1387733      1.973388    0.95
 Alternative hypothesis: rho != 0

dwtest(funcao)

    Durbin-Watson test

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

dwt(funcao)

 lag Autocorrelation D-W Statistic p-value
   1      -0.1387733      1.973388   0.978
 Alternative hypothesis: rho != 0

#Análise das médias dos efeitos
library(FrF2)
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.22524 -0.09830 -0.03607  0.11664  0.19714 

Coefficients:
         Estimate Std. Error t value Pr(>|t|)    
x1       -0.98765    0.14069  -7.020 6.19e-05 ***
x2        0.33735    0.14069   2.398 0.040043 *  
x3       -0.14013    0.23918  -0.586 0.572368    
x1:x2     0.19239    0.04754   4.047 0.002900 ** 
x1:x3     0.40617    0.07174   5.662 0.000309 ***
x2:x3     0.25492    0.07174   3.553 0.006183 ** 
x1:x2:x3 -0.08509    0.02003  -4.247 0.002151 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.178 on 9 degrees of freedom
Multiple R-squared:  0.999, Adjusted R-squared:  0.9982 
F-statistic:  1249 on 7 and 9 DF,  p-value: 1.037e-12

MEPlot(funcao)

# Percebe-se por meio destas análises, quais são os níveis de ajustes dos fatores, visando o aumento da distância em [M]

# Otimização
library(rsm)
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(5,5,4),
      method = c("L-BFGS-B"))

$par
[1] 3 5 4

$value
[1] -5.916073

$counts
function gradient 
       4        4 

$convergence
[1] 0

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

library(NlcOptim)
solnl(X=x0,objfun = obj, lb = c(3,3,2),ub=c(5,5,4))

$par
     [,1]
[1,]    3
[2,]    5
[3,]    4

$fn
        1 
-5.916073 

$counts
     nfval ngval
[1,]    33     6

$lambda
$lambda$lower
          [,1]
[1,] 0.1028581
[2,] 0.0000000
[3,] 0.0000000

$lambda$upper
          [,1]
[1,] 0.0000000
[2,] 0.9130969
[3,] 1.0766090


$grad
           [,1]
[1,]  0.1028581
[2,] -0.9130969
[3,] -1.0766090

$hessian
             [,1]        [,2]        [,3]
[1,] 5.969496e-05   0.1480236  0.01930137
[2,] 1.480236e-01 367.0501881 47.86097263
[3,] 1.930137e-02  47.8609726  6.24159627

library(rsm)
persp(funcao,~x1+x2,zlab = "comprimento",col="blue",
      contours = "colors",theta = 30,phi = 30)

persp(funcao,~x1+x3,zlab = "comprimento",col="blue",
      contours = "colors",theta = 30,phi = 30)

persp(funcao,~x2+x3,zlab = "comprimento",col="blue",
      contours = "colors",theta = 30,phi = 30)

contour(funcao,~x1+x2, image = T )
points(3,5,col="blue",lwd=4) # ponto de ótimo

contour(funcao,~x1+x3, col="purple" )
points(3,4,col="blue",lwd=4) # ponto de ótimo

contour(funcao,~x1+x2, col=rainbow(10) )
points(5,4,col="blue",lwd=4) # ponto de ótimo

#Validação Estatística
DOE

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

$fit
       fit      lwr      upr
1 5.916073 5.426218 6.405928

$se.fit
[1] 0.1233668

$df
[1] 9

$residual.scale
[1] 0.1779652

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

$fit
       fit      lwr      upr
1 5.916073 5.636998 6.195148

$se.fit
[1] 0.1233668

$df
[1] 9

$residual.scale
[1] 0.1779652

# Com aleatorização ou randomização

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

DOE<-rbind(DOE,DOE)
DOE

DOE$Bloco<- as.factor(c(rep(1,8),rep(2,8)))

DOE

set.seed(30);DOE<- DOE[sample(nrow(DOE)),] # Aleatorização

DOE

DOE$y<-(c( 1.83,4.40, 1.83,5.50, 3.66,5.85,3.66,3.63,
          3.44,1.72, 1.94,4.40,5.78,4.20,5.95,3.74))
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 
-0.27063 -0.05281  0.00000  0.05281  0.27063 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) -6.00938    2.66159  -2.258   0.0585 .
x1           0.62500    0.64542   0.968   0.3651  
x2           1.65500    0.64542   2.564   0.0373 *
x3           2.42188    0.84152   2.878   0.0237 *
Bloco2       0.11875    0.09900   1.199   0.2694  
x1:x2       -0.17000    0.15654  -1.086   0.3135  
x1:x3       -0.27688    0.20410  -1.357   0.2170  
x2:x3       -0.33813    0.20410  -1.657   0.1416  
x1:x2:x3     0.07313    0.04950   1.477   0.1831  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.198 on 7 degrees of freedom
Multiple R-squared:  0.9915,    Adjusted R-squared:  0.9818 
F-statistic: 101.9 on 8 and 7 DF,  p-value: 1.497e-06

summary(aov(funcao))

            Df Sum Sq Mean Sq F value   Pr(>F)    
x1           1  0.001   0.001   0.027    0.874    
x2           1 11.239  11.239 286.670 6.14e-07 ***
x3           1 20.498  20.498 522.831 7.75e-08 ***
Bloco        1  0.056   0.056   1.439    0.269    
x1:x2        1  0.039   0.039   0.995    0.352    
x1:x3        1  0.004   0.004   0.100    0.761    
x2:x3        1  0.033   0.033   0.850    0.387    
x1:x2:x3     1  0.086   0.086   2.182    0.183    
Residuals    7  0.274   0.039                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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

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

Residuals:
     Min       1Q   Median       3Q      Max 
-0.39938 -0.06219  0.04750  0.10813  0.26063 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -2.90250    0.24842  -11.68 2.88e-08 ***
x2           0.83812    0.04872   17.20 2.52e-10 ***
x3           1.13188    0.04872   23.23 5.69e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1949 on 13 degrees of freedom
Multiple R-squared:  0.9847,    Adjusted R-squared:  0.9823 
F-statistic: 417.9 on 2 and 13 DF,  p-value: 1.598e-12

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.39176 -0.11794  0.05706  0.12353  0.26824 

Coefficients:
         Estimate Std. Error t value Pr(>|t|)   
x1       -0.77500    0.19128  -4.052  0.00288 **
x2        0.25500    0.19128   1.333  0.21524   
x3        0.63687    0.32517   1.959  0.08183 . 
x1:x2     0.15941    0.06464   2.466  0.03579 * 
x1:x3     0.14313    0.09753   1.467  0.17631   
x2:x3     0.08188    0.09753   0.839  0.42295   
x1:x2:x3 -0.02570    0.02724  -0.944  0.37006   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2419 on 9 degrees of freedom
Multiple R-squared:  0.998, Adjusted R-squared:  0.9965 
F-statistic: 654.8 on 7 and 9 DF,  p-value: 1.884e-11

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

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

Residuals:
     Min       1Q   Median       3Q      Max 
-0.55909 -0.28585  0.08425  0.17886  0.48036 

Coefficients:
      Estimate Std. Error t value Pr(>|t|)   
x1    -0.08162    0.08348  -0.978   0.3475   
x2     0.24005    0.10324   2.325   0.0384 * 
x3     0.37732    0.15165   2.488   0.0285 * 
x2:x3  0.17486    0.04149   4.214   0.0012 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3493 on 12 degrees of freedom
Multiple R-squared:  0.9946,    Adjusted R-squared:  0.9927 
F-statistic: 547.9 on 4 and 12 DF,  p-value: 1.816e-13

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

shapiro.test(funcao$residuals)

    Shapiro-Wilk normality test

data:  funcao$residuals
W = 0.95075, p-value = 0.5017

library(lmtest)
library(car)
residualPlot(funcao)

bptest(funcao)

    studentized Breusch-Pagan test

data:  funcao
BP = 1.2482, df = 3, p-value = 0.7415

durbinWatsonTest(funcao)

 lag Autocorrelation D-W Statistic p-value
   1      -0.2059954      2.334132   0.414
 Alternative hypothesis: rho != 0

dwtest(funcao)

    Durbin-Watson test

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

dwt(funcao)

 lag Autocorrelation D-W Statistic p-value
   1      -0.2059954      2.334132   0.408
 Alternative hypothesis: rho != 0

library(FrF2)
MEPlot(funcao)

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(5,5,4),
      method = c("L-BFGS-B"))

$par
[1] 3 5 4

$value
[1] -5.961862

$counts
function gradient 
       5        5 

$convergence
[1] 0

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

library(NlcOptim)
solnl(X=x0,objfun = obj, lb = c(3,3,2),ub=c(5,5,4))

$par
     [,1]
[1,]    3
[2,]    5
[3,]    4

$fn
        1 
-5.961862 

$counts
     nfval ngval
[1,]    13     4

$lambda
$lambda$lower
           [,1]
[1,] 0.08161962
[2,] 0.00000000
[3,] 0.00000000

$lambda$upper
          [,1]
[1,] 0.0000000
[2,] 0.9394892
[3,] 1.2516129


$grad
            [,1]
[1,]  0.08161962
[2,] -0.93948925
[3,] -1.25161290

$hessian
              [,1]          [,2]          [,3]
[1,]  1.000000e+00  3.748032e-09 -3.058125e-09
[2,]  3.748032e-09  7.383624e-02 -6.914123e-02
[3,] -3.058125e-09 -6.914123e-02  1.087904e-01

DOE

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

$fit
       fit      lwr      upr
1 5.961862 5.559959 6.363764

$se.fit
[1] 0.1844595

$df
[1] 12

$residual.scale
[1] 0.3492838