Resumen ejecutivo

El presente documento desarrolla un ejercicio demostrativo de Diseño de Experimentos (DOE) y Metodología de Superficie de Respuesta (RSM) aplicado a nutrición de camarón blanco (Litopenaeus vannamei). El propósito es evaluar el efecto conjunto de dos factores nutricionales —porcentaje de proteína cruda y porcentaje de lípidos— sobre indicadores productivos relevantes: ganancia de peso, peso final, factor de conversión alimenticia (FCR), supervivencia y tasa específica de crecimiento (SGR).

El estudio se construye con datos simulados exclusivamente con fines académicos y de portafolio técnico. Por tanto, las conclusiones ilustran el procedimiento estadístico y no deben interpretarse como recomendaciones nutricionales comerciales ni como resultados de un ensayo de campo real.

La lógica del análisis es:

diseño experimental → control de calidad de datos → análisis descriptivo → modelado cuadrático → evaluación de falta de ajuste → diagnóstico → superficie de respuesta → optimización individual y multirrespuesta → propuesta de ensayo confirmatorio.

1. Planteamiento del ensayo

1.1 Pregunta de investigación

¿Qué combinación de proteína cruda y lípidos en una dieta experimental permite mejorar el desempeño productivo del camarón, maximizando crecimiento y supervivencia y reduciendo el factor de conversión alimenticia?

1.2 Factores experimentales

Se consideran dos factores controlables:

Factor Nivel bajo (-1) Nivel central (0) Nivel alto (+1)
Proteína cruda (%) 32 36 40
Lípidos (%) 6 8 10

La combinación de tres niveles de proteína y tres niveles de lípidos genera:

\[ 3 \times 3 = 9 \text{ tratamientos} \]

Cada tratamiento dispone de tres réplicas independientes:

\[ 9 \times 3 = 27 \text{ unidades experimentales} \]

En un ensayo real, la unidad experimental sería el tanque, no cada camarón individual. Esta distinción evita pseudorreplicación, debido a que los animales de un mismo tanque comparten dieta, agua, manejo y condiciones ambientales.

1.3 Variables de respuesta

Las principales variables productivas son:

  • Ganancia de peso (g): peso final menos peso inicial.
  • Peso final (g).
  • FCR: alimento consumido / ganancia de biomasa. Valores menores indican mayor eficiencia alimenticia.
  • Supervivencia (%).
  • SGR (%/día): tasa específica de crecimiento.

\[ SGR=100\frac{\ln(W_f)-\ln(W_i)}{t} \]

donde \(W_f\) es el peso final, \(W_i\) el peso inicial y \(t\) el número de días del ensayo.

2. Importación y control de calidad de los datos

archivos <- list.files(
  pattern = "^datos_simulados_rsm_camaron.*\\.csv$",
  ignore.case = TRUE
)

if (length(archivos) == 0) {
  stop("No se encontró el archivo CSV. Coloque el R Markdown y el CSV en la misma carpeta.")
}

# Si existen varias versiones, se utiliza la modificada más recientemente.
info_archivos <- file.info(archivos)
archivo <- rownames(info_archivos)[which.max(info_archivos$mtime)]

datos <- readr::read_delim(
  archivo,
  delim = ";",
  locale = locale(decimal_mark = "."),
  show_col_types = FALSE,
  trim_ws = TRUE
)

cat("Archivo utilizado:", archivo, "\n")
## Archivo utilizado: datos_simulados_rsm_camaron.csv
control <- tibble(
  Indicador = c(
    "Unidades experimentales",
    "Variables",
    "Valores faltantes",
    "Filas duplicadas",
    "Niveles de proteína",
    "Niveles de lípidos"
  ),
  Valor = c(
    nrow(datos),
    ncol(datos),
    sum(is.na(datos)),
    sum(duplicated(datos)),
    n_distinct(datos$proteina_pct),
    n_distinct(datos$lipidos_pct)
  )
)

kable(
  control,
  caption = "Control de calidad y estructura general de la base"
)
Control de calidad y estructura general de la base
Indicador Valor
Unidades experimentales 27
Variables 12
Valores faltantes 0
Filas duplicadas 0
Niveles de proteína 3
Niveles de lípidos 3

La base contiene 27 unidades experimentales y 12 variables. Se identifican 3 niveles de proteína y 3 niveles de lípidos, correspondientes a un diseño factorial \(3\times3\) con replicación.

glimpse(datos)
## Rows: 27
## Columns: 12
## $ orden_corrida     <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1…
## $ replica           <dbl> 2, 3, 3, 1, 1, 1, 2, 1, 3, 2, 3, 3, 2, 2, 1, 1, 1, 2…
## $ proteina_pct      <dbl> 40, 32, 40, 32, 36, 32, 40, 40, 36, 32, 32, 32, 36, …
## $ lipidos_pct       <dbl> 8, 10, 6, 8, 8, 6, 6, 6, 8, 10, 8, 6, 10, 10, 10, 8,…
## $ x1                <dbl> 1, -1, 1, -1, 0, -1, 1, 1, 0, -1, -1, -1, 0, 1, -1, …
## $ x2                <dbl> 0, 1, -1, 0, 0, -1, -1, -1, 0, 1, 0, -1, 1, 1, 1, 0,…
## $ peso_inicial_g    <dbl> 1.935, 1.961, 2.027, 1.964, 1.977, 2.029, 2.046, 1.9…
## $ ganancia_peso_g   <dbl> 18.287, 15.546, 15.533, 15.318, 18.525, 13.680, 16.0…
## $ peso_final_g      <dbl> 20.223, 17.507, 17.560, 17.282, 20.502, 15.709, 18.1…
## $ FCR               <dbl> 1.281, 1.494, 1.437, 1.485, 1.321, 1.591, 1.460, 1.4…
## $ supervivencia_pct <dbl> 93.09, 93.16, 93.89, 94.18, 95.26, 92.53, 92.23, 95.…
## $ SGR_pct_dia       <dbl> 4.190, 3.910, 3.855, 3.884, 4.177, 3.655, 3.894, 3.8…

2.1 Verificación de réplicas

tabla_replicas <- datos %>%
  count(proteina_pct, lipidos_pct, name = "numero_replicas") %>%
  arrange(proteina_pct, lipidos_pct)

kable(
  tabla_replicas,
  caption = "Número de réplicas por combinación experimental"
)
Número de réplicas por combinación experimental
proteina_pct lipidos_pct numero_replicas
32 6 3
32 8 3
32 10 3
36 6 3
36 8 3
36 10 3
40 6 3
40 8 3
40 10 3

Un diseño balanceado debería presentar el mismo número de réplicas en cada combinación. En esta base, el número mínimo de réplicas por tratamiento es 3 y el máximo es 3, confirmándose el balance del diseño.

2.2 Codificación de factores

Para el análisis de superficie de respuesta se emplean variables codificadas:

\[ x_1=\frac{\text{Proteína}-36}{4} \]

\[ x_2=\frac{\text{Lípidos}-8}{2} \]

tabla_codificacion <- datos %>%
  distinct(proteina_pct, lipidos_pct, x1, x2) %>%
  arrange(proteina_pct, lipidos_pct)

kable(
  tabla_codificacion,
  digits = 2,
  caption = "Correspondencia entre niveles reales y variables codificadas"
)
Correspondencia entre niveles reales y variables codificadas
proteina_pct lipidos_pct x1 x2
32 6 -1 -1
32 8 -1 0
32 10 -1 1
36 6 0 -1
36 8 0 0
36 10 0 1
40 6 1 -1
40 8 1 0
40 10 1 1

La codificación centra el experimento en \(x_1=0\) y \(x_2=0\), equivalentes a 36% de proteína y 8% de lípidos. Los niveles extremos se representan mediante -1 y +1, facilitando la interpretación de efectos lineales, interacción y curvatura.

3. Estadística descriptiva por tratamiento

resumen_tratamientos <- datos %>%
  group_by(proteina_pct, lipidos_pct) %>%
  summarise(
    n = n(),
    peso_inicial_medio = mean(peso_inicial_g),
    peso_final_medio = mean(peso_final_g),
    ganancia_media = mean(ganancia_peso_g),
    sd_ganancia = sd(ganancia_peso_g),
    FCR_medio = mean(FCR),
    sd_FCR = sd(FCR),
    supervivencia_media = mean(supervivencia_pct),
    SGR_medio = mean(SGR_pct_dia),
    .groups = "drop"
  )

kable(
  resumen_tratamientos,
  digits = 3,
  caption = "Indicadores productivos medios por tratamiento"
)
Indicadores productivos medios por tratamiento
proteina_pct lipidos_pct n peso_inicial_medio peso_final_medio ganancia_media sd_ganancia FCR_medio sd_FCR supervivencia_media SGR_medio
32 6 3 2.001 15.804 13.803 0.139 1.597 0.006 92.467 3.691
32 8 3 1.977 17.237 15.260 0.341 1.472 0.012 93.823 3.868
32 10 3 2.014 17.073 15.059 0.440 1.482 0.013 93.373 3.818
36 6 3 2.026 18.630 16.604 0.489 1.442 0.009 94.530 3.962
36 8 3 1.951 20.232 18.282 0.224 1.290 0.027 95.003 4.177
36 10 3 1.984 20.054 18.070 0.310 1.307 0.013 94.790 4.131
40 6 3 2.003 17.511 15.508 0.574 1.457 0.018 93.723 3.871
40 8 3 1.964 20.063 18.098 0.248 1.287 0.013 94.350 4.150
40 10 3 1.996 20.821 18.825 0.298 1.303 0.019 94.853 4.187

Interpretación descriptiva. La mayor ganancia media observada se presenta con 40% de proteína y 10% de lípidos, con aproximadamente 18.825 g de ganancia. El menor FCR observado corresponde a 40% de proteína y 8% de lípidos, con un FCR medio de 1.287. La mayor supervivencia media se observa en 36% de proteína y 8% de lípidos, con 95.00%.

Estos resultados descriptivos son útiles, pero todavía no constituyen una optimización. La RSM permitirá estimar qué ocurre entre los niveles ensayados y localizar una región óptima continua.

4. Exploración gráfica

4.1 Ganancia de peso e interacción proteína-lípidos

ggplot(
  resumen_tratamientos,
  aes(
    x = proteina_pct,
    y = ganancia_media,
    group = factor(lipidos_pct),
    linetype = factor(lipidos_pct)
  )
) +
  geom_line(linewidth = 0.9) +
  geom_point(size = 3) +
  scale_x_continuous(breaks = c(32, 36, 40)) +
  labs(
    title = "Interacción entre proteína y lípidos sobre la ganancia de peso",
    x = "Proteína cruda (%)",
    y = "Ganancia media de peso (g)",
    linetype = "Lípidos (%)"
  ) +
  theme_minimal(base_size = 12)

Interpretación. El gráfico permite observar una respuesta no estrictamente lineal. El crecimiento mejora al pasar desde los niveles bajos hacia la región central-alta, pero la magnitud del efecto de la proteína cambia según el nivel de lípidos. La falta de paralelismo perfecto entre las líneas anticipa una posible interacción proteína × lípidos, que será evaluada formalmente mediante el modelo cuadrático.

4.2 Ganancia media con variabilidad experimental

ggplot(
  resumen_tratamientos,
  aes(
    x = factor(proteina_pct),
    y = ganancia_media,
    group = factor(lipidos_pct),
    fill = factor(lipidos_pct)
  )
) +
  geom_col(
    position = position_dodge(width = 0.8),
    width = 0.7
  ) +
  geom_errorbar(
    aes(
      ymin = ganancia_media - sd_ganancia,
      ymax = ganancia_media + sd_ganancia
    ),
    position = position_dodge(width = 0.8),
    width = 0.2
  ) +
  labs(
    title = "Ganancia de peso por combinación nutricional",
    subtitle = "Media ± desviación estándar",
    x = "Proteína cruda (%)",
    y = "Ganancia de peso (g)",
    fill = "Lípidos (%)"
  ) +
  theme_minimal(base_size = 12)

4.3 Factor de conversión alimenticia

ggplot(
  resumen_tratamientos,
  aes(
    x = factor(proteina_pct),
    y = FCR_medio,
    group = factor(lipidos_pct),
    fill = factor(lipidos_pct)
  )
) +
  geom_col(
    position = position_dodge(width = 0.8),
    width = 0.7
  ) +
  labs(
    title = "Factor de conversión alimenticia por tratamiento",
    subtitle = "Valores menores representan mayor eficiencia",
    x = "Proteína cruda (%)",
    y = "FCR medio",
    fill = "Lípidos (%)"
  ) +
  theme_minimal(base_size = 12)

Interpretación. Los tratamientos de la región central-alta muestran una conversión alimenticia más eficiente que la combinación de niveles nutricionales bajos. Sin embargo, elevar simultáneamente ambos nutrientes hasta sus máximos no implica necesariamente el FCR mínimo, lo que justifica buscar un punto de optimización mediante una función cuadrática.

5. Modelo de Superficie de Respuesta

Para cada respuesta se ajusta el modelo cuadrático de segundo orden:

\[ Y=\beta_0+\beta_1x_1+\beta_2x_2+\beta_{12}x_1x_2+ \beta_{11}x_1^2+\beta_{22}x_2^2+\varepsilon \]

donde:

modelo_peso <- lm(
  ganancia_peso_g ~ x1 + x2 + x1:x2 + I(x1^2) + I(x2^2),
  data = datos
)

modelo_fcr <- lm(
  FCR ~ x1 + x2 + x1:x2 + I(x1^2) + I(x2^2),
  data = datos
)

modelo_surv <- lm(
  supervivencia_pct ~ x1 + x2 + x1:x2 + I(x1^2) + I(x2^2),
  data = datos
)

modelo_sgr <- lm(
  SGR_pct_dia ~ x1 + x2 + x1:x2 + I(x1^2) + I(x2^2),
  data = datos
)

# Modelo equivalente con el paquete rsm, útil para análisis canónico.
modelo_peso_rsm <- rsm(
  ganancia_peso_g ~ SO(x1, x2),
  data = datos
)

5.1 Funciones auxiliares para métricas

p_global_modelo <- function(modelo) {
  f <- summary(modelo)$fstatistic
  unname(pf(f[1], f[2], f[3], lower.tail = FALSE))
}

metricas_modelo <- function(modelo, respuesta) {
  s <- summary(modelo)
  tibble(
    Respuesta = respuesta,
    R2 = s$r.squared,
    R2_ajustado = s$adj.r.squared,
    F_global = unname(s$fstatistic[1]),
    p_global = p_global_modelo(modelo)
  )
}

tabla_metricas <- bind_rows(
  metricas_modelo(modelo_peso, "Ganancia de peso"),
  metricas_modelo(modelo_fcr, "FCR"),
  metricas_modelo(modelo_surv, "Supervivencia"),
  metricas_modelo(modelo_sgr, "SGR")
)

kable(
  tabla_metricas,
  digits = 4,
  caption = "Desempeño global de los modelos cuadráticos"
)
Desempeño global de los modelos cuadráticos
Respuesta R2 R2_ajustado F_global p_global
Ganancia de peso 0.9608 0.9514 102.8345 0.0000
FCR 0.9834 0.9794 248.8329 0.0000
Supervivencia 0.5089 0.3919 4.3515 0.0071
SGR 0.9338 0.9180 59.2363 0.0000

5.2 Modelo para ganancia de peso

coef_peso <- as.data.frame(summary(modelo_peso)$coefficients)
coef_peso$Termino <- rownames(coef_peso)
rownames(coef_peso) <- NULL
coef_peso <- coef_peso %>%
  select(Termino, everything())

kable(
  coef_peso,
  digits = 5,
  caption = "Coeficientes del modelo cuadrático para ganancia de peso"
)
Coeficientes del modelo cuadrático para ganancia de peso
Termino Estimate Std. Error t value Pr(>|t|)
(Intercept) 18.25322 0.16481 110.75110 0.00000
x1 1.38478 0.09027 15.34009 0.00000
x2 1.00661 0.09027 11.15089 0.00000
I(x1^2) -1.55967 0.15636 -9.97514 0.00000
I(x2^2) -0.90183 0.15636 -5.76784 0.00001
x1:x2 0.51550 0.11056 4.66263 0.00013 
anova_peso <- as.data.frame(anova(modelo_peso))
anova_peso$Fuente <- rownames(anova_peso)
rownames(anova_peso) <- NULL
anova_peso <- anova_peso %>% select(Fuente, everything())

kable(
  anova_peso,
  digits = 5,
  caption = "ANOVA del modelo de ganancia de peso"
)
ANOVA del modelo de ganancia de peso
Fuente Df Sum Sq Mean Sq F value Pr(>F)
x1 1 34.51697 34.51697 235.31844 0.00000
x2 1 18.23879 18.23879 124.34240 0.00000
I(x1^2) 1 14.59536 14.59536 99.50345 0.00000
I(x2^2) 1 4.87982 4.87982 33.26803 0.00001
x1:x2 1 3.18888 3.18888 21.74012 0.00013
Residuals 21 3.08032 0.14668 NA NA

Interpretación técnica. El modelo explica aproximadamente 96.1% de la variabilidad observada en la ganancia de peso, con un \(R^2\) ajustado de 0.951. La prueba global del modelo presenta \(p=\) <0.001, lo que evidencia una relación estadísticamente relevante entre la formulación nutricional y el crecimiento simulado.

Los efectos lineales de proteína (\(p=\) <0.001) y lípidos (\(p=\) <0.001) son significativos. La interacción proteína × lípidos también resulta significativa (\(p=\) <0.001), indicando que el efecto de un nutriente depende del nivel del otro. Los términos cuadráticos de proteína y lípidos presentan evidencia de curvatura, condición fundamental para que la RSM pueda identificar una región óptima interior.

5.3 Modelo para FCR

coef_fcr <- as.data.frame(summary(modelo_fcr)$coefficients)
coef_fcr$Termino <- rownames(coef_fcr)
rownames(coef_fcr) <- NULL
coef_fcr <- coef_fcr %>%
  select(Termino, everything())

kable(
  coef_fcr,
  digits = 5,
  caption = "Coeficientes del modelo cuadrático para FCR"
)
Coeficientes del modelo cuadrático para FCR
Termino Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.29204 0.00663 194.94709 0.0000
x1 -0.08422 0.00363 -23.20106 0.0000
x2 -0.06728 0.00363 -18.53330 0.0000
I(x1^2) 0.08644 0.00629 13.74857 0.0000
I(x2^2) 0.08161 0.00629 12.97985 0.0000
x1:x2 -0.00983 0.00445 -2.21175 0.0382

Interpretación técnica. El modelo de FCR alcanza un \(R^2\) de 0.983 y un \(R^2\) ajustado de 0.979, mostrando un ajuste muy alto para estos datos simulados. Los coeficientes lineales negativos indican que, dentro de la región inicial del experimento, el incremento de proteína y lípidos reduce el FCR. Los términos cuadráticos positivos muestran posteriormente un cambio de pendiente, generando un mínimo en lugar de una reducción indefinida. La interacción presenta \(p=\) 0.0382.

Desde una perspectiva productiva, este resultado es relevante porque el punto de máximo crecimiento no necesariamente coincide exactamente con el punto de mínimo consumo de alimento por unidad de biomasa producida.

5.4 Modelos para supervivencia y SGR

metricas_secundarias <- bind_rows(
  metricas_modelo(modelo_surv, "Supervivencia"),
  metricas_modelo(modelo_sgr, "SGR")
)

kable(
  metricas_secundarias,
  digits = 4,
  caption = "Modelos de respuestas productivas complementarias"
)
Modelos de respuestas productivas complementarias
Respuesta R2 R2_ajustado F_global p_global
Supervivencia 0.5089 0.3919 4.3515 0.0071
SGR 0.9338 0.9180 59.2363 0.0000

Interpretación. La supervivencia presenta un ajuste más moderado que las variables de crecimiento y eficiencia alimenticia, lo que indica una mayor proporción de variabilidad experimental no explicada por proteína y lípidos. En contraste, el SGR presenta un ajuste elevado y una estructura similar a la ganancia de peso. En un ensayo real sería recomendable registrar variables ambientales —oxígeno disuelto, temperatura, salinidad, pH, amonio y otros parámetros de calidad de agua— para evaluar fuentes adicionales de variabilidad en supervivencia.

6. Prueba de falta de ajuste

La existencia de réplicas permite separar el error residual en error puro y falta de ajuste. Una falta de ajuste no significativa es deseable, porque indica que el modelo cuadrático representa adecuadamente la tendencia de los datos dentro de la región experimental.

prueba_lof <- function(modelo, datos, respuesta) {

  y <- datos[[respuesta]]
  grupo <- interaction(datos$x1, datos$x2, drop = TRUE)
  media_grupo <- ave(y, grupo, FUN = mean)

  ss_error_puro <- sum((y - media_grupo)^2)
  gl_error_puro <- length(y) - nlevels(grupo)

  ss_residual <- sum(residuals(modelo)^2)
  gl_residual <- df.residual(modelo)

  ss_lof <- ss_residual - ss_error_puro
  gl_lof <- gl_residual - gl_error_puro

  cm_error_puro <- ss_error_puro / gl_error_puro
  cm_lof <- ss_lof / gl_lof

  F_lof <- cm_lof / cm_error_puro
  p_lof <- pf(F_lof, gl_lof, gl_error_puro, lower.tail = FALSE)

  tibble(
    Fuente = c("Falta de ajuste", "Error puro"),
    SC = c(ss_lof, ss_error_puro),
    gl = c(gl_lof, gl_error_puro),
    CM = c(cm_lof, cm_error_puro),
    F = c(F_lof, NA),
    p = c(p_lof, NA)
  )
}

lof_peso <- prueba_lof(
  modelo_peso, datos, "ganancia_peso_g"
)

lof_fcr <- prueba_lof(
  modelo_fcr, datos, "FCR"
)

lof_surv <- prueba_lof(
  modelo_surv, datos, "supervivencia_pct"
)

lof_sgr <- prueba_lof(
  modelo_sgr, datos, "SGR_pct_dia"
)
tabla_lof <- tibble(
  Respuesta = c(
    "Ganancia de peso",
    "FCR",
    "Supervivencia",
    "SGR"
  ),
  p_falta_ajuste = c(
    lof_peso$p[1],
    lof_fcr$p[1],
    lof_surv$p[1],
    lof_sgr$p[1]
  )
)

kable(
  tabla_lof,
  digits = 4,
  caption = "Prueba de falta de ajuste de los modelos"
)
Prueba de falta de ajuste de los modelos
Respuesta p_falta_ajuste
Ganancia de peso 0.1969
FCR 0.4514
Supervivencia 0.6322
SGR 0.8023

Interpretación. Para la ganancia de peso, la prueba de falta de ajuste presenta \(p=\) 0.197; para FCR, \(p=\) 0.451. Al ser superiores a 0,05, no existe evidencia suficiente para afirmar que el modelo cuadrático presente una falta de ajuste sistemática en estas respuestas. Esto respalda el uso del modelo de segundo orden dentro de la región experimental evaluada.

7. Diagnóstico de residuos

7.1 Ganancia de peso

diag_peso <- tibble(
  ajustado = fitted(modelo_peso),
  residuo = residuals(modelo_peso),
  residuo_est = rstandard(modelo_peso)
)

ggplot(diag_peso, aes(ajustado, residuo_est)) +
  geom_point(size = 2.7) +
  geom_hline(yintercept = 0, linetype = 2) +
  labs(
    title = "Diagnóstico del modelo: residuos vs valores ajustados",
    x = "Valor ajustado de ganancia de peso",
    y = "Residuo estandarizado"
  ) +
  theme_minimal(base_size = 12)

ggplot(diag_peso, aes(sample = residuo_est)) +
  stat_qq(size = 2.5) +
  stat_qq_line() +
  labs(
    title = "Gráfico Q-Q de residuos - ganancia de peso"
  ) +
  theme_minimal(base_size = 12)

shapiro_peso <- shapiro.test(residuals(modelo_peso))

tibble(
  Prueba = "Shapiro-Wilk",
  W = unname(shapiro_peso$statistic),
  p = shapiro_peso$p.value
) %>%
  kable(
    digits = 4,
    caption = "Normalidad de los residuos del modelo de ganancia"
  )
Normalidad de los residuos del modelo de ganancia
Prueba W p
Shapiro-Wilk 0.9823 0.9107

Interpretación. El gráfico de residuos debe mostrar dispersión aproximadamente aleatoria alrededor de cero, sin una forma de embudo ni patrón curvilíneo evidente. El gráfico Q-Q permite evaluar normalidad aproximada. La prueba de Shapiro-Wilk presenta \(p=\) 0.911; un valor superior a 0,05 es compatible con el supuesto de normalidad de los residuos en este ejercicio.

7.2 FCR

diag_fcr <- tibble(
  ajustado = fitted(modelo_fcr),
  residuo_est = rstandard(modelo_fcr)
)

ggplot(diag_fcr, aes(ajustado, residuo_est)) +
  geom_point(size = 2.7) +
  geom_hline(yintercept = 0, linetype = 2) +
  labs(
    title = "Diagnóstico del modelo FCR",
    x = "FCR ajustado",
    y = "Residuo estandarizado"
  ) +
  theme_minimal(base_size = 12)

8. Superficies de respuesta

rejilla <- tidyr::expand_grid(
  x1 = seq(-1, 1, length.out = 151),
  x2 = seq(-1, 1, length.out = 151)
) %>%
  mutate(
    proteina_pct = 36 + 4*x1,
    lipidos_pct = 8 + 2*x2,
    pred_ganancia = predict(modelo_peso, newdata = .),
    pred_FCR = predict(modelo_fcr, newdata = .),
    pred_supervivencia = predict(modelo_surv, newdata = .),
    pred_SGR = predict(modelo_sgr, newdata = .)
  )

8.1 Mapa de contornos: ganancia de peso

ggplot(
  rejilla,
  aes(
    x = proteina_pct,
    y = lipidos_pct,
    z = pred_ganancia
  )
) +
  geom_contour_filled(bins = 14) +
  geom_point(
    data = datos,
    aes(proteina_pct, lipidos_pct),
    inherit.aes = FALSE,
    size = 1.6
  ) +
  labs(
    title = "Superficie de respuesta de la ganancia de peso",
    subtitle = "Mapa de contornos del modelo cuadrático",
    x = "Proteína cruda (%)",
    y = "Lípidos (%)",
    fill = "Ganancia\nde peso (g)"
  ) +
  theme_minimal(base_size = 12)

Interpretación. Las curvas de nivel representan combinaciones de proteína y lípidos con igual ganancia de peso predicha. La región de mayor respuesta se ubica dentro de la zona central-alta del espacio experimental, en lugar de encontrarse necesariamente en una esquina del diseño. Esto refleja la presencia de curvatura y constituye una de las principales ventajas de la RSM frente a seleccionar simplemente el tratamiento con la media observada más alta.

8.2 Mapa de contornos: FCR

ggplot(
  rejilla,
  aes(
    x = proteina_pct,
    y = lipidos_pct,
    z = pred_FCR
  )
) +
  geom_contour_filled(bins = 14) +
  geom_point(
    data = datos,
    aes(proteina_pct, lipidos_pct),
    inherit.aes = FALSE,
    size = 1.6
  ) +
  labs(
    title = "Superficie de respuesta del FCR",
    subtitle = "En esta respuesta se busca el mínimo",
    x = "Proteína cruda (%)",
    y = "Lípidos (%)",
    fill = "FCR"
  ) +
  theme_minimal(base_size = 12)

8.3 Superficie tridimensional

xseq <- seq(-1, 1, length.out = 45)
yseq <- seq(-1, 1, length.out = 45)

z <- outer(
  xseq,
  yseq,
  Vectorize(function(a, b) {
    predict(
      modelo_peso,
      newdata = data.frame(x1 = a, x2 = b)
    )
  })
)

persp(
  x = 36 + 4*xseq,
  y = 8 + 2*yseq,
  z = z,
  theta = 40,
  phi = 25,
  expand = 0.75,
  ticktype = "detailed",
  xlab = "Proteína (%)",
  ylab = "Lípidos (%)",
  zlab = "Ganancia (g)",
  main = "Superficie de respuesta 3D - ganancia de peso"
)

9. Punto estacionario y optimización individual

Se determina el punto donde la derivada de la superficie respecto de ambos factores es igual a cero. Según los signos de los autovalores de la matriz Hessiana, este punto puede clasificarse como máximo, mínimo o punto silla.

punto_estacionario <- function(modelo, objetivo = "max") {

  b <- coef(modelo)

  H <- matrix(
    c(
      2*b["I(x1^2)"], b["x1:x2"],
      b["x1:x2"], 2*b["I(x2^2)"]
    ),
    nrow = 2,
    byrow = TRUE
  )

  lineales <- c(b["x1"], b["x2"])
  x <- solve(H, -lineales)

  autovalores <- eigen(H, symmetric = TRUE)$values

  tipo <- if (all(autovalores < 0)) {
    "Máximo"
  } else if (all(autovalores > 0)) {
    "Mínimo"
  } else {
    "Punto silla"
  }

  pred <- predict(
    modelo,
    newdata = data.frame(
      x1 = x[1],
      x2 = x[2]
    )
  )

  tibble(
    x1 = x[1],
    x2 = x[2],
    proteina_pct = 36 + 4*x[1],
    lipidos_pct = 8 + 2*x[2],
    respuesta_predicha = unname(pred),
    clasificacion = tipo
  )
}

opt_peso <- punto_estacionario(modelo_peso)
opt_fcr <- punto_estacionario(modelo_fcr)
opt_surv <- punto_estacionario(modelo_surv)
opt_sgr <- punto_estacionario(modelo_sgr)
tabla_optimos <- bind_rows(
  opt_peso %>% mutate(Respuesta = "Ganancia de peso"),
  opt_fcr %>% mutate(Respuesta = "FCR"),
  opt_surv %>% mutate(Respuesta = "Supervivencia"),
  opt_sgr %>% mutate(Respuesta = "SGR")
) %>%
  select(
    Respuesta,
    proteina_pct,
    lipidos_pct,
    respuesta_predicha,
    clasificacion
  )

kable(
  tabla_optimos,
  digits = 3,
  caption = "Puntos estacionarios estimados para cada respuesta"
)
Puntos estacionarios estimados para cada respuesta
Respuesta proteina_pct lipidos_pct respuesta_predicha clasificacion
Ganancia de peso 38.251 9.438 19.005 Máximo
FCR 38.049 8.886 1.256 Mínimo
Supervivencia 37.128 8.914 95.229 Máximo
SGR 38.049 9.040 4.233 Máximo

Interpretación. Para maximizar la ganancia de peso, el modelo estima una combinación cercana a 38.25% de proteína y 9.44% de lípidos, con una ganancia predicha aproximada de 19.00 g. El punto se clasifica como máximo dentro de la superficie cuadrática.

Para FCR, el mínimo estimado se localiza alrededor de 38.05% de proteína y 8.89% de lípidos, con FCR predicho de 1.256.

La cercanía entre las regiones óptimas de crecimiento y eficiencia alimenticia sugiere que es posible buscar una solución de compromiso con buen desempeño simultáneo.

10. Análisis canónico con el paquete rsm

canonical(modelo_peso_rsm)
## $xs
##        x1        x2 
## 0.5627436 0.7189274 
## 
## $eigen
## eigen() decomposition
## $values
## [1] -0.8128729 -1.6486271
## 
## $vectors
##          [,1]       [,2]
## x1 -0.3262564 -0.9452813
## x2 -0.9452813  0.3262564

El análisis canónico es una forma formal de estudiar la geometría de la superficie de respuesta. Los signos de sus raíces o autovalores permiten identificar si la región estacionaria corresponde a un máximo, mínimo o punto silla.

11. Optimización multirrespuesta mediante deseabilidad

Una formulación técnicamente atractiva no debería evaluarse únicamente por el peso. Se plantea una optimización simultánea con tres objetivos:

  1. maximizar ganancia de peso;
  2. minimizar FCR;
  3. maximizar supervivencia.

Se asigna doble ponderación a ganancia de peso y FCR y una ponderación simple a supervivencia.

mayor_mejor <- function(y) {
  (y - min(y)) / (max(y) - min(y))
}

menor_mejor <- function(y) {
  (max(y) - y) / (max(y) - min(y))
}

rejilla_opt <- rejilla %>%
  mutate(
    d_peso = mayor_mejor(pred_ganancia),
    d_fcr = menor_mejor(pred_FCR),
    d_surv = mayor_mejor(pred_supervivencia),

    deseabilidad = (
      d_peso^2 *
      d_fcr^2 *
      d_surv
    )^(1/5)
  )

optimo_multi <- rejilla_opt %>%
  slice_max(
    deseabilidad,
    n = 1,
    with_ties = FALSE
  ) %>%
  select(
    proteina_pct,
    lipidos_pct,
    pred_ganancia,
    pred_FCR,
    pred_supervivencia,
    pred_SGR,
    deseabilidad
  )

kable(
  optimo_multi,
  digits = 3,
  caption = "Óptimo multirrespuesta estimado"
)
Óptimo multirrespuesta estimado
proteina_pct lipidos_pct pred_ganancia pred_FCR pred_supervivencia pred_SGR deseabilidad
37.867 9.067 18.968 1.256 95.193 4.233 0.993 
ggplot(
  rejilla_opt,
  aes(
    x = proteina_pct,
    y = lipidos_pct,
    z = deseabilidad
  )
) +
  geom_contour_filled(bins = 14) +
  geom_point(
    data = optimo_multi,
    aes(proteina_pct, lipidos_pct),
    inherit.aes = FALSE,
    size = 4,
    shape = 8
  ) +
  labs(
    title = "Optimización multirrespuesta por deseabilidad",
    subtitle = "Maximizar crecimiento y supervivencia; minimizar FCR",
    x = "Proteína cruda (%)",
    y = "Lípidos (%)",
    fill = "Deseabilidad"
  ) +
  theme_minimal(base_size = 12)

Interpretación. La optimización simultánea localiza una solución de compromiso alrededor de 37.87% de proteína y 9.07% de lípidos. En ese punto el modelo predice aproximadamente 18.97 g de ganancia de peso, FCR = 1.256, supervivencia de 95.19% y una deseabilidad global de 0.993.

Esta solución es más útil para toma de decisiones que optimizar una sola variable, porque integra simultáneamente crecimiento, eficiencia alimenticia y supervivencia.

12. Comparación entre valor observado y óptimo modelado

comparacion <- tibble(
  Criterio = c(
    "Mayor ganancia observada",
    "Máximo de ganancia según RSM",
    "Mínimo FCR según RSM",
    "Óptimo multirrespuesta"
  ),
  Proteina_pct = c(
    mejor_peso_obs$proteina_pct,
    opt_peso$proteina_pct,
    opt_fcr$proteina_pct,
    optimo_multi$proteina_pct
  ),
  Lipidos_pct = c(
    mejor_peso_obs$lipidos_pct,
    opt_peso$lipidos_pct,
    opt_fcr$lipidos_pct,
    optimo_multi$lipidos_pct
  )
)

kable(
  comparacion,
  digits = 2,
  caption = "Comparación de regiones de interés"
)
Comparación de regiones de interés
Criterio Proteina_pct Lipidos_pct
Mayor ganancia observada 40.00 10.00
Máximo de ganancia según RSM 38.25 9.44
Mínimo FCR según RSM 38.05 8.89
Óptimo multirrespuesta 37.87 9.07

La comparación permite distinguir entre selección del mejor tratamiento ensayado y optimización mediante una superficie continua. Esa diferencia es central en RSM: el mejor tratamiento observado se limita a las combinaciones ejecutadas, mientras que el modelo permite estimar una combinación intermedia potencialmente superior.

13. Implicaciones técnicas para un ensayo real

Los resultados simulados sugieren una región central-alta de proteína y lípidos como zona de interés. Sin embargo, una decisión real de formulación debería incorporar variables adicionales, entre ellas:

Además, si los tanques estuvieran distribuidos espacialmente en diferentes áreas del centro experimental, podría ser conveniente aplicar bloqueo para controlar gradientes ambientales.

14. Propuesta de ensayo confirmatorio

La optimización estadística no debería considerarse el punto final del proceso. La formulación óptima predicha debe validarse experimentalmente.

Un ensayo confirmatorio podría incluir:

  1. Dieta control comercial.
  2. Dieta base o formulación estándar.
  3. Dieta óptima propuesta por la RSM.

La hipótesis de validación sería que la respuesta observada con la dieta óptima se encuentre dentro del intervalo de predicción generado por el modelo y mantenga ventajas productivas frente al control.

En un proyecto de I+D, el flujo completo sería:

\[ \text{Diseño} \rightarrow \text{Ejecución} \rightarrow \text{Captura de datos} \rightarrow \text{Modelado} \rightarrow \text{Optimización} \rightarrow \text{Validación} \rightarrow \text{Escalamiento} \]

15. Conclusiones

  1. El diseño permite evaluar simultáneamente los efectos lineales, la interacción proteína × lípidos y la curvatura de la respuesta.

  2. La ganancia de peso y el FCR presentan una estructura cuadrática marcada en los datos simulados, lo que justifica el uso de Metodología de Superficie de Respuesta.

  3. El modelo de ganancia de peso explica aproximadamente 96.1% de la variabilidad, mientras que el modelo de FCR explica aproximadamente 98.3%.

  4. La falta de ajuste no significativa en las respuestas principales respalda el uso de la forma cuadrática dentro de la región experimental.

  5. El óptimo individual de crecimiento se estima alrededor de 38.25% de proteína y 9.44% de lípidos.

  6. La optimización multirrespuesta propone aproximadamente 37.87% de proteína y 9.07% de lípidos, buscando simultáneamente mayor crecimiento, menor FCR y alta supervivencia.

  7. En una aplicación comercial, la recomendación final debe complementarse con un ensayo confirmatorio y un análisis económico de la formulación.

16. Nota metodológica

Este documento constituye un ejercicio técnico demostrativo basado en datos simulados. Su objetivo es mostrar competencias en diseño experimental, estadística aplicada, R, R Markdown, modelado de superficies de respuesta y comunicación de resultados técnicos en un contexto de investigación nutricional para acuicultura.