Bibliotecas y configuración

  • Configuramos el tema por defecto para todos nuestros gráficos con la función theme_set()
library(tidyverse)
library(ggpubr)
library(plotly)
library(corrplot)
library(corrr)

library(sf)
library(ggspatial)

theme_set(theme_bw())

Datos

# Nombres de variables
nombres <- c(
  "depto",
  "mpio",
  "cultivo",
  "estado",
  "tiempo_estab",
  "topografia",
  "drenaje",
  "riego",
  "fertilizantes",
  "fecha_analisis",
  "ph_agua_suelo",
  "materia_org",
  "fosforo",
  "azufre",
  "acidez",
  "aluminio",
  "calcio",
  "magnesio",
  "potasio",
  "sodio",
  "cice",
  "conductividad",
  "hierro_olsen",
  "cobre",
  "manganeso",
  "zinc",
  "boro",
  "hierro_doble",
  "cobre_doble",
  "manganeso_doble",
  "zinc_doble"
  
)

# Datos de suelos cultivo de cacao
datos_cacao <-
  read_csv("datos/Resultados_de_An_lisis_de_Laboratorio_Suelos_en_Colombia.csv",
           na = "ND") %>%
  select(-c(numfila, Secuencial)) %>%
  set_names(nombres) %>%
  mutate(
    across(c(depto, mpio), str_to_title),
    across(
      c(cultivo, estado, tiempo_estab, topografia, drenaje, riego),
      str_to_sentence
    ),
    across(
      c(
        fosforo,
        calcio,
        magnesio,
        potasio,
        sodio,
        hierro_olsen,
        cobre,
        manganeso,
        zinc,
        cobre_doble,
        manganeso_doble,
        zinc_doble
      ),
      ~ str_replace_all(
        string = .,
        pattern = ",",
        replacement = "."
      )
    ),
    across(
      c(
        fosforo,
        calcio,
        magnesio,
        potasio,
        sodio,
        hierro_olsen,
        cobre,
        manganeso,
        zinc,
        cobre_doble,
        manganeso_doble,
        zinc_doble
      ),
      ~ str_replace_all(
        string = .,
        pattern = "<",
        replacement = ""
      )
    ),
    across(
      c(
        fosforo,
        calcio,
        magnesio,
        potasio,
        sodio,
        hierro_olsen,
        cobre,
        manganeso,
        zinc,
        cobre_doble,
        manganeso_doble,
        zinc_doble
      ),
      as.numeric
    )
  ) %>% 
  select(-c(fecha_analisis, hierro_doble, cobre_doble, manganeso_doble, zinc_doble)) %>% 
  filter(cultivo == "Cacao") %>% 
  filter(ph_agua_suelo <= 10)

datos_cacao %>% head()

Gráficos de dispersión

  • pH vs calcio:
datos_cacao %>% 
  ggplot(aes(x = ph_agua_suelo, y = calcio)) +
  geom_point()

  • Podemos mejorar el gráfico anterior usando logaritmos para los valores del eje Y:
datos_cacao %>% 
  ggplot(aes(x = ph_agua_suelo, y = calcio)) +
  geom_point() +
  scale_y_log10()

  • Podemos agregar una línea de tendencia (modelo lineal) a la relación del pH y el calcio:
grafico1 <-
  datos_cacao %>%
  ggplot(aes(x = ph_agua_suelo, y = calcio)) +
  geom_point() +
  scale_y_log10() +
  geom_smooth(method = "lm")

grafico1

  • Veamos el gráfico anterior de forma interactiva:
ggplotly(grafico1)
  • Podemos controlar que el color de los puntos esté determinado por una variable cualitativa de interés:
grafico2 <-
  datos_cacao %>%
  ggplot(aes(x = ph_agua_suelo, y = calcio, color = estado)) +
  geom_point(alpha = 0.1) +
  scale_y_log10() +
  geom_smooth(method = "lm")

grafico2

  • Podemos tener mayor control del gráfico si lo usamos de forma interactiva:
ggplotly(grafico2)

Gráfico cuantil-cuantil

  • Primera opción con ggplt2:
datos_cacao %>% 
  ggplot(aes(sample = ph_agua_suelo)) +
  geom_qq_line() +
  geom_qq()

datos_cacao %>% 
  ggplot(aes(sample = calcio)) +
  geom_qq_line() +
  geom_qq()

  • Otra opción con ggpubr:
ggqqplot(datos_cacao$ph_agua_suelo)

ggqqplot(datos_cacao$calcio)

Correlaciones

Correlación paramétrica Person

cor(datos_cacao$ph_agua_suelo, datos_cacao$calcio, method = "pearson")
## [1] 0.6756225

Correlación no paramétrica de Spearman

cor(datos_cacao$ph_agua_suelo, datos_cacao$calcio, method = "spearman")
## [1] 0.7760956

Matriz de correlaciones

  • Como no estamos totalmente seguros de que nuestras variables numéricas sigan una distribución normal, vamos a usar correlaciones no paramétricas de Spearman.
datos_cacao %>% 
  select(where(is.numeric)) %>% 
  cor(method = "spearman")
##               ph_agua_suelo materia_org    fosforo azufre acidez aluminio
## ph_agua_suelo     1.0000000 -0.23044439  0.3484182     NA     NA       NA
## materia_org      -0.2304444  1.00000000  0.1198217     NA     NA       NA
## fosforo           0.3484182  0.11982170  1.0000000     NA     NA       NA
## azufre                   NA          NA         NA      1     NA       NA
## acidez                   NA          NA         NA     NA      1       NA
## aluminio                 NA          NA         NA     NA     NA        1
## calcio            0.7760956 -0.07928854  0.3232576     NA     NA       NA
## magnesio          0.6861372 -0.19930419  0.1656682     NA     NA       NA
## potasio                  NA          NA         NA     NA     NA       NA
## sodio             0.2290031 -0.27853424 -0.0260695     NA     NA       NA
## cice              0.5204852 -0.02975566  0.2124675     NA     NA       NA
## conductividad            NA          NA         NA     NA     NA       NA
## hierro_olsen             NA          NA         NA     NA     NA       NA
## cobre                    NA          NA         NA     NA     NA       NA
## manganeso                NA          NA         NA     NA     NA       NA
## zinc                     NA          NA         NA     NA     NA       NA
## boro                     NA          NA         NA     NA     NA       NA
##                    calcio   magnesio potasio      sodio        cice
## ph_agua_suelo  0.77609555  0.6861372      NA  0.2290031  0.52048519
## materia_org   -0.07928854 -0.1993042      NA -0.2785342 -0.02975566
## fosforo        0.32325764  0.1656682      NA -0.0260695  0.21246753
## azufre                 NA         NA      NA         NA          NA
## acidez                 NA         NA      NA         NA          NA
## aluminio               NA         NA      NA         NA          NA
## calcio         1.00000000  0.8575178      NA  0.2618703  0.85477166
## magnesio       0.85751778  1.0000000      NA  0.3260122  0.77706855
## potasio                NA         NA       1         NA          NA
## sodio          0.26187035  0.3260122      NA  1.0000000  0.27572723
## cice           0.85477166  0.7770685      NA  0.2757272  1.00000000
## conductividad          NA         NA      NA         NA          NA
## hierro_olsen           NA         NA      NA         NA          NA
## cobre                  NA         NA      NA         NA          NA
## manganeso              NA         NA      NA         NA          NA
## zinc                   NA         NA      NA         NA          NA
## boro                   NA         NA      NA         NA          NA
##               conductividad hierro_olsen cobre manganeso zinc boro
## ph_agua_suelo            NA           NA    NA        NA   NA   NA
## materia_org              NA           NA    NA        NA   NA   NA
## fosforo                  NA           NA    NA        NA   NA   NA
## azufre                   NA           NA    NA        NA   NA   NA
## acidez                   NA           NA    NA        NA   NA   NA
## aluminio                 NA           NA    NA        NA   NA   NA
## calcio                   NA           NA    NA        NA   NA   NA
## magnesio                 NA           NA    NA        NA   NA   NA
## potasio                  NA           NA    NA        NA   NA   NA
## sodio                    NA           NA    NA        NA   NA   NA
## cice                     NA           NA    NA        NA   NA   NA
## conductividad             1           NA    NA        NA   NA   NA
## hierro_olsen             NA            1    NA        NA   NA   NA
## cobre                    NA           NA     1        NA   NA   NA
## manganeso                NA           NA    NA         1   NA   NA
## zinc                     NA           NA    NA        NA    1   NA
## boro                     NA           NA    NA        NA   NA    1
  • Los valores “NA” en muchos de los pares de variables se debe a que hay datos nulos o ausentes en algunas variables, para ello podemos agregar la opción use = “pairwise.complete.obs”
datos_cacao %>% 
  select(where(is.numeric)) %>% 
  cor(method = "spearman", use = "pairwise.complete.obs")
##               ph_agua_suelo  materia_org     fosforo      azufre      acidez
## ph_agua_suelo  1.0000000000 -0.230444386  0.34841822 -0.30396544 -0.70038649
## materia_org   -0.2304443862  1.000000000  0.11982170  0.28251840  0.02051716
## fosforo        0.3484182193  0.119821704  1.00000000  0.01257518 -0.03923356
## azufre        -0.3039654363  0.282518396  0.01257518  1.00000000  0.04602610
## acidez        -0.7003864876  0.020517165 -0.03923356  0.04602610  1.00000000
## aluminio      -0.7028101188 -0.001360699 -0.04922720  0.04592097  0.99024772
## calcio         0.7760955503 -0.079288537  0.32325764 -0.25783325 -0.30948848
## magnesio       0.6861371704 -0.199304187  0.16566824 -0.25562563 -0.28482922
## potasio        0.3160512140  0.146996260  0.23259160  0.01560293  0.01402571
## sodio          0.2290031112 -0.278534242 -0.02606950 -0.01105171 -0.06589176
## cice           0.5204851874 -0.029755658  0.21246753 -0.16616407  0.36431191
## conductividad  0.1277957941  0.307475956  0.27383162  0.23435181 -0.06423920
## hierro_olsen  -0.7270326028  0.383686315 -0.09815333  0.24264939  0.52848692
## cobre          0.0738065788  0.278897482  0.17717447  0.09594644 -0.05463086
## manganeso     -0.0712990529  0.001240247 -0.11590729 -0.03808012 -0.19061161
## zinc           0.0004949256  0.370980642  0.34240702  0.11985823  0.07559154
## boro           0.1655318172  0.149915399  0.28221224  0.04799044  0.09369194
##                   aluminio      calcio   magnesio     potasio       sodio
## ph_agua_suelo -0.702810119  0.77609555  0.6861372  0.31605121  0.22900311
## materia_org   -0.001360699 -0.07928854 -0.1993042  0.14699626 -0.27853424
## fosforo       -0.049227201  0.32325764  0.1656682  0.23259160 -0.02606950
## azufre         0.045920969 -0.25783325 -0.2556256  0.01560293 -0.01105171
## acidez         0.990247717 -0.30948848 -0.2848292  0.01402571 -0.06589176
## aluminio       1.000000000 -0.33093050 -0.3050850 -0.01657663 -0.06039812
## calcio        -0.330930496  1.00000000  0.8575178  0.52366650  0.26187035
## magnesio      -0.305084992  0.85751778  1.0000000  0.57828386  0.32601216
## potasio       -0.016576632  0.52366650  0.5782839  1.00000000  0.24327762
## sodio         -0.060398124  0.26187035  0.3260122  0.24327762  1.00000000
## cice           0.337459665  0.85477166  0.7770685  0.57584694  0.27572723
## conductividad -0.088750176  0.40536385  0.3317315  0.45723377  0.13207982
## hierro_olsen   0.514250895 -0.47578082 -0.4795432 -0.16715145 -0.22102267
## cobre         -0.072075627  0.22572350  0.2326018  0.19131601 -0.11010921
## manganeso     -0.207066808  0.11003251  0.2215541  0.13887670  0.10920089
## zinc           0.038950674  0.25561072  0.1296750  0.26629655 -0.07927982
## boro           0.072145538  0.26313992  0.2183712  0.33076721  0.03779557
##                      cice conductividad hierro_olsen       cobre    manganeso
## ph_agua_suelo  0.52048519    0.12779579  -0.72703260  0.07380658 -0.071299053
## materia_org   -0.02975566    0.30747596   0.38368632  0.27889748  0.001240247
## fosforo        0.21246753    0.27383162  -0.09815333  0.17717447 -0.115907288
## azufre        -0.16616407    0.23435181   0.24264939  0.09594644 -0.038080117
## acidez         0.36431191   -0.06423920   0.52848692 -0.05463086 -0.190611612
## aluminio       0.33745967   -0.08875018   0.51425090 -0.07207563 -0.207066808
## calcio         0.85477166    0.40536385  -0.47578082  0.22572350  0.110032510
## magnesio       0.77706855    0.33173153  -0.47954324  0.23260183  0.221554070
## potasio        0.57584694    0.45723377  -0.16715145  0.19131601  0.138876702
## sodio          0.27572723    0.13207982  -0.22102267 -0.11010921  0.109200893
## cice           1.00000000    0.41916043  -0.26265305  0.24214105  0.047094988
## conductividad  0.41916043    1.00000000  -0.02295090  0.20858875  0.137934233
## hierro_olsen  -0.26265305   -0.02295090   1.00000000  0.08667420  0.029691996
## cobre          0.24214105    0.20858875   0.08667420  1.00000000  0.149757225
## manganeso      0.04709499    0.13793423   0.02969200  0.14975723  1.000000000
## zinc           0.28498753    0.35045423   0.26559727  0.32518473  0.184524411
## boro           0.28634744    0.30098347   0.01701895  0.04609187 -0.113590574
##                        zinc        boro
## ph_agua_suelo  0.0004949256  0.16553182
## materia_org    0.3709806420  0.14991540
## fosforo        0.3424070205  0.28221224
## azufre         0.1198582288  0.04799044
## acidez         0.0755915370  0.09369194
## aluminio       0.0389506735  0.07214554
## calcio         0.2556107192  0.26313992
## magnesio       0.1296750296  0.21837125
## potasio        0.2662965520  0.33076721
## sodio         -0.0792798162  0.03779557
## cice           0.2849875326  0.28634744
## conductividad  0.3504542308  0.30098347
## hierro_olsen   0.2655972747  0.01701895
## cobre          0.3251847263  0.04609187
## manganeso      0.1845244107 -0.11359057
## zinc           1.0000000000  0.21307846
## boro           0.2130784581  1.00000000

Correlograma

datos_cacao %>% 
  select(where(is.numeric)) %>% 
  cor(method = "spearman", use = "pairwise.complete.obs") %>% 
  corrplot()

  • Podemos mejorar el gráfico anterior de la siguiente manera:
datos_cacao %>% 
  select(where(is.numeric)) %>% 
  cor(method = "spearman", use = "pairwise.complete.obs") %>% 
  corrplot(diag = FALSE,
           type = "lower", #lower o upper
           tl.col = "black",
           tl.srt = 45,
           method = "pie")

  • Otra opción para graficar correlaciones:
correlaciones <- datos_cacao %>% 
  select(where(is.numeric)) %>% 
  correlate(method = "spearman")

network_plot(correlaciones)

Series temporales

Datos TRM

df_trm <- read_csv("datos/TRM_20231022.csv")
df_trm %>% head()
  • ¿Cuál es el tipo de datos?
df_trm %>% 
  glimpse()
## Rows: 7,664
## Columns: 3
## $ VALOR         <dbl> 4761.64, 4781.28, 4802.48, 4836.24, 4815.99, 4818.32, 47…
## $ VIGENCIADESDE <chr> "22/12/2022", "20/12/2022", "17/12/2022", "13/12/2022", …
## $ VIGENCIAHASTA <chr> "22/12/2022", "20/12/2022", "19/12/2022", "13/12/2022", …

Procesando fechas

  • Podemos editar las fechas para que tengan el formato correcto (Date):
df_trm2 <- df_trm %>% 
  mutate(VIGENCIADESDE = dmy(VIGENCIADESDE),
         VIGENCIAHASTA = dmy(VIGENCIAHASTA))

df_trm2 %>% head()
  • Tener las fechas en el formato correcto habilita otras funciones para extraer información como los meses, días del mes, semana, trimestre, etc.
df_trm3 <-
  df_trm2 %>%
  mutate(
    year_es = year(VIGENCIADESDE),
    mes = month(VIGENCIADESDE, label = TRUE, abbr = FALSE),
    dia_mes = day(VIGENCIADESDE),
    dia_semana = weekdays(VIGENCIADESDE),
    semana = week(VIGENCIADESDE),
    semestre = semester(VIGENCIADESDE),
    trimestre = quarter(VIGENCIADESDE)
  )

df_trm3 %>% head()

Gráfico de series temporales

df_trm3 %>% 
  ggplot(aes(x = VIGENCIADESDE, y = VALOR)) +
  geom_line() +
  labs(x = "Fecha", y = "TRM")

  • Podemos agregar una línea de tendencia a la serie temporal:
df_trm3 %>% 
  ggplot(aes(x = VIGENCIADESDE, y = VALOR)) +
  geom_line() +
  geom_smooth() +
  labs(x = "Fecha", y = "TRM")

  • Podemos graficar patrones mensuales con un boxplot:
df_trm3 %>% 
  ggplot(aes(x = mes, y = VALOR)) +
  geom_boxplot()

  • Podríamos tratar de intuir patrones estacionales por año y mes:
df_trm3 %>% 
  group_by(year_es, mes) %>% 
  reframe(promedio = mean(VALOR, na.rm = TRUE)) %>% 
  mutate(year_es = as.factor(year_es)) %>% 
  ggplot(aes(x = mes, y = promedio, color = year_es)) +
  geom_line(aes(group = year_es))

Mapas

Colombia - departamentos

colombia_deptos <- st_read("datos/mapas-colombia/departamentos/MGN_DPTO_POLITICO.shp")
## Reading layer `MGN_DPTO_POLITICO' from data source 
##   `D:\Otros\UdeA\2023-02\estadistica\datos\mapas-colombia\departamentos\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 %>% 
  ggplot() +
  geom_sf()

Colombia - municipios

colombia_mpios <- st_read("datos/mapas-colombia/municipios/MGN_MPIO_POLITICO.shp")
## Reading layer `MGN_MPIO_POLITICO' from data source 
##   `D:\Otros\UdeA\2023-02\estadistica\datos\mapas-colombia\municipios\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 %>% 
  ggplot() +
  geom_sf()

Cacao por departamento

  • Verificamos que los nombres de los departamentos en ambas tablas sean iguales:
unique(datos_cacao$depto)
##  [1] "Meta"               "Boyacá"             "Antioquia"         
##  [4] "Santander"          "Cesar"              "Cundinamarca"      
##  [7] "Valle Del Cauca"    "Casanare"           "Arauca"            
## [10] "Huila"              "Norte De Santander" "Nariño"            
## [13] "Cauca"              "Guaviare"           "Tolima"            
## [16] "Caldas"             "Caquetá"            "Putumayo"          
## [19] "Bolívar"            "Córdoba"            "La Guajira"        
## [22] "Sucre"              "Magdalena"          "Risaralda"         
## [25] "Chocó"              "Bogotá, D.c."       "Atlántico"
unique(colombia_deptos$DPTO_CNMBR)
##  [1] "ANTIOQUIA"                                               
##  [2] "ATLÁNTICO"                                               
##  [3] "BOGOTÁ, D.C."                                            
##  [4] "BOLÍVAR"                                                 
##  [5] "BOYACÁ"                                                  
##  [6] "CALDAS"                                                  
##  [7] "CAQUETÁ"                                                 
##  [8] "CAUCA"                                                   
##  [9] "CESAR"                                                   
## [10] "CÓRDOBA"                                                 
## [11] "CUNDINAMARCA"                                            
## [12] "CHOCÓ"                                                   
## [13] "HUILA"                                                   
## [14] "LA GUAJIRA"                                              
## [15] "MAGDALENA"                                               
## [16] "META"                                                    
## [17] "NARIÑO"                                                  
## [18] "NORTE DE SANTANDER"                                      
## [19] "QUINDIO"                                                 
## [20] "RISARALDA"                                               
## [21] "SANTANDER"                                               
## [22] "SUCRE"                                                   
## [23] "TOLIMA"                                                  
## [24] "VALLE DEL CAUCA"                                         
## [25] "ARAUCA"                                                  
## [26] "CASANARE"                                                
## [27] "PUTUMAYO"                                                
## [28] "ARCHIPIÉLAGO DE SAN ANDRÉS, PROVIDENCIA Y SANTA CATALINA"
## [29] "AMAZONAS"                                                
## [30] "GUAINÍA"                                                 
## [31] "GUAVIARE"                                                
## [32] "VAUPÉS"                                                  
## [33] "VICHADA"
  • Editamos los nombres de los departamentos para que coincidan:
colombia_deptos2 <-
  colombia_deptos %>%
  mutate(DPTO_CNMBR = str_to_title(DPTO_CNMBR))

colombia_deptos2$DPTO_CNMBR %>% unique()
##  [1] "Antioquia"                                               
##  [2] "Atlántico"                                               
##  [3] "Bogotá, D.c."                                            
##  [4] "Bolívar"                                                 
##  [5] "Boyacá"                                                  
##  [6] "Caldas"                                                  
##  [7] "Caquetá"                                                 
##  [8] "Cauca"                                                   
##  [9] "Cesar"                                                   
## [10] "Córdoba"                                                 
## [11] "Cundinamarca"                                            
## [12] "Chocó"                                                   
## [13] "Huila"                                                   
## [14] "La Guajira"                                              
## [15] "Magdalena"                                               
## [16] "Meta"                                                    
## [17] "Nariño"                                                  
## [18] "Norte De Santander"                                      
## [19] "Quindio"                                                 
## [20] "Risaralda"                                               
## [21] "Santander"                                               
## [22] "Sucre"                                                   
## [23] "Tolima"                                                  
## [24] "Valle Del Cauca"                                         
## [25] "Arauca"                                                  
## [26] "Casanare"                                                
## [27] "Putumayo"                                                
## [28] "Archipiélago De San Andrés, Providencia Y Santa Catalina"
## [29] "Amazonas"                                                
## [30] "Guainía"                                                 
## [31] "Guaviare"                                                
## [32] "Vaupés"                                                  
## [33] "Vichada"
  • ¿Coinciden los nombres en ambas tablas?
table(unique(datos_cacao$depto) %in% unique(colombia_deptos2$DPTO_CNMBR))
## 
## TRUE 
##   27
  • Primero debemos resumir a nivel de departamento los datos del cacao. En este caso calcularemos el promedio del pH para cada departamento:
resumen_depto <-
  datos_cacao %>% 
  group_by(depto) %>% 
  reframe(promedio_ph = mean(ph_agua_suelo, na.rm = TRUE))

resumen_depto
  • Unimos la tabla resumida por departamento con la tabla que tiene la información del mapa por departamento:
depto_cacao <-
  colombia_deptos2 %>%
  left_join(resumen_depto, by = c("DPTO_CNMBR" = "depto"))
  • Finalmente graficamos el promedio de pH por departamento:
depto_cacao %>% 
  ggplot(aes(fill = promedio_ph)) +
  geom_sf() +
  scale_fill_viridis_c() +
  annotation_north_arrow(location = "tl") +
  labs(fill = "pH")