Space-Time Prediction of Bike Share Demand in NYC

Introduction

Bike share is a dock-based sharing system, which is also a sustainable and convenient transportation option. The system can provide access to bikes for short-term use, ranging from several minutes to hours. However, the most difficult operational problems is the need to re-balance bikes across the network. As bike share is not useful if a dock has no bikes to pick up, nor if there are no opening docking spaces to deposit. Therefore, the re-balancing is the practice of anticipating bike share demand for all docks at all times and manually redistributing bikes to ensure a bike or a docking space is available when needed.

In this assignment, I’ll dive into the bike share system in New York City (Citi Bike) and forecast space/time demand for bike share pickups. My focused area and time frame is Manhattan, May 2023. By developing an algorithm to predict the trip starts across time and space, I try to en-vision and inform a bike re-balancing plan to ensure supply and enhance efficiency.

Specifically, several mechanisms will be taken in light of expected outcomes. For example, there will be trucks to move bikes between different locations. Also, some incentives will be offered if riders manage to move a bike from place to place.

Data Wrangling

In this assignment,I forecast space/time demand for bike share pickups in New York City and designed an algorithm to inform a bike re-balancing plan. The data I use here is obtained from Citi Bike Trip Data that contains the latest monthly record of Citi Bike usage in New York City.

My model will take a low resolution approach, reducing millions of NYC ride share trips from May through June, 2023, into a 15% subsample and aggregating to hourly intervals and a subset of New York City Census tracts.

Trip Data

Census Data

Several census geography and variables are loaded and wrangled to test generalizability later. Among them, there are Population, Race, Means of Transport, Mean Commute Time, % of taking public transit, etc.

Weather Data

Furthermore, I have incorporated weather data from JFK Airport into my analysis. To be specific, I imported hourly data related to temperature, wind speed, and precipitation. I have also created visualizations to illustrate the trends in temperature and precipitation throughout the duration of my study period.

weather.Panel <- 
  riem_measures(station = "JFK", date_start = "2023-05-01", date_end = "2023-05-31") %>%
  dplyr::select(valid, tmpf, p01i, sknt)%>%
  replace(is.na(.), 0) %>%
    mutate(interval60 = ymd_h(substr(valid,1,13))) %>%
    mutate(week = week(interval60),
           dotw = wday(interval60, label=TRUE)) %>%
    group_by(interval60) %>%
    summarize(Temperature = max(tmpf),
              Precipitation = sum(p01i),
              Wind_Speed = max(sknt)) %>%
    mutate(Temperature = ifelse(Temperature == 0, 42, Temperature))

glimpse(weather.Panel)

## Rows: 721
## Columns: 4
## $ interval60    <dttm> 2023-05-01 00:00:00, 2023-05-01 01:00:00, 2023-05-01 02…
## $ Temperature   <dbl> 56, 54, 52, 50, 50, 50, 51, 55, 57, 59, 60, 59, 58, 59, …
## $ Precipitation <dbl> 3e-02, 1e-04, 0e+00, 0e+00, 0e+00, 0e+00, 0e+00, 0e+00, …
## $ Wind_Speed    <dbl> 14, 19, 21, 16, 19, 19, 18, 17, 15, 18, 15, 18, 20, 17, …

grid.arrange(
  ggplot(weather.Panel, aes(interval60,Precipitation)) + geom_line(color="#FF759F") + 
  labs(title="Percipitation", x="Hour", y="Perecipitation") + plotTheme(),
  ggplot(weather.Panel, aes(interval60,Wind_Speed)) + geom_line(color="#FF759F") + 
    labs(title="Wind Speed", x="Hour", y="Wind Speed") + plotTheme(),
  ggplot(weather.Panel, aes(interval60,Temperature)) + geom_line(color="#FF759F") + 
    labs(title="Temperature", x="Hour", y="Temperature") + plotTheme(),
  top="Weather Data - NYC JFK - May, 2023")

Amenity Data

In addition, I incorporated amenity data into my analysis, as it has the potential to significantly impact the bike re-balancing process. For instance, the proximity of various points of interest such as landmarks, markets, bus stations, and colleges can influence the demand for shared bikes. Therefore, I included these amenity features as variables in my analysis. To further refine the analysis, I employed the k-nearest neighbor (k=1) method to calculate the nearest landmark, market, bus station, and college for each bike station. This helps to account for the influence of these amenities on bike station usage and re-balancing strategies.

Space-Time Panel

I created a data frame panel which has each unique space/time observations. Then the trip counts were summarized by station for each time interval. Also, the socio-economic information and latitude/longitude information were also kept for future data integration and analysis.

length(unique(dat_census$interval60)) * length(unique(dat_census$start_station_id))

## [1] 494760

study.panel <- 
  expand.grid(interval60=unique(dat_census$interval60), 
              start_station_id = unique(dat_census$start_station_id)) %>%
  left_join(., dat_census %>%
              select(start_station_id, start_station_name, Origin.Tract, start_lng, start_lat )%>%
              distinct() %>%
              group_by(start_station_id) %>%
              slice(1)) %>%
  left_join(.,amenity_nn %>%
              select(start_station_id, Landmarks.nn, Markets.nn, Bus_Stations.nn, College.nn)%>%
              distinct() %>%
              group_by(start_station_id) %>%
              slice(1))

## Joining with `by = join_by(start_station_id)`
## Joining with `by = join_by(start_station_id)`

nrow(study.panel)      

## [1] 494760

Time Lags

Moreover, time lag variables about the demand during a given time period are also added. by evaluating the correlations in the lags, we can tell that some of the lags are pretty strong. For example, there’s a Pearson’s R of 0.87 and 0.86 for the lagHour and lag1day respectively, indicating that the demand for an hour ago and a day ago is strongly correlated with the demand now. However, the 12-hour lag exhibits an opposite relationship in terms of demand.

## # A tibble: 6 × 2
##   Variable   correlation
##   <fct>            <dbl>
## 1 lagHour           0.87
## 2 lag2Hours         0.64
## 3 lag3Hours         0.41
## 4 lag4Hours         0.21
## 5 lag12Hours       -0.49
## 6 lag1day           0.85

Exploratory Analysis

Serial Autocorrelation

I begin by examining the time and frequency components of the data. First, from the overall time pattern - there is clearly a daily periodicity and there are lull periods on weekends. Also, the weekend near the 22th of May doesn’t have the same dip in activity due to the extreme weather.

By examining the distribution of trip volume by station for different times of the day, it’s clear that there are a few high volume periods but mostly low volume.

I’ve also analyzed the daily trends in ridership, differentiating between days of the week and distinguishing between weekends and weekdays. The temporal patterns are quite apparent, with weekdays exhibiting the highest trip counts. Furthermore, among weekdays, the peak in trip counts aligns with commuting hours, indicating a strong correlation between ridership and typical commuting times.

The figure below plots out the aggregate number of trips as a function of date. The vertical lines indicate Mondays. From it, the weekly pattern is obvious that there are peaks and troughs within a week and a day.

## `summarise()` has grouped output by 'Legend'. You can override using the
## `.groups` argument.

Further, I plot the ridershare trip count as a function of spatial lags to check the correlation. The analysis reveals a robust correlation between trip initiations and lag features. However, this correlation gradually diminishes as the lag period increases within a single day. Notably, a lag of 12 hours exhibits no statistically significant relationship with trip initiations. However, a lag of 1 day has.

## `geom_smooth()` using formula = 'y ~ x'

Spatial Correlation

Then it comes to spatial correlation, the bike-share ridership also shows spatial autocorrelation. The maps below further show the sum of bike share trips by station and by days of the week, indicating that the majority of trips start in midtown & lower Manhattan.

## `summarise()` has grouped output by 'week', 'Origin.Tract', 'start_lng'. You
## can override using the `.groups` argument.

## `summarise()` has grouped output by 'dotw', 'Origin.Tract', 'start_lng'. You
## can override using the `.groups` argument.

Space/time correlation

Again, the bike share ridership also exhibits strong space/time correlation. To show, I pick up one day in NYC in May, 2023 to visualize the relationship with 15-minute intervals via an animation.

Weather

As for the effect of weather, obviously ridership also varies with precipitation and temperature.

The average number of trips per week seems to trend upward as the temperature increases. This trend is pretty consistent across all panels.

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 24 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 24 rows containing missing values (`geom_point()`).

Modeling and Validation

In this part, I split the data into a training and a test set and created five linear models using the lm funtion. The first models include only temporal controls, but the later ones contain all of our lag information or amenities.

Five linear regressions are further estimated on bike-share train data, each with different fixed effects: 1. reg 1 focuses on just time, including hour fixed effects, day of the week, and Temperature. 2. reg 2 further adds space effects for the across-station differences. 3. reg 3 combines the time and space effects, and also adds more weather effects, such as precipitation. 4. reg 4 takes time lag features into consideration. 5. reg 5 further adds more amenities effects, such as landmarks, markets, colleges, and bus stations.

##Predict for test data

ride.Test.weekNest <- 
  ride.Test %>%
  nest(-week) 

ride.Test.weekNest

## # A tibble: 2 × 2
##    week data                   
##   <dbl> <list>                 
## 1    18 <tibble [95,760 × 34]> 
## 2    19 <tibble [111,720 × 34]>

## # A tibble: 10 × 8
##     week data     Regression      Prediction Observed Absolute_Error   MAE sd_AE
##    <dbl> <list>   <chr>           <list>     <list>   <list>         <dbl> <dbl>
##  1    18 <tibble> ATime_FE        <dbl>      <dbl>    <dbl [95,760]> 0.401 0.651
##  2    19 <tibble> ATime_FE        <dbl>      <dbl>    <dbl>          0.558 0.698
##  3    18 <tibble> BSpace_FE       <dbl>      <dbl>    <dbl [95,760]> 0.398 0.581
##  4    19 <tibble> BSpace_FE       <dbl>      <dbl>    <dbl>          0.475 0.649
##  5    18 <tibble> CTime_Space_FE  <dbl>      <dbl>    <dbl [95,760]> 0.391 0.576
##  6    19 <tibble> CTime_Space_FE  <dbl>      <dbl>    <dbl>          0.479 0.636
##  7    18 <tibble> DTime_Space_FE… <dbl>      <dbl>    <dbl [95,760]> 0.340 0.546
##  8    19 <tibble> DTime_Space_FE… <dbl>      <dbl>    <dbl>          0.400 0.606
##  9    18 <tibble> ETime_Space_FE… <dbl>      <dbl>    <dbl [95,760]> 0.340 0.546
## 10    19 <tibble> ETime_Space_FE… <dbl>      <dbl>    <dbl>          0.400 0.606

Examine Error Metrics for Accuracy

By plotting MAE by model by week, the spatial fixed effects and temporal+spatial fixed effects do a similar job in predicting.However, the model starts to become more accurate when lag effects are taken into consideration, with a MAE under 0.4, which is exactly same as the model which further takes amenities features into consideration.

From the predicting performances plot, the model is getting more accurate and generalizable after applying time, space, lag features. The lag model does a good performance in predicting the bike share in first two weeks of May, although it loses some observations at the beginning of the month and also underestimates a little bit.

week_predictions %>% 
    mutate(interval60 = map(data, pull, interval60),
           start_station_id = map(data, pull, start_station_id)) %>%
    dplyr::select(interval60, start_station_id, Observed, Prediction, Regression) %>%
    unnest() %>%
    gather(Variable, Value, -Regression, -interval60, -start_station_id) %>%
    group_by(Regression, Variable, interval60) %>%
    summarize(Value = sum(Value)) %>%
    ggplot(aes(interval60, Value, colour=Variable)) + 
      geom_line(size = 1.1) + 
      facet_wrap(~Regression, ncol=1) +
      scale_color_manual(values = palette2) +
      labs(title = "Predicted/Observed bike share time series", subtitle = "New York City; A test set of 2 weeks",  x = "Hour", y= "Station Trips") +
      plotTheme()

## `summarise()` has grouped output by 'Regression', 'Variable'. You can override
## using the `.groups` argument.

Based on the Space_Time_Lag(_amenity) model(Reg4/Reg5) which seems to have the best goodness of fit generally, we can observe some spatial patterns from the mean absolute errors by station. Specifically, the stations with higher MAE are aggregated in Midtown & Lower Manhattan.

## `summarise()` has grouped output by 'start_station_id', 'start_lng'. You can
## override using the `.groups` argument.

Space-Time Error Evaluation

By comparing observed versus predicted ridership for various times of day across weekdays and weekends, the presence of space-time errors becomes apparent. The majority of the data points falling below the identity line suggests that the model generally underestimates bike-share usage. This trend is particularly pronounced during the weekday morning rush hours. Furthermore, the slope of the best fit linek varies across different time slots for both weekdays and weekends, indicating fluctuations in the model’s predictive accuracy dependent on the specific time period being analyzed. Therefore, more factors should be taken into consideration.

week_predictions %>% 
    mutate(interval60 = map(data, pull, interval60),
           start_station_id = map(data, pull, start_station_id), 
           start_lat = map(data, pull, start_lat), 
           start_lng = map(data, pull, start_lng),
           dotw = map(data, pull, dotw)) %>%
    select(interval60, start_station_id, start_lng, 
           start_lat, Observed, Prediction, Regression,
           dotw) %>%
    unnest() %>%
  filter(Regression == "DTime_Space_FE_timeLags")%>%
  mutate(weekend = ifelse(dotw %in% c("周日", "周六"), "Weekend", "Weekday"),
         time_of_day = case_when(hour(interval60) < 7 | hour(interval60) > 19 ~ "Overnight",
                                 hour(interval60) >= 7 & hour(interval60) < 10 ~ "AM Rush",
                                 hour(interval60) >= 10 & hour(interval60) < 15 ~ "Mid-Day",
                                 hour(interval60) >= 15 & hour(interval60) <= 19 ~ "PM Rush"))%>%
  ggplot()+
  geom_point(aes(x= Observed, y = Prediction))+
    geom_smooth(aes(x= Observed, y= Prediction), method = "lm", se = FALSE, color = "#FFA7C4" )+
    geom_abline(slope = 1, intercept = 0,color="#622A8C")+
  facet_grid(time_of_day~weekend)+
  labs(title="Observed vs Predicted by the day of the week and hour",
       x="Observed trips", 
       y="Predicted trips")+
  plotTheme()

## `geom_smooth()` using formula = 'y ~ x'

These maps further show that high errors are concentrated in midtown & lower Manhattan - where bike-share riderships are also high.

## `summarise()` has grouped output by 'start_station_id', 'weekend',
## 'time_of_day', 'start_lng'. You can override using the `.groups` argument.

Focusing on the morning commute, the model doesn’t perform well for specific socio-economic groups characterized by higher income, lower public transit usage, and higher white population %.

week_predictions %>% 
    mutate(interval60 = map(data, pull, interval60),
           start_station_id = map(data, pull, start_station_id), 
           start_lat = map(data, pull, start_lat), 
           start_lng = map(data, pull, start_lng),
           dotw = map(data, pull, dotw),
           Percent_Taking_Public_Trans = map(data, pull, Percent_Taking_Public_Trans),
           Med_Inc = map(data, pull, Med_Inc),
           Percent_White = map(data, pull, Percent_White)) %>%
    select(interval60, start_station_id, start_lng, 
           start_lat, Observed, Prediction, Regression,
           dotw, Percent_Taking_Public_Trans, Med_Inc, Percent_White) %>%
    unnest() %>%
  filter(Regression == "DTime_Space_FE_timeLags")%>%
  mutate(weekend = ifelse(dotw %in% c("周日", "周六"), "Weekend", "Weekday"),
         time_of_day = case_when(hour(interval60) < 7 | hour(interval60) > 19 ~ "Overnight",
                                 hour(interval60) >= 7 & hour(interval60) < 10 ~ "AM Rush",
                                 hour(interval60) >= 10 & hour(interval60) < 15 ~ "Mid-Day",
                                 hour(interval60) >= 15 & hour(interval60) <= 19 ~ "PM Rush")) %>%
  filter(time_of_day == "AM Rush") %>%
  group_by(start_station_id, Percent_Taking_Public_Trans, Med_Inc, Percent_White) %>%
  summarize(MAE = mean(abs(Observed-Prediction), na.rm = TRUE))%>%
  gather(-start_station_id, -MAE, key = "variable", value = "value")%>%
  ggplot(.)+
  geom_point(aes(x = value, y = MAE), alpha = 0.4)+
  geom_smooth(aes(x = value, y = MAE,color="#622A8C"), method = "lm", se= FALSE)+
  facet_wrap(~variable, scales = "free")+
  labs(title="Errors as a function of socio-economic variables",
       y="Mean Absolute Error (Trips)")+
  plotTheme()

Cross-validation

To further validate the generalizability of the model, I cross validate the model by time and space. To manage the large number of data points, the model deemed most effective was subjected to a 50-fold cross-validation. This means that the data was divided into 50 parts, with the model being trained on 49 parts and tested on the 1 remaining part, this process being repeated 50 times with different parts each time.

As we can see, the mean absolute error is slightly over 0.37. In terms of the CV goodness of fit metrics,the metrics are clusterd close to the man, suggesting that the model performs consistently and shows good generalizabilty.

## Linear Regression 
## 
## 287712 samples
##     10 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (50 fold) 
## Summary of sample sizes: 281958, 281958, 281957, 281957, 281957, 281958, ... 
## Resampling results:
## 
##   RMSE       Rsquared   MAE      
##   0.6887728  0.3352304  0.3676375
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

Mean Absolute Error	Standard Deviation of MAE
0.3676375	0.0079647

dplyr::select(cvFit$resample, -Resample) %>%
  gather(metric, value) %>%
  left_join(gather(cvFit$results[2:4], metric, mean)) %>%
  ggplot(aes(value)) + 
    geom_histogram(bins=35, fill = "#FFA7C4") +
    facet_wrap(~metric) +
    geom_vline(aes(xintercept = mean), colour = "#B14A90", linetype = 3, size = 1.5) +
    scale_x_continuous(limits = c(0, 1)) +
    labs(x="Goodness of Fit", y="Count", title="CV Goodness of Fit Metrics",
         subtitle = "50 folds, Across-fold mean reprented as dotted lines") +
    plotTheme()

## Joining with `by = join_by(metric)`

## Warning: Removed 6 rows containing missing values (`geom_bar()`).

Conclusion

In this project, the focus was on developing a predictive model for bike rebalancing within New York City’s Citi Bike system using a combination of serial data, spatial characteristics, and time-lag variables. The incorporation of amenity features did not significantly enhance the model’s performance, indicating that the selected variables were already sufficient in capturing the patterns necessary for accurate bike-share ridership predictions. Despite encountering space-time errors, with a tendency to underpredict overall ridership, the model demonstrates effectiveness and provides a solid basis for planning and optimizing Citi Bike’s rebalancing operations.

To further improve the model’s accuarcy and generalizabilty, more deeper analysis into the temporal patterns and spatial distributions of ridership can be taken into consideration.