################################################################################
# Ejemplo del uso del paquete "rsm" para la
# MetodologÃa de superficie de respuesta
require(rsm)
## Loading required package: rsm
Etapa 1. Diseño \(2^{3-1} + 4 \times 0\).
################################################################################
# Generación de la configuración del primer experimento
exper0 <- cube(~ x1 + x2, x3 ~ x1 * x2, n0 = 4,
coding=c(x1 ~ (harina - 1) / 0.1,
x2 ~ (azucar - 0.5) / 0.1,
x3 ~ (mantequilla - 0.25) / 0.1), randomize=FALSE)
exper0
## run.order std.order harina azucar mantequilla
## 1 1 1 0.9 0.4 0.35
## 2 2 2 1.1 0.4 0.15
## 3 3 3 0.9 0.6 0.15
## 4 4 4 1.1 0.6 0.35
## 5 5 5 1.0 0.5 0.25
## 6 6 6 1.0 0.5 0.25
## 7 7 7 1.0 0.5 0.25
## 8 8 8 1.0 0.5 0.25
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (harina - 1)/0.1
## x2 ~ (azucar - 0.5)/0.1
## x3 ~ (mantequilla - 0.25)/0.1
as.data.frame(exper0)
## run.order std.order x1 x2 x3
## 1 1 1 -1 -1 1
## 2 2 2 1 -1 -1
## 3 3 3 -1 1 -1
## 4 4 4 1 1 1
## 5 5 5 0 0 0
## 6 6 6 0 0 0
## 7 7 7 0 0 0
## 8 8 8 0 0 0
set.seed(0)
exper1 <- cube(~ x1 + x2, x3 ~ x1*x2, n0 = 4,
coding=c(x1 ~ (harina - 1) / 0.1,
x2 ~ (azucar - 0.5) / 0.1,
x3 ~ (mantequilla - 0.25) / 0.1))
exper1
## run.order std.order harina azucar mantequilla
## 1 1 6 1.0 0.5 0.25
## 2 2 2 1.1 0.4 0.15
## 3 3 3 0.9 0.6 0.15
## 4 4 5 1.0 0.5 0.25
## 5 5 8 1.0 0.5 0.25
## 6 6 4 1.1 0.6 0.35
## 7 7 7 1.0 0.5 0.25
## 8 8 1 0.9 0.4 0.35
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (harina - 1)/0.1
## x2 ~ (azucar - 0.5)/0.1
## x3 ~ (mantequilla - 0.25)/0.1
# Gráfica de varianza y configuración gráfica
op <- par(no.readonly = TRUE)
par(mfrow = c(1, 2))
varfcn(exper1, ~ FO(x1, x2, x3))
varfcn(exper1, ~ FO(x1, x2, x3), contour = TRUE)
par(op)
# Intento de verificar si podrÃa ajustar un modelo de segundo orden (SO)
# varfcn(exper1,~SO(x1,x2,x3))
# Adicionar un diseño en estrella "rotable"
# djoin(exper1,star(n0=2,alpha="rotatable"))
set.seed(1)
exper2 <- djoin(exper1, star(n0 = 2, alpha = 1.5))
exper2
## run.order std.order harina azucar mantequilla Block
## 1 1 6 1.00 0.50 0.25 1
## 2 2 2 1.10 0.40 0.15 1
## 3 3 3 0.90 0.60 0.15 1
## 4 4 5 1.00 0.50 0.25 1
## 5 5 8 1.00 0.50 0.25 1
## 6 6 4 1.10 0.60 0.35 1
## 7 7 7 1.00 0.50 0.25 1
## 8 8 1 0.90 0.40 0.35 1
## 9 1 5 1.00 0.50 0.10 2
## 10 2 7 1.00 0.50 0.25 2
## 11 3 1 0.85 0.50 0.25 2
## 12 4 3 1.00 0.35 0.25 2
## 13 5 4 1.00 0.65 0.25 2
## 14 6 8 1.00 0.50 0.25 2
## 15 7 6 1.00 0.50 0.40 2
## 16 8 2 1.15 0.50 0.25 2
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (harina - 1)/0.1
## x2 ~ (azucar - 0.5)/0.1
## x3 ~ (mantequilla - 0.25)/0.1
# Generación del diseño en estrella con una longitud de 1.5
op <- par(no.readonly = TRUE)
par(mfrow = c(1, 2))
varfcn(exper2, ~ Block + SO(x1, x2, x3))
varfcn(exper2, ~ Block + SO(x1, x2, x3), contour = TRUE)
par(op)
write.csv2(code2val(exper1, codings = codings(exper1)),
file = "experimento1.csv",
row.names = FALSE)
exper1_r <- read.csv2("experimento1_realizado.csv")
exper1_r
## run.order std.order harina azucar mantequilla puntaje
## 1 1 6 1.0 0.5 0.25 24.7
## 2 2 2 1.1 0.4 0.15 28.9
## 3 3 3 0.9 0.6 0.15 20.2
## 4 4 5 1.0 0.5 0.25 25.0
## 5 5 8 1.0 0.5 0.25 24.7
## 6 6 4 1.1 0.6 0.35 27.1
## 7 7 7 1.0 0.5 0.25 25.5
## 8 8 1 0.9 0.4 0.35 21.5
exper1_r <- coded.data(exper1_r, x1 ~ (harina - 1) / 0.1,
x2 ~ (azucar - 0.5) / 0.1,
x3 ~ (mantequilla - 0.25) / 0.1)
exper1_r
## run.order std.order harina azucar mantequilla puntaje
## 1 1 6 1.0 0.5 0.25 24.7
## 2 2 2 1.1 0.4 0.15 28.9
## 3 3 3 0.9 0.6 0.15 20.2
## 4 4 5 1.0 0.5 0.25 25.0
## 5 5 8 1.0 0.5 0.25 24.7
## 6 6 4 1.1 0.6 0.35 27.1
## 7 7 7 1.0 0.5 0.25 25.5
## 8 8 1 0.9 0.4 0.35 21.5
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (harina - 1)/0.1
## x2 ~ (azucar - 0.5)/0.1
## x3 ~ (mantequilla - 0.25)/0.1
################################################################################
# Modelo 1
m1 <- rsm(puntaje ~ FO(x1, x2, x3), data = exper1_r)
summary(m1)
##
## Call:
## rsm(formula = puntaje ~ FO(x1, x2, x3), data = exper1_r)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 24.70000 0.17963 137.5077 1.678e-08 ***
## x1 3.57500 0.25403 14.0731 0.0001479 ***
## x2 -0.77500 0.25403 -3.0508 0.0379977 *
## x3 -0.12500 0.25403 -0.4921 0.6484543
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9811, Adjusted R-squared: 0.9669
## F-statistic: 69.2 on 3 and 4 DF, p-value: 0.0006658
##
## Analysis of Variance Table
##
## Response: puntaje
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3) 3 53.587 17.8625 69.2010 0.0006658
## Residuals 4 1.033 0.2581
## Lack of fit 1 0.605 0.6050 4.2456 0.1314343
## Pure error 3 0.428 0.1425
##
## Direction of steepest ascent (at radius 1):
## x1 x2 x3
## 0.97672947 -0.21173856 -0.03415138
##
## Corresponding increment in original units:
## harina azucar mantequilla
## 0.097672947 -0.021173856 -0.003415138
(dir1 <- steepest(m1))
## Path of steepest ascent from ridge analysis:
## dist x1 x2 x3 | harina azucar mantequilla | yhat
## 1 0.0 0.000 0.000 0.000 | 1.0000 0.5000 0.2500 | 24.700
## 2 0.5 0.488 -0.106 -0.017 | 1.0488 0.4894 0.2483 | 26.529
## 3 1.0 0.977 -0.212 -0.034 | 1.0977 0.4788 0.2466 | 28.361
## 4 1.5 1.465 -0.318 -0.051 | 1.1465 0.4682 0.2449 | 30.190
## 5 2.0 1.953 -0.423 -0.068 | 1.1953 0.4577 0.2432 | 32.018
## 6 2.5 2.442 -0.529 -0.085 | 1.2442 0.4471 0.2415 | 33.851
## 7 3.0 2.930 -0.635 -0.102 | 1.2930 0.4365 0.2398 | 35.680
## 8 3.5 3.419 -0.741 -0.120 | 1.3419 0.4259 0.2380 | 37.512
## 9 4.0 3.907 -0.847 -0.137 | 1.3907 0.4153 0.2363 | 39.341
## 10 4.5 4.395 -0.953 -0.154 | 1.4395 0.4047 0.2346 | 41.170
## 11 5.0 4.884 -1.059 -0.171 | 1.4884 0.3941 0.2329 | 43.002
Etapa 2. Camino del ascenso máximo.
set.seed(2)
exper2 <- dupe(dir1[2:9, ])
exper2
## run.order std.order dist x1 x2 x3 | harina azucar mantequilla
## 1 1 4 2.0 1.953 -0.423 -0.068 | 1.1953 0.4577 0.2432
## 2 2 1 0.5 0.488 -0.106 -0.017 | 1.0488 0.4894 0.2483
## 3 3 6 3.0 2.930 -0.635 -0.102 | 1.2930 0.4365 0.2398
## 4 4 3 1.5 1.465 -0.318 -0.051 | 1.1465 0.4682 0.2449
## 5 5 2 1.0 0.977 -0.212 -0.034 | 1.0977 0.4788 0.2466
## 6 6 5 2.5 2.442 -0.529 -0.085 | 1.2442 0.4471 0.2415
## 7 7 7 3.5 3.419 -0.741 -0.120 | 1.3419 0.4259 0.2380
## 8 8 8 4.0 3.907 -0.847 -0.137 | 1.3907 0.4153 0.2363
## |.1 yhat
## 1 | 32.018
## 2 | 26.529
## 3 | 35.680
## 4 | 30.190
## 5 | 28.361
## 6 | 33.851
## 7 | 37.512
## 8 | 39.341
write.csv2(exper2, file = "experimento2.csv",
row.names = FALSE)
exper2_r <- read.csv2("experimento2_realizado.csv")
exper2_r
## run.order std.order dist x1 x2 x3 X. harina azucar
## 1 1 4 2.0 1.953 -0.423 -0.068 | 1.1953 0.4577
## 2 2 1 0.5 0.488 -0.106 -0.017 | 1.0488 0.4894
## 3 3 6 3.0 2.930 -0.635 -0.102 | 1.2930 0.4365
## 4 4 3 1.5 1.465 -0.318 -0.051 | 1.1465 0.4682
## 5 5 2 1.0 0.977 -0.212 -0.034 | 1.0977 0.4788
## 6 6 5 2.5 2.442 -0.529 -0.085 | 1.2442 0.4471
## 7 7 7 3.5 3.419 -0.741 -0.120 | 1.3419 0.4259
## 8 8 8 4.0 3.907 -0.847 -0.137 | 1.3907 0.4153
## mantequilla X..1 yhat puntaje
## 1 0.2432 | 32.018 26.6
## 2 0.2483 | 26.529 24.8
## 3 0.2398 | 35.680 27.8
## 4 0.2449 | 30.190 27.5
## 5 0.2466 | 28.361 25.3
## 6 0.2415 | 33.851 26.0
## 7 0.2380 | 37.512 27.3
## 8 0.2363 | 39.341 24.3
exper2_r <- coded.data(exper2_r,
x1 ~ (harina - 1) / 0.1,
x2 ~ (azucar - 0.5) / 0.1,
x3 ~ (mantequilla - 0.25) / 0.1)
exper2_r
## run.order std.order dist harina azucar mantequilla X. harina azucar
## 1 1 4 2.0 1.1953 0.4577 0.2432 | 1.953 -0.423
## 2 2 1 0.5 1.0488 0.4894 0.2483 | 0.488 -0.106
## 3 3 6 3.0 1.2930 0.4365 0.2398 | 2.930 -0.635
## 4 4 3 1.5 1.1465 0.4682 0.2449 | 1.465 -0.318
## 5 5 2 1.0 1.0977 0.4788 0.2466 | 0.977 -0.212
## 6 6 5 2.5 1.2442 0.4471 0.2415 | 2.442 -0.529
## 7 7 7 3.5 1.3419 0.4259 0.2380 | 3.419 -0.741
## 8 8 8 4.0 1.3907 0.4153 0.2363 | 3.907 -0.847
## mantequilla X..1 yhat puntaje
## 1 -0.068 | 32.018 26.6
## 2 -0.017 | 26.529 24.8
## 3 -0.102 | 35.680 27.8
## 4 -0.051 | 30.190 27.5
## 5 -0.034 | 28.361 25.3
## 6 -0.085 | 33.851 26.0
## 7 -0.120 | 37.512 27.3
## 8 -0.137 | 39.341 24.3
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (harina - 1)/0.1
## x2 ~ (azucar - 0.5)/0.1
## x3 ~ (mantequilla - 0.25)/0.1
plot(puntaje ~ dist,data = exper2_r, pch = 19)
################################################################################
# Modelo 2
m2<-lm(puntaje ~ poly(dist, 2), data = exper2_r)
m2
##
## Call:
## lm(formula = puntaje ~ poly(dist, 2), data = exper2_r)
##
## Coefficients:
## (Intercept) poly(dist, 2)1 poly(dist, 2)2
## 26.2000 0.5246 -2.5151
plot(puntaje ~ dist, data = exper2_r, pch=19)
with(exper2_r,{
ord=order(dist)
lines(dist[ord], predict(m2)[ord])
})
Etapa 3. Diseño \(2^{3-1} + 4 \times 0\).
set.seed(4)
exper3 <- dupe(exper1)
codings(exper3) <- c(x1 ~ (harina - 1.25) / 0.1,
x2 ~ (azucar - 0.45) / 0.1,
x3 ~ (mantequilla - 0.25) / 0.1)
exper3
## run.order std.order harina azucar mantequilla
## 1 1 2 1.35 0.35 0.15
## 2 2 3 1.15 0.55 0.15
## 3 3 1 1.15 0.35 0.35
## 4 4 8 1.25 0.45 0.25
## 5 5 6 1.25 0.45 0.25
## 6 6 5 1.25 0.45 0.25
## 7 7 4 1.35 0.55 0.35
## 8 8 7 1.25 0.45 0.25
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (harina - 1.25)/0.1
## x2 ~ (azucar - 0.45)/0.1
## x3 ~ (mantequilla - 0.25)/0.1
write.csv2(code2val(exper3, codings = codings(exper1)),
file = "experimento3.csv",
row.names = FALSE)
exper3_r <- read.csv2("experimento3_realizado.csv")
codings(exper3_r) <- c(x1 ~ (harina - 1.25) / 0.1,
x2 ~ (azucar - 0.45) / 0.1,
x3 ~ (mantequilla - 0.25) / 0.1)
exper3_r
## run.order std.order harina azucar mantequilla puntaje
## 1 1 2 1.35 0.35 0.15 25.3
## 2 2 3 1.15 0.55 0.15 26.0
## 3 3 1 1.15 0.35 0.35 27.3
## 4 4 8 1.25 0.45 0.25 27.2
## 5 5 6 1.25 0.45 0.25 26.2
## 6 6 5 1.25 0.45 0.25 26.6
## 7 7 4 1.35 0.55 0.35 23.7
## 8 8 7 1.25 0.45 0.25 27.8
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (harina - 1.25)/0.1
## x2 ~ (azucar - 0.45)/0.1
## x3 ~ (mantequilla - 0.25)/0.1
################################################################################
# Modelo 3
m3 <- rsm(puntaje ~ FO(x1, x2, x3), data = exper3_r)
summary(m3)
##
## Call:
## rsm(formula = puntaje ~ FO(x1, x2, x3), data = exper3_r)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.26250 0.40509 64.8306 3.391e-07 ***
## x1 -1.07500 0.57289 -1.8765 0.1338
## x2 -0.72500 0.57289 -1.2655 0.2744
## x3 -0.07500 0.57289 -0.1309 0.9022
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.5624, Adjusted R-squared: 0.2341
## F-statistic: 1.713 on 3 and 4 DF, p-value: 0.3015
##
## Analysis of Variance Table
##
## Response: puntaje
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3) 3 6.7475 2.2492 1.7132 0.30145
## Residuals 4 5.2513 1.3128
## Lack of fit 1 3.7813 3.7813 7.7168 0.06911
## Pure error 3 1.4700 0.4900
##
## Direction of steepest ascent (at radius 1):
## x1 x2 x3
## -0.82768868 -0.55820864 -0.05774572
##
## Corresponding increment in original units:
## harina azucar mantequilla
## -0.082768868 -0.055820864 -0.005774572
Etapa 4. Doblar el experimento.
set.seed(5)
exper4 <- foldover(exper3, variable = "x1")
exper4
## run.order std.order harina azucar mantequilla
## 1 1 6 1.25 0.45 0.25
## 2 2 2 1.15 0.35 0.15
## 3 3 5 1.25 0.45 0.25
## 4 4 7 1.25 0.45 0.25
## 5 5 3 1.35 0.55 0.15
## 6 6 1 1.35 0.35 0.35
## 7 7 4 1.15 0.55 0.35
## 8 8 8 1.25 0.45 0.25
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (harina - 1.25)/0.1
## x2 ~ (azucar - 0.45)/0.1
## x3 ~ (mantequilla - 0.25)/0.1
write.csv2(code2val(exper4, codings = codings(exper3)),
file = "experimento4.csv",
row.names = FALSE)
exper4_r <- read.csv2("experimento4_realizado.csv")
codings(exper4_r) <- c(x1 ~ (harina - 1.25) / 0.1,
x2 ~ (azucar - 0.45) / 0.1,
x3 ~ (mantequilla - 0.25) / 0.1)
exper4_r
## run.order std.order harina azucar mantequilla puntaje
## 1 1 6 1.25 0.45 0.25 34.9
## 2 2 2 1.15 0.35 0.15 34.0
## 3 3 5 1.25 0.45 0.25 35.5
## 4 4 7 1.25 0.45 0.25 35.1
## 5 5 3 1.35 0.55 0.15 31.0
## 6 6 1 1.35 0.35 0.35 33.3
## 7 7 4 1.15 0.55 0.35 34.6
## 8 8 8 1.25 0.45 0.25 35.1
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (harina - 1.25)/0.1
## x2 ~ (azucar - 0.45)/0.1
## x3 ~ (mantequilla - 0.25)/0.1
union3_4 <- djoin(exper3_r, exper4_r)
union3_4
## run.order std.order harina azucar mantequilla puntaje Block
## 1 1 2 1.35 0.35 0.15 25.3 1
## 2 2 3 1.15 0.55 0.15 26.0 1
## 3 3 1 1.15 0.35 0.35 27.3 1
## 4 4 8 1.25 0.45 0.25 27.2 1
## 5 5 6 1.25 0.45 0.25 26.2 1
## 6 6 5 1.25 0.45 0.25 26.6 1
## 7 7 4 1.35 0.55 0.35 23.7 1
## 8 8 7 1.25 0.45 0.25 27.8 1
## 9 1 6 1.25 0.45 0.25 34.9 2
## 10 2 2 1.15 0.35 0.15 34.0 2
## 11 3 5 1.25 0.45 0.25 35.5 2
## 12 4 7 1.25 0.45 0.25 35.1 2
## 13 5 3 1.35 0.55 0.15 31.0 2
## 14 6 1 1.35 0.35 0.35 33.3 2
## 15 7 4 1.15 0.55 0.35 34.6 2
## 16 8 8 1.25 0.45 0.25 35.1 2
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (harina - 1.25)/0.1
## x2 ~ (azucar - 0.45)/0.1
## x3 ~ (mantequilla - 0.25)/0.1
################################################################################
# Modelo 4
m4<-rsm(puntaje~Block + FO(x1 ,x2, x3), data = union3_4)
summary(m4)
##
## Call:
## rsm(formula = puntaje ~ Block + FO(x1, x2, x3), data = union3_4)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.26250 0.40329 65.1208 1.388e-15 ***
## Block2 7.92500 0.57034 13.8953 2.543e-08 ***
## x1 -1.07500 0.40329 -2.6656 0.02197 *
## x2 -0.57500 0.40329 -1.4258 0.18169
## x3 0.32500 0.40329 0.8059 0.43739
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9486, Adjusted R-squared: 0.9299
## F-statistic: 50.72 on 4 and 11 DF, p-value: 5.075e-07
##
## Analysis of Variance Table
##
## Response: puntaje
## Df Sum Sq Mean Sq F value Pr(>F)
## Block 1 251.222 251.222 193.0793 2.543e-08
## FO(x1, x2, x3) 3 12.735 4.245 3.2625 0.063177
## Residuals 11 14.313 1.301
## Lack of fit 5 12.652 2.530 9.1464 0.008934
## Pure error 6 1.660 0.277
##
## Direction of steepest ascent (at radius 1):
## x1 x2 x3
## -0.8520282 -0.4557360 0.2575899
##
## Corresponding increment in original units:
## harina azucar mantequilla
## -0.08520282 -0.04557360 0.02575899
Etapa 5. Bloque con estrella.
set.seed(6)
exper5 <- star(exper4, n0=2, alpha = "orthogonal")
exper5
## run.order std.order harina azucar mantequilla
## 1 1 8 1.250000 0.4500000 0.2500000
## 2 2 3 1.250000 0.3085786 0.2500000
## 3 3 6 1.250000 0.4500000 0.3914214
## 4 4 5 1.250000 0.4500000 0.1085786
## 5 5 1 1.108579 0.4500000 0.2500000
## 6 6 4 1.250000 0.5914214 0.2500000
## 7 7 2 1.391421 0.4500000 0.2500000
## 8 8 7 1.250000 0.4500000 0.2500000
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (harina - 1.25)/0.1
## x2 ~ (azucar - 0.45)/0.1
## x3 ~ (mantequilla - 0.25)/0.1
op <- par(no.readonly = TRUE)
par(mfrow = c(1, 2))
union3_4_5 <- djoin(exper3,exper4,exper5)
varfcn(union3_4_5, ~Block + SO(x1, x2, x3))
varfcn(union3_4_5, ~Block + SO(x1, x2, x3), contour=T)
par(op)
write.csv2(code2val(exper5, codings = codings(exper5)),
file = "experimento5.csv",
row.names = FALSE)
exper5_r <- read.csv2("experimento5_realizado.csv")
codings(exper5_r) <- c(x1 ~ (harina - 1.25) / 0.1,
x2 ~ (azucar - 0.45) / 0.1,
x3 ~ (mantequilla - 0.25) / 0.1)
union3_4_5 <- djoin(exper3_r, exper4_r, exper5_r)
################################################################################
# Modelo 5
m5 <- rsm(puntaje ~ Block + SO(x1, x2, x3), data = union3_4_5)
summary(m5)
##
## Call:
## rsm(formula = puntaje ~ Block + SO(x1, x2, x3), data = union3_4_5)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.6996e+01 2.5051e-01 107.7647 < 2.2e-16 ***
## Block2 7.9250e+00 2.9641e-01 26.7366 4.601e-12 ***
## Block3 6.0000e-01 2.9641e-01 2.0242 0.065795 .
## x1 -1.0466e+00 1.7113e-01 -6.1160 5.208e-05 ***
## x2 -7.7224e-01 1.7113e-01 -4.5125 0.000711 ***
## x3 2.5202e-01 1.7113e-01 1.4727 0.166578
## x1:x2 -4.0000e-01 2.0959e-01 -1.9085 0.080537 .
## x1:x3 -1.5000e-01 2.0959e-01 -0.7157 0.487888
## x2:x3 -5.7175e-16 2.0959e-01 0.0000 1.000000
## x1^2 -1.2393e+00 1.9405e-01 -6.3865 3.471e-05 ***
## x2^2 -6.4286e-02 1.9405e-01 -0.3313 0.746139
## x3^2 -1.6429e-01 1.9405e-01 -0.8466 0.413766
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9881, Adjusted R-squared: 0.9772
## F-statistic: 90.77 on 11 and 12 DF, p-value: 8.46e-10
##
## Analysis of Variance Table
##
## Response: puntaje
## Df Sum Sq Mean Sq F value Pr(>F)
## Block 2 311.523 155.762 443.2137 5.678e-12
## FO(x1, x2, x3) 3 21.064 7.021 19.9791 5.849e-05
## TWI(x1, x2, x3) 3 1.460 0.487 1.3848 0.294838
## PQ(x1, x2, x3) 3 16.845 5.615 15.9771 0.000172
## Residuals 12 4.217 0.351
## Lack of fit 5 2.312 0.462 1.6993 0.252578
## Pure error 7 1.905 0.272
##
## Stationary point of response surface:
## x1 x2 x3
## 1.0644952 -9.3180901 0.2810583
##
## Stationary point in original units:
## harina azucar mantequilla
## 1.3564495 -0.4818090 0.2781058
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] -0.03002067 -0.16052168 -1.27731479
##
## $vectors
## [,1] [,2] [,3]
## x1 0.16811860 -0.04985422 0.98450530
## x2 -0.98128332 -0.10360829 0.16232180
## x3 -0.09391048 0.99336795 0.06633959
steepest(m5)
## Path of steepest ascent from ridge analysis:
## dist x1 x2 x3 | harina azucar mantequilla | yhat
## 1 0.0 0.000 0.000 0.000 | 1.2500 0.4500 0.2500 | 26.996
## 2 0.5 -0.227 -0.417 0.156 | 1.2273 0.4083 0.2656 | 27.484
## 3 1.0 -0.235 -0.922 0.307 | 1.2265 0.3578 0.2807 | 27.817
## 4 1.5 -0.189 -1.431 0.408 | 1.2311 0.3069 0.2908 | 28.102
## 5 2.0 -0.126 -1.939 0.473 | 1.2374 0.2561 0.2973 | 28.358
## 6 2.5 -0.055 -2.446 0.514 | 1.2445 0.2054 0.3014 | 28.591
## 7 3.0 0.020 -2.951 0.536 | 1.2520 0.1549 0.3036 | 28.804
## 8 3.5 0.098 -3.456 0.546 | 1.2598 0.1044 0.3046 | 28.999
## 9 4.0 0.178 -3.959 0.546 | 1.2678 0.0541 0.3046 | 29.177
## 10 4.5 0.258 -4.459 0.538 | 1.2758 0.0041 0.3038 | 29.337
## 11 5.0 0.339 -4.961 0.525 | 1.2839 -0.0461 0.3025 | 29.481
Etapa 6. Camino de máximo ascenso.
set.seed(7)
exper6 <- dupe(steepest(m5, dist = (2:9) / 3))
## Path of steepest ascent from ridge analysis:
exper6
## run.order std.order dist x1 x2 x3 | harina azucar
## 1 1 3 1.3333333 -0.207 -1.261 0.379 | 1.2293 0.3239
## 2 2 8 3.0000000 0.020 -2.951 0.536 | 1.2520 0.1549
## 3 3 2 1.0000000 -0.235 -0.922 0.307 | 1.2265 0.3578
## 4 4 7 2.6666667 -0.030 -2.615 0.523 | 1.2470 0.1885
## 5 5 5 2.0000000 -0.126 -1.939 0.473 | 1.2374 0.2561
## 6 6 4 1.6666667 -0.169 -1.600 0.433 | 1.2331 0.2900
## 7 7 6 2.3333333 -0.079 -2.278 0.502 | 1.2421 0.2222
## 8 8 1 0.6666667 -0.241 -0.584 0.212 | 1.2259 0.3916
## mantequilla |.1 yhat
## 1 0.2879 | 28.011
## 2 0.3036 | 28.804
## 3 0.2807 | 27.817
## 4 0.3023 | 28.664
## 5 0.2973 | 28.358
## 6 0.2933 | 28.190
## 7 0.3002 | 28.516
## 8 0.2712 | 27.603
write.csv2(exper6, file = "experimento6.csv",
row.names = FALSE)
exper6_r <- read.csv2("experimento6_realizado.csv")
codings(exper6_r) <- c(x1 ~ (harina - 1.25) / 0.1,
x2 ~ (azucar - 0.45) / 0.1,
x3 ~ (mantequilla - 0.25) / 0.1)
plot(puntaje ~ dist, data = exper6_r)
################################################################################
# Modelo 6
m6 <- lm(puntaje ~ poly(dist, 2), data = exper6_r)
plot(puntaje ~ dist, data = exper6_r)
with(exper6_r,{
ord <- order(dist)
lines(dist[ord], predict(m6)[ord])
})
Etapa 7. Diseño central compuesto. \(2^3\), en estrella \(+2\times0\).
set.seed(8)
exper7 <- ccd(~ x1 + x2 + x3, n0 = c(0, 2),
alpha="orthogonal",
coding=c(x1~(harina - 1.25) / 0.1,
x2~(azucar - 0.3) / 0.1,
x3~(mantequilla - 0.3) / 0.1))
exper7
## run.order std.order harina azucar mantequilla Block
## 1 1 7 1.15 0.4 0.4 1
## 2 2 2 1.35 0.2 0.2 1
## 3 3 6 1.35 0.2 0.4 1
## 4 4 1 1.15 0.2 0.2 1
## 5 5 3 1.15 0.4 0.2 1
## 6 6 8 1.35 0.4 0.4 1
## 7 7 5 1.15 0.2 0.4 1
## 8 8 4 1.35 0.4 0.2 1
## 9 1 4 1.25 0.5 0.3 2
## 10 2 5 1.25 0.3 0.1 2
## 11 3 6 1.25 0.3 0.5 2
## 12 4 7 1.25 0.3 0.3 2
## 13 5 1 1.05 0.3 0.3 2
## 14 6 8 1.25 0.3 0.3 2
## 15 7 3 1.25 0.1 0.3 2
## 16 8 2 1.45 0.3 0.3 2
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (harina - 1.25)/0.1
## x2 ~ (azucar - 0.3)/0.1
## x3 ~ (mantequilla - 0.3)/0.1
write.csv2(code2val(exper7, codings = codings(exper7)),
file = "experimento7.csv",
row.names = TRUE)
exper7_r <- read.csv2("experimento7_realizado.csv")
codings(exper7_r) <- c(x1 ~ (harina - 1.25) / 0.1,
x2 ~ (azucar - 0.3) / 0.1,
x3 ~ (mantequilla - 0.3) / 0.1)
################################################################################
# Modelo 7
m7 <- rsm(puntaje ~ Block + SO(x1, x2, x3), data = exper7_r)
summary(m7)
##
## Call:
## rsm(formula = puntaje ~ Block + SO(x1, x2, x3), data = exper7_r)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 27.20000 0.74128 36.6932 2.831e-07 ***
## Block 0.70000 0.37064 1.8886 0.117568
## x1 -0.90000 0.18532 -4.8564 0.004648 **
## x2 -0.05000 0.18532 -0.2698 0.798093
## x3 0.63750 0.18532 3.4400 0.018436 *
## x1:x2 -0.27500 0.26208 -1.0493 0.342094
## x1:x3 0.10000 0.26208 0.3816 0.718466
## x2:x3 -0.40000 0.26208 -1.5262 0.187476
## x1^2 -1.15000 0.18532 -6.2055 0.001587 **
## x2^2 -0.60000 0.18532 -3.2376 0.023010 *
## x3^2 -0.52500 0.18532 -2.8329 0.036549 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9421, Adjusted R-squared: 0.8262
## F-statistic: 8.131 on 10 and 5 DF, p-value: 0.01602
##
## Analysis of Variance Table
##
## Response: puntaje
## Df Sum Sq Mean Sq F value Pr(>F)
## Block 1 1.9600 1.9600 3.5669 0.117568
## FO(x1, x2, x3) 3 19.5025 6.5008 11.8305 0.010422
## TWI(x1, x2, x3) 3 1.9650 0.6550 1.1920 0.401731
## PQ(x1, x2, x3) 3 21.2550 7.0850 12.8935 0.008653
## Residuals 5 2.7475 0.5495
## Lack of fit 4 1.7675 0.4419 0.4509 0.789337
## Pure error 1 0.9800 0.9800
##
## Stationary point of response surface:
## x1 x2 x3
## -0.3421914 -0.1772769 0.6420873
##
## Stationary point in original units:
## harina azucar mantequilla
## 1.2157809 0.2822723 0.3642087
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] -0.3390336 -0.7534946 -1.1824718
##
## $vectors
## [,1] [,2] [,3]
## x1 0.1562180 -0.1664741 0.973592470
## x2 -0.6513590 0.7236263 0.228246488
## x3 0.7425142 0.6698144 -0.004609058
par(mfrow=c(1,3))
contour(m7, ~ x1 + x2 + x3, at=xs(m7), image=TRUE)
par(mfrow=c(1,1))
par(mfrow = c(2,2))
plot(m7)
par(mfrow = c(1,1))
ajustes <- predict(m7)
residuales <- resid(m7)
set.seed(9)
update(m7, ajustes + sample(residuales, replace = TRUE) ~ .)
##
## Call:
## rsm(formula = ajustes + sample(residuales, replace = TRUE) ~
## Block + FO(x1, x2, x3) + TWI(x1, x2, x3) + PQ(x1, x2, x3),
## data = exper7_r)
##
## Coefficients:
## (Intercept) Block FO(x1, x2, x3)x1
## 27.51250 0.54375 -0.92813
## FO(x1, x2, x3)x2 FO(x1, x2, x3)x3 TWI(x1, x2, x3)x1:x2
## -0.07656 0.56094 -0.49375
## TWI(x1, x2, x3)x1:x3 TWI(x1, x2, x3)x2:x3 PQ(x1, x2, x3)x1^2
## 0.31875 -0.23750 -1.15000
## PQ(x1, x2, x3)x2^2 PQ(x1, x2, x3)x3^2
## -0.51719 -0.60156
set.seed(9)
xs(update(m7, ajustes + sample(residuales, replace = TRUE) ~ .))
## x1 x2 x3
## -0.35437573 0.01010444 0.37035256
boot.crudo <- replicate(500,
xs(update(m7,
ajustes + sample(residuales,
replace = TRUE) ~ .)))
boot1 <- code2val(as.data.frame(t(boot.crudo)),
codings = codings(m7))
# Verificación del último modelo mediante boostrap
par(mfrow = c(1, 3))
plot(azucar ~ harina, data = boot1, col = "gray")
points(1.215, 0.282, col = "red", pch = 7)
plot(mantequilla~harina, data = boot1, col = "gray")
points(1.215, 0.364, col = "red", pch = 7)
plot(mantequilla~azucar, data = boot1, col = "gray")
points(0.282, 0.364, col = "red", pch = 7)
par(mfrow = c(1, 1))
par(mfrow = c(1, 3))
plot(azucar ~ harina, data = boot1,
col = "gray", xlim = c(1, 1.5), ylim = c(0, 0.5), asp = 1)
points(1.215, 0.282, col = "red", pch = 7)
plot(mantequilla ~ harina, data = boot1,
col = "gray", xlim = c(1, 1.5), ylim = c(0, 0.5), asp = 1)
points(1.215, 0.364, col = "red", pch = 7)
plot(mantequilla ~ azucar, data = boot1,
col = "gray", ylim = c(0, 0.5), xlim = c(0, 0.5), asp = 1)
points(0.282, 0.364, col = "red", pch = 7)
par(mfrow = c(1, 1))