MixOptim basic usage
This vignette demonstrates the basic usage of MixOptim functions to perform a basic optimization and plot.
Basic usage
The following code will perform a full range optimization of a mixture of 3 compounds based on two antagonist responses.
library(MixOptim)
#> Loading required package: ggplot2
#> Loading required package: patchwork
#> Loading required package: desirabilitysaPred<-function(x) 17.3359 * x[1] + 30.7765 * x[2] + 20.0501 * x[3] + 0.7506 * x[1] * x[2] + (-6.3443 * x[1] * x[3]) + (-5.9291 * x[2] * x[3]) +(-25.3093 * x[1] * x[2] * x[3])
pvPred<-function(x) 0.884640 * x[1] + 0.789863 * x[2] + 0.825016 * x[3] + (-0.108964 * x[1] * x[2]) + 0.107167 * x[1] * x[3] + (-0.005220 * x[2] * x[3]) + 1.625246 * x[1] * x[2] * x[3]
funcoes <- c(saPred, pvPred)
saD<-dMin(16.8293856, 31.41170555) # from experimental data
pvD<-dMax(0.7796004, 0.9019796) # from experimental data
overallD<-dOverall(saD, pvD)
a <- mixtureOptim(funcoes, overallD, 3, step = 0.01, plot = T)
#> Optimizing values.. 0 %
#> 10 %
#> 20 %
#> 30 %
#> 40 %
#> 50 %
#> 60 %
#> 70 %
#> 80 %
#> 90 %
#> 100 %
desirabilityPlot(funcoes, a$plotData, a$bestComposition, list(saD, pvD), c("min", "max"))
If the user want to get a more accurate composition, it is possible to do an optimization within a specific range (points +- alpha) based on previous values:
b <- mixtureFineOptim(funcoes, overallD, startPoint = a$bestComposition, step = 0.0001, alpha = 0.015)
#> Optimizing values.. 0 %
#> 10 %
#> 20 %
#> 30 %
#> 40 %
#> 50 %
#> 60 %
#> 70 %
#> 80 %
#> 90 %
#> 100 %
b
#> $maxDesirability
#> [1] 0.9603821
#>
#> $bestComposition
#> [1] 0.5630 0.0874 0.3496
desirabilityPlot(funcoes, a$plotData, b$bestComposition, list(saD, pvD), c("min", "max"))
Restricting seach intervals
# data
dados <- read.table(header = T, sep = "\t", text = "
ID TiO2 Vehicle Extender A Extender B Hiding Scrub
1 0.05 0.20 0.30 0.45 7.8953 533.67
2 0.45 0.20 0.30 0.05 32.862 749
3 0.05 0.60 0.30 0.05 3.721 39.5
4 0.05 0.20 0.70 0.05 9.2751 203.25
5 0.25 0.20 0.30 0.25 20.132 555.25
6 0.05 0.40 0.30 0.25 4.7137 51.75
7 0.05 0.20 0.50 0.25 8.3829 342.75
8 0.25 0.40 0.30 0.05 16.245 84.75
9 0.25 0.20 0.50 0.05 22.639 360.75
10 0.05 0.40 0.50 0.05 5.4645 48
11 0.05 0.33 0.43 0.18 5.8882 76
12 0.18 0.20 0.43 0.18 17.256 386.25
13 0.18 0.33 0.30 0.18 12.351 136
14 0.18 0.33 0.43 0.05 14.499 75.5
15 0.10 0.25 0.35 0.30 10.548 325.75
16 0.30 0.25 0.35 0.10 22.096 359
17 0.10 0.45 0.35 0.10 6.2888 40.75
18 0.10 0.25 0.55 0.10 10.629 136.67
19 0.15 0.30 0.40 0.15 11.777 114")
# function with coefficients generated with lm()
hiding<-function(x) 67.748*x[1] + 7.291*x[2] + 11.419*x[3] + 14.578*x[4] - 64.32*x[1]*x[2] + 35.878*x[1]*x[3] - 15.696*x[1]*x[4] - 31.006*x[2]*x[3] - 38.668*x[2]*x[4] - 6.59*x[3]*x[4]
scrub<-function(x) 3937.5*x[1] + 899.3*x[2] + 502*x[3] + 2354.8*x[4] - 8227.2*x[1]*x[2] - 3227.4*x[1]*x[3] - 2447.7*x[1]*x[4] - 2435.3*x[2]*x[3] - 6325.1*x[2]*x[4] - 1050.3*x[3]*x[4]
funcoes2 <- c(hiding, scrub)
des1<-dMax(min(dados$Hiding), max(dados$Hiding))
des2<-dMin(min(dados$Scrub), max(dados$Scrub))
finalD<-dOverall(des1, des2)
# optimize and generate data for plot
teste <- mixtureRangeOptim(funcoes2, finalD, midPoints = c(0.25, 0.4, 0.5, 0.25), alpha = c(0.2, 0.2, 0.2, 0.2), step = 0.01, plot=T)
#> Optimizing values.. 0 %..
#> 10 %..
#> 20 %..
#> 30 %..
#> 40 %..
#> 50 %..
#> 60 %..
#> 70 %..
#> 80 %..
#> 90 %..
#> 100 %..
# find a more accurate value based on previous composition, ingoring the last one and not generating plot data
teste2 <- mixtureRangeOptim(funcoes2, finalD, midPoints = teste$bestComposition, alpha = c(0.01, 0.01, 0.01, 0), step = 0.001, plot=F)
#> Optimizing values.. 0 %..
#> 10 %..
#> 20 %..
#> 30 %..
#> 40 %..
#> 50 %..
#> 60 %..
#> 70 %..
#> 80 %..
#> 90 %..
#> 100 %..
teste2
#> $maxDesirability
#> [1] 0.6474971
#>
#> $bestComposition
#> [1] 0.306 0.344 0.300 0.050
# plot result
desirabilityPlot(funcoes2, teste$plotData, teste2$bestComposition, list(des1, des2), c("max", "min"))
More specific desirability models: full example
# data
dados <- read.table(header = T, dec = ",", sep = "\t", text = "
x1 x2 x3 R1 R2 R3
1 0 0 0,76 8 5
1 0 0 0,75 8 5
0,5 0,5 0 1,4 7 7,5
0,5 0 0,5 0,55 8 10
0 1 0 4,1 4 10
0 1 0 4,4 4 10
0 0,5 0,5 0,9 7 12,5
0 0 1 0,42 9 15
0 0 1 0,4 10 15
0,6667 0,1667 0,1667 0,8 7 7,5
0,1667 0,6667 0,1667 1,7 7 10
0,1667 0,1667 0,6667 0,55 8 12,5
0,3333 0,3333 0,3333 0,8 8 10")
# create models, verify coefficients and significante, and generate functions
lm1 <- lm(data = dados, R1 ~ -1 + x1 + x2 + x3 + x1:x2 + x1:x3 + x2:x3)
summary(lm1)
#>
#> Call:
#> lm(formula = R1 ~ -1 + x1 + x2 + x3 + x1:x2 + x1:x3 + x2:x3,
#> data = dados)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -0.209034 -0.027362 -0.007836 0.068831 0.191668
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> x1 0.7678 0.1012 7.586 0.000128 ***
#> x2 4.2083 0.1012 41.580 1.21e-09 ***
#> x3 0.4274 0.1012 4.222 0.003924 **
#> x1:x2 -4.3273 0.5800 -7.461 0.000142 ***
#> x1:x3 0.3070 0.5800 0.529 0.612922
#> x2:x3 -5.6101 0.5800 -9.672 2.66e-05 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.1458 on 7 degrees of freedom
#> Multiple R-squared: 0.9967, Adjusted R-squared: 0.9939
#> F-statistic: 353.1 on 6 and 7 DF, p-value: 2.524e-08
flm1 <- function(x) 0.7678*x[1] + 4.2083*x[2] + 0.4274*x[3] - 4.3273*x[1]*x[2] + 0.3070*x[1]*x[3] - 5.6101*x[2]*x[3]
lm2 <- lm(data = dados, R2 ~ -1 + x1 + x2 + x3)
summary(lm2)
#>
#> Call:
#> lm(formula = R2 ~ -1 + x1 + x2 + x3, data = dados)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -0.67424 -0.57424 0.02576 0.62576 1.05836
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> x1 7.9742 0.3845 20.74 1.51e-09 ***
#> x2 4.5742 0.3845 11.90 3.17e-07 ***
#> x3 9.3742 0.3845 24.38 3.07e-10 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.656 on 10 degrees of freedom
#> Multiple R-squared: 0.9941, Adjusted R-squared: 0.9923
#> F-statistic: 561.3 on 3 and 10 DF, p-value: 1.936e-11
flm2 <- function(x) 7.9742*x[1] + 4.5742*x[2] + 9.3742*x[3]
lm3 <- lm(data = dados, R3 ~ -1 + x1 + x2 + x3)
summary(lm3)
#>
#> Call:
#> lm(formula = R3 ~ -1 + x1 + x2 + x3, data = dados)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -0.0008461 0.0001539 0.0001539 0.0001539 0.0011539
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> x1 5.000e+00 3.562e-04 14038 <2e-16 ***
#> x2 1.000e+01 3.562e-04 28076 <2e-16 ***
#> x3 1.500e+01 3.562e-04 42114 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.0006076 on 10 degrees of freedom
#> Multiple R-squared: 1, Adjusted R-squared: 1
#> F-statistic: 1.286e+09 on 3 and 10 DF, p-value: < 2.2e-16
flm3 <- function(x) 4.9998461*x[1] + 9.9998461*x[2] + 14.9998461*x[3]
funcoes2 <- c(flm1, flm2, flm3)
# specific range and target for Y1
des1<-dTarget(0.5, 0.6, 0.7)
des2<-dMax(8, max(dados$R2))
des3<-dMin(5, 10)
finalD<-dOverall(des1, des2, des3)
# full range optimization
teste <- mixtureOptim(funcoes2, finalD, 3, step = 0.01, plot = T)
#> Optimizing values.. 0 %
#> 10 %
#> 20 %
#> 30 %
#> 40 %
#> 50 %
#> 60 %
#> 70 %
#> 80 %
#> 90 %
#> 100 %
# plot
desirabilityPlot(funcoes2, teste$plotData, teste$bestComposition, list(des1, des2, des3), c("max", "max", "min"))
# more accurate optimization within +- 0.02 (default value for alpha)
teste2 <- mixtureFineOptim(funcoes2, finalD, teste$bestComposition, step = 0.0001)
#> Optimizing values.. 0 %
#> 10 %
#> 20 %
#> 30 %
#> 40 %
#> 50 %
#> 60 %
#> 70 %
#> 80 %
#> 90 %
#> 100 %
teste2
#> $maxDesirability
#> [1] 0.2581966
#>
#> $bestComposition
#> [1] 0.5657 0.0579 0.3764
# plot
desirabilityPlot(funcoes2, teste$plotData, teste2$bestComposition, list(des1, des2, des3), c("max", "max", "min"))