CI_CATAPULTA_ANEIRSON

library(rsm)
library(lmtest)
library(carData)
library(car)
library(nortest)
library(NlcOptim)
library(FrF2)
library(ggthemes)
library(rio)
library(ROI)
library(MASS)
library(nloptr)
DOE<-bbd(k=~x1+x2+x3,
         block = F,
         n0=5,
         randomize = F,
         coding = list(x1~(M-3)/2,
                       x2~(Mart-3)/1,
                       x3~(fixa-3)/1))
DOE

   run.order std.order M Mart fixa
1          1         1 1    2    3
2          2         2 5    2    3
3          3         3 1    4    3
4          4         4 5    4    3
5          5         5 1    3    2
6          6         6 5    3    2
7          7         7 1    3    4
8          8         8 5    3    4
9          9         9 3    2    2
10        10        10 3    4    2
11        11        11 3    2    4
12        12        12 3    4    4
13        13        13 3    3    3
14        14        14 3    3    3
15        15        15 3    3    3
16        16        16 3    3    3
17        17        17 3    3    3

Data are stored in coded form using these coding formulas ...
x1 ~ (M - 3)/2
x2 ~ (Mart - 3)/1
x3 ~ (fixa - 3)/1

DOE$y1<-c(0.95, 1.35, 2.72, 4.40, 1.21, 1.76, 2.33, 4.82,
          0.59,3.25, 3.47, 6.61, 3.11, 3.33, 3.50,
          3.55, 3.28 )

DOE$y2<-c(1.31, 2.02, 2.86, 4.32, 1.52, 3.24, 3.13, 
          4.31, 1.05, 3.36, 3.21, 5.34, 3.36, 3.15, 
          3.64, 3.45, 3.44  )
DOE

   run.order std.order M Mart fixa   y1   y2
1          1         1 1    2    3 0.95 1.31
2          2         2 5    2    3 1.35 2.02
3          3         3 1    4    3 2.72 2.86
4          4         4 5    4    3 4.40 4.32
5          5         5 1    3    2 1.21 1.52
6          6         6 5    3    2 1.76 3.24
7          7         7 1    3    4 2.33 3.13
8          8         8 5    3    4 4.82 4.31
9          9         9 3    2    2 0.59 1.05
10        10        10 3    4    2 3.25 3.36
11        11        11 3    2    4 3.47 3.21
12        12        12 3    4    4 6.61 5.34
13        13        13 3    3    3 3.11 3.36
14        14        14 3    3    3 3.33 3.15
15        15        15 3    3    3 3.50 3.64
16        16        16 3    3    3 3.55 3.45
17        17        17 3    3    3 3.28 3.44

Data are stored in coded form using these coding formulas ...
x1 ~ (M - 3)/2
x2 ~ (Mart - 3)/1
x3 ~ (fixa - 3)/1

funcao<-rsm(y1~SO(x1,x2,x3),data = DOE)
summary(funcao)

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

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept)  3.35400    0.15956 21.0210 1.387e-07 ***
x1           0.64000    0.12614  5.0738 0.0014410 ** 
x2           1.32750    0.12614 10.5241 1.526e-05 ***
x3           1.30250    0.12614 10.3259 1.731e-05 ***
x1:x2        0.32000    0.17839  1.7938 0.1159236    
x1:x3        0.48500    0.17839  2.7188 0.0298173 *  
x2:x3        0.12000    0.17839  0.6727 0.5227215    
x1^2        -0.97450    0.17387 -5.6047 0.0008119 ***
x2^2        -0.02450    0.17387 -0.1409 0.8919109    
x3^2         0.15050    0.17387  0.8656 0.4153954    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple R-squared:  0.9761,    Adjusted R-squared:  0.9454 
F-statistic: 31.79 on 9 and 7 DF,  p-value: 7.345e-05

Analysis of Variance Table

Response: y1
                Df  Sum Sq Mean Sq F value    Pr(>F)
FO(x1, x2, x3)   3 30.9469 10.3156 81.0413 8.444e-06
TWI(x1, x2, x3)  3  1.4081  0.4694  3.6874  0.070431
PQ(x1, x2, x3)   3  4.0610  1.3537 10.6346  0.005334
Residuals        7  0.8910  0.1273                  
Lack of fit      3  0.7657  0.2552  8.1466  0.035303
Pure error       4  0.1253  0.0313                  

Stationary point of response surface:
        x1         x2         x3 
-0.1566316  8.1411610 -7.3205083 

Stationary point in original units:
        M      Mart      fixa 
 2.686737 11.141161 -4.320508 

Eigenanalysis:
eigen() decomposition
$values
[1]  0.23474831 -0.03809012 -1.04515818

$vectors
         [,1]        [,2]       [,3]
x1 -0.2285521 -0.06143486  0.9715913
x2 -0.3511963 -0.92560278 -0.1411405
x3 -0.9079786  0.37347725 -0.1899727

funcao<-rsm(y1~ FO(x1,x2,x3)+PQ(x1)+TWI(x1,x3),data = DOE)
summary(funcao)

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

Call:
rsm(formula = y1 ~ FO(x1, x2, x3) + PQ(x1) + TWI(x1, x3), data = DOE)

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept)  3.41000    0.12122 28.1305 1.342e-11 ***
x1           0.64000    0.12857  4.9777 0.0004170 ***
x2           1.32750    0.12857 10.3248 5.367e-07 ***
x3           1.30250    0.12857 10.1304 6.495e-07 ***
x1^2        -0.96750    0.17671 -5.4751 0.0001934 ***
x1:x3        0.48500    0.18183  2.6673 0.0218989 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple R-squared:  0.961, Adjusted R-squared:  0.9433 
F-statistic: 54.22 on 5 and 11 DF,  p-value: 2.23e-07

Analysis of Variance Table

Response: y1
               Df  Sum Sq Mean Sq F value    Pr(>F)
FO(x1, x2, x3)  3 30.9469 10.3156 78.0010  1.07e-07
PQ(x1)          1  3.9645  3.9645 29.9771 0.0001934
TWI(x1, x3)     1  0.9409  0.9409  7.1146 0.0218989
Residuals      11  1.4547  0.1322                  
Lack of fit     7  1.3294  0.1899  6.0619 0.0504518
Pure error      4  0.1253  0.0313                  

Stationary point of response surface:
        x1         x2         x3 
 -2.685567   0.000000 -12.034169 

Stationary point in original units:
        M      Mart      fixa 
-2.371134  3.000000 -9.034169 

Eigenanalysis:
eigen() decomposition
$values
[1]  0.05737874  0.00000000 -1.02487874

$vectors
         [,1] [,2]       [,3]
x1 -0.2302556    0  0.9731302
x2  0.0000000   -1  0.0000000
x3 -0.9731302    0 -0.2302556

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

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

persp(funcao,~x2+x3,zlab="Distância",
      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.94154, p-value = 0.3369

#Homocedasticidade
residualPlot(funcao)

bptest(funcao)

    studentized Breusch-Pagan test

data:  funcao
BP = 3.5118, df = 5, p-value = 0.6216

#Autocorrelação Residual
dwt(funcao)

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

 lag Autocorrelation D-W Statistic p-value
   1     -0.08429283      2.001873    0.91
 Alternative hypothesis: rho != 0

dwtest(funcao)

    Durbin-Watson test

data:  funcao
DW = 2.0019, p-value = 0.5338
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.17193, df = 1, p-value = 0.6784

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

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

Residuals:
    Min      1Q  Median      3Q     Max 
-0.2737 -0.1087  0.0175  0.1087  0.2737 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          3.40800    0.11977  28.455 1.70e-08 ***
SO(x1, x2, x3)x1     0.63375    0.09469   6.693 0.000279 ***
SO(x1, x2, x3)x2     1.03625    0.09469  10.944 1.18e-05 ***
SO(x1, x2, x3)x3     0.85250    0.09469   9.004 4.25e-05 ***
SO(x1, x2, x3)x1:x2  0.18750    0.13390   1.400 0.204170    
SO(x1, x2, x3)x1:x3 -0.13500    0.13390  -1.008 0.346941    
SO(x1, x2, x3)x2:x3 -0.04500    0.13390  -0.336 0.746675    
SO(x1, x2, x3)x1^2  -0.48525    0.13051  -3.718 0.007476 ** 
SO(x1, x2, x3)x2^2  -0.29525    0.13051  -2.262 0.058142 .  
SO(x1, x2, x3)x3^2   0.12725    0.13051   0.975 0.362047    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2678 on 7 degrees of freedom
Multiple R-squared:  0.9746,    Adjusted R-squared:  0.942 
F-statistic: 29.89 on 9 and 7 DF,  p-value: 9.025e-05

funcao2<-rsm(y2~ FO(x1,x2,x3)+PQ(x1,x2),data = DOE)
summary(funcao2)

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

Call:
rsm(formula = y2 ~ FO(x1, x2, x3) + PQ(x1, x2), data = DOE)

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept)  3.46158    0.10661 32.4685 2.819e-12 ***
x1           0.63375    0.09486  6.6809 3.464e-05 ***
x2           1.03625    0.09486 10.9240 3.036e-07 ***
x3           0.85250    0.09486  8.9869 2.126e-06 ***
x1^2        -0.47855    0.13057 -3.6650  0.003722 ** 
x2^2        -0.28855    0.13057 -2.2099  0.049229 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple R-squared:   0.96, Adjusted R-squared:  0.9418 
F-statistic:  52.8 on 5 and 11 DF,  p-value: 2.561e-07

Analysis of Variance Table

Response: y2
               Df  Sum Sq Mean Sq F value    Pr(>F)
FO(x1, x2, x3)  3 17.6177  5.8726 81.5778 8.459e-08
PQ(x1, x2)      2  1.3876  0.6938  9.6375  0.003817
Residuals      11  0.7919  0.0720                  
Lack of fit     7  0.6664  0.0952  3.0347  0.150074
Pure error      4  0.1255  0.0314                  

Stationary point of response surface:
       x1        x2        x3 
0.6621529 1.7955996 0.0000000 

Stationary point in original units:
       M     Mart     fixa 
4.324306 4.795600 3.000000 

Eigenanalysis:
eigen() decomposition
$values
[1]  0.0000000 -0.2885526 -0.4785526

$vectors
   [,1] [,2] [,3]
x1    0    0    1
x2    0    1    0
x3    1    0    0

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

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

persp(funcao2,~x2+x3,zlab="Distância",
      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.94864, p-value = 0.4353

residualPlot(funcao2)

bptest(funcao2)

    studentized Breusch-Pagan test

data:  funcao2
BP = 3.4055, df = 5, p-value = 0.6377

dwt(funcao2)

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

 lag Autocorrelation D-W Statistic p-value
   1      0.03699807      1.822462   0.516
 Alternative hypothesis: rho != 0

dwtest(funcao2)

    Durbin-Watson test

data:  funcao2
DW = 1.8225, p-value = 0.2515
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 = 0.033273, df = 1, p-value = 0.8553

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

obj2=function(x){a<- predict(funcao2, 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 = -1,
      upper = 1,method = c("L-BFGS-B"),control =T)

$par
[1] 0.5813953 1.0000000 1.0000000

$value
[1] -6.367035

$counts
function gradient 
      10       10 

$convergence
[1] 0

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

ANE<- optim(par = x0,fn=obj2,lower = -1,
      upper = 1,method = c("L-BFGS-B"),control =T)
ANE

$par
[1] 0.6621529 1.0000000 1.0000000

$value
[1] -5.271596

$counts
function gradient 
       6        6 

$convergence
[1] 0

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

cat("Valor ótimo de x:", ANE$par, "\n")

Valor ótimo de x: 0.6621529 1 1 

cat("Valor ótimo da função objetivo:", ANE$value, "\n")

Valor ótimo da função objetivo: -5.271596 

library(nloptr)


constraint <- function(x) {
  return(as.integer(x)-x)}
x0<-c(1,1,1)

nloptr(x0,eval_f =obj,lb=c(-1,-1,-1),ub=c(1,1,1),
       eval_g_ineq = constraint,
             #eval_jac_g_eq = constraint,
       opts = list("algorithm"="NLOPT_GN_ISRES",
                   xtol_rel=1.0e-8))

Call:
nloptr(x0 = x0, eval_f = obj, lb = c(-1, -1, -1), ub = c(1, 1, 
    1), eval_g_ineq = constraint, opts = list(algorithm = "NLOPT_GN_ISRES", 
    xtol_rel = 1e-08))


Minimization using NLopt version 2.7.1 

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

Number of Iterations....: 100 
Termination conditions:  xtol_rel: 1e-08 
Number of inequality constraints:  3 
Number of equality constraints:    0 
Current value of objective function:  -6.1975 
Current value of controls: 1 1 1

cat("Valor ótimo de x:", ANE$solution, "\n")

Valor ótimo de x: 

cat("Valor ótimo da função objetivo:", ANE$objective, "\n")

Valor ótimo da função objetivo: 

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

$fit
     fit      lwr      upr
1 6.1975 5.138651 7.256349

$se.fit
[1] 0.3149405

$df
[1] 11

$residual.scale
[1] 0.3636619