1. Introducción y Planteamiento del Problema

En los sistemas de producción bovina de doble propósito en Panamá, particularmente en las lecherías Grado C de los distritos de Alanje, Bugaba y Boquerón en la provincia de Chiriquí, la alimentación se fundamenta en el pastoreo extensivo. Estos sistemas presentan grandes retos productivos y bioeconómicos, ya que dependen de pastos de baja calidad o que fluctúan severamente en su biomasa y nutrientes a lo largo del año (Camarena, 2011; González, 2015). La literatura señala que la calidad bromatológica de pasturas naturalizadas disminuye drásticamente en épocas críticas, mientras que pastos introducidos y mejorados como Brachiaria decumbens, Pasto Estrella (Cynodon nlemfuensis) y Mombasa (Panicum maximum) ofrecen distintos y superiores perfiles nutricionales (Guillén Moreno, 2005).

Para mitigar esta caída productiva, Aguilar-Pérez et al. (2019) y Bottini-Luzardo et al. (2016) sugieren la implementación de sistemas forrajeros de alta calidad y la suplementación estratégico-energética, demostrando que mejorar la base forrajera eleva el rendimiento diario de leche. Sin embargo, surge la interrogante de si el efecto de suplementar es igual de eficiente en pastos de baja calidad que en pastos mejorados, como lo evalúan Peniche-González et al. (2014) al indicar que las asociaciones forrajeras optimizan el uso de los suplementos. Además, la respuesta biológica en las lecherías panameñas está fuertemente influenciada por la etapa productiva de las vacas o el número de partos (Camarena, 2011).

Para analizar esto de forma rigurosa, se plantea un Diseño de Bloques Completos al Azar (DBCA) con Arreglo Factorial 4x2. Los bloques aislarán la heterogeneidad de los partos, mientras que el arreglo factorial evaluará el Factor A (Tipo de Forraje: 4 niveles) y el Factor B (Suplementación: 2 niveles).

2. Pregunta de Investigación y Objetivos

Pregunta de Investigación: ¿Cuál es el efecto principal y el efecto de interacción del Tipo de Forraje (Faragua, B. decumbens, Estrella, Mombasa) y el Nivel de Suplementación (Sin Suplemento vs. Con Suplemento) sobre el rendimiento de leche (L/vaca/día) en sistemas de doble propósito Grado C, aislando y controlando la variabilidad biológica del número de partos de las vacas?

Objetivos:

  1. Evaluar el efecto principal de cuatro tipos de forraje en la producción láctea.
  2. Evaluar el efecto principal del nivel de suplementación en la producción láctea.
  3. Determinar si existe un efecto de interacción (sinergia) entre la base forrajera y la suplementación.
  4. Aislar y cuantificar la variabilidad atribuible al número de partos (Bloques).

3. Modelo Estadístico (Arreglo Factorial en Bloques)

El modelo matemático lineal para un diseño factorial \(4 \times 2\) en bloques está dado por:

\[ Y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + \gamma_k + \varepsilon_{ijk} \]

Donde:

  • \(Y_{ijk}\): Producción de leche (L/día) bajo el nivel \(i\) de Forraje, nivel \(j\) de Suplementación en el bloque \(k\).
  • \(\mu\): Media global del experimento.
  • \(\alpha_i\): Efecto fijo del \(i\)-ésimo nivel del Factor A (Forraje: Faragua, B. decumbens, Estrella, Mombasa).
  • \(\beta_j\): Efecto fijo del \(j\)-ésimo nivel del Factor B (Suplementación: Sin Supl, Con Supl).
  • \((\alpha\beta)_{ij}\): Efecto de la interacción Forraje-Suplementación.
  • \(\gamma_k\): Efecto del \(k\)-ésimo Bloque (Partos: 1 a 5).
  • \(\varepsilon_{ijk}\): Error experimental aleatorio, \(\varepsilon_{ijk} \sim N(0, \sigma^2)\).

4. Simulación y Diseño de Campo (Croquis)

Utilizamos 40 vacas alojadas en parcelas individuales de pastoreo, divididas en 5 bloques biológicos según sus partos (8 vacas por bloque). Las combinaciones de los factores A (4) y B (2) conforman los 8 tratamientos experimentales.

4.1. Generación del Croquis Factorial Aleatorizado

set.seed(123) # Semilla para reproducibilidad
niveles_forraje <- c("Faragua", "B.decumbens", "Estrella", "Mombasa")
niveles_supl <- c("Sin_Supl", "Con_Supl")
niveles_bloque <- c("1_Parto", "2_Partos", "3_Partos", "4_Partos", "5_Partos")

# Creamos las 8 combinaciones factoriales
trt_comb <- paste(rep(niveles_forraje, each=2), rep(niveles_supl, times=4), sep="-")

# Generamos el diseño formal en bloques con agricolae
diseno_dbca <- design.rcbd(trt = trt_comb, r = 5, seed = 123)

# Estructuramos la base de datos dividiendo los factores para el análisis
datos_campo <- diseno_dbca$book %>%
  setNames(c("Parcela", "Bloque_Num", "Tratamiento")) %>%
  separate(Tratamiento, into = c("Forraje", "Suplemento"), sep = "-", remove = FALSE) %>%
  mutate(
    Bloque = factor(niveles_bloque[as.numeric(Bloque_Num)], levels = niveles_bloque),
    Forraje = factor(Forraje, levels = niveles_forraje),
    Suplemento = factor(Suplemento, levels = niveles_supl)
  ) %>%
  group_by(Bloque_Num) %>%
  mutate(Posicion = row_number()) %>%
  ungroup()

# Visualización del Croquis
ggplot(datos_campo, aes(x = factor(Posicion), y = forcats::fct_rev(Bloque), fill = Tratamiento)) +
  geom_tile(color = "black", width = 0.95, height = 0.9, size = 0.5) +
  # Etiquetas de texto con tamaño reducido y salto de línea para mejorar estética
  geom_text(aes(label = paste0(Forraje, "\n", Suplemento)), 
            color = "black", fontface = "bold", size = 2.8, lineheight = 0.9) +
  scale_fill_brewer(palette = "Set3") +
  theme_minimal() +
  labs(
    title = "Croquis DBCA Factorial 4x2 en Campo",
    subtitle = "Aleatorización de 8 combinaciones en 5 Bloques (40 Unidades Experimentales)",
    x = "Posición Aleatoria de la Parcela",
    y = "Bloque Biológico (Número de Partos)"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14, family = "sans"),
    plot.subtitle = element_text(hjust = 0.5, size = 12, margin=margin(b=10)),
    axis.title = element_text(face = "bold"),
    axis.text.y = element_text(face = "bold", size = 10),
    axis.text.x = element_blank(),
    panel.grid = element_blank(),
    legend.position = "none"
  )
Figura 1. Croquis de Aleatorización en Campo (DBCA Factorial 4x2)

Figura 1. Croquis de Aleatorización en Campo (DBCA Factorial 4x2)

4.2. Simulación de Rendimiento Lechero Posterior a 60 Días

Se simula la recolección de datos tras 60 días de evaluación, integrando los efectos del forraje y suplemento basados en las estimaciones de Aguilar-Pérez et al. (2019), incluyendo el efecto del bloque (partos).

# Efectos biológicos teóricos (L/día)
efecto_forraje <- c("Faragua" = 4.5, "B.decumbens" = 5.2, "Estrella" = 5.5, "Mombasa" = 6.0)
efecto_supl <- c("Sin_Supl" = 0, "Con_Supl" = 1.2)
efecto_bloque <- c("1_Parto" = -2.0, "2_Partos" = -0.5, "3_Partos" = 1.5, "4_Partos" = 2.0, "5_Partos" = -1.0)

# Aplicamos los efectos al marco de datos con interacción sinérgica
datos_exp <- datos_campo %>%
  mutate(
    Efecto_Inter = case_when(
      Suplemento == "Con_Supl" & Forraje == "Mombasa" ~ 1.5,
      Suplemento == "Con_Supl" & Forraje == "Estrella" ~ 1.0,
      Suplemento == "Con_Supl" & Forraje == "B.decumbens" ~ 0.6,
      Suplemento == "Con_Supl" & Forraje == "Faragua" ~ 0.1,
      TRUE ~ 0
    ),
    Produccion_L_dia = round(
      efecto_forraje[as.character(Forraje)] + 
      efecto_supl[as.character(Suplemento)] + 
      Efecto_Inter + 
      efecto_bloque[as.character(Bloque)] + 
      rnorm(n(), mean = 0, sd = 0.3), 2)
  ) %>%
  select(Parcela, Bloque, Forraje, Suplemento, Produccion_L_dia, Tratamiento)

# Presentación de la tabla de datos
datos_exp %>%
  head(10) %>%
  kbl(caption = "Tabla 1. Muestra de Datos Simulados Recolectados (L/vaca/día)", align = "c") %>%
  kable_classic(full_width = FALSE, html_font = "Arial") %>%
  row_spec(0, bold = TRUE, background = "#ecf0f1")
Tabla 1. Muestra de Datos Simulados Recolectados (L/vaca/día)
Parcela Bloque Forraje Suplemento Produccion_L_dia Tratamiento
101 1_Parto Mombasa Sin_Supl 3.82 Mombasa-Sin_Supl
102 1_Parto Faragua Sin_Supl 2.34 Faragua-Sin_Supl
103 1_Parto Mombasa Con_Supl 6.82 Mombasa-Con_Supl
104 1_Parto B.decumbens Con_Supl 4.87 B.decumbens-Con_Supl
105 1_Parto Estrella Sin_Supl 3.57 Estrella-Sin_Supl
106 1_Parto Faragua Con_Supl 4.13 Faragua-Con_Supl
107 1_Parto Estrella Con_Supl 5.46 Estrella-Con_Supl
108 1_Parto B.decumbens Sin_Supl 2.67 B.decumbens-Sin_Supl
201 2_Partos B.decumbens Con_Supl 6.99 B.decumbens-Con_Supl
202 2_Partos Faragua Con_Supl 5.56 Faragua-Con_Supl

5. Análisis de Varianza (ANOVA Factorial en Bloques)

Se examina el efecto de los factores y su interacción.

# Modelo con interacción (Forraje * Suplemento) ajustando por Bloque
modelo_fact <- aov(Produccion_L_dia ~ Forraje * Suplemento + Bloque, data = datos_exp)

tidy(modelo_fact) %>%
  mutate(
    p.value = ifelse(p.value < 0.001, "< 0.001", as.character(round(p.value, 4))),
    statistic = round(statistic, 2),
    sumsq = round(sumsq, 2),
    meansq = round(meansq, 2)
  ) %>%
  kbl(caption = "Tabla 2. Tabla ANOVA - Diseño Factorial 4x2 en Bloques", align = "c") %>%
  kable_classic(full_width = FALSE, html_font = "Arial") %>%
  row_spec(0, bold = TRUE, background = "#ecf0f1") %>%
  row_spec(1:3, bold = TRUE, background = "#d9edf7") # Resaltar factores e interacción
Tabla 2. Tabla ANOVA - Diseño Factorial 4x2 en Bloques
term df sumsq meansq statistic p.value
Forraje 3 24.64 8.21 114.61 < 0.001
Suplemento 1 44.42 44.42 619.70 < 0.001
Bloque 4 89.53 22.38 312.28 < 0.001
Forraje:Suplemento 3 2.28 0.76 10.59 < 0.001
Residuals 28 2.01 0.07 NA NA

6. Validación de Supuestos

Para confirmar la validez del ANOVA, evaluamos gráficamente y formalmente los supuestos de los residuos.

par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
plot(fitted(modelo_fact), residuals(modelo_fact), pch = 19, col = "#2980b9",
     main = "Residuos vs Ajustados", xlab = "Valores Ajustados", ylab = "Residuos")
abline(h = 0, col = "#c0392b", lty = 2, lwd=2)

qqnorm(residuals(modelo_fact), pch = 19, col = "#2980b9", main = "Q-Q Plot Normal")
qqline(residuals(modelo_fact), col = "#c0392b", lwd = 2)
Figura 2. Diagnóstico Gráfico de los Residuos

Figura 2. Diagnóstico Gráfico de los Residuos 

par(mfrow = c(1, 1))

# Pruebas formales
shapiro_test <- shapiro.test(residuals(modelo_fact))
levene_test <- leveneTest(Produccion_L_dia ~ Tratamiento, data = datos_exp)
dw_test <- dwtest(modelo_fact)

data.frame(
  Supuesto = c("Normalidad (Shapiro-Wilk)", "Homocedasticidad (Levene)", "Independencia (Durbin-Watson)"),
  Estadistico = c(round(shapiro_test$statistic, 4), round(levene_test$`F value`[1], 4), round(dw_test$statistic, 4)),
  Valor_P = c(round(shapiro_test$p.value, 4), round(levene_test$`Pr(>F)`[1], 4), round(dw_test$p.value, 4))
) %>%
  kbl(caption = "Tabla 3. Validación Formal de Supuestos del Modelo", align = "c") %>%
  kable_classic(full_width = FALSE, html_font = "Arial") %>%
  row_spec(0, bold = TRUE, background = "#ecf0f1")
Tabla 3. Validación Formal de Supuestos del Modelo
Supuesto Estadistico Valor_P
W Normalidad (Shapiro-Wilk) 0.9492 0.0710
Homocedasticidad (Levene) 0.1233 0.9961
DW Independencia (Durbin-Watson) 2.4069 0.7561

Se confirma el estricto cumplimiento de los supuestos (\(p > 0.05\)).

7. Pruebas No Paramétricas de Respaldo

A modo de verificación cruzada (útil en muestras pequeñas de campo agronómico), empleamos la prueba de Friedman para bloques.

friedman_res <- friedman.test(Produccion_L_dia ~ Tratamiento | Bloque, data = datos_exp)

data.frame(
  Prueba = "Friedman (DBCA Factorial)",
  Chi_Cuadrado = round(friedman_res$statistic, 4),
  Valor_p = ifelse(friedman_res$p.value < 0.001, "< 0.001", round(friedman_res$p.value, 4))
) %>%
  kbl(caption = "Tabla 4. Prueba No Paramétrica de Confirmación (Friedman)", align = "c") %>%
  kable_classic(full_width = FALSE, html_font = "Arial") %>%
  row_spec(0, bold = TRUE, background = "#ecf0f1")
Tabla 4. Prueba No Paramétrica de Confirmación (Friedman)
Prueba Chi_Cuadrado Valor_p
Friedman chi-squared Friedman (DBCA Factorial) 34.5322 < 0.001

8. Prueba Post-Hoc: Comparación de Tukey e Interacción

Habiéndose detectado una interacción altamente significativa, evaluamos las medias de los 8 tratamientos combinados mediante Tukey HSD y las graficamos.

# Modelo lineal sobre la variable combinada 'Tratamiento' para la post-hoc
modelo_tukey <- aov(Produccion_L_dia ~ Tratamiento + Bloque, data = datos_exp)
tukey_out <- HSD.test(modelo_tukey, "Tratamiento", group = TRUE)

data.frame(
  Tratamiento = rownames(tukey_out$groups),
  Media_L_dia = round(tukey_out$groups$Produccion_L_dia, 2),
  Grupo = tukey_out$groups$groups
) %>%
  kbl(caption = "Tabla 5. Prueba de Rango Múltiple de Tukey (HSD)", align = "c") %>%
  kable_classic(full_width = FALSE, html_font = "Arial") %>%
  row_spec(0, bold = TRUE, background = "#ecf0f1") %>%
  row_spec(which(tukey_out$groups$groups == "a"), bold = TRUE, color = "white", background = "#27ae60")
Tabla 5. Prueba de Rango Múltiple de Tukey (HSD)
Tratamiento Media_L_dia Grupo
Mombasa-Con_Supl 8.69 a
Estrella-Con_Supl 7.67 b
B.decumbens-Con_Supl 7.06 c
Mombasa-Sin_Supl 6.04 d
Faragua-Con_Supl 5.89 de
Estrella-Sin_Supl 5.41 ef
B.decumbens-Sin_Supl 4.90 fg
Faragua-Sin_Supl 4.54 g 
# Gráfico de Interacción Formal (APA)
resumen_inter <- datos_exp %>%
  group_by(Forraje, Suplemento) %>%
  summarise(Media = mean(Produccion_L_dia), .groups = 'drop')

ggplot(resumen_inter, aes(x = Forraje, y = Media, color = Suplemento, group = Suplemento)) +
  geom_line(size = 1.2) +
  geom_point(size = 4, shape = 19) +
  scale_color_manual(values = c("Sin_Supl" = "#c0392b", "Con_Supl" = "#2980b9")) +
  theme_classic() +
  labs(
    title = "Efecto de Interacción: Base Forrajera vs Nivel de Suplementación",
    y = "Producción Promedio Láctea (L/vaca/día)",
    x = "Tipo de Forraje"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 13),
    axis.title = element_text(face = "bold"),
    axis.text = element_text(size = 11, color = "black"),
    legend.title = element_text(face="bold")
  )
Figura 3. Efecto de Interacción (Forraje x Suplemento)

Figura 3. Efecto de Interacción (Forraje x Suplemento)

9. Conclusiones Agropecuarias y Metodológicas

A partir de los resultados inferenciales del modelo factorial simulado, se establecen las siguientes conclusiones, validadas frente a la literatura regional:

  1. La Importancia del Bloqueo (Efecto del Parto): El análisis de varianza (Tabla 2) confirmó que el factor biológico (número de partos) introdujo una varianza altamente significativa. El empleo de un diseño de bloques eliminó exitosamente este ruido experimental, permitiendo observar los efectos puros de la dieta, lo que responde a las advertencias metodológicas sobre la alta heterogeneidad en explotaciones doble propósito en la provincia de Chiriquí (Camarena, 2011; González, 2015).
  2. Sinergia Nutricional (Interacción Forraje x Suplemento): La Figura 3 y el ANOVA evidencian una interacción multiplicativa (sinérgica) significativa. Suplementar vacas sobre forrajes de baja calidad o altamente degradados (como el Faragua) generó un retorno productivo marginal. Sin embargo, administrar suplemento a animales que pastoreaban forrajes mejorados (Estrella y Mombasa) disparó el rendimiento. Esto es congruente con lo reportado por Peniche-González et al. (2014) y Aguilar-Pérez et al. (2019), quienes señalan que una base forrajera rica en digestibilidad permite a los rumiantes un aprovechamiento eficiente de los suplementos energéticos y proteicos.
  3. Recomendación Práctica (Tukey): Según la prueba de Tukey HSD (Tabla 5), la combinación Mombasa-Con_Supl conforma un clúster estadístico superior aislado (Grupo “a”). Los productores panameños de leche Grado C deben concentrar sus esfuerzos de inversión, en primera instancia, en la reconversión de sus praderas y el establecimiento de pastos mejorados como un paso previo, o simultáneo, a la integración de altos niveles de concentrados, ya que esto último incrementa drásticamente los costos de producción sin garantizar el retorno productivo adecuado (Guillén Moreno, 2005).

Referencias

Aguilar-Pérez, C., Solorio-Sánchez, F., Ku-Vera, J., Magaña-Monforte, J., & Santos-Flores, J. (2019). Producción de leche y carne en sistemas silvopastoriles. Bioagrociencias, 12(1), 1-8.

Bottini-Luzardo, M. B., Aguilar-Pérez, C. F., Centurión-Castro, F. G., Solorio-Sánchez, F. J., & Ku-Vera, J. C. (2016). Milk yield and blood urea nitrogen in crossbred cows grazing L. leucocephala in a silvopastoral system in the Mexican tropics. Tropical Grasslands-Forrajes Tropicales, 4(3), 159-167.

Camarena, M. I. (2011). Situación reproductiva y lactacional en fincas lecheras grado C en los distritos de Alanje, Bugaba y Boquerón y su relación con la producción de leche y la condición bioeconómica (Tesis de licenciatura). Universidad de Panamá, David, Chiriquí.

González, M. M. (2015). Evaluación bioeconómica mediante un modelo de simulación (LIFE-SIM) para la producción bovina de leche en la zona baja de la provincia de Chiriquí (Tesis de licenciatura). Universidad de Panamá, David, Chiriquí.

Guillén Moreno, A. J. (2005). Caracterización de los sistemas de alimentación y su influencia sobre la producción y composición láctea en fincas productoras de leche con condiciones tecnológicas grado A en la provincia de Chiriquí (Tesis de maestría). Universidad de Panamá, David, Chiriquí.

Peniche-González, I. N., González-López, Z. U., Aguilar-Pérez, C. F., Ku-Vera, J. C., Ayala-Burgos, A. J., & Solorio-Sánchez, F. J. (2014). Milk production and reproduction of dual-purpose cows with a restricted concentrate allowance and access to an association of Leucaena leucocephala and Cynodon nlemfuensis. Journal of Applied Animal Research, 42(3), 345-351.