---
toc: true
toc-location: right
title: "Análisis de datos"
---

```{r carga}
#| echo: FALSE
#| warning: FALSE
#| eval: TRUE

library(fontawesome)
library(tidyverse)
library(hrbrthemes)
library(ggplot2)
library("readxl")
library("plotrix")
library(ggpubr)
library(rstatix)
library(Skillings.Mack)
library(agricolae)
library(MASS)
library("gridExtra") 
library(extrafont)
loadfonts(device = "win")
library(showtext)
library(AgroR)
library(ggpmisc)
library(dplyr)
library(tidyr)
library(gridExtra)
library(forcats)  # Para reordenar factores

font_add_google("Montserrat", "montserrat")




```

# Análisis de datos

A continuación se presentan los principales resultados del proyecto.

## Análisis previo

::: justify
Dado que las distintas especies de maleza que se pueden encontrar en una parcela de cultivo no se distribuyen de manera aleatoria, e incluso existe una fuerte evidencia de agrupamiento [@riepma1963], antes de someter los resultados a un análisis estadístico convencional, se se realizó el siguiente análisis para determinar la frecuencia de cada una de las especies y su distribución dentro de la parcela experimental.
:::

### Nivel de infestación de maleza

::: justify
Este análisis proporciona un resumen visual de la frecuencia inicial de las distintas especies de maleza en las subparcelas y es fundamental para entender la distribución y densidad de las especies de malezas presentes.
:::

#### Conteo previo

::: justify
El siguiente gráfico muestra la frecuencia de aparición de las distintas especies de maleza encontradas en las 72 subparcelas en el conteo previo a la aplicación herbicida. De esta manera podemos determinar cuales especies de maleza se presentan con la suficiente frecuencia en nuestra parcela experimental como para ser sujetas a un análisis estadístico que nos permita realizar una comparación entre tratamientos. Se puede observar que solamente las especies de *C. murale*, *A. palmeri*, *C. acutifolia* y el complejo de maleza se presentaron con suficiente frecuencia para ser analizados (@fig-previo).
:::

```{r analisisprevio}
#| echo: TRUE
#| warning: FALSE
#| eval: TRUE
#| code-fold: true
#| code-summary: "Mostrar el código."


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

data_previo <- filter(data, data$conteo == "Previo")


# Función para preparar los datos y crear el gráfico
create_species_plot <- function(species_name, data, total_parcelas) {
  # Filtrar los datos por especie
  species_data <- data %>% filter(especie == species_name)
  
  # Contar ocurrencias por parcela
  species_counts <- species_data %>%
    group_by(parcela) %>%
    summarise(Count = n()) %>%
    ungroup()
  
  # Crear una matriz vacía de 12x6
  matrix_data <- matrix(0, nrow = 12, ncol = 6, dimnames = list(1:12, LETTERS[1:6]))
  
  # Rellenar la matriz con los datos contados
  for (i in 1:nrow(species_counts)) {
    parcela <- species_counts$parcela[i]
    row_num <- as.numeric(substr(parcela, 2, 3)) # Corregir para extraer correctamente el número de fila (p.ej., 'A01' -> 1, 'A12' -> 12)
    col_num <- match(substr(parcela, 1, 1), LETTERS) # Extraer letra de columna de la parcela (p.ej., 'A01' -> A)
    matrix_data[row_num, col_num] <- species_counts$Count[i]
  }
  
  # Convertir la matriz en un dataframe largo para ggplot2
  df_long <- as.data.frame(as.table(matrix_data))
  colnames(df_long) <- c("Row", "Column", "Count")
  df_long$Count <- as.numeric(as.character(df_long$Count))
  
  # Invertir el orden de las filas
  df_long$Row <- factor(df_long$Row, levels = 1:12)
  
  # Calcular el total de ocurrencias y el porcentaje
  total_occurrences <- sum(df_long$Count)
  percentage <- (total_occurrences / total_parcelas) * 100
  
  # Crear el gráfico
  p <- ggplot(df_long, aes(x = Column, y = Row)) +
    geom_tile(aes(fill = Count > 0), color = "black") + # Ajustar la escala de color
    scale_fill_manual(values = c("FALSE" = "white", "TRUE" = "green")) + # Solo pintar verde las celdas con valores > 0
    geom_text(aes(label = ifelse(Count > 0, Count, "")), vjust = 0.5, color = "black") + # Mostrar solo los valores > 0
    labs(
      title = paste0(species_name, " (", total_occurrences, "/", total_parcelas, ", ", round(percentage, 2), "%)"),
      x = "Columnas",
      y = "Filas"
    ) +
    theme_minimal() +
    theme(
      text = element_text(family = "montserrat"),
      axis.text = element_text(size = 10, color = "#0a231fff"),
      axis.title = element_text(size = 12, face = "bold", color = "#ae8f5fff"),
      plot.title = element_text(size = 13, face = "bold", color = "#ae8f5fff"),
      legend.position = "none"
    ) +
    coord_fixed(ratio = 7.5 / 8)  # Ajustar la relación de aspecto
  
  return(p)
}

# Crear gráficos para las primeras 3 especies como ejemplo
species_list <- unique(data$especie)
total_parcelas <- length(unique(data$parcela)) # Corregir el total de parcelas basado en el número de parcelas únicas


# Graficar especies individuales directamente
InicialComplejo <- create_species_plot(species_list[4], data_previo, total_parcelas) # Complejo
InicialChualC <- create_species_plot(species_list[2], data_previo, total_parcelas) # Chual C
InicialQuelite <- create_species_plot(species_list[1], data_previo, total_parcelas) # Quelite
InicialTomatillo <- create_species_plot(species_list[3], data_previo, total_parcelas) # Tomatillo
InicialBorraja <- create_species_plot(species_list[5], data_previo, total_parcelas) # Borraja
InicialMalva <- create_species_plot(species_list[6], data_previo, total_parcelas) # Malva
InicialChualB <- create_species_plot(species_list[7], data_previo, total_parcelas) # Chual B
InicialColaR <- create_species_plot(species_list[8], data_previo, total_parcelas) # Cola de Rata
InicialTrebol <- create_species_plot(species_list[9], data_previo, total_parcelas) # Trebol
InicialVerdolaga <- create_species_plot(species_list[10], data_previo, total_parcelas) # Verdolaga
InicialMalvaP <- create_species_plot(species_list[11], data_previo, total_parcelas) # Malva P


library(patchwork)


InicialCombinacion <- InicialComplejo + InicialChualC + 
  InicialQuelite + 
  InicialTomatillo +
  InicialMalva + 
  InicialChualB +
  InicialBorraja +
  InicialColaR +
  InicialVerdolaga + 
  InicialTrebol + 
  InicialMalvaP +
  plot_layout(ncol = 4) 



#ggsave("./Plots/InicialCombinacion.png", plot = InicialCombinacion, width=10, height=14)



```

![Presencia de las distintas especies de maleza en el conteo previo a la aplicación de herbicida.](plots/InicialCombinacion.png){#fig-previo group="gallery"}

#### Desarrollo durante el experimento

::: justify
Este análisis rastrea el progreso de la infestación de maleza a lo largo del tiempo en las parcelas experimentales. Incluye datos recolectados en diferentes momentos del experimento. Los cambios en la densidad y crecimiento de las malezas, representados en gráficos temporales, indican cómo los tratamientos afectan la infestación a lo largo del tiempo. Es útil para identificar tendencias y patrones de crecimiento.

La @fig-dev_Complejo_TS muestra el comportamiento del complejo de maleza en el testigo sin aplicación. Si observamos los valores del número de plantas (a la izquierda del guión en cada recuadro), podemos notar que la variación desde el conteo previo hasta los 28 dda se comportó de manera aleatoria, presentándose casos en los que el número de plantas de maleza fue mayor y en otros casos menor, con un promedio de 37 en el conteo previo y 30 en el conteo a los 28 dda, es decir, 19% menos. Por otra parte, el cambio en la altura máxima, durante el mismo periodo, presentó una tendencia mucho más coherente, aumentando desde 17 cm de promedio en el conteo previo, hasta 59 cm promedio en el conteo final; que representa un incremento del 253%.
:::

![Desarrollo de la presencia del complejo de maleza durante el experimento en las parcelas del testigo sin aplicación.](plots/dev_Complejo_TS.png){#fig-dev_Complejo_TS group="gallery"}

::: justify
La @fig-dev_Complejo_TR y @fig-dev_Complejo_T1 ejemplifican el comportamiento observado en las parcelas con aplicación herbicida. En este caso podemos observar que el número de plantas de maleza disminuyo moderadamente, y que la altura máxima se comportó de manera aleatoria. En las parcelas con aplicación del testigo regional, se contaron, en promedio, 35 plantas en el conteo previo y 13 plantas en el conteo a los 28 dda, una disminución del 62%. Si utilizamos la ecuación de @abbott1925, y ajustamos de acuerdo al 19% de disminución en el número de plantas del testigo sin aplicación, esto representa solamente un 53% de control. Por otro lado, la altura máxima promedio de las plantas en el testigo regional incrementó 3 cm (16%).
:::

![Desarrollo de la presencia del complejo de maleza durante el experimento en las parcelas del testigo regional.](plots/dev_Complejo_TR.png){#fig-dev_Complejo_TR group="gallery"}

![Desarrollo de la presencia del complejo de maleza durante el experimento en las parcelas del tratamiento 1.](plots/dev_Complejo_T1.png){#fig-dev_Complejo_T1 group="gallery"}

::: justify
De lo anterior, podemos suponer que, bajo las condiciones del experimento, el efecto de la aplicación de los tratamientos herbicidas, respecto al testigo sin aplicación, se expresó de manera más evidente en el cambio en altura máxima que en el cambio en el número de plantas.
:::

### Correlación de variables evaluadas y rendimiento

::: justify
Para determinar si las variables evaluadas tuvieron un efecto sobre el rendimiento de trigo, y en que magnitud lo hicieron, se realizaron los siguientes análisis de correlación. El objetivo de este análisis es validar que las variables medidas son las correctas para explicar la variabilidad en el rendimiento. Esto es crucial, ya que demuestra que las variables de control evaluadas, el cambio en el número de plantas y la altura máxima, son suficientes para explicar las diferencias en el rendimiento observadas entre las subparcelas. Si estas variables explican una parte significativa de la variabilidad en el rendimiento, se puede proceder con confianza a los análisis comparativos entre tratamientos, sabiendo que se están utilizando las variables adecuadas para evaluar el impacto de los tratamientos sobre el rendimiento y el control de maleza.
:::

::: {.callout-note collapse="true"}
En el contexto de análisis estadístico, es importante entender las diferencias entre el coeficiente de correlación de Pearson (r) y el coeficiente de determinación (R cuadrada).

1.  Coeficiente de Correlación de Pearson (r):

-   El coeficiente de correlación de Pearson, denotado como r, es una medida de la fuerza y la dirección de la relación lineal entre dos variables cuantitativas.

-   r puede tomar valores entre -1 y 1. Los rangos y sus interpretaciones generales son los siguientes:

    -   r=1: Correlación positiva perfecta. Las variables aumentan o disminuyen juntas de manera lineal perfecta.
    -   0.7≤r\<1: Correlación positiva fuerte.
    -   0.4≤r\<0.7: Correlación positiva moderada.
    -   0.2≤r\<0.4: Correlación positiva débil.
    -   0≤r\<0.2: Correlación positiva muy débil o inexistente.
    -   r=0: No hay correlación lineal entre las variables.
    -   −0.2≤r\<0: Correlación negativa muy débil o inexistente.
    -   −0.4≤r\<−0.2: Correlación negativa débil.
    -   −0.7≤r\<−0.4: Correlación negativa moderada.
    -   −1≤r\<−0.7: Correlación negativa fuerte.
    -   r=−1: Correlación negativa perfecta. Una variable aumenta mientras la otra disminuye de manera lineal perfecta.

Es importante recordar que el coeficiente de correlación de Pearson solo mide relaciones lineales y no implica causalidad entre las variables.

-   El coeficiente de correlación de Pearson no tiene unidades y no está limitado a la relación lineal; también puede capturar relaciones no lineales en algunos casos.

2.  Coeficiente de Determinación (R cuadrada):

-   El coeficiente de determinación, denotado como R2, es una medida que indica cuánta variación en la variable de respuesta (variable dependiente) es explicada por el modelo de regresión lineal.

-   R2 se calcula como el cuadrado del coeficiente de correlación de Pearson (r), por lo que R2=r2.

-   R2 varía de 0 a 1, donde:

    -   R2 = 0 indica que el modelo de regresión no explica ninguna variabilidad en la variable dependiente.
    -   R2 = 1 indica que el modelo de regresión explica toda la variabilidad en la variable dependiente.

-   R2 se interpreta como la proporción de la variabilidad total de la variable dependiente que es explicada por el modelo de regresión lineal.

En resumen:

-   El coeficiente de correlación de Pearson (r) mide la fuerza y la dirección de la relación lineal entre dos variables.
-   El coeficiente de determinación (R cuadrada) es una medida de cuánta variabilidad en la variable dependiente es explicada por el modelo de regresión lineal.
:::

#### Cambio en el número de plantas *vs* Rendimiento

::: justify
Los datos del cambio en el número de plantas *vs* rendimiento (kg/ha) mostraron una correlación significativa (*P* = 0.011), con coeficiente de correlación de Pearson de -0.299, que se considera una correlación negativa débil; y un coeficiente de determinación de 0.089, lo cual indica que aproximadamente 9% de la variabilidad del rendimiento es explicada por el cambio en el número de plantas.
:::

```{r CorrelacionRendimientoNoplantas, dev = "svg", dev.args=list(bg="transparent")}
#| label: fig-CorrVarRendimientoControl
#| fig-cap: "Correlación entre el "
#| echo: TRUE
#| warning: FALSE
#| eval: TRUE
#| lightbox:
#|   group: gallery
#| code-fold: true
#| code-summary: "Mostrar el código."

correlacionRendimiento = read_excel(path = "csv/correlacionRendimiento.xlsx")


correlacionRendimiento_ggControl <- ggplot(correlacionRendimiento, aes(cambioNoPlantas, Rendimiento)) + 
  geom_point() +
  geom_smooth(method = 'lm', formula = y ~ x, se = TRUE, fill = "blue", color = "blue", fullrange = TRUE, na.rm = TRUE, alpha = 0.10, linetype = "dashed") +
  theme_classic() + 
  labs(title = expression("Rendimiento vs"), subtitle = "Cambio en número de plantas") + 
  xlab("Cambio en número de plantas") + 
  ylab("Rendimiento") + 
  theme(
    text = element_text(family = "montserrat"), 
    axis.text = element_text(size = 9), 
    axis.title = element_text(size = 9, face = "bold"), 
    plot.title = element_text(size = 9, face = "bold")
  ) +
  stat_regline_equation(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~~")), 
                        formula = y ~ x, label.x = -60, label.y = 1500, output.type = "expression") + 
  stat_cor(method = "pearson", label.x = -60, label.y = 1000, p.accuracy = 0.001, r.accuracy = 0.001)

print(correlacionRendimiento_ggControl)




```

#### Cambio en altura máxima *vs* rendimiento (kg/ha)

::: justify
Los datos del cambio en la altura máxima *vs* rendimiento (kg/ha) mostraron una correlación significativa (*P* \< 0.001), con coeficiente de correlación de Pearson de -0.424, que se considera una correlación negativa moderada; y un coeficiente de determinación de 0.18, lo cual indica que aproximadamente 18% de la variabilidad del rendimiento es explicada por el cambio en la altura máxima.
:::

```{r CorrRendimientoCrecimiento, dev = "svg", dev.args=list(bg="transparent")}
#| label: fig-CorrVarRendimientoCrecimiento
#| fig-cap: "Correlación entre el "
#| echo: TRUE
#| warning: FALSE
#| eval: TRUE
#| lightbox:
#|   group: gallery
#| code-fold: true
#| code-summary: "Mostrar el código."

correlacionRendimiento = read_excel(path = "csv/correlacionRendimiento.xlsx")


correlacionRendimiento_ggCrecimiento <- ggplot(correlacionRendimiento, aes(cambioAltura, Rendimiento)) + 
  geom_point() +
  geom_smooth(method = 'lm', formula = y ~ x, se = TRUE, fill = "blue", color = "blue", fullrange = TRUE, na.rm = TRUE, alpha = 0.10, linetype = "dashed") +
  theme_classic() + 
  labs(title = expression("Rendimiento vs"), subtitle = "Crecimiento") + 
  xlab("Crecimiento") + 
  ylab("Rendimiento") + 
  theme(
    text = element_text(family = "montserrat"), 
    axis.text = element_text(size = 9), 
    axis.title = element_text(size = 9, face = "bold"), 
    plot.title = element_text(size = 9, face = "bold")
  ) +
  stat_regline_equation(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~~")), 
                        formula = y ~ x, label.x = 0, label.y = 1500, output.type = "expression") + 
  stat_cor(method = "pearson", label.x = 0, label.y = 1000, p.accuracy = 0.001, r.accuracy = 0.001)

print(correlacionRendimiento_ggCrecimiento)




```

#### Predicción del modelo de regresión vs rendimiento real

::: justify
La ecuación de regresión del modelo anterior nos permite calcular el rendimiento (*y*) conociendo el cambio en la altura máxima de las plantas de maleza (*x*). En la @tbl-cuadroCorrelacionRendimiento se muestra la comparación entre el rendimiento estimado por el modelo y el rendimiento promedio real de las 72 subparcelas y de cada uno de los tratamientos, en donde se puede apreciar un error promedio menor al 4%.
:::

```{r cuadroCorrelacionRendimiento, dev = "svg", dev.args=list(bg="transparent")}
#| label: tbl-cuadroCorrelacionRendimiento
#| tbl-cap: "Correlación entre el "
#| echo: TRUE
#| warning: FALSE
#| eval: TRUE
#| code-fold: true
#| code-summary: "Mostrar el código."
#| 


# Cargar paquetes necesarios
library(knitr)
library(kableExtra)

# Datos
tratamientos <- c("TODOS", "TS", "TR", "T4", "T3", "T2", "T1")
cambioAltura <- c(11, 42, 3, 7, -1, 4, 9)
rendimientoEstimado <- c(4739, 3986, 4931, 4825, 5018, 4897, 4777)
rendimientoReal <- c(4739, 3738, 5117, 5058, 4672, 4921, 4929)
diferenciaKg <- c(0, 248, -186, -233, 346, -24, -152)
diferenciaPorcentaje <- c(0.00, 6.64, -3.63, -4.61, 7.41, -0.48, -3.09)

# Crear el data frame
tabla <- data.frame(
  Tratamiento = tratamientos,
  cambioAltura = cambioAltura,
  rendimientoEstimado = rendimientoEstimado,
  rendimientoReal = rendimientoReal,
  diferenciaKg = diferenciaKg,
  diferenciaPorcentaje = diferenciaPorcentaje
)

# Calcular y agregar el error promedio
errorPromedio <- 3.69
errorPromedioFila <- data.frame(
  Tratamiento = "Error promedio (%)",
  cambioAltura = "",
  rendimientoEstimado = "",
  rendimientoReal = "",
  diferenciaKg = "",
  diferenciaPorcentaje = errorPromedio
)

# Agregar fila de error promedio
tabla <- rbind(tabla, errorPromedioFila)

# Crear y mostrar la tabla
kable(tabla, format = "html", caption = "Resultados del experimento") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive")) %>%
  row_spec(nrow(tabla), bold = TRUE, color = "white", background = "#0a231fff")



```

### Análisis de covarianza

::: justify
En este apartado, se llevó a cabo un análisis de covarianza (ANCOVA) para investigar si la combinación de dos variables, el cambio en el número de plantas y el cambio en la altura máxima, proporciona un modelo que explique mejor la variabilidad en el rendimiento (kg/ha) en comparación con el uso exclusivo de la altura máxima.

Los datos de la combinación de las dos variables vs rendimiento (kg/ha) mostraron una correlación significativa (P \< 0.001), con coeficiente de correlación de Pearson de -0.458, que se considera una correlación negativa moderada; y un coeficiente de determinación de 0.21, lo cual indica que aproximadamente 21% de la variabilidad del rendimiento es explicada por el cambio en la altura máxima.

Los datos de la combinación de las dos variables vs rendimiento (kg/ha) mostraron una correlación significativa (P \< 0.001), con coeficiente de correlación de Pearson de -0.458, que se considera una correlación negativa moderada; y un coeficiente de determinación de 0.21, lo cual indica que aproximadamente 21% de la variabilidad del rendimiento es explicada por la combinación del cambio en el número de plantas y el cambio en la altura máxima.
:::

```{r CovarianzaRendimiento, dev = "svg", dev.args=list(bg="transparent")}
#| label: fig-covarianzaCrecimiento
#| fig-cap: "Correlación entre el "
#| echo: TRUE
#| warning: FALSE
#| eval: TRUE
#| lightbox:
#|   group: gallery
#| code-fold: true
#| code-summary: "Mostrar el código."


# Ajustar el modelo de ANCOVA
modelo_ancova <- aov(Rendimiento ~ cambioNoPlantas + cambioAltura, data = correlacionRendimiento)

# Resumen del modelo
summary(modelo_ancova)

# Mostrar los coeficientes del modelo
coeficientes <- coef(modelo_ancova)

# Mostrar los coeficientes
print(coeficientes)

# Ajustar el modelo de ANCOVA
modelo_ancova2 <- lm(Rendimiento ~ cambioNoPlantas + cambioAltura, data = correlacionRendimiento)

# Resumen del modelo
resumen_modelo <- summary(modelo_ancova2)

# Extraer R^2
R2 <- resumen_modelo$r.squared

# Calcular R
R <- sqrt(R2)

# Mostrar los valores
cat("R^2:", R2, "\n")
cat("R:", R, "\n")


```

```{r graficoCovarianzaRendimiento, dev = "svg", dev.args=list(bg="transparent")}
#| label: fig-graficoCovarianzaCrecimiento
#| fig-cap: "Correlación entre el "
#| echo: TRUE
#| warning: FALSE
#| eval: TRUE
#| lightbox:
#|   group: gallery
#| code-fold: true
#| code-summary: "Mostrar el código."

# Instalar y cargar la librería plotly
#install.packages("htmlwidgets")
#install.packages("crosstalk")
#install.packages("promises")
#install.packages("later")
#install.packages("plotly")
library(htmlwidgets)
library(crosstalk)
library(promises)
library(later)
library(plotly)

# Ajustar el modelo de regresión lineal múltiple
modelo_ancova <- lm(Rendimiento ~ cambioNoPlantas + cambioAltura, data = correlacionRendimiento)

# Crear una cuadrícula de puntos para el plano de regresión
grid <- expand.grid(
  cambioNoPlantas = seq(min(correlacionRendimiento$cambioNoPlantas), max(correlacionRendimiento$cambioNoPlantas), length.out = 30),
  cambioAltura = seq(min(correlacionRendimiento$cambioAltura), max(correlacionRendimiento$cambioAltura), length.out = 30)
)

# Predecir los valores de Rendimiento para la cuadrícula
grid$Rendimiento <- predict(modelo_ancova, newdata = grid)

# Convertir las predicciones a una matriz con las dimensiones correctas
z_matrix <- matrix(grid$Rendimiento, nrow = length(unique(grid$cambioNoPlantas)), ncol = length(unique(grid$cambioAltura)))

# Crear el gráfico 3D con el plano de regresión
p4 <- plot_ly() %>%
  add_markers(data = correlacionRendimiento, x = ~cambioNoPlantas, y = ~cambioAltura, z = ~Rendimiento, 
              marker = list(size = 3), name = "Datos") %>%
  add_surface(x = ~unique(grid$cambioNoPlantas), y = ~unique(grid$cambioAltura), z = ~z_matrix,
              showscale = FALSE, opacity = 0.5, name = "Plano de Regresión") %>%
  layout(title = "Gráfico 3D con Plano de Regresión",
         scene = list(
           xaxis = list(title = "Cambio en No. Plantas"),
           yaxis = list(title = "Cambio en Altura"),
           zaxis = list(title = "Rendimiento"),
           camera = list(
             eye = list(x = 1.75, y = 1.75, z = 1.0) # Ajustar estos valores para cambiar el ángulo
           )
         ))

# Mostrar el gráfico 3D
p4




```

```{r graficoDosCovarianzaRendimiento, dev = "svg", dev.args=list(bg="transparent")}
#| label: fig-graficoDosCovarianzaCrecimiento
#| fig-cap: "Correlación entre el "
#| echo: TRUE
#| warning: FALSE
#| eval: TRUE
#| lightbox:
#|   group: gallery
#| code-fold: true
#| code-summary: "Mostrar el código."


library(scatterplot3d)
# Crear el gráfico 3D
scatter <- scatterplot3d(x = correlacionRendimiento$cambioNoPlantas, y = correlacionRendimiento$cambioAltura, z = correlacionRendimiento$Rendimiento,
                         main = "Gráfico 3D con Plano de Regresión",
                         xlab = "Cambio en No. Plantas", ylab = "Cambio en Altura (cm)", zlab = "Rendimiento (t/ha)",
                         pch = 16, highlight.3d = TRUE, angle = 55, 
                            label.tick.marks = T, 
                            z.ticklabs = c("0", "1", "2", "3", "4", "5", "6", "7"))

# Ajustar el modelo de regresión lineal múltiple
modelo_ancova <- lm(Rendimiento ~ cambioNoPlantas + cambioAltura, data = correlacionRendimiento)

# Crear una cuadrícula de puntos para el plano de regresión
grid <- expand.grid(
  cambioNoPlantas = seq(min(correlacionRendimiento$cambioNoPlantas), max(correlacionRendimiento$cambioNoPlantas), length.out = 30),
  cambioAltura = seq(min(correlacionRendimiento$cambioAltura), max(correlacionRendimiento$cambioAltura), length.out = 30)
)

# Predecir los valores de Rendimiento para la cuadrícula
grid$Rendimiento <- predict(modelo_ancova, newdata = grid)

# Verificar algunos valores predichos
#head(grid)

# Graficar el plano de regresión
scatter$plane3d(modelo_ancova, draw_polygon = TRUE, polygon_args = list(border = "blue"))


```

#### Predicción del modelo de regresión vs rendimiento real

::: justify
En la @tbl-cuadroCovarianzaRendimiento se muestra la comparación entre el rendimiento estimado por el modelo combinado y el rendimiento promedio real de las 72 subparcelas. El error promedio fue de 3.55%, sin mucha diferencia respecto al 3.69% obtenido al utilizar solamente la altura máxima. Esto sugiere que ambos modelos son muy similares en términos de su capacidad para explicar la variabilidad del rendimiento.
:::

```{r cuadroCovarianzaRendimiento, dev = "svg", dev.args=list(bg="transparent")}
#| label: tbl-cuadroCovarianzaRendimiento
#| tbl-cap: "Correlación entre el "
#| echo: TRUE
#| warning: FALSE
#| eval: TRUE
#| code-fold: true
#| code-summary: "Mostrar el código."


# Cargar paquetes necesarios
library(knitr)
library(kableExtra)

# Datos
tratamientos <- c("TODOS", "TS", "TR", "T4", "T3", "T2", "T1")
cambioNoPlantas <- c(-20, -8, -22, -17, -26, -23, -24)
cambioAltura <- c(11, 42, 3, 7, -1, 4, 9)
rendimientoEstimado <- c(4739, 3927, 4934, 4773, 5060, 4914, 4824)
rendimientoReal <- c(4739, 3738, 5117, 5058, 4672, 4921, 4929)
errorKg <- c(-0.28, 189.49, -182.74, -285.18, 388.05, -6.70, -105.02)
errorPorcentaje <- c(-0.01, 5.07, -3.57, -5.64, 8.31, -0.14, -2.13)

# Crear el data frame
tabla <- data.frame(
  Tratamiento = tratamientos,
  cambioNoPlantas = cambioNoPlantas,
  cambioAltura = cambioAltura,
  rendimientoEstimado = rendimientoEstimado,
  rendimientoReal = rendimientoReal,
  errorKg = errorKg,
  errorPorcentaje = errorPorcentaje
)

# Calcular y agregar el error promedio
errorPromedio <- mean(abs(tabla$errorPorcentaje))
errorPromedioFila <- data.frame(
  Tratamiento = "Error promedio (%)",
  cambioNoPlantas = "",
  cambioAltura = "",
  rendimientoEstimado = "",
  rendimientoReal = "",
  errorKg = "",
  errorPorcentaje = round(errorPromedio, 2)
)

# Agregar fila de error promedio
tabla <- rbind(tabla, errorPromedioFila)


# Crear y mostrar la tabla
kable(tabla, format = "html", caption = "Resultados del experimento") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive")) %>%
  row_spec(nrow(tabla), bold = TRUE, color = "white", background = "#0a231fff")


```

::: justify
Dado que los modelos son comparables en su capacidad predictiva, y considerando que es un principio estadístico común optar por el modelo más sencillo cuando los modelos son similares en desempeño, se podría preferir el modelo que utiliza únicamente la altura máxima. Este enfoque no solo simplifica la interpretación y la aplicación práctica, sino que también reduce la posibilidad de sobreajuste.
:::

## Comparaciones entre el método tradicional y el método basado en imagenes (BI).

::: justify
Los datos los conteos por el método tradicional y basado en imágenes de *C. murale* mostraron una correlación significativa (P \< 0.001), con coeficiente de correlación de Pearson de 0.703, que se considera una correlación positiva fuerte; y un coeficiente de determinación de 0.86 (@fig-Cmurale).
:::

```{r Cmurale, dev = "svg", dev.args=list(bg="transparent")}
#| label: fig-Cmurale
#| fig-cap: "Correlación entre el número de plantas de C. murale contadas mediante el método tradicional y el método basado en imágenes."
#| echo: TRUE
#| warning: FALSE
#| eval: TRUE
#| lightbox:
#|   group: gallery
#| code-fold: true
#| code-summary: "Mostrar el código."

NOPLANTAS0DDA = read_excel(path = "csv/previoManualvsBI.xlsx")


Cmurale.NOPLANTAS.0dda <- ggplot(NOPLANTAS0DDA, aes(CmuraleBI, CmuraleManual)) + 
  geom_point() +
  geom_smooth(method = 'lm', formula = y ~ x - 1, se = TRUE, fill = "blue", color = "blue", fullrange = TRUE, na.rm = TRUE, alpha = 0.10, linetype = "dashed") + 
  geom_abline(intercept = 0, slope = 1, color = "black") +
  theme_classic() + 
  labs(title = expression("Número de plantas ("*italic("C. murale")*")"), subtitle = "Conteo previo") + 
  xlab("Método BI") + 
  ylab("Método tradicional") + 
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 5, 10, 15, 20, 25, 30, 35)) + 
  scale_y_continuous(expand = c(0, 0), breaks = c(0, 5, 10, 15, 20, 25, 30, 35)) +
  coord_cartesian(xlim = c(0, 35), ylim = c(0, 35)) + 
  theme(
    text = element_text(family = "montserrat"), 
    axis.text = element_text(size = 9), 
    axis.title = element_text(size = 9, face = "bold"), 
    plot.title = element_text(size = 9, face = "bold"),
    aspect.ratio = 1
  ) +
  stat_regline_equation(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~~")), 
                        formula = y ~ x - 1, label.x = 3, label.y = 33, output.type = "expression") + 
  stat_cor(method = "pearson", label.x = 3, label.y = 30, p.accuracy = 0.001, r.accuracy = 0.001) 



print(Cmurale.NOPLANTAS.0dda)




```

::: justify
Los datos los conteos por el método tradicional y basado en imágenes de *A. palmeri* mostraron una correlación significativa (P \< 0.001), con coeficiente de correlación de Pearson de 0.833, que se considera una correlación positiva fuerte; y un coeficiente de determinación de 0.88 (@fig-Apalmeri).
:::

```{r Apalmeri, dev = "svg", dev.args=list(bg="transparent")}
#| label: fig-Apalmeri
#| fig-cap: "Correlación entre el número de plantas de A. palmeri contadas mediante el método tradicional y el método basado en imágenes."
#| echo: TRUE
#| warning: FALSE
#| eval: TRUE
#| lightbox:
#|   group: gallery
#| code-fold: true
#| code-summary: "Mostrar el código."


Apalmeri.NOPLANTAS.0dda <- ggplot(NOPLANTAS0DDA, aes(ApalmeriBI, ApalmeriManual)) + 
  geom_point() +
  geom_smooth(method = 'lm', formula = y ~ x - 1, se = TRUE, fill = "blue", color = "blue", fullrange = TRUE, na.rm = TRUE, alpha = 0.10, linetype = "dashed") + 
  geom_abline(intercept = 0, slope = 1, color = "black") +
  theme_classic() + 
  labs(title = expression("Número de plantas ("*italic("A. palmeri")*")"), subtitle = "Conteo previo") + 
  xlab("Método BI") + 
  ylab("Método tradicional") + 
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 10, 20, 30, 40, 50, 60, 65)) + 
  scale_y_continuous(expand = c(0, 0), breaks = c(0, 10, 20, 30, 40, 50, 60, 65)) +
  coord_cartesian(xlim = c(0, 65), ylim = c(0, 65)) + 
  theme(
    text = element_text(family = "montserrat"), 
    axis.text = element_text(size = 9), 
    axis.title = element_text(size = 9, face = "bold"), 
    plot.title = element_text(size = 9, face = "bold"),
    aspect.ratio = 1
  ) +
  stat_regline_equation(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~~")), 
                        formula = y ~ x - 1, label.x = 3, label.y = 58, output.type = "expression") + stat_cor(method = "pearson", label.x = 3, label.y = 55, p.accuracy = 0.001, r.accuracy = 0.001)

Apalmeri.NOPLANTAS.0dda



```

::: justify
Los datos los conteos por el método tradicional y basado en imágenes de *P. acutifolia* mostraron una correlación significativa (P \< 0.001), con coeficiente de correlación de Pearson de 0.774, que se considera una correlación positiva fuerte; y un coeficiente de determinación de 0.85 (@fig-Pacutifolia).
:::

```{r Pacutifolia, dev = "svg", dev.args=list(bg="transparent")}
#| label: fig-Pacutifolia
#| fig-cap: "Correlación entre el número de plantas de P. acutifolia contadas mediante el método tradicional y el método basado en imágenes."
#| echo: TRUE
#| warning: FALSE
#| eval: TRUE
#| lightbox:
#|   group: gallery
#| code-fold: true
#| code-summary: "Mostrar el código."

Pacutifolia.NOPLANTAS.0dda <- ggplot(NOPLANTAS0DDA, aes(PacutifoliaBI, PacutifoliaManual)) + 
  geom_point() +
  geom_smooth(method = 'lm', formula = y ~ x - 1, se = TRUE, fill = "blue", color = "blue", fullrange = TRUE, na.rm = TRUE, alpha = 0.10, linetype = "dashed") + 
  geom_abline(intercept = 0, slope = 1, color = "black") +
  theme_classic() + 
  labs(title = expression("Número de plantas ("*italic("P. acutifolia")*")"), subtitle = "Conteo previo") + 
  xlab("Método BI") + 
  ylab("Método tradicional") + 
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 5, 10, 15)) + 
  scale_y_continuous(expand = c(0, 0), breaks = c(0, 5, 10, 15)) +
  coord_cartesian(xlim = c(0, 15), ylim = c(0, 15)) + 
  theme(
    text = element_text(family = "montserrat"), 
    axis.text = element_text(size = 9), 
    axis.title = element_text(size = 9, face = "bold"), 
    plot.title = element_text(size = 9, face = "bold"),
    aspect.ratio = 1
  ) +
  stat_regline_equation(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~~")), 
                        formula = y ~ x - 1, label.x = 1, label.y = 13, output.type = "expression") + stat_cor(method = "pearson", label.x = 1, label.y = 12, p.accuracy = 0.001, r.accuracy = 0.001)

Pacutifolia.NOPLANTAS.0dda



```

::: justify
Los datos los conteos por el método tradicional y basado en imágenes del complejo de maleza mostraron una correlación significativa (P \< 0.001), con coeficiente de correlación de Pearson de 0.823, que se considera una correlación positiva fuerte; y un coeficiente de determinación de 0.94 (@fig-Complejo).
:::

```{r Complejo, dev = "svg", dev.args=list(bg="transparent")}
#| label: fig-Complejo
#| fig-cap: "Correlación entre el número de plantas del complejo de maleza contadas mediante el método tradicional y el método basado en imágenes."
#| echo: TRUE
#| warning: FALSE
#| eval: TRUE
#| lightbox:
#|   group: gallery
#| code-fold: true
#| code-summary: "Mostrar el código."

Complejo.NOPLANTAS.0dda <- ggplot(NOPLANTAS0DDA, aes(ComplejoBI, ComplejoManual)) + 
  geom_point() +
  geom_smooth(method = 'lm', formula = y ~ x - 1, se = TRUE, fill = "blue", color = "blue", fullrange = TRUE, na.rm = TRUE, alpha = 0.10, linetype = "dashed") + 
  geom_abline(intercept = 0, slope = 1, color = "black") +
  theme_classic() + 
  labs(title = expression("Número de plantas ("*italic("Complejo de maleza")*")"), subtitle = "Conteo previo") + 
  xlab("Método BI") + 
  ylab("Método tradicional") + 
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 85)) + 
  scale_y_continuous(expand = c(0, 0), breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 85)) +
  coord_cartesian(xlim = c(0, 85), ylim = c(0, 85)) + 
  theme(
    text = element_text(family = "montserrat"), 
    axis.text = element_text(size = 9), 
    axis.title = element_text(size = 9, face = "bold"), 
    plot.title = element_text(size = 9, face = "bold"),
    aspect.ratio = 1
  ) +
  stat_regline_equation(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~~")), 
                        formula = y ~ x - 1, label.x = 3, label.y = 78, output.type = "expression") + stat_cor(method = "pearson", label.x = 3, label.y = 73, p.accuracy = 0.001, r.accuracy = 0.001)

Complejo.NOPLANTAS.0dda

```

## Control de maleza

### Conteo previo (0 dda)

#### Descriptivo

```{r noPlantas0ddaDescriptivo}



# 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(noPlantas0dda, na.rm = TRUE),
            mediana = median(noPlantas0dda, na.rm = TRUE),
            sd = sd(noPlantas0dda, na.rm = TRUE),
            varianza = var(noPlantas0dda, na.rm = TRUE),
            min = min(noPlantas0dda, na.rm = TRUE),
            max = max(noPlantas0dda, na.rm = TRUE))


print(summary_stats)


```

#### Datos sin transformar

##### Análisis de Varianza

```{r noPlantas0ddaANOVA}


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

# Mostrar los resultados del ANOVA
print(anova_results)


```

##### Análisis de normalidad de los residuos utilizando Shapiro-Wilk

```{r shapironoPlantas0dda}


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



```

##### Comparaciones múltiples

```{r noPlantas0ddaTUKEY}

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

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


```

#### Transformación de Box-Cox

```{r noPlantas0ddaBoxcox}
#| echo: TRUE
#| warning: FALSE
#| eval: TRUE
#| lightbox:
#|   group: gallery2


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

# 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)

# Aplicar la transformación de Box-Cox
boxcox_results <- boxcox(noPlantas0dda ~ 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))

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

# Ajustar el modelo con los datos transformados
model_transformed <- aov(transformed_noPlantas0dda ~ 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_noPlantas0dda ~ Tratamiento, data = data)

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

# 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)

```

### Conteo final (28 dda)

#### Descriptivo

```{r noPlantas28ddaDescriptivo}



# 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(noPlantas28dda, na.rm = TRUE),
            mediana = median(noPlantas28dda, na.rm = TRUE),
            sd = sd(noPlantas28dda, na.rm = TRUE),
            varianza = var(noPlantas28dda, na.rm = TRUE),
            min = min(noPlantas28dda, na.rm = TRUE),
            max = max(noPlantas28dda, na.rm = TRUE))


print(summary_stats)


```

#### Datos sin transformar

##### Análisis de Varianza

```{r noPlantas28ddaANOVA}


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

# Mostrar los resultados del ANOVA
print(anova_results)


```

##### Análisis de normalidad de los residuos utilizando Shapiro-Wilk

```{r shapironoPlantas28dda}


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



```

##### Comparaciones múltiples

```{r noPlantas28ddaTUKEY}

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

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


```

#### Transformación de Box-Cox

```{r noPlantas28ddaBoxcox}
#| echo: TRUE
#| warning: FALSE
#| eval: TRUE
#| lightbox:
#|   group: gallery2


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

# 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$noPlantas28dda <- data$noPlantas28dda + 1

# Aplicar la transformación de Box-Cox
boxcox_results <- boxcox(noPlantas28dda ~ 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))

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

# Ajustar el modelo con los datos transformados
model_transformed <- aov(transformed_noPlantas28dda ~ 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_noPlantas28dda ~ Tratamiento, data = data)

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

# 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)

```

#### Análisis no paramétrico

```{r noPlantas28ddaKruskal}

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


```

##### Test de Dunn

```{r noPlantas28ddaDunn}


library(FSA)
dunn_test <- dunnTest(noPlantas28dda ~ Tratamiento , data = data, method = "bonferroni")
print(dunn_test)


```

### Cambio en número de plantas (28 dda)

#### Descriptivo

```{r cambioNoPlantasDescriptivo}



# 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(cambioNoPlantas, na.rm = TRUE),
            mediana = median(cambioNoPlantas, na.rm = TRUE),
            sd = sd(cambioNoPlantas, na.rm = TRUE),
            varianza = var(cambioNoPlantas, na.rm = TRUE),
            min = min(cambioNoPlantas, na.rm = TRUE),
            max = max(cambioNoPlantas, na.rm = TRUE))


print(summary_stats)


```

#### Datos sin transformar

##### Análisis de Varianza

```{r cambioNoPlantasANOVA}


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

# Mostrar los resultados del ANOVA
print(anova_results)


```

##### Análisis de normalidad de los residuos utilizando Shapiro-Wilk

```{r shapirocambioNoPlantas}


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



```

##### Comparaciones múltiples

```{r cambioNoPlantasTUKEY}

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

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


```

### Altura máxima (0 dda)

#### Descriptivo

```{r hMax0ddaDescriptivo}



# 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(hMax0dda, na.rm = TRUE),
            mediana = median(hMax0dda, na.rm = TRUE),
            sd = sd(hMax0dda, na.rm = TRUE),
            varianza = var(hMax0dda, na.rm = TRUE),
            min = min(hMax0dda, na.rm = TRUE),
            max = max(hMax0dda, na.rm = TRUE))


print(summary_stats)


```

#### Datos sin transformar

##### Análisis de Varianza

```{r hMax0ddaANOVA}


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

# Mostrar los resultados del ANOVA
print(anova_results)


```

##### Análisis de normalidad de los residuos utilizando Shapiro-Wilk

```{r shapirohMax0dda}


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



```

##### Comparaciones múltiples

```{r hMax0ddaTUKEY}

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

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


```

### Altura máxima (28 dda)

#### Descriptivo

```{r hMax28ddaDescriptivo}



# 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(hMax28dda, na.rm = TRUE),
            mediana = median(hMax28dda, na.rm = TRUE),
            sd = sd(hMax28dda, na.rm = TRUE),
            varianza = var(hMax28dda, na.rm = TRUE),
            min = min(hMax28dda, na.rm = TRUE),
            max = max(hMax28dda, na.rm = TRUE))


print(summary_stats)


```

#### Datos sin transformar

##### Análisis de Varianza

```{r hMax28ddaANOVA}


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

# Mostrar los resultados del ANOVA
print(anova_results)


```

##### Análisis de normalidad de los residuos utilizando Shapiro-Wilk

```{r shapirohMax28dda}


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



```

#### Análisis de homogeneidad de varianzas utilizando la prueba de Levene

```{r levenehMax28dda}


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


```

#### Análisis de independencia de los residuos

```{r residuoshMax28dda}


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


```

##### Comparaciones múltiples

```{r hMax28ddaTUKEY}

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

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


```

### Cambio en Altura máxima (28 dda)

#### Descriptivo

```{r cambioAlturaDescriptivo}



# 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(cambioAltura, na.rm = TRUE),
            mediana = median(cambioAltura, na.rm = TRUE),
            sd = sd(cambioAltura, na.rm = TRUE),
            varianza = var(cambioAltura, na.rm = TRUE),
            min = min(cambioAltura, na.rm = TRUE),
            max = max(cambioAltura, na.rm = TRUE))


print(summary_stats)


```

#### Datos sin transformar

##### Análisis de Varianza

```{r cambioAlturaANOVA}


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

# Mostrar los resultados del ANOVA
print(anova_results)


```

##### Análisis de normalidad de los residuos utilizando Shapiro-Wilk

```{r shapirocambioAltura}


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



```

#### Análisis de homogeneidad de varianzas utilizando la prueba de Levene

```{r levenecambioAltura}


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


```

#### Análisis de independencia de los residuos

```{r residuoscambioAltura}


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


```

##### Comparaciones múltiples

```{r cambioAlturaTUKEY}

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

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


```

## Rendimiento

### Descriptivo

```{r RendimientoDescriptivo}



# 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))


print(summary_stats)


```

### Datos sin transformar

#### Análisis de Varianza

```{r RendimientoANOVA}


# 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)


```

#### Análisis de normalidad de los residuos utilizando Shapiro-Wilk

```{r shapiroRendimiento}


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



```

#### Análisis de homogeneidad de varianzas utilizando la prueba de Levene

```{r leveneRendimiento}


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


```

#### Análisis de independencia de los residuos

```{r residuosRendimiento}


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


```

#### Comparaciones múltiples

```{r RendimientoTUKEY}

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

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


```

### Transformación de Box-Cox

```{r Boxcox}

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

# 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))

# 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)
print(shapiro_test_transformed)
print(levene_test_transformed)

# 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)

```

### Análisis no paramétrico

```{r RendimientoKruskal}

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


```

#### Test de Dunn

```{r RendimientoDunn}


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


```

### Boxplot de rendimientos por tratamiento

```{r RendimientoBoxplot}
#| label: fig-RendimientoBoxplot
#| fig-cap: "RendimientoBoxplot."
#| echo: TRUE
#| warning: FALSE
#| eval: TRUE
#| lightbox:
#|   group: gallery
#| code-fold: true
#| code-summary: "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)






```

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

```{r RendimientoBarras}
#| label: fig-RendimientoBarras
#| fig-cap: "RendimientoBarras."
#| echo: TRUE
#| warning: FALSE
#| eval: TRUE
#| lightbox:
#|   group: gallery
#| code-fold: true
#| code-summary: "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)



```

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

```{r RendimientoBarrasPorcentaje}
#| label: fig-RendimientoBarrasPorcentaje
#| fig-cap: "RendimientoBarrasPorcentaje."
#| echo: TRUE
#| warning: FALSE
#| eval: TRUE
#| lightbox:
#|   group: gallery
#| code-fold: true
#| code-summary: "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)



```

## Comparación respecto a la producción histórica

![Pérdidas de producción debido a diversos factores](images/factoresProduccion2.png){#fig-perdidas group="gallery"}

------------------------------------------------------------------------

::: columns
::: {.column width="49%"}

------------------------------------------------------------------------
