Objetivo Analizar un Diseño
Central Compuesto Rotable para tres factores codificados
x1, x2 y x3, ajustar un modelo de
segundo orden para la respuesta Y, evaluar 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.
Los datos corresponden a un experimento para optimizar el crecimiento de cristales en función de tres variables codificadas:
\[ x_1, \quad x_2, \quad x_3 \]
Se buscan valores altos de la respuesta Y, por lo que el objetivo por defecto del análisis es:
\[ \max(Y) \]
El diseño indicado es un Diseño Central Compuesto Rotable, compuesto por:
-1 y
1.±α sobre
cada eje.(0, 0, 0) para estimar
error puro y curvatura.Para tres factores, el valor teórico de rotabilidad es:
\[ \alpha = (2^k)^{1/4} = (2^3)^{1/4} = 1.68179 \]
archivo <- params$archivo_excel
hoja <- params$hoja
respuesta_param <- params$respuesta
# Si el Excel está en la misma carpeta del RMarkdown, se leerá automáticamente.
# Si no se encuentra, se usa una base interna de respaldo construida con los datos del archivo entregado.
if (file.exists(archivo)) {
datos_raw <- readxl::read_excel(archivo, sheet = hoja)
} else {
datos_raw <- tibble::tribble(
~name, ~run.no.in.std.order, ~run.no, ~run.no.std.rp, ~Block.ccd, ~X1, ~X2, ~X3, ~Y,
"C1.1", "C1.1", 1, "C1.1", 1, -1, -1, -1, 66,
"C1.2", "C1.2", 2, "C1.2", 1, 1, -1, -1, 70,
"C1.3", "C1.3", 3, "C1.3", 1, -1, 1, -1, 78,
"C1.4", "C1.4", 4, "C1.4", 1, 1, 1, -1, 60,
"C1.5", "C1.5", 5, "C1.5", 1, -1, -1, 1, 80,
"C1.6", "C1.6", 6, "C1.6", 1, 1, -1, 1, 70,
"C1.7", "C1.7", 7, "C1.7", 1, -1, 1, 1, 100,
"C1.8", "C1.8", 8, "C1.8", 1, 1, 1, 1, 75,
"C1.9", "C1.9", 9, "C1.9", 1, 0, 0, 0, 113,
"C1.10", "C1.10", 10, "C1.10", 1, 0, 0, 0, 100,
"C1.11", "C1.11", 11, "C1.11", 1, 0, 0, 0, 118,
"C1.12", "C1.12", 12, "C1.12", 1, 0, 0, 0, 88,
"C1.13", "C1.13", 13, "C1.13", 1, 0, 0, 0, 100,
"C1.14", "C1.14", 14, "C1.14", 1, 0, 0, 0, 85,
"S2.1", "S2.1", 15, "S2.1", 2, -1.68179283050743, 0, 0, 100,
"S2.2", "S2.2", 16, "S2.2", 2, 1.68179283050743, 0, 0, 80,
"S2.3", "S2.3", 17, "S2.3", 2, 0, -1.68179283050743, 0, 68,
"S2.4", "S2.4", 18, "S2.4", 2, 0, 1.68179283050743, 0, 63,
"S2.5", "S2.5", 19, "S2.5", 2, 0, 0, -1.68179283050743, 65,
"S2.6", "S2.6", 20, "S2.6", 2, 0, 0, 1.68179283050743, 82
)
}
nombres_originales <- names(datos_raw)
# Identificación automática de la respuesta.
y_idx <- which(tolower(names(datos_raw)) == tolower(respuesta_param))
if (length(y_idx) == 0) {
y_idx <- which(tolower(names(datos_raw)) %in% c("y", "respuesta", "response"))
}
if (length(y_idx) == 0) stop("No se encontró la columna de respuesta. Revise params$respuesta.")
y_idx <- y_idx[1]
# En este archivo, los tres factores codificados son las tres columnas inmediatamente anteriores a Y.
x_idx <- (y_idx - 3):(y_idx - 1)
if (min(x_idx) < 1) stop("No fue posible identificar tres columnas de factores antes de la respuesta Y.")
# Identificación del bloque, si existe.
block_idx <- grep("block|bloque", names(datos_raw), ignore.case = TRUE)
if (length(block_idx) == 0) {
bloque_tmp <- rep(1, nrow(datos_raw))
} else {
bloque_tmp <- datos_raw[[block_idx[1]]]
}
# Identificación del número de corrida, si existe.
run_idx <- grep("run.no$|run_no$|corrida|orden", names(datos_raw), ignore.case = TRUE)
if (length(run_idx) == 0) {
corrida_tmp <- seq_len(nrow(datos_raw))
} else {
corrida_tmp <- datos_raw[[run_idx[1]]]
if (!is.numeric(corrida_tmp)) corrida_tmp <- seq_len(nrow(datos_raw))
}
datos <- tibble::tibble(
corrida = as.integer(corrida_tmp),
block = factor(bloque_tmp),
x1 = as.numeric(datos_raw[[x_idx[1]]]),
x2 = as.numeric(datos_raw[[x_idx[2]]]),
x3 = as.numeric(datos_raw[[x_idx[3]]]),
y = as.numeric(datos_raw[[y_idx]])
)
datos <- datos |>
dplyr::mutate(
tratamiento = interaction(round(x1, 8), round(x2, 8), round(x3, 8), drop = TRUE),
tipo_punto = dplyr::case_when(
abs(x1) < 1e-8 & abs(x2) < 1e-8 & abs(x3) < 1e-8 ~ "Central",
rowSums(abs(cbind(x1, x2, x3)) > 1 + 1e-6) == 1 ~ "Axial/estrella",
TRUE ~ "Factorial"
),
tipo_punto = factor(tipo_punto, levels = c("Factorial", "Axial/estrella", "Central"))
) |>
dplyr::arrange(corrida)
alpha <- params$alpha
alpha_detectado <- max(abs(unlist(datos[, c("x1", "x2", "x3")], use.names = FALSE)))
DT::datatable(
datos,
caption = "Matriz experimental del Diseño Central Compuesto Rotable",
options = list(pageLength = 20, scrollX = TRUE)
)resumen_diseno <- tibble::tibble(
Componente = c(
"Número de factores",
"Número total de corridas",
"Puntos factoriales",
"Puntos axiales o estrella",
"Puntos centrales",
"Bloques experimentales",
"Alpha teórico rotable, k = 3",
"Alpha detectado en la matriz"
),
Valor = c(
3,
nrow(datos),
sum(datos$tipo_punto == "Factorial"),
sum(datos$tipo_punto == "Axial/estrella"),
sum(datos$tipo_punto == "Central"),
nlevels(datos$block),
round((2^3)^(1/4), 6),
round(alpha_detectado, 6)
)
)
DT::datatable(resumen_diseno, options = list(dom = "t"), rownames = FALSE)Respuesta al literal a). El diseño empleado es un
Diseño Central Compuesto Rotable para tres factores. La
matriz contiene puntos factoriales, axiales y centrales. El valor axial
detectado coincide con el valor rotable esperado para
k = 3, aproximadamente α = 1.68179.
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", options = list(dom = "t"))p_corrida <- ggplot(datos, aes(x = factor(corrida), y = y, fill = tipo_punto,
text = paste0(
"Corrida: ", corrida,
"<br>Tipo: ", tipo_punto,
"<br>Bloque: ", block,
"<br>x1: ", round(x1, 4),
"<br>x2: ", round(x2, 4),
"<br>x3: ", round(x3, 4),
"<br>Y: ", y
))) +
geom_col() +
labs(
title = "Respuesta observada por corrida experimental",
x = "Corrida",
y = "Respuesta Y",
fill = "Tipo de punto"
) +
theme_minimal(base_size = 13)
plotly::ggplotly(p_corrida, tooltip = "text")plotly::plot_ly(
datos,
x = ~x1, y = ~x2, z = ~x3,
color = ~y,
colors = paleta_viridis,
type = "scatter3d",
mode = "markers",
marker = list(size = 6),
text = ~paste0(
"Corrida: ", corrida,
"<br>Tipo: ", tipo_punto,
"<br>Y: ", y
),
hoverinfo = "text"
) |>
plotly::layout(
title = "Distribución 3D de los puntos experimentales",
scene = list(
xaxis = list(title = "x1"),
yaxis = list(title = "x2"),
zaxis = list(title = "x3")
)
)El modelo cuadrático completo en variables codificadas 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 \]
Si el diseño tiene bloques y params$usar_bloque = TRUE,
se incluye el bloque como efecto no controlable para depurar
variabilidad experimental. La optimización se realiza únicamente sobre
x1, x2 y x3.
usar_bloque_en_modelo <- isTRUE(params$usar_bloque) && nlevels(datos$block) > 1
formula_so <- y ~ x1 + x2 + x3 + I(x1^2) + I(x2^2) + I(x3^2) + x1:x2 + x1:x3 + x2:x3
formula_modelo <- if (usar_bloque_en_modelo) {
y ~ block + x1 + x2 + x3 + I(x1^2) + I(x2^2) + I(x3^2) + x1:x2 + x1:x3 + x2:x3
} else {
formula_so
}
modelo_lm <- lm(formula_modelo, data = datos)
modelo_sin_bloque <- lm(formula_so, 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 del modelo cuadrático ajustado",
options = list(pageLength = 15, scrollX = TRUE)
)get_coef <- function(nombre, modelo = modelo_lm) {
cc <- coef(modelo)
if (nombre %in% names(cc)) as.numeric(cc[[nombre]]) else 0
}
fmt <- function(x) format(round(x, 4), nsmall = 4, trim = TRUE)
fmt_signed <- function(x) ifelse(x >= 0, paste0(" + ", fmt(x)), paste0(" - ", fmt(abs(x))))
ecuacion <- paste0(
"Ŷ = ", fmt(get_coef("(Intercept)")),
fmt_signed(get_coef("x1")), "x1",
fmt_signed(get_coef("x2")), "x2",
fmt_signed(get_coef("x3")), "x3",
fmt_signed(get_coef("I(x1^2)")), "x1²",
fmt_signed(get_coef("I(x2^2)")), "x2²",
fmt_signed(get_coef("I(x3^2)")), "x3²",
fmt_signed(get_coef("x1:x2")), "x1x2",
fmt_signed(get_coef("x1:x3")), "x1x3",
fmt_signed(get_coef("x2:x3")), "x2x3"
)
htmltools::HTML(paste0("<div class='info-box'><b>Modelo ajustado en escala codificada:</b><br><span style='font-size:18px;'>", ecuacion, "</span></div>"))anova_modelo <- as.data.frame(anova(modelo_lm))
anova_modelo$Fuente <- rownames(anova_modelo)
rownames(anova_modelo) <- NULL
anova_modelo <- anova_modelo |>
dplyr::select(Fuente, dplyr::everything()) |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5)))
DT::datatable(anova_modelo, caption = "ANOVA del modelo ajustado", options = list(pageLength = 20, scrollX = TRUE))glance_modelo <- broom::glance(modelo_lm) |>
dplyr::select(r.squared, adj.r.squared, sigma, statistic, p.value, df, df.residual, AIC, BIC) |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5)))
DT::datatable(glance_modelo, caption = "Métricas globales del modelo", options = list(dom = "t", scrollX = TRUE))componentes <- broom::tidy(modelo_lm) |>
dplyr::filter(term != "(Intercept)") |>
dplyr::mutate(
decision_5pct = dplyr::case_when(
p.value < 0.01 ~ "Altamente significativo",
p.value < 0.05 ~ "Significativo",
p.value < 0.10 ~ "Marginal",
TRUE ~ "No significativo"
),
dplyr::across(where(is.numeric), ~ round(.x, 5))
) |>
dplyr::arrange(p.value)
DT::datatable(componentes, caption = "Componentes ordenados por importancia estadística", options = list(pageLength = 15, scrollX = TRUE))# Pruebas parciales tipo drop1 para valorar la contribución de cada término.
# Debe interpretarse con criterio jerárquico: no conviene retirar términos lineales si aparecen en interacciones o cuadrados relevantes.
drop1_tab <- as.data.frame(drop1(modelo_lm, test = "F"))
drop1_tab$Termino <- rownames(drop1_tab)
rownames(drop1_tab) <- NULL
drop1_tab <- drop1_tab |>
dplyr::select(Termino, dplyr::everything()) |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5)))
DT::datatable(drop1_tab, caption = "Evaluación parcial de términos del modelo", options = list(pageLength = 15, scrollX = TRUE))Respuesta al literal b). El modelo de segundo orden se ajusta mediante mínimos cuadrados. Los componentes más importantes se identifican con la tabla de coeficientes, el ANOVA y las pruebas parciales. Para conservar la jerarquía del modelo de superficie de respuesta, la interpretación no debe basarse solo en eliminar términos por valor-p; también debe considerarse la estructura cuadrática completa.
La prueba de falta de ajuste compara el modelo cuadrático contra un modelo con medias por tratamiento. El error puro se estima a partir de corridas repetidas en puntos idénticos, especialmente en el centro del diseño.
modelo_puro <- lm(y ~ tratamiento, data = datos)
tabla_falta_ajuste <- tryCatch({
as.data.frame(anova(modelo_sin_bloque, modelo_puro))
}, error = function(e) NULL)
if (!is.null(tabla_falta_ajuste)) {
tabla_falta_ajuste$Modelo <- rownames(tabla_falta_ajuste)
rownames(tabla_falta_ajuste) <- NULL
tabla_falta_ajuste <- tabla_falta_ajuste |>
dplyr::select(Modelo, dplyr::everything()) |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5)))
DT::datatable(tabla_falta_ajuste, caption = "Prueba de falta de ajuste usando error puro", options = list(dom = "t", scrollX = TRUE))
} else {
htmltools::HTML("<div class='warn-box'>No fue posible calcular la prueba de falta de ajuste con la estructura actual de datos.</div>")
}if (!is.null(tabla_falta_ajuste) && "Pr(>F)" %in% names(tabla_falta_ajuste)) {
p_lof <- tabla_falta_ajuste$`Pr(>F)`[2]
if (!is.na(p_lof) && p_lof > 0.05) {
cat("<div class='ok-box'><b>Interpretación:</b> con un nivel de significancia de 5%, no se detecta falta de ajuste significativa. Esto favorece la adecuación del modelo cuadrático.</div>")
} else if (!is.na(p_lof)) {
cat("<div class='bad-box'><b>Interpretación:</b> se detecta falta de ajuste significativa. El modelo cuadrático podría no describir completamente la superficie o podrían existir observaciones influyentes.</div>")
}
}diagnostico <- tibble::tibble(
corrida = datos$corrida,
block = datos$block,
tipo_punto = datos$tipo_punto,
y_observado = datos$y,
y_ajustado = fitted(modelo_lm),
residual = residuals(modelo_lm),
residual_estandar = rstandard(modelo_lm),
leverage = hatvalues(modelo_lm),
cook = cooks.distance(modelo_lm)
)
DT::datatable(
diagnostico |> dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5))),
caption = "Datos de diagnóstico del modelo",
options = list(pageLength = 20, scrollX = TRUE)
)test_dw <- lmtest::dwtest(modelo_lm)
dw_tab <- tibble::tibble(
Prueba = "Durbin-Watson",
Estadistico = unname(test_dw$statistic),
p_value = test_dw$p.value,
Interpretacion = ifelse(test_dw$p.value > 0.05,
"No se detecta autocorrelación significativa de primer orden",
"Existe evidencia de autocorrelación de primer orden")
) |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5)))
DT::datatable(dw_tab, caption = "Prueba de independencia de residuales", options = list(dom = "t"), rownames = FALSE)# Levene requiere grupos. Aquí se agrupan los residuales por cuartiles del valor ajustado.
# Esto evalúa si la variabilidad residual cambia a través de la escala de respuesta predicha.
qs <- unique(quantile(diagnostico$y_ajustado, probs = seq(0, 1, 0.25), na.rm = TRUE))
if (length(qs) >= 3) {
diagnostico$grupo_levene <- cut(diagnostico$y_ajustado, breaks = qs, include.lowest = TRUE)
test_levene <- car::leveneTest(residual ~ grupo_levene, data = diagnostico, center = median)
levene_tab <- as.data.frame(test_levene)
levene_tab$Fuente <- rownames(levene_tab)
rownames(levene_tab) <- NULL
levene_tab <- levene_tab |>
dplyr::select(Fuente, dplyr::everything()) |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5)))
DT::datatable(levene_tab, caption = "Test de Levene sobre residuales agrupados por valores ajustados", options = list(dom = "t"), rownames = FALSE)
} else {
htmltools::HTML("<div class='warn-box'>No hay suficientes grupos distintos para aplicar Levene con cuartiles del valor ajustado.</div>")
}if (nlevels(diagnostico$tipo_punto) >= 2) {
test_levene_tipo <- car::leveneTest(residual ~ tipo_punto, data = diagnostico, center = median)
levene_tipo_tab <- as.data.frame(test_levene_tipo)
levene_tipo_tab$Fuente <- rownames(levene_tipo_tab)
rownames(levene_tipo_tab) <- NULL
levene_tipo_tab <- levene_tipo_tab |>
dplyr::select(Fuente, dplyr::everything()) |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5)))
DT::datatable(levene_tipo_tab, caption = "Test de Levene por tipo de punto experimental", options = list(dom = "t"), rownames = FALSE)
}res_std <- as.numeric(scale(residuals(modelo_lm)))
if (sd(res_std, na.rm = TRUE) > 0) {
test_ks <- suppressWarnings(ks.test(res_std, "pnorm"))
ks_tab <- tibble::tibble(
Prueba = "Kolmogorov-Smirnov sobre residuales estandarizados",
Estadistico = unname(test_ks$statistic),
p_value = test_ks$p.value,
Interpretacion = ifelse(test_ks$p.value > 0.05,
"No se rechaza normalidad al 5%",
"Se rechaza normalidad al 5%")
) |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 5)))
DT::datatable(ks_tab, caption = "Prueba de normalidad de Kolmogorov-Smirnov", options = list(dom = "t"), rownames = FALSE)
} else {
htmltools::HTML("<div class='warn-box'>No es posible aplicar Kolmogorov-Smirnov porque la desviación estándar residual es cero.</div>")
}Nota técnica. La prueba Kolmogorov-Smirnov aplicada a residuales estandarizados es una aproximación, porque la media y la desviación estándar se estiman con los propios datos. Por eso se complementa con el gráfico Q-Q y las demás gráficas de diagnóstico.
p_res_fit <- ggplot(diagnostico, aes(x = y_ajustado, y = residual_estandar,
text = paste0(
"Corrida: ", corrida,
"<br>Y observado: ", round(y_observado, 4),
"<br>Y ajustado: ", round(y_ajustado, 4),
"<br>Residual estándar: ", round(residual_estandar, 4)
))) +
geom_point(size = 3) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_smooth(method = "loess", se = FALSE) +
labs(
title = "Residuales estandarizados vs valores ajustados",
x = "Valores ajustados",
y = "Residual estandarizado"
) +
theme_minimal(base_size = 13)
plotly::ggplotly(p_res_fit, tooltip = "text")p_res_orden <- ggplot(diagnostico, aes(x = corrida, y = residual_estandar,
text = paste0(
"Corrida: ", corrida,
"<br>Residual estándar: ", round(residual_estandar, 4)
))) +
geom_line() +
geom_point(size = 3) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(
title = "Residuales estandarizados en orden de corrida",
x = "Corrida",
y = "Residual estandarizado"
) +
theme_minimal(base_size = 13)
plotly::ggplotly(p_res_orden, tooltip = "text")qq <- tibble::tibble(
teorico = qnorm(ppoints(length(res_std))),
observado = sort(res_std)
)
q_obs <- quantile(qq$observado, probs = c(0.25, 0.75), na.rm = TRUE)
q_teo <- qnorm(c(0.25, 0.75))
pendiente <- as.numeric(diff(q_obs) / diff(q_teo))
intercepto <- as.numeric(q_obs[1] - pendiente * q_teo[1])
x_linea <- range(qq$teorico, na.rm = TRUE)
y_linea <- intercepto + pendiente * x_linea
plotly::plot_ly(qq, x = ~teorico, y = ~observado, type = "scatter", mode = "markers",
marker = list(size = 7),
text = ~paste0("Teórico: ", round(teorico, 4), "<br>Observado: ", round(observado, 4)),
hoverinfo = "text") |>
plotly::add_lines(x = x_linea, y = y_linea, inherit = FALSE,
line = list(dash = "dash"), name = "Línea de referencia") |>
plotly::layout(
title = "Gráfico Q-Q de residuales estandarizados",
xaxis = list(title = "Cuantiles teóricos normales"),
yaxis = list(title = "Cuantiles observados")
)p_hist <- ggplot(diagnostico, aes(x = residual_estandar)) +
geom_histogram(aes(y = after_stat(density)), bins = 8, alpha = 0.75) +
geom_density(linewidth = 1) +
labs(
title = "Distribución de residuales estandarizados",
x = "Residual estandarizado",
y = "Densidad"
) +
theme_minimal(base_size = 13)
plotly::ggplotly(p_hist)p_cook <- ggplot(diagnostico, aes(x = corrida, y = cook,
text = paste0(
"Corrida: ", corrida,
"<br>Cook: ", round(cook, 5),
"<br>Leverage: ", round(leverage, 5)
))) +
geom_col() +
geom_hline(yintercept = 4 / nrow(datos), linetype = "dashed") +
labs(
title = "Distancia de Cook por corrida",
subtitle = "La línea punteada corresponde al criterio 4/n",
x = "Corrida",
y = "Distancia de Cook"
) +
theme_minimal(base_size = 13)
plotly::ggplotly(p_cook, tooltip = "text")Respuesta al literal c). La adecuación del modelo debe concluirse integrando tres evidencias: significancia global del modelo, falta de ajuste no significativa y residuales sin patrones severos. Si el modelo tiene buen R² ajustado, falta de ajuste no significativa, residuales aproximadamente normales, varianza homogénea e independencia, entonces describe adecuadamente la superficie experimental.
Las gráficas de contorno muestran la respuesta predicha para dos
factores, manteniendo fijo el tercer factor. Primero se presentan cortes
en el centro del diseño, es decir, con el tercer factor igual a
0.
bloque_ref <- levels(datos$block)[1]
predecir_y <- function(x) {
nd <- data.frame(x1 = x[1], x2 = x[2], x3 = x[3])
if (usar_bloque_en_modelo) nd$block <- factor(bloque_ref, levels = levels(datos$block))
as.numeric(predict(modelo_lm, newdata = nd))
}
crear_grid <- function(var_x, var_y, fijo = list(), n = 80) {
sec <- seq(-alpha, alpha, length.out = n)
base <- expand.grid(a = sec, b = sec)
nd <- data.frame(
x1 = rep(0, nrow(base)),
x2 = rep(0, nrow(base)),
x3 = rep(0, nrow(base))
)
nd[[var_x]] <- base$a
nd[[var_y]] <- base$b
if (length(fijo) > 0) {
for (nm in names(fijo)) nd[[nm]] <- fijo[[nm]]
}
if (usar_bloque_en_modelo) nd$block <- factor(bloque_ref, levels = levels(datos$block))
base$yhat <- as.numeric(predict(modelo_lm, newdata = nd))
z <- matrix(base$yhat, nrow = n, ncol = n, byrow = FALSE)
list(xs = sec, ys = sec, z = z, data = cbind(nd, yhat = base$yhat))
}
plot_contorno <- function(g, titulo, xlab, ylab) {
plotly::plot_ly(
x = g$xs,
y = g$ys,
z = t(g$z),
type = "contour",
colorscale = escala_viridis_plotly,
contours = list(showlabels = TRUE)
) |>
plotly::layout(
title = titulo,
xaxis = list(title = xlab),
yaxis = list(title = ylab)
)
}
plot_superficie <- function(g, titulo, xlab, ylab) {
plotly::plot_ly(
x = g$xs,
y = g$ys,
z = t(g$z),
type = "surface",
colorscale = escala_viridis_plotly
) |>
plotly::layout(
title = titulo,
scene = list(
xaxis = list(title = xlab),
yaxis = list(title = ylab),
zaxis = list(title = "Y predicha")
)
)
}La parte cuadrática del modelo puede escribirse como:
\[ \hat{Y} = \beta_0 + \mathbf{b}'\mathbf{x} + \mathbf{x}'\mathbf{B}\mathbf{x} \]
Donde:
\[ \mathbf{x} = (x_1, x_2, x_3)' \]
El punto estacionario se calcula con:
\[ \mathbf{x}_s = -\frac{1}{2}\mathbf{B}^{-1}\mathbf{b} \]
Luego se clasifican los valores propios de \(\mathbf{B}\):
b_lineal <- c(
x1 = get_coef("x1"),
x2 = get_coef("x2"),
x3 = get_coef("x3")
)
B <- matrix(0, nrow = 3, ncol = 3)
rownames(B) <- colnames(B) <- c("x1", "x2", "x3")
B[1, 1] <- get_coef("I(x1^2)")
B[2, 2] <- get_coef("I(x2^2)")
B[3, 3] <- get_coef("I(x3^2)")
B[1, 2] <- B[2, 1] <- get_coef("x1:x2") / 2
B[1, 3] <- B[3, 1] <- get_coef("x1:x3") / 2
B[2, 3] <- B[3, 2] <- get_coef("x2:x3") / 2
x_est <- tryCatch(as.numeric(-0.5 * solve(B, b_lineal)), error = function(e) rep(NA_real_, 3))
names(x_est) <- c("x1", "x2", "x3")
y_est <- if (all(is.finite(x_est))) predecir_y(x_est) else NA_real_
val_propios <- tryCatch(eigen(B)$values, error = function(e) rep(NA_real_, 3))
clasificacion <- dplyr::case_when(
all(is.finite(val_propios)) && all(val_propios < 0) ~ "Máximo local",
all(is.finite(val_propios)) && all(val_propios > 0) ~ "Mínimo local",
all(is.finite(val_propios)) ~ "Punto silla",
TRUE ~ "No clasificable"
)
dentro_region <- if (all(is.finite(x_est))) all(abs(x_est) <= alpha + 1e-8) else FALSE
tabla_estacionario <- tibble::tibble(
Medida = c("x1 estacionario", "x2 estacionario", "x3 estacionario", "Y predicha", "Clasificación", "Dentro de [-alpha, alpha]"),
Valor = c(round(x_est[1], 6), round(x_est[2], 6), round(x_est[3], 6), round(y_est, 6), clasificacion, dentro_region)
)
DT::datatable(tabla_estacionario, caption = "Punto estacionario del modelo cuadrático", options = list(dom = "t"), rownames = FALSE)B_tab <- as.data.frame(round(B, 6))
B_tab$Fila <- rownames(B_tab)
B_tab <- B_tab |> dplyr::select(Fila, dplyr::everything())
vp_tab <- tibble::tibble(
Valor_propio = paste0("lambda_", seq_along(val_propios)),
Valor = round(val_propios, 6)
)
DT::datatable(B_tab, caption = "Matriz cuadrática B", options = list(dom = "t"), rownames = FALSE)Aunque exista un punto estacionario, también se calcula el mejor punto dentro de la región explorada por el diseño. Se usa la restricción:
\[ -\alpha \leq x_i \leq \alpha, \quad i = 1,2,3 \]
objetivo <- tolower(params$objetivo)
if (!objetivo %in% c("max", "min")) stop("params$objetivo debe ser 'max' o 'min'.")
fn_optim <- function(z) {
val <- predecir_y(z)
if (objetivo == "max") -val else val
}
# Se usan varios puntos iniciales para reducir el riesgo de depender de un único arranque.
iniciales <- rbind(
c(0, 0, 0),
c(-1, -1, -1), c(1, -1, -1), c(-1, 1, -1), c(1, 1, -1),
c(-1, -1, 1), c(1, -1, 1), c(-1, 1, 1), c(1, 1, 1),
c(-alpha, 0, 0), c(alpha, 0, 0),
c(0, -alpha, 0), c(0, alpha, 0),
c(0, 0, -alpha), c(0, 0, alpha)
)
opts <- purrr::map(seq_len(nrow(iniciales)), function(i) {
optim(
par = iniciales[i, ],
fn = fn_optim,
method = "L-BFGS-B",
lower = rep(-alpha, 3),
upper = rep(alpha, 3)
)
})
y_opts <- purrr::map_dbl(opts, ~ predecir_y(.x$par))
idx_best <- if (objetivo == "max") which.max(y_opts) else which.min(y_opts)
optimo <- opts[[idx_best]]
x_opt <- optimo$par
y_opt <- predecir_y(x_opt)
tabla_opt <- tibble::tibble(
Criterio = c("Objetivo", "x1 óptimo", "x2 óptimo", "x3 óptimo", "Y predicha en el óptimo", "Convergencia optim"),
Valor = c(objetivo, round(x_opt[1], 6), round(x_opt[2], 6), round(x_opt[3], 6), round(y_opt, 6), optimo$convergence)
)
DT::datatable(tabla_opt, caption = "Óptimo acotado dentro de la región experimental", options = list(dom = "t"), rownames = FALSE)opt_df <- tibble::tibble(
x1 = x_opt[1],
x2 = x_opt[2],
x3 = x_opt[3],
y = y_opt,
tipo = "Óptimo acotado"
)
plotly::plot_ly() |>
plotly::add_markers(
data = datos,
x = ~x1, y = ~x2, z = ~x3,
color = ~y,
colors = paleta_viridis,
marker = list(size = 5),
text = ~paste0("Corrida: ", corrida, "<br>Y: ", y),
hoverinfo = "text",
name = "Puntos experimentales"
) |>
plotly::add_markers(
data = opt_df,
x = ~x1, y = ~x2, z = ~x3,
marker = list(size = 9, symbol = "diamond"),
text = ~paste0("Óptimo<br>x1: ", round(x1, 4), "<br>x2: ", round(x2, 4), "<br>x3: ", round(x3, 4), "<br>Y predicha: ", round(y, 4)),
hoverinfo = "text",
name = "Óptimo"
) |>
plotly::layout(
title = "Ubicación del óptimo acotado en el espacio experimental",
scene = list(
xaxis = list(title = "x1"),
yaxis = list(title = "x2"),
zaxis = list(title = "x3")
)
)La técnica de escalonamiento usa el gradiente del modelo para avanzar desde el centro del diseño hacia valores crecientes o decrecientes de la respuesta.
Para maximizar:
\[ \mathbf{d} = \frac{\nabla \hat{Y}}{\|\nabla \hat{Y}\|} \]
Para minimizar:
\[ \mathbf{d} = -\frac{\nabla \hat{Y}}{\|\nabla \hat{Y}\|} \]
En el modelo cuadrático:
\[ \nabla \hat{Y} = \mathbf{b} + 2\mathbf{B}\mathbf{x} \]
gradiente <- function(x) {
as.numeric(b_lineal + 2 * B %*% x)
}
generar_ruta <- function(x0 = c(0, 0, 0), paso = params$tamano_paso, n_pasos = params$n_pasos, objetivo = params$objetivo) {
sentido <- ifelse(tolower(objetivo) == "max", 1, -1)
actual <- as.numeric(x0)
ruta <- tibble::tibble(
paso = 0,
x1 = actual[1],
x2 = actual[2],
x3 = actual[3],
y_pred = predecir_y(actual),
norma_gradiente = sqrt(sum(gradiente(actual)^2))
)
for (i in seq_len(n_pasos)) {
g <- gradiente(actual)
norma <- sqrt(sum(g^2))
if (!is.finite(norma) || norma < 1e-10) break
direccion <- sentido * g / norma
# Cálculo del máximo avance posible antes de salir de [-alpha, alpha].
t_lim <- rep(Inf, 3)
for (j in 1:3) {
if (direccion[j] > 0) t_lim[j] <- (alpha - actual[j]) / direccion[j]
if (direccion[j] < 0) t_lim[j] <- (-alpha - actual[j]) / direccion[j]
}
avance <- min(paso, t_lim, na.rm = TRUE)
if (!is.finite(avance) || avance <= 1e-10) break
siguiente <- actual + avance * direccion
siguiente[abs(siguiente) > alpha] <- sign(siguiente[abs(siguiente) > alpha]) * alpha
ruta <- dplyr::bind_rows(
ruta,
tibble::tibble(
paso = i,
x1 = siguiente[1],
x2 = siguiente[2],
x3 = siguiente[3],
y_pred = predecir_y(siguiente),
norma_gradiente = sqrt(sum(gradiente(siguiente)^2))
)
)
actual <- siguiente
if (any(abs(actual) >= alpha - 1e-8)) break
}
ruta
}
ruta_escalonamiento <- generar_ruta(objetivo = objetivo)
DT::datatable(
ruta_escalonamiento |> dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 6))),
caption = paste0("Ruta de escalonamiento ", ifelse(objetivo == "max", "ascendente", "descendente")),
options = list(pageLength = 20, scrollX = TRUE)
)plotly::plot_ly() |>
plotly::add_markers(
data = datos,
x = ~x1, y = ~x2, z = ~x3,
color = ~y,
colors = paleta_viridis,
marker = list(size = 5),
name = "Puntos experimentales",
text = ~paste0("Corrida: ", corrida, "<br>Y: ", y),
hoverinfo = "text"
) |>
plotly::add_trace(
data = ruta_escalonamiento,
x = ~x1, y = ~x2, z = ~x3,
type = "scatter3d",
mode = "lines+markers",
marker = list(size = 5),
line = list(width = 5),
name = "Ruta de escalonamiento",
text = ~paste0("Paso: ", paso, "<br>Y predicha: ", round(y_pred, 4)),
hoverinfo = "text"
) |>
plotly::layout(
title = paste0("Ruta de escalonamiento ", ifelse(objetivo == "max", "ascendente", "descendente")),
scene = list(
xaxis = list(title = "x1"),
yaxis = list(title = "x2"),
zaxis = list(title = "x3")
)
)p_ruta_y <- ggplot(ruta_escalonamiento, aes(x = paso, y = y_pred,
text = paste0(
"Paso: ", paso,
"<br>x1: ", round(x1, 4),
"<br>x2: ", round(x2, 4),
"<br>x3: ", round(x3, 4),
"<br>Y predicha: ", round(y_pred, 4)
))) +
geom_line() +
geom_point(size = 3) +
labs(
title = paste0("Evolución de Y predicha durante el escalonamiento ", ifelse(objetivo == "max", "ascendente", "descendente")),
x = "Paso",
y = "Y predicha"
) +
theme_minimal(base_size = 13)
plotly::ggplotly(p_ruta_y, tooltip = "text")Estrategia aplicada:
x1, x2 y
x3, con región aproximada [-α, α].b y la matriz cuadrática B;
luego se obtiene xs = -0.5 B^{-1}b.B indican si el punto estacionario es máximo, mínimo o
punto silla.-α ≤ xi ≤ α.(0,0,0) siguiendo el gradiente del modelo para maximizar o
en dirección contraria para minimizar.mejor_ruta <- if (objetivo == "max") {
ruta_escalonamiento[which.max(ruta_escalonamiento$y_pred), ]
} else {
ruta_escalonamiento[which.min(ruta_escalonamiento$y_pred), ]
}
cat("<div class='ok-box'>")Conclusión:
El análisis corresponde a un Diseño Central Compuesto Rotable con tres factores.
cat("El modelo de segundo orden permite estimar la curvatura de la superficie y ubicar un punto candidato de operación. ")El modelo de segundo orden permite estimar la curvatura de la superficie y ubicar un punto candidato de operación.
Punto estacionario: Máximo local , con
cat("x1 = ", round(x_est[1], 4), ", x2 = ", round(x_est[2], 4), ", x3 = ", round(x_est[3], 4), ", Y predicha = ", round(y_est, 4), ". ")x1 = -1.1899 , x2 = 0.3359 , x3 = 0.5166 , Y predicha = 106.0005 .
cat("<b>Óptimo acotado por optimización:</b> x1 = ", round(x_opt[1], 4), ", x2 = ", round(x_opt[2], 4), ", x3 = ", round(x_opt[3], 4), ", Y predicha = ", round(y_opt, 4), ". ")Óptimo acotado por optimización: x1 = -1.1899 , x2 = 0.3359 , x3 = 0.5166 , Y predicha = 106.0005 .
cat("<b>Mejor punto de la ruta de escalonamiento:</b> paso ", mejor_ruta$paso, ", x1 = ", round(mejor_ruta$x1, 4), ", x2 = ", round(mejor_ruta$x2, 4), ", x3 = ", round(mejor_ruta$x3, 4), ", Y predicha = ", round(mejor_ruta$y_pred, 4), ".")Mejor punto de la ruta de escalonamiento: paso 9 , x1 = -1.1789 , x2 = 0.333 , x3 = 0.514 , Y predicha = 106.0001 .
cat("Se recomienda realizar una <b>corrida confirmatoria</b> en el punto óptimo sugerido para verificar experimentalmente la respuesta.")Se recomienda realizar una corrida confirmatoria en el punto óptimo sugerido para verificar experimentalmente la respuesta.
Para maximizar la respuesta, dejar en el encabezado:
Para minimizar la respuesta, cambiar a:
En este archivo los factores se analizan en escala codificada porque
la matriz entregada contiene x1, x2 y
x3. Si se conocen los valores reales de cada variable, se
puede convertir el óptimo usando:
\[ X_i = X_{0i} + \Delta_i x_i \]
Donde \(X_{0i}\) es el valor central real del factor y \(\Delta_i\) es el semirango experimental.