This analysis identifies SEPTA bus stops that are under- or over-performing relative to their predicted ridership based on service characteristics, demographics, and accessibility. By mapping the residuals from our spatial predictive models, we reveal neighborhoods with unmet transit demand and underutilized services, providing actionable recommendations for SEPTA bus service planning and investment priorities.
Key Findings:
Research Question: Where should SEPTA prioritize investments to maximize ridership and serve unmet transit demand?
Transit agencies operate under persistent budget, fleet, and staffing constraints. As a result, they must make deliberate choices about where to add service, upgrade stops, or scale back underperforming routes. Traditional performance metrics such as raw ridership counts or simple boardings-per-hour thresholds identify busy or quiet locations, but they do not answer a critical question: Is a stop performing better or worse than expected given its context?
A low-ridership stop in a low-density, auto-oriented neighborhood may be performing as expected, while a similarly low-ridership stop in a dense, transit-dependent corridor could signal unmet demand or barriers to access.
To move beyond raw counts, we use predictive modeling to estimate “expected” ridership at each stop based on service characteristics, neighborhood demographics, and accessibility. We then compare these predictions to actual ridership. Stops that overperform (higher ridership than predicted) and those that underperform (lower ridership than predicted) reveal where the network is over- or under-delivering relative to its potential.
This residuals-based approach allows us to:
Distinguish between structurally low-demand locations and places where transit is underperforming.
Identify corridors with latent demand that could benefit from targeted investments.
Highlight stops where ridership is strong despite barriers, suggesting supportive local conditions that can be replicated elsewhere.
For SEPTA, this analysis supports:
The remainder of the project builds a modeling and spatial analysis framework to quantify these patterns and translate them into planning guidance for SEPTA’s bus network.
We integrate four main categories of data: 1. Stop-level ridership 2. Transit service characteristics 3. Neighborhood socio-demographics 4. Accessibility and spatial context All spatial data are ultimately projected to EPSG:2272 (PA State Plane South) to ensure consistent distance and area calculations in the Philadelphia region.
# Set options
options(scipen = 999)
options(tigris_use_cache = TRUE)
# Packages
if(!require(pacman)){install.packages("pacman"); library(pacman, quietly = T)}
p_load(tidyverse, lubridate, sf, tigris, tidycensus, riem, viridis, gridExtra, knitr, kableExtra, tidytransit, stringr, spdep, corrplot, car, MASS, conflicted, caret, lmtest, hexbin, webshot2)
conflict_prefer("select", "dplyr")
conflict_prefer("filter", "dplyr")# Define consistent plot themes
plotTheme <- theme(
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 10),
plot.caption = element_text(size = 8),
axis.text.x = element_text(size = 10, angle = 45, hjust = 1),
axis.text.y = element_text(size = 10),
axis.title = element_text(size = 11, face = "bold"),
panel.background = element_blank(),
panel.grid.major = element_line(colour = "#D0D0D0", linewidth = 0.2),
panel.grid.minor = element_blank(),
axis.ticks = element_blank(),
legend.position = "right"
)
mapTheme <- theme(
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 10),
plot.caption = element_text(size = 8),
axis.line = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.title = element_blank(),
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_line(colour = 'transparent'),
panel.grid.minor = element_blank(),
legend.position = "right",
plot.margin = margin(1, 1, 1, 1, 'cm'),
legend.key.height = unit(1, "cm"),
legend.key.width = unit(0.2, "cm")
)We use SEPTA’s Fall 2024 bus ridership summary, which represents a recent but relatively stable service period that predates the most recent service reductions. This allows us to model ridership patterns under a consistent service regime.
We first rename and retain only the variables needed for analysis and then construct a set of derived ridership metrics:
Source: SEPTA Open Data Portal
bus_ridership_raw <- st_read("./data/septa_fall_2024_bus_summary.geojson", quiet = TRUE)# preview
glimpse(bus_ridership_raw)Rows: 18,290
Columns: 16
$ FID <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
$ Mode <chr> "Bus", "Bus", "Bus", "Bus", "Bus", "Bus", "Bus", "Bus", "Bus…
$ Route <chr> "18", "18", "18", "18", "18", "18", "18", "18", "18", "18", …
$ Direction <chr> "Eastbound", "Eastbound", "Eastbound", "Eastbound", "Eastbou…
$ Route_1 <chr> "370", "371", "372", "373", "375", "467", "15421", "15424", …
$ Stop <chr> "Upsal St & Mansfield Av", "Chew Av & Chelten Av", "Olney Av…
$ Lat <dbl> 40.06365, 40.04550, 40.03903, 40.03887, 40.05443, 40.07590, …
$ Lon <dbl> -75.16982, -75.16500, -75.14496, -75.14416, -75.09078, -75.0…
$ Sign_Up <chr> "Fall 2024", "Fall 2024", "Fall 2024", "Fall 2024", "Fall 20…
$ Wk_Ons <int> 18, 166, 668, 692, 9, 12, 228, 19, 36, 56, 48, 19, 24, 53, 8…
$ Wk_Offs <int> 19, 175, 1163, 144, 44, 47, 0, 0, 4, 31, 11, 3, 5, 19, 11, 8…
$ Sat_Ons <int> 13, 146, 419, 515, 4, 3, 181, 10, 19, 30, 18, 7, 9, 25, 22, …
$ Sat_Offs <int> 9, 87, 682, 85, 24, 31, 0, 2, 5, 15, 9, 6, 8, 9, 1, 0, 15, 9…
$ Sun_Ons <int> 12, 92, 266, 347, 4, 4, 126, 6, 16, 16, 16, 9, 10, 17, 12, 4…
$ Sun_Offs <int> 4, 64, 455, 72, 20, 21, 0, 2, 5, 12, 4, 3, 2, 6, 5, 0, 12, 5…
$ geometry <POINT [°]> POINT (-75.16982 40.06365), POINT (-75.165 40.0455), P…# clean data
bus_ridership <- bus_ridership_raw %>%
select(- Sign_Up, - Mode) %>%
rename(fid = FID,
route = Route,
direction = Direction,
stop_id = Route_1,
stop_name = Stop,
latitude = Lat,
longitude = Lon,
wk_ons = Wk_Ons,
wk_offs = Wk_Offs,
sat_ons = Sat_Ons,
sat_offs = Sat_Offs,
sun_ons = Sun_Ons,
sun_offs = Sun_Offs)
#aggregate ridership metrics and compute average daily boardings
bus_ridership <- bus_ridership %>%
mutate(total_boardings = wk_ons + sat_ons + sun_ons,
total_alightings = wk_offs + sat_offs + sun_offs,
total_ridership = total_boardings + total_alightings,
avg_daily_boardings = (wk_ons * 5 + sat_ons + sun_ons) / 7,
avg_daily_alightings = (wk_offs * 5 + sat_offs + sun_offs) / 7,
avg_weekday_boardings = wk_ons / 5,
weekend_share = (sat_ons + sun_ons) / (total_boardings + 0.001),
# categories based on average weekday boardings
ridership_category = case_when(
avg_weekday_boardings == 0 ~ "No Activity",
avg_weekday_boardings < 5 ~ "Very Low (< 5/day)",
avg_weekday_boardings < 20 ~ "Low (5-20/day)",
avg_weekday_boardings < 50 ~ "Medium (20-50/day)",
avg_weekday_boardings < 100 ~ "High (50-100/day)",
TRUE ~ "Very High (100+/day)"),
log_wk_boardings = log(wk_ons + 1),
log_total_boardings = log(total_boardings + 1))
glimpse(bus_ridership)Rows: 18,290
Columns: 24
$ fid <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 1…
$ route <chr> "18", "18", "18", "18", "18", "18", "18", "18", …
$ direction <chr> "Eastbound", "Eastbound", "Eastbound", "Eastboun…
$ stop_id <chr> "370", "371", "372", "373", "375", "467", "15421…
$ stop_name <chr> "Upsal St & Mansfield Av", "Chew Av & Chelten Av…
$ latitude <dbl> 40.06365, 40.04550, 40.03903, 40.03887, 40.05443…
$ longitude <dbl> -75.16982, -75.16500, -75.14496, -75.14416, -75.…
$ wk_ons <int> 18, 166, 668, 692, 9, 12, 228, 19, 36, 56, 48, 1…
$ wk_offs <int> 19, 175, 1163, 144, 44, 47, 0, 0, 4, 31, 11, 3, …
$ sat_ons <int> 13, 146, 419, 515, 4, 3, 181, 10, 19, 30, 18, 7,…
$ sat_offs <int> 9, 87, 682, 85, 24, 31, 0, 2, 5, 15, 9, 6, 8, 9,…
$ sun_ons <int> 12, 92, 266, 347, 4, 4, 126, 6, 16, 16, 16, 9, 1…
$ sun_offs <int> 4, 64, 455, 72, 20, 21, 0, 2, 5, 12, 4, 3, 2, 6,…
$ geometry <POINT [°]> POINT (-75.16982 40.06365), POINT (-75.165…
$ total_boardings <int> 43, 404, 1353, 1554, 17, 19, 535, 35, 71, 102, 8…
$ total_alightings <int> 32, 326, 2300, 301, 88, 99, 0, 4, 14, 58, 24, 12…
$ total_ridership <int> 75, 730, 3653, 1855, 105, 118, 535, 39, 85, 160,…
$ avg_daily_boardings <dbl> 16.428571, 152.571429, 575.000000, 617.428571, 7…
$ avg_daily_alightings <dbl> 15.4285714, 146.5714286, 993.1428571, 125.285714…
$ avg_weekday_boardings <dbl> 3.6, 33.2, 133.6, 138.4, 1.8, 2.4, 45.6, 3.8, 7.…
$ weekend_share <dbl> 0.5813818, 0.5891075, 0.5062820, 0.5546972, 0.47…
$ ridership_category <chr> "Very Low (< 5/day)", "Medium (20-50/day)", "Ver…
$ log_wk_boardings <dbl> 2.944439, 5.117994, 6.505784, 6.541030, 2.302585…
$ log_total_boardings <dbl> 3.784190, 6.003887, 7.210818, 7.349231, 2.890372…We verified that there were no negative or missing ridership values and retained all rows. We then clipped the dataset to stops within the City of Philadelphia using Census place boundaries and projected the data to EPSG:2272.
# download Philadelphia boundary
philly_boundary <- places(state = "PA", year = 2022) %>%
filter(NAME == "Philadelphia") %>%
st_transform(4326)
# retain only boundary
philly_boundary <- philly_boundary %>%
st_geometry() %>%
st_as_sf()
# clip stops to Philadelphia boundary
ridership_philly <- st_intersection(bus_ridership, philly_boundary)
# project to EPSG:2272 (PA State Plane South)
ridership_philly <- ridership_philly %>%
st_transform(2272)
glimpse(ridership_philly)Rows: 10,898
Columns: 24
$ fid <int> 1, 2, 3, 4, 5, 6, 8, 9, 10, 11, 12, 13, 14, 15, …
$ route <chr> "18", "18", "18", "18", "18", "18", "18", "18", …
$ direction <chr> "Eastbound", "Eastbound", "Eastbound", "Eastboun…
$ stop_id <chr> "370", "371", "372", "373", "375", "467", "15424…
$ stop_name <chr> "Upsal St & Mansfield Av", "Chew Av & Chelten Av…
$ latitude <dbl> 40.06365, 40.04550, 40.03903, 40.03887, 40.05443…
$ longitude <dbl> -75.16982, -75.16500, -75.14496, -75.14416, -75.…
$ wk_ons <int> 18, 166, 668, 692, 9, 12, 19, 36, 56, 48, 19, 24…
$ wk_offs <int> 19, 175, 1163, 144, 44, 47, 0, 4, 31, 11, 3, 5, …
$ sat_ons <int> 13, 146, 419, 515, 4, 3, 10, 19, 30, 18, 7, 9, 2…
$ sat_offs <int> 9, 87, 682, 85, 24, 31, 2, 5, 15, 9, 6, 8, 9, 1,…
$ sun_ons <int> 12, 92, 266, 347, 4, 4, 6, 16, 16, 16, 9, 10, 17…
$ sun_offs <int> 4, 64, 455, 72, 20, 21, 2, 5, 12, 4, 3, 2, 6, 5,…
$ total_boardings <int> 43, 404, 1353, 1554, 17, 19, 35, 71, 102, 82, 35…
$ total_alightings <int> 32, 326, 2300, 301, 88, 99, 4, 14, 58, 24, 12, 1…
$ total_ridership <int> 75, 730, 3653, 1855, 105, 118, 39, 85, 160, 106,…
$ avg_daily_boardings <dbl> 16.428571, 152.571429, 575.000000, 617.428571, 7…
$ avg_daily_alightings <dbl> 15.4285714, 146.5714286, 993.1428571, 125.285714…
$ avg_weekday_boardings <dbl> 3.6, 33.2, 133.6, 138.4, 1.8, 2.4, 3.8, 7.2, 11.…
$ weekend_share <dbl> 0.5813818, 0.5891075, 0.5062820, 0.5546972, 0.47…
$ ridership_category <chr> "Very Low (< 5/day)", "Medium (20-50/day)", "Ver…
$ log_wk_boardings <dbl> 2.944439, 5.117994, 6.505784, 6.541030, 2.302585…
$ log_total_boardings <dbl> 3.784190, 6.003887, 7.210818, 7.349231, 2.890372…
$ geometry <POINT [US_survey_foot]> POINT (2690583 276595.2), POI…# summary statistics
summary_stats <- ridership_philly %>%
st_drop_geometry() %>%
summarise(
n_stops = n(),
n_routes = n_distinct(route),
total_weekday_boardings = sum(wk_ons, na.rm = TRUE),
mean_weekday_boardings = mean(wk_ons, na.rm = TRUE),
median_weekday_boardings = median(wk_ons, na.rm = TRUE),
max_weekday_boardings = max(wk_ons, na.rm = TRUE),
pct_zero_ridership = sum(wk_ons == 0) / n() * 100
)
cat("\nSummary Statistics:\n")
Summary Statistics:table(summary_stats) %>%
as.data.frame() %>%
print() n_stops n_routes total_weekday_boardings mean_weekday_boardings
1 10898 94 370845 34.028720866214
median_weekday_boardings max_weekday_boardings pct_zero_ridership Freq
1 13 2566 11.6902183886952 1# Ridership by category
cat("\nStops by ridership category:\n")
Stops by ridership category:table(ridership_philly$ridership_category) %>%
as.data.frame() %>%
rename(Category = Var1, Count = Freq) %>%
mutate(Percent = round(Count / sum(Count) * 100, 1)) %>%
print() Category Count Percent
1 High (50-100/day) 105 1.0
2 Low (5-20/day) 2913 26.7
3 Medium (20-50/day) 650 6.0
4 No Activity 1274 11.7
5 Very High (100+/day) 51 0.5
6 Very Low (< 5/day) 5905 54.2The summary shows the distribution of weekday boardings and the share of stops in each ridership category.
# Map of stops by ridership category
p1 <- ggplot(ridership_philly) +
geom_sf(aes(color = ridership_category, size = wk_ons), alpha = 0.6) +
scale_color_viridis_d(option = "plasma", name = "Ridership") +
scale_size_continuous(range = c(0.5, 5), guide = "none") +
labs(title = "Bus Stop Ridership Across Philadelphia") +
theme_void()
p1
A map of stops by category provides an initial spatial sense of where activity is concentrated across the city.
Transit service characteristics come from SEPTA’s GTFS feed and describe how much service reaches each stop and how well-connected it is within the broader network.
Headways: The average headway (in minutes) of the route with shortest average headway at the stop. We calculated the minimum average headways for 5 timeframes
headwayMorn corresponds to the minimum average during the morning rush hour (4 - 6 am).
headwayRushAM corresponds to the minimum average during the morning rush hour (6:30 - :830 am).
headwayMidday corresponds to the minimum average during the morning rush hour (noon - 2 pm).
headwayRushPM corresponds to the minimum average during the morning rush hour (5 - 7:30 pm).
headwayNight corresponds to the minimum average during the morning rush hour (10 - 11:59 pm).
Connectivity: Important variables of connectivity were the number of routes at each stop as well as accessibility of subway/trolley and regional rail stations.
distRail: the distance (in miles) between bus stop X and the nearest subway or trolley stop
distRR: the distance (in miles) between bus stop X and the nearest regional rail stop
railConnection: If the nearest rail or trolley station is within .1 miles of bus station X
rrConnection: If the nearest rail or trolley station is within .1 miles of bus station X
Source: OpenDataPhilly, SEPTA Open Data Portal
# SEPTA bus stations GTFS
septaBus <- read_gtfs("data/google_bus.zip")To get the the number of routes at each station, we first created a column of lists of routes at each station (which is crucial later). Then, we simply counted how many commas were in the list, with there being n + 1 routes for every n commas in the list.
#Count routes by number of commas in string
stopsRoutes <- septaBus$route_stops %>%
group_by(stop_id) %>%
summarise(routes = list(unique(route_id)),
numRoutes = 1 + str_count(routes, ","))Using Zhanchao’s [work](https://github.com/zyang91/septa-cut-analysis/blob/main/scripts/analysis.r) as an example, we found station-level headways at 5 given times: morning, A.M. rush hour, mid-day, P.M. rush hour, and night.
#Get active service ids
activeIds <- septaBus$calendar_dates %>%
distinct(service_id) %>%
pull(service_id)
#Get headways of routes in sec
freqMorn <- septaBus %>%
get_route_frequency(
start_time = "04:00:00",
end_time = "06:00:00",
service_ids = activeIds
)
freqPeakAM <- septaBus %>%
get_route_frequency(
start_time = "06:30:00",
end_time = "08:30:00",
service_ids = activeIds
)
freqMidday <- septaBus %>%
get_route_frequency(
start_time = "12:00:00",
end_time = "14:00:00",
service_ids = activeIds
)
freqPeakPM <- septaBus %>%
get_route_frequency(
start_time = "17:00:00",
end_time = "19:30:00",
service_ids = activeIds
)
freqNight <- septaBus %>%
get_route_frequency(
start_time = "22:00:00",
end_time = "23:59:59",
service_ids = activeIds
)
#headways in minutes
freqMorn <- freqMorn %>%
mutate(avgHeadwayMin = mean_headways / 60)
freqPeakAM <- freqPeakAM %>%
mutate(avgHeadwayMin = mean_headways / 60)
freqPeakPM <- freqPeakPM %>%
mutate(avgHeadwayMin = mean_headways / 60)
freqMidday <- freqMidday %>%
mutate(avgHeadwayMin = mean_headways / 60)
freqNight <- freqNight %>%
mutate(avgHeadwayMin = mean_headways / 60)
#Make a combined headways by time for each route
routeHeadways <- freqPeakAM %>% select(route_id, rushAM = avgHeadwayMin)
routeHeadways <-routeHeadways %>%
left_join(freqPeakPM %>% select(route_id, rushPM = avgHeadwayMin),
by = "route_id")
routeHeadways <-routeHeadways %>%
left_join(freqMidday %>% select(route_id, midday = avgHeadwayMin),
by = "route_id")
routeHeadways <-routeHeadways %>%
left_join(freqNight %>% select(route_id, night = avgHeadwayMin),
by = "route_id")
routeHeadways <-routeHeadways %>%
left_join(freqMorn %>% select(route_id, morn = avgHeadwayMin),
by = "route_id")However, this results in the average headway for each route. We need the average headway at each station. For stations with multiple bus routes, we decided the minimum of the two would illustrate the connectivity more accurately.
#prepare datagrame for new coumns
stopsRoutes$headwayRushAM <- NA
stopsRoutes$headwayRushPM <- NA
stopsRoutes$headwayMidday <- NA
stopsRoutes$headwayNight <- NA
stopsRoutes$headwayMorn <- NA
#Initalize var and ensure all headways will be less than
minAM <- 9999
minPM <- 9999
minMidday <- 9999
minNight <- 9999
minMorn <- 9999
#iterate through all stops
for(i in 1:nrow(stopsRoutes)){
#reset to ensure new iteration's headways will compare against themselves
minAM <- 9999
minPM <- 9999
minMidday <- 9999
minNight <- 9999
minMorn <- 9999
#print(paste("i: ",i))
#iterate through all routes at stop
for(j in 1:length(stopsRoutes$routes[[i]])){
#print(j)
#check if route in route_id
if(stopsRoutes$routes[[i]][j] %in% routeHeadways$route_id){
#find min headway for time periods
if(minAM > routeHeadways$rushAM[routeHeadways$route_id == stopsRoutes$routes[[i]][j]] &
!is.na(routeHeadways$rushAM[routeHeadways$route_id == stopsRoutes$routes[[i]][j]])){
minAM <- routeHeadways$rushAM[routeHeadways$route_id == stopsRoutes$routes[[i]][j]]
}
if(minPM > routeHeadways$rushPM[routeHeadways$route_id == stopsRoutes$routes[[i]][j]] &
!is.na(routeHeadways$rushPM[routeHeadways$route_id == stopsRoutes$routes[[i]][j]])){
minPM <- routeHeadways$rushPM[routeHeadways$route_id == stopsRoutes$routes[[i]][j]]
}
if(minMidday > routeHeadways$midday[routeHeadways$route_id == stopsRoutes$routes[[i]][j]] &
!is.na(routeHeadways$midday[routeHeadways$route_id == stopsRoutes$routes[[i]][j]])){
minMidday <- routeHeadways$midday[routeHeadways$route_id == stopsRoutes$routes[[i]][j]]
}
if(minNight > routeHeadways$night[routeHeadways$route_id == stopsRoutes$routes[[i]][j]] &
!is.na(routeHeadways$night[routeHeadways$route_id == stopsRoutes$routes[[i]][j]])){
minNight <- routeHeadways$night[routeHeadways$route_id == stopsRoutes$routes[[i]][j]]
}
if(minMorn > routeHeadways$morn[routeHeadways$route_id == stopsRoutes$routes[[i]][j]] &
!is.na(routeHeadways$morn[routeHeadways$route_id == stopsRoutes$routes[[i]][j]])){
minMorn <- routeHeadways$morn[routeHeadways$route_id == stopsRoutes$routes[[i]][j]]
}
}
}
#Set min timezone headway
stopsRoutes$headwayRushAM[i] <- minAM
stopsRoutes$headwayRushPM[i] <- minPM
stopsRoutes$headwayMidday[i] <- minMidday
stopsRoutes$headwayNight[i] <- minNight
stopsRoutes$headwayMorn[i] <- minMorn
}
#Switch 9999s for NAs
stopsRoutes <- stopsRoutes %>%
mutate(
headwayRushPM = case_when(headwayRushPM == 9999 ~ NA,
TRUE ~ headwayRushPM),
headwayRushAM = case_when(headwayRushAM == 9999 ~ NA,
TRUE ~ headwayRushAM),
headwayMorn = case_when(headwayMorn == 9999 ~ NA,
TRUE ~ headwayMorn),
headwayNight = case_when(headwayNight == 9999 ~ NA,
TRUE ~ headwayNight),
headwayMidday = case_when(headwayMidday == 9999 ~ NA,
TRUE ~ headwayMidday),
)First, we downloaded the station locations for trolley, subway, and regional rail stations.Next, we prepared the bus stop data for spatial analytics in the context of the other stations.
# bus, rail, & trolley stations
allStations <- read_sf("data/Transit_Stops/Transit_Stops_(Spring_2025).shp")
#Regional rail stations
rrStations <- read_sf("data/Regional_Rail_Stations/Regional_Rail_Stations.shp")
#Rail and trolly line abbr
railTrolley <- c("T1", "T2", "T3", "T4", "T5",
"G2", "D1", "D2", "MFL", "BSL",
"M1")
#trim down to strictly rail/trolley stations
railStations <- allStations %>%
filter(LineAbbr %in% railTrolley)
#Get locations for bus stops
#GTFS for consistency
busStops <- septaBus$stops %>%
select(stop_id, stop_name, wheelchair_boarding, stop_lon, stop_lat)
#combining headways with wheelchairs and coords
stopsRoutes2 <- stopsRoutes %>% left_join(busStops, by = "stop_id")
stopsRoutes2 <- stopsRoutes2 %>%
filter(!is.na(stop_lon) & !is.na(stop_lat))
#get the SFs onto same projection, must go through WGS first
stopsRoutes2 <- st_as_sf(stopsRoutes2,
coords = c("stop_lon","stop_lat"),
crs = 4326)
stopsRoutes2 <- st_transform(stopsRoutes2, 3857)Now, with the data consistent, we found the distance to the nearest station (rail/trolley and regional) in miles.
#add column of distance to rail in miles, with a boolean if it is close
#add equivalent columns for regional rail
stopsRoutes2 <- stopsRoutes2 %>%
mutate(
distRail = apply(st_distance(geometry,railStations$geometry),
1,
min) * 0.00062137,
railConnection = distRail < .1,
distRR = apply(st_distance(stopsRoutes2,
rrStations),
1,min)* 0.00062137,
rrConnection = distRR < .1
)We examined the created data to see how distributions and patterns may affect future analysis.
#Create a pivot
pivotHeadways <- stopsRoutes2 %>%
select(stop_id, headwayMorn, headwayRushAM,
headwayMidday, headwayRushPM,
headwayNight, rrConnection, railConnection) %>%
pivot_longer(headwayMorn:headwayNight, names_to = "time",
values_to = "avgHeadway") %>%
mutate(time = case_when(
time == "headwayRushAM" ~ "Morning Rush Hour",
time == "headwayRushPM" ~ "Evening Rush Hour",
time == "headwayMidday" ~ "Mid-day",
time == "headwayNight" ~ "Night",
time == "headwayMorn" ~ "Morning"
)) %>%
filter(!is.na(avgHeadway))
#Histogram of all station headways
ggplot(data = pivotHeadways, aes(x = avgHeadway,
y = ..density..)) +
geom_histogram(bins = 50, fill = "#9ACD32", alpha = .7) +
geom_density() +
facet_wrap(vars(time), ncol = 3) +
xlab("Minimum Average Headway (min)") + ylab("Density") +
ggtitle("Minimum Average Headway For Bus Stations")+
plotTheme
#of stations where regional rail stops
ggplot(data = pivotHeadways[pivotHeadways$rrConnection,], aes(x = avgHeadway,
y = ..density..)) +
geom_histogram(bins = 50, fill = "#9ACD32", alpha = .7) +
geom_density() +
facet_wrap(vars(time), ncol = 3) +
xlab("Minimum Average Headway (min)") + ylab("Density") +
ggtitle("Minimum Average Headway For Bus Stations with Regional Rail\nConnections")+
plotTheme
#of stations where heavy rail/trolley stops
ggplot(data = pivotHeadways[pivotHeadways$railConnection,], aes(x = avgHeadway,
y = ..density..)) +
geom_histogram(bins = 50, fill = "#9ACD32", alpha = .7) +
geom_density() +
facet_wrap(vars(time), ncol = 3) +
xlab("Minimum Average Headway (min)") + ylab("Density") +
ggtitle("Minimum Average Headway For Bus Stations with Heavy rail\nor Trolly Connections")+
plotTheme
#Histogram of distance to nearest rail/trolley station
ggplot(data = stopsRoutes2, aes(x = distRail,
y = ..density..)) +
geom_histogram(bins = 50, fill = "#9ACD32", alpha = .7) +
xlab("Distance (miles)") + ylab("Density") +
ggtitle("Distance to Nearest Subway/Trolley Station")+
plotTheme
#to nearest regional rail station
ggplot(data = stopsRoutes2, aes(x = distRR,
y = ..density..)) +
geom_histogram() +
xlab("Distance (miles)") + ylab("Density") +
ggtitle("Distance to Nearest Regional Rail Station") +
plotTheme
Both distances to nearest regional rail and subway/trolley stations are right-skewed, with a majority of stations relatively close to 0 and creating an increasingly small tail as distance increases. The headway distributions are shaped a little more complexly. Both the evening and morning rush hour have a general right skew, but not as steadily as the distance variables; the mid-day headways are ever so slightly right-skewed, but the pattern is also close to be constant; finally the morning and evening headway distribution is almost completely constant. The right-skewed distributions may require a logarithmic transformation for normalization purposes.
We retrieved demographic and socioeconomic data from the American Community Survey (ACS) 5-year estimates for 2022 at the census tract level. Our variable selection included factors potentially associated with transit ridership.
Demographics & Economics: - Total population (pop_tot), median age (med_age), median household income (med_hh_inc), and employment status (employment_status)
Commuting Behavior: - Total workers by means of transportation (means_of_trans_tot), workers with no vehicle access (means_of_trans_no_veh), solo drivers (means_of_trans_car_alone), and public transit users (means_of_trans_publictransit)
Housing Tenure: - Total occupied units (tenure_total), owner-occupied (owner_occupied), and renter-occupied (renter_occupied) housing units
We retrieved data for all Pennsylvania census tracts (FIPS code 42) and transformed the geometry to Pennsylvania State Plane South (EPSG:2272), which uses feet as the unit of measurement and is appropriate for accurate distance calculations in the Philadelphia region.
Source: American Community Survey (ACS) 2022 5-year estimate
#census_api_key("YOUR_Key_Here", overwrite = TRUE, install = TRUE)census_api_key(Sys.getenv("MY_CENSUS_API_KEY"))
vars <- c("pop_tot" = "B01001_001",
"med_hh_inc" = "B19013_001",
"med_age" = "B01002_001",
"means_of_trans_tot" = "B08141_001",
"means_of_trans_no_veh" = "B08141_002", # Going to work with no access to a vehicle
"means_of_trans_car_alone" = "B08141_006", # Going to work alone in car
"means_of_trans_publictransit" = "B08141_016", # Going to work via public transit
"employment_status" = "B23025_001",
"tenure_total" = "B25003_001", # Occupied housing units
"owner_occupied" = "B25003_002",
"renter_occupied" = "B25003_003")
fips_pa <- 42
# Retrieve data from the ACS for 2022
acs_data_pa <- get_acs(geography = "tract",
variable = vars,
year = 2022,
state = fips_pa,
output = "wide",
geometry = T,
progress_bar = F) %>%
st_transform(2272)We cleaned the raw ACS data by creating separate columns for tracts, counties, and states, then filtered to retain only Philadelphia County tracts. Afterwards, we engineered several features from the raw census variables to better capture neighborhood characteristics relevant to transit ridership:
Density & Scale: - Population density (people per acre)
Commuting Patterns & Transit Dependency: - Percent of workers with no vehicle access (pct_no_vehicle) - Percent driving alone (pct_car_alone) - Percent using public transit (pct_public_transit) - Percent of population that commutes to work (pct_commuters)
Demographic Indicators: - Binary flag for “young” areas (median age < 35 years) - Percent renter-occupied (pct_renter) and owner-occupied (pct_owner) housing units
The percentage-based features normalize raw counts by tract, allowing for more meaningful comparisons across neighborhoods with varying populations.
# Separate NAME column into its constituent parts
acs_data_pa <- acs_data_pa %>%
separate(NAME,
into = c("TRACT", "COUNTY", "STATE"),
sep = "; ",
remove = T) %>%
mutate(TRACT = parse_number(TRACT),
COUNTY = sub(x = COUNTY, " County", ""))
# Filter out Philadelphia tracts
acs_data_phl <- acs_data_pa %>%
filter(COUNTY == "Philadelphia")
# Creating a new dataframe with engineered features, to keep the original census derived data clean.
acs_data_phl_features <- acs_data_phl %>%
mutate(
pop_density = pop_totE / (st_area(geometry) %>% as.numeric() / 43560), # per acre
# Percentages based on the transportation to work metric.
pct_no_vehicle = (means_of_trans_no_vehE / means_of_trans_totE) * 100,
pct_car_alone = (means_of_trans_car_aloneE / means_of_trans_totE) * 100,
pct_public_transit = (means_of_trans_publictransitE / means_of_trans_totE) * 100,
# Originally had a cutoff of 44 (max age of Millenials) but lowered the threshold to cover skewing that could occur with medians.
is_young_area = ifelse(med_ageE < 35, 1, 0),
# How many people in a tract are actually going to a workplace or are working (means of trans includes WFH)
pct_commuters = (means_of_trans_totE / pop_totE) * 100,
# To see if tracts have higher renter rates or home ownership rates
pct_renter = (renter_occupiedE / tenure_totalE) * 100,
pct_owner = (owner_occupiedE / tenure_totalE) * 100
)We examined the distributions of our engineered features to identify potential skewness, outliers, and patterns that may require transformation or special handling during model development.
hist_data <- acs_data_phl_features %>%
st_drop_geometry() %>%
select(pop_density, pct_no_vehicle, pct_car_alone, pct_public_transit,
pct_commuters, pct_renter, pct_owner, med_ageE, med_hh_incE, pop_totE) %>%
pivot_longer(everything(), names_to = "variable", values_to = "value")
ggplot(hist_data, aes(x = value)) +
geom_histogram(bins = 30, fill = "steelblue", color = "white") +
facet_wrap(~variable, scales = "free") +
theme_minimal() +
labs(title = "Distributions of Numeric Features Engineered from Census Data",
x = "Value", y = "Count")
Here, we see patterns that must be accounted for during model development. Percent of workers with no vehicle access (pct_no_vehicle) and median household income (med_hh_inc) both show right skewness, with most tracts at lower values and a long tail of high values. Population density (pop_density) appears relatively symmetric but contains potential outliers. These distribution characteristics may require transformations (e.g. a log transform) or handling of outliers during model development.
We incorporated worker-weighted transit job accessibility from the University of Minnesota’s data (DRUM), which reports the average number of jobs reachable from each tract by walk/roll plus transit within travel-time thresholds from 5 to 60 minutes.
access_jobs_raw <- read.csv("data/UMN_Accessibility_Pennsylvania/Pennsylvania_42_transit_census_tract_2022.csv")From this table we extractd, for each census tract:
jobs_15 – average jobs reachable within 15 minutes by transitjobs_30 – average jobs reachable within 30 minutes by transitWe reshaped the data from long (one row per tract–threshold pair) to wide form with separate columns for jobs_15 and jobs_30, enabling a straightforward join to ACS tracts and, subsequently, to stops.
access_jobs <- access_jobs_raw %>%
select (Census.ID , Threshold, Weighted_average_total_jobs) %>%
rename(census_tract = Census.ID)
access_jobs <- access_jobs %>%
filter(Threshold %in% c(15, 30)) %>%
pivot_wider(
id_cols = census_tract,
names_from = Threshold,
values_from = Weighted_average_total_jobs,
names_prefix = "jobs_"
)To capture the monocentric structure of Philadelphia, we compute distance from each bus stop to City Hall, treated as the central business district (CBD):
distance_cbd).This measure allowed us to test whether ridership systematically declines with increasing distance from the core.
# city Hall coordinates in WGS84
city_hall_wgs <- st_sfc(
st_point(c(-75.1636, 39.9526)),
crs = 4326)
# transform City Hall to EPSG:2272
city_hall <- st_transform(city_hall_wgs, 2272)
ridership_philly <- st_transform(ridership_philly , 2272)
# distance to CBD in miles
ridership_philly <- ridership_philly %>%
mutate(
distance_cbd_m = st_distance(geometry, city_hall) %>%
as.numeric(), # in feet (for EPSG:2272)
distance_cbd = distance_cbd_m / 5280 # convert to miles
)We used neighborhood boundaries from OpenDataPhilly to assign each stop a neighborhood name (neigh_name). After transforming neighborhood polygons to match the ridership CRS, we performed a spatial join so that each stop inherits the neighborhood in which it falls. This variable served both as a descriptive unit for mapping and as the basis for neighborhood-level clustering in Section 4.
neighborhoods <- st_read("data/philadelphia_neighborhoods/philadelphia_neighborhoods.shp")Reading layer `philadelphia_neighborhoods' from data source
`/Users/sujankakumanu/Documents/musa-5080-portfolio/assignments/final/data/philadelphia_neighborhoods/philadelphia_neighborhoods.shp'
using driver `ESRI Shapefile'
Simple feature collection with 159 features and 5 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -75.28026 ymin: 39.86701 xmax: -74.95576 ymax: 40.13799
Geodetic CRS: WGS 84neighborhoods <- st_transform(neighborhoods, st_crs(ridership_philly))
# join neighborhood names
neighborhoods_clean <- neighborhoods %>%
select(neigh_name = NAME, geometry)
ridership_philly <- ridership_philly %>%
st_join(neighborhoods_clean)Our exploratory analysis examined the distribution of ridership, its spatial patterns, and its relationships with key predictors. This helped us diagnose transformations, identify spatial clustering, and anticipate potential multicollinearity before formal model building.
Average weekday boardings are highly right-skewed: most stops had relatively modest activity, while a small number of locations served very high volumes of riders. A histogram of avg_weekday_boardings made this skew clear, and summary statistics confirmed the long upper tail.
Applying a log transformation (log(avg_weekday_boardings + 1)) yielded a distribution that was much closer to symmetric and suitable for linear modeling. We therefore used log-transformed ridership as the dependent variable in subsequent regression models.
# Ridership distribution
p1 <- ggplot(ridership_philly, aes(x = avg_weekday_boardings)) +
geom_histogram(bins = 50, fill = "#440154FF", alpha = 0.7) +
geom_vline(aes(xintercept = median(avg_weekday_boardings)),
color = "#FDE725FF", linetype = "dashed", size = 1) +
labs(title = "Distribution of Average Weekday Boardings",
subtitle = "Right Skew",
x = "Average Weekday Boardings )",
y = "Count of Stations") +
theme_minimal()
# Log-transformed ridership
p2 <- ggplot(ridership_philly, aes(x = log(avg_weekday_boardings + 1))) +
geom_histogram(bins = 50, fill = "#440154FF", alpha = 0.7) +
geom_vline(aes(xintercept = median(log(avg_weekday_boardings + 1))),
color = "#FDE725FF", linetype = "dashed", size = 1) +
labs(title = "Distribution of Log(Average Weekday Boardings)",
subtitle = "More normal after log transformation",
x = "Log(Average Weekday Boardings + 1)",
y = "Count of Stations") +
theme_minimal()
p1
ggsave("presentation-assets/eda_p1.png",
width = 10, height = 6, dpi = 300)
p2
ggsave("presentation-assets/eda_p2.png",
width = 10, height = 6, dpi = 300)
# Summary statistics
summary(ridership_philly$avg_weekday_boardings) %>%
as.matrix() %>%
t() %>%
as.data.frame() %>%
kable(caption = "Ridership Summary Statistics", digits = 1) %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Min. | 1st Qu. | Median | Mean | 3rd Qu. | Max. |
|---|---|---|---|---|---|
| 0 | 0.6 | 2.6 | 6.8 | 7.2 | 513.2 |
Mapping stops with point size and color scaled to avg_weekday_boardings revealed the macro-level geography of bus ridership in Philadelphia. High-ridership stops cluster along major corridors such as arterial streets, key transfer hubs, and downtown-adjacent neighborhoods, while low-ridership stops are more dispersed along peripheral and lower-density segments.
This map provided important context for later residual analysis. Many of the highest-demand locations already align with expected transit corridors, but there are also pockets of strong ridership in outlying neighborhoods that merit closer attention.
# Map of ridership by stop
ggplot() +
geom_sf(data = ridership_philly,
aes(size = avg_weekday_boardings, color = avg_weekday_boardings),
alpha = 0.6) +
scale_color_viridis(option = "plasma", name = "Avg Weekday\nRidership",
trans = "log10") +
scale_size_continuous(name = "Avg Weekday\nRidership",
range = c(0.5, 2),
trans = "log10") +
labs(title = "SEPTA Bus Ridership Across the Network",
subtitle = "Size and color show average weekday boardings (log scale)") +
theme_void() +
theme(plot.title = element_text(face = "bold", size = 14),
legend.position = "right")
To identify spatial clusters of high and low ridership, we conducted a Local Indicators of Spatial Association (LISA) analysis using Local Moran’s I. This technique reveals whether bus stops with similar ridership levels cluster together geographically.
Methodology:
We log-transformed average weekday boardings to address the right-skewed distribution identified in Section 3.1. For each bus stop, we identified its 5 nearest neighbors and created a spatial weights matrix using row-standardized weights. We then calculated Local Moran’s I statistics and their significance levels (p < 0.05).
Cluster Classification:
The resulting map identifies corridors of high ridership (High-High) and areas with consistently low ridership (Low-Low), providing spatial context for understanding where ridership patterns deviate from our model predictions in later sections.
# Log transforming boardings because above, we saw that it's skewed
ridership_philly$log_boardings <- log(ridership_philly$avg_weekday_boardings + 1)
coords <- st_coordinates(ridership_philly)
# Set-up for LISA, with 5 nearest neighbors and a weight matrix.
knn_nb <- knearneigh(coords, k = 5)
nb <- knn2nb(knn_nb)
lw <- nb2listw(nb, style = "W")
# Local Moran's I
lisa <- localmoran(ridership_philly$log_boardings, lw)
# Adding results back to the ridership data
ridership_philly$Ii <- lisa[, "Ii"]
ridership_philly$z_Ii <- lisa[, "Z.Ii"]
ridership_philly$p_Ii <- lisa[, "Pr(z != E(Ii))"]mean_val <- mean(ridership_philly$log_boardings, na.rm = TRUE)
ridership_philly$quad_sig <- "Not significant"
sig_level <- 0.05
# This loop determines how to categorize values on the map. We only care about significant values.
# Then, we compare each value to the mean to determine if it's high or low. Then we can compare to the Moran's I value to see what the neighbor values are.
for (i in 1:nrow(ridership_philly)) {
if (ridership_philly$p_Ii[i] < sig_level) {
if (ridership_philly$log_boardings[i] > mean_val & ridership_philly$Ii[i] > 0) {
ridership_philly$quad_sig[i] <- "High-High"
} else if (ridership_philly$log_boardings[i] < mean_val & ridership_philly$Ii[i] > 0) {
ridership_philly$quad_sig[i] <- "Low-Low"
} else if (ridership_philly$log_boardings[i] > mean_val & ridership_philly$Ii[i] < 0) {
ridership_philly$quad_sig[i] <- "High-Low"
} else if (ridership_philly$log_boardings[i] < mean_val & ridership_philly$Ii[i] < 0) {
ridership_philly$quad_sig[i] <- "Low-High"
}
}
}ggplot(ridership_philly) +
geom_sf(aes(color = quad_sig), size = 2, alpha = 0.8) +
scale_color_manual(values = c(
"High-High" = "red",
"Low-Low" = "blue",
"High-Low" = "lightpink",
"Low-High" = "lightblue",
"Not significant" = "grey70"
)) +
labs(
title = "LISA (Local Moran's I) - Bus Ridership Hotspots/Coldspots",
color = "Cluster Type"
) +
theme_void() +
theme(
legend.position = "right",
plot.title = element_text(hjust = 0.5, size = 16)
)
Before feature engineering and model development, we conducted correlation analysis to examine relationships between potential predictors and our dependent variable. This preliminary analysis identifies strong predictors and flags potential multicolinearity concerns, which will be formally assessed using variance inflation factors (VIF) during the modeling phase.
# Join all predictor datasets
# Spatial join with census data with features
ridership_full <- st_join(
ridership_philly,
acs_data_phl_features %>%
select(GEOID, pop_density, pct_no_vehicle, pct_car_alone,
pct_public_transit, pct_commuters, pct_renter, pct_owner,
med_ageE, med_hh_incE)
) %>%
st_drop_geometry() %>%
# Joining with the route and stop level features
left_join(stopsRoutes2 %>% st_drop_geometry() %>%
select(stop_id, numRoutes, headwayMorn, headwayRushAM,
headwayMidday, headwayRushPM, headwayNight,
distRail, distRR),
by = "stop_id") %>%
left_join(access_jobs %>% mutate(GEOID = as.character(census_tract)),
by = "GEOID")
# Only selecting the predictors from the joined data above, to make the correlation easier
cor_data <- ridership_full %>%
select(avg_weekday_boardings, numRoutes, headwayMorn, headwayRushAM,
headwayMidday, headwayRushPM, headwayNight, distRail, distRR,
pop_density, pct_no_vehicle, pct_car_alone, pct_public_transit,
pct_commuters, pct_renter, pct_owner, med_ageE, med_hh_incE,
jobs_15, jobs_30, distance_cbd) %>%
na.omit()
# Calculate and plot correlation matrix
cor_matrix <- cor(cor_data)
corrplot(cor_matrix,
method = "color",
tl.srt = 45,
tl.col = "black",
addCoef.col = "black",
number.cex = 0.8, # Making the correlation numbers a bit smaller
col = colorRampPalette(c("#6D9EC1", "white", "#E46726"))(200),
title = "Correlation Matrix of Predictors",
mar = c(0, 0, 2, 0)) # margins for readability
What defines “strong” correlation varies by context, particularly in transportation planning where built environment and behavioral factors can have varying levels of association. In addition, we will be doing further feature engineering in Section 4 that will adjust these associations. For a more easily interpretable view of potential multicollinearity, we can summarize variables with correlations where |r| > 0.7.
# Identify highly correlated predictor pairs
high_cor_pairs <- cor_matrix %>%
as.data.frame() %>%
rownames_to_column("var1") %>%
pivot_longer(-var1, names_to = "var2", values_to = "correlation") %>%
filter(var1 < var2, # Remove duplicate pairs
var1 != "avg_weekday_boardings", # Exclude outcome variable
var2 != "avg_weekday_boardings",
abs(correlation) > 0.7) %>% # Flag high correlations
arrange(desc(abs(correlation)))
# Display problematic pairs
kable(high_cor_pairs,
caption = "Highly Correlated Predictor Pairs (|r| > 0.7)",
digits = 3) %>%
kable_styling(bootstrap_options = c("striped", "hover"))| var1 | var2 | correlation |
|---|---|---|
| pct_owner | pct_renter | -1.000 |
| headwayMidday | headwayRushPM | 0.866 |
| distance_cbd | distRail | 0.854 |
| distance_cbd | jobs_30 | -0.770 |
| jobs_30 | pct_car_alone | -0.747 |
| pct_car_alone | pct_no_vehicle | -0.747 |
| headwayRushAM | headwayRushPM | 0.746 |
| distance_cbd | pct_car_alone | 0.742 |
| headwayMidday | headwayRushAM | 0.732 |
| headwayNight | headwayRushPM | 0.730 |
| jobs_15 | jobs_30 | 0.721 |
| headwayNight | headwayRushAM | 0.713 |
| med_hh_incE | pct_commuters | 0.712 |
Below, we present some scatterplots of select predictors with respect to our dependent variable of avg_weekday_boardings. We give particular attention to predictors that are candidates for interaction terms, such as Job Accessibility in 30 mins, and Distance to Rail Stations.
p2 <- ggplot(ridership_full, aes(x = pct_no_vehicle, y = avg_weekday_boardings)) +
geom_point(alpha = 0.4, color = "gray") +
geom_smooth(method = "lm", se = TRUE, color = "maroon") +
labs(title = "% No Vehicle Access vs Ridership",
x = "% Workers with No Vehicle", y = "Avg Weekday Boardings") +
plotTheme
p3 <- ggplot(ridership_full, aes(x = jobs_30, y = avg_weekday_boardings)) +
geom_point(alpha = 0.4, color = "gray") +
geom_smooth(method = "lm", se = TRUE, color = "maroon") +
labs(title = "Job Accessibility vs Ridership",
x = "Jobs Reachable in 30 min", y = "Avg Weekday Boardings") +
plotTheme
p4 <- ggplot(ridership_full, aes(x = distRail, y = avg_weekday_boardings)) +
geom_point(alpha = 0.4, color = "gray") +
geom_smooth(method = "lm", se = TRUE, color = "maroon") +
labs(title = "Distance to Rail vs Ridership",
x = "Distance to Rail Station (miles)", y = "Avg Weekday Boardings") +
plotTheme
p5 <- ggplot(ridership_full, aes(x = pct_renter, y = avg_weekday_boardings)) +
geom_point(alpha = 0.4, color = "gray") +
geom_smooth(method = "lm", se = TRUE, color = "maroon") +
labs(title = "% Renter-Occupied vs Ridership",
x = "% Renter-Occupied Housing", y = "Avg Weekday Boardings") +
plotTheme
# Arrange plots in grid
grid.arrange(p2, p3, p4, p5, ncol = 2)
As you can see, our variables present very weak correlations with the dependent variable (something that was also shown in the correlation matrix above). We suspect this will impact the results of our OLS models below, and lend to non-linear models performing better in our analysis.
EDA revealed three key themes that shaped our modeling choices:
Ridership and several predictors are strongly right-skewed We addressed this by using log transformations for ridership, vehicle access, income, and selected headway measures to better satisfy linear model assumptions.
There is clear spatial clustering in ridership High–high and low–low clusters suggested that unobserved neighborhood and corridor factors influence ridership. We incorporated this spatial structure through neighborhood clusters, distance measures, and later residual mapping.
Some predictors are highly correlated We handled potential multicollinearity by avoiding redundant variables, constructing interaction terms, and checking variance inflation factors (VIF) in our models.
These insights motivated the feature engineering strategies described next.
In prior sections, we identified some predictors that have skewed distributions. Others show high correlations with each other or have conventionally accepted associations, suggesting opportunities for combining variables or creating composite measures. These are all examples of predictors where some transformations, such as interaction terms and log transformation, can create more useful predictors in our final model. This work is done below.
Below, we log-transform predictors that were previously discussed to have right-skewed distributions. This is done to make the distribution more normal, aiding in model assumptions - especially for OLS.
pct_no_vehicle)med_hh_inc)headwayRushAM)headwayRushPM)distRail)ridership_full <- ridership_full %>%
mutate(
log_pct_no_vehicle = log(pct_no_vehicle + 1),
log_med_hh_inc = log(med_hh_incE + 1),
log_headwayRushAM = log(headwayRushAM + 1),
log_headwayRushPM = log(headwayRushPM + 1),
log_distRail = log(distRail + 1),
log_avg_weekday_boardings = log(avg_weekday_boardings + 1)
)We create four interaction terms to test whether the effect of one predictor depends on the level of another. First, we create a new predictor called proximity_rail, which is the inverse of distance from bus stops to rail stations. This allows us to create interactions between it and predictors that have a positive correlation with stop ridership, ensuring consistent directional interpretation.
numRoutes × pct_no_vehicle: Tests whether having more route options affects ridership more where less workers have access to vehicles.jobs_30 × pct_no_vehicle: Tests whether job accessibility within 30 minutes matters more where less workers have access to vehicles.proximity_rail × jobs_30: Tests whether job accessibility has a stronger ridership effect at stops closer to rail stations.proximity_rail × pct_renter: Tests whether rail station proximity influences ridership more in renter-heavy areas.ridership_full <- ridership_full %>%
mutate(
proximity_rail = 1 / (distRail + 0.1), # 0.1 to avoid div by 0
routes_x_noveh = numRoutes * pct_no_vehicle,
jobs_x_noveh = jobs_30 * pct_no_vehicle,
railprox_x_jobs = proximity_rail * jobs_30,
railprox_x_renter = proximity_rail * pct_renter
)To capture neighborhood-level context beyond individual census or stop attributes, we created a categorical neighborhood cluster variable using k-means clustering. We selected three features that represent distinct dimensions of neighborhood socioeconomic and demand characteristics: average weekday boardings, median household income, and population density. These variables were combined into a single dataset neigh_for_clust and rows with missing values were removed to ensure model stability.
Before running the algorithm, the numeric variables were standardized using scale() so that income, density, and ridership contribute equally regardless of their units or magnitude. We used k = 5 clusters, to ensure meaningful variation across Philadelphia neighborhoods. K-means was run with 25 random starts to improve convergence reliability.
The resulting cluster assignments were then joined back to the neighborhoods shapefile to create a spatial layer of neighborhood types, and also merged into the main ridership_full dataset so each bus stop inherits the cluster of the neighborhood it belongs to. This cluster variable captures broader neighborhood context patterns (e.g., high-density/high-ridership, low-income/low-ridership areas) and helps the predictive model incorporate spatial and socioeconomic structure not captured by single variables alone.
# Build neighborhood-level table (one row per neigh_name)
neigh_for_clust <- ridership_full %>%
group_by(neigh_name) %>%
summarise(
avg_weekday_boardings = mean(avg_weekday_boardings, na.rm = TRUE),
med_hh_incE = median(med_hh_incE, na.rm = TRUE),
pop_density = mean(pop_density, na.rm = TRUE),
.groups = "drop"
) %>%
na.omit()
set.seed(5080)
k_value <- 5
km_neigh <- kmeans(
scale(neigh_for_clust %>% select(-neigh_name)),
centers = k_value,
nstart = 25
)
# attach cluster labels to neighborhood summary table
neigh_clusters <- neigh_for_clust %>%
mutate(
neighborhood_cluster = factor(km_neigh$cluster)
)
# join clusters to neighborhood sf
neighborhoods_clusters <- neighborhoods_clean %>%
left_join(neigh_clusters, by = "neigh_name")
# Attach cluster to stops
ridership_full <- ridership_full %>%
left_join(
neigh_clusters %>% select(neigh_name, neighborhood_cluster),
by = "neigh_name"
)Stops map colored by neighborhood cluster
#attaching cluster to ridership_philly because it has the stop geometry
ridership_philly <- ridership_philly %>%
left_join(
neigh_clusters %>% select(neigh_name, neighborhood_cluster),
by = "neigh_name"
) %>%
filter(!is.na(neighborhood_cluster))
neighborhoods_clusters <- st_transform(
neighborhoods_clusters,
st_crs(ridership_philly)
)
ggplot() +
geom_sf(data = neighborhoods_clusters,
fill = "gray95", color = "white", linewidth = 0.2) +
geom_sf(data = ridership_philly,
aes(color = neighborhood_cluster),
size = 0.7, alpha = 0.7) +
scale_color_brewer(palette = "Set1", name = "Neighborhood\nCluster") +
labs(
title = "Bus Stops by Neighborhood Cluster",
) +
theme_void() +
theme(
legend.position = "right",
plot.title = element_text(face = "bold")
)
Our model building structure focused on increasing complexity with each level. Below is our process:
flowchart LR A(Local Census Data) --> B(Station Features) B --> C(Accessibilty & Built Environment) C --> D(Interaction Terms) D --> E(Neighborhood Clusters)
This model only takes into account surrounding census tract data for each bus stop. After several iterations (involving removing multicolinear and insignificant variables), the three census variables we chose were:
Population density
The log of percent of workers with no vehicle
The log of median household income
#Preparing for 10-fold CV test
ctrl <- trainControl(method = "cv", number = 10)
#Building first simple model: census data
mod1 <- train(log_avg_weekday_boardings ~ pop_density + log_pct_no_vehicle + log_med_hh_inc, data = ridership_full, method = "lm",
trControl = ctrl,
na.action = na.exclude)
summary(mod1)
Call:
lm(formula = .outcome ~ ., data = dat)
Residuals:
Min 1Q Median 3Q Max
-2.3074 -0.7366 -0.1247 0.6151 4.6540
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.0896909 0.2505008 12.33 <0.0000000000000002 ***
pop_density 0.0091368 0.0006647 13.75 <0.0000000000000002 ***
log_pct_no_vehicle 0.2317117 0.0135893 17.05 <0.0000000000000002 ***
log_med_hh_inc -0.2314745 0.0215829 -10.72 <0.0000000000000002 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9824 on 10003 degrees of freedom
Multiple R-squared: 0.1212, Adjusted R-squared: 0.121
F-statistic: 459.9 on 3 and 10003 DF, p-value: < 0.00000000000000022This model solely bases itself on surrounding census tract data for each station. This resulted in an \(R^2\) of 0.121 and a residual standard error of 0.982. The table below describes how each variable relates to the model.
| Variable | Coefficient | P-Score |
|---|---|---|
| Intercept | 3.090 | < .001 |
| Population Density | 0.009 | < .001 |
| Log(Percent of Workers with No Vehicle) | 0.232 | < .001 |
| Log(Median Household Income) | -0.231 | < .001 |
Stations qualities pertain to the qualities the features of individual bus stations. After exploring multicolinearity and levels of significance, the additional variables added are:
The minimum average headway during morning rush hour at the station (in minutes)
The number of routes associated with the station
The minimum average headway at night (in minutes)
mod2 <- train(log_avg_weekday_boardings ~ pop_density + log_pct_no_vehicle + log_med_hh_inc + log_headwayRushAM + headwayNight + numRoutes,
data = ridership_full, method = "lm",
trControl = ctrl,
na.action = na.exclude)
summary(mod2)
Call:
lm(formula = .outcome ~ ., data = dat)
Residuals:
Min 1Q Median 3Q Max
-2.3581 -0.6803 -0.0951 0.5971 4.1872
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.6051460 0.2808941 19.955 < 0.0000000000000002 ***
pop_density 0.0081927 0.0006552 12.504 < 0.0000000000000002 ***
log_pct_no_vehicle 0.1562244 0.0136782 11.421 < 0.0000000000000002 ***
log_med_hh_inc -0.2131326 0.0214679 -9.928 < 0.0000000000000002 ***
log_headwayRushAM -0.6916743 0.0445837 -15.514 < 0.0000000000000002 ***
headwayNight -0.0041473 0.0007276 -5.700 0.000000012320859 ***
numRoutes -0.0475404 0.0064616 -7.357 0.000000000000203 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9509 on 9736 degrees of freedom
Multiple R-squared: 0.1748, Adjusted R-squared: 0.1743
F-statistic: 343.6 on 6 and 9736 DF, p-value: < 0.00000000000000022Adding the station features increased the \(R^2\) to 0.174. Furthermore, it lowered the residual standard error from 0.982 to 0.951 The variables relate to the model as follows:
| Variable | Coefficient | P-Score |
|---|---|---|
| Intercept | 5.605 | < .001 |
| Population Density | 0.008 | < .001 |
| Log(Percent of Workers with No Vehicle) | 0.156 | < .001 |
| Log(Median Household Income) | -0.213 | < .001 |
| Log(Average Headway at Morning Rush Hour) | -0.692 | < .001 |
| Average Headway at Night | -.004 | < .001 |
| Number of Routes | -.047 | < .001 |
Accessibility factors focus on the ability for potential riders to reach the station. After exploring multicolinearity and levels of significance, the additional variables added are:
The distance to the nearest subway/trolley station (in miles)
The number of jobs within 30 minutes of car-less travel (in thousands)
mod3 <- train(log_avg_weekday_boardings ~ pop_density + log_pct_no_vehicle + log_med_hh_inc + log_headwayRushAM + headwayMorn + numRoutes + distRail + I(jobs_30/1000), data = ridership_full, method = "lm",
trControl = ctrl,
na.action = na.exclude)
summary(mod3)
Call:
lm(formula = .outcome ~ ., data = dat)
Residuals:
Min 1Q Median 3Q Max
-2.2374 -0.6749 -0.0969 0.5845 4.0574
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.2396832 0.3142627 19.855 < 0.0000000000000002 ***
pop_density 0.0063788 0.0006845 9.319 < 0.0000000000000002 ***
log_pct_no_vehicle 0.0585610 0.0192104 3.048 0.00231 **
log_med_hh_inc -0.2702899 0.0247857 -10.905 < 0.0000000000000002 ***
log_headwayRushAM -0.5596686 0.0420555 -13.308 < 0.0000000000000002 ***
headwayMorn -0.0049220 0.0004651 -10.583 < 0.0000000000000002 ***
numRoutes -0.0504721 0.0065796 -7.671 0.0000000000000187 ***
distRail -0.0308129 0.0057431 -5.365 0.0000000826940779 ***
`I(jobs_30/1000)` 0.0006411 0.0001085 5.907 0.0000000036004048 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9427 on 9849 degrees of freedom
Multiple R-squared: 0.1895, Adjusted R-squared: 0.1888
F-statistic: 287.8 on 8 and 9849 DF, p-value: < 0.00000000000000022Adding the station features increased the \(R^2\) to 0.189. Furthermore, it lowered the residual standard error to 0.943 The variables relate to the model as follows:
| Variable | Coefficient | P-Score |
|---|---|---|
| Intercept | 6.240 | < .001 |
| Population Density | 0.006 | < .001 |
| Log(Percent of Workers with No Vehicle) | 0.059 | ~ .002 |
| Log(Median Household Income) | -0.270 | < .001 |
| Log(Average Headway at Morning Rush Hour) | -0.560 | < .001 |
| Average Headway at Night | -.005 | < .001 |
| Number of Routes | -.050 | < .001 |
| Distance to Nearest Subway or Trolley Station | -.031 | < .001 |
| Jobs within 30 minutes of Car-less Travel (in thousands) | 6.41 * 10^(-4) | < .001 |
Some variables amplify the impact of others. To better illustrate these connections, we added the following interaction terms:
The number of routes * percentage of workers with no vehicle
The number of jobs * percentage of workers with no vehicle
Proximity to the nearest subway/trolley station * renting percentage of local census tract
With the addition of interaction variables, some the following variables were removed due to multicolinearity:
Number of jobs within 30 minutes of car-less travel
Number of routes using the station
mod4 <- train(log_avg_weekday_boardings ~ log_pct_no_vehicle + log_med_hh_inc +
pop_density + log_headwayRushAM + headwayMorn+ distRail + routes_x_noveh + jobs_x_noveh + railprox_x_renter, data = ridership_full, method = "lm",
trControl = ctrl,
na.action = na.exclude)
summary(mod4)
Call:
lm(formula = .outcome ~ ., data = dat)
Residuals:
Min 1Q Median 3Q Max
-2.7381 -0.6722 -0.0919 0.5989 3.9036
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.699972127525 0.301436128140 18.909 < 0.0000000000000002
log_pct_no_vehicle 0.097381350180 0.020756655298 4.692 0.00000274764777542
log_med_hh_inc -0.230209565812 0.022526581020 -10.219 < 0.0000000000000002
pop_density 0.007016923967 0.000682590901 10.280 < 0.0000000000000002
log_headwayRushAM -0.580846743208 0.041685016708 -13.934 < 0.0000000000000002
headwayMorn -0.004821650365 0.000459408744 -10.495 < 0.0000000000000002
distRail -0.030822239839 0.005832017109 -5.285 0.00000012838998866
routes_x_noveh -0.002512884931 0.000235515557 -10.670 < 0.0000000000000002
jobs_x_noveh 0.000000012073 0.000000003123 3.866 0.000111
railprox_x_renter 0.000789533451 0.000101385386 7.787 0.00000000000000753
(Intercept) ***
log_pct_no_vehicle ***
log_med_hh_inc ***
pop_density ***
log_headwayRushAM ***
headwayMorn ***
distRail ***
routes_x_noveh ***
jobs_x_noveh ***
railprox_x_renter ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9384 on 9848 degrees of freedom
Multiple R-squared: 0.1969, Adjusted R-squared: 0.1961
F-statistic: 268.2 on 9 and 9848 DF, p-value: < 0.00000000000000022Adding the station features increased the \(R^2\) to 0.196. Furthermore, it lowered the residual standard error to 0.938 The variables relate to the model as follows:
| Variable | Coefficient | P-Score |
|---|---|---|
| Intercept | 5.700 | < .001 |
| Population Density | 0.007 | < .001 |
| Log(Percent of Workers with No Vehicle) | 0.097 | < .001 |
| Log(Median Household Income) | -0.230 | < .001 |
| Log(Average Headway at Morning Rush Hour) | -0.581 | < .001 |
| Average Headway at Night | -.005 | < .001 |
| Distance to Nearest Subway or Trolley Station | -.031 | < .001 |
| Number of Routes * Vehicle-less Worker Percentage | -.003 | < .001 |
| Number of Jobs within 30 Minutes * Vehicle-less Worker Percentage | 1.211* 10^(-8) | < .001 |
| Proximity to Subway/Trolley Station * Renting Percentage | 0.001 | < .001 |
Some variables amplify the impact of others. To better illustrate these connections, we added the following interaction terms:
With the addition of interaction variables, some the following variables were removed due to lack of significance:
#All categories above 10, analysis safe
print(ridership_full %>%
count(neighborhood_cluster) %>%
arrange(n)) neighborhood_cluster n
1 <NA> 62
2 2 995
3 5 1031
4 4 1741
5 1 2770
6 3 4299#As.factor used to make neighborhood cluster categorical
mod5 <- train(log_avg_weekday_boardings ~ log_pct_no_vehicle +
log_med_hh_inc + pop_density + log_headwayRushAM + headwayMorn +
routes_x_noveh + jobs_x_noveh + railprox_x_renter +
as.factor(neighborhood_cluster),
data = ridership_full, method = "lm",
trControl = ctrl,
na.action = na.exclude)
summary(mod5)
Call:
lm(formula = .outcome ~ ., data = dat)
Residuals:
Min 1Q Median 3Q Max
-2.8970 -0.6606 -0.0856 0.5803 3.8585
Coefficients:
Estimate Std. Error t value
(Intercept) 4.153216299736 0.377170359929 11.012
log_pct_no_vehicle 0.058230990130 0.019191637877 3.034
log_med_hh_inc -0.111225129388 0.030517274661 -3.645
pop_density 0.003495040539 0.000744823394 4.692
log_headwayRushAM -0.576685668852 0.041189938568 -14.001
headwayMorn -0.004767795038 0.000456476269 -10.445
routes_x_noveh -0.002747983904 0.000239239412 -11.486
jobs_x_noveh 0.000000010025 0.000000003162 3.170
railprox_x_renter 0.000806885658 0.000099408544 8.117
`as.factor(neighborhood_cluster)2` 0.683799974025 0.054408591803 12.568
`as.factor(neighborhood_cluster)3` 0.456374854958 0.033510336885 13.619
`as.factor(neighborhood_cluster)4` 0.605082957306 0.046182722925 13.102
`as.factor(neighborhood_cluster)5` 0.326468738344 0.046026340676 7.093
Pr(>|t|)
(Intercept) < 0.0000000000000002 ***
log_pct_no_vehicle 0.002418 **
log_med_hh_inc 0.000269 ***
pop_density 0.000002736027981430 ***
log_headwayRushAM < 0.0000000000000002 ***
headwayMorn < 0.0000000000000002 ***
routes_x_noveh < 0.0000000000000002 ***
jobs_x_noveh 0.001529 **
railprox_x_renter 0.000000000000000536 ***
`as.factor(neighborhood_cluster)2` < 0.0000000000000002 ***
`as.factor(neighborhood_cluster)3` < 0.0000000000000002 ***
`as.factor(neighborhood_cluster)4` < 0.0000000000000002 ***
`as.factor(neighborhood_cluster)5` 0.000000000001401856 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9283 on 9845 degrees of freedom
Multiple R-squared: 0.2143, Adjusted R-squared: 0.2134
F-statistic: 223.8 on 12 and 9845 DF, p-value: < 0.00000000000000022Adding the station features increased the \(R^2\) to 0.213. Furthermore, it lowered the residual standard error to 0.928. The variables relate to the model as follows:
| Variable | Coefficient | P-Score |
|---|---|---|
| Intercept | 4.153 | < .001 |
| Population Density | 0.003 | < .001 |
| Log(Percent of Workers with No Vehicle) | 0.058 | .002 |
| Log(Median Household Income) | -0.111 | < .001 |
| Log(Average Headway at Morning Rush Hour) | -0.577 | < .001 |
| Average Headway at Night | -.005 | < .001 |
| Number of Routes * Vehicle-less Worker Percentage | -.003 | < .001 |
| Number of Jobs within 30 Minutes * Vehicle-less Worker Percentage | 1.002* 10^(-8) | .002 |
| Proximity to Subway/Trolley Station * Renting Percentage | 0.001 | < .001 |
| Neighborhood Cluster 2 | 0.684 | < .001 |
| Neighborhood Cluster 3 | 0.456 | < .001 |
| Neighborhood Cluster 4 | 0.605 | < .001 |
| Neighborhood Cluster 5 | 0.326 | < .001 |
Note: Generalized VIF^(1/(2DF)) is substituted for VIF in this model due to the categorical nature of the “Neighborhood Clusters” variable
Some of the assumptions do not work with the data format of train(), must make lm() copies to use to test the relationships.
asmpt1 <- lm(log_avg_weekday_boardings ~ pop_density + log_pct_no_vehicle + log_med_hh_inc, data = ridership_full)
asmpt2 <- lm(log_avg_weekday_boardings ~ pop_density + log_pct_no_vehicle + log_med_hh_inc, data = ridership_full)
asmpt3 <- lm(log_avg_weekday_boardings ~ pop_density + log_pct_no_vehicle + log_med_hh_inc + log_headwayRushAM + headwayMorn + numRoutes + distRail + I(jobs_30/1000), data = ridership_full)
asmpt4 <- lm(log_avg_weekday_boardings ~ log_pct_no_vehicle + log_med_hh_inc +
pop_density + log_headwayRushAM + headwayMorn+ distRail + routes_x_noveh + jobs_x_noveh + railprox_x_renter, data = ridership_full)
asmpt5 <- lm(log_avg_weekday_boardings ~ log_pct_no_vehicle +
log_med_hh_inc + pop_density + log_headwayRushAM + headwayMorn +
routes_x_noveh + jobs_x_noveh + railprox_x_renter +
as.factor(neighborhood_cluster),
data = ridership_full)We used a residual plot to check the linearity of our models.
#Linearity test for model 1
diaMod1 <- data.frame(residuals = residuals(mod1),
fitted = fitted(mod1))
ggplot(diaMod1, aes(x = fitted, y = residuals)) +
geom_point() +
geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
labs(title = "Model 1 Residual Plot", x = "Fitted Values", y = "Residuals") +
plotTheme
#for model 2
diaMod2 <- data.frame(residuals = residuals(mod2),
fitted = fitted(mod2))
ggplot(diaMod2, aes(x = fitted, y = residuals)) +
geom_point() +
geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
labs(title = "Model 2 Residual Plot", x = "Fitted Values", y = "Residuals") +
plotTheme
#for model 3
diaMod3 <- data.frame(residuals = residuals(mod3),
fitted = fitted(mod3))
ggplot(diaMod3, aes(x = fitted, y = residuals)) +
geom_point() +
geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
labs(title = "Model 3 Residual Plot", x = "Fitted Values", y = "Residuals") +
plotTheme
#for model 4
diaMod4 <- data.frame(residuals = residuals(mod4),
fitted = fitted(mod4))
ggplot(diaMod4, aes(x = fitted, y = residuals)) +
geom_point() +
geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
labs(title = "Model 4 Residual Plot", x = "Fitted Values", y = "Residuals") +
plotTheme
#for model 5
diaMod5 <- data.frame(residuals = residuals(mod5),
fitted = fitted(mod5))
ggplot(diaMod5, aes(x = fitted, y = residuals)) +
geom_point() +
geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
labs(title = "Model 5 Residual Plot", x = "Fitted Values", y = "Residuals") +
plotTheme
Relationship is not linear in regards to negative residuals across all five OLS models. Such a display violates the assumptions of OLS regression modeling and illustrates that OLS is most likely not the method we should use.
The Breusch-Pagan tests the homoscedasticity of the model, with p-values above 0.05 acceptable. Values below 0.05 demonstrates heteroscedasticity, meaning the variance of the variables is not constant.
bptest(asmpt1)
studentized Breusch-Pagan test
data: asmpt1
BP = 99.558, df = 3, p-value < 0.00000000000000022bptest(asmpt2)
studentized Breusch-Pagan test
data: asmpt2
BP = 99.558, df = 3, p-value < 0.00000000000000022bptest(asmpt3)
studentized Breusch-Pagan test
data: asmpt3
BP = 395.94, df = 8, p-value < 0.00000000000000022bptest(asmpt4)
studentized Breusch-Pagan test
data: asmpt4
BP = 546.28, df = 9, p-value < 0.00000000000000022bptest(asmpt5)
studentized Breusch-Pagan test
data: asmpt5
BP = 565.99, df = 12, p-value < 0.00000000000000022OLS modeling in this context fails too, meaining our prediction intervals based upon these models are not the most trustworthy.
VIF testing explores the level of correlation between variables, with higher values relating to more correlation. A VIF value over 5 is concerning and over 10 is dangerous.
vif(asmpt1) pop_density log_pct_no_vehicle log_med_hh_inc
1.319542 1.573044 1.237306 vif(asmpt2) pop_density log_pct_no_vehicle log_med_hh_inc
1.319542 1.573044 1.237306 vif(asmpt3) pop_density log_pct_no_vehicle log_med_hh_inc log_headwayRushAM
1.481148 3.340368 1.743951 1.768063
headwayMorn numRoutes distRail I(jobs_30/1000)
1.699740 1.140079 2.330049 2.532604 vif(asmpt4)log_pct_no_vehicle log_med_hh_inc pop_density log_headwayRushAM
3.935158 1.453617 1.486348 1.752829
headwayMorn distRail routes_x_noveh jobs_x_noveh
1.673428 2.424581 2.060528 3.654173
railprox_x_renter
2.153705 #Note that with the categorical variable in model 5 means that GVIF^(1/(2*DF)) is the stand in for VIF
vif(asmpt5) GVIF Df GVIF^(1/(2*Df))
log_pct_no_vehicle 3.437801 1 1.854131
log_med_hh_inc 2.726215 1 1.651125
pop_density 1.808486 1 1.344800
log_headwayRushAM 1.748924 1 1.322469
headwayMorn 1.688318 1 1.299353
routes_x_noveh 2.172771 1 1.474032
jobs_x_noveh 3.829587 1 1.956933
railprox_x_renter 2.115886 1 1.454608
as.factor(neighborhood_cluster) 8.030983 4 1.297466No variables are above 5 in any model, meaning that they are not correlated enough to diminish results.
Cook’s distance analyses demonstrate how much influence each data point has over the OLS model. It does so by measuring how much fitted values would change without each data point. Typically, a point with a cook’s distance over 4/n (n being the number of data points, between 10,885 and 10,894 depending on the model).
Looking at the results of the cook’s distance analysis, the OLS model has an high number of influential points. (constantly above 500). However, in the context of the scale of the data set, the share
Comparing the results of our OLS Models, we looked at both RMSE, MAE, and R squared to see how each iteration of the model performed. Furthermore, we included the exponentiated RMSE and MAE (expRMSE and expMAE) in order to translate the units of measurement to from log(number of boardings) to number of boardings.
?exp
results <- data.frame(
Model = c("Census Data", "+ Station Features",
"+ Accessibility",
"+ Interaction Terms",
"+ Neighborhood Clustering"),
RMSE = c(mod1$results$RMSE,
mod2$results$RMSE,
mod3$results$RMSE,
mod4$results$RMSE,
mod5$results$RMSE),
MAE = c(mod1$results$MAE,
mod2$results$MAE,
mod3$results$MAE,
mod4$results$MAE,
mod5$results$MAE),
expRMSE = c(exp(mod1$results$RMSE),
exp(mod2$results$RMSE),
exp(mod3$results$RMSE),
exp(mod4$results$RMSE),
exp(mod5$results$RMSE)),
expMAE = c(exp(mod1$results$MAE),
exp(mod2$results$MAE),
exp(mod3$results$MAE),
exp(mod4$results$MAE),
exp(mod5$results$MAE)),
rSquared = c(mod1$results$Rsquared,
mod2$results$Rsquared,
mod3$results$Rsquared,
mod4$results$Rsquared,
mod5$results$Rsquared)
)
kable(results)| Model | RMSE | MAE | expRMSE | expMAE | rSquared |
|---|---|---|---|---|---|
| Census Data | 0.9823305 | 0.7898383 | 2.670673 | 2.203040 | 0.1214881 |
| + Station Features | 0.9513893 | 0.7589638 | 2.589304 | 2.136062 | 0.1738714 |
| + Accessibility | 0.9431210 | 0.7513378 | 2.567983 | 2.119834 | 0.1883762 |
| + Interaction Terms | 0.9387835 | 0.7502112 | 2.556869 | 2.117447 | 0.1960932 |
| + Neighborhood Clustering | 0.9290225 | 0.7395651 | 2.532033 | 2.095024 | 0.2125934 |
First, in using count models such as Poisson and Negative Binomial, we run into an issue: our dependent variable is an average so some values show up as decimals. In order to account for this, we can round all the averages. While there are some trade offs with this approach, including introducing more zeros to an already zero-inflated distribution and adding a small amount of error in estimates, we believe this preserves some interpretability and comparability with our OLS models.
Below, we first do this rounding before fitting a Poisson Regression model.
In terms of picking predictors, we went with the predictors from OLS Model 5 (all including neighborhood effects). When running Poisson, pop_density was found to be insignificant and removed.
ridership_full <- ridership_full %>%
mutate(count_weekday_boardings = as.integer(round(avg_weekday_boardings)))
# Fit Poisson regression
model_poisson <- glm(
count_weekday_boardings ~ log_pct_no_vehicle + log_med_hh_inc + log_headwayRushAM + headwayMorn +
routes_x_noveh + jobs_x_noveh + railprox_x_renter + as.factor(neighborhood_cluster),
data = ridership_full,
family = poisson(link = "log")
)
# View results
summary(model_poisson)
Call:
glm(formula = count_weekday_boardings ~ log_pct_no_vehicle +
log_med_hh_inc + log_headwayRushAM + headwayMorn + routes_x_noveh +
jobs_x_noveh + railprox_x_renter + as.factor(neighborhood_cluster),
family = poisson(link = "log"), data = ridership_full)
Coefficients:
Estimate Std. Error z value
(Intercept) 2.424734680451 0.149531344565 16.216
log_pct_no_vehicle 0.192649704189 0.009031895261 21.330
log_med_hh_inc 0.065825654515 0.012124111405 5.429
log_headwayRushAM -0.543844395852 0.016888860918 -32.201
headwayMorn -0.011879054351 0.000195224188 -60.848
routes_x_noveh -0.003048092804 0.000087454990 -34.853
jobs_x_noveh -0.000000025610 0.000000001082 -23.673
railprox_x_renter 0.002287005556 0.000026666858 85.762
as.factor(neighborhood_cluster)2 0.958201129355 0.019975196918 47.970
as.factor(neighborhood_cluster)3 0.606919925170 0.016863811521 35.989
as.factor(neighborhood_cluster)4 0.993237736697 0.019676060747 50.480
as.factor(neighborhood_cluster)5 0.517983157920 0.020641761112 25.094
Pr(>|z|)
(Intercept) < 0.0000000000000002 ***
log_pct_no_vehicle < 0.0000000000000002 ***
log_med_hh_inc 0.0000000566 ***
log_headwayRushAM < 0.0000000000000002 ***
headwayMorn < 0.0000000000000002 ***
routes_x_noveh < 0.0000000000000002 ***
jobs_x_noveh < 0.0000000000000002 ***
railprox_x_renter < 0.0000000000000002 ***
as.factor(neighborhood_cluster)2 < 0.0000000000000002 ***
as.factor(neighborhood_cluster)3 < 0.0000000000000002 ***
as.factor(neighborhood_cluster)4 < 0.0000000000000002 ***
as.factor(neighborhood_cluster)5 < 0.0000000000000002 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 138103 on 9857 degrees of freedom
Residual deviance: 106235 on 9846 degrees of freedom
(1040 observations deleted due to missingness)
AIC: 133126
Number of Fisher Scoring iterations: 6A key assumption of the Poisson model is that the variance equals the mean. When the variance is substantially larger than the mean, the data are considered overdispersed. This condition biases the model’s standard error estimates. We assess the level of overdispersion by calculating the dispersion parameter from the Pearson residuals.
# Fit Poisson model
# Calculate dispersion parameter
dispersion <- sum(residuals(model_poisson, type = "pearson")^2) /
model_poisson$df.residual
cat("Dispersion parameter:", round(dispersion, 3), "\n")Dispersion parameter: 19.818 # Rule of thumb:
# < 1.5: Poisson OK
# 1.5 - 3: Mild overdispersion, NegBin recommended
# > 3: Serious overdispersion, NegBin essential
# Example output:
# Dispersion parameter: 4.523
# Warning: Serious overdispersion detected!A dispersion parameter of 19.818 is incredibly high. This shows severe overdispersion; the actual variance in ridership is over 19 times the mean assumed by the Poisson model. Because this makes the Poisson model’s standard errors inaccurate, it is not adequate. We will use the Negative Binomial model next, as it handles this excess variability. Both pop_density and log_hh_med_inc were removed as predictors for being insignificant and not having an effect on final results.
model_nb <- glm.nb(
count_weekday_boardings ~ log_pct_no_vehicle + log_headwayRushAM + headwayMorn +
routes_x_noveh + jobs_x_noveh + railprox_x_renter + as.factor(neighborhood_cluster),
data = ridership_full
)
summary(model_nb)
Call:
glm.nb(formula = count_weekday_boardings ~ log_pct_no_vehicle +
log_headwayRushAM + headwayMorn + routes_x_noveh + jobs_x_noveh +
railprox_x_renter + as.factor(neighborhood_cluster), data = ridership_full,
init.theta = 0.7102090634, link = log)
Coefficients:
Estimate Std. Error z value
(Intercept) 3.707101080332 0.171660420357 21.596
log_pct_no_vehicle 0.164581210930 0.024534794402 6.708
log_headwayRushAM -0.768096341631 0.056174516427 -13.673
headwayMorn -0.008893493231 0.000622807427 -14.280
routes_x_noveh -0.002257856648 0.000317404748 -7.113
jobs_x_noveh -0.000000019325 0.000000004164 -4.641
railprox_x_renter 0.002095006857 0.000130188271 16.092
as.factor(neighborhood_cluster)2 0.927999259613 0.063760596532 14.554
as.factor(neighborhood_cluster)3 0.645615110481 0.042302176705 15.262
as.factor(neighborhood_cluster)4 0.906239662663 0.055682302752 16.275
as.factor(neighborhood_cluster)5 0.521331984988 0.056133518395 9.287
Pr(>|z|)
(Intercept) < 0.0000000000000002 ***
log_pct_no_vehicle 0.00000000001972 ***
log_headwayRushAM < 0.0000000000000002 ***
headwayMorn < 0.0000000000000002 ***
routes_x_noveh 0.00000000000113 ***
jobs_x_noveh 0.00000345894885 ***
railprox_x_renter < 0.0000000000000002 ***
as.factor(neighborhood_cluster)2 < 0.0000000000000002 ***
as.factor(neighborhood_cluster)3 < 0.0000000000000002 ***
as.factor(neighborhood_cluster)4 < 0.0000000000000002 ***
as.factor(neighborhood_cluster)5 < 0.0000000000000002 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(0.7102) family taken to be 1)
Null deviance: 14007 on 10001 degrees of freedom
Residual deviance: 11085 on 9991 degrees of freedom
(896 observations deleted due to missingness)
AIC: 56239
Number of Fisher Scoring iterations: 1
Theta: 0.7102
Std. Err.: 0.0116
2 x log-likelihood: -56214.6530 cat("\n")Below, we conduct Spatial Cross-Validation (CV), which is necessary because standard random CV is unreliable for spatial data. Randomly splitting stops mixes geographically correlated observations across training and testing sets, leading to spatial leakage. This leakage allows the model to “peek” at adjacent information, resulting in overly optimistic accuracy estimates. By using the LOGO method based on our neighborhood_cluster groups, we ensure the Negative Binomial model is tested on unobserved, independent spatial regions, providing a better measure of how well the model generalizes across our diverse clusters.
# Get unique clusters
clusters <- ridership_full %>%
filter(!is.na(neighborhood_cluster)) %>%
distinct(neighborhood_cluster) %>%
pull()
cat("Running LOGO Cross-Validation by Neighborhood Cluster...\n")Running LOGO Cross-Validation by Neighborhood Cluster...cat("Total Clusters (Folds):", length(clusters), "\n\n")Total Clusters (Folds): 5 # Initialize results storage
cv_results_cluster <- tibble() # Store cluster-level metrics
cv_predictions_cluster <- tibble() # Store all stop-level predictions (here, you will see how large errors can be at stop level. Can inform future work)
# Loop through each cluster
for (cluster in clusters) {
# Print progress
cat("Processing cluster:", cluster, "...\n")
# Split data: Leave-One-Group-Out (LOGO)
train_data <- ridership_full %>%
filter(neighborhood_cluster != cluster,
!is.na(count_weekday_boardings))
test_data <- ridership_full %>%
filter(neighborhood_cluster == cluster,
!is.na(count_weekday_boardings))
# --- Fit Model ---
model_cv <- glm.nb(
count_weekday_boardings ~ log_pct_no_vehicle + log_med_hh_inc +
pop_density + numRoutes + log_headwayRushAM +
routes_x_noveh + jobs_x_noveh + railprox_x_renter,
data = train_data,
na.action = na.exclude
)
# --- Predict and Calculate Metrics ---
# Predict
test_data <- test_data %>%
mutate(
prediction = suppressWarnings(
predict(model_cv, newdata = ., type = "response")
),
error = count_weekday_boardings - prediction # Actual - Predicted
)
# Calculate derived metrics using the error
test_data <- test_data %>%
mutate(
pct_error_stop = (error / (count_weekday_boardings + 1)) * 100 # +1 to avoid div by zero (okay if we rounded down to 0...we normalize to % anyways)
)
# Store stop-level predictions
cv_predictions_cluster <- bind_rows(
cv_predictions_cluster,
test_data %>%
select(stop_id, neigh_name, neighborhood_cluster,
count_weekday_boardings, prediction, error, pct_error_stop)
)
# Calculate cluster-level metrics
mae <- mean(abs(test_data$error), na.rm = TRUE)
rmse <- sqrt(mean(test_data$error^2, na.rm = TRUE))
mean_actual <- mean(test_data$count_weekday_boardings, na.rm = TRUE)
mean_predicted <- mean(test_data$prediction, na.rm = TRUE)
mean_error <- mean(test_data$error, na.rm = TRUE) # Signed error (Actual - Predicted)
# Store cluster summary
cv_results_cluster <- bind_rows(
cv_results_cluster,
tibble(
fold = cluster,
test_cluster = cluster,
n_test = nrow(test_data),
mae = mae,
rmse = rmse,
mean_actual = mean_actual,
mean_predicted = mean_predicted,
mean_error = mean_error,
pct_error_mae = (mae / mean_actual) * 100,
signed_pct_error = (mean_error / mean_actual) * 100
)
)
}Processing cluster: 3 ...
Processing cluster: 4 ...
Processing cluster: 1 ...
Processing cluster: 2 ...
Processing cluster: 5 ...cat("\nCross-Validation Complete\n")
Cross-Validation Completecat("Successful folds:", sum(!is.na(cv_results_cluster$mae)), "/", length(clusters), "\n")Successful folds: 5 / 5 cat("--- Overall Mean CV Performance ---\n")--- Overall Mean CV Performance ---cat("Mean MAE:", round(mean(cv_results_cluster$mae, na.rm = TRUE), 2), "boardings\n")Mean MAE: 7.53 boardingscat("Mean RMSE:", round(mean(cv_results_cluster$rmse, na.rm = TRUE), 2), "boardings\n")Mean RMSE: 14.92 boardingscat("Mean Signed Error:", round(mean(cv_results_cluster$mean_error, na.rm = TRUE), 2), "boardings\n")Mean Signed Error: 0.18 boardingsBecause ridership varies greatly across the city, raw metrics like MAE and RMSE are poor indicators for our goal of finding under/over-performing areas. Absolute errors are difficult to interpret across different scales. Thus, we focus on the cluster-level percentage errors (MAE % Error and Signed % Error) below. These normalized results offer a more accurate and directly comparable overview of errors and success across diverse neighborhood clusters.
neg_bin_cv_kable <- cv_results_cluster %>%
# Select and rename columns for clarity in the final report
select(
`Test Cluster` = test_cluster,
`Stops (N)` = n_test,
`MAE (Boardings)` = mae,
`RMSE (Boardings)` = rmse,
`Mean Error` = mean_error,
`MAE % Error` = pct_error_mae,
`Signed % Error` = signed_pct_error
) %>%
# Order by cluster for consistent presentation
arrange(`Test Cluster`) %>%
# Create the kable object
kable(
caption = "Negative Binomial Model Performance: Leave-One-Cluster-Out Cross-Validation",
digits = 2,
align = 'c'
) %>%
# Apply styling
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = F
) %>%
# Format the header for better readability
add_header_above(c(" " = 2, "Absolute Error" = 2, "Signed Error" = 1, "Percentage Error" = 2))
# Saving for presentation
save_kable(
neg_bin_cv_kable,
file = "presentation-assets/negative_binomial_cv_results.png",
zoom = 2
)
neg_bin_cv_kable| Test Cluster | Stops (N) | MAE (Boardings) | RMSE (Boardings) | Mean Error | MAE % Error | Signed % Error |
|---|---|---|---|---|---|---|
| 1 | 2770 | 4.22 | 6.48 | -1.93 | 168.42 | -76.90 |
| 2 | 995 | 8.64 | 15.65 | 2.06 | 83.65 | 20.00 |
| 3 | 4299 | 6.17 | 11.58 | -0.08 | 96.05 | -1.31 |
| 4 | 1741 | 12.00 | 30.27 | 2.01 | 92.43 | 15.52 |
| 5 | 1031 | 6.60 | 10.64 | -1.17 | 106.33 | -18.82 |
These results are discussed in the following section.
Because of the broken assumptions of OLS, as discussed in Section 5.1, and the overdispersion parameter after Poisson Regression, as discussed in Section 5.2 and 5.3, we present the Negative Binomial Model as our final recommendation.
The Spatial Cross-Validation by cluster reveals significant differences in model performance and predictive accuracy across the distinct neighborhood groups. The overall prediction error, measured by the MAE % Error, varies widely, from 68.42% to 106.33% across the five clusters. The raw MAE (measured in boardings) is an absolute error metric that is poorly suited for comparing performance across clusters because it fails to account for the magnitude of a cluster’s average ridership, thus hiding the true relative predictive success.
Thus, more critically, the Signed % Error identifies prediction tendencies essential for understanding service by cluster type. Cluster 2 shows the highest positive residual (+20.00% signed error), indicating the model consistently under-predicted ridership in this group; this highlights stops likely experiencing unmet demand. Conversely, Cluster 1 shows the most negative residual (−76.90% signed error), meaning the model heavily over-predicted ridership in that cluster, suggesting those stops, despite having reasons for high ridership, are facing some sort of barriers that our model does not capture.
The cluster-level percentage error informs high-level policy, but a granular understanding requires mapping the stop-level residuals. Since hundereds of individual, overlapping stop points make direct interpretation difficult, we employ a hexbin map. This technique effectively aggregates nearby stop residuals into visible bins, providing a far more interpretable view of the spatial % error intensity.
Data prep and plotting of this data is done below.
stop_geometry <- ridership_philly %>%
st_transform(2272) %>%
select(stop_id, geometry) %>%
distinct(stop_id, .keep_all = TRUE)
stop_error_map_data <- stop_geometry %>%
# Join the non-spatial error data
left_join(
cv_predictions_cluster %>%
select(stop_id, pct_error_stop, error),
by = "stop_id"
) %>%
filter(!is.na(error))map_data_for_hex <- stop_error_map_data %>%
st_coordinates() %>%
as.data.frame() %>%
mutate(pct_error_stop = stop_error_map_data$pct_error_stop)
p_hexbin_error <- ggplot(map_data_for_hex, aes(x = X, y = Y, z = pct_error_stop)) +
stat_summary_hex(
fun = mean,
bins = 75
) +
scale_fill_gradient2(
low = "maroon", # Red: Negative % Error (Over-predicted)
mid = "gray95", # Near Zero
high = "steelblue", # Blue: Positive % Error (Under-predicted / Unmet Demand)
midpoint = 0,
limits = c(-50, 50),
oob = scales::squish,
name = "Avg. Residual Error\n(Actual - Predicted, as %)"
) +
labs(
title = "Hexbin Map of Under- vs. Over-Performance (Avg. Percentage Error)",
subtitle = "Positive Residuals (Blue Hexagons) = Under Prediction (i.e there is more ridership than predicted)\n
Negative Residuals (Red Hexagons) = Over Prediction (i.e there is less ridership than predicted)",
x = NULL, y = NULL
) +
theme_void() +
theme(
plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
legend.position = "right"
)
ggsave(
filename = "presentation-assets/p_hexbin_error.png",
plot = p_hexbin_error,
width = 10,
height = 8,
dpi = 300
)
p_hexbin_error
In the hexbin map in the previous section, the clearest and most critical finding is the pervasive presence of negative residuals (over-prediction) across numerous geographic areas. This systematic pattern indicates that our model, even after accounting for known demand factors—such as neighborhood effects and transit dependence (e.g., workers without vehicle access) — still expects higher ridership than is actually observed.
This consistent trend in observed ridership (where Actual < Predicted) should be interpreted by policymakers as reason to investigate and eliminate uncaptured barriers. We interpret these results as potential of high ridership, yet ridership is being suppressed by unseen factors. Potential barriers in these hexbin areas warranting investigation include:
The hexbin map also identifies clusters (in North and West Philadelphia in particular) defined by positive residuals (under-prediction), indicating areas where observed ridership exceeds the model’s prediction (Actual > Predicted). This finding reveals unmet transit demand across the network. This trend means that the existing service could be operating under-capacity. This trend of unmet demand should be interpreted by policymakers as reason to increase service support. Potential interventions include:
An important note is our model’s solely examines racial and class aspects of bus ridership implicitly through the percent of workers with no vehicle and neighborhood clustering variables. Furthermore, bus ridership data inherently has racial and class components baked in, as riding the bus is demarcated as “lower-class” or “gross” by those who can afford to avoid taking it. Class is highly tied to race within Philadelphia, as Black Philadelphians are more likely to be lower-income and, thus, rely on the bus for access to the city. These residents are more likely to face barriers to the bus network, such as the cost of a ticket and — due to the lower-income communities being pushed from economic centers like Market Street — longer travel times.
What We Couldn’t Include (due to time/data constraints):
Data Quality Issues: - Ridership from 2024 may not reflect current patterns - Census data lags real-time demographic changes - GTFS shows scheduled service, not actual service delivered
Model Assumptions:
A key assumption in our chosen mode — negative binomial — is that the variance is greater than the mean.
Limitations in our Process:
We were very limited by the short turn around of this report. There are impactful variables and processes we excluded, did not engineer, or did not find due to the temporal limitations. Furthermore, we are limited to techniques that we are knowledgeble about. It is possible there is a better prediction technique in existence that were are unaware of and, thus, did not use.
If we had more time/data:
Concluding the analysis, the Negative Binomial Model was selected as the final model, using key variables such as transit dependence (percent of workers with no access to vehicles), morning rush hour headways, and neighborhood clusters. The model’s signed residuals were then utilized to diagnose where stops are under- or over-performing relative to expected ridership.
We found that the model systematically overpredicts across most of the Philadelphia bus network. This finding, particularly visible in the hexbin analysis, indicates a widespread failure to capture localized ridership barriers. This suggests that SEPTA prioritize investigation to identify and remove any unseen barriers to access.
Conversely, the model significantly underpredicted ridership in specific corridors, notably parts of North and West Philadelphia. This suggests the existence of unmet transit demand in these areas. For these locations, SEPTA should immediately pursue service improvements and capital investment to better support ridership and capture demand.
For SEPTA: The recent SEPTA funding crisis is a stark reminder of the limited budget and resources the agency is working with. Thus, SEPTA could use this model to understand where communities are most reliant on bus services and, therefore, where they may want to allocate more resources towards. SEPTA should also aim to avoid a “transit death spiral”, where reduced service leads to reduced ridership which then leads back to service reductions. Rather, SEPTA should use areas where ridership is lower than conventionally expected as reason to investigate barriers and increase customer satisfaction.
For Riders: Our model’s over-prediction of boardings at most stations across the city illustrates a clear disconnect between potential riders and SEPTA. Riders should organize around the flaws and barriers of the system in order to make it more inclusive. Furthermore, communities where our model under-predicted should advocate for increased service, as their stations are more likely to be under-served by current service.
Acknowledgments: This analysis uses data from SEPTA, the US Census Bureau, UMN Data Repository, and OpenDataPhilly contributors.