Rebalancing is a common and important problem facing bike share systems around the globe. Below, we consider the value of a range of linear models incorporating spatiotemporal data to predict future rideshare demand for the Indego bikeshare network in Philadelphia. We find that spatiotemporal features have high predictive power for this use case and enable us to construct a model worth deploying. Although the predictions are still hampered by non-spatially random error, indicating the possibility of further improving the model, we recommend deploying the model as a useful tool in ride share balancing.
Bike sharing systems have emerged as a popular, eco-friendly transportation alternative in cities worldwide. These systems provide a network of bicycles and docking stations, allowing users to pick up a bike from one location and drop it off at another, facilitating short, point-to-point trips. However, one significant challenge in these systems is maintaining a balanced distribution of bikes across the network. This is where bike re-balancing becomes crucial.
Due to daily commuting patterns and other socio-economic factors, certain areas in a city may experience a surplus of bikes while others face a shortage. For instance, residential areas might see a high number of bikes in the morning as people commute to work, leaving these stations depleted by evening. Conversely, business districts may have an excess of bikes during the day, which need to be redistributed for the evening commute home. Without effective re-balancing, the utility and efficiency of the bike share system are significantly reduced, leading to customer dissatisfaction and a decline in usage.
Below, we propose and evaluate a rebalancing strategy for Indego, the bike share system in Philadelphia, based on a spatiotemporal predictive model. Launched in 2015, Indego has 1,400 bikes distributed across 130 stations, mostly concentrated in Center City, West Philadelphia, and South Philadelphia.1 Using Indego data from the 4th quarter of 2022, as well as supplemental weather data, we propose a predictive model optimized for overnight rebalancing with trucks, which aims to minimize expected customer dissatisfaction the next day.2 Based on our predictions, an integer program could be used to optimize redistribution routes.3
For our analysis, we rely on three datasets: 1) Indego station and ridership data from Q4 of 2022, 2) weather data collected at Philadelphia International Airport during those same dates and accessed via the riem package in R, and 3) American Community Survey data from 2021.
3.1.2 Exploratory Analysis
Examining our bike share data, we notice clear temporal trends: peaks during certain hours of the day and a longer-time decline in trips across the quarter.
Show the code
ggplot(dat_census %>%group_by(interval60) %>%tally())+geom_line(aes(x = interval60, y = n))+labs(title="Bike share trips per Hour",subtitle ="Philadelphia, Oct. - Dec., 2022",x="Date", y="Number of trips")+ plotTheme
Breaking these trends down, we notice peaks at certain times of day and certain days of the week. As may be expected, morning and afternoon rush hours are associated with the highest ridership volumes during the week. Overall ridership is lower on the weekend, and peaks during the midday instead of morning or afternoon.
Show the code
dat_census %>%mutate(hour =hour(interval60)) %>%ggplot(aes(x = hour, color =as.factor(dotw)))+geom_freqpoly(binwidth =1)+labs(title="Bike share trips in Philadelphia, by day of the week, Oct. - Dec., 2022",x="Hour", y="Trip Counts",color ="Day of the Week")+ plotTheme
Show the code
ggplot(dat_census %>%mutate(hour =hour(interval60),weekend =ifelse(dotw %in% weekend, "Weekend", "Weekday")))+geom_freqpoly(aes(x = hour, color = weekend), binwidth =1)+labs(title="Bike share trips in Philadelphia, Oct. - Dec., 2022, \nweekend vs. weekday",x="Hour", y="Trip Counts",color ="Weekend?")+ plotTheme
Although the average number of trips per hour per station remains around 1.3 no matter the time of day, a higher proportion of stations experience ridership demands exceeding this during the AM and PM rush as compared to overnight or during the middle of the day.
Show the code
dat_census %>%mutate(time_of_day =case_when(hour(interval60) <7|hour(interval60) >18~"Overnight",hour(interval60) >=7&hour(interval60) <10~"AM Rush",hour(interval60) >=10&hour(interval60) <15~"Mid-Day",hour(interval60) >=15&hour(interval60) <=18~"PM Rush"))%>%group_by(interval60, start_station, time_of_day) %>%tally()%>%group_by(start_station, time_of_day)%>%summarize(mean_trips =mean(n))%>%ggplot()+geom_density(aes(mean_trips, color = time_of_day, fill = time_of_day), alpha =0.3)+labs(title="Mean Number of Hourly Trips Per Station. Philadelphia, Oct. - Dec., 2022",x="Number of trips", y="Density",color ="Time of Day",fill ="Time of Day")+#scale_x_log10()facet_wrap(~time_of_day)+ plotTheme
Additionally, these trends are born out spatially, as we observe that the highest ridership per station clusters in the downtown and University City areas during the afternoon rush hour.
To build our model, we incorporate various predictors, including weather and spatial and temporal lag. Our main contribution in terms of feature engineering is to add temporal lags. This is important, as we find that they are strongly correlated with our dependent variable, Trip_Count.
Correlations between Engineered Features and Dependent Variable
Variable
correlation
lagHour
0.90
lag2Hours
0.75
lag3Hours
0.58
lag4Hours
0.42
lag12Hours
-0.31
lag1day
0.83
Finally, we note clear spatial and temporal autocorrelation in our data: this animated map, for example, indicates that data cluster at certain times of the day in certain places.
For our analysis, we use a range of OLS regressions. Given the spatial autocorrelation noted above, OLS may not be the ideal approach; a spatial lag or geographically weighted regression, for example, might improve performance and generalizability. Overall, however, for the purposes of this task, OLS is sufficient.
Using 100-fold cross validation, we find that the MAE for our best model is 0.52`–not bad. This suggests that, for a given station at a given hour, our model will, on average, correctly predict demand to within one bike.
Of the six models we consider, we find that those incorporating both time and space features, as well as the lag of various time features, perform the best. Interesting, temporal features seem more important than spatial features; incorporating time lags into the model yielded significant improvement, while incorporating spatial lag resulted in only a marginal further improvement. The best models are those that incorporate time lag; the inclusion of holiday lags, spatial lag, and socioeconomic variables yield little to no improvement. It’s likely that these features are in some ways already implicit in the pure spatial and temporal features of the dataset; given the highly spatially-autocorrelated distribution of income, for instance, this variable is already baked into the location data.
Show the code
week_predictions %>% dplyr::select(week, Regression, MAE) %>%gather(Variable, MAE, -Regression, -week) %>%ggplot(aes(week, MAE)) +geom_bar(aes(fill = Regression), position ="dodge", stat="identity") +scale_fill_manual(values = palette7) +labs(title ="Mean Absolute Errors by model specification and week") + plotTheme +theme(legend.position ="bottom")
One potential challenge posed by our models is that they consistently underpredict rather than overpredict. When trying to meet rideshare demand, this poses a problem; if Indego consistently delivers too many bikes, they will never miss out on ridership, but if they consistently underprepare, they will lose potential riders and even possibly cost themselves long-term subscribers.
Show the code
week_predictions %>%mutate(interval60 =map(data, pull, interval60),start_station =map(data, pull, start_station)) %>% dplyr::select(interval60, start_station, Observed, Prediction, Regression) %>%unnest() %>%gather(Variable, Value, -Regression, -interval60, -start_station) %>%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) +labs(title ="Predicted/Observed bike share time series", subtitle ="Philadelphia; A test set of 3 Months", x ="Hour", y="Station Trips") + plotTheme
We also observe that the errors in our model are not spatially random, indicating that we have not fully accounted for spatial fixed effects in our model. In particular, they cluster in areas of high demand. Consistent underprediction in areas of high demand undermines the viability of our model, as these are precisely the areas in which we would hope to have the most accurate predictions.
If we plot observed vs. predicted for different times of day during the week and weekend, some patterns begin to emerge.
Underprediction: Across almost all plots, the predicted values tend to be lower than the observed values, especially as the number of observed trips increases. This is indicated by the data points lying above the 45-degree line (which represents perfect prediction). In other words, our model struggles more with the higher-demand stations.
Time-of-Day Effects: There might be different patterns of underprediction based on the time of day. For example, during the AM Rush and PM Rush, the difference between observed and predicted seems larger than during midday or overnight. This could suggest that the model does not capture the peak travel times as accurately.
Weekday vs. Weekend: There could be a difference in the prediction quality between weekdays and weekends. The weekend plots might show more scatter away from the line than the weekdays, indicating a higher variance in the accuracy of predictions during weekends.
Model Limitations: The scatter and the deviation from the line could indicate that the model is missing some factors that influence the number of trips, especially during specific times like rush hours or weekends when travel patterns could be significantly different from the usual.
Outliers: There are some data points that are far away from the majority of the data, especially in the higher range of observed trips. These outliers suggest that there are extreme cases where the model fails to predict the correct number of trips by a large margin.
In order to improve our model, we could attempt to account for some of these issues with further feature engineering. For example, we might add binary variables to mark peak hours, attempt to factor in special events, or even consider switching to different models. One simple upgrade here could be using a penalized model such as a ridge, lasso, or elastic net to account for what is effectively overfitting. Other approaches, such as a generalized additive model (GAM) or a tree-based model (e.g., a Random Forest) might also help.
Lastly, we find that our model does not generalize perfectly; errors are higher in wealthier and whiter neighborhoods, and lower in neighborhoods that are more reliant on public transportation (in reality, these are likely all the same neighborhoods, as these variables are highly colinear). This introduces an interesting optimization problem: do we prefer to optimize for an equitable distribution of error, whereby we seek to avoid burdening low-income, minority neighborhoods with disproportionate rates of error, or do we prefer to minimize error in high-demand areas (which are predominately white and wealthy)?
All that said, the predictive power of our model is good. With an MAE around 0.5, we can expect to be fairly close (within roughly 1 bike) to the accurate demand for bikes at a given station at any given hour. The major challenge to address is underprediction at high-demand stations; this is both a question of model improvements and rebalancing priorities. Not only do we want to maximize model accuracy, but we also need to consider at which stations we most care about accuracy. This shades into a cost-benefit analysis.
Further improvement of the model depends on two factors: improved feature engineering and more computing power. Accounting for features like rush hour with binary variables, for example, might help. The best thing, however, would likely be to basically throw as much data into the model as possible. Simply, training on more historical data would allow for a more sophisticated model that could, for example, distinguish between different holidays that might prompt different levels of commuting (Veterans’ Day versus Christmas Eve, for instance), or simply learn from more historic data. Likewise, more sophisticated models (random forest, Lasso regression, etc.) might increase model accuracy. One key component would be switching from an OLS regression to one that explicitly incorporated spatial dependencies, such as spatial lag regression or geographically-weighted regression, in addition to incorporating temporal features. (This, however, is beyond the scope of this assignment and likely would require more computing power than we have access to). In truth, though, the best approach is probably the best machine learning approach in general: throw as much data as possible into the model. Training on the roughly 7 years of available Indego data, rather than the current three months’ worth, would probably do more than anything else to increase model accuracy and generalizability–but would also be prohibitively computational intensive, at least in the context of this assignment.
6 Conclusion
We find that our model is effective in predicting rideshare demand and would be useful in informing rebalancing efforts for Indego. Overall, the model’s accuracy is high enough to be useful–certainly in comparison to the naive approach of responding to demand as it emerges. That said, the model still suffers from non-random spatial distribution of errors. Specifically, it does not generalize well to high-demand areas, which is a substantial issue given that these are exactly the areas where the impact of underprediction is biggest. Thus, while the model is still a solid predictive tool worth implementing, there are important opportunities to improve it by better accounting for spatial fixed effects.
