Visualización de datos

Evaluaciones agropecuarias

Bibliotecas

library(readxl)
library(tidyverse)
library(janitor)
library(ggthemes)
library(plotly)

Datos EVAs

datos <- read_excel("datos/Base agrícola 2019 - 2023.xlsx",
                    skip = 6) |> 
  clean_names()

datos |> head()

Datos embalses

df_embalses <- 
  read_csv("datos/PorcVoluUtilDiar.csv")

df_embalses |> head()

Cantidades y proporciones

Registros por año

• Primero calculamos el total de registros por año: primera capa para ggplot2 –> Datos.

conteo_ano <-
  datos |> 
  count(ano, name = "total")

conteo_ano

• Agreguemos las dos capas adicionales para obtener el gráfico:

conteo_ano |> 
  ggplot(mapping = aes(x = ano, y = total)) +
  geom_col()

• También es posibles declarar el gráfico de la siguiente manera:

ggplot(data = conteo_ano,
       mapping = aes(x = ano, y = total)) +
  geom_col()

• Ahora vamos a agregar título, subtítulo, etiquetas a los ejes “x” y “y”, editar colores de las barras y cambiar el tema del gráfico:

conteo_ano |> 
  ggplot(mapping = aes(x = ano, y = total)) +
  geom_col(color = "#657C6A", fill = "#BB3E00", width = 0.75, alpha = 0.5) +
  labs(
    x = "Año",
    y = "Registros (n)",
    title = "Número de evaluaciones agropecuarias por año",
    subtitle = "Años 2019 a 2024 en Colombia",
    caption = "Datos obtenidos desde la UPRA"
  ) +
  theme_minimal()

• Ahora agregaremos texto encima de las barras: geom_text()

conteo_ano |> 
  ggplot(mapping = aes(x = ano, y = total)) +
  geom_col(color = "#657C6A", fill = "#BB3E00", width = 0.75, alpha = 0.5) +
  geom_text(aes(label = total), color = "#BB3E00", size = 5) +
  labs(
    x = "Año",
    y = "Registros (n)",
    title = "Número de evaluaciones agropecuarias por año",
    subtitle = "Años 2019 a 2024 en Colombia",
    caption = "Datos obtenidos desde la UPRA"
  ) +
  theme_minimal() 

• Otra opción para agregar texto es: geom_label()

conteo_ano |> 
  ggplot(mapping = aes(x = ano, y = total)) +
  geom_col(color = "#657C6A", fill = "#BB3E00", width = 0.75, alpha = 0.5) +
  geom_label(aes(label = total), color = "white", fill = "#BB3E00", size = 3.5) +
  labs(
    x = "Año",
    y = "Registros (n)",
    title = "Número de evaluaciones agropecuarias por año",
    subtitle = "Años 2019 a 2024 en Colombia",
    caption = "Datos obtenidos desde la UPRA"
  ) +
  theme_minimal() 

Registros por departamento

• Primera opción:

datos |> 
  count(departamento, name = "total") |> 
  ggplot(aes(x = departamento, y = total)) +
  geom_col()

• Una opción mejorada podría ser invertir o rotar los ejes: coord_flip()

datos |> 
  count(departamento, name = "total") |> 
  ggplot(aes(x = reorder(departamento, total), y = total)) +
  geom_col() +
  coord_flip()

• Otra opción: rotando las etiquetas del eje “x”. En los documentos de Quarto podemos controlar el ancho y el alto de la ventana en donde se incorporan los gráficos:
  • El ancho de las figuras máximo podrá ser configurado en 9
  • La altura de las figuras podrá ser configurada en cualquier número positivo

datos |> 
  count(departamento, name = "total") |> 
  ggplot(aes(x = reorder(departamento, total), y = total)) +
  geom_col() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

Diagrama de barras anidado (2 variables)

• ¿Para el cultivo de café en cada departamento cuál es la proporción de registros que superan la media nacional del rendimiento?

df_resumen_cafe <-
  datos |>
  filter(cultivo == "Café") |>
  mutate(mayor_promedio = if_else(
    condition = rendimiento_t_ha > mean(rendimiento_t_ha),
    true = "Sí",
    false = "No"
  )) |> 
  count(departamento, mayor_promedio, name = "total")

df_resumen_cafe

• Queremos graficar el departamento y la variable “mayor_promedio” —> Primera opción:

df_resumen_cafe |> 
  ggplot(aes(x = departamento, y = total, color = mayor_promedio, fill = mayor_promedio)) +
  geom_col() +
  coord_flip()

• Segunda opción: poner una barra al lado de la otra en cada departamento:

df_resumen_cafe |> 
  ggplot(aes(x = departamento, y = total, color = mayor_promedio, fill = mayor_promedio)) +
  geom_col(position = "dodge")  +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

• Tercera opción: relativar el total —> Usar proporciones

df_resumen_cafe |> 
  ggplot(aes(x = departamento, y = total, color = mayor_promedio, fill = mayor_promedio)) +
  geom_col(position = "fill")  +
  coord_flip() +
  labs(
    x = "",
    y = "Proporción",
    color = "¿Mayor al promedio nacional?",
    fill = "¿Mayor al promedio nacional?"
  ) +
  theme_minimal() +
  theme(legend.position = "top") 

• El mismo gráfico anterior editando color (scale_color_manual) y fill (scale_fill_manual):

df_resumen_cafe |> 
  ggplot(aes(x = departamento, y = total, color = mayor_promedio, fill = mayor_promedio)) +
  geom_col(position = "fill", alpha = 0.85)  +
  coord_flip() +
  labs(
    x = "",
    y = "Proporción",
    color = "¿Mayor al promedio nacional?",
    fill = "¿Mayor al promedio nacional?",
    title = "Cultivo de café"
  ) +
  scale_color_manual(values = c("#213448", "#94B4C1")) +
  scale_fill_manual(values = c("#213448", "#94B4C1")) +
  theme_minimal() +
  theme(legend.position = "top") 

Diagrama de barras anidado (3 variables)

• ¿Por año para el cultivo de maíz en cada departamento cuál es la proporción de registros que superan la media nacional del rendimiento?

df_resumen_maiz <-
  datos |>
  filter(cultivo == "Maíz") |>
  mutate(mayor_promedio = if_else(
    condition = rendimiento_t_ha > mean(rendimiento_t_ha),
    true = "Sí",
    false = "No"
  )) |> 
  count(ano, departamento, mayor_promedio, name = "total")

df_resumen_maiz

• Graficamos con facetas (facet_wrap):

df_resumen_maiz |> 
  ggplot(aes(x = departamento, y = total, color = mayor_promedio, fill = mayor_promedio)) +
  facet_wrap(facets = ~ ano) +
  geom_col(position = "fill", alpha = 0.85)  +
  coord_flip() +
  labs(
    x = "",
    y = "Proporción",
    color = "¿Mayor al promedio nacional?",
    fill = "¿Mayor al promedio nacional?",
    title = "Cultivo de maíz"
  ) +
  scale_color_manual(values = c("#213448", "#94B4C1")) +
  scale_fill_manual(values = c("#213448", "#94B4C1")) +
  theme_minimal() +
  theme(legend.position = "top") 

Diagrama de barras anidado (4 variables)

• ¿Qué proporción de evaluaciones agropecuarias por año, departamento y desagregación superan la media nacional de rendimieento?
• Nota: para que la comparación sea más justa calcularemos la media nacional para cada desagregación y luego sí generamos la variable “mayor_promedio”.

df_resumen_tomate <-
  datos |>
  filter(cultivo == "Tomate") |>
  group_by(desagregacion_cultivo) |> 
  mutate(media_nacional = mean(rendimiento_t_ha)) |> 
  ungroup() |>
  mutate(mayor_promedio = if_else(
    condition = rendimiento_t_ha > media_nacional,
    true = "Sí",
    false = "No"
  )) |>
  count(ano,
        departamento,
        desagregacion_cultivo,
        mayor_promedio,
        name = "total")

df_resumen_tomate

• Ahora graficamos con facetas: facet_wrap() - Opción 1

df_resumen_tomate |> 
  ggplot(aes(x = total, y = departamento, color = mayor_promedio, fill = mayor_promedio)) +
  facet_wrap(~ ano ~ desagregacion_cultivo, ncol = 2) +
  geom_col(position = "fill", alpha = 0.85)  +
  labs(
    x = "",
    y = "Proporción",
    color = "¿Mayor al promedio nacional?",
    fill = "¿Mayor al promedio nacional?",
    title = "Cultivo de tomate",
    subtitle = "Comparativo de desagregaciones"
  ) +
  scale_color_manual(values = c("#FE7743", "#273F4F")) +
  scale_fill_manual(values = c("#FE7743", "#273F4F")) +
  theme_minimal() +
  theme(legend.position = "top") 

• Ahora graficamos con facetas: facet_wrap() - Opción 2

df_resumen_tomate |> 
  ggplot(aes(x = total, y = departamento, color = mayor_promedio, fill = mayor_promedio)) +
  facet_wrap(~ desagregacion_cultivo  ~ ano, ncol = 5) +
  geom_col(position = "fill", alpha = 0.85)  +
  labs(
    x = "",
    y = "Proporción",
    color = "¿Mayor al promedio nacional?",
    fill = "¿Mayor al promedio nacional?",
    title = "Cultivo de tomate",
    subtitle = "Comparativo de desagregaciones"
  ) +
  scale_color_manual(values = c("#FE7743", "#273F4F")) +
  scale_fill_manual(values = c("#FE7743", "#273F4F")) +
  theme_minimal() +
  theme(legend.position = "top") 

• Ahora graficamos con facetas: facet_grid() - Opción 3

df_resumen_tomate |> 
  ggplot(aes(x = total, y = departamento, color = mayor_promedio, fill = mayor_promedio)) +
  facet_grid(~ ano ~ desagregacion_cultivo) +
  geom_col(position = "fill", alpha = 0.85)  +
  labs(
    x = "",
    y = "Proporción",
    color = "¿Mayor al promedio nacional?",
    fill = "¿Mayor al promedio nacional?",
    title = "Cultivo de tomate",
    subtitle = "Comparativo de desagregaciones"
  ) +
  scale_color_manual(values = c("#FE7743", "#273F4F")) +
  scale_fill_manual(values = c("#FE7743", "#273F4F")) +
  theme_minimal() +
  theme(legend.position = "top") 

• Ahora graficamos con facetas: facet_grid() - Opción 4

df_resumen_tomate |> 
  ggplot(aes(x = total, y = departamento, color = mayor_promedio, fill = mayor_promedio)) +
  facet_grid(~ desagregacion_cultivo ~ ano) +
  geom_col(position = "fill", alpha = 0.85)  +
  labs(
    x = "",
    y = "Proporción",
    color = "¿Mayor al promedio nacional?",
    fill = "¿Mayor al promedio nacional?",
    title = "Cultivo de tomate",
    subtitle = "Comparativo de desagregaciones"
  ) +
  scale_color_manual(values = c("#FE7743", "#273F4F")) +
  scale_fill_manual(values = c("#FE7743", "#273F4F")) +
  theme_minimal() +
  theme(legend.position = "top") 

Interactividad con plotly

• La función que utilizamos para conseguir la interactividad es “ggplotly”.
• La primera opción es incorporar el gráfico dentro de esta función:

ggplotly(
  df_resumen_cafe |>
    ggplot(
      aes(
        x = departamento,
        y = total,
        color = mayor_promedio,
        fill = mayor_promedio
      )
    ) +
    geom_col(position = "fill", alpha = 0.85)  +
    coord_flip() +
    labs(
      x = "",
      y = "Proporción",
      color = "¿Mayor al promedio nacional?",
      fill = "¿Mayor al promedio nacional?",
      title = "Cultivo de café"
    ) +
    scale_color_manual(values = c("#213448", "#94B4C1")) +
    scale_fill_manual(values = c("#213448", "#94B4C1")) +
    theme_minimal() +
    theme(legend.position = "top")
)

• La segunda opción es primero crear el gráfico, asignarlo (guardarlo) y luego usar la función “ggplotly”:

# Creación y asignación del gráfico
grafico1 <-
  df_resumen_cafe |>
    ggplot(
      aes(
        x = departamento,
        y = total,
        color = mayor_promedio,
        fill = mayor_promedio
      )
    ) +
    geom_col(position = "fill", alpha = 0.85)  +
    coord_flip() +
    labs(
      x = "",
      y = "Proporción",
      color = "¿Mayor al promedio nacional?",
      fill = "¿Mayor al promedio nacional?",
      title = "Cultivo de café"
    ) +
    scale_color_manual(values = c("#213448", "#94B4C1")) +
    scale_fill_manual(values = c("#213448", "#94B4C1")) +
    theme_minimal() +
    theme(legend.position = "top")

# Uso de ggplotly
ggplotly(grafico1)

Cantidad - Promedio

• ¿Cuál es el top 10 de municipios en promedio de rendimiento para el cultivo de papa criolla?

df_papa_criolla <-
  datos |> 
  filter(desagregacion_cultivo == "Papa criolla")

df_papa_criolla |> 
  group_by(municipio) |> 
  reframe(promedio = mean(rendimiento_t_ha, na.rm = TRUE)) |> 
  arrange(desc(promedio)) |> 
  slice(1:10) |> 
  ggplot(aes(x = reorder(municipio, promedio), y = promedio)) +
  geom_col()

• ¿Cuál es el promedio de rendimiento por departamento y año? Este gráfico es errado!

df_papa_criolla |> 
  group_by(departamento, ano) |> 
  reframe(promedio = mean(rendimiento_t_ha, na.rm = TRUE)) |> 
  mutate(ano = as.factor(ano)) |>
  ggplot(aes(x = departamento, y = promedio, color = ano, fill = ano)) +
  geom_col()

• Dado que estamos graficando el promedio no es posible totalizarlo, es decir, necesitamos separar las barras de cada año en cada departamento:

df_papa_criolla |> 
  group_by(departamento, ano) |> 
  reframe(promedio = mean(rendimiento_t_ha, na.rm = TRUE)) |> 
  mutate(ano = as.factor(ano)) |>
  ggplot(aes(x = departamento, y = promedio, color = ano, fill = ano)) +
  geom_col(position = "dodge")

• También podría intercambiar “color-fill” con “x”: además, editamos los nueve colores.

paleta_9 <-
  c(
    "#a6cee3",
    "#1f78b4",
    "#b2df8a",
    "#33a02c",
    "#fb9a99",
    "#e31a1c",
    "#fdbf6f",
    "#ff7f00",
    "#cab2d6"
  )

df_papa_criolla |> 
  group_by(departamento, ano) |> 
  reframe(promedio = mean(rendimiento_t_ha, na.rm = TRUE)) |> 
  mutate(ano = as.factor(ano)) |>
  ggplot(aes(x = ano, y = promedio, color = departamento, fill = departamento)) +
  geom_col(position = "dodge") +
  scale_color_manual(values = paleta_9) +
  scale_fill_manual(values = paleta_9)

Cantidad - Desviación estándar

df_papa_criolla |> 
  group_by(departamento, ano) |> 
  reframe(desviacion = sd(rendimiento_t_ha, na.rm = TRUE)) |> 
  mutate(ano = as.factor(ano)) |>
  ggplot(aes(x = ano, y = desviacion, color = departamento, fill = departamento)) +
  geom_col(position = "dodge") +
  scale_color_manual(values = paleta_9) +
  scale_fill_manual(values = paleta_9)

Cantidad - Sumatoria

• ¿En cuál semestre hay mayor participación del volumen de papa criolla producida entre el 2019 y 2023?

df_papa_criolla |> 
  group_by(ano, periodo) |> 
  reframe(total_pdn = sum(produccion_t, na.rm = TRUE)) |> 
  ggplot(aes(x = ano, y = total_pdn, color = periodo, fill = periodo)) +
  geom_col()

• En el gráfico anterior es necesario unificar los semestres antes de graficar:

df_papa_criolla |> 
  group_by(ano, periodo) |> 
  reframe(total_pdn = sum(produccion_t, na.rm = TRUE)) |> 
  mutate(semestre = if_else(
    condition = str_detect(periodo, "A"),
    true = "A",
    false = "B"
  )) |> 
  ggplot(aes(x = ano, y = total_pdn, color = semestre, fill = semestre)) +
  geom_col(position = "fill")

• El patrón de comportamiento anterior se replica en los departamentos?

df_papa_criolla |> 
  group_by(ano, periodo, departamento) |> 
  reframe(total_pdn = sum(produccion_t, na.rm = TRUE)) |> 
  mutate(semestre = if_else(
    condition = str_detect(periodo, "A"),
    true = "A",
    false = "B"
  )) |> 
  ggplot(aes(x = ano, y = total_pdn, color = semestre, fill = semestre)) +
  facet_wrap(~departamento) +
  geom_col(position = "fill")

Uso de esquisse

• Primero, instalar la biblioteca “esquisse”
• Segundo, ejecutar la siguiente función: “esquisser(viewer =”browser”)”

library(esquisse)

esquisser(viewer = "browser")

library(dplyr)
library(ggplot2)

df_papa_criolla %>%
 filter(!(departamento %in% "Caldas")) %>%
 ggplot() +
 aes(x = departamento) +
 geom_bar(fill = "#AB18CB") +
 labs(x = "Departamento", 
 y = "Total (n)", title = "Conteo por departamento", subtitle = "Papa criolla (2019-2023)", caption = "Colombia") +
 coord_flip() +
 theme_minimal()

Distribuciones

• Haremos los ejemplos con el cultivo de tomate:

df_tomate <-
  datos |> 
  filter(cultivo == "Tomate") |> 
  filter(rendimiento_t_ha > 0)

Individuales

Histogramas

df_tomate |> 
  ggplot(aes(x = rendimiento_t_ha)) +
  geom_histogram()

• Mejorar el gráfico anterior: Regla de Sturges —> https://es.wikipedia.org/wiki/Regla_de_Sturges

\[1 + log_2(n) = 1 + log_2(5007) = 13.28\]

df_tomate |>
  ggplot(aes(x = rendimiento_t_ha)) +
  geom_histogram(
    color = "#309898",
    fill = "#F4631E",
    alpha = 0.75,
    bins = 13
  )

• Aplicamos logaritmos al eje x:

df_tomate |>
  ggplot(aes(x = rendimiento_t_ha)) +
  geom_histogram(
    color = "#309898",
    fill = "#F4631E",
    alpha = 0.75,
    bins = 13
  ) +
  scale_x_log10()

Densidad

df_tomate |> 
  ggplot(aes(x = rendimiento_t_ha)) +
  geom_density()

• Mejorar el gráfico anterior:

df_tomate |>
  ggplot(aes(x = rendimiento_t_ha)) +
  geom_density(
    color = "#309898",
    fill = "#F4631E",
    alpha = 0.65
  )

• Agregamos el logaritmo en el eje x:

df_tomate |>
  ggplot(aes(x = rendimiento_t_ha)) +
  geom_density(
    color = "#309898",
    fill = "#F4631E",
    alpha = 0.65
  ) +
  scale_x_log10()

Boxplot

df_tomate |>
  ggplot(aes(x = "", y = rendimiento_t_ha)) +
  geom_boxplot(
    width = 0.4,
    color = "#129990",
    outlier.color = "red",
    fill = "#90D1CA",
    alpha = 0.65
  ) 

• Con transformación logarítmica:

df_tomate |>
  ggplot(aes(x = "", y = rendimiento_t_ha)) +
  geom_boxplot(
    width = 0.4,
    color = "#129990",
    outlier.color = "red",
    fill = "#90D1CA",
    alpha = 0.65
  ) +
  scale_y_log10()

Violín

df_tomate |>
  ggplot(aes(x = "", y = rendimiento_t_ha)) +
  geom_violin(
    color = "#129990",
    fill = "#90D1CA",
    alpha = 0.65
  ) 

• Con transformación logarítmica:

df_tomate |> 
  ggplot(aes(x = "", y = rendimiento_t_ha)) +
  geom_violin(
    color = "#129990",
    fill = "#90D1CA",
    alpha = 0.65
  ) +
  scale_y_log10()

• Violín + Boxplot:

df_tomate |>
  ggplot(aes(x = "", y = rendimiento_t_ha)) +
  geom_violin(color = "#129990",
              fill = "#90D1CA",
              alpha = 0.65) +
  geom_boxplot(
    width = 0.15,
    color = "#096B68",
    fill = "#096B68",
    alpha = 0.25,
    outlier.color = "red"
  )

Agrupadas

Histogramas

• Comparemos la distribución del rendimiento en las dos desagregaciones:

df_tomate |> 
  ggplot(aes(x = rendimiento_t_ha, color = desagregacion_cultivo)) +
  geom_histogram(fill = "white") +
  scale_x_log10()

• Ahora configuremos color y fill:

df_tomate |> 
  ggplot(aes(x = rendimiento_t_ha,
             color = desagregacion_cultivo,
             fill = desagregacion_cultivo)) +
  geom_histogram(alpha = 0.75) +
  scale_x_log10()

• Ahora vamos a utilizar facetas:

media_rto_tomate <-
  df_tomate |> 
  pull(rendimiento_t_ha) |> 
  mean()

df_tomate |> 
  ggplot(aes(x = rendimiento_t_ha,
             color = desagregacion_cultivo,
             fill = desagregacion_cultivo)) +
  geom_histogram(alpha = 0.75, show.legend = FALSE) +
  scale_x_log10() +
  facet_wrap(~desagregacion_cultivo, ncol = 1, scales = "free_y") +
  geom_vline(xintercept = media_rto_tomate, color = "red", linetype = 2)

• Ahora agregamos el departamento al gráfico anterior:

df_tomate |> 
  filter(!departamento %in% c("Amazonas", "Vichada", "Arauca",
                              "Atlántico", "Bolívar", "Chocó",
                              "Córdoba", "La Guajira", "Magdalena",
                              "Sucre")) |> 
  ggplot(aes(x = rendimiento_t_ha,
             color = desagregacion_cultivo,
             fill = desagregacion_cultivo)) +
  geom_histogram(alpha = 0.75, show.legend = FALSE) +
  scale_x_log10() +
  facet_wrap(~departamento+desagregacion_cultivo, ncol = 2, scales = "free_y")

Densidad

• Comparemos la distribución del rendimiento en las dos desagregaciones:

df_tomate |> 
  ggplot(aes(x = rendimiento_t_ha, color = desagregacion_cultivo)) +
  geom_density(fill = "white", alpha = 0.5) +
  scale_x_log10()

• Ahora configuremos color y fill:

df_tomate |> 
  ggplot(aes(x = rendimiento_t_ha,
             color = desagregacion_cultivo,
             fill = desagregacion_cultivo)) +
  geom_density(alpha = 0.75) +
  scale_x_log10()

• Ahora vamos a utilizar facetas:

media_rto_tomate <-
  df_tomate |> 
  pull(rendimiento_t_ha) |> 
  mean()

df_tomate |> 
  ggplot(aes(x = rendimiento_t_ha,
             color = desagregacion_cultivo,
             fill = desagregacion_cultivo)) +
  geom_density(alpha = 0.75, show.legend = FALSE) +
  scale_x_log10() +
  facet_wrap(~desagregacion_cultivo, ncol = 1, scales = "free_y") +
  geom_vline(xintercept = media_rto_tomate, color = "red", linetype = 2)

• Ahora agregamos el departamento al gráfico anterior:

df_tomate |> 
  filter(!departamento %in% c("Amazonas", "Vichada", "Arauca",
                              "Atlántico", "Bolívar", "Chocó",
                              "Córdoba", "La Guajira", "Magdalena",
                              "Sucre")) |> 
  ggplot(aes(x = rendimiento_t_ha,
             color = desagregacion_cultivo,
             fill = desagregacion_cultivo)) +
  geom_density(alpha = 0.75, show.legend = FALSE) +
  scale_x_log10() +
  facet_wrap(~departamento+desagregacion_cultivo, ncol = 2, scales = "free_y")

Boxplot

df_tomate |> 
  ggplot(aes(x = desagregacion_cultivo,
             y = rendimiento_t_ha)) +
  geom_boxplot(color = "#2A4759", fill = "#F79B72")

• El mismo gráfico anterior con logaritmos:

df_tomate |> 
  ggplot(aes(x = desagregacion_cultivo,
             y = rendimiento_t_ha)) +
  geom_boxplot(color = "#2A4759", fill = "#F79B72") +
  scale_y_log10()

• El mismo gráfico anterior con logaritmos y violín:

df_tomate |> 
  ggplot(aes(x = desagregacion_cultivo,
             y = rendimiento_t_ha)) +
  geom_violin() +
  geom_boxplot(color = "#2A4759", fill = "#F79B72", width = 0.25) +
  scale_y_log10()

• Comparando departamentos y desagregación:

df_tomate |> 
  filter(!departamento %in% c("Amazonas", "Vichada", "Arauca",
                              "Atlántico", "Bolívar", "Chocó",
                              "Córdoba", "La Guajira", "Magdalena",
                              "Sucre")) |> 
  ggplot(aes(x = departamento,
             y = rendimiento_t_ha,
             color = desagregacion_cultivo)) +
  geom_boxplot() +
  scale_y_log10() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top")

• Ahora agregamos una cuarta variable al gráfico:

df_tomate |> 
  filter(!departamento %in% c("Amazonas", "Vichada", "Arauca",
                              "Atlántico", "Bolívar", "Chocó",
                              "Córdoba", "La Guajira", "Magdalena",
                              "Sucre")) |> 
  ggplot(aes(x = departamento,
             y = rendimiento_t_ha,
             color = desagregacion_cultivo)) +
  facet_wrap(~ano, ncol = 1) +
  geom_boxplot() +
  scale_y_log10() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top")

Relaciones x vs y

X: tiempo vs Y

df_embalses |> 
  filter(Name == "CHUZA") |> 
  ggplot(aes(x = Date, y = Value)) +
  geom_line()

• Agregamos una línea de tendencia:

df_embalses |> 
  filter(Name == "CHUZA") |> 
  ggplot(aes(x = Date, y = Value)) +
  geom_line() +
  geom_smooth()

• Ahora graficamos todos los embalses (usamos la versión interactiva):

ggplotly(df_embalses |>
           ggplot(aes(
             x = Date, y = Value, color = Name
           )) +
           geom_line())

• El mismo gráfico pero usano áreas:

df_embalses |> 
  filter(Name %in%  c("CHUZA", "AGREGADO BOGOTA")) |> 
  ggplot(aes(x = Date, y = Value, fill = Name)) +
  geom_area()

• ¿Cómo es la distribución mensual en el embalse Chuza?

df_embalses |> 
  filter(Name == "CHUZA") |> 
  mutate(mes = month(Date, label = TRUE, abbr = FALSE)) |> 
  ggplot(aes(x = mes, y = Value)) +
  geom_boxplot()

X: numérica vs Y: numérica

Diagrama de puntos

• ¿Existe alguna relación entre el área perdida y el rendimiento?

df_yuca <-
  datos |> 
  filter(cultivo == "Yuca") |> 
  filter(area_cosechada_ha > 0) |> 
  filter(area_sembrada_ha >= area_cosechada_ha) |> 
  mutate(area_perdida = area_sembrada_ha - area_cosechada_ha,
         porcen_perdida = area_perdida / area_sembrada_ha)

df_yuca |> 
  ggplot(aes(x = porcen_perdida, y = rendimiento_t_ha)) +
  geom_point()

• Ahora graficamos de forma interativa:

ggplotly(
  df_yuca |> 
  ggplot(aes(x = porcen_perdida, y = rendimiento_t_ha)) +
  geom_point()
)

• Ahora agregamos logaritmos en ambos ejes:

df_yuca |> 
  ggplot(aes(x = porcen_perdida, y = rendimiento_t_ha)) +
  geom_point() +
  scale_x_log10() +
  scale_y_log10()

• Ahora agregamos un suavizamiento de tipo lineal (rojo) y no lineal (azul):

df_yuca |> 
  ggplot(aes(x = porcen_perdida, y = rendimiento_t_ha)) +
  geom_point(color = "gray50") +
  scale_x_log10() +
  scale_y_log10() +
  geom_smooth(method = "lm", color = "red") +
  geom_smooth(color = "blue")

• Ahora agregamos facetas para ver el comportamiento en cada año:

df_yuca |> 
  ggplot(aes(x = porcen_perdida, y = rendimiento_t_ha)) +
  facet_wrap(~ano) +
  geom_point(color = "gray50") +
  scale_x_log10() +
  scale_y_log10() +
  geom_smooth(method = "lm", color = "red") +
  geom_smooth(color = "blue")

• ¿Cuál es la relación que existe entre el área sembrada y el rendimiento?

df_yuca |> 
  ggplot(aes(x = area_sembrada_ha, y = rendimiento_t_ha)) +
  geom_point() +
  scale_y_log10() +
  scale_x_log10() +
  geom_smooth(method = "lm", color = "red") +
  geom_smooth(color = "blue")

• Ahora vamos a graficar la relación existente entre el nivel de dos embalses:

df_embalses |> 
  filter(Name %in% c("CHUZA", "RIOGRANDE2")) |> 
  pivot_wider(names_from = Name,
              values_from = Value) |> 
  ggplot(aes(x = CHUZA, y = RIOGRANDE2)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red") +
  geom_smooth(color = "blue")

• Ahora quitamos los puntos y dejamos solo las líneas del suavizamiento:

ggplotly(
  df_embalses |> 
  filter(Name %in% c("CHUZA", "RIOGRANDE2")) |> 
  pivot_wider(names_from = Name,
              values_from = Value) |> 
  ggplot(aes(x = CHUZA, y = RIOGRANDE2)) +
  geom_smooth(method = "lm", color = "red") +
  geom_smooth(color = "blue")
)

• Perfil de rezagos:

df_chuza_rezagos <-
  df_embalses |> 
  filter(Name %in% c("CHUZA")) |> 
  mutate(rezago1 = lag(Value, n = 1),
         rezago7 = lag(Value, n = 7),
         rezago30 = lag(Value, n = 30),
         rezago90 = lag(Value, n = 90),
         rezago180 = lag(Value, n = 180),
         rezago365 = lag(Value, n = 365))

df_chuza_rezagos |> 
  ggplot(aes(x = rezago1, y = Value)) +
  geom_point()

• Rezago 7:

df_chuza_rezagos |> 
  ggplot(aes(x = rezago7, y = Value)) +
  geom_point() +
  geom_smooth(meth = "lm")

• Rezago 30:

df_chuza_rezagos |> 
  ggplot(aes(x = rezago30, y = Value)) +
  geom_point() +
  geom_smooth(method = "lm")

• Rezago 90:

df_chuza_rezagos |> 
  ggplot(aes(x = rezago90, y = Value)) +
  geom_point() +
  geom_smooth(method = "lm") +
  geom_smooth(color = "red")

• Rezago 180:

df_chuza_rezagos |> 
  ggplot(aes(x = rezago180, y = Value)) +
  geom_point() +
  geom_smooth(method = "lm") +
  geom_smooth(color = "red")

• Rezago 365:

df_chuza_rezagos |> 
  ggplot(aes(x = rezago365, y = Value)) +
  geom_point() +
  geom_smooth(method = "lm") +
  geom_smooth(color = "red")

• Todos los rezagos en un mismo gráfico:

df_chuza_rezagos |>
  pivot_longer(cols = rezago1:rezago365, values_to = "rezago") |>
  mutate(name = factor(
    name,
    levels = c(
      "rezago1",
      "rezago7",
      "rezago30",
      "rezago90",
      "rezago180",
      "rezago365"
    )
  )) |>
  ggplot(aes(x = rezago, y = Value)) +
  facet_wrap( ~ name) +
  geom_point() +
  geom_smooth(method = "lm") +
  geom_smooth(color = "red")

Densidad 2D

df_yuca |> 
  ggplot(aes(x = porcen_perdida, y = rendimiento_t_ha)) +
  geom_bin_2d() +
  scale_x_log10() +
  scale_y_log10() +
  geom_smooth(method = "lm", color = "red") +
  geom_smooth(color = "blue")

Matriz de correlaciones

• En este caso queremos relacionar el nivel de todos los embalses.

df_embalses_ancho <-
  df_embalses |> 
  pivot_wider(names_from = Name,
              values_from = Value)

df_embalses_ancho

corrplot

• El primer paso es calcular la matriz de correlaciones.

mtx_cor1 <-
  df_embalses_ancho |> 
  select(where(is.numeric)) |> 
  cor(use = "pairwise.complete.obs")

mtx_cor1 |> head()

                AGREGADO BOGOTA ALTOANCHICAYA      AMANI    BETANIA     CALIMA1
AGREGADO BOGOTA      1.00000000  -0.105853184 0.04548885 0.05313706 0.248772126
ALTOANCHICAYA       -0.10585318   1.000000000 0.10083206 0.10618523 0.008379681
AMANI                0.04548885   0.100832058 1.00000000 0.01802802 0.638348426
BETANIA              0.05313706   0.106185231 0.01802802 1.00000000 0.113169192
CALIMA1              0.24877213   0.008379681 0.63834843 0.11316919 1.000000000
CHUZA               -0.12312142   0.140558590 0.10272721 0.05856395 0.210925440
                      CHUZA   ESMERALDA     GUAVIO   MIRAFLORES      MUNA
AGREGADO BOGOTA -0.12312142  0.43810069 0.35097490  0.441335542 0.2169859
ALTOANCHICAYA    0.14055859  0.05849095 0.09355349  0.028559484 0.1941361
AMANI            0.10272721  0.33933476 0.37151014  0.606414433 0.1123771
BETANIA          0.05856395 -0.02966879 0.06579716 -0.007968756 0.2672809
CALIMA1          0.21092544  0.22801642 0.38234049  0.648877898 0.2304914
CHUZA            1.00000000  0.57366196 0.69027732  0.359097272 0.1048208
                     PENOL      PLAYAS    PORCE II   PORCE III       PRADO
AGREGADO BOGOTA 0.19131187 -0.10522977 -0.15738029 -0.08260161 -0.12356294
ALTOANCHICAYA   0.09010404  0.15362189  0.12866931  0.11896751  0.12739429
AMANI           0.67681615  0.36653419  0.37172095  0.25907644  0.49303076
BETANIA         0.07741600  0.02425467  0.05320275  0.07128835  0.26110821
CALIMA1         0.71753130  0.40784527  0.16766161  0.15065918  0.63155955
CHUZA           0.36863896  0.19645910  0.07382421  0.10676375 -0.04835912
                    PUNCHINA RIOGRANDE2   SALVAJINA SAN LORENZO   TRONERAS
AGREGADO BOGOTA -0.082918196 0.30748849  0.05111805  0.21820352 -0.1388754
ALTOANCHICAYA    0.166095691 0.09794493  0.03054739  0.09221607  0.1895705
AMANI            0.252537732 0.75529761  0.29878376  0.69847895  0.3255307
BETANIA          0.001769815 0.12406871  0.18151679  0.01120944  0.1222360
CALIMA1          0.130732601 0.62240825  0.40582078  0.45731070  0.2439332
CHUZA            0.171145666 0.20545636 -0.08758519  0.37546900  0.3164525
                     URRA1   TOPOCORO   EL QUIMBO   ITUANGO
AGREGADO BOGOTA 0.30581672 0.25184288 0.287782514 0.1054153
ALTOANCHICAYA   0.07017185 0.09601921 0.003257804 0.1705238
AMANI           0.63661494 0.83301216 0.375878584 0.3239157
BETANIA         0.03466484 0.08668939 0.179215561 0.3403021
CALIMA1         0.41519563 0.69511287 0.464293658 0.3955245
CHUZA           0.29185387 0.23719174 0.461513334 0.4277711

• Después graficamos el correlograma:

library(corrplot)
mtx_cor1 |> 
  corrplot()

• Ahora mejoramos el gráfico:

mtx_cor1 |> 
  corrplot(
    diag = FALSE,
    type = "lower",
    tl.cex = 0.5,
    tl.col = "black",
    tl.srt = 45,
    method = "pie",
    order = "hclust"
  )

Gráfico de redes

• Necesitamos calcular la matriz de correlaciones con la función correlate():

library(corrr)

mtx_cor2 <-
  df_embalses_ancho |> 
  select(where(is.numeric)) |> 
  correlate()

mtx_cor2

• Ahora realizamos el gráfico:

mtx_cor2 |> 
  network_plot()

Incertidumbre

Bandas confianza

• Lo primero que hacemos es resumir la serie de forma mensual.

library(tsibble)

df_betania <-
  df_embalses |> 
  filter(Name == "BETANIA") |> 
  mutate(fecha_mes = yearmonth(Date)) |> 
  group_by(fecha_mes) |> 
  reframe(promedio = mean(Value, na.rm = TRUE),
          mediana = median(Value, na.rm = TRUE),
          p25 = quantile(Value, probs = 0.25, na.rm = TRUE),
          p75 = quantile(Value, probs = 0.75, na.rm = TRUE),
          desviacion = sd(Value, na.rm = TRUE))

df_betania

• Podemos graficar la incertidumbre con la mediana y los percentiles:

df_betania |> 
  ggplot(aes(x = fecha_mes, y = mediana, ymin = p25, ymax = p75)) +
  geom_line() +
  geom_ribbon()

• Mejoramos el gráfico anterior:

df_betania |> 
  ggplot(aes(x = fecha_mes, y = mediana, ymin = p25, ymax = p75)) +
  geom_ribbon(fill = "#B2D8CE") +
  geom_line(color = "#5459AC")

• Ahora usemos el promedio y la desviación estándar:

df_betania |> 
  ggplot(aes(x = fecha_mes,
             y = promedio,
             ymin = promedio - desviacion,
             ymax = promedio + desviacion)) +
  geom_ribbon(fill = "#B2D8CE") +
  geom_line(color = "#5459AC")

• Ahora vamos a representar +1, -1, +2 y -2 desviaciones estándar:

grafico_incertidumbre <-
  df_betania |>
  ggplot(aes(x = fecha_mes, y = promedio)) +
  geom_ribbon(aes(
    ymin = promedio - (desviacion * 2),
    ymax = promedio + (desviacion * 2)
  ),
  fill = "#5459AC",
  alpha = 0.5) +
  geom_ribbon(
    aes(ymin = promedio - desviacion, ymax = promedio + desviacion),
    fill = "#52357B",
    alpha = 0.65
  ) +
  geom_line(color = "#648DB3")

grafico_incertidumbre

• Ahora le damos interactividad:

ggplotly(grafico_incertidumbre)

Barras de error

Barras

• Primero conformamos los datos:

df_limites_cafe <-
  datos |> 
  filter(cultivo == "Café") |> 
  group_by(ano) |> 
  reframe(promedio = mean(rendimiento_t_ha),
          desviacion = sd(rendimiento_t_ha)) |> 
  mutate(lim_inferior = promedio - desviacion,
         lim_superior = promedio + desviacion)

df_limites_cafe

• Ahora graficamos las barras de error:

df_limites_cafe |> 
  ggplot(aes(x = ano, y = promedio, ymin = lim_inferior, ymax = lim_superior)) +
  geom_col(fill = "#533B4D") +
  geom_errorbar(width = 0.15, color = "#F564A9")

Puntos + Línea de error

df_limites_cafe |> 
  ggplot(aes(x = ano, y = promedio, ymin = lim_inferior, ymax = lim_superior)) +
  geom_errorbar(width = 0.15, color = "#F564A9") +
  geom_point(color = "#533B4D")

• Otra opción es con pointrange:

df_limites_cafe |> 
  ggplot(aes(x = ano, y = promedio, ymin = lim_inferior, ymax = lim_superior)) +
  geom_pointrange()