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.

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 R² reflejará la bondad del ajuste de un modelo a la variable que pretende explicar.

Por otro lado, el R² 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 R² y R² 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 R² 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"))

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 R² 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 ⁴:

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

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.