Apéndice I — Análisis de rendimiento – Métodos basados en imágenes como alternativa para la evaluación de maleza en cultivo de trigo en el Sur de Sonora

Para determinar el rendimiento se muestreo manualmente 2 metros cuadrados del centro de cada una de las 72 subparcelas.

I.0.1 Descriptivo

# Leer los datos desde el archivo Excel
data <- read_excel("csv/correlacionRendimiento.xlsx")

# Mostrar una vista previa de los datos
head(data)

# A tibble: 6 × 28
  cambioNoPlantas cambioAltura Tratamiento bloque Parcela Rendimiento
            <dbl>        <dbl> <chr>        <dbl> <chr>         <dbl>
1             -34           15 T1               1 A1             4130
2              -6            0 T1               4 A4             5270
3             -37           20 T1               7 B8             6050
4             -40           25 T1              10 B10            5700
5             -33            0 T1               2 C3             5450
6             -13           15 T1               5 C6             4960
# ℹ 22 more variables: noPlantas0dda <dbl>, noPlantas14dda <dbl>,
#   noPlantas19dda <dbl>, noPlantas28dda <dbl>, hMax0dda <dbl>,
#   hMax14dda <dbl>, hMax19dda <dbl>, hMax28dda <dbl>, Columna <chr>,
#   Fila <dbl>, tiempo0dda <dbl>, tiempo14dda <dbl>, tiempo19dda <dbl>,
#   tiempo28dda <dbl>, tiempo0ddaBI <dbl>, cober0dda <dbl>, cober14dda <dbl>,
#   cober19dda <dbl>, cober28dda <dbl>, cloro14dda <dbl>, cloro19dda <dbl>,
#   cloro28dda <dbl>

# Convertir la columna 'Tratamiento' a factor
data$Tratamiento <- as.factor(data$Tratamiento)

# Análisis descriptivo
summary_stats <- data %>% 
  group_by(Tratamiento) %>% 
  summarise(n = n(),
            media = mean(Rendimiento, na.rm = TRUE),
            mediana = median(Rendimiento, na.rm = TRUE),
            sd = sd(Rendimiento, na.rm = TRUE),
            varianza = var(Rendimiento, na.rm = TRUE),
            min = min(Rendimiento, na.rm = TRUE),
            max = max(Rendimiento, na.rm = TRUE))


print(summary_stats)

# A tibble: 6 × 8
  Tratamiento     n media mediana    sd varianza   min   max
  <fct>       <int> <dbl>   <dbl> <dbl>    <dbl> <dbl> <dbl>
1 T1             12 4929.    5025  727.  528572.  3820  6050
2 T2             12 4921.    4855  452.  203855.  4315  5870
3 T3             12 4672.    4395 1206. 1453852.  2420  6550
4 T4             12 5058.    5005  788.  621397.  3720  6350
5 TR             12 5117.    5220  851.  723479.  3030  6320
6 TS             12 3738.    4205 1383. 1912002.   290  5170

I.0.2 Datos sin transformar

I.0.2.1 Análisis de Varianza

# Análisis estadístico
# Análisis de Varianza (ANOVA) considerando el efecto de Tratamiento y bloque
model <- aov(Rendimiento ~ Tratamiento + bloque, data = data)
anova_results <- Anova(model)

# Mostrar los resultados del ANOVA
print(anova_results)

Anova Table (Type II tests)

Response: Rendimiento
              Sum Sq Df F value   Pr(>F)   
Tratamiento 15858243  5  3.6503 0.005679 **
bloque       3397657  1  3.9104 0.052232 . 
Residuals   56477066 65                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

I.0.2.2 Análisis de normalidad de los residuos utilizando Shapiro-Wilk

shapiro_test <- shapiro.test(residuals(model))
print(shapiro_test)


    Shapiro-Wilk normality test

data:  residuals(model)
W = 0.9449, p-value = 0.003365

I.0.2.3 Análisis de homogeneidad de varianzas utilizando la prueba de Levene

levene_test <- car::leveneTest(Rendimiento ~ Tratamiento, data = data)
print(levene_test)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  5  1.1804 0.3283
      66

I.0.2.4 Análisis de independencia de los residuos

plot(residuals(model) ~ fitted(model))
abline(h = 0, col = "red")

I.0.2.5 Comparaciones múltiples

# Comparaciones múltiples usando el test de Tukey
tukey_results <- HSD.test(model, "Tratamiento")

# Mostrar los resultados del test de Tukey
print(tukey_results)

$statistics
   MSerror Df     Mean       CV      MSD
  868877.9 65 4739.097 19.66907 1117.439

$parameters
   test      name.t ntr StudentizedRange alpha
  Tukey Tratamiento   6         4.152742  0.05

$means
   Rendimiento       std  r       se  Min  Max    Q25  Q50    Q75
T1    4929.167  727.0296 12 269.0845 3820 6050 4182.5 5025 5462.5
T2    4921.250  451.5032 12 269.0845 4315 5870 4565.0 4855 5235.0
T3    4671.667 1205.7577 12 269.0845 2420 6550 4022.5 4395 5555.0
T4    5058.333  788.2874 12 269.0845 3720 6350 4572.5 5005 5567.5
TR    5116.667  850.5756 12 269.0845 3030 6320 4725.0 5220 5640.0
TS    3737.500 1382.7517 12 269.0845  290 5170 3530.0 4205 4557.5

$comparison
NULL

$groups
   Rendimiento groups
TR    5116.667      a
T4    5058.333      a
T1    4929.167      a
T2    4921.250      a
T3    4671.667     ab
TS    3737.500      b

attr(,"class")
[1] "group"

I.0.3 Transformación de Box-Cox

# Cargar las librerías necesarias
library(MASS)
library(car)

Warning: package 'car' was built under R version 4.2.3

Loading required package: carData


Attaching package: 'car'

The following object is masked from 'package:dplyr':

    recode

The following object is masked from 'package:purrr':

    some

# Leer los datos desde el archivo Excel
data <- read_excel("csv/correlacionRendimiento.xlsx")

# Convertir la columna 'Tratamiento' a factor
data$Tratamiento <- as.factor(data$Tratamiento)
data$bloque <- as.factor(data$bloque)


# Aplicar la transformación de Box-Cox
boxcox_results <- boxcox(Rendimiento ~ Tratamiento + bloque, data = data, lambda = seq(-2, 2, 0.1))

# Encontrar el mejor lambda
best_lambda <- boxcox_results$x[which.max(boxcox_results$y)]
print(paste("Mejor lambda:", best_lambda))

[1] "Mejor lambda: 2"

# Transformar los datos usando el mejor lambda
if (best_lambda == 0) {
  data$transformed_Rendimiento <- log(data$Rendimiento)
} else {
  data$transformed_Rendimiento <- (data$Rendimiento^best_lambda - 1) / best_lambda
}

# Ajustar el modelo con los datos transformados
model_transformed <- aov(transformed_Rendimiento ~ Tratamiento + bloque, data = data)

# Verificación de supuestos con los datos transformados
shapiro_test_transformed <- shapiro.test(residuals(model_transformed))
levene_test_transformed <- car::leveneTest(transformed_Rendimiento ~ Tratamiento, data = data)

# Resultados del modelo y verificación de supuestos
anova_results_transformed <- Anova(model_transformed)
print(anova_results_transformed)

Anova Table (Type II tests)

Response: transformed_Rendimiento
                Sum Sq Df F value   Pr(>F)   
Tratamiento 2.4395e+14  5  3.9999 0.003633 **
bloque      4.0746e+14 11  3.0367 0.003109 **
Residuals   6.7088e+14 55                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

print(shapiro_test_transformed)


    Shapiro-Wilk normality test

data:  residuals(model_transformed)
W = 0.98134, p-value = 0.3625

print(levene_test_transformed)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  5   0.993 0.4289
      66

# Análisis de independencia de los residuos transformados
plot(residuals(model_transformed) ~ fitted(model_transformed))
abline(h = 0, col = "red")

# Comparaciones múltiples usando el test de Tukey en datos transformados
tukey_results_transformed <- HSD.test(model_transformed, "Tratamiento")
print(tukey_results_transformed)

$statistics
       MSerror Df     Mean       CV     MSD
  1.219784e+13 55 11755444 29.70998 4209793

$parameters
   test      name.t ntr StudentizedRange alpha
  Tukey Tratamiento   6         4.175514  0.05

$means
   transformed_Rendimiento     std  r      se       Min      Max      Q25
T1                12390604 3541710 12 1008209 7296199.5 18301250  8747112
T2                12202784 2274303 12 1008209 9309612.0 17228450 10419650
T3                11578583 5728882 12 1008209 2928199.5 21451250  8101737
T4                13078175 3960655 12 1008209 6919199.5 20161250 10455987
TR                13421733 4018910 12 1008209 4590449.5 19971200 11162850
TS                 7860787 3956766 12 1008209   42049.5 13364450  6232850
        Q50      Q75
T1 12626325 14919687
T2 11786125 13702650
T3  9658025 15448850
T4 12525125 15500637
TR 13653000 15904799
TS  8853025 10388787

$comparison
NULL

$groups
   transformed_Rendimiento groups
TR                13421733      a
T4                13078175      a
T1                12390604      a
T2                12202784      a
T3                11578583     ab
TS                 7860787      b

attr(,"class")
[1] "group"

I.0.4 Resumen gráfico

Código

library(ggdist)
library(ggtext)

# Leer los datos desde el archivo Excel
data <- read_excel("csv/correlacionRendimiento.xlsx")

# Mostrar una vista previa de los datos
#head(data)

# Convertir la columna 'Tratamiento' a factor
data$Tratamiento <- as.factor(data$Tratamiento)

# Análisis descriptivo
summary_stats <- data %>% 
  group_by(Tratamiento) %>% 
  summarise(n = n(),
            media = mean(Rendimiento, na.rm = TRUE),
            mediana = median(Rendimiento, na.rm = TRUE),
            sd = sd(Rendimiento, na.rm = TRUE),
            varianza = var(Rendimiento, na.rm = TRUE),
            min = min(Rendimiento, na.rm = TRUE),
            max = max(Rendimiento, na.rm = TRUE))



# Preparar los datos del test de Tukey
tukey_groups <- tukey_results_transformed$groups
tukey_groups$Tratamiento <- rownames(tukey_groups)
rownames(tukey_groups) <- NULL

# Unir los resultados del test de Tukey con los summary_stats
summary_stats <- summary_stats %>%
  left_join(tukey_groups, by = "Tratamiento")

# Ordenar el eje Y de manera descendente
data$Tratamiento <- factor(data$Tratamiento, levels = summary_stats$Tratamiento[order(summary_stats$media)])

graficoRendimientoTUKEY <- data %>%
  ggplot(aes(x=Rendimiento, y = Tratamiento, fill = Tratamiento)) +
  geom_boxplot(width = 0.1) +
  geom_dots(
    side = 'bottom',
    position = position_nudge(y = -0.075),
    height = 0.55
  ) +
  stat_slab(
    position = position_nudge(y = 0.075),
    height = 0.55,
    trim = FALSE,
    alpha = 0.6
  ) +
  labs(
    title= "",
    y = element_blank(),
    x = "Rendimiento (kg/ha)"
  ) +
  theme_minimal(
    base_size = 20)+
  theme(
    text = element_text(family = "montserrat"),
    panel.grid.minor = element_blank(),
    legend.position = 'none'
  )+
  geom_text(data = summary_stats, aes(x = 7200, y = Tratamiento, label = paste0(round(media, 2), "  ", groups)), 
            hjust = -0.2, size = 3, color = "#0a231fff", fontface = "bold", hjust = 1)

graficoRendimientoTUKEY

#ggsave("./Plots/graficoRendimientoTUKEY.png", plot = graficoRendimientoTUKEY, width=1920, height=1080, units = "px")

Figura I.1

I.0.5 Análisis no paramétrico

# Test de Kruskal-Wallis
kruskal_test <- kruskal.test(Rendimiento ~ Tratamiento, data = data)
print(kruskal_test)


    Kruskal-Wallis rank sum test

data:  Rendimiento by Tratamiento
Kruskal-Wallis chi-squared = 13.388, df = 5, p-value = 0.02

I.0.5.1 Test de Dunn

library(FSA)

Warning: package 'FSA' was built under R version 4.2.3

Registered S3 methods overwritten by 'FSA':
  method       from
  confint.boot car 
  hist.boot    car

## FSA v0.9.5. See citation('FSA') if used in publication.
## Run fishR() for related website and fishR('IFAR') for related book.


Attaching package: 'FSA'

The following object is masked from 'package:car':

    bootCase

dunn_test <- dunnTest(Rendimiento ~ Tratamiento , data = data, method = "bonferroni")
print(dunn_test)

Dunn (1964) Kruskal-Wallis multiple comparison

  p-values adjusted with the Bonferroni method.

   Comparison          Z     P.unadj      P.adj
1     T1 - T2  0.1365546 0.891382871 1.00000000
2     T1 - T3  0.8485892 0.396109891 1.00000000
3     T2 - T3  0.7120347 0.476443315 1.00000000
4     T1 - T4 -0.3560173 0.721827588 1.00000000
5     T2 - T4 -0.4925719 0.622315092 1.00000000
6     T3 - T4 -1.2046066 0.228355216 1.00000000
7     T1 - TR -0.6730191 0.500935165 1.00000000
8     T2 - TR -0.8095737 0.418185256 1.00000000
9     T3 - TR -1.5216083 0.128107254 1.00000000
10    T4 - TR -0.3170017 0.751242289 1.00000000
11    T1 - TS  2.5311369 0.011369346 0.17054020
12    T2 - TS  2.3945823 0.016639308 0.24958963
13    T3 - TS  1.6825476 0.092462693 1.00000000
14    T4 - TS  2.8871542 0.003887436 0.05831154
15    TR - TS  3.2041560 0.001354591 0.02031886

I.0.6 Boxplot de rendimientos por tratamiento

Mostrar el código.

# Extraer los grupos de Tukey
tukey_groups <- as.data.frame(tukey_results_transformed$groups)


# Añadir la columna 'Tratamiento' a 'tukey_groups' para mapearlo en el gráfico
tukey_groups$Tratamiento <- rownames(tukey_groups)

# Ordenar el eje Y por rendimiento medio (de mayor a menor)
data$Tratamiento <- fct_reorder(data$Tratamiento, -tukey_groups$transformed_Rendimiento[match(data$Tratamiento, tukey_groups$Tratamiento)])

# Parámetros para las líneas verticales
lineas <- data.frame(
  x = c(6600, 5750, 7070),
  label = c("Historico", "2024", "2023")
)


# Paleta de colores personalizada
colores <- c("TR" = "#AE8F5F", "T4" = "#00AA00", "T1" = "#007ACC", "T2" = "#FF8C00", "T3" = "#CC3333", "TS" = "#9370DB")



# Visualización gráfica con las letras de Tukey


# Boxplot de rendimientos por tratamiento con ejes invertidos y letras de Tukey
boxplot_tratamiento <- ggplot(data, aes(y = Tratamiento, x = Rendimiento, fill = Tratamiento)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Boxplot de Rendimiento por Tratamiento", x = "Rendimiento (t/ha)", y = "Tratamiento") +
  theme(legend.position = "right") + # Colocar la leyenda a la derecha
  geom_text(data = tukey_groups, aes(x = max(data$Rendimiento) + 700, y = Tratamiento, label = groups), 
            hjust = 0, color = "black", size = 5, position = position_dodge(width = 0.9))+ # Colocar la leyenda a la derecha
  scale_x_continuous(expand = c(0, 0), breaks = seq(0, max(data$Rendimiento), by = 1000)) +
  coord_cartesian(xlim = c(0, max(data$Rendimiento) + 1000), ylim = c(0.5, 6.5))+
  geom_vline(data = lineas, aes(xintercept = x), linetype = "dashed", color = "blue", size = 1) +
  annotate("text", x = lineas$x, y = 6.4, label = lineas$label, vjust = -0.5, color = "black", size = 3)+
  scale_fill_manual(values = colores)




# Mostrar el boxplot de tratamientos con ejes invertidos y letras de Tukey
print(boxplot_tratamiento)

Figura I.2: RendimientoBoxplot.

I.0.7 Gráfico de barras para el promedio de rendimiento por tratamiento

Mostrar el código.

# Gráfico de barras para el promedio de rendimiento por tratamiento

barplot_tratamiento <- data %>%
  group_by(Tratamiento) %>%
  summarise(media_rendimiento = mean(Rendimiento, na.rm = TRUE)) %>%
  mutate(Tratamiento = fct_reorder(Tratamiento, -media_rendimiento)) %>%  # Reordenar factores de mayor a menor
  ggplot(aes(x = media_rendimiento, y = Tratamiento, fill = Tratamiento)) +
  geom_bar(stat = "identity") +
  theme_minimal() +
  labs(title = "Promedio de Rendimiento por Tratamiento", x = "Rendimiento promedio (kg/ha)", y = "Tratamiento") +
  theme(legend.position = "right") + # Colocar la leyenda a la derecha
  scale_x_continuous(expand = c(0, 0), breaks = seq(0, max(data$Rendimiento), by = 1000)) +
  coord_cartesian(xlim = c(0, max(data$Rendimiento) + 1000), ylim = c(0.5, 6.5)) + # Ajustar el límite del eje X y Y
  geom_text(aes(x = media_rendimiento - 300, y = Tratamiento, label = tukey_groups$groups[match(Tratamiento, tukey_groups$Tratamiento)]), 
            hjust = 1, color = "black", size = 4) + # Colocar las letras Tukey dentro de las barras, cerca del eje Y
  geom_vline(data = lineas, aes(xintercept = x), linetype = "dashed", color = "blue", size = 1) +
  annotate("text", x = lineas$x, y = 6.4, label = lineas$label, vjust = -0.5, color = "black", size = 3)+
  scale_fill_manual(values = colores)



# Mostrar el gráfico de barras de tratamientos con letras de Tukey
print(barplot_tratamiento)

Figura I.3: RendimientoBarras.

I.0.8 Gráfico de barras para el promedio de rendimiento por tratamiento (Porcentaje respecto al testigo regional)

Mostrar el código.

# Calcular el rendimiento promedio por tratamiento
promedios <- data %>%
  group_by(Tratamiento) %>%
  summarise(media_rendimiento = mean(Rendimiento, na.rm = TRUE))

# Obtener el rendimiento promedio del tratamiento TR
rendimiento_tr <- promedios %>% filter(Tratamiento == "TR") %>% pull(media_rendimiento)

# Calcular el porcentaje relativo al tratamiento TR
promedios <- promedios %>%
  mutate(porcentaje_relativo = (media_rendimiento / rendimiento_tr) * 100)

# Ordenar el eje Y por porcentaje relativo (de mayor a menor)
promedios$Tratamiento <- fct_reorder(promedios$Tratamiento, -promedios$porcentaje_relativo)

# Crear el gráfico de barras
barplot_porcentaje <- ggplot(promedios, aes(x = porcentaje_relativo, y = Tratamiento, fill = Tratamiento)) +
  geom_bar(stat = "identity") +
  theme_minimal() +
  labs(title = "Rendimiento Relativo por Tratamiento (Referencia TR = 100%)", 
       x = "Porcentaje relativo a TR (%)", y = "Tratamiento") +
  theme(legend.position = "none") + # Ocultar la leyenda ya que no es necesaria
  geom_text(aes(label = round(porcentaje_relativo, 1)), 
            hjust = -0.1, color = "black", size = 4) + # Mostrar el porcentaje en las barras
  scale_x_continuous(expand = c(0, 0), limits = c(0, max(promedios$porcentaje_relativo) * 1.1)) +
  scale_fill_manual(values = colores)

# Mostrar el gráfico de barras
print(barplot_porcentaje)

Figura I.4: RendimientoBarrasPorcentaje.

I.0.9 Comparación respecto a la producción histórica

Figura I.5: Pérdidas de producción debido a diversos factores