Predicción de demanda de servicios urbanos con open data + Facebook Prophet

De todos los datasets que publica el portal de Open Data de Buenos Aires, mi favorito es sin dudas el que contiene los reclamos registrados por el Sistema Único de Atención Ciudadana (SUACI). El SUACI, también llamado BA 147, equivale a lo que en otras latitudes se conoce como servicio 311. El 311 es el número telefónico, complementado por un servicio web y en general una app también, al que los ciudadanos recurren para realizar reclamos al gobierno de la ciudad. En contraste con el servicio 911, el 311 (o 147 en Buenos Aires) se utiliza para reportar problemas que no involucran urgencias de salud o seguridad. Por ejemplo, si la cuadra de uno aparece llena de basura después de un evento multitudinario en las cercanías, se llama al 147 o se usa la app para que la ciudad envíe una cuadrilla de limpieza. En cambio, si una persona sufre un infarto en la vía pública, se llama al 911.

Si esto les resulta un poco confuso, no se preocupen, nos es culpa nuestra; hay un ligero exceso de nombres distintos para cosas similares. Los nombres SUACI y BA 147 coexisten porque -creo- SUACI registra las solicitudes al 147 pero también reclamos enviados a la ciudad por otros medios. En cuanto a 311 vs. 147, la popularidad de 311 como número para reclamos aún muy lejos del archifamoso 911. Por eso muchas ciudades en el mundo, BA incluida, coinciden en usar el 911 para emergencias pero varían en el número reservado para reclamos cotidianos.

Lo interesante de los servicios tipo 311 es que, cuando sus registros se comparten con el público, permiten hacer estudios sobre muchas facetas de la ciudad. Por ejemplo, identificar edificios peligrosos por su deterioro, o probar que las “fronteras” entre comunidades distintas dentro de la ciudad generan más reclamos que áreas homogéneas.

Un área que creo de particular interés es la de predecir demanda a futuro de servicios urbanos, para ayudar a planificar la asignación de los siempre limitados recursos públicos. Analizando la cantidad de reclamos que la población realiza a lo largo del tiempo, podemos detectar tendencias y pronosticar demanda futura, así como predecir los momentos y lugares de “tranquilidad” (que requieren menos recursos) así como aquellos que generan picos (donde se van a necesitar refuerzos).

A todo ésto, el área de I+D de Facebook liberó recientemente sus algoritmos de modelado y predicción para datos seriados en el tiempo, bajo el nombre de Prophet. La razón por la cual Prophet me resultó llamativo de inmediato es que hace muy fácil incorporar el efecto de días atípicos en un modelo predictivo para procesos que se desarrollan a lo largo del tiempo. En palabras de urbanista: podemos pronosticar la demanda de servicios urbanos usando registros históricos, generando un modelo que toma en cuenta el efecto de los diversos días feriados.

Allá vamos!

Obteniendo los datos de la ciudad

El primer paso es descargar los datos que necesitamos. Visitamos Buenos Aires Data, buscamos “Sistema Único de Atención Ciudadana”, y llegamos a la página de descarga que nos interesa - aquí. Tenemos disponible un link para descargar un archivo comprimido que contiene los datos que buscamos. Podemos dejar que R se encargue de acceder al sitio, descargar y descomprimir por nosotros:

# Donde vive nuestra data?
url <- "https://data.buenosaires.gob.ar/api/datasets/rJg_9jlR5/download"

# Donde queremos guardarla
destino <- '/home/havb/data/gcba/suaci'

# Creamos un archivo temporal para dscargar el zip con todos los datasets
temp <- tempfile()

# Descargamos -puede tomar unos cuantos minutos
download.file(url, temp)

# Des-zipeamos
unzip(temp, exdir = destino, junkpaths = TRUE)

# Eliminamos el archivo temporal
unlink(temp)

Terminado el trámite, revisamos nuestro botín:

list.files(destino)

## [1] "sistema-unico-de-atencion-ciudadana-2011.csv"
## [2] "sistema-unico-de-atencion-ciudadana-2012.csv"
## [3] "sistema-unico-de-atencion-ciudadana-2013.csv"
## [4] "sistema-unico-de-atencion-ciudadana-2014.csv"
## [5] "sistema-unico-de-atencion-ciudadana-2015.csv"
## [6] "sistema-unico-de-atencion-ciudadana-2016.csv"

He aquí el primer problema. En lugar de darnos un dataset con todo el historial de registros, nos dan una pila archivos que contienen la data separada por año. Cuando uno analiza datos a lo largo del tiempo, lo natural es tenerlos todos juntos, y luego segmentarlos -por año, por mes, o por otro período arbitrario- cuando es necesario. Y desde ya que cuando uno analiza tendencias para predecir, como estamos haciendo ahora, necesita todos los datos juntos. En fin, no es tan grave… solo es cuestión de pegar los datasets uno después del otro.

library(tidyverse)
library(lubridate)

# Obtener la lista de archivos .csv que descargamos 

archivos <- list.files(path = destino, pattern = "\\.csv", full.names = TRUE)

# A procesar!

suaci <-  archivos %>% 
  # leer cada archivo, guardar el resultado en un sólo dataframe
  map_df(read_csv2) %>% 
  # determinar la fecha a partir unieendo los campos FECHA_INGRESO, HORA_INGRESO, 
  # interpretando el resultado en formato "día/mes/año hora:min:sec [AM|PM]" 
  mutate(fecha = parse_date_time(paste(FECHA_INGRESO, HORA_INGRESO), "d/m/Y IMS p", tz = "America/Argentina/Buenos_Aires"))

Ahora, exploremos los datos unificados. Rapidito, un top ten de los reclamos más frecuentes:

suaci %>% 
  group_by(CONCEPTO, RUBRO) %>% 
  summarise(total = n()) %>% 
  arrange(desc(total)) %>%
  head(10)

## # A tibble: 10 x 3
## # Groups:   CONCEPTO [10]
##    CONCEPTO                                           RUBRO          total
##    <chr>                                              <chr>          <int>
##  1 SOLICITUD DE PARTIDAS                              REGISTRO CIV… 818151
##  2 SOLICITUD DE PAGO VOLUNTARIO DE INFRACCIONES       TRANSPORTE Y… 358298
##  3 PERSONAS SIN TECHO EVALUACION                      ATENCION SOC… 138741
##  4 LUMINARIAS APAGADAS                                ALUMBRADO     115234
##  5 SOLICITUD DE RETIRO DE RESIDUOS VOLUMINOSOS        SANEAMIENTO …  86978
##  6 SOLICITUD DE RETIRO DE RESTOS DE OBRAS O DEMOLICI… SANEAMIENTO …  81570
##  7 RETIRO DE ESCOMBROS (RESIDUOS ARIDOS)              SANEAMIENTO …  75983
##  8 LUMINARIA: APAGADA                                 ALUMBRADO      63788
##  9 VEHICULOS ABANDONADOS EN VIA PUBLICA LEY 342       TRANSPORTE Y…  47773
## 10 SOLICITUD DE REPOSICION O CAMBIO DE UBICACION DE … SANEAMIENTO …  47709

Hay algunos problemas. Categorías como “LUMINARIAS APAGADAS” y “LUMINARIA: APAGADA” son obvias referencias al mismo reclamo, registradas bajo categorías distintas. Podemos dejarlo pasar, porque no afecta a nuestros fines.

Lo que vamos a modelar son las solicitudes al departamento de saneamiento urbano: reclamos de la ciudadanía para que la ciudad retire residuos voluminosos, escombros de obra, ramas podadas, etc. Aquí vale la pena insistir en la riqueza del dataset. Entre las categorías que aparecen con solo mirar las principales, aparecen la atención a personas sin techo, y los problemas de iluminación pública. Hay, tantas, tantas cosas que se pueden hacer mediante el análisis espacial y temporal de la data!

Aislamos los reclamos que nos interesan en este momento:

solicitudes_saneamiento <- suaci %>%
  filter(grepl("RESIDUOS|RESTOS|ESCOMBROS ", CONCEPTO) == TRUE & 
           RUBRO == "SANEAMIENTO URBANO") %>% 
  group_by(dia = format(fecha, "%Y-%m-%d")) %>% 
  summarise(total = n())

Y ahora los graficamos para ver que pinta tienen:

library(hrbrthemes)

ggplot(solicitudes_saneamiento, aes(as.Date(dia), total)) + 
  geom_line(color = "orange") +
  scale_x_date(date_labels = "%Y-%m") +
  theme_ipsum() +
  labs(y = "solicitudes", x = "fecha",
       title="Demanda diaria de servicios de saneamiento urbano",
       subtitle = "Ciudad de Buenos Aires: 2011 - 2016",
       caption = "Fuente: https://data.buenosaires.gob.ar/")

A primera vista, es obvio que los datos están sesgados, ya que la explosión en el número de solicitudes a partir de mediados del 2012 es atribuible a la forma en que la ciudad empezó a tomar nota de los reclamos, más que a un furor ciudadano por solicitar limpieza. Cuando hagamos nuestro modelo, vamos a descartar los registros de 2011 y 2012 porque tenemos claro que no son representativos.

El siguiente ingrediente que necesitamos es una lista de los feriados públicos en la Argentina durante el período analizado.

Cuando los datos están dispersos: la hora del scraping

Buscar una lista oficial de feriados resultó en frustración. Si bien el gobierno nacional publica una lista de los feriados vigentes para el año en curso, no existe un archivo para descargar con las fechas exactas de los feriados en años anteriores.

En los sitios web de varios diarios locales encontramos un historial de feriados, pero dispersos en distintas páginas web, y sin opción para descargarlos en un archivo de texto. Vamos a tener que “scrapear” la data. Cómo nos hacen laburar, che. Si tan solo alguien nos diera una API para toda información de consulta permanente, no tendríamos que hacer estas cosas!

ira

Nuestro proveedor de fechas de feriados será La Nación, que publica un bonito calendario de feriados oficiales, con la posibilidad de consultar los de años anteriores.

library(rvest)

# Una función que hace scraping de todos los feriados publicados en lanacion.com.ar para un año dado

scrape_feriados_del_anio <- function(year) {
  
  # Donde está la data de los feriados
  baseurl <- "http://servicios.lanacion.com.ar/feriados/"
  
  get_feriados_mes <- function(month, feriados_page) {
    
    # Los feriados aparecen en dos categorías: "inamovible" y "trasladable"
    
    inamovibles <- feriados_page %>%
      map(html_nodes, "li.inamovible") %>%
      map(html_text) %>%
      .[[month]] %>%
      {if (!is_empty(.)) paste(year, month, ., sep = "/")}
  
    trasladables  <- feriados_page %>%
      map(html_nodes, "li.trasladable") %>%
      map(html_text) %>% 
      .[[month]] %>%
      {if (!is_empty(.)) paste(year, month, ., sep = "/")}
    
    todos <- c(inamovibles, trasladables) 
    
    if (!is.null(todos)) return(todos)
    
    }
  
  feriados_page <- read_html(paste0(baseurl, year)) %>% 
    html_nodes(".bloque")
  
  feriados_del_anio <- map(1:12, get_feriados_mes, feriados_page) %>% 
    unlist()
  
  return(feriados_del_anio)
}
  

# Descargamos los feriados de 2011 a 2016,  
feriados_2013_2016 <- map(2013:2016, scrape_feriados_del_anio) %>% 
  #los unimos en un unico vector
  unlist %>% 
  # los definimos como fecha
  ymd %>% 
  # Los ordenamos cronológicamente (necesario porque el scraper los trae como texto 
  # en orden alfabetico)
  sort

Tenemos nuestra lista de feriados?

feriados_2013_2016

##  [1] "2013-01-01" "2013-01-31" "2013-02-11" "2013-02-12" "2013-02-20"
##  [6] "2013-03-24" "2013-03-29" "2013-04-01" "2013-04-02" "2013-05-01"
## [11] "2013-05-25" "2013-06-20" "2013-06-21" "2013-07-09" "2013-08-19"
## [16] "2013-10-14" "2013-11-25" "2013-12-08" "2013-12-25" "2014-01-01"
## [21] "2014-03-03" "2014-03-04" "2014-03-24" "2014-04-02" "2014-04-18"
## [26] "2014-05-01" "2014-05-02" "2014-05-25" "2014-06-20" "2014-07-09"
## [31] "2014-08-18" "2014-10-13" "2014-11-24" "2014-12-08" "2014-12-25"
## [36] "2014-12-26" "2015-01-01" "2015-02-16" "2015-02-17" "2015-03-23"
## [41] "2015-03-24" "2015-04-02" "2015-04-03" "2015-05-01" "2015-05-25"
## [46] "2015-06-20" "2015-07-09" "2015-08-17" "2015-10-12" "2015-11-27"
## [51] "2015-12-07" "2015-12-08" "2015-12-25" "2016-01-01" "2016-02-08"
## [56] "2016-02-09" "2016-03-24" "2016-03-25" "2016-04-02" "2016-05-01"
## [61] "2016-05-25" "2016-06-17" "2016-06-20" "2016-07-08" "2016-07-09"
## [66] "2016-08-15" "2016-10-10" "2016-11-28" "2016-12-08" "2016-12-09"
## [71] "2016-12-25"

Oh si.

Modelando y prediciendo

Con todos los ingredientes a mano, es hora de hacer vaticinios. Como suele pasar cuando uno trabaja con datos, hacer el modelo es la parte más sencilla… la mayor parte del tiempo la empleamos en reunir y limpiar los datos!

Creamos un modelo:

library(prophet)

modelo <- solicitudes_saneamiento %>% 
  filter(year(as.Date(dia)) > 2012) %>% 
  transmute(ds = as.Date(dia), y = total) %>%  
  prophet(holidays = data.frame(holiday = "feriado", ds = feriados_2013_2016))

## Initial log joint probability = -116.558
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance

Y predecimos la demanda para el año siguiente (todo el 2017):

forecast <- predict(modelo, make_future_dataframe(modelo, periods = 365))

De paso, hacemos dos preguntas rápidas. Cuantas solicitudes diarias recibe la ciudad, y que efecto tiene un día feriado en el nivel de demanda?

Durante el período 2013-2016, el área de saneamiento urbano de la ciudad recibió un promedio de 394 reclamos diarios, con un máximo 1075. El día en el que menos reclamos se registraron sólo hubo 3.

modelo$history$y %>% summary

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     3.0   121.0   392.0   394.1   603.5  1075.0

Según nuestro modelo, el efecto de que un feriado caiga en un día particular del mes es una reducción de 401 reclamos. En un día típico, esto haría que prácticamente no haya reclamos.

forecast %>% 
  select(ds, feriado) %>% 
  filter(abs(feriado) > 0) %>% 
  head()

##           ds   feriado
## 1 2013-01-01 -401.1635
## 2 2013-01-31 -401.1635
## 3 2013-02-11 -401.1635
## 4 2013-02-12 -401.1635
## 5 2013-02-20 -401.1635
## 6 2013-03-24 -401.1635

Visualizando por separado la tendencias general y las periódicas (por día de la semana y día del mes) notamos varios efectos.

La tendencia general fue una baja del nivel de demanda desde el 2014 hasta el 2016, revertida luego. Se evidencia una fuerte suba de allí en más. Valdría la pena discernir si esto se debe a que la data sub-representa los reclamos del 2014-2016, o si efectivamente ocurrió una baja de demanda en esos años. Si ésto último fuera el caso, a que podría deberse?
El día de mayor actividad es el Lunes, a continuación de los días más tranquilo, los del fin de semana. Está claro que en Domingo la gente no reclama mucho, pero al día siguiente si. En el resto de los días la demanda es pareja
Tal como habíamos comprobado revisando los números, los días feriados generan una caída de unos 400 reclamos
Los meses de mayor actividad son los de la segunda mitad del año. Durante las vacaciones de invierno se observa una baja de la demanda, leve en comparación a la de fin de año…. a partir de mediados de Diciembre, la demanda cae en picada, y se mantiene mínima en Enero.

prophet_plot_components(modelo, forecast)

Y por último, trazamos la predicción de nivel de demanda para lo que queda del año:

plot(modelo, forecast) + theme_ipsum() +
  labs(y = "solicitudes", x = "fecha",
       title="Pronóstico: demanda diaria de servicios de saneamiento urbano",
       subtitle = "Ciudad de Buenos Aires: 2013 - 2016, extendido a final del 2017")

Yo diría que, a no ser que esté operando con capacidad de sobra, Saneamiento Urbano va a necesitar más recursos!

Para la próxima: hacer análisis espacial además del temporal.

prophet R rvest