1 Introducción

Partiendo de los datos de movimientos de BiciMAD realizamos dos modelos, uno para predicir la demanda por estación y hora del día y otro para predecir la demanda por estación y día en los próximos meses.

En esta ocasión, no hemos considerado los datos de marzo a junio de 2020 ya que la crisis del Coronavirus ha cambiado mucho el comportamiento normal de los usuarios de BiciMAD. En la presentación de la serie hemos visto como durante la segunda quincena de marzo y hasta finales de abril, no ha habido uso de BiciMAD. En mayo y junio todavía existían bastantes restricciones horarias en cuanto a la movilidad de la población. Todas estas circunstancias nos han llevado a prescindir de ese periodo temporal.

2 Preparación de los datos

# Librerías utilizadas durante todo el proyecto
library(ndjson)
library(stringr)
library(chron)
library(dplyr)
library(rjson)
library(geosphere)
library(caret)
library(ranger)
library(data.table)
library(ggplot2)
library(lubridate)
library(Hmisc)
library(utils)
# Directorio de trabajo
setwd("C:/Users/rodsp/Desktop/DataScience/DATASETS/TFM/MOVIMIENTOS TOTAL")

# Leer datos
movimientos <- as.data.frame(fread("movimientos2019-jun2020.csv"))
summary(movimientos)
##    _id.$oid         user_day_code        user_type        ageRange    
##  Length:5062186     Length:5062186     Min.   :1.000   Min.   :1.000  
##  Class :character   Class :character   1st Qu.:1.000   1st Qu.:4.000  
##  Mode  :character   Mode  :character   Median :1.000   Median :4.000  
##                                        Mean   :1.148   Mean   :4.239  
##                                        3rd Qu.:1.000   3rd Qu.:5.000  
##                                        Max.   :3.000   Max.   :6.000  
##   idplug_base    idunplug_base   idplug_station   idunplug_station
##  Min.   : 1.00   Min.   : 1.00   Min.   :  1.00   Min.   :  1.0   
##  1st Qu.: 5.00   1st Qu.: 5.00   1st Qu.: 51.00   1st Qu.: 51.0   
##  Median :12.00   Median :12.00   Median : 96.00   Median : 96.0   
##  Mean   :12.13   Mean   :12.11   Mean   : 99.84   Mean   : 99.9   
##  3rd Qu.:19.00   3rd Qu.:19.00   3rd Qu.:149.00   3rd Qu.:149.0   
##  Max.   :30.00   Max.   :30.00   Max.   :219.00   Max.   :219.0   
##  unplug_hourTime     travel_time         zip_code            fe_year    
##  Length:5062186     Min.   :       0   Length:5062186     Min.   :2019  
##  Class :character   1st Qu.:     462   Class :character   1st Qu.:2019  
##  Mode  :character   Median :     700   Mode  :character   Median :2019  
##                     Mean   :    1491                      Mean   :2019  
##                     3rd Qu.:    1071                      3rd Qu.:2019  
##                     Max.   :19272889                      Max.   :2020  
##     fe_month         fe_day         fe_wday         fe_hour     
##  Min.   : 1.00   Min.   : 1.00   Min.   :1.000   Min.   : 0.00  
##  1st Qu.: 3.00   1st Qu.: 8.00   1st Qu.:3.000   1st Qu.: 9.00  
##  Median : 5.00   Median :14.00   Median :4.000   Median :14.00  
##  Mean   : 5.83   Mean   :15.08   Mean   :4.066   Mean   :13.67  
##  3rd Qu.: 9.00   3rd Qu.:23.00   3rd Qu.:6.000   3rd Qu.:18.00  
##  Max.   :12.00   Max.   :31.00   Max.   :7.000   Max.   :23.00  
##   fe_season        
##  Length:5062186    
##  Class :character  
##  Mode  :character  
##                    
##                    
## 
# Nos quedamos solo con los datos hasta febrero 2020
# Conviertimos a fecha_hms
movimientos$unplug_hourTime <- ymd_hms(movimientos$unplug_hourTime)
movimientos$fecha <- as.Date(movimientos$unplug_hourTime, format="%d/%m/%Y")
movimientos <- subset(movimientos, movimientos$fecha < '2020-03-01') 

# Eliminamos la variable unplug_hourTime ya que da la misma info que fecha
movimientos$unplug_hourTime <- NULL

2.1 Exploración y limpieza de datos de movimientos

A continuación, exploramos y limpiamos las variables de este dataset creado con los datos de movimientos de BiciMAD desde enero de 2019 a febrero de 2020.

2.1.1 ids

Las variables de id no son relevantes ya que son códigos generamos por BiciMAD para cada viaje, pero no nos proporciona información de los usuarios ni de las estaciones, por lo que las eliminamos.

# Cambiamos el nombre a _id.$oid ya que da problemas al tratarla
names(movimientos)[1] <- 'id1'

# Eliminamos 
movimientos$id1 <- NULL
movimientos$user_day_code <- NULL

2.1.2 user_type

Para nuestra predicción hemos optado por eliminar los usuarios tipo 3 ya que no son viajes propiamente dichos sino movimientos que realizan los empleados de BiciMAD para el mantimiento de las bicis.

movimientos_sint3 <- subset(movimientos, movimientos$user_type != 3) 

2.1.3 ageRange ,zip_code y idunplug_base

Estas variables relativas al usuario tampoco las consideramos relevantes para el modelo, pues lo que estamos prediciendo son número de viajes por estación a cada hora del día sin saber previamente la información de los usuarios para cada uno de esos viajes.

movimientos_sint3$ageRange <- NULL
movimientos_sint3$zip_code <- NULL
movimientos_sint3$idunplug_base<-NULL

2.1.4 fe_season

Recategorizamos esta variable como factor y cambiamos por su frecuencia.

movimientos_sint3$fe_season <- as.factor(movimientos_sint3$fe_season)

# Cambiar por su frecuencia
movimientos_sint3 <- as.data.table(movimientos_sint3)
movimientos_sint3[ , fe2_season := .N , by = .(fe_season)]   

# Eliminamos la original
movimientos_sint3$fe_season <- NULL

2.1.5 travel_time

Consideramos que los viajes de menos de 60 sg son viajes erroneos, por ejemplo, porque la bici no funciona correctamente y se vuelve a dejar en la misma estación y los viajes superiores a 7200 sg, por ejemplo, porque no se ha contabilizado bien el tiempo (suele pasar y tenemos experiencia en ello). Se crea un nuevo dataframe llamado movimientos_sint3_time.

Por lo tanto, ahora tenemos viajes con los siguientes valores:

  • El viaje más corto es de 61 sg.
  • El viaje más largo es de 7199 sg, es decir, unas 2 horas
  • La duración media de los viajes es de 841 sg, es decir, 14 minutos

Después de filtrar la duración de viajes con la que trabajaremos, eliminamos la variable travel time de nuestro data frame, pues es información de usuarios.

# Eliminación de viajes cortos y largos
movimientos_sint3_time <- subset(movimientos_sint3, movimientos_sint3$travel_time > 60) 
movimientos_sint3_time <- subset(movimientos_sint3_time, movimientos_sint3_time$travel_time < 7200) 

# Visualización
summary(movimientos_sint3_time$travel_time)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    61.0   454.0   667.0   841.5   971.0  7199.0
colnames(movimientos_sint3_time)
##  [1] "user_type"        "idplug_base"      "idplug_station"   "idunplug_station"
##  [5] "travel_time"      "fe_year"          "fe_month"         "fe_day"          
##  [9] "fe_wday"          "fe_hour"          "fecha"            "fe2_season"
movimientos_sint3_time$travel_time<-NULL  

2.1.6 Variables destino

Como nuestro objetivo es predecir el número de viajes que salen de la estación de origen, eliminamos las variables de estación de destino.

movimientos_sint3_time$idplug_base <- NULL
movimientos_sint3_time$idplug_station <- NULL
movimientos_sint3_time$user_type<-NULL

2.2 Datos de clima

El servicio BiciMAD funciona en el centro de la ciudad de Madrid. Dado que es un servicio al aire libre, consideramos relevante estudiar si la metereologia del día influye en su uso (parece que sí según los datos de meses y estaciones vistos en el EDA).

Descargamos estos datos de la Agencia Estatal de Metereología (Aemet) que también cuenta con una base de datos pública a través de su API. Hemos descargado los datos de clima de los meses que nos interesan para la estación de Retiro ya que es la más centrica y acorde al area de influencia de BiciMAD.

Para la lectura de datos utilizamos la función read.csv2 ya que los decimales los configura con una coma en lugar de punto, con esta función asume que esas comas son puntos y trata la variable como numérica (sino lo hace como character).

# Directorio de trabajo
setwd("C:/Users/rodsp/Desktop/DataScience/DATASETS/TFM/MOVIMIENTOS TOTAL")

# Lectura de datos
clima <- read.csv2("aemet_ene2019-feb2020.csv") 

head(clima)
str(clima)
## 'data.frame':    425 obs. of  19 variables:
##  $ fecha      : chr  "01/01/2019" "02/01/2019" "03/01/2019" "04/01/2019" ...
##  $ indicativo : int  3195 3195 3195 3195 3195 3195 3195 3195 3195 3195 ...
##  $ nombre     : chr  "MADRID, RETIRO" "MADRID, RETIRO" "MADRID, RETIRO" "MADRID, RETIRO" ...
##  $ provincia  : chr  "MADRID" "MADRID" "MADRID" "MADRID" ...
##  $ altitud    : int  667 667 667 667 667 667 667 667 667 667 ...
##  $ tmed       : num  6.6 6 7.3 5.9 5.4 6.6 7.1 7.1 7.6 4.7 ...
##  $ prec       : chr  "0" "0" "0" "0" ...
##  $ tmin       : num  1.4 0.4 2.8 0.4 -0.1 -0.5 0.5 1.2 2 0.6 ...
##  $ horatmin   : chr  "7:30" "6:30" "7:00" "8:05" ...
##  $ tmax       : num  11.7 11.7 11.8 11.4 10.8 13.8 13.7 13 13.2 8.8 ...
##  $ horatmax   : chr  "15:10" "15:15" "15:20" "15:10" ...
##  $ dir        : int  6 6 5 6 9 2 8 34 3 3 ...
##  $ velmedia   : num  0.3 0.6 0.6 0.3 0.6 1.1 0.6 0.8 1.4 2.8 ...
##  $ racha      : num  5.3 7.8 7.8 3.6 3.1 3.3 5 3.6 9.7 10 ...
##  $ horaracha  : chr  "0:00" "23:20" "0:10" "5:40" ...
##  $ presMax    : num  954 952 951 953 955 ...
##  $ horaPresMax: chr  "9" "0" "10" "Varias" ...
##  $ presMin    : num  952 948 949 950 952 ...
##  $ horaPresMin: chr  "5" "16" "5" "5" ...

2.2.1 Limpieza datos clima

Primero, quitamos las variables no relevantes para el ejercicio: el indicativo, el nombre de la estación, la provincia y la altitud ya que son siempre iguales para la estación Retiro.

Dirección y racha del viento tienen valores nulos que imputamos por aleatorio.

Para poder unir esta dataframe al anterior de movimiento lo hacemos con la variable fecha, para ello, debe tener el mismo formato en ambos dataframes.

# Quitamos ciudad, provincia etc y las de hora
clima_final <- clima[ , -c(2:5, 9, 11, 15, 16, 17, 19)]

# Imputamos dir y racha
clima_final$dir <- impute(clima_final$dir, "random")
clima_final$racha <- impute(clima_final$racha, "random")

# Cambiamos a fecha la fecha del clima
clima_final$fecha <- as.Date(clima_final$fecha, format="%d/%m/%Y")

2.2.2 Feature engineering con variables de clima

  • Con las precipitaciones creamos una nueva variable binaria: - 1: llovió - 0: no llovió

  • La dirección del viento la modificamos a intervalos abiertos por la derecha según Rosa del viento (según este método se podrían hacer más o menos intervalos, hemos cogido el intermedio): - Norte: 0-45º - Noreste: 45-90º - Este: 90-135º - Sureste: 135-180º - Sur: 180-225º - Suroeste: 225-270º - Oeste: 270-335º - Noroeste: 335-360º En nuestro caso solo hay valores hasta el 99 por lo que solo habrá viento Norte, Noreste y Este que hemos llamado respectivamente 1, 2 y 3

  • La velocidad del viento la modificamos en tramos (está medido en m/s): - Flojo: 0 - 5.4 que serán los valores 1 - Fresco: 5.5 - 17.1 que serán los valores 2 - Temporal > 17.2 que serán los valores 3 (aunque no habrá ninguno)

Estas ideas sobre la dirección y velocidad del viento han sido extraidas de esta página: https://es.windfinder.com/wind/windspeed.htm

clima_final$fe_prec_bin <- factor(replace(clima_final$prec, which(clima_final$prec != 0), 1))

clima_final$fe_dir <- cut(clima_final$dir, breaks = c(-1, 45, 90, 135), labels = c("1", "2", "3")) # no ponemos más intervalos porque hemos visto en el summary que el máximo es 99

clima_final$fe_velmedia <- cut(clima_final$velmedia, breaks = c(-1, 5.4, 17.1, 50), labels = c("1", "2", "3"))

2.2.3 Unión al dataframe de movimientos

Estas variables las unimos al dataframe anterior de movimientos a través de la variable fecha, creando un nuevo dataframe que llamaremos mov_con_clima.

#Unión de ambos df
mov_con_clima <-merge(movimientos_sint3_time, clima_final, by = "fecha")

head(mov_con_clima)

2.3 Datos de las estaciones

Cargamos los datos de las estaciones que también tratamos anteriormente, en este caso, utilizaremos las variables de localización de la estación, su id para poder unirlo al otro dataframe y el número de bases con las que cuenta.

Como estamos prediciendo el número de viajes desde la estación origen, consideramos estos valores para ella (idunplug_station).

Este conjunto de datos nos permite añadir más información a nuestro data frame, especialmente localizar en qué parte de la ciudad se situa cada estación.

2.3.1 Lectura y limpieza de datos

# Directorio de trabajo
setwd("C:/Users/rodsp/Desktop/DataScience/DATASETS/TFM/MOVIMIENTOS TOTAL")

# Leer datos
estaciones <- as.data.frame(fread("estaciones.csv"))
summary(estaciones)
##     activate     name           reservations_count     light      
##  Min.   :1   Length:214         Min.   :0          Min.   :0.000  
##  1st Qu.:1   Class :character   1st Qu.:0          1st Qu.:2.000  
##  Median :1   Mode  :character   Median :0          Median :2.000  
##  Mean   :1                      Mean   :0          Mean   :1.925  
##  3rd Qu.:1                      3rd Qu.:0          3rd Qu.:2.000  
##  Max.   :1                      Max.   :0          Max.   :3.000  
##   total_bases     free_bases       number            longitude     
##  Min.   :12.0   Min.   : 0.00   Length:214         Min.   :-3.725  
##  1st Qu.:24.0   1st Qu.:10.00   Class :character   1st Qu.:-3.704  
##  Median :24.0   Median :12.00   Mode  :character   Median :-3.695  
##  Mean   :23.9   Mean   :11.63                      Mean   :-3.694  
##  3rd Qu.:24.0   3rd Qu.:13.00                      3rd Qu.:-3.683  
##  Max.   :30.0   Max.   :24.00                      Max.   :-3.656  
##   no_available        address             latitude       dock_bikes   
##  Min.   :0.000000   Length:214         Min.   :40.39   Min.   : 0.00  
##  1st Qu.:0.000000   Class :character   1st Qu.:40.41   1st Qu.:10.00  
##  Median :0.000000   Mode  :character   Median :40.42   Median :12.00  
##  Mean   :0.004673                      Mean   :40.42   Mean   :11.37  
##  3rd Qu.:0.000000                      3rd Qu.:40.43   3rd Qu.:13.00  
##  Max.   :1.000000                      Max.   :40.47   Max.   :19.00  
##        id        
##  Min.   :  1.00  
##  1st Qu.: 57.25  
##  Median :112.50  
##  Mean   :111.60  
##  3rd Qu.:165.75  
##  Max.   :219.00
# Nos quedamos con las variables que nos interesan y creamos el df
estaciones_origen <- estaciones[ , c(5,8,11,13)]
# Cambiamos los nombres
names(estaciones_origen)[4] <- "idunplug_station"
names(estaciones_origen)[1] <- "total_bases_origen"
names(estaciones_origen)[2] <- "longitude_origen"
names(estaciones_origen)[3] <- "latitude_origen"

2.3.2 Feature engineering con datos estaciones

Habitualmente encontramos las coordenadas geográficas representadas por Latitud y Longitud separadas por una coma, por ejemplo: 40.4168961º, -3.7024255º. Creamos una nueva variable que así lo represente. Después, las cambiamos por su frecuencia para que entren correctamente en el modelo.

# Creamos coordenadas
estaciones_origen$coordenadas_origen <- paste(estaciones_origen$latitude_origen, "º", ", ", estaciones_origen$longitude_origen, "º", sep = "")

# Cambiamos por frecuencia. Antes nos aseguramos que es data.table
estaciones_origen <- as.data.table(estaciones_origen)

estaciones_origen[ , fe_coordenadas_origen := .N , by = .(coordenadas_origen)] 

# Eliminamos original
estaciones_origen$coordenadas_origen <- NULL

También creamos una variable que situe cada estación según lo lejos que esté del centro. Tomamos como referencia el Kilometro cero de la Puerta del Sol y calculamos la distancia en km que hay desde ese punto a cada estación.

estaciones_origen <- estaciones_origen%>% rowwise() %>% 
  mutate(dist_centro_origen = distHaversine(c(latitude_origen, longitude_origen), c(40.41664,-3.7059988))/1000)

head(estaciones_origen)

Y unimos al dataframe con el que veníamos trabajando haciendo un merge tomando como referencia la estación de desanclaje.

mov_con_clima_estaciones<-merge(as.data.frame(mov_con_clima), as.data.frame(estaciones_origen), by= "idunplug_station")

2.4 Creación variable target

Nuestro objetivo es predecir el número de viajes que se realizarán cada hora del día por estación. Para ello, creamos nuestra variable objetivo o target con el número de viajes que se realizan en cada estación en cada hora del día y le unimos el df que hemos trabajado.

# Contamos los viajes por estación por día y hora
target <- mov_con_clima_estaciones%>%count(idunplug_station, fecha, fe_hour, sort=TRUE)

# Nos quedamos con las filas únicas.
mov_con_clima_estaciones<-distinct(mov_con_clima_estaciones)


# Y unimos el df original 

mov_con_clima_estaciones_target <-merge(target, mov_con_clima_estaciones, by = c("idunplug_station", "fecha", "fe_hour"))

# Cambiamos nombre a la target
names(mov_con_clima_estaciones_target)[4] <- "target"

# Visualización df final con la target
head(mov_con_clima_estaciones_target)
#summary(mov_con_clima_estaciones_target)
sum(is.na(mov_con_clima_estaciones_target))
## [1] 0
str(mov_con_clima_estaciones_target)
## 'data.frame':    1203308 obs. of  25 variables:
##  $ idunplug_station     : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ fecha                : Date, format: "2019-01-01" "2019-01-01" ...
##  $ fe_hour              : int  12 15 16 17 18 19 20 22 23 9 ...
##  $ target               : int  1 1 2 2 9 1 1 1 2 1 ...
##  $ fe_year              : int  2019 2019 2019 2019 2019 2019 2019 2019 2019 2019 ...
##  $ fe_month             : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ fe_day               : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ fe_wday              : int  3 3 3 3 3 3 3 3 3 3 ...
##  $ fe2_season           : int  1582753 1582753 1582753 1582753 1582753 1582753 1582753 1582753 1582753 1582753 ...
##  $ tmed                 : num  6.6 6.6 6.6 6.6 6.6 6.6 6.6 6.6 6.6 6.6 ...
##  $ prec                 : chr  "0" "0" "0" "0" ...
##  $ tmin                 : num  1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 ...
##  $ tmax                 : num  11.7 11.7 11.7 11.7 11.7 11.7 11.7 11.7 11.7 11.7 ...
##  $ dir                  : int  6 6 6 6 6 6 6 6 6 6 ...
##  $ velmedia             : num  0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 ...
##  $ racha                : num  5.3 5.3 5.3 5.3 5.3 5.3 5.3 5.3 5.3 5.3 ...
##  $ presMin              : num  952 952 952 952 952 ...
##  $ fe_prec_bin          : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ fe_dir               : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ fe_velmedia          : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ total_bases_origen   : int  30 30 30 30 30 30 30 30 30 30 ...
##  $ longitude_origen     : num  -3.7 -3.7 -3.7 -3.7 -3.7 ...
##  $ latitude_origen      : num  40.4 40.4 40.4 40.4 40.4 ...
##  $ fe_coordenadas_origen: int  1 1 1 1 1 1 1 1 1 1 ...
##  $ dist_centro_origen   : num  0.473 0.473 0.473 0.473 0.473 ...

A continuación hacemos el mismo proceso, pero para crear una nueva variable target que usaremos en el segundo modelo. En este caso el número de viajes por estación y día.

target2<-mov_con_clima_estaciones_target%>%group_by(idunplug_station,fecha)%>%summarise(n=sum(target)) # Agrupamos por fecha.

mov_con_clima_estaciones2<-mov_con_clima_estaciones[,-7] # Quitamos la variable "fe_hour"
mov_con_clima_estaciones2<-distinct(mov_con_clima_estaciones2) # Filtramos solo por las filas únicas.

mov_dia <-merge(target2, mov_con_clima_estaciones2, by = c("idunplug_station", "fecha")) # Juntamos la variable target con el data frame.
names(mov_dia)[3]<-"target"

3 Modelo

Ahora que tenemos nuestro data frame pulido y con nuevas variables que consideramos importantes, procedemos a hacer un modelo predictivo del número de viajes según la estación, el día y la hora.

3.1 Modelo por hora.

3.1.1 Separación train y test

Hacemos una partición de 70% train y 30% test.

set.seed(1234)
trainIndex <- createDataPartition(mov_con_clima_estaciones_target$target, p=0.7, list=FALSE)
data_train <- mov_con_clima_estaciones_target[trainIndex,]
data_test  <- mov_con_clima_estaciones_target[-trainIndex,]

3.1.2 Modelo random forest.

Vamos a hacer un modelo de random forest, utilizando la implementación de la librería ranger.

modelo_ranger1 <- ranger( 
                      formula    = target~., 
                      data       = data_train,
                      num.trees= 100,
                      importance = 'impurity',          
                      verbose    = TRUE)
## Growing trees.. Progress: 22%. Estimated remaining time: 1 minute, 49 seconds.
## Growing trees.. Progress: 47%. Estimated remaining time: 1 minute, 9 seconds.
## Growing trees.. Progress: 71%. Estimated remaining time: 38 seconds.
## Growing trees.. Progress: 96%. Estimated remaining time: 5 seconds.
modelo_ranger1
## Ranger result
## 
## Call:
##  ranger(formula = target ~ ., data = data_train, num.trees = 100,      importance = "impurity", verbose = TRUE) 
## 
## Type:                             Regression 
## Number of trees:                  100 
## Sample size:                      842317 
## Number of independent variables:  24 
## Mtry:                             4 
## Target node size:                 5 
## Variable importance mode:         impurity 
## Splitrule:                        variance 
## OOB prediction error (MSE):       4.921925 
## R squared (OOB):                  0.4466778

Este modelo ranger de tipo regresión, con 100 árboles, nos da un R squared (OOB) de 0.44. Pero este resultado es sobre el conjunto train. Debemos probar el modelo con el conjunto test, para comprobar la precisión del modelo.

Analizando la importancia que le da el modelo a las variables, destacan la hora y la estación de origen. Las variables de geolocalización también tienen bastante relevancia, así como las climáticas.

rang_ranger1 <- as.data.frame(importance(modelo_ranger1))
rang_ranger1$vars <- rownames(rang_ranger1)
names(rang_ranger1)[1]<- c('Importance')
imp_var<-rang_ranger1 %>%
    arrange(desc(Importance))

imp_var_top<-head(imp_var,15)

ggplot(imp_var_top, aes(x=reorder(vars,-Importance), y=Importance))+
  geom_bar(stat="identity",color="darkgreen", fill="lightgreen")+
  theme(axis.text.x=element_text(size=12, angle=45, color="blue",face="bold.italic"))+
  labs(x = "Variables",y = "Importancia")

Predicción sobre test:

pred_mod1 <- predict(modelo_ranger1 , data_test)
# Coeficiente de determinación de la predicción TEST
R2 <- 1 - (sum((data_test$target-pred_mod1$predictions)^2)/sum((data_test$target-mean(data_test$target))^2))

# Coeficientee de determinación ajustado.
n<-length(pred_mod1$predictions)
k<-pred_mod1$num.independent.variables

R2A<- 1-((n-1)/(n-k-1))*(1-R2)

R2
## [1] 0.4514743
R2A
## [1] 0.4514378

El R2 ajustado en test es 0.45. Un resultado muy parecido al R squared OOB del modelo. No tenemos un R2 muy alto, pero sí un modelo que nos complementará la información que obtengamos en el siguiente.

Este modelo nos permitirá adelantarnos a la demanda de las estaciones de Bicimad para cada hora y día del año. Sabiendo cuales van a ser los picos de demanda durante el día, se optimizarían las estaciones ya que los usuarios nunca encontrarían una estación vacía.

3.2 Modelo por día.

En este caso, la partición la hacemos sobre nuestro data frame de número de viajes por día.

set.seed(1234)
trainIndex2 <- createDataPartition(mov_dia$target, p=0.7, list=FALSE)
data_train2 <- mov_dia[trainIndex2,]
data_test2  <- mov_dia[-trainIndex2,]

Hacemos un modelo random forest (ranger), pero en este caso con 200 árboles. Tenemos menos datos por lo que el esfuerzo computacional es menor.

modelo_ranger2 <- ranger( 
                      formula    = target~., 
                      data       = data_train2,
                      num.trees= 200,
                      importance = 'impurity',          
                      verbose    = TRUE)
modelo_ranger2
## Ranger result
## 
## Call:
##  ranger(formula = target ~ ., data = data_train2, num.trees = 200,      importance = "impurity", verbose = TRUE) 
## 
## Type:                             Regression 
## Number of trees:                  200 
## Sample size:                      52842 
## Number of independent variables:  23 
## Mtry:                             4 
## Target node size:                 5 
## Variable importance mode:         impurity 
## Splitrule:                        variance 
## OOB prediction error (MSE):       240.8185 
## R squared (OOB):                  0.8151926

Predicción sobre test:

pred_mod2 <- predict(modelo_ranger2 , data_test2)

Vemos que la importancia de las variables no es muy diferente al modelo anterior excepto por la mayor relevancia de la longitud y la distancia de las estaciones al centro. Tampoco está la hora de desanclaje porque no entraba en el modelo (predicción por día).

imp2<-as.data.frame(modelo_ranger2$variable.importance)

imp2<-cbind(rownames(imp2), imp2$`modelo_ranger2$variable.importance`)
imp2<-as.data.frame(imp2)
imp2$V2<-as.numeric(imp2$V2)
imp2<-imp2[order(-imp2$V2),]
colnames(imp2)[1:2]<-c("Variable", "Importancia")
imp2
# Coeficiente de determinación de la predicción TEST
R2 <- 1 - (sum((data_test2$target-pred_mod2$predictions)^2)/sum((data_test2$target-mean(data_test2$target))^2))

# Coeficientee de determinación ajustado.
n<-length(pred_mod2$predictions)
k<-pred_mod2$num.independent.variables

R2A<- 1-((n-1)/(n-k-1))*(1-R2)

R2
## [1] 0.8150124
R2A
## [1] 0.8148243

El coeficiente de determinación del segundo modelo es muy similar al R squared OOB, 0,81. Tenemos una buena medida de ajuste.

Con este modelo no podemos hacer una predicción tan detallada de la demanda como hicimos en el primero (demanda por hora), en cambio su R2 es mucho más alto. Este modelo permitiría a las empresas reubicar las bicicletas para que haya siempre suficiente stock en las estaciones, y así satisfacer la demanda del día.