Dissertação_DOE_ANE

library(rsm)
library(lmtest)
library(carData)
library(car)
library(nortest)
library(NlcOptim)
library(FrF2)
library(ggthemes)
library(rio)
library(ROI)
library(MASS)
library(nloptr)

library(readxl)

DOE<-read_excel("kaio2023.xlsx")

#======================= Função y1 ======================================
funcao<-rsm(log(y1)~SO(x1,x2,x3,x4),data = DOE)
summary(funcao)

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

               Estimate  Std. Error t value Pr(>|t|)  
(Intercept)  9.78955416  4.09681031  2.3896  0.03416 *
x1          -0.07834421  0.23726327 -0.3302  0.74694  
x2           0.01770484  0.05462103  0.3241  0.75141  
x3          -7.39651588  9.02236357 -0.8198  0.42831  
x4           0.72032059  0.36524762  1.9721  0.07209 .
x1:x2        0.00013723  0.00168297  0.0815  0.93636  
x1:x3        0.57778829  0.25244583  2.2888  0.04102 *
x1:x4       -0.01667476  0.01262229 -1.3211  0.21112  
x2:x3        0.03545788  0.06731889  0.5267  0.60799  
x2:x4       -0.00102694  0.00336594 -0.3051  0.76552  
x3:x4       -0.32267558  0.50489166 -0.6391  0.53477  
x1^2        -0.00351166  0.00546561 -0.6425  0.53263  
x2^2        -0.00038459  0.00038867 -0.9895  0.34194  
x3^2        -4.35477226  8.74498006 -0.4980  0.62750  
x4^2        -0.03785521  0.02186245 -1.7315  0.10896  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple R-squared:  0.5315,    Adjusted R-squared:  -0.01516 
F-statistic: 0.9723 on 14 and 12 DF,  p-value: 0.5255

Analysis of Variance Table

Response: log(y1)
                    Df  Sum Sq  Mean Sq F value Pr(>F)
FO(x1, x2, x3, x4)   4 0.10658 0.026646  0.6533 0.6356
TWI(x1, x2, x3, x4)  6 0.31688 0.052813  1.2949 0.3302
PQ(x1, x2, x3, x4)   4 0.13171 0.032928  0.8073 0.5438
Residuals           12 0.48944 0.040786               
Lack of fit         10 0.46256 0.046256  3.4424 0.2460
Pure error           2 0.02687 0.013437               

Stationary point of response surface:
        x1         x2         x3         x4 
18.4635538 42.4294158  0.4368292  3.0104023 

Eigenanalysis:
eigen() decomposition
$values
[1]  0.0223225928 -0.0004071608 -0.0385706135 -4.3798685396

$vectors
          [,1]          [,2]        [,3]         [,4]
x1  0.93854096 -0.0706245437 -0.33140536  0.065753555
x2  0.06837813  0.9974866972 -0.01812372  0.004031261
x3  0.07438735 -0.0008331845  0.01300111 -0.997144322
x4 -0.33004821  0.0056363534 -0.94322478 -0.036924525

funcao<-lm(log(y1)~ -1+ SO(x1,x2,x3,x4),data = DOE)
summary(funcao)

Call:
lm.default(formula = log(y1) ~ -1 + SO(x1, x2, x3, x4), data = DOE)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.28474 -0.14459 -0.02912  0.11370  0.28880 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)  
SO(x1, x2, x3, x4)x1     3.917e-01  1.549e-01   2.529   0.0252 *
SO(x1, x2, x3, x4)x2     1.013e-01  4.897e-02   2.068   0.0592 .
SO(x1, x2, x3, x4)x3     9.316e+00  6.653e+00   1.400   0.1848  
SO(x1, x2, x3, x4)x4     1.034e+00  3.979e-01   2.598   0.0221 *
SO(x1, x2, x3, x4)x1:x2 -1.201e-03  1.852e-03  -0.648   0.5281  
SO(x1, x2, x3, x4)x1:x3  3.102e-01  2.641e-01   1.175   0.2612  
SO(x1, x2, x3, x4)x1:x4 -2.169e-02  1.453e-02  -1.493   0.1592  
SO(x1, x2, x3, x4)x2:x3 -1.212e-02  7.506e-02  -0.161   0.8742  
SO(x1, x2, x3, x4)x2:x4 -1.919e-03  3.904e-03  -0.492   0.6313  
SO(x1, x2, x3, x4)x3:x4 -5.011e-01  5.828e-01  -0.860   0.4055  
SO(x1, x2, x3, x4)x1^2  -1.144e-02  5.069e-03  -2.257   0.0418 *
SO(x1, x2, x3, x4)x2^2  -7.998e-04  4.058e-04  -1.971   0.0704 .
SO(x1, x2, x3, x4)x3^2  -1.578e+01  8.547e+00  -1.846   0.0878 .
SO(x1, x2, x3, x4)x4^2  -5.619e-02  2.389e-02  -2.352   0.0351 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2357 on 13 degrees of freedom
Multiple R-squared:  0.9997,    Adjusted R-squared:  0.9993 
F-statistic:  2662 on 14 and 13 DF,  p-value: < 2.2e-16

funcao<- lm(log(y1)~ -1 + FO(x1,x2,x3,x4)+PQ(x1,x2,x3,x4),data = DOE)
summary(funcao)

Call:
lm.default(formula = log(y1) ~ -1 + FO(x1, x2, x3, x4) + PQ(x1, 
    x2, x3, x4), data = DOE)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.37389 -0.10701  0.04037  0.13657  0.41008 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)   
FO(x1, x2, x3, x4)x1    4.502e-01  1.247e-01   3.611  0.00186 **
FO(x1, x2, x3, x4)x2    7.607e-02  3.145e-02   2.419  0.02578 * 
FO(x1, x2, x3, x4)x3    1.457e+01  5.272e+00   2.763  0.01238 * 
FO(x1, x2, x3, x4)x4    3.772e-01  1.336e-01   2.824  0.01085 * 
PQ(x1, x2, x3, x4)x1^2 -1.292e-02  3.491e-03  -3.700  0.00152 **
PQ(x1, x2, x3, x4)x2^2 -8.771e-04  3.482e-04  -2.519  0.02089 * 
PQ(x1, x2, x3, x4)x3^2 -1.791e+01  6.609e+00  -2.709  0.01391 * 
PQ(x1, x2, x3, x4)x4^2 -5.961e-02  2.166e-02  -2.751  0.01270 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2276 on 19 degrees of freedom
Multiple R-squared:  0.9995,    Adjusted R-squared:  0.9993 
F-statistic:  4997 on 8 and 19 DF,  p-value: < 2.2e-16

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

#par(mfrow=c(1,2))
persp(funcao,~x1+x2,zlab="y1",
      col = rainbow(250),
      contours = ("colors"))

persp(funcao,~x1+x3,zlab="y1",
      col = rainbow(350),
      contours = ("colors"))

persp(funcao,~x1+x4,zlab="y1",
      col = rainbow(150),
      contours = ("colors"))

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

persp(funcao,~x2+x4,zlab="y1",
      col = rainbow(450),
      contours = ("colors"))

persp(funcao,~x3+x4, zlab="y1",
      col = rainbow(450),
      contours = ("colors"))

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

legend("topright", legend = c("Real","Estimada"),
       col=c("turquoise","blue"),lwd=5)

hist(funcao$residuals,lwd=3,col="yellow",
     main="Histograma dos Resíduos",xlab = "Resíduos")

#Normalidade
shapiro.test(funcao$residuals)

    Shapiro-Wilk normality test

data:  funcao$residuals
W = 0.96992, p-value = 0.5997

#Homocedasticidade
residualPlot(funcao)

bptest(funcao)

    studentized Breusch-Pagan test

data:  funcao
BP = 8.595, df = 7, p-value = 0.2831

#Autocorrelação Residual
dwt(funcao)

 lag Autocorrelation D-W Statistic p-value
   1       0.1328287      1.721409   0.412
 Alternative hypothesis: rho != 0

dwtest(funcao)

    Durbin-Watson test

data:  funcao
DW = 1.7214, p-value = 0.2182
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 = 0.62908, df = 1, p-value = 0.4277

#================== Grupo 2 ====================================================

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

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

               Estimate  Std. Error t value  Pr(>|t|)    
(Intercept)  4.5737e+00  1.9856e+00  2.3034 0.0399516 *  
x1          -2.8633e-02  1.1500e-01 -0.2490 0.8075807    
x2           1.9077e-02  2.6474e-02  0.7206 0.4849486    
x3          -1.4559e+01  4.3730e+00 -3.3293 0.0060063 ** 
x4           6.9160e-01  1.7703e-01  3.9067 0.0020846 ** 
x1:x2        4.9076e-04  8.1570e-04  0.6016 0.5586063    
x1:x3       -2.2292e-01  1.2236e-01 -1.8219 0.0934701 .  
x1:x4       -1.1146e-02  6.1178e-03 -1.8219 0.0934701 .  
x2:x3        3.2971e-02  3.2628e-02  1.0105 0.3321938    
x2:x4       -4.2569e-03  1.6314e-03 -2.6093 0.0228296 *  
x3:x4       -7.8576e-01  2.4471e-01 -3.2110 0.0074786 ** 
x1^2         2.9847e-03  2.6491e-03  1.1267 0.2818933    
x2^2        -3.0915e-04  1.8838e-04 -1.6411 0.1267062    
x3^2         2.1630e+01  4.2385e+00  5.1031 0.0002605 ***
x4^2        -4.8079e-03  1.0596e-02 -0.4537 0.6581133    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple R-squared:  0.9224,    Adjusted R-squared:  0.8319 
F-statistic: 10.19 on 14 and 12 DF,  p-value: 0.0001334

Analysis of Variance Table

Response: log(y2)
                    Df  Sum Sq  Mean Sq F value    Pr(>F)
FO(x1, x2, x3, x4)   4 0.73314 0.183285 19.1294  3.84e-05
TWI(x1, x2, x3, x4)  6 0.24088 0.040147  4.1902 0.0167034
PQ(x1, x2, x3, x4)   4 0.39323 0.098307 10.2603 0.0007566
Residuals           12 0.11498 0.009581                  
Lack of fit         10 0.09005 0.009005  0.7227 0.7052373
Pure error           2 0.02492 0.012461                  

Stationary point of response surface:
         x1          x2          x3          x4 
13.47558118 60.30306499  0.36139555  0.07541881 

Eigenanalysis:
eigen() decomposition
$values
[1] 21.6373488007  0.0058586264 -0.0003357832 -0.0153731802

$vectors
            [,1]         [,2]         [,3]       [,4]
x1 -0.0051464291  0.900534819 -0.194045000 0.38904639
x2  0.0007634794  0.166935096  0.980604468 0.10269844
x3  0.9998217568 -0.002779319 -0.002246847 0.01853865
x4 -0.0181489919 -0.401449860 -0.027502146 0.91528807

funcao2<-lm(log(y2)~ -1+SO(x1,x2,x3,x4),data = DOE)
summary(funcao2)

Call:
lm.default(formula = log(y2) ~ -1 + SO(x1, x2, x3, x4), data = DOE)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.19793 -0.04579  0.01175  0.04093  0.16298 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
SO(x1, x2, x3, x4)x1     0.1909667  0.0741922   2.574 0.023120 *  
SO(x1, x2, x3, x4)x2     0.0581172  0.0234643   2.477 0.027775 *  
SO(x1, x2, x3, x4)x3    -6.7507094  3.1874571  -2.118 0.054027 .  
SO(x1, x2, x3, x4)x4     0.8379983  0.1906325   4.396 0.000723 ***
SO(x1, x2, x3, x4)x1:x2 -0.0001344  0.0008875  -0.151 0.881930    
SO(x1, x2, x3, x4)x1:x3 -0.3479607  0.1265161  -2.750 0.016528 *  
SO(x1, x2, x3, x4)x1:x4 -0.0134906  0.0069602  -1.938 0.074617 .  
SO(x1, x2, x3, x4)x2:x3  0.0107418  0.0359611   0.299 0.769885    
SO(x1, x2, x3, x4)x2:x4 -0.0046737  0.0018707  -2.498 0.026669 *  
SO(x1, x2, x3, x4)x3:x4 -0.8691201  0.2792367  -3.112 0.008246 ** 
SO(x1, x2, x3, x4)x1^2  -0.0007201  0.0024286  -0.297 0.771513    
SO(x1, x2, x3, x4)x2^2  -0.0005031  0.0001944  -2.588 0.022515 *  
SO(x1, x2, x3, x4)x3^2  16.2923237  4.0948148   3.979 0.001574 ** 
SO(x1, x2, x3, x4)x4^2  -0.0133754  0.0114478  -1.168 0.263632    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1129 on 13 degrees of freedom
Multiple R-squared:  0.9988,    Adjusted R-squared:  0.9976 
F-statistic: 787.2 on 14 and 13 DF,  p-value: < 2.2e-16

funcao2<-lm(log(y2)~ - 1+ FO(x1,x2,x3,x4)+PQ(x2,x3)+TWI(x1,x4)+TWI(x1,x3)
+TWI(x3,x4)+TWI(x2,x4), data = DOE)
summary(funcao2)

Call:
lm.default(formula = log(y2) ~ -1 + FO(x1, x2, x3, x4) + PQ(x2, 
    x3) + TWI(x1, x4) + TWI(x1, x3) + TWI(x3, x4) + TWI(x2, x4), 
    data = DOE)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.162001 -0.046483 -0.005419  0.060789  0.163302 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
FO(x1, x2, x3, x4)x1  0.1786620  0.0372888   4.791 0.000170 ***
FO(x1, x2, x3, x4)x2  0.0567220  0.0159865   3.548 0.002472 ** 
FO(x1, x2, x3, x4)x3 -5.8962530  2.0774739  -2.838 0.011354 *  
FO(x1, x2, x3, x4)x4  0.7899106  0.1620835   4.873 0.000143 ***
PQ(x2, x3)x2^2       -0.0004614  0.0001713  -2.694 0.015373 *  
PQ(x2, x3)x3^2       16.9026002  3.5055312   4.822 0.000159 ***
TWI(x1, x4)          -0.0142844  0.0062949  -2.269 0.036563 *  
TWI(x1, x3)          -0.3902983  0.0938392  -4.159 0.000657 ***
TWI(x3, x4)          -0.8973451  0.2537110  -3.537 0.002533 ** 
TWI(x2, x4)          -0.0048148  0.0017122  -2.812 0.011998 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1043 on 17 degrees of freedom
Multiple R-squared:  0.9987,    Adjusted R-squared:  0.9979 
F-statistic:  1292 on 10 and 17 DF,  p-value: < 2.2e-16

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

persp(funcao2,~x1+x2,zlab="y2",
      col = rainbow(250),
      contours = ("colors"))

persp(funcao2,~x1+x3,zlab="y2",
      col = rainbow(250),
      contours = ("colors"))

persp(funcao2,~x2+x3,zlab="y2",
      col = rainbow(250),
      contours = ("colors"))

persp(funcao2,~x3+x4,zlab="y2",
      col = rainbow(250),
      contours = ("colors"))

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

legend("topright", legend = c("Real","Estimada"),
       col=c("turquoise","blue"),lwd=5)

hist(funcao2$residuals,lwd=3,col="yellow",main="Histograma dos Resíduos",
     xlab="Resíduos")

shapiro.test(funcao2$residuals)

    Shapiro-Wilk normality test

data:  funcao2$residuals
W = 0.98362, p-value = 0.9333

residualPlot(funcao2)

bptest(funcao2)

    studentized Breusch-Pagan test

data:  funcao2
BP = 7.0044, df = 9, p-value = 0.6367

dwt(funcao2)

 lag Autocorrelation D-W Statistic p-value
   1       0.1637541      1.650433    0.29
 Alternative hypothesis: rho != 0

dwtest(funcao2)

    Durbin-Watson test

data:  funcao2
DW = 1.6504, p-value = 0.153
alternative hypothesis: true autocorrelation is greater than 0

bgtest(funcao2)

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

data:  funcao2
LM test = 2.0418, df = 1, p-value = 0.153

#================== Grupo 3 ====================================================

funcao3<-lm(y3~ SO(x1,x2,x3,x4),data = DOE)
summary(funcao3)

Call:
lm.default(formula = y3 ~ SO(x1, x2, x3, x4), data = DOE)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.3750 -0.6042  0.1667  1.0833  1.8750 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)  
(Intercept)              43.708333  37.381855   1.169   0.2650  
SO(x1, x2, x3, x4)x1     -3.291667   2.164938  -1.520   0.1543  
SO(x1, x2, x3, x4)x2     -0.069444   0.498396  -0.139   0.8915  
SO(x1, x2, x3, x4)x3    -42.916667  82.325679  -0.521   0.6116  
SO(x1, x2, x3, x4)x4      6.125000   3.332747   1.838   0.0910 .
SO(x1, x2, x3, x4)x1:x2  -0.008333   0.015356  -0.543   0.5973  
SO(x1, x2, x3, x4)x1:x3   5.000000   2.303473   2.171   0.0507 .
SO(x1, x2, x3, x4)x1:x4  -0.250000   0.115174  -2.171   0.0507 .
SO(x1, x2, x3, x4)x2:x3  -0.333333   0.614260  -0.543   0.5973  
SO(x1, x2, x3, x4)x2:x4   0.058333   0.030713   1.899   0.0818 .
SO(x1, x2, x3, x4)x3:x4   1.250000   4.606947   0.271   0.7907  
SO(x1, x2, x3, x4)x1^2    0.072917   0.049872   1.462   0.1694  
SO(x1, x2, x3, x4)x2^2    0.001296   0.003546   0.366   0.7211  
SO(x1, x2, x3, x4)x3^2  -45.833333  79.794658  -0.574   0.5763  
SO(x1, x2, x3, x4)x4^2   -0.583333   0.199487  -2.924   0.0127 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.843 on 12 degrees of freedom
Multiple R-squared:  0.8267,    Adjusted R-squared:  0.6246 
F-statistic:  4.09 on 14 and 12 DF,  p-value: 0.009624

funcao3<-lm(log(y3)~ -1+SO(x1,x2,x3,x4),data = DOE)
summary(funcao3)

Call:
lm.default(formula = log(y3) ~ -1 + SO(x1, x2, x3, x4), data = DOE)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.198481 -0.088862  0.007507  0.076826  0.158937 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)   
SO(x1, x2, x3, x4)x1    -6.504e-03  1.010e-01  -0.064  0.94961   
SO(x1, x2, x3, x4)x2     2.238e-02  3.193e-02   0.701  0.49564   
SO(x1, x2, x3, x4)x3     5.935e+00  4.337e+00   1.369  0.19434   
SO(x1, x2, x3, x4)x4     6.439e-01  2.594e-01   2.482  0.02749 * 
SO(x1, x2, x3, x4)x1:x2 -1.067e-03  1.208e-03  -0.884  0.39298   
SO(x1, x2, x3, x4)x1:x3  2.143e-01  1.721e-01   1.245  0.23526   
SO(x1, x2, x3, x4)x1:x4 -2.172e-02  9.471e-03  -2.294  0.03911 * 
SO(x1, x2, x3, x4)x2:x3 -4.122e-02  4.893e-02  -0.842  0.41480   
SO(x1, x2, x3, x4)x2:x4  4.202e-03  2.545e-03   1.651  0.12272   
SO(x1, x2, x3, x4)x3:x4  1.882e-03  3.800e-01   0.005  0.99612   
SO(x1, x2, x3, x4)x1^2   1.419e-03  3.304e-03   0.430  0.67457   
SO(x1, x2, x3, x4)x2^2  -3.663e-05  2.645e-04  -0.138  0.89198   
SO(x1, x2, x3, x4)x3^2  -1.001e+01  5.572e+00  -1.796  0.09573 . 
SO(x1, x2, x3, x4)x4^2  -5.655e-02  1.558e-02  -3.631  0.00305 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1537 on 13 degrees of freedom
Multiple R-squared:  0.9984,    Adjusted R-squared:  0.9967 
F-statistic:   584 on 14 and 13 DF,  p-value: 8.027e-16

funcao3<-rsm(log(y3)~ FO(x1,x2,x4)+PQ(x4)+
TWI(x2,x4)+TWI(x1,x4),data = DOE)
summary(funcao3)

Near-stationary-ridge situation detected -- stationary point altered
 Change 'threshold' if this is not what you intend

Call:
rsm(formula = log(y3) ~ FO(x1, x2, x4) + PQ(x4) + TWI(x2, x4) + 
    TWI(x1, x4), data = DOE)

              Estimate Std. Error t value Pr(>|t|)   
(Intercept)  1.4374725  0.6379394  2.2533  0.03562 * 
x1           0.0751201  0.0289471  2.5951  0.01731 * 
x2          -0.0178812  0.0077192 -2.3164  0.03126 * 
x4           0.5603044  0.2130472  2.6300  0.01605 * 
x4^2        -0.0516056  0.0139505 -3.6992  0.00142 **
x2:x4        0.0045714  0.0024013  1.9037  0.07144 . 
x1:x4       -0.0196440  0.0090050 -2.1814  0.04126 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple R-squared:  0.719, Adjusted R-squared:  0.6347 
F-statistic:  8.53 on 6 and 20 DF,  p-value: 0.0001123

Analysis of Variance Table

Response: log(y3)
               Df  Sum Sq  Mean Sq F value    Pr(>F)
FO(x1, x2, x4)  3 0.60436 0.201453  9.7043 0.0003681
PQ(x4)          1 0.28407 0.284069 13.6840 0.0014199
TWI(x2, x4)     1 0.07523 0.075231  3.6240 0.0714424
TWI(x1, x4)     1 0.09879 0.098787  4.7587 0.0412565
Residuals      20 0.41518 0.020759                  
Lack of fit    12 0.33098 0.027582  2.6204 0.0894042
Pure error      8 0.08420 0.010526                  

Stationary point of response surface:
       x1        x2        x4 
 7.975413 -1.855968  3.828562 

Eigenanalysis:
eigen() decomposition
$values
[1]  0.001900641  0.000000000 -0.053506289

$vectors
         [,1]          [,2]       [,3]
x1  0.9571241  2.266549e-01  0.1803914
x2 -0.2227335  9.739751e-01 -0.0419791
x4 -0.1852115 -1.110223e-16  0.9826987

funcao3<-lm(log(y3)~ -1+FO(x1,x2,x3,x4)+
             PQ(x3,x4)+TWI(x1,x4)+TWI(x2,x4),data = DOE)
summary(funcao3)

Call:
lm.default(formula = log(y3) ~ -1 + FO(x1, x2, x3, x4) + PQ(x3, 
    x4) + TWI(x1, x4) + TWI(x2, x4), data = DOE)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.21501 -0.09701  0.01168  0.09727  0.17383 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
FO(x1, x2, x3, x4)x1  0.077709   0.024980   3.111 0.005755 ** 
FO(x1, x2, x3, x4)x2 -0.017421   0.007153  -2.435 0.024899 *  
FO(x1, x2, x3, x4)x3  6.959486   2.559484   2.719 0.013616 *  
FO(x1, x2, x3, x4)x4  0.607929   0.169624   3.584 0.001979 ** 
PQ(x3, x4)x3^2       -8.805172   3.264975  -2.697 0.014286 *  
PQ(x3, x4)x4^2       -0.056249   0.013438  -4.186 0.000501 ***
TWI(x1, x4)          -0.020396   0.007908  -2.579 0.018386 *  
TWI(x2, x4)           0.004438   0.002240   1.981 0.062287 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1404 on 19 degrees of freedom
Multiple R-squared:  0.9981,    Adjusted R-squared:  0.9972 
F-statistic:  1223 on 8 and 19 DF,  p-value: < 2.2e-16

funcao3<-lm((y3)~ -1+ FO(x1,x3,x4)+PQ(x4)+TWI(x1,x3)+TWI(x1,x4),
            data = DOE)
summary(funcao3)

Call:
lm.default(formula = (y3) ~ -1 + FO(x1, x3, x4) + PQ(x4) + TWI(x1, 
    x3) + TWI(x1, x4), data = DOE)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.8216 -0.8888  0.1061  1.2105  4.2496 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
FO(x1, x3, x4)x1   0.3448     0.1570   2.196 0.039463 *  
FO(x1, x3, x4)x3 -36.2300    17.1591  -2.111 0.046895 *  
FO(x1, x3, x4)x4  10.6816     2.3920   4.466 0.000213 ***
PQ(x4)            -0.6513     0.1929  -3.377 0.002848 ** 
TWI(x1, x3)        2.0156     0.9902   2.036 0.054605 .  
TWI(x1, x4)       -0.3060     0.1172  -2.611 0.016321 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.997 on 21 degrees of freedom
Multiple R-squared:  0.9863,    Adjusted R-squared:  0.9824 
F-statistic: 251.5 on 6 and 21 DF,  p-value: < 2.2e-16

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

persp(funcao3,~x1+x3,zlab="y3",
      col = rainbow(250),
      contours = ("colors"))

persp(funcao3,~x1+x4,zlab="y3",
      col = rainbow(250),
      contours = ("colors"))

persp(funcao3,~x3+x4,zlab="y3",
      col = rainbow(250),
      contours = ("colors"))

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

legend("topright", legend = c("Real","Estimada"),
       col=c("turquoise","blue"),lwd=5)

hist(funcao3$residuals,lwd=3,col="yellow",main="Histograma dos Resíduos",
     xlab="Resíduos")

shapiro.test(funcao3$residuals)

    Shapiro-Wilk normality test

data:  funcao3$residuals
W = 0.96008, p-value = 0.3712

residualPlot(funcao3)

bptest(funcao3)

    studentized Breusch-Pagan test

data:  funcao3
BP = 4.7593, df = 5, p-value = 0.446

dwt(funcao3)

 lag Autocorrelation D-W Statistic p-value
   1      0.01888434      1.862135    0.66
 Alternative hypothesis: rho != 0

dwtest(funcao3)

    Durbin-Watson test

data:  funcao3
DW = 1.8621, p-value = 0.3312
alternative hypothesis: true autocorrelation is greater than 0

bgtest(funcao3,order = 1)

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

data:  funcao3
LM test = 0.013641, df = 1, p-value = 0.907

#=================== Otimização========================================
obj=function(x){a<- predict(funcao, newdata = data.frame(x1=x[1],
                                                         x2=x[2],
                                                         x3=x[3],
                                                         x4=x[4]))
return(-a)}

obj2=function(x){a<- predict(funcao2, newdata = data.frame(x1=x[1],
                                                         x2=x[2],
                                                         x3=x[3],
                                                         x4=x[4]))
return(-a)}

obj3=function(x){a<- predict(funcao3, newdata = data.frame(x1=x[1],
                                                           x2=x[2],
                                                           x3=x[3],
                                                           x4=x[4]))
return(-a)}

x0<-c(14,30,.3,1)
optim(par = x0,fn=obj,lower = c(14,30,.3,1),
      upper = c(22,60,.5,5),method = c("L-BFGS-B"),control =T)

$par
[1] 17.4253753 43.3637777  0.4067494  3.1639912

$value
[1] -9.130589

$counts
function gradient 
      39       39 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

optim(par = x0,fn=obj2,lower = c(14,30,.3,1),
      upper = c(22,60,.5,5),method = c("L-BFGS-B"),control =T)

$par
[1] 14.00000 35.38084  0.30000  5.00000

$value
[1] -2.795558

$counts
function gradient 
      26       26 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

optim(par = x0,fn=obj3,lower = c(14,30,.3,1),
      upper = c(22,60,.5,5),method = c("L-BFGS-B"),control =T)

$par
[1] 14.000000 30.000000  0.300000  4.911891

$value
[1] -18.1378

$counts
function gradient 
       7        7 

$convergence
[1] 0

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

#============== Intervalo de predição===================================
# Valores ótimos
predict(funcao, newdata = data.frame(x1=16
,x2=44,x3=.333333333,x4=3.16),
        se.fit = T,interval = "pred")

$fit
       fit      lwr     upr
1 9.007509 8.494298 9.52072

$se.fit
[1] 0.09121297

$df
[1] 19

$residual.scale
[1] 0.2276041

predict(funcao2, newdata = data.frame(x1=16
                                      ,x2=44,x3=.333333333,x4=3.16),
        se.fit = T,interval = "pred")

$fit
       fit      lwr      upr
1 2.451433 2.215513 2.687353

$se.fit
[1] 0.04027101

$df
[1] 17

$residual.scale
[1] 0.1043169

predict(funcao3, newdata = data.frame(x1=16
                                      ,x2=44,x3=.333333333,x4=3.16),
        se.fit = T,interval = "pred")

$fit
      fit      lwr      upr
1 15.9716 11.57548 20.36773

$se.fit
[1] 0.6932738

$df
[1] 21

$residual.scale
[1] 1.996997

library(nloptr)

obj=function(x){a<- predict(funcao, newdata = data.frame(x1=x[1],
                                                         x2=x[2],
                                                         x3=x[3],
                                                         x4=x[4]))

return(-a)}

constraint <- function(x) {
  return(as.integer(x)-x[1])}

#constraint2 <- function(x) {
#  return(as.integer(x)-x[2])}
#integer_constraints <- c("x1", "x2")

x0<-c(16,30,.3,1)

nloptr(x0,eval_f =obj,lb=c(14,30,.3,1),ub=c(22,60,.5,5),
      # 
      eval_g_ineq = constraint,
             # xint = variable_types,
       #eval_jac_g_eq = constraint,
       opts = list("algorithm"="NLOPT_GN_ISRES",
                   xtol_rel=1.0e-100,maxeval = 100000,
                   xint = constraint))

Call:
nloptr(x0 = x0, eval_f = obj, lb = c(14, 30, 0.3, 1), ub = c(22, 
    60, 0.5, 5), eval_g_ineq = constraint, opts = list(algorithm = "NLOPT_GN_ISRES", 
    xtol_rel = 1e-100, maxeval = 1e+05, xint = constraint))


Minimization using NLopt version 2.7.1 

NLopt solver status: 5 ( NLOPT_MAXEVAL_REACHED: Optimization stopped because 
maxeval (above) was reached. )

Number of Iterations....: 100000 
Termination conditions:  xtol_rel: 1e-100   maxeval: 1e+05 
Number of inequality constraints:  4 
Number of equality constraints:    0 
Current value of objective function:  -8.72609953406255 
Current value of controls: 22 31 0.4067437 3.163782