R RSM_DOE_ANEIRSON

library(rsm)
library(lmtest)

Loading required package: zoo

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

library(car)

Loading required package: carData

library(FrF2)

Loading required package: DoE.base
Loading required package: grid
Loading required package: conf.design
Registered S3 method overwritten by 'DoE.base':
  method           from       
  factorize.factor conf.design

Attaching package: ‘DoE.base’

The following objects are masked from ‘package:stats’:

    aov, lm

The following object is masked from ‘package:graphics’:

    plot.design

The following object is masked from ‘package:base’:

    lengths

library(nortest)
library(ROI)

ROI: R Optimization Infrastructure
Registered solver plugins: nlminb, glpk.
Default solver: auto.

library(nloptr)
library(rio)

The following rio suggested packages are not installed: ‘arrow’, ‘hexView’, ‘fst’, ‘pzfx’, ‘rmatio’, ‘readODS’, ‘qs’
Use 'install_formats()' to install them

Attaching package: ‘rio’

The following object is masked from ‘package:conf.design’:

    factorize

DOE<-bbd(k=~ x1+x2+x3,
             randomize = F,
             block = F,
             n0=3,
         coding = list(x1~(A-0)/1,
                       x2~(T-1000)/100,
                        x3~(Carb-20)/10))
DOE$y<-c(2.49,3.46,4.08,4.40,2.87,
         3.99,4.62,4.92,1.83,1.16,
         2.39,2.25,2.07,2.64,2.82)
DOE$y2<-c(6.92,317,5.20,110.37,3.89,
          6.89,194.96,544.56,5.92,4.08,
          467.93,183.51,7.48,57.33,32.4050)
DOE

   run.order std.order  A    T Carb    y      y2
1          1         1 -1  900   20 2.49   6.920
2          2         2  1  900   20 3.46 317.000
3          3         3 -1 1100   20 4.08   5.200
4          4         4  1 1100   20 4.40 110.370
5          5         5 -1 1000   10 2.87   3.890
6          6         6  1 1000   10 3.99   6.890
7          7         7 -1 1000   30 4.62 194.960
8          8         8  1 1000   30 4.92 544.560
9          9         9  0  900   10 1.83   5.920
10        10        10  0 1100   10 1.16   4.080
11        11        11  0  900   30 2.39 467.930
12        12        12  0 1100   30 2.25 183.510
13        13        13  0 1000   20 2.07   7.480
14        14        14  0 1000   20 2.64  57.330
15        15        15  0 1000   20 2.82  32.405

Data are stored in coded form using these coding formulas ...
x1 ~ (A - 0)/1
x2 ~ (T - 1000)/100
x3 ~ (Carb - 20)/10

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

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

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept)  2.51000    0.34984  7.1748 0.0008182 ***
x1           0.33875    0.21423  1.5812 0.1746636    
x2           0.21500    0.21423  1.0036 0.3616406    
x3           0.54125    0.21423  2.5265 0.0527551 .  
x1:x2       -0.16250    0.30297 -0.5364 0.6147022    
x1:x3       -0.20500    0.30297 -0.6766 0.5286646    
x2:x3        0.13250    0.30297  0.4373 0.6801081    
x1^2         1.64500    0.31534  5.2166 0.0034194 ** 
x2^2        -0.54750    0.31534 -1.7362 0.1430375    
x3^2        -0.05500    0.31534 -0.1744 0.8683807    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple R-squared:  0.8954,    Adjusted R-squared:  0.7072 
F-statistic: 4.758 on 9 and 5 DF,  p-value: 0.05031

Analysis of Variance Table

Response: y
                Df  Sum Sq Mean Sq F value  Pr(>F)
FO(x1, x2, x3)   3  3.6314  1.2105  3.2969 0.11592
TWI(x1, x2, x3)  3  0.3440  0.1147  0.3123 0.81653
PQ(x1, x2, x3)   3 11.7454  3.9151 10.6634 0.01301
Residuals        5  1.8358  0.3672                
Lack of fit      3  1.5292  0.5097  3.3250 0.23975
Pure error       2  0.3066  0.1533                

Stationary point of response surface:
       x1        x2        x3 
0.2724875 0.8075658 5.3853867 

Stationary point in original units:
           A            T         Carb 
   0.2724875 1080.7565779   73.8538671 

Eigenanalysis:
eigen() decomposition
$values
[1]  1.65444822 -0.05362437 -0.55832385

$vectors
          [,1]        [,2]        [,3]
x1  0.99737086 -0.06559914 -0.03079145
x2 -0.03864649 -0.12205694 -0.99177041
x3 -0.06130098 -0.99035289  0.12427121

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

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

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

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept)  2.47615    0.23702 10.4472 2.483e-06 ***
x1           0.33875    0.17444  1.9419 0.0840478 .  
x2           0.21500    0.17444  1.2325 0.2489804    
x3           0.54125    0.17444  3.1028 0.0126652 *  
x1^2         1.64923    0.25601  6.4421 0.0001192 ***
x2^2        -0.54327    0.25601 -2.1221 0.0628241 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple R-squared:  0.8752,    Adjusted R-squared:  0.8059 
F-statistic: 12.62 on 5 and 9 DF,  p-value: 0.0007543

Analysis of Variance Table

Response: y
               Df  Sum Sq Mean Sq F value    Pr(>F)
FO(x1, x2, x3)  3  3.6314  1.2105  4.9725 0.0264474
PQ(x1, x2)      2 11.7342  5.8671 24.1016 0.0002431
Residuals       9  2.1909  0.2434                  
Lack of fit     7  1.8843  0.2692  1.7559 0.4100087
Pure error      2  0.3066  0.1533                  

Stationary point of response surface:
        x1         x2         x3 
-0.1026994  0.1978761  0.0000000 

Stationary point in original units:
           A            T         Carb 
  -0.1026994 1019.7876106   20.0000000 

Eigenanalysis:
eigen() decomposition
$values
[1]  1.6492308  0.0000000 -0.5432692

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

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

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

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

residualPlot(funcao)

hist(funcao$residuals)

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="red",lwd=3)

legend("topright",  legend = c("Real","Estimada"),
       col = c("turquoise","red"),lwd=2)

shapiro.test(funcao$residuals)

    Shapiro-Wilk normality test

data:  funcao$residuals
W = 0.92752, p-value = 0.2505

bptest(funcao)

    studentized Breusch-Pagan test

data:  funcao
BP = 9.7392, df = 5, p-value = 0.08297

dwtest(funcao)

    Durbin-Watson test

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

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.12841      2.070564   0.842
 Alternative hypothesis: rho != 0

bgtest(funcao,order = 1)

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

data:  funcao
LM test = 0.39314, df = 1, p-value = 0.5307

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

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

Residuals:
    Min      1Q  Median      3Q     Max 
-20.651  -2.494   1.907  18.157  57.330 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
SO(x1, x2, x3)x1       95.98      11.72   8.189 0.000179 ***
SO(x1, x2, x3)x2      -61.83      11.72  -5.275 0.001875 ** 
SO(x1, x2, x3)x3      171.27      11.72  14.612 6.45e-06 ***
SO(x1, x2, x3)x1:x2   -51.23      16.58  -3.090 0.021378 *  
SO(x1, x2, x3)x1:x3    86.65      16.58   5.227 0.001962 ** 
SO(x1, x2, x3)x2:x3   -70.64      16.58  -4.262 0.005311 ** 
SO(x1, x2, x3)x1^2     66.04      14.36   4.601 0.003690 ** 
SO(x1, x2, x3)x2^2     43.83      14.36   3.053 0.022424 *  
SO(x1, x2, x3)x3^2    121.53      14.36   8.466 0.000148 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 33.15 on 6 degrees of freedom
Multiple R-squared:  0.9906,    Adjusted R-squared:  0.9766 
F-statistic: 70.55 on 9 and 6 DF,  p-value: 2.146e-05

persp(funcao2,~x1+x2,col = rainbow(280),
      contours = "colors", zlab = "Super.Espec")

persp(funcao2,~x1+x3,col = rainbow(280),
      contours = "colors", zlab = "Super.Espec")

persp(funcao2,~x2+x3,col = rainbow(280),
      contours = "colors", zlab = "Super.Espec")

#===================== Resíduos ================================================
residualPlot(funcao2)

hist(funcao2$residuals)

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="red",lwd=3)

legend("topright",  legend = c("Real","Estimada"),
       col = c("turquoise","red"),lwd=2)

shapiro.test(funcao$residuals)

    Shapiro-Wilk normality test

data:  funcao$residuals
W = 0.92752, p-value = 0.2505

bptest(funcao)

    studentized Breusch-Pagan test

data:  funcao
BP = 9.7392, df = 5, p-value = 0.08297

dwtest(funcao)

    Durbin-Watson test

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

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.12841      2.070564   0.834
 Alternative hypothesis: rho != 0

bgtest(funcao,order = 1)

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

data:  funcao
LM test = 0.39314, df = 1, p-value = 0.5307

#================== Otimização =================================================

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)}

constraint <- function(x) {
  return(as.integer(x) - x[1])} # Defini que x1 deve ser inteiro

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

$par
[1] 1.0000000 0.1978761 1.0000000

$value
[1] -5.026656

$counts
function gradient 
       5        5 

$convergence
[1] 0

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

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

$par
[1]  1 -1  1

$value
[1] -769.0062

$counts
function gradient 
       2        2 

$convergence
[1] 0

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

ANE<- nloptr(x0,eval_f =  obj2,lb=c(-1,-1,-1),ub=c(1,1,1),
             eval_g_ineq = constraint, 
       opts = list("algorithm"="NLOPT_LN_COBYLA",
                   "xtol_rel"=1.0e-8))



ANE

Call:

nloptr(x0 = x0, eval_f = obj2, lb = c(-1, -1, -1), ub = c(1, 
    1, 1), eval_g_ineq = constraint, opts = list(algorithm = "NLOPT_LN_COBYLA", 
    xtol_rel = 1e-08))


Minimization using NLopt version 2.7.1 

NLopt solver status: 4 ( NLOPT_XTOL_REACHED: Optimization stopped 
because xtol_rel or xtol_abs (above) was reached. )

Number of Iterations....: 85 
Termination conditions:  xtol_rel: 1e-08 
Number of inequality constraints:  3 
Number of equality constraints:    0 
Optimal value of objective function:  -769.00625 
Optimal value of controls: 1 -1 1

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

Valor ótimo de x: 1 -1 1 

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

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

#================== Desirability ===============================================

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

if(a > 5){
  return(1)}

if(a<= 5 && a>= 4){
  return((a-4)/(5-4))}
return(-a)}

obj2=function(x){b<- predict(funcao2, newdata = data.frame(x1=x[1],
                                                           x2=x[2],
                                                           x3=x[3]))
if(b < 500 ){
  return(0)}

if(b > 580){
  return(1)}

if(b<=580 && b>= 500){
  return((b-500)/(580-500))}
return(-b)}


#cons_eq= function(x){g<- x[1]+x[2]+x[3] -1
#return(list(ceq=g, c=NULL))


#}
x0<-x0<-c(0,0,0) # ponto de partida
#Desirability function

d = function(x){
  
  -(obj(x)^0.5*obj2(x)^0.5)
}

d

function(x){
  
  -(obj(x)^0.5*obj2(x)^0.5)
}

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


x0<-c(1,    -0.228600183,   1)

optim(par = x0,fn=d,lower = -1,upper = 1,
      method = c("L-BFGS-B"))

$par
[1]  1.0000000 -0.2001506  1.0000000

$value
[1] -0.9698395

$counts
function gradient 
      61       61 

$convergence
[1] 52

$message
[1] "ERROR: ABNORMAL_TERMINATION_IN_LNSRCH"

nloptr(x0,eval_f =d,lb=c(-1,-1,-1),ub=c(1,1,1),
       eval_g_ineq = constraint, 
       opts = list("algorithm"="NLOPT_LN_COBYLA",
                   xtol_rel=1.0e-8))

Call:

nloptr(x0 = x0, eval_f = d, lb = c(-1, -1, -1), ub = c(1, 1, 
    1), eval_g_ineq = constraint, opts = list(algorithm = "NLOPT_LN_COBYLA", 
    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:  -0.969841241490465 
Current value of controls: 1 -0.2001425 0.9999997

# "algorithm"="NLOPT_LN_COBYLA"

#============== Intervalo de Confiança da Predição =============================
predict(funcao, newdata = data.frame(x1=1,x2=-0.200182988225692,x3=1),
        se.fit = T,interval = "pred")

$fit
       fit      lwr      upr
1 4.940575 3.602641 6.278508

$se.fit
[1] 0.3261448

$df
[1] 9

$residual.scale
[1] 0.493389

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

$fit
       fit     lwr     upr
1 580.0085 472.224 687.793

$se.fit
[1] 29.00341

$df
[1] 6

$residual.scale
[1] 33.15318