Analysis of San Francisco Bay Area Bike Sharing Program

1 Introduction

Introduction of SF Bay Area Bike Share: http://www.bayareabikeshare.com/

Data used: Year 1 Data, from Mar 2014 - Aug 2014, trips data ONLY.

One of the most desirable wishes of the riders who uses this bike-sharing program is that they are able to get a bike at the nearest station from their current location and then able to dock at their destination bike station upon arriving. However, this is unlikely the case for some because of the irregularities in people movement across the days, timings of the day (think peak-hours traffic) and the starting and end location (most people would be heading towards central area).

This is a headache for the managers in this bike-program as they seek to keep their riders happy. What are some of the actions that they could take to ensure this problem is mitigated? One of the current practices is to assign employees to move the bikes from areas with least demand at certain timings to other areas with much higher demand. However, how would a manager know which station has highest demand, least demand, what time to move? To add on, the days of the week and who are using it also would matter. In fact, there are also limited employees to shift the bikes around on a certain day.

Hence, this submission seeks to reveal some insights related to this problem such as - what is the distribution of the number of rides across the dates, days and time. What are the most common to and fro paths? Eventually, from these exploratory data analysis, we will be able to understand the data better and then create an appropriate prediction model and tool to help manager understand the situtation much better! This tool will help to inform the manager in advance, which day would have the highest increment and fall in number of bikes. Hence, the manager will then be able to station his employees across those stations and transport the bikes to balance the situation.

After conducting a Google search on the stations, I realized that the data consists of rides and stations across different cities. Hence, for this analysis, I would focus on the trips that start or end stations within San Francisco only (i.e. includes trips that start from SF stations and end at non-SF stations, trips that start from non-SF stations and end at SF stations)

Several hypothesis that could be uncovered:

• Are there more riders on the weekdays or weekends?

• Are there more customers or subscribers using the service?

• How does the distribution of trips look like (across time, days of week, dates and certain stations)?

• Which are the stations that tend to be the most ‘active’?

2 Exploratory Data Analysis

To start off, let’s load the necessary packages, the CSV data file provided and the stations within San Francisco.

** Please place the tripdata file in the working directory for the code to run **

library(ggplot2)
library(lubridate)
library(plyr)
library(dplyr)
library(forecast)
tripdata <- read.csv("./tripdata.csv")
sf <- c("South Van Ness at Market", "Market at 10th", "San Francisco City Hall", 
        "Golden Gate at Polk", 
        "Civic Center BART (7th at Market)", "Powell Street BART", 
        "Powell at Post (Union Square)", "5th at Howard", "Market at 4th", "Post at Kearny", 
        "Commercial at Montgomery", "Washington at Kearny", "Grant Avenue at Columbus Avenue", 
        "Embarcadero at Sansome", "Embarcadero at Vallejo", "Broadway St at Battery St", 
        "Davis at Jackson", "Clay at Battery", "Harry Bridges Plaza (Ferry Building)", 
        "Steuart at Market", "Beale at Market", "Mechanics Plaza (Market at Battery)", 
        "Embarcadero at Folsom", "Temporary Transbay Terminal (Howard at Beale)", 
        "Market at Sansome", "Spear at Folsom", "Howard at 2nd", "Embarcadero at Bryant", 
        "Yerba Buena Center of the Arts (3rd @ Howard)", "2nd at Folsom", "2nd at South Park", 
        "2nd at Townsend", "San Francisco Caltrain (Townsend at 4th)", "Townsend at 7th",
        "San Francisco Caltrain 2 (330 Townsend)")

2.1 Removal of trips within non-SF cities

### Create 2 columns that tag the start and end stations on whether they are SF stations ###
tripdata$go <- ifelse(tripdata$Start.Station %in% sf, "SF", "NotSF")
tripdata$back <- ifelse(tripdata$End.Station %in% sf, "SF", "NotSF")
tripdata$goback <- paste(tripdata$go, tripdata$back, sep = " to ")

Having done so, let’s have a look at the distributions of the trips of inter-cities and intra-cities stations within this data.

### Group up by the *goback* variable, before tallying the results up and ordering
tripdata %>% group_by(goback) %>% tally() %>% arrange(desc(n))

Interesting! There are 19031 trips made in non-SF cities within the data. We have also found people who have in fact travelled across cities (14 such cases), must have took them some time which I wonder how long, but that will be for another time :) In addition, this works out well in our favour because 14 is a negligible number, thus we can remove them without affecting the numbers much!

Now, let’s retain only those trips that are within “SF to SF”.

tripdata <- filter(tripdata, goback == 'SF to SF')
### Also to remove the tagged columns that were added as they are no longer required ###
tripdata <- tripdata[,-c(12,13,14)]

Right now, our dataset contains of trips that are only within SF city, great! Now we can move on to perform some simple exploratory analysis to see what’s been going on.

2.2 Further cleaning of dataset

### Convert the Start.Date and End.Date columns into Date-time class ###
StartDate <- strptime(tripdata$Start.Date, "%m/%d/%Y %H:%M")
tripdata$Start.Date <- as.POSIXct(StartDate) 
EndDate <- strptime(tripdata$End.Date, "%m/%d/%Y %H:%M")
tripdata$End.Date <- as.POSIXct(EndDate) 
### Remove the unnecessary stored values ###
remove(StartDate,EndDate)
### Creates a copy of the original dataset ###
original <- tripdata
### Creates several useful columns (breakdown start/end date by Month, Day, Day of Week for analysis ###
tripdata <- mutate(tripdata, sdate = date(Start.Date), smonth = month(Start.Date,label = TRUE), sday = day(Start.Date), swday = wday(Start.Date,label = TRUE), shr = hour(Start.Date), edate = date(Start.Date), emonth = month(End.Date,label = TRUE), eday = day(End.Date), ewday = wday(End.Date,label = TRUE), ehr = hour(End.Date))
### To turn the Start.Date and End.Date into character format to avoid conflicting issues in tallying ###
tripdata$Start.Date <- as.character(tripdata$Start.Date)
tripdata$End.Date <- as.character(tripdata$End.Date)

2.3 Distribution of rides by days of week, time, dates and subscriber type

With the cleaning up done, let’s now tally up the number of rides across the weekdays and see what’s the distribution like to answer our hypothesis!

### Tally up the number of trips for certain weekday ###
cwd <- ddply(tripdata, .(swday), tally)
cwd$wkday <- ifelse(cwd$swday %in% c("Mon","Tues","Wed","Thurs","Fri"), "Weekday","Weekend")
ggplot(cwd, aes(x = swday, y = n, fill = wkday)) + geom_bar(stat='identity') + labs(title = "Ridership over days of week", x = "Days of Week", y = "Count")

The result is as expected, the number of rides occurred mainly during weekdays and fell to a lower number during weekends. Hence, we would expect the availability of bikes problem to occur more during weekdays.

This brings me to the point, could the number of non-subscribers drive up the usage on weekdays? As much as we want to cater to both groups of customers, however, if there is a large number of non-subscribers using the bikes, this would cause much displeasure to our subscribers. Hence, we may want to increase the prices for customers usage such that we can allow more subscribers to use instead.

Let’s find out how does the proportion of customers v.s. subscribers look like in the distribution.

cwdrider <- ddply(tripdata, .(swday,Subscriber.Type), tally)
ggplot(cwdrider, aes(x = swday, y = n, fill = Subscriber.Type)) + geom_bar(stat='identity') + labs(title = "Ridership over days of week by Subscriber Type", x = "Days of Week", y = "Count")

Looks like the customers didn’t create high demand during the weekdays… However, it would be too quick to simply dismiss them as not a cause. Could they utilize the service mostly during the peak hours? If so, this could cause a problem as well!

Now let’s move on to examining the distribution of ridership over the hour of the day. We shall construct similar code but in the hours context.

### Just some cleaning up to do ###
remove(cwd,cwdrider)
###
hrrider <- ddply(tripdata, .(shr,Subscriber.Type), tally)
ggplot(hrrider, aes(x = shr, y = n, fill = Subscriber.Type)) + geom_bar(stat='identity') + labs(title = "Ridership over time of day by Subscriber Type", x = "Time of day (hr)", y = "Count")

Whew! Looks like customers are not causing too much of a problem then it seems… In addition, we are also able to see the distribution of trips over the time of the day. For curiosity’s sake, let’s combine all these factors together and see how does ridership vary over the days of week over hour of the day, broken down by subscriber type!

whrider <- ddply(tripdata, .(shr,swday,Subscriber.Type), tally)
ggplot(whrider, aes(x = shr, y = n, fill = Subscriber.Type)) + facet_grid(. ~ swday) + geom_bar(stat='identity') + labs(title = "Ridership over WDay over Time by Subscriber Type", x = "Time", y = "Count")

3 quick points that we can infer from the graphs:

• The distribution of trips over time during weekdays is approximately the same across all 5 days, having peaks at 7 - 9 a.m. and 4 - 6 p.m.

• The distribution of trips over the weekend is, different from weekdays but similar across both days!

• The number of trips taken by customers are roughly more than the number of trips taken by subscribers on weekends, which is a good sign :)

It is also intuitive to have a look at how does the ridership change over time within the given period of March to August, and the effects of weekends and holidays have on them!

Extracted several holiday dates celebrated in the U.S. in 2014

• 11 May 2014 (Sun) - Mother’s Day

• 26 May 2014 (Mon) - Memorial Day

• 15 Jun 2014 (Sun) - Father’s Day

• 04 Jul 2014 (Fri) - Independence Day

### Some cleaning up as usual ###
remove(hrrider,whrider)
### Convert the initial start date back to POSICXlt class, and renamed it ###
tripdata$Start.DatewTime <- strptime(tripdata$Start.Date, "%Y-%m-%d %H:%M:%S")
tripdata$End.DatewTime <- strptime(tripdata$End.Date, "%Y-%m-%d %H:%M:%S")
### Extracted the date out and turned into character class ###
tripdata$Start.Date <- as.character(date(tripdata$Start.DatewTime))
tripdata$End.Date <- as.character(date(tripdata$End.DatewTime))
### Convert the DateWTime back to character so that we can tally w/o issues ###
tripdata$Start.DatewTime <- as.character(tripdata$Start.DatewTime)
tripdata$End.DatewTime <- as.character(tripdata$End.DatewTime)
### To add columns for holidays and weekends ###
tripdata <- mutate(tripdata, wkend = (swday %in% c("Sat","Sun")),
                   hol = (Start.Date %in% c("2014-05-11","2014-05-26","2014-06-15",
                                                  "2014-07-04")))
### Tally up the number of trips by start date ###
tripdate <- ddply(tripdata, .(Start.Date, wkend, hol), tally)
### Convert into date format
tripdate$Start.Date <- date(tripdate$Start.Date)
### Plot the graph with red lines for weekends, blue lines for holidays
ggplot(tripdate, aes(Start.Date, n)) + geom_line() + geom_area(aes(y = wkend*max(n)), fill ="red", alpha = 0.30) + geom_area(aes(y = hol*max(n)), fill = "blue", alpha = 0.30) +
    labs(title = "No. of trips from Mar - Aug", x = "Date", y = "Count") + theme_minimal()

When plotted across the dates from 1st Mar 14 - 31st Aug 14, there is a fall in ridership on every weekend, this finding is consistent with our previous findings. Yet, there also seems to be a slight drop for the days with holidays, although it seems minor/negligible.

Also, we note that there is a general increment trend as well which is a good sign! This could very well be due to the transition of season, from spring to summer, thus the increase in number of trips?

In addition, it seems like the distribution of the trips follow a general trend with seasonality effect, which we could actually predict using the Holts-Winter model..

After knowing how does the ridership number change over several variables, now we wonder if the ridership displays a homogenous distribution across the paths taken? This is most unlikely to be the case because we should expect several stations that are nearer to the centre to be more crowded/congested during peakhours!

We shall now investigate by looking at the top 5 starting and ending locations

### General cleaning ###
remove(tripdate)
### 
startpaths <- ddply(tripdata, .(Start.Station), tally) %>% arrange(desc(n))
head(startpaths)

endpaths <- ddply(tripdata, .(End.Station), tally) %>% arrange(desc(n))
head(endpaths)

The top 6 locations that appear in the Start Station also appear in the top 6 ending locations. This is not surprising because it is very likely that these 6 stations are places where people usually head to work and return from. This is in fact a good sign, because it would be worrying if these 2 results do not match. This would imply that riders are heading to a location, which they would very likely not return from, or vice versa. This would result in a huge imblanace in the number of bikes at each location.

However, the hypothesis still remains! If everyone is travelling to, for example, San Francisco Caltrain (Townsend at 4th), then we will be expecting that station to be full quickly!

Let’s see if this is the case for the top 6 stations and bottom 6 stations across time.

startdata <- tripdata[tripdata$Start.Station %in% startpaths$Start.Station[1:6],] %>%
    ddply(.(Start.Station,shr, wkend), tally) 
startdata$wkend <- ifelse(startdata$wkend == TRUE, 'Weekend', 'Weekday')
ggplot(startdata, aes(x = shr, y = n, colour = wkend)) + facet_wrap( ~ Start.Station, ncol = 2) + geom_line(aes(group = wkend)) + geom_point(aes(group = wkend)) + labs(title = "Distribution of trips starting from each station across time by weekday/weekend", x = "Time (hr)", y = "Count")

We notice that each of the 6 stations have very different distributions across time and different magnitude as well.

Let’s examine the case for the top 6 ending stations with highest number of trips.

enddata <- tripdata[tripdata$End.Station %in% endpaths$End.Station[1:6],] %>%
    ddply(.(End.Station,shr, wkend), tally) 
enddata$wkend <- ifelse(enddata$wkend == TRUE, 'Weekend', 'Weekday')
ggplot(enddata, aes(x = shr, y = n, colour = wkend)) + facet_wrap( ~ End.Station, ncol = 2) + geom_line(aes(group = wkend)) + geom_point(aes(group = wkend)) + labs(title = "Distribution of trips starting from each station across time by weekday/weekend", x = "Time (hr)", y = "Count")

Based on the graphs displayed, this does confirm our hypothesis that these tend to be places when people would usually go/leave in the morning and then go/leave at night! This is evident from the reflection of the first graph about the y axis, which resembles the distribution of the second graph.

From here, we could in fact see that San Francisco Caltrain (Townsend at 4th), San Francisco Caltrain 2 (330 Townsend) and Harry Bridges Plaza are places where people either go home to or to take on another form of transport (most likely the commuter train) in the evening to some where else (high number of departures in the evening)! On the other end, the other 3 stations are places where many people are most likely to work at, hence the high number of arrivals in the morning and high number of departures in the evening.

Hence, we can then see how several factors such as time of the day, day of the week, date, subscriber type can then affect the general demand and supply across the stations.

3 Creation of the prediction tool

Next step here, we will attempt to model the predicted net change of number of bikes at the San Francisco Caltrain (Townsend at 4th) station using the data provided. This will be done by several steps.

1. Extracts out the trips that either start or end at the station.

2. Separate this dataset into 2 - one for trips starting at that station, one for trips ending at that station.

3. Extracts only the date information out of the Start.DatewTime and End.DatewTime columns.

4. Aggregate the number of trips by day.

5. Create a single column of all the dates running from 1st Mar 14 to 31st Aug 14.

6. Merge these datasets together by date.

7. Trips that start at that station signifies number of bikes leaving the station, trips that end at that station signifies the number of bikes leaving the station. Thus, the number of bikes ending at station deducted by trips starting at station, would give us the net change (i.e. a number of positive 8 at a particular date, would represent that there is a net increase of 8 bikes at that date)

8. Plot the net change of bikes at the trip over date.

9. Fit into the Holt-Winter’s additive model (found to be the best)

• The tool will show the top 3 highest and top 3 lowest number of incoming bikes for each day. It isn’t feasible to show by day by hour, because it is afterall an estimate. If we were to follow the tool strictly, then we could have been easily caught offguard if the prediction isn’t right. *

### General cleaning up as usual ###
remove(startdata,enddata,startpaths,endpaths)
### Extracts out the data that start or end at San Francisco Caltrain (Townsend at 4th) ###
sfct <- tripdata[tripdata$Start.Station %in% "San Francisco Caltrain (Townsend at 4th)" | tripdata$End.Station %in% "San Francisco Caltrain (Townsend at 4th)",]
### Splits the data set into 2, one that starts at that station, another one that ends at the station ###
startdt <- select(sfct, Start.Station, Start.DatewTime) %>% filter(Start.Station == "San Francisco Caltrain (Townsend at 4th)")
enddt <- select(sfct, End.Station, End.DatewTime) %>% filter(End.Station == "San Francisco Caltrain (Townsend at 4th)")
### Extracts only the date and hour out ###
startdt$Start.DatewTime <- as.POSIXct(format(as.POSIXct(startdt$Start.DatewTime), "%Y-%m-%d"))
enddt$End.DatewTime <- as.POSIXct(format(as.POSIXct(enddt$End.DatewTime), "%Y-%m-%d"))
### Aggregate the number of trips by date and hour ###
startdt <- aggregate(Start.Station ~ Start.DatewTime, data = startdt, FUN = length)
enddt <- aggregate(End.Station ~ End.DatewTime, data = enddt, FUN = length)
### Create a column of sequential date and time from 1st Mar 14 to 31st Aug 14 ###
datetime <- seq(from = as.POSIXct("2014-03-01"), to = as.POSIXct("2014-08-31"), by = "day")
datetime <- data.frame(datetime)
### Merge all the subsequent columns together, replace NA values with 0 ###
combdt <- merge(datetime, startdt, by.x = "datetime", by.y = "Start.DatewTime", all.x = TRUE)
combdt$Start.Station[is.na(combdt$Start.Station)] <- 0
combdt <- merge(combdt, enddt, by.x = "datetime", by.y = "End.DatewTime", all.x = TRUE)
combdt$End.Station[is.na(combdt$End.Station)] <- 0
### The start.station column signifies number of bikes leaving the station, the end.station column signifies the number of bikes entering the station. To get the net change, we form a new column of Net = End.Station - Start.Station ###
combdt <- mutate(combdt, Net = End.Station - Start.Station)
### How the data looks like for head and tail ###
head(combdt)

tail(combdt)

### Plot the data to see the distribution of net change in number of bikes across date and hour ###
ggplot(combdt, aes(x = datetime, y = Net)) + geom_line() + labs(title = "Net change of bikes over date", x = "Net change of bikes", y = "Date")

What we have noticed here is that the number of trips definitely follows a pattern, due to the change in usage from weekdays to weekend. Hence, we tend to notice 4 peak periods in each month (the 4 set of weekdays in each month). In addition, it also seems that there is a larger variation of net change as time progresses.

With these findings in mind, we can now definitely attempt to plot Holt-Winters model to predict what could the net change be in the upcoming month!

combdt <- select(combdt, datetime, Net)
combdts <- ts(combdt$Net, frequency = 7, start = c(1,3))
### Fit an additive Holt Winter's Additive model ###
fitHW <- HoltWinters(combdts, seasonal = "additive")
plot(fitHW, lwd = 2, col = "blue", main = "Holt Winter Additive")

The original plot is in blue, the plotted graph is in red. Looks like the fitted model seems to be doing decent, and much better at the later points! Now let’s attempt to predict and plot the next 2 weeks net change in bikes.

### Fitted line is in red ###
fa <- forecast(fitHW, 14)
plot(fa, lwd = 2, col = "blue", main = "HW Additive Prediction")

accuracy(fa)

                     ME     RMSE      MAE  MPE MAPE      MASE       ACF1
Training set -0.6762827 10.93165 8.317779 -Inf  Inf 0.8831715 0.09600309

Box.test(residuals(fitHW), lag = 7, type = "Ljung")

    Box-Ljung test

data:  residuals(fitHW)
X-squared = 5.6411, df = 7, p-value = 0.5822

The predicted model seems to work well because the observed p-value is higher than 0.05 which is a good indicator!

Now let’s attempt to then use the Holt-Winter’s additive model on all the other stations to predict the net change in bikes for the next 2 weeks (i.e. 1st Sep 14 - 14th Sep 14). After which, the results will be stored into a table and then used for the tool.

remove(startdt, enddt, sfct)
### Create an empty data frame to store all the values ###
final <- data.frame()
for (i in 1:length(sf)) {
    station <- sf[i]
    dt <- tripdata[tripdata$Start.Station == station | tripdata$End.Station == station,]
    startdt <- select(dt, Start.Station, Start.DatewTime) %>% filter(Start.Station == station)
    enddt <- select(dt, End.Station, End.DatewTime) %>% filter(End.Station == station)
    ### Extracts only the date and hour out ###
    startdt$Start.DatewTime <- as.POSIXct(format(as.POSIXct(startdt$Start.DatewTime),"%Y-%m-%d"))
    enddt$End.DatewTime <- as.POSIXct(format(as.POSIXct(enddt$End.DatewTime),"%Y-%m-%d"))
    ### Aggregate the number of trips by date and hour ###
    startdt <- aggregate(Start.Station ~ Start.DatewTime, data = startdt, FUN = length)
    enddt <- aggregate(End.Station ~ End.DatewTime, data = enddt, FUN = length)
    ### Create a column of sequential date and time from 1st Mar 14 to 31st Aug 14 ###
    datetime <- seq(from = as.POSIXct("2014-03-01"), to = as.POSIXct("2014-08-31"), by = "day")
    datetime <- data.frame(datetime)
    ### Merge all the subsequent columns together, replace NA values with 0 ###
    combdt <- merge(datetime, startdt, by.x = "datetime", by.y = "Start.DatewTime", all.x = TRUE)
    combdt$Start.Station[is.na(combdt$Start.Station)] <- 0
    combdt <- merge(combdt, enddt, by.x = "datetime", by.y = "End.DatewTime", all.x = TRUE)
    combdt$End.Station[is.na(combdt$End.Station)] <- 0
    ### The start.station column signifies number of bikes leaving the station, the end.station column signifies the number of bikes entering the station. To get the net change, we form a new column of Net = End.Station - Start.Station ###
    combdt <- mutate(combdt, Net = End.Station - Start.Station)
    
    combdt <- select(combdt, datetime, Net)
    combdts <- ts(combdt$Net, frequency = 7, start = c(1,3))
    ### Fit an additive Holt Winter's Additive model ###
    fitHW <- HoltWinters(combdts, seasonal = "additive")
    
    ### Predict the values for the next 14 days ###
    fa <- forecast(fitHW, 14)
    
    predicted <- as.character(as.numeric(fa$mean))
    datetime <- seq(from = as.POSIXct("2014-09-01"), to = as.POSIXct("2014-09-14"), by = "day")
    datetime <- as.character(datetime)
    station_name <- rep(station, 14)
    
    overall <- data.frame(cbind(station_name, datetime, predicted))
    final <- rbind(final,overall)
}
final$station_name <- as.character(final$station_name)
final$datetime <- as.character(final$datetime)
final$predicted <- as.numeric(as.character(final$predicted))
write.csv(final, file = "./final.csv")

Now that we have the final table, that consists of the station name, the date and the predicted change in number of bikes at that particular station and date, we are ready to use this CSV file for our tool!

What this tool would then help to achieve, is giving the manager the ability to know in advance, for a particular day, which station would have the highest number of incoming and highest number of outgoing bikes. The manager will then be able to assign and allocate employees to these stations and balance the bikes in a more efficient manner!

Numbers in positive will represent a net increase in number of bikes. Numbers in negative will represent a net decrease in number of bikes.

E.g. On September 12th 2014, We can expect a huge number of bikes entering San Francisco Caltrain (Townsend at 4th) at 46.81, on the other hand, we notice that there will be a large number of bikes leaving Embarcadero at Bryant. Hence, we should have most number of employees at these 2 stations to facilitate the balance of number of bikes.

4 Limitations of tool

This tool would have been better enhanced if the number of bikes was first available at 1st March 14, 00:00:00 and the number of slots available at each station. This would then be able to give a more accurate indicator of whether this particular bike station would get filled up soon.

E.g. if there is 25 bikes in the station on this day, and it is predicted to have a net change of 25 bikes the next day. This would certainly be a problem if the number of bike slots is 50, because this would indicate that the bike station would be filled up by the end of next day if the prediction is right! However, this wouldn’t be the case if the number of bike slots is 75 instead.

In addition, the prediction model could have been better optimized to fit the data. This would then result in a more accurate number.

Lastly, this tool could also be more accurate if we were given additional inputs such as the temperature, humidity level, etc. because weather conditions can definitely affect the number of trips as well.