1 Resumen

Este informe compara la producción registrada de maíz y arroz en las provincias de Herrera y Los Santos durante tres años agrícolas. Se presentan dos estructuras de análisis: un diseño completamente aleatorizado (DCA), considerando las cuatro combinaciones provincia–cultivo como tratamientos, y un diseño de bloques completos al azar (DBCA), utilizando el año agrícola como bloque.

También se describe la precipitación anual y se explora su asociación con la producción mediante correlaciones y regresión lineal simple. Debido a que los datos son registros históricos agregados y no provienen de un experimento con asignación aleatoria de tratamientos, los diseños se emplean como estructuras estadísticas de comparación; sus resultados no deben interpretarse como evidencia causal de un experimento aleatorizado.

2 Objetivos

  • Describir la producción agrícola y la precipitación anual por provincia y año.
  • Comparar las cuatro combinaciones de provincia y cultivo mediante un ANOVA bajo la estructura DCA.
  • Incorporar el año agrícola como bloque en el modelo DBCA.
  • Revisar los supuestos de normalidad de residuos y homogeneidad de varianzas.
  • Explorar la relación entre precipitación y producción para cada cultivo.

3 Datos y metodología

3.1 Base de datos

Los valores de producción se expresan en quintales y la precipitación en milímetros. En la tabla original, el tercer año aparece escrito como 2016/2015; para mantener la secuencia cronológica se registra aquí como 2016-2017. Esta corrección de etiqueta debe confirmarse con la fuente original.

datos <- data.frame(
  Anio = c("2014-2015", "2015-2016", "2016-2017",
           "2014-2015", "2015-2016", "2016-2017"),
  Provincia = c(rep("Herrera", 3), rep("Los Santos", 3)),
  Precipitacion = c(940.5, 860.6, 1371.3, 704.0, 644.3, 1232.1),
  Maiz = c(182924, 440500, 470000, 1378443, 1302800, 1701100),
  Arroz = c(152922, 66412, 116953, 702391, 648100, 936960)
)

datos$Anio <- factor(datos$Anio,
                     levels = c("2014-2015", "2015-2016", "2016-2017"))
datos$Provincia <- factor(datos$Provincia,
                          levels = c("Herrera", "Los Santos"))

datos_largos <- datos |>
  pivot_longer(
    cols = c(Maiz, Arroz),
    names_to = "Cultivo",
    values_to = "Produccion"
  ) |>
  mutate(
    Cultivo = factor(Cultivo, levels = c("Maiz", "Arroz")),
    Tratamiento = interaction(Provincia, Cultivo, sep = " – ", lex.order = TRUE)
  )

formato_tabla(datos, 2)
Anio Provincia Precipitacion Maiz Arroz
2014-2015 Herrera 940.5 182924 152922
2015-2016 Herrera 860.6 440500 66412
2016-2017 Herrera 1371.3 470000 116953
2014-2015 Los Santos 704.0 1378443 702391
2015-2016 Los Santos 644.3 1302800 648100
2016-2017 Los Santos 1232.1 1701100 936960

3.2 Tratamientos y bloques

Las cuatro combinaciones utilizadas como tratamientos son:

  1. Herrera – Maíz.
  2. Herrera – Arroz.
  3. Los Santos – Maíz.
  4. Los Santos – Arroz.

En el DBCA, cada año agrícola constituye un bloque y contiene una observación para cada una de las cuatro combinaciones. Por tanto, hay cuatro tratamientos, tres bloques y doce observaciones en el conjunto largo.

estructura <- datos_largos |>
  count(Anio, Tratamiento, name = "Observaciones")

formato_tabla(estructura)
Anio Tratamiento Observaciones
2014-2015 Herrera – Maiz 1
2014-2015 Herrera – Arroz 1
2014-2015 Los Santos – Maiz 1
2014-2015 Los Santos – Arroz 1
2015-2016 Herrera – Maiz 1
2015-2016 Herrera – Arroz 1
2015-2016 Los Santos – Maiz 1
2015-2016 Los Santos – Arroz 1
2016-2017 Herrera – Maiz 1
2016-2017 Herrera – Arroz 1
2016-2017 Los Santos – Maiz 1
2016-2017 Los Santos – Arroz 1

3.3 Estadística descriptiva

Se calculan el número de observaciones, media, mediana, desviación estándar, mínimo y máximo para cada combinación de provincia y cultivo.

descriptivos <- datos_largos |>
  group_by(Provincia, Cultivo) |>
  summarise(
    n = n(),
    Media = mean(Produccion),
    Mediana = median(Produccion),
    Desviacion = sd(Produccion),
    Minimo = min(Produccion),
    Maximo = max(Produccion),
    .groups = "drop"
  )

formato_tabla(descriptivos, 2)
Provincia Cultivo n Media Mediana Desviacion Minimo Maximo
Herrera Maiz 3 364474.7 440500 157917.85 182924 470000
Herrera Arroz 3 112095.7 116953 43459.06 66412 152922
Los Santos Maiz 3 1460781.0 1378443 211531.04 1302800 1701100
Los Santos Arroz 3 762483.7 702391 153519.94 648100 936960

3.4 Visualización de la producción

El gráfico compara la producción por año, cultivo y provincia. Se utilizan escalas verticales independientes por provincia para facilitar la lectura, ya que los volúmenes registrados son muy diferentes.

ggplot(datos_largos,
       aes(x = Anio, y = Produccion, fill = Cultivo)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.7) +
  facet_wrap(~ Provincia, scales = "free_y") +
  scale_fill_manual(values = c("Maiz" = "#e0a52b", "Arroz" = "#3c8d70"),
                    labels = c("Maíz", "Arroz")) +
  scale_y_continuous(labels = scales::label_comma()) +
  labs(
    title = "Producción agrícola por provincia y cultivo",
    subtitle = "Años agrícolas 2014–2015 a 2016–2017",
    x = "Año agrícola", y = "Producción (quintales)", fill = "Cultivo"
  ) +
  tema_informe +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))
Producción registrada de maíz y arroz por año agrícola y provincia.

Producción registrada de maíz y arroz por año agrícola y provincia.

4 Diseño completamente aleatorizado (DCA)

4.1 Planteamiento del modelo

Se ajusta un ANOVA de un factor con cuatro niveles de tratamiento. El modelo es:

\[ Y_{ij} = \mu + \tau_i + \varepsilon_{ij} \]

donde \(Y_{ij}\) representa la producción observada, \(\mu\) es la media general, \(\tau_i\) es el efecto de la combinación provincia–cultivo y \(\varepsilon_{ij}\) es el error aleatorio.

Las hipótesis son:

  • \(H_0\): las medias de producción de los cuatro tratamientos son iguales.
  • \(H_1\): al menos una media de producción difiere.
modelo_DCA <- aov(Produccion ~ Tratamiento, data = datos_largos)
anova_DCA <- anova(modelo_DCA)
formato_tabla(data.frame(Fuente = rownames(anova_DCA), anova_DCA,
                         check.names = FALSE), 4)
Fuente Df Sum Sq Mean Sq F value Pr(>F)
Tratamiento 3 3.115177e+12 1.038392e+12 43.6572 0
Residuals 8 1.902810e+11 2.378512e+10 NA NA

4.2 Modelo factorial complementario

Como análisis complementario, se separan los efectos de provincia, cultivo y su interacción. Este modelo permite examinar si la diferencia entre cultivos depende de la provincia.

modelo_DCA_factorial <- aov(Produccion ~ Provincia * Cultivo,
                            data = datos_largos)
anova_factorial_DCA <- anova(modelo_DCA_factorial)
formato_tabla(data.frame(Fuente = rownames(anova_factorial_DCA),
                         anova_factorial_DCA,
                         check.names = FALSE), 4)
Fuente Df Sum Sq Mean Sq F value Pr(>F)
Provincia 1 2.288206e+12 2.288206e+12 96.2032 0.0000
Cultivo 1 6.778391e+11 6.778391e+11 28.4984 0.0007
Provincia:Cultivo 1 1.491324e+11 1.491324e+11 6.2700 0.0367
Residuals 8 1.902810e+11 2.378512e+10 NA NA

4.3 Verificación de supuestos del DCA

La normalidad se examina sobre los residuos mediante Shapiro–Wilk y la homogeneidad de varianzas mediante Levene. Con solo tres observaciones por tratamiento, estas pruebas tienen poca potencia; por ello, deben interpretarse junto con los gráficos diagnósticos y no como una garantía de cumplimiento.

residuos_DCA <- residuals(modelo_DCA)

shapiro_DCA <- shapiro.test(residuos_DCA)
levene_DCA <- car::leveneTest(Produccion ~ Tratamiento,
                              data = datos_largos)

shapiro_DCA

    Shapiro-Wilk normality test

data:  residuos_DCA
W = 0.9677, p-value = 0.8853
levene_DCA
par(mfrow = c(2, 2))
plot(modelo_DCA)
Diagnóstico gráfico de los residuos del modelo DCA.

Diagnóstico gráfico de los residuos del modelo DCA. 

par(mfrow = c(1, 1))

4.4 Comparaciones múltiples y alternativa no paramétrica

Tukey permite comparar pares de medias bajo el modelo ANOVA. Kruskal–Wallis se presenta como alternativa no paramétrica para comparar los grupos, aunque tampoco incorpora bloques.

TukeyHSD(modelo_DCA)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Produccion ~ Tratamiento, data = datos_largos)

$Tratamiento
                                          diff          lwr       upr     p adj
Herrera – Arroz-Herrera – Maiz       -252379.0  -655630.633  150872.6 0.2626137
Los Santos – Maiz-Herrera – Maiz     1096306.3   693054.701 1499558.0 0.0001095
Los Santos – Arroz-Herrera – Maiz     398009.0    -5242.633  801260.6 0.0530183
Los Santos – Maiz-Herrera – Arroz    1348685.3   945433.701 1751937.0 0.0000237
Los Santos – Arroz-Herrera – Arroz    650388.0   247136.367 1053639.6 0.0037887
Los Santos – Arroz-Los Santos – Maiz -698297.3 -1101548.966 -295045.7 0.0024234
kruskal.test(Produccion ~ Tratamiento, data = datos_largos)

    Kruskal-Wallis rank sum test

data:  Produccion by Tratamiento
Kruskal-Wallis chi-squared = 10.385, df = 3, p-value = 0.01556

5 Diseño de bloques completos al azar (DBCA)

5.1 Planteamiento del modelo

Se incorpora el año agrícola como bloque para controlar diferencias entre años. El modelo aditivo es:

\[ Y_{ij} = \mu + \tau_i + \beta_j + \varepsilon_{ij} \]

donde \(\tau_i\) representa el efecto del tratamiento y \(\beta_j\) el efecto del bloque (año agrícola).

  • \(H_0\) para tratamientos: las cuatro medias de producción son iguales, una vez considerado el bloque.
  • \(H_1\): al menos una media de tratamiento difiere.
  • Para bloques, se evalúa si existen diferencias entre años agrícolas.
modelo_DBCA <- aov(Produccion ~ Anio + Tratamiento,
                   data = datos_largos)
anova_DBCA <- anova(modelo_DBCA)

formato_tabla(data.frame(Fuente = rownames(anova_DBCA), anova_DBCA,
                         check.names = FALSE), 4)
Fuente Df Sum Sq Mean Sq F value Pr(>F)
Anio 2 1.036410e+11 5.182048e+10 3.5887 0.0944
Tratamiento 3 3.115177e+12 1.038392e+12 71.9108 0.0000
Residuals 6 8.664002e+10 1.444000e+10 NA NA

5.2 Modelo factorial con bloques

Se presenta también el modelo que incorpora año, provincia, cultivo e interacción provincia–cultivo. En este conjunto balanceado, el efecto de tratamiento corresponde a la combinación de provincia y cultivo.

modelo_DBCA_factorial <- aov(Produccion ~ Anio + Provincia * Cultivo,
                             data = datos_largos)
anova_factorial_DBCA <- anova(modelo_DBCA_factorial)

formato_tabla(data.frame(Fuente = rownames(anova_factorial_DBCA),
                         anova_factorial_DBCA,
                         check.names = FALSE), 4)
Fuente Df Sum Sq Mean Sq F value Pr(>F)
Anio 2 1.036410e+11 5.182048e+10 3.5887 0.0944
Provincia 1 2.288206e+12 2.288206e+12 158.4630 0.0000
Cultivo 1 6.778391e+11 6.778391e+11 46.9418 0.0005
Provincia:Cultivo 1 1.491324e+11 1.491324e+11 10.3277 0.0183
Residuals 6 8.664002e+10 1.444000e+10 NA NA

5.3 Verificación de supuestos del DBCA

Se revisa la normalidad de los residuos y la homogeneidad de varianzas entre tratamientos. La prueba de Levene no evalúa el efecto de bloque; se utiliza aquí como diagnóstico de dispersión entre tratamientos.

residuos_DBCA <- residuals(modelo_DBCA)

shapiro_DBCA <- shapiro.test(residuos_DBCA)
levene_DBCA <- car::leveneTest(Produccion ~ Tratamiento,
                               data = datos_largos)

shapiro_DBCA

    Shapiro-Wilk normality test

data:  residuos_DBCA
W = 0.94537, p-value = 0.5706
levene_DBCA
par(mfrow = c(2, 2))
plot(modelo_DBCA)
Diagnóstico gráfico de los residuos del modelo DBCA.

Diagnóstico gráfico de los residuos del modelo DBCA. 

par(mfrow = c(1, 1))

5.4 Comparaciones múltiples y prueba de Friedman

Tukey se aplica a las medias ajustadas por el modelo aditivo. Friedman constituye una alternativa no paramétrica para un diseño en bloques, siempre que exista una observación por tratamiento dentro de cada bloque.

TukeyHSD(modelo_DBCA, which = "Tratamiento")
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Produccion ~ Anio + Tratamiento, data = datos_largos)

$Tratamiento
                                          diff         lwr        upr     p adj
Herrera – Arroz-Herrera – Maiz       -252379.0  -592026.86   87268.86 0.1439977
Los Santos – Maiz-Herrera – Maiz     1096306.3   756658.47 1435954.20 0.0001289
Los Santos – Arroz-Herrera – Maiz     398009.0    58361.14  737656.86 0.0256780
Los Santos – Maiz-Herrera – Arroz    1348685.3  1009037.47 1688333.20 0.0000387
Los Santos – Arroz-Herrera – Arroz    650388.0   310740.14  990035.86 0.0023237
Los Santos – Arroz-Los Santos – Maiz -698297.3 -1037945.20 -358649.47 0.0015907
friedman.test(Produccion ~ Tratamiento | Anio,
              data = datos_largos)

    Friedman rank sum test

data:  Produccion and Tratamiento and Anio
Friedman chi-squared = 9, df = 3, p-value = 0.02929

6 Precipitación y producción agrícola

6.1 Comportamiento de la precipitación

formato_tabla(
  datos |>
    select(Anio, Provincia, Precipitacion) |>
    arrange(Provincia, Anio),
  2
)
Anio Provincia Precipitacion
2014-2015 Herrera 940.5
2015-2016 Herrera 860.6
2016-2017 Herrera 1371.3
2014-2015 Los Santos 704.0
2015-2016 Los Santos 644.3
2016-2017 Los Santos 1232.1 
ggplot(datos,
       aes(x = Anio, y = Precipitacion,
           group = Provincia, color = Provincia)) +
  geom_line(linewidth = 1.1) +
  geom_point(size = 3) +
  scale_color_manual(values = c("Herrera" = "#276678",
                                "Los Santos" = "#d17a22")) +
  labs(
    title = "Precipitación anual por provincia",
    x = "Año agrícola", y = "Precipitación (mm)", color = "Provincia"
  ) +
  tema_informe +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))
Variación de la precipitación anual registrada en Herrera y Los Santos.

Variación de la precipitación anual registrada en Herrera y Los Santos.

6.2 Correlación y regresión exploratoria

La precipitación se repite para los dos cultivos de una misma provincia y año. Para evitar duplicar artificialmente las observaciones, las correlaciones se calculan por cultivo utilizando las seis combinaciones provincia–año disponibles.

Con solo seis observaciones por cultivo, las correlaciones y regresiones son exploratorias. No permiten establecer causalidad ni estimar con precisión el efecto de la precipitación sobre la producción.

resultados_cor <- lapply(levels(datos_largos$Cultivo), function(cultivo_actual) {
  tmp <- filter(datos_largos, Cultivo == cultivo_actual)

  pearson <- cor.test(tmp$Precipitacion, tmp$Produccion,
                      method = "pearson")
  spearman <- cor.test(tmp$Precipitacion, tmp$Produccion,
                       method = "spearman", exact = FALSE)

  data.frame(
    Cultivo = cultivo_actual,
    Pearson_r = unname(pearson$estimate),
    Pearson_p = pearson$p.value,
    Spearman_rho = unname(spearman$estimate),
    Spearman_p = spearman$p.value
  )
})

resultados_cor <- bind_rows(resultados_cor)
formato_tabla(resultados_cor, 4)
Cultivo Pearson_r Pearson_p Spearman_rho Spearman_p
Maiz -0.1650 0.7548 -0.0857 0.8717
Arroz -0.1639 0.7563 -0.1429 0.7872 
modelos_clima <- lapply(levels(datos_largos$Cultivo), function(cultivo_actual) {
  tmp <- filter(datos_largos, Cultivo == cultivo_actual)
  lm(Produccion ~ Precipitacion, data = tmp)
})
names(modelos_clima) <- levels(datos_largos$Cultivo)

summary(modelos_clima$Maiz)

Call:
lm(formula = Produccion ~ Precipitacion, data = tmp)

Residuals:
      1       2       3       4       5       6 
-736207 -507025 -296040  375268  278410  885593 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)
(Intercept)   1253350.8  1056348.3   1.186    0.301
Precipitacion    -355.4     1062.2  -0.335    0.755

Residual standard error: 687300 on 4 degrees of freedom
Multiple R-squared:  0.02722,   Adjusted R-squared:  -0.216 
F-statistic: 0.1119 on 1 and 4 DF,  p-value: 0.7548
summary(modelos_clima$Arroz)

Call:
lm(formula = Produccion ~ Precipitacion, data = tmp)

Residuals:
      1       2       3       4       5       6 
-288206 -391472 -233826  211664  144853  556987 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)   638371.6   627651.5   1.017    0.367
Precipitacion   -209.7      631.1  -0.332    0.756

Residual standard error: 408400 on 4 degrees of freedom
Multiple R-squared:  0.02687,   Adjusted R-squared:  -0.2164 
F-statistic: 0.1104 on 1 and 4 DF,  p-value: 0.7563
ggplot(datos_largos,
       aes(x = Precipitacion, y = Produccion, color = Provincia)) +
  geom_point(size = 3) +
  geom_smooth(method = "lm", se = FALSE, linewidth = 0.9) +
  facet_wrap(~ Cultivo, scales = "free_y",
             labeller = as_labeller(c(Maiz = "Maíz", Arroz = "Arroz"))) +
  scale_y_continuous(labels = scales::label_comma()) +
  labs(
    title = "Precipitación y producción agrícola",
    subtitle = "Seis observaciones provincia–año por cultivo",
    x = "Precipitación anual (mm)",
    y = "Producción (quintales)",
    color = "Provincia"
  ) +
  tema_informe
Relación exploratoria entre precipitación anual y producción, separada por cultivo.

Relación exploratoria entre precipitación anual y producción, separada por cultivo.

7 Comparación de los modelos

Se comparan los cuadrados medios del error de ambos modelos. Una reducción del cuadrado medio del error al incorporar bloques puede indicar que el año explica parte de la variabilidad, pero no basta por sí sola para declarar que el DBCA es superior. También deben considerarse los grados de libertad, los supuestos y la naturaleza observacional de los datos.

MSE_DCA <- deviance(modelo_DCA) / df.residual(modelo_DCA)
MSE_DBCA <- deviance(modelo_DBCA) / df.residual(modelo_DBCA)

comparacion <- data.frame(
  Modelo = c("DCA", "DBCA"),
  Grados_libertad_residuales = c(df.residual(modelo_DCA),
                                 df.residual(modelo_DBCA)),
  Cuadrado_medio_error = c(MSE_DCA, MSE_DBCA),
  Desviacion_residual = c(sigma(modelo_DCA), sigma(modelo_DBCA))
)

formato_tabla(comparacion, 3)
Modelo Grados_libertad_residuales Cuadrado_medio_error Desviacion_residual
DCA 8 23785122435 154224.3
DBCA 6 14440004153 120166.6

8 Conclusiones y limitaciones

Las conclusiones deben redactarse después de revisar los valores de F, los valores p, los intervalos y las comparaciones múltiples generadas al ejecutar el documento. No se anticipa aquí qué tratamiento es estadísticamente diferente.

Aspectos que deben considerarse al interpretar los resultados:

  • Producción y rendimiento: los valores suministrados están expresados en quintales totales. Si no se dispone del área cosechada, no deben denominarse rendimiento por hectárea.
  • Registros históricos: las observaciones no son unidades experimentales asignadas aleatoriamente. Por tanto, los análisis DCA y DBCA son aproximaciones de comparación estadística y no prueban causalidad.
  • Tamaño de muestra: existen tres años y seis registros provincia–año para cada cultivo. La capacidad de detectar diferencias y validar supuestos es limitada.
  • Precipitación: la asociación con producción puede estar influida por superficie sembrada, tecnología, variedades, manejo agrícola, plagas y otros factores no incluidos.
  • DBCA: con tres bloques, la estimación del efecto del año y del error residual tiene pocos grados de libertad.

9 Archivos exportables

El siguiente bloque permite guardar las tablas principales en formato CSV dentro del directorio de trabajo de R.

write.csv(datos_largos, "Base_Datos_Agricola.csv", row.names = FALSE)
write.csv(descriptivos, "Estadistica_Descriptiva.csv", row.names = FALSE)
write.csv(as.data.frame(anova_DCA), "Resultados_ANOVA_DCA.csv")
write.csv(as.data.frame(anova_DBCA), "Resultados_ANOVA_DBCA.csv")
write.csv(resultados_cor, "Correlaciones_Precipitacion_Produccion.csv",
          row.names = FALSE)

cat("Análisis finalizado. Las tablas CSV se guardaron en:\n")
Análisis finalizado. Las tablas CSV se guardaron en:
getwd()
[1] "/cloud/project"

Fin del informe.