Introducción

La Inferencia estadística de procesos en modelos de superficie de respuesta (RSM, por sus siglas en inglés) de segundo orden se refiere a la aplicación de técnicas estadísticas para analizar, interpretar y realizar conclusiones sobre el comportamiento de un sistema o proceso modelado a través de un polinomio de segundo grado. Los RSM son útiles en la optimización de procesos y diseño de experimentos.

Fundamentos de RSM de Segundo Orden

1. Modelo de Segundo Orden

El modelo general de RSM de segundo orden se expresa como:

\[Y = \beta_0 + \sum_{i=1}^k \beta_i X_i + \sum_{i=1}^k \beta_{ii} X_i^2 + \sum_{i<j} \beta_{ij} X_i X_j + \epsilon,\]

donde:

  • \(Y\): variable respuesta.
  • \(X_i\): variables independientes o factores.
  • \(\beta_0, \beta_i, \beta_{ii}, \beta_{ij}\): coeficientes del modelo.
  • \(\epsilon\): término de error.

2. Componentes Claves

  • Términos lineales (\(\beta_i X_i\)): Efecto principal de los factores.
  • Términos cuadráticos (\(\beta_{ii} X_i^2\)): Efecto de curvatura de cada factor.
  • Términos de interacción (\(\beta_{ij} X_i X_j\)): Efecto combinado de dos factores.

3. Propósito

  • Optimizar \(Y\) encontrando los valores óptimos de \(X_i\).
  • Identificar la forma de la superficie de respuesta (e.g., convexa, cóncava, con máximos o mínimos locales).
  • Comprender las relaciones entre los factores y la respuesta.

Inferencia Estadística en RSM de Segundo Orden

La inferencia estadística en este contexto implica el uso de análisis para evaluar la significancia y la calidad del modelo.

1. Estimación de Parámetros

Los coeficientes (\(\beta\)) del modelo se estiman comúnmente mediante mínimos cuadrados ordinarios (OLS). Esto genera valores estimados para los parámetros basados en los datos experimentales.

2. Pruebas de Hipótesis

  • Significancia de los coeficientes:
    • Se evalúa mediante pruebas \(t\) para cada coeficiente, con la hipótesis nula \(H_0: \beta_i = 0\).
  • Ajuste global del modelo:
    • Se usa la prueba \(F\) para evaluar si el modelo explica significativamente la variabilidad de \(Y\).

3. Análisis de Varianza (ANOVA)

Permite descomponer la variabilidad total en componentes debido al modelo y al error, ayudando a evaluar la calidad del ajuste.

4. Bondad de Ajuste

  • \(R^2\): Proporción de la variabilidad explicada por el modelo.
  • \(R^2\) ajustado: Considera el número de parámetros del modelo y la cantidad de datos.

5. Validación del Modelo

  • Análisis de residuales para verificar supuestos de normalidad, homocedasticidad e independencia.
  • Gráficos de diagnóstico.

6. Optimización

Una vez validado el modelo, se puede buscar:

  • Máximos o mínimos de la superficie.
  • Puntos de estacionariedad analizando las derivadas parciales.

Ejercicios

Ejercicio 1: Optimización de un Proceso de Producción

Una fábrica desea optimizar la calidad de un producto (\(Y\)) en función de dos factores controlables:

El objetivo es ajustar un modelo de segundo orden que relacione \(Y\) con \(X_1\) y \(X_2\) y determinar las condiciones óptimas para maximizar \(Y\).

En lugar de usar los valores tradicionales de \(\pm 1.414\) (que corresponden a \(\sqrt{2}\)), usamos valores distintos para los puntos axiales. Vamos a usar los puntos axiales a \(\pm 2\), lo que representa una mayor distancia respecto al centro, para captar la curvatura de la superficie de respuesta en una escala diferente.

Efectos principales

Corresponde a los puntos experimentales con \(X_1\) y \(X_2\) a los niveles -1, 0 y 1.

Réplicas del centro

Se añaden más puntos en el centro de la matriz experimental, con \(X_1 = 0\) y \(X_2 = 0\), y repeticiones para evaluar la variabilidad.

Puntos axiales

Los puntos en los que \(X_1\) y \(X_2\) toman valores más alejados del centro, específicamente a \(\pm 2\).

Datos Experimentales

X1 X2 Y
-1 -1 4.749
1 -1 6.359
-1 1 1.194
1 1 5.657
-1.414 0 2.408
1.414 0 6.650
0 -1.414 6.170
0 1.414 2.533
0 0 2.859
0 0 3.163
0 0 2.861
0 0 2.860

Pasos para Resolver en Excel

1. Ajustar el Modelo Cuadrático

  • Agregar columnas en la tabla para los términos necesarios del modelo:
    • \(X_1^2\): Cuadrado de \(X_1\).
    • \(X_2^2\): Cuadrado de \(X_2\).
    • \(X_1 \times X_2\): Producto de \(X_1\) y \(X_2\).

2. Usar Regresión Múltiple

  • Seleccionar la herramienta Análisis de Datos en Excel.
  • Elegir Regresión y configurar:
    • Rango de entrada Y: Seleccionar la columna \(Y\).
    • Rango de entrada X: Seleccionar las columnas \(X_1\), \(X_2\), \(X_1^2\), \(X_2^2\), \(X_1 \times X_2\).
    • Marcar la casilla de Rótulos si has incluido encabezados en las selecciones.
  • Generar el análisis.

3. Interpretación de Resultados

  • Coeficientes estimados (\(\beta_0, \beta_1, \beta_2, \beta_{ii}, \beta_{12}\)): Observar los valores en la tabla de coeficientes.
  • Significancia estadística: Evaluar los valores \(p\) (valores menores a 0.05 indican significancia).
  • Calidad del modelo:
    • Examinar \(R^2\) y \(R^2\) ajustado para evaluar la proporción de variabilidad explicada por el modelo.

4. Predecir y Optimizar

  • Utilizar los coeficientes ajustados para predecir \(Y\) mediante la fórmula:

\[ Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_{11} X_1^2 + \beta_{22} X_2^2 + \beta_{12} X_1 X_2 \]

  • Encontrar los valores de \(X_1\) y \(X_2\) que maximizan \(Y\).

Resolución en R

Cargar los datos

data <- data.frame(
  X1 = c(-1, 1, -1, 1, -1.414, 1.414, 0, 0, 0, 0, 0, 0),
  X2 = c(-1, -1, 1, 1, 0, 0, -1.414, 1.414, 0, 0, 0, 0),
  Y = c(4.749, 6.359, 1.194, 5.657, 2.408, 6.65, 6.17, 2.533, 2.859, 3.163, 2.861, 2.86)
)

Ajustar el modelo cuadrático

model <- lm(Y ~ X1 + X2 + I(X1^2) + I(X2^2) + X1:X2, data = data)

# Ver el resumen del modelo
summary(model)
## 
## Call:
## lm(formula = Y ~ X1 + X2 + I(X1^2) + I(X2^2) + X1:X2, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.18137 -0.07687 -0.05609  0.12779  0.22726 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.93574    0.08483  34.605 3.88e-08 ***
## X1           1.50913    0.05999  25.156 2.60e-07 ***
## X2          -1.17514    0.05999 -19.588 1.15e-06 ***
## I(X1^2)      0.80913    0.06708  12.062 1.97e-05 ***
## I(X2^2)      0.72036    0.06708  10.738 3.85e-05 ***
## X1:X2        0.71325    0.08483   8.408 0.000154 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1697 on 6 degrees of freedom
## Multiple R-squared:  0.9954, Adjusted R-squared:  0.9916 
## F-statistic:   261 on 5 and 6 DF,  p-value: 6.26e-07

Interpretación:

El modelo es adecuado, con un excelente ajuste (alto R^2) y es estadísticamente significativo (p-value muy bajo). Esto sugiere que el modelo de segundo orden está describiendo correctamente la relación entre las variables predictoras y la respuesta.

1. Prueba F para validar el modelo

H0: El modelo es insignificante (los efectos no tienen impacto sobre Y)

anova(model)

Interpretación:

El modelo muestra que tanto la Temperatura (X1) como el Tiempo de mezclado (X2) tienen efectos significativos en la calidad del producto, con un ajuste sólido del modelo.

2. Prueba t para validar los coeficientes (validez de los betas)

Para cada coeficiente, vemos si podemos rechazar H0: beta = 0 (los coeficientes no son significativos)

summary(model)
## 
## Call:
## lm(formula = Y ~ X1 + X2 + I(X1^2) + I(X2^2) + X1:X2, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.18137 -0.07687 -0.05609  0.12779  0.22726 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.93574    0.08483  34.605 3.88e-08 ***
## X1           1.50913    0.05999  25.156 2.60e-07 ***
## X2          -1.17514    0.05999 -19.588 1.15e-06 ***
## I(X1^2)      0.80913    0.06708  12.062 1.97e-05 ***
## I(X2^2)      0.72036    0.06708  10.738 3.85e-05 ***
## X1:X2        0.71325    0.08483   8.408 0.000154 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1697 on 6 degrees of freedom
## Multiple R-squared:  0.9954, Adjusted R-squared:  0.9916 
## F-statistic:   261 on 5 and 6 DF,  p-value: 6.26e-07

Interpretación:

La Temperatura (X1) y el Tiempo de mezclado (X2) tienen un impacto significativo en la calidad del producto. Además, ambos factores muestran efectos no lineales.

3. Prueba de la falta de ajuste (Lack of Fit)

Comprobamos si el modelo cuadrático se ajusta bien a los datos

model_complete <- lm(Y ~ poly(X1, 2) + poly(X2, 2), data = data)

# Realizar la comparación de modelos usando ANOVA
anova(model, model_complete)

Interpretación:

El modelo 1 (con los términos lineales, cuadráticos y de interacción) ya es un buen ajuste para los datos. El modelo 2 (polinómico) no ofrece una mejora sustancial, ya que la prueba de “Lack of Fit” muestra que no hay evidencia suficiente para sugerir que el modelo 1 no esté capturando toda la variabilidad de \(Y\).

El valor \(p\) bajo (menor que 0.05) en la prueba indica que el modelo más flexible (modelo 2) no mejora significativamente el ajuste. Por lo tanto, podemos concluir que el modelo 1 es adecuado y no presenta una falta de ajuste significativa.

4. Prueba de curvatura

Ho No hay curvatura significativa en el modelo
H1 Hay curvatura significativa en el modelo

La prueba de curvatura evalúa si el modelo cuadrático tiene una curvatura significativa. Para esto, podemos usar los términos cuadráticos e interactivos en el modelo y observar si son significativos.

# Modelo lineal (solo términos lineales)
model_linear <- lm(Y ~ X1 + X2, data = data)

# Modelo no lineal (con términos cuadráticos y de interacción)
model_nonlinear <- lm(Y ~ X1 + X2 + I(X1^2) + I(X2^2) + X1:X2, data = data)

# Comparar los dos modelos usando ANOVA
curvature_test <- anova(model_linear, model_nonlinear)
curvature_test

Interpretación

El modelo con términos cuadráticos y de interacción es altamente significativo y proporciona un ajuste excelente a los datos. La curvatura está presente y es importante para modelar correctamente la relación entre \(X_1\), \(X_2\) y \(Y\).

El modelo lineal por sí solo no es suficiente para capturar toda la variabilidad de \(Y\), como lo demuestra la prueba de curvatura.

5. Graficar la superficie de respuesta con plotly

# Generar una cuadrícula para predecir
x1_seq <- seq(-2.5, 2.5, length.out = 50)
x2_seq <- seq(-2.5, 2.5, length.out = 50)
grid <- expand.grid(X1 = x1_seq, X2 = x2_seq)
grid$Y <- predict(model, newdata = grid)

Graficar la superficie de respuesta

library(plotly)
## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
plot_ly(
  data = grid,
  x = ~X1,
  y = ~X2,
  z = ~Y,
  type = "mesh3d"
) %>%
  layout(
    scene = list(
      xaxis = list(title = "X1 (Temperatura)"),
      yaxis = list(title = "X2 (Tiempo de Mezclado)"),
      zaxis = list(title = "Calidad (Y)")
    )
  )
# Cargar librerías necesarias
library(ggplot2)

# Definir el modelo (suponiendo que ya lo has ajustado)
model <- lm(Y ~ X1 + X2 + I(X1^2) + I(X2^2) + X1:X2, data = data)

# Crear una malla de valores para X1 y X2 con mayor resolución
x1_seq <- seq(-2, 2, length.out = 200)  # Aumentamos la resolución
x2_seq <- seq(-2, 2, length.out = 200)  # Aumentamos la resolución
grid <- expand.grid(X1 = x1_seq, X2 = x2_seq)  # Generar la cuadrícula

# Predecir valores de Y para cada punto en la malla usando el modelo
grid$Y <- predict(model, newdata = grid)

# Encontrar el punto óptimo utilizando optimización
opt <- optim(
  par = c(0, 0),  # Valores iniciales (0, 0)
  fn = function(x) {
    -predict(model, newdata = data.frame(X1 = x[1], X2 = x[2]))  # Negamos para maximizar
  },
  control = list(fnscale = -1)  # Maximizar la función
)

# Los valores óptimos
opt_x1 <- opt$par[1]
opt_x2 <- opt$par[2]
opt_Y <- -opt$value  # El valor máximo de Y

# Graficar las curvas de nivel y añadir el punto óptimo
ggplot(grid, aes(x = X1, y = X2, z = Y)) +
  geom_contour_filled(aes(fill = ..level..), color = "white") +  # Curvas de nivel rellenas
  geom_point(aes(x = opt_x1, y = opt_x2), color = "red", size = 4, shape = 17) +  # Añadir punto óptimo
  labs(title = "Curvas de Nivel con Punto Óptimo",
       x = "Temperatura (X1)", 
       y = "Tiempo de Mezclado (X2)", 
       fill = "Calidad (Y)") +
  theme_minimal() +
  annotate("text", x = opt_x1, y = opt_x2, label = paste("Óptimo\n(", round(opt_x1, 2), ",", round(opt_x2, 2), ")", sep = ""), 
           color = "red", size = 5, vjust = -1, hjust = 0.5) +  # Añadir texto con la ubicación del óptimo
  coord_cartesian(xlim = c(-2, 2), ylim = c(-2, 2)) +  # Ajustar los límites de los ejes
  theme(plot.margin = margin(1, 1, 1, 1, "cm"))  # Ampliar el margen de la gráfica
## Warning: The dot-dot notation (`..level..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(level)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

6. Optimización para maximizar Y

Optimización para encontrar el valor óptimo de X1 y X2 que maximicen Y

Definimos la función del modelo

modelo_funcion <- function(par) {
  X1 <- par[1]  # Primer parámetro (X1)
  X2 <- par[2]  # Segundo parámetro (X2)
  
  # Ecuación del modelo
  Y <- 2.93574 + 1.50913*X1 - 1.17514*X2 + 0.80913*X1^2 + 0.72036*X2^2 + 0.71325*X1*X2
  
  return(Y)  # Retornamos el valor de Y
}

Resultado óptimo

resultado_optimo <- optim(c(0, 0), modelo_funcion, method = "BFGS")

# Resultado óptimo
X1_optimo <- resultado_optimo$par[1]  # Valor óptimo de X1
X2_optimo <- resultado_optimo$par[2]  # Valor óptimo de X2
Y_optimo <- modelo_funcion(resultado_optimo$par)  # Valor de Y en el óptimo

Resultados óptimos

cat("El valor óptimo de X1 es:", X1_optimo, "\n")
## El valor óptimo de X1 es: -1.652684
cat("El valor óptimo de X2 es:", X2_optimo, "\n")
## El valor óptimo de X2 es: 1.633851
cat("El valor mínimo de Y en ese punto es:", Y_optimo, "\n")
## El valor mínimo de Y en ese punto es: 0.7286832
library(numDeriv)
# Calcular la matriz Hessiana en el punto óptimo (segunda derivada)
hessiana <- hessian(func = modelo_funcion, x = resultado_optimo$par)

# Imprimir la matriz Hessiana
cat("\nLa matriz Hessiana en el punto óptimo es:\n")
## 
## La matriz Hessiana en el punto óptimo es:
print(hessiana)
##         [,1]    [,2]
## [1,] 1.61826 0.71325
## [2,] 0.71325 1.44072
# Verificar la naturaleza de la matriz Hessiana (positiva definida o negativa definida)
eigenvalues <- eigen(hessiana)$values

# Ver los valores propios de la matriz Hessiana
cat("\nValores propios de la matriz Hessiana:", eigenvalues, "\n")
## 
## Valores propios de la matriz Hessiana: 2.248243 0.8107371
# Caracterización del punto óptimo
if (all(eigenvalues > 0)) {
  cat("\nEl punto es un mínimo local (matriz Hessiana positiva definida).")
} else if (all(eigenvalues < 0)) {
  cat("\nEl punto es un máximo local (matriz Hessiana negativa definida).")
} else {
  cat("\nEl punto es un punto silla (matriz Hessiana con valores propios de signo mixto).")
}
## 
## El punto es un mínimo local (matriz Hessiana positiva definida).

Intervalos de confianza para los parámetros

confint(model)
##                  2.5 %     97.5 %
## (Intercept)  2.7281599  3.1433253
## X1           1.3623321  1.6559206
## X2          -1.3219364 -1.0283480
## I(X1^2)      0.6449912  0.9732777
## I(X2^2)      0.5562144  0.8845009
## X1:X2        0.5056673  0.9208327

El coeficiente de \(X_1\) está en el intervalo \([1.489, 1.529]\), lo que significa que podemos estar un 95% seguros de que el valor real del parámetro \(\beta_1\) está dentro de ese rango.

Propagación de incertidumbre al punto óptimo

Método de Monte Carlo: Este es el más directo, donde se simula los coeficientes dentro de los intervalos de confianza y luego se calcula el punto óptimo para cada simulación.

Monte Carlo

modelo_funcion <- function(X1, X2, coef_simulados) {
  # coef_simulados es un vector con los coeficientes simulados
  beta_0 <- coef_simulados[1]
  beta_1 <- coef_simulados[2]
  beta_2 <- coef_simulados[3]
  beta_3 <- coef_simulados[4]
  beta_4 <- coef_simulados[5]
  beta_5 <- coef_simulados[6]
  
  # Calculamos el valor de Y en función de X1 y X2
  Y <- beta_0 + beta_1 * X1 + beta_2 * X2 + beta_3 * X1^2 + beta_4 * X2^2 + beta_5 * X1 * X2
  
  return(Y)
}
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:plotly':
## 
##     select
# Número de simulaciones
n_sim <- 1000

# Simulación de los coeficientes dentro de sus intervalos de confianza
simulaciones <- replicate(n_sim, {
  # Generar un conjunto de coeficientes aleatorios dentro de los intervalos de confianza
  coef_simulados <- mvrnorm(1, coef(model), vcov(model))
  
  # Encontrar el punto óptimo con los coeficientes simulados
  optim_result <- optim(c(0, 0), function(par) {
    # Aquí, par[1] es X1 y par[2] es X2
    modelo_funcion(par[1], par[2], coef_simulados)
  })
  
  # Guardar el punto óptimo encontrado
  optim_result$par
})

# Calcular el intervalo de confianza para el punto óptimo
IC_optimo <- apply(simulaciones, 1, function(x) quantile(x, c(0.025, 0.975)))

# Mostrar los intervalos de confianza para el punto óptimo
IC_optimo
##            [,1]     [,2]
## 2.5%  -2.411839 1.183632
## 97.5% -1.245943 2.458000

Interpetación:

El intervalo de confianza del 95% para el valor óptimo del tiempo de mezclado es entre 1.227 y 2.545.

Esto significa que, en el 95% de las simulaciones, el valor óptimo del tiempo de mezclado se encuentra dentro de este rango. Es un intervalo positivo, lo cual tiene más sentido en un contexto práctico.

Conclusion

Tras realizar la optimización del modelo, el punto óptimo de \(X_1 = X_1^*\) y \(X_2 = X_2^*\) se encuentra en el valor de \(Y = Y^*\). La matriz Hessiana evaluada en este punto tiene todos los valores propios positivos, lo que indica que el punto es un mínimo local. La gráfica 3D muestra que este punto es efectivamente el más bajo en el espacio de las variables \(X_1\) y \(X_2\).

Por lo tanto, podemos concluir que hemos encontrado el punto óptimo para optimizar la calidad del producto en función de la temperatura y el tiempo de mezclado.

Referencias