Introducción

En este proyecto se analizará la demanda de pasajeros en el Subte de Buenos Aires durante el mes de abril 2019. Esta información se obtendrá de una base de datos que contiene ingresos horarios por molinetes 2 en la red de subterráneos 3.

A continuación, según DGEyC (2018), se definen los tipos de pasajeros que componen la demanda a estudiar:

A su vez, se buscará una relación entre la información de pasajeros totales y los datos horarios recopilados por las estaciones meteorológicas pertenecientes al Servicio Meteorológico Nacional (SMN). Esta base de datos contiene características climáticas como temperatura (TEMP), presión barométrica (PNM), presencia de precipitaciones (PP), intensidad (FF) y dirección del viento (DD), y humedad (HUM).

En vista de que el presente estudio está inspirado en el trabajo realizado por Evanoff (2015), surge la siguente pregunta de investigación: ¿viaja más gente en el Metro cuando llueve?

A través de las técnicas estadísticas aplicadas, Evanoff (2015) logró confirmar que la demanda de pasajeros del Metro de Nueva York crece, conforme se intensifican las lluvias.

Ahora bien, ¿esta situación se sostiene para el Subte de Buenos Aires? En caso contrario, ¿cuál sería la causa de que los usuarios no elijan este modo de transporte?

Presentación de los datos

  1. Base Molinetes por estación por hora - días hábiles de abril 2019. Fuente: SBASE.
  2. Datos meteorológicos horarios - días hábiles de abril 2019. Fuente: SMN.
  3. Registro de precipitaciones horario - días hábiles de abril 2019. Fuente: Twitter SMN-OCBA (@SMN_OCBA).
  4. Información censal por radio. Fuente primaria: INDEC.
  5. Reclamos a través del Sistema Único de Atención Ciudadana - año 2019. Fuente primaria: Jefatura de Gabinete de Ministros.
  6. Estaciones y líneas de Subte. Fuente primaria: SBASE.

Como primer paso, se eliminarán los archivos del ambiente de trabajo de R y se seteará a las variables cuantitativas sin notación científica:

rm(list=ls())
options(scipen=999)

A continuación, se cargarán los paquetes y se activarán las librerías útiles:

install.packages("tidyverse", repos = "https://cran.rstudio.com", dependencies = TRUE)
install.packages("scales", repos = "https://cran.rstudio.com", dependencies = TRUE)
install.packages("dplyr", repos = "https://cran.rstudio.com", dependencies = TRUE)
install.packages("GGally", repos = "https://cran.rstudio.com", dependencies = TRUE)
install.packages("gridExtra", repos = "https://cran.rstudio.com", dependencies = TRUE)
install.packages("corrplot", repos = "https://cran.rstudio.com", dependencies = TRUE)
install.packages("stargazer", repos = "https://cran.rstudio.com", dependencies = TRUE)
install.packages("viridis", repos = "https://cran.rstudio.com", dependencies = TRUE)
install.packages("readxl", repos = "https://cran.rstudio.com", dependencies = TRUE)
install.packages("caTools", repos = "https://cran.rstudio.com", dependencies = TRUE)
install.packages("lubridate", repos = "https://cran.rstudio.com", dependencies = TRUE)
install.packages("sf", repos = "https://cran.rstudio.com", dependencies = TRUE)
install.packages("knitr", repos = "https://cran.rstudio.com", dependencies = TRUE)
install.packages("png", repos = "https://cran.rstudio.com", dependencies = TRUE)
install.packages("sp", repos = "https://cran.rstudio.com", dependencies = TRUE)
install.packages("latticeExtra", repos = "https://cran.rstudio.com", dependencies = TRUE)
library(tidyverse)
library(scales)
library(dplyr)
library(GGally)
library(gridExtra)
library(corrplot) 
library(stargazer)
library(viridis)
library(readxl)
library(caTools)
library(lubridate)
library(sf)
library(knitr)
library(png)
library(sp)
library(latticeExtra)

Manejo de los datos

Como el objetivo principal de esta investigación es analizar la intensidad de la relación entre el uso del Subte y la presencia de precipitaciones, se cargarán al ambiente de trabajo la base Molinetes y de Datos Horarios de temperatura correspondientes a los días hábiles del mes de abril 2019.

Previo a este paso, se deben chequear en la aplicación TweetDeck los tweets horarios realizados por la cuenta @SMN_OCBA (perteneciente al Observatorio Ciudad de Buenos Aires) durante el período en análisis, que mencionen la palabra “precipitación”, “lluvia”, “llovizna”, o las abreviaturas “precip” o “ppcion”, y registrar en la base de Datos Horarios un 1 cuando se presente el evento, y 0 en caso contrario.

Para el armado de esta base, se pegarán las filas correspondientes a cada día, unas debajo de las otras, en orden calendario:

molinetes_crudo <- read_excel("turnstiles_crudo_Apr19.xlsx")

datohorario_0401 <- read_excel("datohorario20190401.xlsx")
datohorario_0403 <- read_excel("datohorario20190403.xlsx")
datohorario_0404 <- read_excel("datohorario20190404.xlsx")
datohorario_0405 <- read_excel("datohorario20190405.xlsx")
datohorario_0408 <- read_excel("datohorario20190408.xlsx")
datohorario_0409 <- read_excel("datohorario20190409.xlsx")
datohorario_0410 <- read_excel("datohorario20190410.xlsx")
datohorario_0411 <- read_excel("datohorario20190411.xlsx")
datohorario_0412 <- read_excel("datohorario20190412.xlsx")
datohorario_0415 <- read_excel("datohorario20190415.xlsx")
datohorario_0416 <- read_excel("datohorario20190416.xlsx")
datohorario_0417 <- read_excel("datohorario20190417.xlsx")
datohorario_0422 <- read_excel("datohorario20190422.xlsx")
datohorario_0423 <- read_excel("datohorario20190423.xlsx")
datohorario_0424 <- read_excel("datohorario20190424.xlsx")
datohorario_0425 <- read_excel("datohorario20190425.xlsx")
datohorario_0426 <- read_excel("datohorario20190426.xlsx")
datohorario_0429 <- read_excel("datohorario20190429.xlsx")

dato_horario <- bind_rows(datohorario_0401, datohorario_0403, datohorario_0404, datohorario_0405, datohorario_0408, datohorario_0409, datohorario_0410, datohorario_0411, datohorario_0412, datohorario_0415, datohorario_0416, datohorario_0417, datohorario_0422, datohorario_0423, datohorario_0424, datohorario_0425, datohorario_0426, datohorario_0429)

Con el fin de emprolijar los datos para que sean más fáciles de usar, en primer lugar, se incluirá en la base Molinetes una columna con el día de la semana correspondiente a cada fecha. En segundo lugar, se creará una primary key para unir ambas bases y, por último, se filtrará a la base de Datos Horarios por estación meteorológica “Buenos Aires”:

molinetes_crudo <- molinetes_crudo %>%
  mutate(dia_semana = wday(fecha1, label = TRUE),
         fecha_horario = paste(fecha1, hora, sep = " ")) 

dato_horario <- dato_horario %>%
  filter(NOMBRE == "BUENOS AIRES") %>%
  mutate(fecha_horario = paste(FECHA, HORA, sep = " "),
         PP = as.numeric(PP))

molinetes <- molinetes_crudo %>% 
  left_join(dato_horario) %>%
  select(-fecha_horario, -FECHA, -HORA, -NOMBRE, -pax_pases_pagos, -pax_franq)
## Joining, by = "fecha_horario"

Análisis Exploratorio

Correlaciones entre variables

De acuerdo con Novales (2010), la correlación es la fuerza, dirección y proporcionalidad de una relación lineal entre dos variables estadísticas. Dos variables cuantitativas están correlacionadas si al disminuir los valores de una lo hacen también los de la otra y viceversa.

Para el conjunto de datos del presente proyecto, la matriz de correlación será bastante difícil de interpretar dado que, luego del procedimiento de data wrangling, la base Molinetes contiene siete variables cuantitativas.

Con el fin de resumir esta información, mediante la función corrplot, se obtendrá una visualización gráfica de dicha matriz:

numerico <- molinetes[,sapply(molinetes, is.numeric)]
matriz_corr <- cor(numerico)
corrplot(matriz_corr, method = "square", order="hclust", type="lower", 
         tl.cex=0.75, add = FALSE, tl.pos="lower") 
corrplot(matriz_corr, method = "square", order="hclust", type="upper", 
         tl.cex=0.75, add = TRUE, tl.pos="upper")

Como era de esperarse, las variables que están fuertemente asociadas entre sí son humedad y presencia de precipitaciones, y cantidad de pasajeros pagos con totales.

Pasajeros totales durante días hábiles

El siguiente análisis estará orientado a explicar la dinámica de la demanda de pasajeros sin tener en cuenta los factores meteorológicos. Para eso, se creará una tabla intermedia con la información de pasajeros totales por estación:

rank_estaciones <- molinetes %>%
    group_by(nombre_estacion) %>%
    summarise(ingresos = sum(pax_TOTAL)) %>%
    arrange(-ingresos)

Una vez que se cuenta con la información de pasajeros totales por estación, se graficará el top 25 de estaciones durante el mes en estudio:

##Pasajeros totales por estación - Top 25
ggplot(rank_estaciones[1:25,], aes(x = reorder(nombre_estacion, -ingresos), ingresos)) +
  geom_bar(stat = 'identity', fill = "#FFDB6D", color = "#C4961A") +
  theme(axis.text.x = element_text(angle = 90)) +
  scale_y_continuous() + 
  xlab("Estación") + ylab("Ingresos totales") +
  ggtitle("Pasajeros totales por estación - Top 25 estaciones. \nSubte de Buenos Aires. Días hábiles abril 2019.") 

Como resultado del gráfico, la cantidad de pasajeros totales en la estación Constitución de la Línea C está muy por encima de la estación inmediata posterior Congreso de Tucumán de la Línea D. Esto es así porque, en línea con BID y SBASE (2015), Constitución está ubicada aledaña a la terminal ferroviaria del General Roca (FCGR) que conecta CABA con el conurbano bonaerense sur.

En esta línea, es posible verificar que en el top 25, figuran la mayoría de las estaciones de la red con las que se puede trasbordar a un ferrocarril: Federico Lacroze, Retiro, Rosas, San Pedrito y Palermo, de ahí la importancia de un sistema intermodal de transporte.

Pasajeros totales durante las horas del día

A continuación, se creará una tabla intermedia con la información de pasajeros totales, seleccionando las variables hora, línea y estación:

horas <- molinetes %>%
  select(hora, nombre_estacion, linea, pax_TOTAL) %>%
  group_by(hora, linea, nombre_estacion) %>%
  summarise(total_ingresos = sum(pax_TOTAL)) %>%
  arrange(hora, desc(total_ingresos))

Una vez que se cuenta con la información de pasajeros totales por hora, se graficará la demanda de pasajeros según hora del día por línea y por estación:

##Pasajeros totales por hora por línea
ggplot(horas, aes(x = hora, y = total_ingresos, fill = linea)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle=45)) +
  scale_y_continuous() +
  scale_x_continuous(limit=c(5,23)) +
  scale_fill_manual(values = c("deepskyblue", "darkred", "dodgerblue4", "forestgreen", "darkorchid4", "gold")) +
  guides(fill = FALSE) +
  facet_grid(~linea) +
  xlab("Hora") + ylab("Ingresos totales") + 
  ggtitle("Pasajeros totales por hora por línea \nSubte de Buenos Aires. Días hábiles abril 2019.") 

##Pasajeros totales por hora por estación
ggplot(horas, aes(hora, total_ingresos, fill = as.factor(hora))) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle=45)) +
  guides(fill = FALSE) +
  scale_y_continuous() +
  scale_x_continuous(limit=c(4,23)) +
  scale_fill_viridis(discrete = TRUE, option = "magma") +
  guides(fill = FALSE) +
  facet_wrap(~nombre_estacion) +
  xlab("Hora") + ylab("Ingresos totales") + 
  ggtitle("Pasajeros totales por estación por hora. \nSubte de Buenos Aires. Días hábiles abril 2019.") 

En concreto, es notoria la diferencia de pasajeros totales entre las horas pico (8-10hs. y 18-20hs.) y las horas valle en el total de la red. Por un lado, B y D son las líneas que transportan mayor cantidad de pasajeros debido a que su traza pasa por debajo de las avenidas Corrientes y Cabildo - Santa Fe respectivamente, las cuales cuentan con importantes centros comerciales (BID y SBASE, 2015).

En contraposición, las líneas que transportan menor cantidad de pasajeros son H y E. De acuerdo con enelSubte.com (2018), la H es una de las líneas que más ha crecido en materia de pasajeros y es debido a su carácter transversal, que combina con el resto de las líneas de la red en el corredor Pueyrredón - Jujuy, a excepción de la C.

Por otro lado, pese a que Constitución es la estación con mayor cantidad de pasajeros totales, esta situación se acentúa en la hora pico de la mañana por sobre la hora pico de la tarde, de donde se infiere que la mayor cantidad de pasajeros que utiliza el Subte, reside en zona sur y trabaja en CABA. Todas estas observaciones se relacionan también con que las estaciones Catedral de la Línea D y Leandro N. Alem de la Línea B tienen su pico en el horario de las 18-20hs. Sería interesante poder validar esta hipótesis mediante los datos de operaciones con Sistema Único de Boleto Electrónico (SUBE) publicados por el Ministerio de Transporte, pero excede el presente estudio.

Impacto de las precipitaciones en el servicio Subte

Ahora bien, para avanzar sobre el objetivo general de este proyecto, será preciso crear dos tablas intermedias: una con la cantidad de validaciones por molinete por día y otra con la cantidad de validaciones totales durante las horas con presencia de precipitaciones:

validaciones <- molinetes %>%
  group_by(fecha1) %>%
  summarise (validaciones = sum(pax_pagos))
kable(head(validaciones, n = 10),
      format = "pandoc",
      col.names = c("Fecha", "Validaciones"),
      align = c("c"))
Fecha Validaciones
2019-04-01 1232971
2019-04-03 1293758
2019-04-04 1302598
2019-04-05 1312744
2019-04-08 1259626
2019-04-09 1271521
2019-04-10 1279026
2019-04-11 1350940
2019-04-12 1289296
2019-04-15 1253173
precipitaciones <- molinetes %>%
  group_by(fecha1, hora) %>%
  filter(PP == 1) %>%
  summarise(Pasajeros = sum(pax_pagos)) %>%
  arrange(fecha1, hora)
kable(head(precipitaciones, n = 10),
      format = "pandoc",
      col.names = c("Fecha", "Hora", "Pasajeros"),
      align = c("c"))
Fecha Hora Pasajeros
2019-04-04 13 83202
2019-04-04 14 77687
2019-04-04 15 81691
2019-04-04 16 89718
2019-04-10 22 18721
2019-04-10 23 1614
2019-04-11 7 90331
2019-04-11 8 118897
2019-04-11 9 97159
2019-04-26 6 32957

Luego, se graficará el registro de precipitaciones por hora en el total de la red:

##Histograma de pasajeros totales por hora
ggplot(molinetes, aes(x=pax_TOTAL, fill = PP == 1)) + 
  geom_histogram(alpha=1, binwidth=100, color="black") + 
  coord_cartesian(xlim=c(0,3500), ylim=c(0, 3000)) + 
  scale_y_continuous(breaks=seq(0, 3500, by=500)) +
  scale_fill_manual(name = "¿Precipitaciones?",
                    labels = c("No", "Sí"),
                    values = c("#FC4E07", "#00AFBB")) +
  xlab("Pasajeros totales por hora") + 
  ylab("Frecuencia") + 
  ggtitle("Histograma de pasajeros totales: \nPresencia vs ausencia de lluvia. \nSubte de Buenos Aires. Días hábiles abril 2019.") 

En función de la información meteorológica recabada, 4/18 días tuvieron presencia de precipitaciones, siendo el viernes 26/04 el día que presentó lluvias que duraron casi todo el día.

Sin embargo, en materia de pasajeros, para el 26/04 hubo un total de 1.187.048 validaciones a diferencia del día jueves 11/04 que hubo un total de 1.350.940, siendo el día hábil del mes con mayor cantidad de validaciones.

Modelado estadístico - Regresión lineal

Variables ficticias

Al realizar un análisis de regresión, las variables independientes pueden ser tanto variables cuantitativas como variables de carácter cualitativo, llamadas dummies (Novales, 2010). Dichas variables indican la presencia (=1) o ausencia de un atributo (=0).

Para el conjunto de datos en estudio, se convertirán en factores a las variables estación, día de la semana y hora:

estacion_factor <- factor(molinetes$nombre_estacion)
dia_semana_factor <- factor(molinetes$dia_semana)
hora_factor <- factor(molinetes$hora)
f1 <- model.matrix(~ estacion_factor -1)
f2 <- model.matrix(~ dia_semana_factor -1)
f3 <- model.matrix(~ hora_factor -1)

Variables de interacción

En este proyecto, el modelo de regresión se utilizará con el interés de estimar la relación de una o más variables independientes con la variable dependiente.

En esta línea, Novales (2010) define a la interacción como la interferencia que una o varias variables pueden tener en la asociación con otras. Con este fin, se crearán las siguientes variables de interacción:

  • Estación / Día
  • Estación / Hora
  • Día / Hora
estacion_dia <- paste0(as.character(molinetes$nombre_estacion), '_' , as.character(molinetes$dia_semana))
estacion_hora <- paste0(as.character(molinetes$nombre_estacion), '_' , as.character(molinetes$hora))
dia_hora <- paste0(as.character(molinetes$dia_semana), '_' , as.character(molinetes$hora))
f4 <- model.matrix(~factor(estacion_dia) -1)
f5 <- model.matrix(~factor(estacion_hora) -1)
f6 <- model.matrix(~factor(dia_hora) -1)
dummy <- cbind(f1,f2,f3,f4,f5,f6)

Medidas de varianza explicada

En esta investigación, se utilizarán dos medidas de varianza explicada. Por un lado, el coeficiente de correlación R2 reflejará la bondad del ajuste de un modelo a la variable que pretende explicar.

Por otro lado, el R2 ajustado se utilizará para ver el grado de intensidad o efectividad que tienen las variables independientes en explicar la variable dependiente (Novales, 2010).

A continuación, se construirán funciones de R2 y R2 ajustado, en función de sus expresiones matemáticas:

r_squared <- function(y, pred) {
  RSS <- sum((y - pred )^2)
  TSS <- sum((y - mean(y))^2)
  r_sq <- 1 - (RSS / TSS)
}

adj_r_squared <- function(r_sq, num_examples, num_features) {
  adj_r <- 1 - (1 - r_sq) * ((num_examples - 1) / (num_examples - num_features - 1))
  return(adj_r)
}

Entrenamiento del modelo

En el siguiente punto, se dividirá el total de observaciones de la población en dos muestras a partir de una semilla aleatoria: 70% entrenamiento y 30% prueba.

set.seed(100)
split <- sample.split(molinetes$pax_TOTAL, 0.7)
train_df <- subset(molinetes, split == TRUE)
test_df <- subset(molinetes, split == FALSE)

Selección de modelos lineales

En cuanto a la selección del modelo lineal óptimo, se construirán cuatro modelos lineales en función de las observaciones de la muestra de entrenamiento, y se evaluarán potenciales variables predictoras:

  1. Incorporar solamente la variable estación.
  2. Al modelo 1, sumarle la variable hora.
  3. Al modelo 2, sumarle la variable día de la semana.
  4. Hacer interactuar a las variables del modelo 3 entre sí, y sumar las variables referidas a indicadores meteorológicos.
train_predictores <- train_df %>%
  select(pax_TOTAL, hora, dia_semana, nombre_estacion, TEMP, HUM, PP)

modelo1 <- lm(pax_TOTAL ~ nombre_estacion, data = train_predictores)
modelo2 <- lm(pax_TOTAL ~ nombre_estacion + hora, data = train_predictores)
modelo3 <- lm(pax_TOTAL ~ nombre_estacion + hora + dia_semana, data = train_predictores)
modelo4 <- lm(pax_TOTAL ~ PP + HUM + TEMP + nombre_estacion*hora + hora:dia_semana + nombre_estacion:dia_semana, data = train_predictores)

stargazer(modelo1, modelo2, 
          title = "Resultados de la regresión",
          dep.var.caption = "Pasajeros totales",
          header = FALSE, 
          single.row = TRUE, 
          type = "text", 
          font.size = "tiny", 
          digits.extra = 4, 
          digits = 4, 
          keep="NULL", 
          column.labels = c("Modelo 1", "Modelo 2"))
## 
## Resultados de la regresión
## =============================================================================
##                                         Pasajeros totales                    
##                     ---------------------------------------------------------
##                                             pax_TOTAL                        
##                               Modelo 1                     Modelo 2          
##                                 (1)                          (2)             
## =============================================================================
## Observations                   20,866                       20,866           
## R2                             0.4604                       0.4675           
## Adjusted R2                    0.4582                       0.4653           
## Residual Std. Error    841.5500 (df = 20779)        836.0219 (df = 20778)    
## F Statistic         206.1940*** (df = 86; 20779) 209.7090*** (df = 87; 20778)
## =============================================================================
## Note:                                             *p<0.1; **p<0.05; ***p<0.01
stargazer(modelo3, modelo4, 
          title = "Resultados de la regresión (continuación)",
          dep.var.caption = "Pasajeros totales",
          header = FALSE, 
          single.row = TRUE, 
          type = "text", 
          font.size = "tiny", 
          digits.extra = 4, 
          digits = 4, 
          keep="NULL", 
          column.labels = c("Modelo 3", "Modelo 4"))
## 
## Resultados de la regresión (continuación)
## =============================================================================
##                                         Pasajeros totales                    
##                     ---------------------------------------------------------
##                                             pax_TOTAL                        
##                               Modelo 3                     Modelo 4          
##                                 (1)                          (2)             
## =============================================================================
## Observations                   20,866                       20,866           
## R2                             0.4677                       0.6998           
## Adjusted R2                    0.4653                       0.6901           
## Residual Std. Error    836.0140 (df = 20774)        636.4271 (df = 20216)    
## F Statistic         200.5430*** (df = 91; 20774) 72.6058*** (df = 649; 20216)
## =============================================================================
## Note:                                             *p<0.1; **p<0.05; ***p<0.01

Bondad de ajuste del Modelo 4

Para evaluar la bondad de ajuste de un modelo, según Novales (2010), un valor de R2 cercano a 0, debe interpretarse como que no existe un buen ajuste lineal, lo que no excluye que existan otros tipos de relaciones funcionales.

Además, acorde con lo dicho por Novales (2010), un coeficiente de correlación muy cercano a 1 o a −1, no debe interpretarse como que existe una relación causa - efecto importante entre las dos variables. La relación podría deberse al efecto de otras variables no incluidas en el estudio, lo que se llama relación espuria.

Para eso, se analizarán las medidas de varianza explicada del Modelo 4, que contiene los indicadores meteorológicos:

train_pred <- predict(modelo4, train_df)
train_y <- train_df$pax_TOTAL
train_r_2 <- r_squared(train_y, train_pred)
train_adj_r_2 <- adj_r_squared(train_r_2, dim(train_df)[1], length(coef(modelo4))-1)

test_pred <- predict(modelo4, newdata = test_df)
test_y <- test_df$pax_TOTAL
test_r_2 <- r_squared(test_y, test_pred)

df_r2 <- data.frame(train_r_2, train_adj_r_2, test_r_2)
kable(df_r2, 
      digits = 4,
      format = "pandoc",
      col.names = c("Train R<sup>2</sup>", "Train R<sup>2</sup> ajustado", "Test R<sup>2</sup>"),
      align = c("c"))
Train R2 Train R2 ajustado Test R2
0.6998 0.6901 0.4245
summary(modelo4)$coef[2] ##Coeficiente que acompaña variable PP
## [1] 431.4085

Como era de esperarse, el Modelo 4 es el modelo que mejor representa la realidad, ya que, para el conjunto de datos seleccionado, tiene una bondad de ajuste del 69%, lo que implica una correlación positiva moderada. Por lo tanto, de acuerdo con la tabla Resultados de la Regresión, los modelos 1, 2 y 3, tienen R2 por debajo de 0,5, lo que implica una relación débil entre las variables explicativas y explicadas (Novales, 2010).

En resumen, en el Modelo 4, el coeficiente de la variable “precipitaciones” es 431. Esto significa que, ante la presencia de lluvias, habrá 431 pasajeros totales más, a diferencia de lo que se verifica en los datos crudos (ver Histograma de pasajeros totales).

Residuos del set de entrenamiento

pred <- predict(modelo4)

p1 <- ggplot() + geom_histogram(aes(modelo4$residuals), binwidth = 200) +
  xlab('Residuos') + ylab('Frecuencias') +
  ggtitle('Histograma de residuos')
p2 <- ggplot() + geom_point(aes(pred, modelo4$residuals)) +
  geom_hline(yintercept = 0, linetype = 'dotted', size = 1.2, color = 'red') + 
  xlab('Predicciones') + ylab('Residuos') + 
  ggtitle('Residuos vs predicciones')
p3 <- ggplot() + geom_point(aes(pred, train_df$pax_TOTAL)) + 
  geom_abline(slope = 1, linetype = 'dotted', color = 'red', size = 1.2) + 
  xlab("Predicciones") + ylab("Pasajeros totales por hora") + 
  ggtitle("Predicciones vs pasajeros \ntotales por hora")
p4 <- ggplot() + geom_point(aes(1:nrow(train_df), modelo4$residuals), data = train_df) +
  geom_hline(yintercept = 0, linetype = 'dotted', size = 1.2, color = 'red') + 
  xlab('Índice') + ylab('Residuos') + 
  ggtitle('Residuos en orden \nsecuencial por hora \npor estación')
grid.arrange(p1, p2, p3, p4, ncol = 2)

Residuos del set de test

test_resid <- test_y - test_pred

p5 <- ggplot() + geom_histogram(aes(test_resid), binwidth = 200) +
  xlab('Residuos') + ylab('Frecuencias') + xlim(c(-20000, 20000)) +
  ggtitle('Histograma de residuos')
p6 <- ggplot() + geom_point(aes(test_pred, test_resid)) +
  geom_hline(yintercept = 0, linetype = 'dotted', size = 1.2, color = 'red') + 
  xlab('Predicciones') + ylab('Residuos') + 
  ggtitle('Residuos vs predicciones')
p7 <- ggplot() + geom_point(aes(test_pred, test_y)) + 
  geom_abline(slope = 1, linetype = 'dotted', color = 'red', size = 1.2) + 
  xlab("Predicciones") + ylab("Pasajeros totales por hora") + 
  ggtitle("Predicciones vs pasajeros \ntotales por hora")
p8 <- ggplot() + geom_point(aes(1:nrow(test_df), test_resid), data = test_df) +
  geom_hline(yintercept = 0, linetype = 'dotted', size = 1.2, color = 'red') + 
  xlab('Índice') + ylab('Residuos') + 
  ggtitle('Residuos en orden \nsecuencial por hora \npor estación')
grid.arrange(p5, p6, p7, p8, ncol = 2)

Los gráficos presentados anteriormente muestran que los residuos presentan una distribución aproximadamente normal con una compresión cercana a 0, que se desplega en los extremos superiores. A su vez, nos indican que la predicción resultaría insuficiente para los niveles más altos de pasajeros.

No obstante, hasta este momento el análisis desarrollado no es condición suficiente para ratificar la hipótesis de que hay mayor cantidad de personas que viajarán en el Subte de Buenos Aires cuando llueve a comparación de cuando no llueve.

A continuación, se analizará una posible causa que refute este fenómeno: ¿será que los usuarios no eligen el Subte como medio de transporte cuando llueve porque hay filtraciones en las estaciones? Esta situación puede darse por múltiples razones: lluvias intensas, rotura de cañerías u otros problemas provenientes de la superficie.

Con este fin, se utilizará la base Radios elaborada por INDEC para el Censo Nacional de Población, Hogares y Viviendas (CNPHV) 2010 y la base de reclamos a través del Sistema Único de Atención al Ciudadano (SUACI) del año 2019, disponible en el portal Buenos Aires Data:

radios <- st_read("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/informacion-censal-por-radio/CABA_rc.geojson")
## Reading layer `CABA_rc' from data source `http://cdn.buenosaires.gob.ar/datosabiertos/datasets/informacion-censal-por-radio/CABA_rc.geojson' using driver `GeoJSON'
## Simple feature collection with 3554 features and 8 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -58.53092 ymin: -34.70574 xmax: -58.33455 ymax: -34.528
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
names(radios)
## [1] "RADIO_ID"    "BARRIO"      "COMUNA"      "POBLACION"   "VIVIENDAS"  
## [6] "HOGARES"     "HOGARES_NBI" "AREA_KM2"    "geometry"
suaci <- read.csv("suaci_2019.csv", 
                  encoding="UTF-8", sep=";")
names(suaci)
## [1] "categoria"        "subcategoria"     "concepto"        
## [4] "tipo_prestacion"  "fecha_ingreso"    "domicilio_cgpc"  
## [7] "domicilio_barrio" "lat"              "long"

Una vez que se cargaron en R las bases Radios y SUACI, se crearán tablas intermedias por comunas para localizar su información en un mapa de la Ciudad.

Para emprolijar los datos, se creará una primary key para unir ambas bases y se filtrará la base SUACI por la subcategoría de reclamos “Subtes”:

comunas <- radios %>% 
  mutate (COMUNA = paste("COMUNA", COMUNA)) %>%
  group_by(COMUNA) %>%
  summarise(POBLACION = sum(POBLACION),
            VIVIENDAS = sum(VIVIENDAS),
            HOGARES = sum(HOGARES),
            HOGARES_NBI = sum(HOGARES_NBI),
            AREA_KM2 = sum(AREA_KM2))

suaci_subte <- suaci %>% 
  mutate(COMUNA = as.character(domicilio_cgpc)) %>%
  filter(subcategoria == "SUBTES") %>%
  group_by(COMUNA)

comunas_filt <- comunas %>% 
  left_join(suaci_subte)
## Joining, by = "COMUNA"

Información geográfica

Como el próximo objetivo de esta investigación es localizar los reclamos por filtraciones en las estaciones de Subte, se cargará la información geográfica de líneas y estaciones de Subte, también disponibles en el portal de datos de Gobierno Abierto de GCBA 4:

subte_lineas <- st_read("subte_lineas-1.geojson")
## Reading layer `subte_lineas-1' from data source `G:\Archivos personales\Laboral\CDyPP\subte_lineas-1.geojson' using driver `GeoJSON'
## Simple feature collection with 81 features and 2 fields
## geometry type:  MULTILINESTRING
## dimension:      XY
## bbox:           xmin: -58.48639 ymin: -34.64331 xmax: -58.36993 ymax: -34.55564
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
subte_estaciones <- st_read("subte_estaciones.geojson")
## Reading layer `subte_estaciones' from data source `G:\Archivos personales\Laboral\CDyPP\subte_estaciones.geojson' using driver `GeoJSON'
## Simple feature collection with 87 features and 3 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -58.48639 ymin: -34.64331 xmax: -58.36993 ymax: -34.55564
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

Reclamos por filtraciones

Ahora bien, con la base de reclamos por comuna filtrada por subcategoría “Subtes”, se creará un subconjunto de los datos con la cantidad, por comuna, de reclamos por filtraciones:

filtraciones_en_subte <- comunas_filt %>% 
  filter(concepto == "FILTRACIONES - SUBTE") %>% 
  group_by(COMUNA) %>%
  summarise(RECLAMOS = n())

Por último, se creará un mapa temático que refleje el total de reclamos por filtraciones por comuna:

ggplot() +
  geom_sf(data = radios) +
  geom_sf(data = filtraciones_en_subte, aes(fill = RECLAMOS)) +
  geom_sf(data = subte_lineas, color = "yellow") +
  geom_sf(data = subte_estaciones, color = "black") +
  scale_fill_viridis(option = "magma") +
  labs(title = "Reclamos por filtraciones en el Subte de Buenos Aires",
       subtitle = "Registros de atención al ciudadano. \nCiudad Autónoma de Buenos Aires. Año 2019",
       fill = "Reclamos")

Conclusiones

Durante el primer trimestre del año 2019, se registraron 38 reclamos por filtraciones, de donde se infiere que los usuarios consideran que su presencia afecta la seguridad de las estaciones.

En el mapa presentado anteriormente, se puede observar que casi todas las comunas que son atravesadas por las trazas de la red contaron con reclamos, en especial la Comuna 1 que está compuesta por los barrios del Área Central de la Ciudad, donde está ubicado el Microcentro y se concentran la mayoría de los empleos.

A partir de este gráfico se puede inferir que las estaciones Congreso de Tucumán, Juramento y José Hernández; San Pedrito, Carabobo y Puan; Plaza de los Virreyes, Medalla Milagrosa y Emilio Mitre; Hospitales, Parque Patricios y Caseros, de las líneas D, A, E y H respectivamente, no poseen filtraciones. No obstante, esto no quiere decir que todo el resto de las estaciones de la red sí posean, sino que por lo menos hay una estación por comuna con este inconveniente.

Como conclusión final del trabajo, se puede decir que los usuarios escogen algún medio de transporte alternativo que no se vea afectado ante la probabilidad de precipitaciones.

Referencias bibliográficas

  • Banco Interamericano de Desarrollo (BID) & Subterráneos de Buenos Aires Sociedad del Estado (SBASE). 2015. Plan Estratégico y Técnico para la Expansión de la Red de Subtes de Buenos Aires.
  • Buenos Aires Data. 2019. Datasets. Disponible en: https://data.buenosaires.gob.ar/ Acceso en 31 de mayo de 2019.
  • Dirección General de Estadísticas y Censos (DGEyC). 2018. Anuario Estadístico 2017. Ministerio de Economía y Finanzas. Gobierno de la Ciudad de Buenos Aires. - enelSubte.com. 2018. La línea H ya superó a la línea E en cantidad de pasajeros transportados. Disponible en: https://enelsubte.com/noticias/la-linea-h-ya-supero-a-la-linea-e-en-cantidad-de-pasajeros-transportados/ Acceso en 25 de junio de 2019.
  • Evanoff, M. 2015. Analyze the New York Subway System. Disponible en: https://rpubs.com/mevanoff24/107775. Acceso en 31 de mayo de 2019.
  • Novales, A. 2010. Análisis de Regresión. Departamento de Economía Cuantitativa. Universidad Complutense.
  • Servicio Meteorológico Nacional (SMN). 2019. Datos Horarios, Datos Abiertos. Disponible en: https://www.smn.gob.ar/descarga-de-datos. Acceso en 14 de junio de 2019.
  • TweetDeck. 2019. Tweet content matching @SMN_OCBA from 01 Apr 2019 to 29 Apr 2019. Disponible en: https://tweetdeck.twitter.com/ Acceso en 14 de junio de 2019.

  1. Economista (FCE-UBA) y casi Magister en Generación y Análisis de Información Estadística (UNTREF). Trabajo como Analista de Demanda en la Gerencia de Planeamiento de Subterráneos de Buenos Aires S.E. (SBASE), aunque este informe lo realizo como condición para la aprobar el Curso de Ciencia de Datos y Políticas Públicas. Contacto:

  2. La base Molinetes está disponible en el portal de datos de Gobierno Abierto del GCBA con información desglosada por molinete cada 15 minutos. Para fines de este análisis, se generó una nueva base con datos horarios desglosados por estación.

  3. No se cuenta con información proveniente de los validadores del Premetro.

  4. Hasta el momento, en el portal no se encuentran las trazas definitivas, las tres estaciones inauguradas el día 03/06/19 de la Línea E, ni la estación Facultad de Derecho de la Línea H inaugurada el 17/05/18. Los archivos con formato GEOJSON fueron editados en la plataforma http://geojson.io/ para incorporar la estación Facultad de Derecho y la traza definitiva de la Línea H. No se incorporaron las estaciones Correo Central, Catalinas y Retiro de la Línea E debido a que la base SUACI cuenta con información hasta el 31/03/19 inclusive.