Objetivo Ajustar un modelo de segundo orden para la respuesta Y, diagnosticar los supuestos del modelo, construir gráficas de contorno y superficies de respuesta, y ubicar el punto óptimo mediante análisis estacionario, optimización acotada y escalonamiento ascendente o descendente.
Se trabajó con un diseño Box-Behnken con tres
factores codificados en los niveles -1, 0 y
1:
| Factor | Símbolo | Nivel -1 | Nivel 0 | Nivel 1 |
|---|---|---|---|---|
| Temperatura, °C | \(x_1\) | 120 | 145 | 170 |
| Contenido de humedad, % | \(x_2\) | 25 | 35 | 45 |
| Velocidad de tornillo | \(x_3\) | 800 | 900 | 1000 |
La transformación usada para regresar de escala codificada a escala real es:
\[ T = 145 + 25x_1, \qquad H = 35 + 10x_2, \qquad V = 900 + 100x_3 \]
archivo <- params$archivo_excel
# Si el Excel está en la misma carpeta del RMarkdown, se leerá automáticamente.
# Si no se encuentra, el script usa una base interna de respaldo con los datos del archivo cargado.
if (file.exists(archivo)) {
datos_raw <- readxl::read_excel(archivo, sheet = 1) |> janitor::clean_names()
} else {
datos_raw <- tibble::tribble(
~name, ~run_no_in_std_order, ~run_no, ~run_no_std_rp, ~temperatura, ~porc_humedad, ~vel_tornillo, ~y,
1, 1, 1, 1, -1, -1, 0, 3.7,
2, 2, 2, 2, 1, -1, 0, 2.0,
3, 3, 3, 3, -1, 1, 0, 33.4,
4, 4, 4, 4, 1, 1, 0, 4.4,
5, 5, 5, 5, -1, 0, -1, 10.0,
6, 6, 6, 6, 1, 0, -1, 4.5,
7, 7, 7, 7, -1, 0, 1, 6.0,
8, 8, 8, 8, 1, 0, 1, 3.3,
9, 9, 9, 9, 0, -1, -1, 3.9,
10, 10, 10, 10, 0, 1, -1, 31.0,
11, 11, 11, 11, 0, -1, 1, 3.7,
12, 12, 12, 12, 0, 1, 1, 21.0,
13, 13, 13, 13, 0, 0, 0, 3.5,
14, 14, 14, 14, 0, 0, 0, 4.5,
15, 15, 15, 15, 0, 0, 0, 6.4
)
}
# Normalización de nombres para trabajar cómodo en R.
datos <- datos_raw |>
dplyr::rename(
x1 = temperatura,
x2 = porc_humedad,
x3 = vel_tornillo,
y = y
) |>
dplyr::mutate(
dplyr::across(c(x1, x2, x3, y), as.numeric),
temperatura_c = 145 + 25 * x1,
humedad_pct = 35 + 10 * x2,
velocidad_tornillo = 900 + 100 * x3,
tratamiento = interaction(x1, x2, x3, drop = TRUE)
)
if (!"run_no" %in% names(datos)) datos$run_no <- seq_len(nrow(datos))
datos <- datos |> dplyr::arrange(run_no)
DT::datatable(
datos,
caption = "Matriz experimental Box-Behnken con factores codificados, valores reales y respuesta Y",
options = list(pageLength = 15, scrollX = TRUE)
)resumen_y <- datos |>
dplyr::summarise(
n = dplyr::n(),
media = mean(y),
desviacion = sd(y),
minimo = min(y),
q1 = quantile(y, 0.25),
mediana = median(y),
q3 = quantile(y, 0.75),
maximo = max(y)
)
DT::datatable(round(resumen_y, 4), caption = "Resumen descriptivo de la respuesta Y")p_inicial <- ggplot(datos, aes(x = factor(run_no), y = y, text = paste0(
"Corrida: ", run_no,
"<br>x1: ", x1,
"<br>x2: ", x2,
"<br>x3: ", x3,
"<br>Y: ", y
))) +
geom_col(fill = "#0b5cad") +
labs(
title = "Respuesta observada por corrida experimental",
x = "Número de corrida",
y = "Respuesta Y"
) +
theme_minimal(base_size = 13)
plotly::ggplotly(p_inicial, tooltip = "text")El modelo cuadrático ajustado es:
\[ Y = \beta_0 + \beta_1x_1 + \beta_2x_2 + \beta_3x_3 + \beta_{11}x_1^2 + \beta_{22}x_2^2 + \beta_{33}x_3^2 + \beta_{12}x_1x_2 + \beta_{13}x_1x_3 + \beta_{23}x_2x_3 + \varepsilon \]
modelo_lm <- lm(
y ~ x1 + x2 + x3 + I(x1^2) + I(x2^2) + I(x3^2) + x1:x2 + x1:x3 + x2:x3,
data = datos
)
modelo_rsm <- rsm::rsm(y ~ SO(x1, x2, x3), data = datos)
coeficientes <- broom::tidy(modelo_lm, conf.int = TRUE) |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5)))
DT::datatable(
coeficientes,
caption = "Coeficientes estimados del modelo de segundo orden",
options = list(pageLength = 12, scrollX = TRUE)
)metricas_modelo <- broom::glance(modelo_lm) |>
dplyr::select(r.squared, adj.r.squared, sigma, statistic, p.value, df, logLik, AIC, BIC, deviance, df.residual) |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5)))
DT::datatable(metricas_modelo, caption = "Métricas globales del modelo")anova_modelo <- anova(modelo_lm) |>
as.data.frame() |>
tibble::rownames_to_column("Fuente") |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5)))
DT::datatable(anova_modelo, caption = "ANOVA secuencial del modelo cuadrático")La falta de ajuste se estima comparando el modelo cuadrático contra un modelo saturado por combinación de tratamientos. Esta prueba solo tiene sentido cuando existen repeticiones; en este diseño hay repeticiones en el punto central.
if (any(table(datos$tratamiento) > 1)) {
modelo_puro_error <- lm(y ~ tratamiento, data = datos)
falta_ajuste <- anova(modelo_lm, modelo_puro_error) |>
as.data.frame() |>
tibble::rownames_to_column("Modelo") |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5)))
DT::datatable(falta_ajuste, caption = "Prueba de falta de ajuste frente a error puro")
} else {
htmltools::HTML("<div class='warn-box'>No hay repeticiones suficientes para separar falta de ajuste y error puro.</div>")
}En diseños pequeños como este, las pruebas formales tienen baja potencia. Por eso se reportan los test solicitados y se complementan con gráficas de diagnóstico. La interpretación debe integrar ambas partes.
diag_df <- broom::augment(modelo_lm) |>
dplyr::bind_cols(datos |> dplyr::select(run_no, x1, x2, x3, temperatura_c, humedad_pct, velocidad_tornillo)) |>
dplyr::mutate(
resid_estandar = rstandard(modelo_lm),
resid_student = rstudent(modelo_lm),
cooks_d = cooks.distance(modelo_lm)
)
DT::datatable(
diag_df |> dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5))),
caption = "Tabla de diagnóstico: ajustados, residuales, residuales estandarizados y distancia de Cook",
options = list(pageLength = 15, scrollX = TRUE)
)# El modelo se ajustó con los datos ordenados por run_no.
# Durbin-Watson evalúa autocorrelación serial de primer orden en los residuales.
test_dw <- lmtest::dwtest(modelo_lm)
dw_tabla <- tibble::tibble(
prueba = "Durbin-Watson",
estadistico = unname(test_dw$statistic),
p_valor = test_dw$p.value,
interpretacion = dplyr::case_when(
test_dw$p.value < 0.05 ~ "Evidencia de autocorrelación serial en residuales",
TRUE ~ "No se detecta autocorrelación serial significativa"
)
) |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5)))
DT::datatable(dw_tabla, caption = "Test de Durbin-Watson")Para aplicar Levene a residuales en un diseño con pocos datos, se agrupan los residuales por terciles de valores ajustados. Esta es una estrategia práctica para evaluar si la dispersión cambia según el nivel de respuesta esperado.
breaks_levene <- unique(quantile(diag_df$.fitted, probs = seq(0, 1, length.out = 4), na.rm = TRUE))
if (length(breaks_levene) >= 3) {
diag_df <- diag_df |>
dplyr::mutate(grupo_levene = cut(.fitted, breaks = breaks_levene, include.lowest = TRUE))
test_levene <- car::leveneTest(.resid ~ grupo_levene, data = diag_df, center = median)
levene_tabla <- as.data.frame(test_levene) |>
tibble::rownames_to_column("Fuente") |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5)))
DT::datatable(levene_tabla, caption = "Test de Levene sobre residuales agrupados por valores ajustados")
} else {
htmltools::HTML("<div class='warn-box'>No se pudieron crear grupos suficientes para aplicar Levene.</div>")
}resid_std <- as.numeric(scale(residuals(modelo_lm)))
test_ks <- ks.test(resid_std, "pnorm")
test_shapiro <- shapiro.test(residuals(modelo_lm))
tabla_normalidad <- tibble::tibble(
prueba = c("Kolmogorov-Smirnov con residuales estandarizados", "Shapiro-Wilk complementario"),
estadistico = c(unname(test_ks$statistic), unname(test_shapiro$statistic)),
p_valor = c(test_ks$p.value, test_shapiro$p.value),
interpretacion = dplyr::case_when(
p_valor < 0.05 ~ "Se rechaza normalidad al 5%",
TRUE ~ "No se rechaza normalidad al 5%"
)
) |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5)))
DT::datatable(tabla_normalidad, caption = "Pruebas de normalidad de residuales")p_resid_fit <- ggplot(diag_df, aes(x = .fitted, y = .resid, text = paste0(
"Corrida: ", run_no,
"<br>Ajustado: ", round(.fitted, 3),
"<br>Residual: ", round(.resid, 3),
"<br>Cook D: ", round(cooks_d, 4)
))) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray40") +
geom_point(size = 3, color = "#0b5cad") +
labs(title = "Residuales vs valores ajustados", x = "Valor ajustado", y = "Residual") +
theme_minimal(base_size = 13)
plotly::ggplotly(p_resid_fit, tooltip = "text")# Q-Q plot interactivo construido directamente con plotly.
# Se evita ggplotly(stat_qq_line), porque en algunas versiones genera el error:
# "objeto 'x' no encontrado" al convertir la línea de referencia.
qq_df <- tibble::tibble(
teorico = stats::qnorm(stats::ppoints(length(diag_df$.std.resid))),
observado = sort(diag_df$.std.resid)
)
# Línea de referencia del Q-Q plot usando los cuartiles 25% y 75%.
q_teo <- stats::qnorm(c(0.25, 0.75))
q_obs <- stats::quantile(qq_df$observado, probs = c(0.25, 0.75), names = FALSE)
pendiente_qq <- diff(q_obs) / diff(q_teo)
intercepto_qq <- q_obs[1] - pendiente_qq * q_teo[1]
rango_qq <- range(qq_df$teorico)
linea_qq <- tibble::tibble(
teorico = rango_qq,
observado = intercepto_qq + pendiente_qq * rango_qq
)
plotly::plot_ly(
data = qq_df,
x = ~teorico,
y = ~observado,
type = "scatter",
mode = "markers",
marker = list(color = "#0b5cad", size = 8),
text = ~paste0(
"Cuantil teórico: ", round(teorico, 3),
"<br>Cuantil observado: ", round(observado, 3)
),
hoverinfo = "text",
name = "Residuales"
) |>
plotly::add_lines(
data = linea_qq,
x = ~teorico,
y = ~observado,
line = list(color = "#d99000", width = 3),
hoverinfo = "skip",
name = "Referencia"
) |>
plotly::layout(
title = "Gráfico Q-Q de residuales estandarizados",
xaxis = list(title = "Cuantil teórico"),
yaxis = list(title = "Cuantil observado")
)p_scale <- ggplot(diag_df, aes(x = .fitted, y = sqrt(abs(resid_estandar)), text = paste0(
"Corrida: ", run_no,
"<br>Ajustado: ", round(.fitted, 3),
"<br>|Residual est.|^0.5: ", round(sqrt(abs(resid_estandar)), 3)
))) +
geom_point(size = 3, color = "#0b5cad") +
geom_smooth(method = "loess", se = FALSE, color = "#d99000") +
labs(title = "Scale-Location", x = "Valor ajustado", y = expression(sqrt('|Residual estandarizado|'))) +
theme_minimal(base_size = 13)
plotly::ggplotly(p_scale, tooltip = "text")p_run <- ggplot(diag_df, aes(x = run_no, y = .resid, text = paste0(
"Corrida: ", run_no,
"<br>Residual: ", round(.resid, 3)
))) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray40") +
geom_line(color = "#0b5cad") +
geom_point(size = 3, color = "#0b5cad") +
labs(title = "Residuales en orden de corrida", x = "Orden de corrida", y = "Residual") +
theme_minimal(base_size = 13)
plotly::ggplotly(p_run, tooltip = "text")p_cook <- ggplot(diag_df, aes(x = run_no, y = cooks_d, text = paste0(
"Corrida: ", run_no,
"<br>Cook D: ", round(cooks_d, 4)
))) +
geom_col(fill = "#0b5cad") +
geom_hline(yintercept = 4 / nrow(datos), linetype = "dashed", color = "#d99000") +
labs(title = "Distancia de Cook", x = "Orden de corrida", y = "Cook's D") +
theme_minimal(base_size = 13)
plotly::ggplotly(p_cook, tooltip = "text")etiquetas <- list(
x1 = "Temperatura codificada, x1 (-1=120, 0=145, 1=170 °C)",
x2 = "Humedad codificada, x2 (-1=25, 0=35, 1=45 %)",
x3 = "Velocidad codificada, x3 (-1=800, 0=900, 1=1000)"
)
pred_fun <- function(z) {
nuevo <- data.frame(x1 = z[1], x2 = z[2], x3 = z[3])
as.numeric(predict(modelo_lm, newdata = nuevo))
}
codificar_a_real <- function(x1, x2, x3) {
tibble::tibble(
temperatura_c = 145 + 25 * x1,
humedad_pct = 35 + 10 * x2,
velocidad_tornillo = 900 + 100 * x3
)
}
crear_malla <- function(var_x, var_y, fijo_val = 0, n = 90) {
secuencia <- seq(-1, 1, length.out = n)
malla <- expand.grid(a = secuencia, b = secuencia)
names(malla) <- c(var_x, var_y)
faltante <- setdiff(c("x1", "x2", "x3"), c(var_x, var_y))
malla[[faltante]] <- fijo_val
malla <- malla |> dplyr::select(x1, x2, x3)
malla$yhat <- as.numeric(predict(modelo_lm, newdata = malla))
malla
}
plot_contorno <- function(var_x, var_y, fijo_val = 0) {
faltante <- setdiff(c("x1", "x2", "x3"), c(var_x, var_y))
malla <- crear_malla(var_x, var_y, fijo_val = fijo_val, n = 100)
p <- ggplot(malla, aes(x = .data[[var_x]], y = .data[[var_y]], z = yhat)) +
geom_raster(aes(fill = yhat), interpolate = TRUE) +
geom_contour(color = "white", alpha = 0.75, linewidth = 0.35) +
geom_point(
data = datos,
aes(x = .data[[var_x]], y = .data[[var_y]]),
inherit.aes = FALSE,
color = "black",
size = 2
) +
scale_fill_viridis_c(name = "Y predicha", option = "C") +
labs(
title = paste0("Contorno: ", var_x, " vs ", var_y),
subtitle = paste0(faltante, " fijo en ", fijo_val, " escala codificada"),
x = etiquetas[[var_x]],
y = etiquetas[[var_y]]
) +
theme_minimal(base_size = 13)
plotly::ggplotly(p)
}
plot_superficie <- function(var_x, var_y, fijo_val = 0, n = 70) {
faltante <- setdiff(c("x1", "x2", "x3"), c(var_x, var_y))
secuencia <- seq(-1, 1, length.out = n)
malla <- crear_malla(var_x, var_y, fijo_val = fijo_val, n = n)
# expand.grid deja var_x cambiando más rápido; se transpone para la convención de plotly.
zmat <- t(matrix(malla$yhat, nrow = n, ncol = n, byrow = FALSE))
plotly::plot_ly(x = secuencia, y = secuencia, z = zmat) |>
plotly::add_surface(colorscale = "Viridis") |>
plotly::layout(
title = paste0("Superficie de respuesta: ", var_x, " vs ", var_y, " | ", faltante, " = ", fijo_val),
scene = list(
xaxis = list(title = etiquetas[[var_x]]),
yaxis = list(title = etiquetas[[var_y]]),
zaxis = list(title = "Y predicha")
)
)
}La estrategia usada para hallar el óptimo combina tres enfoques:
Para el modelo:
\[ \hat{Y} = \beta_0 + b'x + x'Bx \]
el punto estacionario se obtiene con:
\[ x_s = -\frac{1}{2}B^{-1}b \]
beta <- coef(modelo_lm)
b <- c(
x1 = beta[["x1"]],
x2 = beta[["x2"]],
x3 = beta[["x3"]]
)
B <- matrix(0, nrow = 3, ncol = 3, dimnames = list(c("x1", "x2", "x3"), c("x1", "x2", "x3")))
B["x1", "x1"] <- beta[["I(x1^2)"]]
B["x2", "x2"] <- beta[["I(x2^2)"]]
B["x3", "x3"] <- beta[["I(x3^2)"]]
B["x1", "x2"] <- B["x2", "x1"] <- beta[["x1:x2"]] / 2
B["x1", "x3"] <- B["x3", "x1"] <- beta[["x1:x3"]] / 2
B["x2", "x3"] <- B["x3", "x2"] <- beta[["x2:x3"]] / 2
x_estacionario <- tryCatch(
-0.5 * solve(B, b),
error = function(e) rep(NA_real_, 3)
)
hessiano <- 2 * B
autovalores <- eigen(hessiano)$values
clasificacion <- dplyr::case_when(
all(is.na(x_estacionario)) ~ "No calculable: matriz B singular",
all(autovalores > 0) ~ "Mínimo local",
all(autovalores < 0) ~ "Máximo local",
TRUE ~ "Punto silla"
)
punto_est <- tibble::tibble(
x1 = x_estacionario[1],
x2 = x_estacionario[2],
x3 = x_estacionario[3],
y_pred = if (any(is.na(x_estacionario))) NA_real_ else pred_fun(x_estacionario),
dentro_region = all(x_estacionario >= -1 & x_estacionario <= 1),
clasificacion = clasificacion
) |>
dplyr::bind_cols(codificar_a_real(x_estacionario[1], x_estacionario[2], x_estacionario[3])) |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5)))
auto_tabla <- tibble::tibble(
componente = paste0("lambda_", seq_along(autovalores)),
autovalor_hessiano = round(autovalores, 5)
)
DT::datatable(punto_est, caption = "Punto estacionario del modelo cuadrático")DT::datatable(auto_tabla, caption = "Autovalores del Hessiano para clasificar el punto estacionario")# Resultado canónico desde el paquete rsm.
# Es útil para validar el análisis estacionario calculado manualmente.
rsm::canonical(modelo_rsm)## $xs
## x1 x2 x3
## -0.09157068 -0.66683781 0.06866606
##
## $eigen
## eigen() decomposition
## $values
## [1] 8.920897 2.332509 -2.590905
##
## $vectors
## [,1] [,2] [,3]
## x1 0.3127667 -0.07444555 0.94690805
## x2 -0.9291750 0.18279137 0.32128040
## x3 0.1970045 0.98032911 0.01200189optim_multistart <- function(maximizar = TRUE) {
iniciales <- as.matrix(expand.grid(x1 = c(-1, 0, 1), x2 = c(-1, 0, 1), x3 = c(-1, 0, 1)))
resultados <- purrr::map(seq_len(nrow(iniciales)), function(i) {
optim(
par = iniciales[i, ],
fn = function(z) if (maximizar) -pred_fun(z) else pred_fun(z),
method = "L-BFGS-B",
lower = c(-1, -1, -1),
upper = c(1, 1, 1)
)
})
valores <- purrr::map_dbl(resultados, ~ pred_fun(.x$par))
idx <- if (maximizar) which.max(valores) else which.min(valores)
z <- resultados[[idx]]$par
tibble::tibble(
objetivo = ifelse(maximizar, "Máximo acotado", "Mínimo acotado"),
x1 = z[1],
x2 = z[2],
x3 = z[3],
y_pred = valores[idx]
) |>
dplyr::bind_cols(codificar_a_real(z[1], z[2], z[3]))
}
opt_max <- optim_multistart(maximizar = TRUE)
opt_min <- optim_multistart(maximizar = FALSE)
optimos <- dplyr::bind_rows(opt_max, opt_min) |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5)))
DT::datatable(optimos, caption = "Máximo y mínimo predichos dentro de la región experimental codificada [-1,1]")objetivo_usuario <- tolower(params$objetivo)
opt_elegido <- if (objetivo_usuario == "min") opt_min else opt_max
nuevo_opt <- data.frame(
x1 = opt_elegido$x1,
x2 = opt_elegido$x2,
x3 = opt_elegido$x3
)
pred_intervalo <- as.data.frame(predict(modelo_lm, newdata = nuevo_opt, interval = "prediction", level = 0.95)) |>
dplyr::bind_cols(opt_elegido) |>
dplyr::select(objetivo, x1, x2, x3, temperatura_c, humedad_pct, velocidad_tornillo, fit, lwr, upr) |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5)))
DT::datatable(pred_intervalo, caption = "Punto óptimo elegido e intervalo de predicción al 95%")El escalonamiento usa la dirección del gradiente del modelo:
\[ \nabla \hat{Y}(x) = b + 2Bx \]
Para maximizar, se avanza en dirección del gradiente. Para minimizar, se avanza en la dirección contraria.
grad_fun <- function(z) {
as.numeric(b + 2 * B %*% z)
}
camino_escalonado <- function(objetivo = c("max", "min"), inicio = c(0, 0, 0), paso = 0.10, n_pasos = 30) {
objetivo <- match.arg(objetivo)
z <- inicio
salida <- list()
for (k in 0:n_pasos) {
salida[[k + 1]] <- tibble::tibble(
paso = k,
objetivo = ifelse(objetivo == "max", "Ascendente", "Descendente"),
x1 = z[1],
x2 = z[2],
x3 = z[3],
y_pred = pred_fun(z)
)
g <- grad_fun(z)
if (objetivo == "min") g <- -g
norma <- sqrt(sum(g^2))
if (norma < 1e-10) break
z_nuevo <- z + paso * g / norma
z_nuevo <- pmin(pmax(z_nuevo, -1), 1)
if (all(abs(z_nuevo - z) < 1e-8)) break
z <- z_nuevo
}
dplyr::bind_rows(salida) |>
dplyr::bind_cols(codificar_a_real(
dplyr::bind_rows(salida)$x1,
dplyr::bind_rows(salida)$x2,
dplyr::bind_rows(salida)$x3
))
}
camino_max <- camino_escalonado("max", paso = 0.10, n_pasos = 30)
camino_min <- camino_escalonado("min", paso = 0.10, n_pasos = 30)
camino_total <- dplyr::bind_rows(camino_max, camino_min) |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5)))
DT::datatable(
camino_total,
caption = "Camino de escalonamiento ascendente y descendente desde el centro del diseño",
options = list(pageLength = 20, scrollX = TRUE)
)p_escalonamiento <- ggplot(camino_total, aes(x = paso, y = y_pred, color = objetivo, text = paste0(
"Paso: ", paso,
"<br>Tipo: ", objetivo,
"<br>x1: ", x1,
"<br>x2: ", x2,
"<br>x3: ", x3,
"<br>Temperatura: ", temperatura_c,
"<br>Humedad: ", humedad_pct,
"<br>Velocidad: ", velocidad_tornillo,
"<br>Y predicha: ", y_pred
))) +
geom_line(linewidth = 1) +
geom_point(size = 2.8) +
labs(
title = "Evolución de Y predicha durante el escalonamiento",
x = "Paso secuencial",
y = "Y predicha",
color = "Camino"
) +
theme_minimal(base_size = 13)
plotly::ggplotly(p_escalonamiento, tooltip = "text")if (objetivo_usuario == "min") {
recomendacion <- opt_min
texto_obj <- "minimizar"
} else {
recomendacion <- opt_max
texto_obj <- "maximizar"
}
rec_tabla <- recomendacion |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5)))
htmltools::HTML(paste0(
"<div class='ok-box'>",
"<b>Recomendación:</b> Si el objetivo es <b>", texto_obj, "</b> la respuesta Y, usar como punto candidato el óptimo acotado dentro de la región experimental. ",
"Este punto debe validarse con una corrida confirmatoria porque el óptimo proviene del modelo ajustado y no de una observación directa.",
"</div>"
))Paso 1. Se trabajó con variables codificadas \((-1,0,1)\) porque el diseño Box-Behnken se interpreta mejor en escala codificada y porque las unidades originales tienen magnitudes diferentes.
Paso 2. Se ajustó un modelo cuadrático completo con efectos lineales, términos de curvatura e interacciones. Este modelo permite representar superficies con máximos, mínimos o puntos silla.
Paso 3. Se validó el modelo mediante ANOVA, métricas de ajuste y diagnósticos de residuales. Se aplicaron Durbin-Watson, Levene y Kolmogorov-Smirnov, junto con gráficas de residuales, Q-Q, Scale-Location, orden de corrida y distancia de Cook.
Paso 4. Se calculó el punto estacionario resolviendo \(x_s=-\frac{1}{2}B^{-1}b\). Luego se clasificó con los autovalores del Hessiano. Si los autovalores tienen signos mezclados, el punto no es máximo ni mínimo, sino punto silla.
Paso 5. Como el óptimo útil debe estar dentro del dominio experimental, se aplicó optimización acotada en \([-1,1]^3\). Así se obtuvo el máximo y el mínimo predicho dentro de los niveles estudiados.
Paso 6. Finalmente, se construyó un camino de escalonamiento. Para maximizar se avanza en la dirección del gradiente; para minimizar se avanza en dirección opuesta. Este camino sirve para proponer corridas confirmatorias secuenciales hasta llegar a la mejor zona experimental.