In bitsandbricks/dataviz: Lecciones de visualización de datos

library(learnr)
library(dplyr)
library(ggplot2)
library(ggmap)
library(lubridate)

delitos <- dataviz::delitos
delitos$fecha <- ymd(delitos$fecha)

set.seed("99")
muestra_de_fechas <- delitos %>% 
    sample_n(5) %>% 
    pull(fecha)

bbox <- make_bbox(delitos$longitud, delitos$latitud)
CABA <- get_stamenmap(bbox = bbox, maptype = "toner-lite", zoom = 12)

knitr::opts_chunk$set(message = FALSE, warning = FALSE, comment = "")

Trabajando con fechas

Se hace cada vez más fácil registrar eventos con gran precisión temporal además de espacial y temporal; es decir, datos que incluyen fecha y hora, o coordenadas que los ubican en un sitio determinado. Para entender datasets con datos en gran volumen que poseen atributos de posición y tiempo, es útil visualizar el ritmo en el que ocurren (diario, mensual, anual, etc) y la forma en la que se distribuyen en el espacio.

Para practicar, trabajaremos con un dataset publicado por la Ciudad Autónoma de Buenos Aires, con delitos registrados durante el 2020. Los datos fueron publicados por en el Mapa del delito de la Ciudad (https://mapa.seguridadciudad.gob.ar/).

knitr::include_graphics('https://live.staticflickr.com/3410/3346781444_48266ac902_b.jpg')

subtipo: "Siniestro Vial"

Las primeras filas del dataset lucen asi:

head(delitos, 50)

Teniendo en cuenta que el dataset es del 2020, el Primer Año de la Pandemia, nos podríamos preguntar si se ve algún descenso brusco en robos, hurtos, homicidios o siniestros viales con el inicio de la cuarentena. Para responder eso, vamos a tener que lidiar con datos de tipo "fecha". Esto podría ser bastante engorroso, pero por suerte podemos usar una herramienta que simplifica las cosas.

La fecha como clase de variable

La fecha es un tipo de dato que puede ser expresado de muchas maneras, dependiendo de que nos interese tener en cuenta: el día de la semana al que corresponde, el mes, el año, etc. El paquete lubridate hace fácil tomar una variable que contiene fechas en cualquier formato (por ejemplo "20/07/2020") para extraer el atributo relacionado que deseemos (como su día, "lunes", o su mes, "julio").

Para empezar, convertimos el campo "fecha" al tipo de dato especializado, que se llama... fecha (date). Aquí tenemos que prestar atención al formato en que aparecen los datos de la columna, en general algo como "2018-07-21" (mes, día y año) o "2018-07-21 12:14:24" (mes, día, año y hora, minutos, segundos). Con nuestros datos se da el primer caso, por lo cual la función para convertir ese campo en fecha es ymd(); para el segundo caso, seria ymd_hms()

library(lubridate)

delitos$fecha <- ymd(delitos$fecha)

Repasemos algunas de los nuevos trucos que podemos hacer con el tiempo. Tomemos cinco fechas elegidas al azar:

muestra_de_fechas

Mediante las funciones disponibles en lubridate, podemos extraer:

El día de la semana al que corresponde cada fecha:

wday(muestra_de_fechas)

wday(muestra_de_fechas, label = TRUE)

El mes:

month(muestra_de_fechas)

month(muestra_de_fechas, label = TRUE)

El año:

year(muestra_de_fechas)

Y varias opciones más, que se pueden repasar en https://cran.r-project.org/web/packages/lubridate/vignettes/lubridate.html

Con lo visto hasta aquí, tenemos suficiente para mostrar patrones temporales en los datos.

Empecemos por un gráfico de barras (geom_bar())con la cantidad de eventos registrados por mes. Para que aparezca un conteo por mes del año, asignaremos al eje de las $x$ el mes al que corresponde cada valor de la columna "fecha", o sea month(fecha, label = TRUE):

ggplot(radios) +
  geom___(___)

ggplot(delitos) + 
    geom_bar(aes(x = month(___)))

ggplot(delitos) + 
    geom_bar(aes(x = month(fecha, label = TRUE)))

(Usamos el parámetro label = TRUE para obtener el nombre del mes en lugar de su número -"dic" en lugar de 12. Prueben realizar el gráfico sin ese parámetro en la llamada a month() para ver que pasa.)

En el gráfico, se ve una reducción drástica a partir de la cuarentena decretada en abril. Y también un incremento gradual y sostenido a partir de allí, que de todos modos no llega a los niveles pre-pandemia de enero y febrero.

Para ver la composición interna de los conteos mensuales, cuantos casos corresponden a cada categoría, podemos hacer un gráfico de "barras apiladas" como vimos en la clase 3. La sintaxis es igual que antes, pero esta vez asignamos la variable tipo"_ al color de relleno de las barras determinado por el parámetro fill.

ggplot(delitos) + 
    geom_bar(aes(x = month(fecha, label = TRUE), ___))

ggplot(delitos) + 
    geom_bar(aes(x = month(fecha, label = TRUE), fill = ___))

ggplot(delitos) + 
    geom_bar(aes(x = month(fecha, label = TRUE), fill = tipo))

Las barras apiladas son prolijas, pero pueden hacer difícil evaluar la evolución de categorías individuales. Recordemos que para mostrar los subconjuntos en barras independientes, una al lado de la otra, podemos usar el parámetro position = "dodge". Inténtenlo agregando el parámetro al ejercicio anterior para ver como queda.

Ahora comparemos la cantidad de eventos registrados, por tipo, para cada día de la semana. Basta con usar la función que extrae el día de la semana, wday(), en lugar de month(). Lo demás es idéntico, incluyendo el uso de label = TRUE para que obtener el nombre del día -"lun"-, en lugar de su posición en la semana -"2", porque para lubridate las semanas empiezan el domingo-.

ggplot(delitos) + 
    geom_bar(aes(x = ___(fecha, label = TRUE), fill = tipo))

ggplot(delitos) + 
    geom_bar(aes(x = wday(fecha, label = TRUE), fill = tipo))

Como era de esperar, durante los fines de semana se observa una menor cantidad de eventos, aunque quizás no en la categoría Homicidio, que es difícil de discernir por su relativa escasez. Para solucionar el problema podríamos filtrar los datos antes de visualizarlos, como hemos hecho antes con la función filter(), mostrando sólo la categoría de interés.

Otra opción, que probaremos ahora, es una visualización en "facetas". Se trata de generar una visualización que muestra distintos aspectos de los datos en paneles separados. Con un ejemplo vamos a dejarlo mas claro. Para realizar un gráfico en facetas, basta con sumar una línea con la función facet_wrap(vars(x, y, ...)), dónde "x", "y", etc son los nombres de las variables cuyas categorías recibirán paneles distintos. Intentémoslo con el último gráfico, agregando "tipo" como variable a facetar:

ggplot(delitos) + 
    geom_bar(aes(x = wday(fecha, label = TRUE), fill = tipo)) +
    facet_wrap(vars(___))

ggplot(delitos) + 
    geom_bar(aes(x = wday(fecha, label = TRUE), fill = tipo)) +
    facet_wrap(vars(tipo))

Obtuvimos barras separadas, pero los homicidios siguen difíciles de distinguir. Se debe a que por defecto facet_wrap() mantiene a escala todos los paneles, de manera que se puedan comparar cantidades de forma directa. La desventaja es que se pierde legibilidad de categorías con valores ínfimos en comparación con otras. Por eso disponemos del parámetro scales, que permite graficar los datos con escala libre, vía scales = "free". Probemos:

ggplot(delitos) + 
    geom_bar(aes(x = wday(fecha, label = TRUE), fill = tipo)) +
    facet_wrap(vars(___), scales = ___)

ggplot(delitos) + 
    geom_bar(aes(x = wday(fecha, label = TRUE), fill = tipo)) +
    facet_wrap(vars(tipo), scales = "free")

Ahora queda más claro que los homicidios siguen un patrón diario distinto al de las demás categorías. También que, en términos relativos, las lesiones son las que mas se reducen durante los fines de semana.

También podemos evaluar el ritmo según la hora del día. ¿Cómo se haría con nuestro dataset?

Mirando al espacio

Pasemos ahora al análisis espacial de nuestros datos. Para facilitar la visualización volveremos a usar el paquete ggmap, que aporta funciones que facilitan la creación de mapas.

library(ggmap)

Obteniendo un mapa base

Como vimos en la clase previa, para obtener un mapa de fondo o "mapa base" necesitamos obtener una bounding box de nuestros datos, que luego pasamos a get_stamenmap().

La diferencia con nuestra experiencia anterior es que antes trabajamos con dataframes espaciales, que contenían geometrías georreferencias: polígonos, líneas, puntos. Pero este dataframe no tiene nada de eso. La buena noticia es que igual podemos hacer mapas, porque contiene dos columnas clave para ello: las de "longitud" y "latitud", que permiten ubicar puntos sobre la faz de la tierra. Y de puntos se trata, porque en general sólo podremos trabajar con columnas de coordenadas cuando los datos espaciales representen puntos, pero no líneas y polígonos. Estas geoemtrías de mayor complejidad se manejan de forma cómodo usando dataframes espaciales como los que exploramos en la clase anterior.

Así como contamos con st_bbox() para obtener la bounding box de dataframes espaciales, cuando tenemos las coordenadas en un par de columnas recurrimos a make_bbox(),:

bbox <- make_bbox(delitos$longitud, delitos$latitud)

bbox

En base a la "bounding box" solicitamos nuestro mapa base:

CABA <- get_stamenmap(bbox = bbox, maptype = "toner-lite", zoom = 12)

Para verlo:

ggmap(CABA)

De coordenadas al mapa

De aquí en más podemos suporponer nuestros datos en distintas capas, con la misma sintaxis que conocemos de ggplot. Para mapear las ubicaciones de los delitos en el dataset, usamos geom_point() y los campos de longitud y latitud para los ejes $x$ e $y$:

ggmap(CABA) +
    geom____(data = delitos, aes(___))

ggmap(CABA) +
    geom_point(data = delitos, aes(x = ___, y = ___))

ggmap(CABA) +
    geom_point(data = delitos, aes(x = longitud, y = latitud))

Aquí nos topamos con un problema habitual al trabajar con grandes volúmenes de datos. Hay tantos puntos proyectados sobre el mapa, que se hace imposible interpretar dónde existen más o menos. Hacemos algunos ajustes: - un color que resalte más contra el mapa base, y que no se confunda con él - un tamaño de punto más pequeño -y aplicación de una ligera transparencia

Todo ello vía los atributos "color", "size" y "alpha". ¿Cuál es el valor ideal para cada uno? En general, no queda otra que recurrir a la prueba y error para encontrar la receta justa. Probemos con color = "orange", size = 0.1 y alpha = 0.1:

ggmap(CABA) +
    geom_point(data = delitos, aes(x = longitud, y = latitud),
               ___)

ggmap(CABA) +
    geom_point(data = delitos, aes(x = longitud, y = latitud),
               color = ___, size = ___, alpha = ___)

ggmap(CABA) +
    geom_point(data = delitos, aes(x = longitud, y = latitud),
               color = "orange", size = 0.1, alpha = 0.1)

Ahora si aparecen ciertos patrones, por ejemplo la mayor frecuencia de casos de casos cerca de las principales de circulación de la ciudad. Aún así, se hace difícil identificar de un golpe de vista las "zonas calientes", los puntos de máxima concentración.

Mapas de densidad

Una solución práctica para el problema de la cantidad de puntos es una técnica llamada "binning": dividir el espacio en una grilla de celdas, contar cuantos puntos caen dentro de cada una, y visualizar las cantidades agregadas. En el mundo ggplot esto se lleva a cabo con geom_bind2d().

ggmap(CABA) +
    ___(data = delitos, aes(x = longitud, y = latitud))

ggmap(CABA) +
    geom_bin2d(data = delitos, aes(x = longitud, y = latitud))

Ahora si, resaltan las áreas de mayor concentración de incidentes. Se puede mejorar un poco el gráfico usando:

una mayor cantidad de celdas para aumentar la resolución
una escala de colores diseñada para ayudar a detectar diferencias por tonalidad, como Viridis.

la cantidad de celdas se define con el parámetro "bins", por ejemplo bins = 100. La escala de color Viridis, como ya habíamos visto, se agrega sumando una llamada a scale_fill_viridis_c() -porque aquí la data es continua, si fuera discreta usaríamos scale_fill_viridis_d().

ggmap(CABA) +
    geom_bin2d(data = delitos, aes(x = longitud, y = latitud), ___) +
    ___

ggmap(CABA) +
    geom_bin2d(data = delitos, aes(x = longitud, y = latitud), bins = ___) +
    scale_fill_viridis_c()

ggmap(CABA) +
    geom_bin2d(data = delitos, aes(x = longitud, y = latitud), bins = 100) +
    scale_fill_viridis_c()

Una alternativa al binning es la llamada kernel density estimation, muy utilizada en aplicaciones de análisis espacial para estimar la intensidad de una determinada variable en cualquier punto del área analizada, incluso en aquellos para los cuales no hay observaciones. La idea es asumir que los valores observados corresponden a una distribución continua sobre el espacio, y determinar cual es la más probable en base a los puntos donde existen datos. Podemos visualizar esta distribución estimada geom_density2d_filled así:

ggmap(CABA) +
    ___(data = delitos, aes(x = longitud, y = latitud), alpha = 0.5)

ggmap(CABA) +
    geom_density2d_filled(data = delitos, aes(x = longitud, y = latitud), alpha = 0.5)

Nótese que aplicamos transparencia usando el parámetro alpha = 0.5. ¿Por qué? ¿Qué pasa si lo quitamos?

Visualizando multiples categorías

Hasta aquí hemos analizado la distribución espacial de eventos en su totalidad, sin diferenciar su tipo. Veamos ahora las diferencias por categoría. Podemos reintentar el mapa de puntos, esta vez diferenciándolos por color. Recuperamos el código que usamos antes para mostrar puntos, y esta vez asignamos la columna "tipo" al atributo estético color:

ggmap(CABA) +
    geom_point(data = delitos, 
               aes(x = longitud, y = latitud, ___),
               size = 0.1, alpha = 0.1)

ggmap(CABA) +
    geom_point(data = delitos, 
               aes(x = longitud, y = latitud, color = ___),
               size = 0.1, alpha = 0.1)

ggmap(CABA) +
    geom_point(data = delitos, 
               aes(x = longitud, y = latitud, color = tipo),
               size = 0.1, alpha = 0.1)

Aquí tenemos dos problemas:

La leyenda ("tipo_delito") es difícil de leer, dado que muestra los puntos tal como los definimos: pequeños y con mucha transparencia. Esos atributos son útiles en el mapa, donde tenemos cientos de miles de puntos, pero muy poco prácticos para la leyenda, donde sólo hay un minúsculo punto por etiqueta.
Los puntos sobre el mapa se superponen en tal medida que es difícil identificar patrones espaciales distintos según su categoría.

El primer problema se resuelve fijando "a mano" los atributos de la leyenda, asi:

ggmap(CABA) +
    geom_point(data = delitos,
               aes(x = longitud, y = latitud, color = tipo),
               size = 0.1, alpha = 0.1) +
    guides(color = guide_legend(override.aes = list(size = 1, alpha = 1)))

El segundo, usando facetado para mostrar en su propio mapa a cada categoría:

ggmap(CABA) +
    geom_point(data = delitos,
               aes(x = longitud, y = latitud, color = tipo),
               size = 0.1, alpha = 0.1) +
    guides(color = guide_legend(override.aes = list(size = 1, alpha = 1))) +
    facet_wrap(vars(tipo))

El facetado ayuda a que no se nos mezclen los colores, y hace evidente cuales categorías son mas frecuentes que otras. Pero con nuestra abundancia de puntos no ayuda encontrar los sitios de alta concentración, y hace que se pierdan de vista los casos de la categoría poco frecuente (homicidios).

Para hacer las diferencias aún mas nítidas, podemos facetar una estimación de densidad en lugar de puntos. ¿Cómo lo haríamos?

ggmap(CABA) +
    ___(data = delitos, aes(x = longitud, y = latitud), alpha = 0.5) +
    facet_wrap(vars(tipo))

ggmap(CABA) +
    geom_density2d_filled(data = delitos, aes(x = longitud, y = latitud), alpha = 0.5) +
    facet_wrap(vars(tipo))

Combinando espacio y tiempo

El facetado también nos permite visualizar el cambio de posición a través del tiempo.

Por ejemplo, podemos comparar cierto tipo delito (hurto sin violencia) mostrando dónde ocurre en cada día de la semana.

Primero activamos el paquete dplyr para acceder a su función de filtrado de datos,

library(dplyr)

Y luego mostramos sólo las filas del dataframe donde la columna tipo contiene "Homicidio", - en forma de puntos en el mapa (geom_point()), - con el "subtipo" de homicidio representado por el color de los puntos - y un facetado por día de la semana (facet_wrap(vars(wday(fecha, label = TRUE))))

ggmap(CABA) +
    ___(data = filter(delitos, ___),
               aes(x = longitud, y = latitud, ___), alpha = .5) +
    ___(vars(wday(fecha, label = TRUE)))

ggmap(CABA) +
    geom_point(data = filter(delitos, tipo == ___),
               aes(x = longitud, y = latitud, color = ___), alpha = .5) +
    facet_wrap(vars(wday(fecha, label = TRUE)))

ggmap(CABA) +
    geom_point(data = filter(delitos, tipo == "Homicidio"),
               aes(x = longitud, y = latitud, color = subtipo), alpha = .5) +
    facet_wrap(vars(wday(fecha, label = TRUE)))

Vale aclarar que el poco elegante facet_wrap(vars(wday(fecha, label = TRUE))) podría cambiarse por un más legible facet_wrap(vars(dia_semana)) si como paso previo agregamos al dataframe la columna "dia_semana", guardando allí el valor obtenido con wday().

Volviendo al tiempo y el espacio, también podemos concentrarnos en un tipo de delito en particular, y evaluar en que zonas se concentra de acuerdo a la hora del día. Nuestra data ya tiene la hora del día declarada en una columna, "franja". Si no fuera así, y tuviéramos la hora como parte de la fecha (estilo "2020-09-18 14:00:00") podríamos obtenerla con ayuda de hour() que funciona de forma similar a las ya vistas month() y wday().

Entonces, mostremos sólo las filas del dataframe donde la columna tipo contiene "Hurto (sin violencia)", - en forma de mapa de densidad (geom_density2d_filled()), - y un facetado por hora del día (facet_wrap(vars(franja)))

ggmap(CABA) +
    ___(data = filter(delitos, ___),
               aes(x = longitud, y = latitud), alpha = .5) +
    ___

ggmap(CABA) +
    geom_density2d_filled(data = filter(delitos, tipo == ___),
               aes(x = longitud, y = latitud), alpha = .5) +
    facet_wrap(___)

ggmap(CABA) +
    geom_density2d_filled(data = filter(delitos, tipo == "Hurto (sin violencia)"),
               aes(x = longitud, y = latitud), alpha = .5) +
    facet_wrap(vars(franja))

En el resultado se puede ver como los hurtos se concentran nítidamente en las áreas de mayor actividad comercial durante el día (lo que los porteños llaman "el centro"), sobre todo desde el mediodía hasta 4 o 5 de la tarde, cuando pierde intensidad y se dispersa en la dirección de las principales avenidas de la ciudad.

Para terminar, pulimos la visualización - filtrando las filas que registran la franja horaria como desconocida (!is.na(franja)) - retirando la leyenda, ya que nos interesa mostrar como se mueve la densidad a lo largo del día más que las cantidades - eligiendo la cantidad de filas en la que se distribuirán las facetas (nrow = 4) - agregando título, subtítulo, y nota al pie con fuente - eligiendo un tema apropiado

ggmap(CABA) +
    geom_density2d_filled(data = filter(delitos, !is.na(franja), tipo == "Hurto (sin violencia)"),
               aes(x = longitud, y = latitud), alpha = .5) +
    guides(fill = FALSE) +
    facet_wrap(vars(franja), nrow = 4) +
    labs(title = "Ciudad de Buenos Aires: concentración espacial de hurtos",
         subtitle = "según hora del día, durante el año 2020",
         caption = "fuente: https://mapa.seguridadciudad.gob.ar") +
    theme_void()