Bibliotecas y configuración

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

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

suelos %>% head()

Funciones para manejo de datos

  • count
  • filter
  • select
  • mutate
  • group_by + reframe (summarise)

Funciones para estadística descriptiva

Función Descripción Tipo de variable
mean() Calcular promedio cuantitativa
weighted.mean() Calcular promedio ponderado Cunatitativa
median() Calcular mediana Cuantitativa
sd() Calcular desviación estándar cuantitativa
var() Calcular la varianza Cuantitativa
range() Calcular el rango Cuantitativa
IQR() Calcular rango intercuartílico Cuantitativa
quantile() Calcular cuartiles, deciles y percentiles Cuantitativa
min() Valor mínimo Cuantitativa
max() Valor máximo Cuantitativa

Conteos por cultivo

suelos %>% 
  count(cultivo, sort = TRUE)

Cultivo de Cacao

Datos Cacao

datos_cacao <- suelos %>% 
  filter(cultivo == "Cacao")

datos_cacao %>% head()

Conteos por departamento

datos_cacao %>% 
  count(depto, sort = TRUE) %>% 
  mutate(porcentaje = (n / sum(n)) * 100)

Conteos por cultivo y estado del cultivo

datos_cacao %>% 
  count(estado) %>% 
  mutate(porcentaje = (n / sum(n)) * 100)

Conteos por cultivo y topografía

datos_cacao %>% 
  count(topografia, sort = TRUE) %>% 
  mutate(porcentaje = (n / sum(n)) * 100)

Promedio de Aluminio

mean(datos_cacao$aluminio, na.rm = TRUE)
## [1] 2.043409

Resumen descriptivo general

  • ¿Hay algún error en el resumen descriptivo?
datos_cacao %>% 
  reframe(
    promedio_ph = mean(ph_agua_suelo, na.rm = TRUE),
    mediana_ph = median(ph_agua_suelo, na.rm = TRUE),
    desv_ph = sd(ph_agua_suelo, na.rm = TRUE),
    minimo_ph = min(ph_agua_suelo, na.rm = TRUE),
    maximo_ph = max(ph_agua_suelo, na.rm = TRUE)
  ) %>% 
  mutate(coef_var_ph = (desv_ph / promedio_ph) * 100)
  • Obtengamos los cuartiles para el pH:
datos_cacao %>% 
  pull(ph_agua_suelo) %>% 
  quantile()
##     0%    25%    50%    75%   100% 
##   3.69   4.79   5.31   6.09 380.00
# Este código es equivalente al anterior
#quantile(datos_cacao$ph_agua_suelo)
  • Obtengamos los deciles del pH:
datos_cacao %>% 
  pull(ph_agua_suelo) %>% 
  quantile(probs = seq(0, 1, 0.1))
##     0%    10%    20%    30%    40%    50%    60%    70%    80%    90%   100% 
##   3.69   4.47   4.69   4.88   5.08   5.31   5.57   5.92   6.26   6.72 380.00
  • Obtengamos los percentiles del pH:
datos_cacao %>% 
  pull(ph_agua_suelo) %>% 
  quantile(probs = seq(0, 1, 0.01))
##       0%       1%       2%       3%       4%       5%       6%       7% 
##   3.6900   4.0500   4.1600   4.2300   4.2900   4.3200   4.3500   4.3800 
##       8%       9%      10%      11%      12%      13%      14%      15% 
##   4.4100   4.4400   4.4700   4.4900   4.5200   4.5400   4.5600   4.5900 
##      16%      17%      18%      19%      20%      21%      22%      23% 
##   4.6100   4.6300   4.6500   4.6700   4.6900   4.7055   4.7200   4.7500 
##      24%      25%      26%      27%      28%      29%      30%      31% 
##   4.7700   4.7900   4.8000   4.8200   4.8400   4.8600   4.8800   4.9000 
##      32%      33%      34%      35%      36%      37%      38%      39% 
##   4.9100   4.9300   4.9500   4.9700   4.9900   5.0100   5.0300   5.0600 
##      40%      41%      42%      43%      44%      45%      46%      47% 
##   5.0800   5.1000   5.1200   5.1500   5.1700   5.1900   5.2200   5.2400 
##      48%      49%      50%      51%      52%      53%      54%      55% 
##   5.2600   5.2900   5.3100   5.3400   5.3600   5.3900   5.4100   5.4400 
##      56%      57%      58%      59%      60%      61%      62%      63% 
##   5.4600   5.4900   5.5190   5.5400   5.5700   5.6100   5.6400   5.6800 
##      64%      65%      66%      67%      68%      69%      70%      71% 
##   5.7100   5.7400   5.7730   5.8100   5.8400   5.8900   5.9200   5.9600 
##      72%      73%      74%      75%      76%      77%      78%      79% 
##   5.9800   6.0200   6.0500   6.0900   6.1200   6.1500   6.1890   6.2200 
##      80%      81%      82%      83%      84%      85%      86%      87% 
##   6.2600   6.3100   6.3500   6.3900   6.4400   6.4800   6.5200   6.5800 
##      88%      89%      90%      91%      92%      93%      94%      95% 
##   6.6300   6.6600   6.7200   6.7805   6.8460   6.9100   7.0000   7.0900 
##      96%      97%      98%      99%     100% 
##   7.2180   7.3500   7.4800   7.6490 380.0000

Filtrando datos

  • ¿Cuántos registros superan el valor de 10 de pH?
datos_cacao %>% 
  filter(ph_agua_suelo > 10)
  • Filtramos datos con pH menor o igual a 10:
datos_cacao2 <- datos_cacao %>% 
  filter(ph_agua_suelo <= 10)

datos_cacao2 %>% head()
  • Calculamos nuevamente el resumen descriptivo:
datos_cacao2 %>% 
  reframe(
    promedio_ph = mean(ph_agua_suelo, na.rm = TRUE),
    mediana_ph = median(ph_agua_suelo, na.rm = TRUE),
    desv_ph = sd(ph_agua_suelo, na.rm = TRUE),
    minimo_ph = min(ph_agua_suelo, na.rm = TRUE),
    maximo_ph = max(ph_agua_suelo, na.rm = TRUE)
  ) %>% 
  mutate(coef_var_ph = (desv_ph / promedio_ph) * 100)

Moda pH

moda <- function(x) {
  ux = unique(x)
  tab = tabulate(match(x, ux))
  ux[tab == max(tab)]
}

moda(datos_cacao2$ph_agua_suelo)
## [1] 4.96 4.86 4.84

Resumen por departamento

  • Un resumen descriptivo compuesto de varias métricas estadísticas:
datos_cacao2 %>% 
  group_by(depto) %>% 
  reframe(
    promedio_ph = mean(ph_agua_suelo, na.rm = TRUE),
    mediana_ph = median(ph_agua_suelo, na.rm = TRUE),
    desv_ph = sd(ph_agua_suelo, na.rm = TRUE),
    minimo_ph = min(ph_agua_suelo, na.rm = TRUE),
    maximo_ph = max(ph_agua_suelo, na.rm = TRUE),
    rango_interq = IQR(ph_agua_suelo, na.rm = TRUE),
    percentil5 = quantile(ph_agua_suelo, probs = 0.05, na.rm = TRUE),
    percentil95 = quantile(ph_agua_suelo, probs = 0.95, na.rm = TRUE),
    coef_asimetria = skewness(ph_agua_suelo, na.rm = TRUE),
    coef_curtosis = kurtosis(ph_agua_suelo, na.rm = TRUE)
  ) %>% 
  mutate(coef_var_ph = (desv_ph / promedio_ph) * 100) %>% 
  arrange(desc(promedio_ph))

Resumen por departamento y estado del cultivo

datos_cacao2 %>% 
  group_by(depto, estado) %>% 
  reframe(
    promedio_ph = mean(ph_agua_suelo, na.rm = TRUE),
    mediana_ph = median(ph_agua_suelo, na.rm = TRUE),
    desv_ph = sd(ph_agua_suelo, na.rm = TRUE),
    minimo_ph = min(ph_agua_suelo, na.rm = TRUE),
    maximo_ph = max(ph_agua_suelo, na.rm = TRUE)
  ) %>% 
  mutate(coef_var_ph = (desv_ph / promedio_ph) * 100)

Promedio ponderado

animales <- c(20, 150, 23, 24, 2)
peso <- c(110.1, 120.3, 112.5, 110.8, 95.6)

# Promedio aritmético
mean(peso)
## [1] 109.86
# Promedio pondedaro
weighted.mean(x = peso, w = animales)
## [1] 117.2826

Visualizaciones

Gráficos con ggplot2

Capas ggplot2

Estructura básica

data %>% 
  ggplot(mapping = aes(x = ..., y = ...)) +
  geom_...

Cantidades (conteos)

  • Podemos graficar el conteo de registros por departamento:
datos_cacao2 %>% 
  count(depto) %>% 
  ggplot(aes(x = depto, y = n)) +
  geom_col()

  • Podemos mejorar el gráfico anterior girando las coordenadas y reordenando los departamentos por el número de datos:
datos_cacao2 %>% 
  count(depto) %>% 
  ggplot(aes(x = reorder(depto, n), y = n)) +
  geom_col() +
  coord_flip()

  • Ejemplo general de edición de color, etiquetas y tema del gráfico:
datos_cacao2 %>% 
  count(depto) %>% 
  ggplot(aes(x = reorder(depto, n), y = n)) +
  geom_col(color = "#46a079", fill = "orchid") +
  coord_flip() +
  labs(title = "Título",
       subtitle = "Subtítulo",
       x = "Eje x",
       y = "Eje y",
       caption = "Nota en la parte inferior del gráfico")

  • Agreguemos etiquetas a los ejes, cambiemos el color de las barras y editemos el tema del gráfico:
grafico1 <-
  datos_cacao2 %>%
  count(depto) %>%
  ggplot(aes(x = reorder(depto, n), y = n)) +
  geom_col(color = "black", fill = "dodgerblue2") +
  coord_flip() +
  labs(
    x = "Departamento",
    y = "Total (n)",
    title = "Distribución de registros por departamento",
    subtitle = "Cultivo de Cacao"
  ) +
  theme_bw()

grafico1

  • Visualizando de forma interactiva el gráfico anterior:
ggplotly(grafico1)
  • Otra forma de representar visualmente los conteos o frecuencias absolutas, es a través del gráfico llamado treemap:
conteos_depto <-
  datos_cacao2 %>%
  count(depto)

library(treemap)
treemap(conteos_depto,
        index = "depto",
        vSize = "n",
        title = "")

  • Al gráfico anterior lo podemos agregar más variables (además del departamento):
conteos_depto_estado <-
  datos_cacao2 %>%
  count(depto, estado)

treemap(conteos_depto_estado,
        index = c("depto", "estado"),
        vSize = "n",
        title = "")

Cantidades (promedio)

  • Promedio de pH por departamento:
datos_cacao2 %>% 
  group_by(depto) %>% 
  reframe(promedio_ph = mean(ph_agua_suelo, na.rm = TRUE)) %>% 
  ggplot(aes(x = reorder(depto, promedio_ph), y = promedio_ph)) +
  geom_col() +
  coord_flip() +
  labs(x = "Departamento", y = "pH")

  • Promedio de pH por departamento y estado del cultivo:
grafico2 <-
  datos_cacao2 %>%
  group_by(depto, estado) %>%
  reframe(promedio_ph = mean(ph_agua_suelo, na.rm = TRUE)) %>%
  ggplot(aes(
    x = reorder(depto, promedio_ph),
    y = promedio_ph,
    fill = estado
  )) +
  geom_col(position = "dodge") +
  labs(x = "Departamento", y = "pH", fill = "Estado:") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top")

grafico2

  • El mismo gráfico anterior agregando interactividad:
ggplotly(grafico2)
  • Otra forma de representar visualmente el mismo resultado anterior:
datos_cacao2 %>%
  group_by(depto, estado) %>%
  reframe(promedio_ph = mean(ph_agua_suelo, na.rm = TRUE)) %>%
  ggplot(aes(x = depto,
             y = estado,
             color = promedio_ph)) +
  geom_point(size = 5) +
  scale_color_viridis_c() +
  labs(x = "Departamento", y = "Estado", color = "pH") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.key.width = unit(2, "cm"))

Proporciones

  • Primera opción:
datos_cacao2 %>%
  count(depto, estado) %>%
  ggplot(aes(x = depto, y = n, fill = estado)) +
  geom_col()

  • Una mejor opción:
datos_cacao2 %>%
  count(depto, estado) %>%
  ggplot(aes(x = depto, y = n, fill = estado)) +
  geom_col(position = "fill") +
  coord_flip() +
  labs(x = "Departamento",
       y = "Proporción",
       fill = "Estado")

Distribuciones (histogramas)

  • pH:
datos_cacao2 %>% 
  ggplot(aes(x = ph_agua_suelo)) +
  geom_histogram(color = "black", fill = "blue")

  • pH por estado del cultivo:
datos_cacao2 %>% 
  ggplot(aes(x = ph_agua_suelo, fill = estado)) +
  geom_histogram(color = "white")

  • ¿Podemos separar los gráficos?
datos_cacao2 %>% 
  ggplot(aes(x = ph_agua_suelo, fill = estado)) +
  facet_wrap(~estado, scales = "free") +
  geom_histogram(color = "white") +
  theme(legend.position = "none")

Distribuciones (densidades)

  • pH:
datos_cacao2 %>% 
  ggplot(aes(x = ph_agua_suelo)) +
  geom_density(color = "black")

  • pH por estado del cultivo:
datos_cacao2 %>% 
  ggplot(aes(x = ph_agua_suelo, fill = estado)) +
  geom_density(color = "white", alpha = 0.5)

  • ¿Podemos separar los gráficos?
datos_cacao2 %>% 
  ggplot(aes(x = ph_agua_suelo, fill = estado)) +
  facet_wrap(~estado, scales = "free") +
  geom_density()

  • ¿Podemos incluir todas las variables numéricas en un sólo gráfico?
datos_cacao2 %>% 
  select(where(is.numeric)) %>% 
  pivot_longer(cols = everything()) %>% 
  ggplot(aes(x = value)) +
  facet_wrap(~name, scales = "free") +
  geom_density() +
  scale_x_log10()

Distribuciones (boxplots)

  • pH por departamento:
datos_cacao2 %>% 
  ggplot(aes(x = depto, y = ph_agua_suelo)) +
  geom_boxplot() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

  • pH por departamento y estado de cultivo:
datos_cacao2 %>% 
  ggplot(aes(x = depto, y = ph_agua_suelo, fill = estado)) +
  geom_boxplot() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top")

Incertidumbre

datos_cacao2 %>% 
  group_by(depto) %>% 
  reframe(promedio_ph = mean(ph_agua_suelo, na.rm = TRUE),
          desv_ph = sd(ph_agua_suelo, na.rm = TRUE)) %>% 
  ggplot(aes(x = depto, y = promedio_ph,
             ymin = promedio_ph - desv_ph,
             ymax = promedio_ph + desv_ph)) +
  geom_point() +
  geom_errorbar(width = 0.25)  +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

  • El mismo gráfico anterior puede ser obtenido con la función geom_pointrange:
datos_cacao2 %>% 
  group_by(depto) %>% 
  reframe(promedio_ph = mean(ph_agua_suelo, na.rm = TRUE),
          desv_ph = sd(ph_agua_suelo, na.rm = TRUE)) %>% 
  ggplot(aes(x = depto, y = promedio_ph,
             ymin = promedio_ph - desv_ph,
             ymax = promedio_ph + desv_ph)) +
  geom_pointrange() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

  • Podemos usar barras para representar el promedio y agregar las líneas que representen la desviación estándar:
datos_cacao2 %>% 
  group_by(depto) %>% 
  reframe(promedio_ph = mean(ph_agua_suelo, na.rm = TRUE),
          desv_ph = sd(ph_agua_suelo, na.rm = TRUE)) %>% 
  ggplot(aes(x = depto, y = promedio_ph,
             ymin = promedio_ph - desv_ph,
             ymax = promedio_ph + desv_ph)) +
  geom_col(fill = "gray50") +
  geom_point() +
  geom_errorbar(width = 0.2) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

  • Podemos agregar las métricas generales (promedio, mediana y desviación) al gráfico de incertidumbre:
promedio_ph_general <- mean(datos_cacao2$ph_agua_suelo, na.rm = TRUE)
mediana_ph_general <- median(datos_cacao2$ph_agua_suelo, na.rm = TRUE)
desviacion_ph_general <- sd(datos_cacao2$ph_agua_suelo, na.rm = TRUE)

datos_cacao2 %>% 
  group_by(depto) %>% 
  reframe(promedio_ph = mean(ph_agua_suelo, na.rm = TRUE),
          desv_ph = sd(ph_agua_suelo, na.rm = TRUE)) %>% 
  ggplot(aes(x = depto, y = promedio_ph,
             ymin = promedio_ph - desv_ph,
             ymax = promedio_ph + desv_ph)) +
  geom_point() +
  geom_errorbar(width = 0.25)  +
  geom_hline(yintercept = promedio_ph_general, linetype = 2, color = "red") +
  geom_hline(yintercept = mediana_ph_general, linetype = 2, color = "forestgreen") +
  geom_hline(yintercept = promedio_ph_general - desviacion_ph_general,
             linetype = 2, color = "blue") +
  geom_hline(yintercept = promedio_ph_general + desviacion_ph_general,
             linetype = 2, color = "blue") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = "Departamento",
       y = "pH",
       caption = "Línea roja: promedio\nLínea verde: mediana\nLíneas azules: más 1 y menos 1 DE")