Visualizaciones

Estadística

Bibliotecas

library(tidyverse)
library(janitor)
library(readxl)

library(ggsci) # colores preestablecidos
library(plotly) # Visualización interactiva

Datos

Nivel de embalses

df_embalses <- read_csv("datos-ejemplos/PorcVoluUtilDiar.csv")
df_embalses |> head()

Evaluaciones agropecuarias

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

df_evas |> head()

Créditos agropecuarios

df_creditos <- read_csv("datos-ejemplos/Colocaciones_de_Cr_dito_Sector_Agropecuario_-_2021-_2024_20240910.csv") |> 
  clean_names()

df_creditos |> head()

Cantidades

Créditos agropecuarios por año

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

• Mejorando la estética del gráfico:

df_creditos |> 
  count(ano, name = "total") |> 
  ggplot(aes(x = ano, y = total)) +
  geom_col(color = "cornflowerblue", fill = "#F87A53", alpha = 0.5) +
  labs(x = "Año",
       y = "Total (n)",
       title = "Créditos agropecuarios en Colombia",
       subtitle = "Años 2021 a 2024",
       caption = "*2024: información incompleta") +
  theme_bw()

• Agregando texto a las barras - componente global:

df_creditos |> 
  count(ano, name = "total") |> 
  ggplot(aes(x = ano, y = total, label = total)) +
  geom_col(color = "cornflowerblue", fill = "#F87A53", alpha = 0.5) +
  geom_text() +
  labs(x = "Año",
       y = "Total (n)",
       title = "Créditos agropecuarios en Colombia",
       subtitle = "Años 2021 a 2024",
       caption = "*2024: información incompleta") +
  theme_bw()

• Agregando texto a las barras - componente específico:

df_creditos |> 
  count(ano, name = "total") |> 
  ggplot(aes(x = ano, y = total)) +
  geom_col(color = "cornflowerblue", fill = "#F87A53", alpha = 0.5) +
  geom_label(aes(label = total), color = "white", fill = "#3B1E54") +
  labs(x = "Año",
       y = "Total (n)",
       title = "Créditos agropecuarios en Colombia",
       subtitle = "Años 2021 a 2024",
       caption = "*2024: información incompleta") +
  theme_bw()

EVAs por departamento

df_evas |> 
  count(departamento, name = "total") |> 
  ggplot(aes(x = departamento, y = total)) +
  geom_col(color = "cornflowerblue", fill = "#F87A53", alpha = 0.5) +
  geom_label(aes(label = total), color = "white", fill = "#3B1E54") +
  labs(x = "Año",
       y = "Total (n)",
       title = "Evaluaciones agropecuarias por departamento",
       subtitle = "Colombia 2019 a 2023") +
  theme_bw()

• Aplicando mejoras al gráfico:

df_evas |> 
  count(departamento, name = "total") |> 
  ggplot(aes(x = departamento, y = total)) +
  geom_col(color = "cornflowerblue", fill = "#F87A53", alpha = 0.5) +
  geom_label(aes(label = total), color = "white", fill = "#3B1E54", size = 1.7) +
  labs(x = "Año",
       y = "Total (n)",
       title = "Evaluaciones agropecuarias por departamento",
       subtitle = "Colombia 2019 a 2023") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

• Rotando las coordenas del gráfico:

df_evas |> 
  count(departamento, name = "total") |> 
  ggplot(aes(x = reorder(departamento, total), y = total)) +
  geom_col(color = "cornflowerblue", fill = "#F87A53", alpha = 0.5) +
  geom_label(aes(label = total), color = "white", fill = "#3B1E54", size = 1.7) +
  labs(x = "Año",
       y = "Total (n)",
       title = "Evaluaciones agropecuarias por departamento",
       subtitle = "Colombia 2019 a 2023") +
  theme_bw() +
  coord_flip()

Mediana por embalse

df_embalses |> 
  group_by(Name) |> 
  reframe(mediana = median(Value)) |> 
  ggplot(aes(x = reorder(Name, mediana), y = mediana)) +
  geom_point(color = "red") +
  coord_flip() +
  labs(x = "")

• Agregando la mediana a nivel de Colombia:

mediana_general_vol <- df_embalses$Value |> median()

df_embalses |> 
  group_by(Name) |> 
  reframe(mediana = median(Value)) |> 
  ggplot(aes(x = reorder(Name, mediana), y = mediana)) +
  geom_point(color = "red") +
  labs(x = "",
       caption = "Línea azul: mediana general") +
  geom_hline(yintercept = mediana_general_vol,
             color = "blue",
             linetype = 2) +
  coord_flip()

Facetas: créditos por año y departamento

df_creditos |> 
  count(departamento_inversion, ano, name = "total") |> 
  ggplot(aes(x = ano, y = total)) +
  facet_wrap(~ departamento_inversion, ncol = 3, scales = "free_y") +
  geom_col()

Agregando una tercera variable (1/2)

• ¿Cómo es la distribución de créditos por tipo de productor, género y año?

df_creditos |> 
  count(ano, genero, tipo_productor) |> 
  ggplot(aes(x = ano, y = n, color = genero, fill = genero)) +
  facet_wrap(~ tipo_productor, scales = "free_y") +
  geom_col(alpha = 0.5)

• Mejoramos el gráfico anterior:

df_creditos |> 
  count(ano, genero, tipo_productor) |> 
  ggplot(aes(x = ano, y = n, color = genero, fill = genero)) +
  facet_wrap(~ tipo_productor, scales = "free_y") +
  geom_col(alpha = 0.5, position = "dodge")

• Ajustamos colores de forma manual:

df_creditos |> 
  count(ano, genero, tipo_productor) |> 
  ggplot(aes(x = ano, y = n, color = genero, fill = genero)) +
  facet_wrap(~ tipo_productor, scales = "free_y") +
  geom_col(alpha = 0.5, position = "dodge") +
  scale_color_manual(values = c("orangered", "magenta2", "gray40")) +
  scale_fill_manual(values = c("orangered", "magenta2", "gray40"))

• Ajustamos detalles estéticos del gráfico anterior:

df_creditos |> 
  count(ano, genero, tipo_productor) |> 
  ggplot(aes(x = ano, y = n, color = genero, fill = genero)) +
  facet_wrap(~ tipo_productor, scales = "free_y") +
  geom_col(alpha = 0.5, position = "dodge") +
  scale_color_npg() +
  scale_fill_npg() +
  theme(legend.position = "top") + #left top right bottom
  labs(x = "Año",
       y = "Total (n)",
       title = "Créditos agropecuarios en Colombia por producto y género",
       subtitle = "Años 2021 a 2024",
       caption = "*2024: información incompleta",
       color = "Género",
       fill = "Género")

Agregando una tercera variable (2/2)

• ¿Cómo es la distribución del valor de la inversión en créditos por tipo de productor, género y año?

df_creditos |> 
  group_by(ano, genero, tipo_productor) |> 
  reframe(total = sum(valor_inversion, na.rm = TRUE)) |> 
  ggplot(aes(x = ano, y = total, color = genero, fill = genero)) +
  facet_wrap(~ tipo_productor, scales = "free_y") +
  geom_col(alpha = 0.5, position = "dodge")

• ¿Cómo es la distribución del valor promedio de la inversión en créditos por tipo de productor, género y año?

df_creditos |> 
  group_by(ano, genero, tipo_productor) |> 
  reframe(total = sum(valor_inversion, na.rm = TRUE),
          creditos = n()) |> 
  mutate(promedio = total / creditos) |> 
  ggplot(aes(x = ano, y = promedio, color = genero, fill = genero)) +
  facet_wrap(~ tipo_productor, scales = "free_y") +
  geom_col(alpha = 0.5, position = "dodge")

Gráfico interactivo

ggplotly(
  df_creditos |> 
  count(ano, name = "total") |> 
  ggplot(aes(x = ano, y = total)) +
  geom_col()
)

Proporciones

Distribución de créditos agropecuarios

• ¿Cómo es la distribución de créditos por tipo de productor, género y año?

df_creditos |> 
  count(ano, genero, tipo_productor) |> 
  ggplot(aes(x = ano, y = n, color = genero, fill = genero)) +
  facet_wrap(~ tipo_productor, scales = "free_y") +
  geom_col(alpha = 0.5, position = "fill") +
  coord_flip()

• ¿Cómo es la distribución de créditos por departamento, tipo de productor, género y año?

df_creditos |> 
  count(ano, genero, tipo_productor, departamento_inversion) |> 
  ggplot(aes(x = ano, y = n, color = genero, fill = genero)) +
  facet_wrap(~ departamento_inversion ~ tipo_productor, ncol = 4) +
  geom_col(position = "fill", alpha = 0.5) +
  theme(legend.position = "top") +
  coord_flip()

• El mismo gráfico anterior cambiando el tipo de producto por género:

df_creditos |> 
  count(ano, genero, tipo_productor, departamento_inversion) |> 
  ggplot(aes(x = ano, y = n, color = tipo_productor, fill = tipo_productor)) +
  facet_wrap(~ departamento_inversion ~ genero, ncol = 3) +
  geom_col(position = "fill", alpha = 0.5) +
  theme(legend.position = "top") +
  coord_flip()

• El mismo gráfico anterior utilizando otro tipo de faceta:

df_creditos |> 
  count(ano, genero, tipo_productor, departamento_inversion) |> 
  ggplot(aes(x = ano, y = n, color = genero, fill = genero)) +
  facet_grid(~ departamento_inversion ~ tipo_productor) +
  geom_col(position = "fill", alpha = 0.5) +
  theme(legend.position = "top") +
  coord_flip()

EVAs

• ¿Cómo es la proporción de evaluaciones agropecuarias por departamento y año que han superado la media histórica del rendimiento en café?

# Filtramos el cultivo de interés
df_eva_cultivo <-
  df_evas |> 
  filter(cultivo == "Café")

# Calculamos primero la media histórica del rendimiento
media_rto_cultivo <-
  df_eva_cultivo$rendimiento_t_ha |> mean()

# Finalmente extraemos la variable "resultado" que indica si 
# el rendimiento de esa evaluación agropecuaria es superior o inferior al
# promedio histórico
df_resumen_rto <-
  df_eva_cultivo |> 
  mutate(
    resultado = if_else(
      condition = rendimiento_t_ha > media_rto_cultivo,
      true = "Superior",
      false = "Inferior"
    )
  )

df_resumen_rto |> 
  count(departamento, ano, resultado) |> 
  ggplot(aes(x = ano, y = n, color = resultado, fill = resultado)) +
  facet_wrap(~ departamento, ncol = 4) +
  geom_col(position = "fill", alpha = 0.5) +
  coord_flip() +
  theme(legend.position = "top")

Embalses

• ¿Qué proporción del tiempo los embalses superan el 50% del nivel en cada año?

df_embalses2 <-
  df_embalses |>
  mutate(
    resultado = if_else(
      condition = Value > 0.5,
      true = "Superior",
      false = "Inferior"
    ),
    year_es = year(Date)
  )

df_embalses2 |> 
  count(year_es, resultado) |> 
  mutate(year_es = as.factor(year_es)) |> 
  ggplot(aes(x = year_es, y = n, fill = resultado)) +
  geom_col(position = "fill") +
  coord_flip()

• Ahora observamos el mismo gráfico a nivel de embalse:

df_embalses2 |> 
  count(year_es, resultado, Name) |> 
  mutate(year_es = as.factor(year_es)) |> 
  ggplot(aes(x = year_es, y = n, fill = resultado)) +
  geom_col(position = "fill") +
  coord_flip() +
  facet_wrap(~ Name) +
  theme(legend.position = "top")

Gráfico interactivo

ggplotly(
 df_creditos |> 
  count(ano, genero, tipo_productor) |> 
  ggplot(aes(x = ano, y = n, color = genero, fill = genero)) +
  facet_wrap(~ tipo_productor, scales = "free_y") +
  geom_col(alpha = 0.5, position = "fill") +
  coord_flip() 
)

Distribuciones

Individuales - Histogramas

• ¿Cómo es la distribución del nivel de los embalses en Colombia?

df_embalses |> 
  ggplot(aes(x = Value)) +
  geom_histogram(color = "red",
                 fill = "gray70",
                 bins = 30)

• ¿La distribución observada a nivel nacional sigue el mismo patrón en los embalses?

df_embalses |> 
  ggplot(aes(x = Value)) +
  geom_histogram(color = "red",
                 fill = "gray70",
                 bins = 30) +
  facet_wrap(~ Name, ncol = 1)

Individuales - Densidades

df_embalses |>
  ggplot(aes(x = Value)) +
  geom_density(color = "red",
               fill = "gray70")

Individuales - Boxplot

df_embalses |>
  ggplot(aes(x = Value)) +
  geom_boxplot(color = "red", fill = "gray70")

Comparativas - Boxplot

df_embalses |>
  ggplot(aes(x = Name, y = Value)) +
  geom_boxplot(color = "red", fill = "gray70") +
  coord_flip()

Comparativas - Densidades

• Una primera opción es graficar todas las densidades en un solo gráfico:

df_embalses |> 
  ggplot(aes(x = Value, color = Name)) +
  geom_density()

• Agregando interactividad al gráfico anterior:

ggplotly(
  df_embalses |> 
  ggplot(aes(x = Value, color = Name)) +
  geom_density()
)

• Como son muchas etiquetas (embalses) optamos por utilizar facetas:

df_embalses |> 
  ggplot(aes(x = Value, color = Name)) +
  geom_density(show.legend = FALSE) +
  facet_wrap(~Name, ncol = 4, scales = "free_y")

Relaciones X-Y

Diagrama de dispersión

• ¿Existe alguna asocición entre el área perdida y el rendimiento para el cultivo de cacao en Colombia?

# Primero calculamos el área perdida y la expresamos en porcentaje
df_eva_cacao <-
  df_evas |> 
  filter(cultivo == "Cacao") |> 
  mutate(area_perdida_porc = (area_sembrada_ha - area_cosechada_ha) / area_sembrada_ha) |> 
  filter(area_perdida_porc > 0)

df_eva_cacao |> 
  ggplot(aes(x = area_perdida_porc, y = rendimiento_t_ha)) +
  geom_point() +
  geom_smooth(method = "lm")

• ¿Hay alguna asociación entre el área que se siembra y el rendimiento?

df_eva_cacao |> 
  ggplot(aes(x = area_sembrada_ha, y = rendimiento_t_ha)) +
  geom_point() +
  scale_x_log10() +
  geom_smooth(method = "lm")

• Ahora calculamos la correlación de estas dos variables:

cor(x = df_eva_cacao$area_sembrada_ha,
    y = df_eva_cacao$rendimiento_t_ha,
    use = "pairwise.complete.obs")

[1] 0.02916119

• Ahora calculamos la misma correlación de estas dos variables pero aplicando la transformación logarítmica sobre ambas:

cor(x = log1p(df_eva_cacao$area_sembrada_ha),
    y = log1p(df_eva_cacao$rendimiento_t_ha),
    use = "pairwise.complete.obs")

[1] 0.209491

Diagrama de dispersión (tiempo)

• Primero grafiquemos un solo embalse:

df_embalses |> 
  filter(Name == "AGREGADO BOGOTA") |>
  ggplot(aes(x = Date, y = Value)) +
  geom_line() +
  geom_smooth(color = "red")

• ¿Existe alguna tendencia a nivel nacional en los niveles de los embalses?

df_embalses |> 
  ggplot(aes(x = Date, y = Value)) +
  geom_bin_2d() +
  geom_smooth(color = "red")

Correlaciones (1/2)

• Cálculo de correlaciones:

matriz_cor1 <-
  df_embalses |>
  pivot_wider(names_from = Name, values_from = Value) |>
  select(where(is.numeric)) |>
  cor(use = "pairwise.complete.obs")

• Gráfico inicial:

library(corrplot)
matriz_cor1 |> 
  corrplot()

• Mejorando el gráfico:

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

Correlaciones (2/2)

library(corrr)

mtx_cor2 <-
  df_embalses |>
  pivot_wider(names_from = Name, values_from = Value) |>
  select(where(is.numeric)) |> 
  correlate()

mtx_cor2 |> network_plot()

Incertidumbre

Barras de errores

df_eva_cacao |> 
  group_by(departamento) |> 
  reframe(promedio = mean(rendimiento_t_ha, na.rm = TRUE),
          desviacion = sd(rendimiento_t_ha, na.rm = TRUE)) |> 
  ggplot(aes(x = departamento, y = promedio,
             ymin = promedio - desviacion,
             ymax = promedio + desviacion)) +
  geom_point() +
  geom_errorbar() +
  coord_flip()

• Otra opción:

df_eva_cacao |> 
  group_by(departamento) |> 
  reframe(promedio = mean(rendimiento_t_ha, na.rm = TRUE),
          desviacion = sd(rendimiento_t_ha, na.rm = TRUE)) |> 
  ggplot(aes(x = departamento, y = promedio,
             ymin = promedio - desviacion,
             ymax = promedio + desviacion)) +
  geom_pointrange() +
  coord_flip()

Mapas

library(geodata)
library(sf)
library(ggspatial)

Tipos de datos

Extracción de mapas con R

• Colombia:

colombia_pais <- gadm(country = "COL", level = 0, path = "datos-ejemplos/")
colombia_pais |> plot()

• Colombia - Departamentos:

colombia_deptos <- gadm(country = "COL", level = 1, path = "datos-ejemplos/")
colombia_deptos |> plot()

• Colombia - Municipios:

colombia_mpios <- gadm(country = "COL", level = 2, path = "datos-ejemplos/")
colombia_mpios |> plot()

Lectura de mapas en R

mapa_colombia <- readRDS("datos-ejemplos/gadm/gadm41_COL_0_pk.rds")
mapa_colombia |> plot()

Colombia: departamentos .shp

colombia_deptos_shp <- 
  st_read("datos-ejemplos/MGN2021_DPTO_POLITICO/MGN_DPTO_POLITICO.shp") |> 
  mutate(DPTO_CCDGO = as.numeric(DPTO_CCDGO))

Reading layer `MGN_DPTO_POLITICO' from data source 
  `D:\Otros\UdeA\2024-02\01-estadistica\estadistica-202402\datos-ejemplos\MGN2021_DPTO_POLITICO\MGN_DPTO_POLITICO.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 33 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.73562 ymin: -4.229406 xmax: -66.84722 ymax: 13.39473
Geodetic CRS:  MAGNA-SIRGAS

colombia_deptos_shp |>
  ggplot() +
  geom_sf()

Colombia: municipios .shp

colombia_mpios_shp <- st_read("datos-ejemplos/MGN2021_MPIO_POLITICO/MGN_MPIO_POLITICO.shp")

Reading layer `MGN_MPIO_POLITICO' from data source 
  `D:\Otros\UdeA\2024-02\01-estadistica\estadistica-202402\datos-ejemplos\MGN2021_MPIO_POLITICO\MGN_MPIO_POLITICO.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1121 features and 12 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.73562 ymin: -4.229406 xmax: -66.84722 ymax: 13.39473
Geodetic CRS:  MAGNA-SIRGAS

colombia_mpios_shp |>
  ggplot() +
  geom_sf()

• Puedo aplicar un filtro tal cual como lo hemos hecho previamente:

colombia_mpios_shp |> 
  filter(DPTO_CNMBR == "ANTIOQUIA") |> 
  ggplot() +
  geom_sf()

EVAs por departamento - Maíz

• Filtramos el maíz y cambiamos el nombre de la variable codigo_dane_departamento

df_eva_maiz <-
  df_evas |> 
  filter(cultivo == "Maíz") |> 
  mutate(codigo_dane_departamento = as.numeric(codigo_dane_departamento)) |> 
  rename(DPTO_CCDGO = codigo_dane_departamento)

• Resumimos la información por departamento:

df_resumen_maiz <-
  df_eva_maiz |>
  group_by(DPTO_CCDGO) |>
  reframe(promedio_rto = mean(rendimiento_t_ha, na.rm = TRUE))

df_resumen_maiz

• Unimos el resumen con el mapa:

df_maiz_mapa <-
  colombia_deptos_shp |> 
  left_join(df_resumen_maiz, by = "DPTO_CCDGO")

df_maiz_mapa

• Graficamos el mapa con el rendimiento promedio del maíz por departamento:

df_maiz_mapa |> 
  ggplot(mapping = aes(fill = promedio_rto)) +
  geom_sf() +
  scale_fill_viridis_c() +
  theme_void()  +
  annotation_north_arrow(location = "tl") +
  labs(title = "Rendimiento del maíz en Colombia",
       subtitle = "Departamental",
       fill = "Rendimiento (t/ha)")

Créditos agropecuarios por departamento

• Filtramos el año 2023 y utilizamos para el ejemplo Capital de trabajo Microcrédito rural (165000)

df_filtro_creditos <-
  df_creditos |> 
  filter(ano == 2023) |> 
  filter(id_rubro == 165000)

df_filtro_creditos |> head()

• Ahora graficamos los puntos (latitud y longitud) en un mapa estático:

colombia_deptos_shp |>
  ggplot() +
  geom_sf() +
  geom_point(data = df_filtro_creditos,
             mapping = aes(x = longitud, y = latitud, color = valor_inversion)) +
  scale_color_viridis_c()  +
  theme_void()  +
  annotation_north_arrow(location = "tl") +
  labs(title = "Créditos agropecuarios en Colombia - Año 2023",
       subtitle = "Capital de trabajo microcrédito rural",
       color = "Valor (COP)")

Imagen raster

library(terra)
elevacion <- rast("datos-ejemplos/wc2.1_2.5m_elev/wc2.1_2.5m_elev.tif")
elevacion |> plot()

Extracción de información desde raster

elevacion_colombia <-
  elevacion |> 
  mask(colombia_deptos_shp) |> 
  crop(colombia_deptos_shp)

elevacion_colombia |> plot()

Referencias

• Paletas de colores en R