Se ajustan modelos cuadráticos de segundo orden para las tres respuestas del proceso de calcificación de cubos de tomate:
\[ Y = \beta_0 + \beta_1x_1 + \beta_2x_2 + \beta_3x_3 + \beta_{12}x_1x_2 + \beta_{13}x_1x_3 + \beta_{23}x_2x_3 + \beta_{11}x_1^2 + \beta_{22}x_2^2 + \beta_{33}x_3^2 + \varepsilon \]
Los factores codificados son:
Las metas de optimización son:
datos <- tibble::tribble(
~x1, ~x2, ~x3, ~Y1, ~Y2, ~Y3,
1, 1, 0, 3720.0, 23.04, 3.86,
1, -1, 0, 5578.8, 24.81, 3.75,
-1, 1, 0, 390.4, 15.40, 4.35,
-1, -1, 0, 248.1, 12.85, 4.20,
1, 0, 1, 7490.2, 26.48, 3.61,
1, 0, -1, 1842.9, 20.59, 4.00,
-1, 0, 1, 253.5, 14.81, 4.33,
-1, 0, -1, 152.0, 10.10, 4.35,
0, 1, 1, 2890.3, 20.79, 4.01,
0, 1, -1, 1162.5, 21.57, 4.33,
0, -1, 1, 1698.1, 22.85, 3.69,
0, -1, -1, 804.9, 17.75, 3.92,
0, 0, 0, 1505.9, 23.53, 3.85,
0, 0, 0, 1274.3, 20.00, 4.13,
0, 0, 0, 1660.3, 24.12, 3.77
) %>%
mutate(
conc_CaCl2 = 0.75 + 0.70*x1,
temperatura = 50 + 15*x2,
tiempo = 2.0 + 1.5*x3,
corrida = row_number()
)
DT::datatable(
datos,
rownames = FALSE,
caption = htmltools::tags$caption(
style = "caption-side: top; text-align: left; font-weight: bold;",
"Tabla 1. Diseño experimental codificado y variables reales del proceso."
),
options = list(pageLength = 15, scrollX = TRUE)
) %>%
DT::formatRound(columns = c("Y1", "Y2", "Y3", "conc_CaCl2", "temperatura", "tiempo"), digits = 3)formula_cuadratica <- as.formula(
y ~ x1 + x2 + x3 + x1:x2 + x1:x3 + x2:x3 + I(x1^2) + I(x2^2) + I(x3^2)
)
ajustar_modelo <- function(respuesta) {
f <- update(formula_cuadratica, as.formula(paste(respuesta, "~ .")))
lm(f, data = datos)
}
modelos <- list(
Y1 = ajustar_modelo("Y1"),
Y2 = ajustar_modelo("Y2"),
Y3 = ajustar_modelo("Y3")
)
nombres_respuesta <- c(
Y1 = "Contenido total de calcio, mg/g",
Y2 = "Firmeza, N/g",
Y3 = "pH"
)
metas_respuesta <- tibble::tribble(
~respuesta, ~nombre, ~objetivo, ~limite_inferior, ~limite_superior, ~meta_texto,
"Y1", "Contenido total de calcio, mg/g", "minimizar", 700, 800, "menor que 800; transición deseable 700-800",
"Y2", "Firmeza, N/g", "maximizar", 20.0, 20.5, "mayor que 20; transición deseable 20.0-20.5",
"Y3", "pH", "minimizar", 3.92, 3.95, "menor que 3.95; transición deseable 3.92-3.95"
)
coeficientes_modelos <- purrr::imap_dfr(
modelos,
~ broom::tidy(.x) %>%
mutate(respuesta = .y, .before = 1)
)
DT::datatable(
coeficientes_modelos,
rownames = FALSE,
caption = htmltools::tags$caption(
style = "caption-side: top; text-align: left; font-weight: bold;",
"Tabla 2. Coeficientes estimados de los modelos cuadráticos."
),
options = list(pageLength = 30, scrollX = TRUE)
) %>%
DT::formatRound(columns = c("estimate", "std.error", "statistic", "p.value"), digits = 5)El error estándar de cada modelo se reporta como \(s = \sqrt{MSE}\), es decir, el error residual estándar del modelo ajustado.
calidad_modelos <- purrr::imap_dfr(
modelos,
~ broom::glance(.x) %>%
transmute(
respuesta = .y,
R2 = r.squared,
R2_ajustado = adj.r.squared,
sigma_error_estandar = sigma,
F_global = statistic,
p_valor_global = p.value,
gl_residuales = df.residual,
AIC = AIC,
BIC = BIC
)
)
DT::datatable(
calidad_modelos,
rownames = FALSE,
caption = htmltools::tags$caption(
style = "caption-side: top; text-align: left; font-weight: bold;",
"Tabla 3. Calidad de ajuste y error estándar residual de cada modelo."
),
options = list(pageLength = 5, scrollX = TRUE)
) %>%
DT::formatRound(columns = c("R2", "R2_ajustado", "sigma_error_estandar", "F_global", "p_valor_global", "AIC", "BIC"), digits = 5)El test de Kolmogorov-Smirnov se aplica sobre residuales estandarizados. Con muestras pequeñas y parámetros estimados, debe interpretarse como diagnóstico complementario. En aplicaciones formales también suele revisarse el gráfico Q-Q y, opcionalmente, Shapiro-Wilk.
crear_grupos_fitted <- function(mod, k = 3) {
q <- quantile(fitted(mod), probs = seq(0, 1, length.out = k + 1), na.rm = TRUE)
q <- unique(q)
if (length(q) <= 2) {
return(factor(ifelse(fitted(mod) <= median(fitted(mod)), "bajo", "alto")))
}
cut(fitted(mod), breaks = q, include.lowest = TRUE)
}
ks_diagnostico <- function(mod) {
r <- as.numeric(scale(residuals(mod)))
out <- suppressWarnings(ks.test(r, "pnorm", exact = FALSE))
tibble(estadistico = unname(out$statistic), p_valor = out$p.value)
}
levene_diagnostico <- function(mod) {
grupo <- crear_grupos_fitted(mod, k = 3)
out <- car::leveneTest(residuals(mod), grupo, center = median)
tibble(
estadistico = out$`F value`[1],
p_valor = out$`Pr(>F)`[1]
)
}
dw_diagnostico <- function(mod) {
out <- lmtest::dwtest(mod)
tibble(
estadistico = unname(out$statistic),
p_valor = out$p.value
)
}
diagnosticos <- purrr::imap_dfr(
modelos,
function(mod, resp) {
bind_rows(
ks_diagnostico(mod) %>% mutate(respuesta = resp, supuesto = "Normalidad", prueba = "Kolmogorov-Smirnov"),
levene_diagnostico(mod) %>% mutate(respuesta = resp, supuesto = "Homocedasticidad", prueba = "Levene sobre grupos de valores ajustados"),
dw_diagnostico(mod) %>% mutate(respuesta = resp, supuesto = "Independencia", prueba = "Durbin-Watson")
) %>%
mutate(
decision_5pct = if_else(p_valor < 0.05, "Rechazar H0 al 5%", "No rechazar H0 al 5%")
) %>%
select(respuesta, supuesto, prueba, estadistico, p_valor, decision_5pct)
}
)
DT::datatable(
diagnosticos,
rownames = FALSE,
caption = htmltools::tags$caption(
style = "caption-side: top; text-align: left; font-weight: bold;",
"Tabla 4. Diagnósticos de normalidad, homocedasticidad e independencia de residuales."
),
options = list(pageLength = 10, scrollX = TRUE)
) %>%
DT::formatRound(columns = c("estadistico", "p_valor"), digits = 5)plot_diagnostico_residuales <- function(mod, resp) {
df <- tibble(
ajustado = fitted(mod),
residual = residuals(mod),
residual_estandar = rstandard(mod),
corrida = datos$corrida
)
p1 <- plot_ly(
df, x = ~ajustado, y = ~residual_estandar,
type = "scatter", mode = "markers+text",
text = ~corrida, textposition = "top center",
hovertemplate = paste(
"Corrida: %{text}<br>",
"Ajustado: %{x:.3f}<br>",
"Residual est.: %{y:.3f}<extra></extra>"
)
) %>%
layout(
title = paste(resp, "- residuales estandarizados vs ajustados"),
xaxis = list(title = "Valores ajustados"),
yaxis = list(title = "Residual estandarizado")
)
p2 <- plot_ly(
df, x = ~corrida, y = ~residual_estandar,
type = "scatter", mode = "lines+markers",
hovertemplate = paste(
"Corrida: %{x}<br>",
"Residual est.: %{y:.3f}<extra></extra>"
)
) %>%
layout(
title = paste(resp, "- residuales estandarizados en orden experimental"),
xaxis = list(title = "Corrida"),
yaxis = list(title = "Residual estandarizado")
)
subplot(p1, p2, nrows = 1, margin = 0.06, titleX = TRUE, titleY = TRUE)
}
plot_diagnostico_residuales(modelos$Y1, "Y1")plot_diagnostico_residuales(modelos$Y2, "Y2")plot_diagnostico_residuales(modelos$Y3, "Y3")anova_modelo <- function(mod, resp) {
as.data.frame(anova(mod)) %>%
tibble::rownames_to_column("fuente") %>%
as_tibble() %>%
mutate(respuesta = resp, .before = 1)
}
anova_todos <- purrr::imap_dfr(modelos, anova_modelo)
DT::datatable(
anova_todos,
rownames = FALSE,
caption = htmltools::tags$caption(
style = "caption-side: top; text-align: left; font-weight: bold;",
"Tabla 5. ANOVA secuencial de los modelos cuadráticos."
),
options = list(pageLength = 30, scrollX = TRUE)
) %>%
DT::formatRound(columns = c("Df", "Sum Sq", "Mean Sq", "F value", "Pr(>F)"), digits = 5)falta_ajuste <- function(mod, resp) {
yvar <- all.vars(formula(mod))[1]
df <- datos %>%
mutate(.y = .data[[yvar]])
sse_res <- sum(residuals(mod)^2)
df_res <- df.residual(mod)
puro <- df %>%
group_by(x1, x2, x3) %>%
summarise(
n = n(),
ss_puro = sum((.y - mean(.y))^2),
.groups = "drop"
)
sse_puro <- sum(puro$ss_puro)
df_puro <- sum(puro$n - 1)
sse_lof <- sse_res - sse_puro
df_lof <- df_res - df_puro
ms_lof <- sse_lof / df_lof
ms_puro <- sse_puro / df_puro
f_lof <- ms_lof / ms_puro
p_lof <- pf(f_lof, df_lof, df_puro, lower.tail = FALSE)
tibble(
respuesta = resp,
fuente = c("Falta de ajuste", "Error puro", "Error residual"),
gl = c(df_lof, df_puro, df_res),
SC = c(sse_lof, sse_puro, sse_res),
CM = c(ms_lof, ms_puro, sse_res / df_res),
F = c(f_lof, NA_real_, NA_real_),
p_valor = c(p_lof, NA_real_, NA_real_)
)
}
lof_todos <- purrr::imap_dfr(modelos, falta_ajuste)
DT::datatable(
lof_todos,
rownames = FALSE,
caption = htmltools::tags$caption(
style = "caption-side: top; text-align: left; font-weight: bold;",
"Tabla 6. Falta de ajuste contra error puro."
),
options = list(pageLength = 10, scrollX = TRUE)
) %>%
DT::formatRound(columns = c("SC", "CM", "F", "p_valor"), digits = 5)Para cada respuesta se optimiza el modelo dentro de la región experimental codificada \([-1,1]^3\). Para \(Y_1\) y \(Y_3\) se minimiza; para \(Y_2\) se maximiza.
decodificar <- function(x1, x2, x3) {
tibble(
conc_CaCl2 = 0.75 + 0.70*x1,
temperatura = 50 + 15*x2,
tiempo = 2.0 + 1.5*x3
)
}
predecir_respuestas <- function(x1, x2, x3) {
nd <- tibble(x1 = x1, x2 = x2, x3 = x3)
tibble(
Y1_pred = as.numeric(predict(modelos$Y1, newdata = nd)),
Y2_pred = as.numeric(predict(modelos$Y2, newdata = nd)),
Y3_pred = as.numeric(predict(modelos$Y3, newdata = nd))
)
}
optimizar_individual <- function(resp) {
mod <- modelos[[resp]]
objetivo <- metas_respuesta %>% filter(respuesta == resp) %>% pull(objetivo)
f_obj <- function(x) {
nd <- tibble(x1 = x[1], x2 = x[2], x3 = x[3])
yhat <- as.numeric(predict(mod, newdata = nd))
if (objetivo == "maximizar") -yhat else yhat
}
inicios <- expand.grid(
x1 = seq(-1, 1, length.out = 5),
x2 = seq(-1, 1, length.out = 5),
x3 = seq(-1, 1, length.out = 5)
)
soluciones <- purrr::pmap_dfr(
inicios,
function(x1, x2, x3) {
out <- optim(
par = c(x1, x2, x3),
fn = f_obj,
method = "L-BFGS-B",
lower = c(-1, -1, -1),
upper = c(1, 1, 1)
)
tibble(x1 = out$par[1], x2 = out$par[2], x3 = out$par[3], valor_objetivo = out$value)
}
) %>%
distinct(round(x1, 8), round(x2, 8), round(x3, 8), .keep_all = TRUE)
mejor <- if (objetivo == "maximizar") {
soluciones %>% mutate(valor_respuesta = -valor_objetivo) %>% arrange(desc(valor_respuesta)) %>% slice(1)
} else {
soluciones %>% mutate(valor_respuesta = valor_objetivo) %>% arrange(valor_respuesta) %>% slice(1)
}
bind_cols(
tibble(respuesta = resp, objetivo = objetivo),
mejor %>% select(x1, x2, x3, valor_respuesta),
decodificar(mejor$x1, mejor$x2, mejor$x3),
predecir_respuestas(mejor$x1, mejor$x2, mejor$x3)
)
}
optimos_individuales <- purrr::map_dfr(names(modelos), optimizar_individual) %>%
mutate(
nombre_respuesta = nombres_respuesta[respuesta],
cumple_meta_individual = case_when(
respuesta == "Y1" ~ if_else(valor_respuesta < 800, "Sí: Y1 < 800", "No"),
respuesta == "Y2" ~ if_else(valor_respuesta > 20, "Sí: Y2 > 20", "No"),
respuesta == "Y3" ~ if_else(valor_respuesta < 3.95, "Sí: Y3 < 3.95", "No")
)
) %>%
select(respuesta, nombre_respuesta, objetivo, x1, x2, x3, conc_CaCl2, temperatura, tiempo,
valor_respuesta, Y1_pred, Y2_pred, Y3_pred, cumple_meta_individual)
DT::datatable(
optimos_individuales,
rownames = FALSE,
caption = htmltools::tags$caption(
style = "caption-side: top; text-align: left; font-weight: bold;",
"Tabla 7. Óptimos individuales por variable de respuesta."
),
options = list(pageLength = 5, scrollX = TRUE)
) %>%
DT::formatRound(columns = c("x1", "x2", "x3", "conc_CaCl2", "temperatura", "tiempo", "valor_respuesta", "Y1_pred", "Y2_pred", "Y3_pred"), digits = 4)En cada gráfica se muestran la superficie de respuesta y el gráfico de contorno correspondiente. Para representar el óptimo individual en dos dimensiones, se grafican \(x_1\) y \(x_3\), dejando fijo \(x_2\) en el valor óptimo individual de cada respuesta.
generar_malla <- function(var_x = "x1", var_y = "x3", fijo = list(x2 = 0), n = 75) {
sec <- seq(-1, 1, length.out = n)
grid <- expand.grid(a = sec, b = sec)
names(grid) <- c(var_x, var_y)
for (nm in c("x1", "x2", "x3")) {
if (!nm %in% names(grid)) grid[[nm]] <- fijo[[nm]]
}
grid %>% select(x1, x2, x3)
}
plot_superficie_contorno <- function(resp, var_x = "x1", var_y = "x3") {
opt <- optimos_individuales %>% filter(respuesta == resp) %>% slice(1)
var_fija <- setdiff(c("x1", "x2", "x3"), c(var_x, var_y))
fijo <- as.list(opt[[var_fija]])
names(fijo) <- var_fija
grid <- generar_malla(var_x = var_x, var_y = var_y, fijo = fijo, n = 80) %>%
mutate(yhat = as.numeric(predict(modelos[[resp]], newdata = .)))
xs <- sort(unique(grid[[var_x]]))
ys <- sort(unique(grid[[var_y]]))
zmat <- matrix(grid$yhat, nrow = length(xs), ncol = length(ys))
opt_z <- as.numeric(predict(modelos[[resp]], newdata = opt %>% select(x1, x2, x3)))
p_surface <- plot_ly(x = xs, y = ys, z = zmat, type = "surface") %>%
add_markers(
x = opt[[var_x]], y = opt[[var_y]], z = opt_z,
marker = list(size = 5),
name = "Óptimo individual",
hovertemplate = paste0(
"Óptimo<br>", var_x, ": %{x:.3f}<br>", var_y, ": %{y:.3f}<br>", resp, ": %{z:.3f}<extra></extra>"
)
) %>%
layout(
title = paste0(resp, " - superficie"),
scene = list(
xaxis = list(title = var_x),
yaxis = list(title = var_y),
zaxis = list(title = paste0(resp, " predicho"))
)
)
p_contour <- plot_ly(x = xs, y = ys, z = zmat, type = "contour", contours = list(showlabels = TRUE)) %>%
add_markers(
x = opt[[var_x]], y = opt[[var_y]],
name = "Óptimo individual",
hovertemplate = paste0(
"Óptimo<br>", var_x, ": %{x:.3f}<br>", var_y, ": %{y:.3f}<extra></extra>"
)
) %>%
layout(
title = paste0(resp, " - contorno"),
xaxis = list(title = var_x),
yaxis = list(title = var_y)
)
subplot(p_surface, p_contour, nrows = 1, widths = c(0.56, 0.44), margin = 0.05, titleX = TRUE, titleY = TRUE) %>%
layout(title = paste0(nombres_respuesta[resp], "; ", var_fija, " fijo en ", round(opt[[var_fija]], 3)))
}
plot_superficie_contorno("Y1", var_x = "x1", var_y = "x3")plot_superficie_contorno("Y2", var_x = "x1", var_y = "x3")plot_superficie_contorno("Y3", var_x = "x1", var_y = "x3")El método gráfico se implementa como una búsqueda sobre la región experimental. Se identifica la región factible donde se cumplen simultáneamente:
\[ \hat{Y}_1 < 800, \qquad \hat{Y}_2 > 20, \qquad \hat{Y}_3 < 3.95. \]
Como el método gráfico normalmente entrega una región factible y no un único punto matemático, se reporta como solución representativa el punto factible más cercano al centro de los rangos óptimos \((750, 20.25, 3.935)\), usando distancia normalizada por el ancho de cada intervalo.
crear_malla3d <- function(n = 81) {
expand.grid(
x1 = seq(-1, 1, length.out = n),
x2 = seq(-1, 1, length.out = n),
x3 = seq(-1, 1, length.out = n)
) %>% as_tibble()
}
malla3d <- crear_malla3d(n = 81) %>%
mutate(
Y1_pred = as.numeric(predict(modelos$Y1, newdata = .)),
Y2_pred = as.numeric(predict(modelos$Y2, newdata = .)),
Y3_pred = as.numeric(predict(modelos$Y3, newdata = .)),
factible = Y1_pred < 800 & Y2_pred > 20 & Y3_pred < 3.95,
distancia_normalizada = ((Y1_pred - 750)/100)^2 + ((Y2_pred - 20.25)/0.5)^2 + ((Y3_pred - 3.935)/0.03)^2
)
solucion_grafica_grid <- malla3d %>%
filter(factible) %>%
arrange(distancia_normalizada) %>%
slice(1)
# Refinamiento continuo con optimización restringida por penalización suave.
objetivo_grafico <- function(x) {
pred <- predecir_respuestas(x[1], x[2], x[3])
dist <- ((pred$Y1_pred - 750)/100)^2 + ((pred$Y2_pred - 20.25)/0.5)^2 + ((pred$Y3_pred - 3.935)/0.03)^2
pen <- 0
pen <- pen + 1e5 * max(0, pred$Y1_pred - 800)^2
pen <- pen + 1e5 * max(0, 20 - pred$Y2_pred)^2
pen <- pen + 1e5 * max(0, pred$Y3_pred - 3.95)^2
dist + pen
}
out_grafico <- optim(
par = as.numeric(solucion_grafica_grid %>% select(x1, x2, x3)),
fn = objetivo_grafico,
method = "L-BFGS-B",
lower = c(-1, -1, -1),
upper = c(1, 1, 1)
)
solucion_grafica <- tibble(x1 = out_grafico$par[1], x2 = out_grafico$par[2], x3 = out_grafico$par[3]) %>%
bind_cols(decodificar(.$x1, .$x2, .$x3)) %>%
bind_cols(predecir_respuestas(.$x1, .$x2, .$x3)) %>%
mutate(
cumple_Y1 = if_else(Y1_pred < 800, "Sí", "No"),
cumple_Y2 = if_else(Y2_pred > 20, "Sí", "No"),
cumple_Y3 = if_else(Y3_pred < 3.95, "Sí", "No"),
interpretacion = paste0(
"Solución factible por superposición de contornos: Y1 = ", round(Y1_pred, 2),
" mg/g, Y2 = ", round(Y2_pred, 2),
" N/g y Y3 = ", round(Y3_pred, 3), "."
)
)
DT::datatable(
solucion_grafica,
rownames = FALSE,
caption = htmltools::tags$caption(
style = "caption-side: top; text-align: left; font-weight: bold;",
"Tabla 8. Solución simultánea obtenida por el método gráfico."
),
options = list(pageLength = 5, scrollX = TRUE)
) %>%
DT::formatRound(columns = c("x1", "x2", "x3", "conc_CaCl2", "temperatura", "tiempo", "Y1_pred", "Y2_pred", "Y3_pred"), digits = 4)Para visualizar la solución gráfica se muestran los contornos en el plano \(x_1\)-\(x_3\), fijando \(x_2\) en el valor de la solución gráfica.
plot_overlay_factible <- function(x2_fijo = solucion_grafica$x2[1]) {
sec <- seq(-1, 1, length.out = 120)
grid <- expand.grid(x1 = sec, x3 = sec) %>%
as_tibble() %>%
mutate(x2 = x2_fijo) %>%
mutate(
Y1_pred = as.numeric(predict(modelos$Y1, newdata = data.frame(x1 = x1, x2 = x2, x3 = x3))),
Y2_pred = as.numeric(predict(modelos$Y2, newdata = data.frame(x1 = x1, x2 = x2, x3 = x3))),
Y3_pred = as.numeric(predict(modelos$Y3, newdata = data.frame(x1 = x1, x2 = x2, x3 = x3))),
factible = Y1_pred < 800 & Y2_pred > 20 & Y3_pred < 3.95
)
z_factible <- matrix(as.numeric(grid$factible), nrow = length(sec), ncol = length(sec))
z_y1 <- matrix(grid$Y1_pred, nrow = length(sec), ncol = length(sec))
z_y2 <- matrix(grid$Y2_pred, nrow = length(sec), ncol = length(sec))
z_y3 <- matrix(grid$Y3_pred, nrow = length(sec), ncol = length(sec))
plot_ly() %>%
add_heatmap(x = sec, y = sec, z = z_factible, showscale = FALSE, name = "Región factible") %>%
add_contour(x = sec, y = sec, z = z_y1, contours = list(start = 800, end = 800, size = 1, showlabels = TRUE), name = "Y1 = 800") %>%
add_contour(x = sec, y = sec, z = z_y2, contours = list(start = 20, end = 20, size = 1, showlabels = TRUE), name = "Y2 = 20") %>%
add_contour(x = sec, y = sec, z = z_y3, contours = list(start = 3.95, end = 3.95, size = 1, showlabels = TRUE), name = "Y3 = 3.95") %>%
add_markers(
x = solucion_grafica$x1, y = solucion_grafica$x3,
name = "Solución gráfica",
hovertemplate = "Solución gráfica<br>x1: %{x:.3f}<br>x3: %{y:.3f}<extra></extra>"
) %>%
layout(
title = paste0("Región factible por superposición de contornos, x2 fijo = ", round(x2_fijo, 3)),
xaxis = list(title = "x1: concentración de CaCl2 codificada"),
yaxis = list(title = "x3: tiempo codificado")
)
}
plot_overlay_factible()Se emplean funciones de deseabilidad tipo Derringer-Suich. Para respuestas a minimizar:
\[ d(y)=\begin{cases} 1, & y \le L \\ \left(\frac{U-y}{U-L}\right)^s, & L < y < U \\ 0, & y \ge U \end{cases} \]
Para respuestas a maximizar:
\[ d(y)=\begin{cases} 0, & y \le L \\ \left(\frac{y-L}{U-L}\right)^s, & L < y < U \\ 1, & y \ge U \end{cases} \]
La deseabilidad global es:
\[ D = (d_1d_2d_3)^{1/3}. \]
d_menor <- function(y, L, U, s = 1) {
ifelse(y <= L, 1,
ifelse(y >= U, 0, ((U - y)/(U - L))^s))
}
d_mayor <- function(y, L, U, s = 1) {
ifelse(y <= L, 0,
ifelse(y >= U, 1, ((y - L)/(U - L))^s))
}
calcular_deseabilidad <- function(x1, x2, x3) {
pred <- predecir_respuestas(x1, x2, x3)
d1 <- d_menor(pred$Y1_pred, L = 700, U = 800, s = 1)
d2 <- d_mayor(pred$Y2_pred, L = 20.0, U = 20.5, s = 1)
d3 <- d_menor(pred$Y3_pred, L = 3.92, U = 3.95, s = 1)
D <- (d1*d2*d3)^(1/3)
tibble(
Y1_pred = pred$Y1_pred,
Y2_pred = pred$Y2_pred,
Y3_pred = pred$Y3_pred,
d1 = d1,
d2 = d2,
d3 = d3,
D_global = D
)
}
objetivo_deseabilidad <- function(x) {
-calcular_deseabilidad(x[1], x[2], x[3])$D_global
}
# Búsqueda inicial en malla para evitar quedar atrapado en zonas con D = 0.
malla_d <- crear_malla3d(n = 81) %>%
bind_cols(calcular_deseabilidad(.$x1, .$x2, .$x3)) %>%
arrange(desc(D_global))
inicios_d <- malla_d %>%
slice(1:30) %>%
select(x1, x2, x3)
soluciones_d <- purrr::pmap_dfr(
inicios_d,
function(x1, x2, x3) {
out <- optim(
par = c(x1, x2, x3),
fn = objetivo_deseabilidad,
method = "L-BFGS-B",
lower = c(-1, -1, -1),
upper = c(1, 1, 1)
)
tibble(x1 = out$par[1], x2 = out$par[2], x3 = out$par[3], D_global = -out$value)
}
) %>%
distinct(round(x1, 8), round(x2, 8), round(x3, 8), .keep_all = TRUE) %>%
arrange(desc(D_global))
opt_d <- soluciones_d %>% slice(1)
optimo_global <- opt_d %>%
select(x1, x2, x3) %>%
bind_cols(decodificar(.$x1, .$x2, .$x3)) %>%
bind_cols(calcular_deseabilidad(.$x1, .$x2, .$x3)) %>%
mutate(
cumple_Y1 = if_else(Y1_pred < 800, "Sí: Y1 < 800", "No"),
cumple_Y2 = if_else(Y2_pred > 20, "Sí: Y2 > 20", "No"),
cumple_Y3 = if_else(Y3_pred < 3.95, "Sí: Y3 < 3.95", "No"),
interpretacion = paste0(
"Con estos niveles se obtiene una deseabilidad global D = ", round(D_global, 4),
". Las respuestas predichas son: Y1 = ", round(Y1_pred, 2), " mg/g, ",
"Y2 = ", round(Y2_pred, 2), " N/g y Y3 = ", round(Y3_pred, 3), "."
)
)
DT::datatable(
optimo_global,
rownames = FALSE,
caption = htmltools::tags$caption(
style = "caption-side: top; text-align: left; font-weight: bold;",
"Tabla 9. Óptimo global mediante función de deseabilidad."
),
options = list(pageLength = 5, scrollX = TRUE)
) %>%
DT::formatRound(columns = c("x1", "x2", "x3", "conc_CaCl2", "temperatura", "tiempo", "Y1_pred", "Y2_pred", "Y3_pred", "d1", "d2", "d3", "D_global"), digits = 5)La superficie se muestra en el plano \(x_1\)-\(x_3\), fijando \(x_2\) en el óptimo global de deseabilidad.
plot_deseabilidad <- function(x2_fijo = optimo_global$x2[1]) {
sec <- seq(-1, 1, length.out = 100)
grid <- expand.grid(x1 = sec, x3 = sec) %>%
as_tibble() %>%
mutate(x2 = x2_fijo) %>%
bind_cols(calcular_deseabilidad(.$x1, .$x2, .$x3))
zD <- matrix(grid$D_global, nrow = length(sec), ncol = length(sec))
p_surface <- plot_ly(x = sec, y = sec, z = zD, type = "surface") %>%
add_markers(
x = optimo_global$x1, y = optimo_global$x3, z = optimo_global$D_global,
name = "Óptimo global",
hovertemplate = "Óptimo global<br>x1: %{x:.3f}<br>x3: %{y:.3f}<br>D: %{z:.4f}<extra></extra>"
) %>%
layout(
title = "Deseabilidad global - superficie",
scene = list(
xaxis = list(title = "x1"),
yaxis = list(title = "x3"),
zaxis = list(title = "D global")
)
)
p_contour <- plot_ly(x = sec, y = sec, z = zD, type = "contour", contours = list(showlabels = TRUE)) %>%
add_markers(
x = optimo_global$x1, y = optimo_global$x3,
name = "Óptimo global",
hovertemplate = "Óptimo global<br>x1: %{x:.3f}<br>x3: %{y:.3f}<extra></extra>"
) %>%
layout(
title = "Deseabilidad global - contorno",
xaxis = list(title = "x1"),
yaxis = list(title = "x3")
)
subplot(p_surface, p_contour, nrows = 1, widths = c(0.56, 0.44), margin = 0.05, titleX = TRUE, titleY = TRUE) %>%
layout(title = paste0("Deseabilidad global; x2 fijo = ", round(x2_fijo, 3)))
}
plot_deseabilidad()comparacion <- bind_rows(
solucion_grafica %>%
transmute(
metodo = "Gráfico / superposición de contornos",
x1, x2, x3, conc_CaCl2, temperatura, tiempo,
Y1_pred, Y2_pred, Y3_pred,
D_global = calcular_deseabilidad(x1, x2, x3)$D_global,
comentario = interpretacion
),
optimo_global %>%
transmute(
metodo = "Función de deseabilidad",
x1, x2, x3, conc_CaCl2, temperatura, tiempo,
Y1_pred, Y2_pred, Y3_pred,
D_global,
comentario = interpretacion
)
)
DT::datatable(
comparacion,
rownames = FALSE,
caption = htmltools::tags$caption(
style = "caption-side: top; text-align: left; font-weight: bold;",
"Tabla 10. Comparación de soluciones simultáneas."
),
options = list(pageLength = 5, scrollX = TRUE)
) %>%
DT::formatRound(columns = c("x1", "x2", "x3", "conc_CaCl2", "temperatura", "tiempo", "Y1_pred", "Y2_pred", "Y3_pred", "D_global"), digits = 5)Interpretación general. Los óptimos individuales pueden ubicarse en diferentes zonas de la región experimental; por eso no necesariamente son compatibles entre sí. La solución gráfica identifica una región factible por superposición de restricciones, mientras que la deseabilidad traduce cada meta a una escala común de 0 a 1 y selecciona el punto que maximiza el compromiso global. Si la deseabilidad global es baja, ello indica conflicto entre respuestas o rangos deseados estrechos; en ese caso conviene revisar los pesos, exponentes de deseabilidad o los límites tecnológicos usados para cada respuesta. También debe revisarse si existen predicciones físicamente imposibles, especialmente en modelos cuadráticos con óptimos en frontera; en tal caso la interpretación debe limitarse a la región experimental y validarse con corridas confirmatorias.
# Ejecutar manualmente si desea guardar tablas en archivos CSV.
write.csv(coeficientes_modelos, "coeficientes_modelos.csv", row.names = FALSE)
write.csv(calidad_modelos, "calidad_modelos.csv", row.names = FALSE)
write.csv(diagnosticos, "diagnosticos_supuestos.csv", row.names = FALSE)
write.csv(anova_todos, "anova_modelos.csv", row.names = FALSE)
write.csv(optimos_individuales, "optimos_individuales.csv", row.names = FALSE)
write.csv(optimo_global, "optimo_global_deseabilidad.csv", row.names = FALSE)
write.csv(comparacion, "comparacion_soluciones.csv", row.names = FALSE)