Instead of a single “grand model” that averages out the nuances of New York City, we present three focused analyses. Each “story” targets a specific question about how weather, events, and station identity shape the daily rhythm of the subway.

Baseline Rhythms

Before diving into anomalies, we must understand the normal heartbeat of the city. Ridership follows a distinct weekly pattern, with sharp peaks on weekdays and a relaxed, single-hump profile on weekends.

# Ridership over time
p_riders_over_time <- station_day_aug %>%
  distinct(date, day_ridership_system, dow, is_holiday) %>%
  filter(is_holiday == "No") %>%
  ggplot(aes(x = date, y = day_ridership_system / 1e6, color = dow)) +
  geom_point(size = 0.8, alpha = 0.6) +
  geom_smooth(method = "loess", span = 0.1, linewidth = 1, se = FALSE) +
  scale_y_continuous(labels = comma) +
  labs(x = NULL, y = "System Ridership (Millions)", color = NULL) +
  theme(legend.position = "none")

# Boxplot
p_dow_riders <- station_day_aug %>%
  distinct(date, day_ridership_system, dow, is_holiday) %>%
  filter(is_holiday == "No") %>%
  ggplot(aes(x = dow, y = day_ridership_system / 1e6, fill = dow)) +
  geom_boxplot(alpha = 0.8, outlier.size = 0.5) +
  labs(x = NULL, y = NULL) +
  theme(axis.text.y = element_blank(), legend.position = "none")

(p_riders_over_time + p_dow_riders) +
  plot_layout(widths = c(3, 1)) +
  plot_annotation(title = "System Ridership: The Weekly Pulse")

Story 1: The “Goldilocks” Zone

Question: Is there a “perfect” temperature for subway ridership?

New Yorkers are famously resilient, but even they have limits. We hypothesize that ridership peaks in a comfortable “Goldilocks zone” of mild temperatures—not too hot, not too cold. First, let’s examine the raw data to see weather effects before formal modeling.

# Temperature
p_temp <- station_day_aug %>%
  filter(is_holiday == "No") %>%
  distinct(date, tmax, dow, day_ridership_system) %>%
  ggplot(aes(x = tmax, y = day_ridership_system / 1e6)) +
  geom_point(alpha = 0.3, size = 1) +
  geom_smooth(method = "loess", span = 1, color = okabe_ito$orange, se = FALSE) +
  facet_grid(. ~ dow) +
  labs(title = "Ridership vs Max Temp", x = "Temp (°F)", y = "Ridership (M)")

# Precip
p_prcp <- station_day_aug %>%
  filter(is_holiday == "No", prcp > 0, snow == 0) %>%
  distinct(date, prcp, pct_change_system) %>%
  ggplot(aes(x = prcp, y = pct_change_system)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", color = okabe_ito$blue, se = FALSE) +
  scale_x_log10() +
  labs(title = "Rain Impact", x = "Rain (in, log)", y = "Log % Change")

# Snow
p_snow <- station_day_aug %>%
  filter(is_holiday == "No", snow > 0) %>%
  distinct(date, snow, pct_change_system) %>%
  ggplot(aes(x = snow, y = pct_change_system)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", color = okabe_ito$green, se = FALSE) +
  scale_x_log10() +
  labs(title = "Snow Impact", x = "Snow (in, log)", y = NULL)

p_temp / (p_prcp + p_snow)

We often assume a linear relationship between weather and behavior, but human comfort is non-linear. We can explore this with a Generalized Additive Model (GAM), which allows for non-linear relationships and controls for other factors.

The plot below shows the predicted ridership based on temperature, holding other factors constant.

\[ \begin{aligned} \text{Ridership} \sim & \; s(\text{Temperature}) \\ & + s(\text{Precipitation}) \\ & + s(\text{Snowfall}) \\ & + \text{Day of Week} + \text{Month} \end{aligned} \]

# Aggregating to System Level
system_daily <- station_day %>%
  group_by(date) %>%
  summarize(
    total_ridership = sum(entries_day, na.rm = TRUE),
    tmax = first(tmax),
    prcp = first(prcp),
    snow = first(snow),
    dow = first(dow),
    month = first(month),
    is_holiday = first(is_holiday),
    .groups = "drop"
  ) %>%
  filter(total_ridership > 100000) # Remove potential bad data days

# Calculate average ridership for each day of week
dow_avg <- system_daily %>%
  filter(is_holiday == "No") %>%
  group_by(dow) %>%
  summarize(avg_ridership = mean(total_ridership), .groups = "drop")

# Normalize ridership by day of week average
system_daily <- system_daily %>%
  left_join(dow_avg, by = "dow") %>%
  mutate(normalized_ridership = total_ridership / avg_ridership)

# Fit GAM on normalized ridership with temperature and precipitation interaction
gam_model <- gam(normalized_ridership ~ te(tmax, prcp, k = c(8, 5)),
  data = system_daily %>% filter(is_holiday == "No")
)

# Create prediction grid
temp_range <- seq(20, 95, length.out = 50)
prcp_range <- seq(0, 2, length.out = 50)

pred_grid <- expand.grid(
  tmax = temp_range,
  prcp = prcp_range
)

# Predict normalized ridership
pred_grid$predicted <- predict(gam_model, newdata = pred_grid)

# Calculate percent change from average (1.0 for normalized ridership)
# Average normalized ridership = 1.0 by definition
pred_grid$pct_change <- (pred_grid$predicted - 1) / 1 * 100

# Find the optimal temperature for subtitle
dry_preds <- pred_grid %>% filter(prcp == 0)
optimal_idx <- which.max(dry_preds$predicted)
optimal_temp <- dry_preds$tmax[optimal_idx]

# Heatmap showing combined effects
ggplot(pred_grid, aes(x = tmax, y = prcp, fill = pct_change)) +
  geom_tile() +
  scale_fill_gradientn(
    colors = c("#08306b", "#2171b5", "#6baed6", "#c6dbef", "white", "#fcbba1", "#fb6a4a", "#cb181d"),
    values = scales::rescale(c(-35, -20, -12, -5, 0, 2, 5, 10)),
    name = "% Change\nfrom Average",
    limits = c(-35, 10),
    breaks = c(10, 0, -10, -20, -30)
  ) +
  scale_x_continuous(breaks = seq(20, 90, by = 10)) +
  scale_y_continuous(breaks = seq(0, 2, by = 0.5)) +
  geom_contour(aes(z = pct_change), color = "gray40", alpha = 0.7, bins = 10) +
  labs(
    title = "The Weather-Ridership Surface: Temperature × Precipitation",
    subtitle = paste0("% change from average daily ridership (optimal at ", round(optimal_temp), "°F, dry day)"),
    x = "Max Daily Temperature (°F)",
    y = "Precipitation (inches)"
  ) +
  theme(
    panel.grid = element_blank(),
    legend.position = "right"
  )

Finding: This heatmap reveals how temperature and precipitation jointly affect ridership relative to the average. The red “sweet spot” sits between 60-75°F with little to no rain—confirming the Goldilocks hypothesis. Key patterns:

  • Optimal Zone (red, bottom center): Ridership peaks 4-6% above average on mild, dry days around 67°F. This is when walking to stations is most comfortable.
  • Cold + Wet (upper left, deep blue): The worst conditions, with ridership dropping 20-35% below average as commuters avoid frigid, rainy walks to stations.
  • Hot days (bottom right): Even dry days above 85°F see ridership 3-5% below average as extreme heat discourages travel.
  • Rain effect (moving up): Each inch of precipitation causes roughly 5-10% ridership decline, regardless of temperature.

The contour lines show that precipitation has a stronger marginal effect than temperature—a half-inch of rain causes more ridership loss than a 20°F temperature swing within the moderate range.

Story 2: The “Blast Radius” of Mega-Events

Question: Do events like the NYC Marathon affect the whole city, or just the route?

Major events transform New York’s geography for a day, but how far does the impact extend? We analyzed two iconic 2024 events: the NYC Marathon (Nov 3) and the Pride March (June 30). For each station, we compared event-day ridership against a baseline of typical Sundays (median of all non-holiday Sundays in surrounding months), then mapped the percent change geographically.

library(sf)

# Load station locations
stations <- station_day %>%
  distinct(station_complex_id, station_complex) %>%
  left_join(read_csv("data/stations/mta_station_complexes.csv", show_col_types = FALSE) %>%
    select(station_complex_id, latitude, longitude),
  by = "station_complex_id"
  )

# Define events
events_list <- list(
  "NYC Marathon" = as.Date("2024-11-03"),
  "Pride March" = as.Date("2024-06-30")
)

# Function to calculate impact
calc_event_impact <- function(evt_date, name) {
  # Use broader baseline: same day of week across multiple months for robustness
  evt_month <- month(evt_date)
  nearby_months <- c(evt_month - 1, evt_month, evt_month + 1)
  nearby_months <- ((nearby_months - 1) %% 12) + 1  # Handle year wrap

  baseline_dates <- station_day %>%
    filter(
      year(date) == year(evt_date),
      month(date) %in% nearby_months,
      dow == wday(evt_date, label = TRUE),
      is_holiday == "No",
      date != evt_date
    ) %>%
    pull(date) %>%
    unique()

  # Get baseline for each station (median across typical Sundays)
  baseline <- station_day %>%
    filter(date %in% baseline_dates) %>%
    group_by(station_complex_id) %>%
    summarize(baseline_entries = median(entries_day, na.rm = TRUE), .groups = "drop")

  # Get event day ridership
  event_day_data <- station_day %>%
    filter(date == evt_date) %>%
    select(station_complex_id, event_entries = entries_day)

  # Merge and calc pct change
  baseline %>%
    left_join(event_day_data, by = "station_complex_id") %>%
    mutate(
      pct_change = (event_entries - baseline_entries) / baseline_entries,
      # Cap for visualization
      pct_change_cap = pmin(pmax(pct_change, -0.5), 1.0),
      event_name = name
    )
}
# Generate data for both events
marathon_impact <- calc_event_impact(events_list[["NYC Marathon"]], "NYC Marathon (Nov 3)")
pride_impact <- calc_event_impact(events_list[["Pride March"]], "Pride March (June 30)")

plot_data <- bind_rows(marathon_impact, pride_impact) %>%
  left_join(stations, by = "station_complex_id") %>%
  filter(!is.na(latitude))

# Load NYC neighborhood tabulation areas for detailed boundaries
# Data preloaded by data/download_nyc_geo.R
nyc_nta <- st_read("data/geo/nyc_nta.geojson", quiet = TRUE)

# Filter to Manhattan neighborhoods
if ("borough" %in% names(nyc_nta)) {
  manhattan_sf <- nyc_nta %>% filter(borough == "Manhattan")
} else if ("boro_name" %in% names(nyc_nta)) {
  manhattan_sf <- nyc_nta %>% filter(boro_name == "Manhattan")
} else {
  manhattan_sf <- nyc_nta
}

# Plot with detailed Manhattan neighborhoods
ggplot() +
  # Add Manhattan neighborhoods as background with subtle borders
  geom_sf(data = manhattan_sf, fill = "gray92", color = "gray75", linewidth = 0.2) +
  # Add station points
  geom_point(
    data = plot_data,
    aes(x = longitude, y = latitude, color = pct_change_cap),
    alpha = 0.85, size = 3
  ) +
  scale_color_gradient2(
    low = okabe_ito$blue, mid = "white", high = okabe_ito$vermillion, midpoint = 0,
    limits = c(-0.5, 1.0),
    labels = percent_format(), name = "Ridership Change\nvs. Baseline"
  ) +
  facet_wrap(~event_name) +
  coord_sf(
    crs = 4326,
    xlim = c(-74.02, -73.91),
    ylim = c(40.70, 40.88)
  ) +
  labs(
    title = "The Blast Radius: Event Impact on Station Ridership",
    subtitle = "Red = increased ridership (+100%), White = normal, Blue = decreased (-50%)",
    x = NULL, y = NULL
  ) +
  theme_minimal() +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank(),
    legend.position = "bottom",
    strip.text = element_text(face = "bold", size = 12)
  )

Finding: The “Blast Radius” is real and geographically distinct.

  • Marathon: Ridership surges at stations near the finish line (Columbus Circle, 59 St) and along the route through the Upper East Side and Harlem. Some stations along closed streets show decreases as access becomes difficult.
  • Pride: The impact is highly concentrated in the West Village and Chelsea, with ridership spiking >50% at key stations like Christopher St-Sheridan Sq and 14 St. The effect is tightly localized—Midtown stations are largely unaffected.

These patterns suggest that event planning could benefit from station-specific service adjustments rather than system-wide changes.

Story 3: Station Personalities

Question: Do all subway stations behave the same way?

Not all subway stations serve the same purpose. Some are dominated by 9-to-5 commuters, others by evening leisure travelers, and some by late-night revelers. We used K-means clustering on hourly ridership profiles to classify stations into distinct “personality types.”

library(sf)

# Load station locations (need this here for map)
stations <- station_day %>%
  distinct(station_complex_id, station_complex) %>%
  left_join(read_csv("data/stations/mta_station_complexes.csv", show_col_types = FALSE) %>%
    select(station_complex_id, latitude, longitude),
  by = "station_complex_id"
  )

# Load NYC neighborhoods for map background
# Data preloaded by data/download_nyc_geo.R
nyc_nta <- st_read("data/geo/nyc_nta.geojson", quiet = TRUE)

if ("borough" %in% names(nyc_nta)) {
  manhattan_sf <- nyc_nta %>% filter(borough == "Manhattan")
} else if ("boro_name" %in% names(nyc_nta)) {
  manhattan_sf <- nyc_nta %>% filter(boro_name == "Manhattan")
} else {
  manhattan_sf <- nyc_nta
}

# 1. Create Hourly Profiles (Normalize by total daily rides)
hourly_profiles <- station_hour %>%
  filter(hour >= 5) %>% # Focus on active hours
  group_by(station_complex_id, hour) %>%
  summarize(avg_ridership = median(total_ridership, na.rm = TRUE), .groups = "drop_last") %>%
  mutate(profile = avg_ridership / sum(avg_ridership)) %>%
  ungroup() %>%
  select(station_complex_id, hour, profile) %>%
  pivot_wider(names_from = hour, values_from = profile, values_fill = 0)

# 2. K-Means Clustering
set.seed(123)
# Scale not needed as rows sum to 1 (compositional), but let's just use raw profile
kmeans_res <- kmeans(hourly_profiles %>% select(-station_complex_id), centers = 2)

station_clusters <- hourly_profiles %>%
  select(station_complex_id) %>%
  mutate(cluster = factor(kmeans_res$cluster)) %>%
  left_join(stations %>% select(station_complex_id, station_complex), by = "station_complex_id")

# Label clusters based on peak time (simple heuristic)
# (We'd inspect this manually in a real workflow, here we infer)
# Let's visualize the centroids to label them
centroids <- kmeans_res$centers %>%
  as.data.frame() %>%
  mutate(cluster = factor(1:nrow(kmeans_res$centers))) %>%
  pivot_longer(cols = -cluster, names_to = "hour", values_to = "pct") %>%
  mutate(hour = as.numeric(hour))

# Auto-label based on profile shape
# "Residential" stations have strong morning peak (commuters leaving home)
# "Commercial" stations have evening peak (workers leaving offices)
cluster_labels <- centroids %>%
  group_by(cluster) %>%
  summarize(
    morning_peak = sum(pct[hour %in% 7:9]),
    evening_peak = sum(pct[hour %in% 17:19]),
    .groups = "drop"
  ) %>%
  mutate(label = if_else(
    morning_peak > 0.15, "Residential", "Commercial"
  )) %>%
  select(cluster, label)

station_clusters <- station_clusters %>%
  left_join(cluster_labels, by = "cluster")

# Count stations by type
station_counts <- station_clusters %>%
  count(label) %>%
  filter(!is.na(label))

# Plot profiles
p_profiles <- ggplot(centroids %>% left_join(cluster_labels, by = "cluster") %>% arrange(label, hour),
  aes(x = hour, y = pct, color = label, group = label)
) +
  geom_line(linewidth = 1.2) +
  scale_color_manual(
    values = c("Residential" = okabe_ito$vermillion, "Commercial" = okabe_ito$blue)
  ) +
  labs(title = "Station Personality Profiles", x = "Hour of Day", y = "% of Daily Rides", color = "Station Type")

p_profiles

Finding: K-means clustering reveals two clear station archetypes based on commute flow direction:

  • Residential (43 stations): Sharp morning peak around 8 AM as commuters leave their neighborhoods to head to work. These are origin stations where people begin their daily commute.

  • Commercial (112 stations): Gradual buildup throughout the day with a pronounced evening peak (5-7 PM) as workers leave offices and head home. The activity stays elevated later into the evening as people also leave restaurants, theaters, and shops.

# Map of station types
station_map_data <- station_clusters %>%
  left_join(stations, by = c("station_complex_id", "station_complex")) %>%
  filter(!is.na(latitude) & !is.na(label))

p_map <- ggplot() +
  geom_sf(data = manhattan_sf, fill = "gray92", color = "gray75", linewidth = 0.2) +
  geom_point(
    data = station_map_data,
    aes(x = longitude, y = latitude, color = label),
    alpha = 0.8, size = 2.5
  ) +
  scale_color_manual(
    values = c("Residential" = okabe_ito$vermillion, "Commercial" = okabe_ito$blue),
    name = "Station Type"
  ) +
  coord_sf(
    crs = 4326,
    xlim = c(-74.02, -73.91),
    ylim = c(40.70, 40.88)
  ) +
  labs(
    title = "Geographic Distribution of Station Types",
    subtitle = paste0(
      "Workplace: ", station_counts$n[station_counts$label == "Workplace"], " stations | ",
      "Leisure/Commercial: ", station_counts$n[station_counts$label == "Leisure/Commercial"], " stations"
    ),
    x = NULL, y = NULL
  ) +
  theme_minimal() +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank(),
    legend.position = "bottom"
  )

p_map

Geographic Distribution: The map confirms our interpretation of commute flow. Commercial stations (blue) dominate Midtown and Lower Manhattan—Grand Central, Bryant Park, Rockefeller Center, and the Financial District—where workers exit the subway in the evening to head home. Residential stations (red) are concentrated in the Upper West Side, Upper East Side, and other neighborhood areas where commuters enter the subway each morning.

The few red spots in Midtown are notable exceptions that prove the rule: Penn Station and Times Square show morning-peak patterns not because they’re residential, but because they’re major transit hubs—commuters from New Jersey and Long Island arrive here in the morning to transfer to other lines or walk to nearby offices. Roosevelt Island is an actual residential neighborhood. These stations function as “entry points” into Manhattan, behaving like residential stations despite their Midtown location.

Story 4: Rain Sensitivity by Hour

Question: When does rain hurt ridership the most?

Beyond station types, we can examine how rain sensitivity varies throughout the day. This analysis reveals important patterns about when commuters are most likely to change their behavior in response to weather.

# Get precipitation data from station_day
daily_prcp <- station_day %>%
  select(station_complex_id, date, prcp) %>%
  distinct()

# Calculate avg drop in ridership on rainy days (>0.5 in) for each station by hour
rain_effect_hourly <- station_hour %>%
  filter(hour >= 5) %>%
  left_join(daily_prcp, by = c("station_complex_id", "date")) %>%
  mutate(is_rainy = if_else(prcp >= 0.5, "Rain", "Dry")) %>%
  group_by(station_complex_id, hour, is_rainy) %>%
  summarize(avg_ridership = median(total_ridership, na.rm = TRUE), .groups = "drop") %>%
  pivot_wider(names_from = is_rainy, values_from = avg_ridership) %>%
  filter(!is.na(Rain) & !is.na(Dry)) %>%
  mutate(rain_change_pct = (Rain - Dry) / Dry)

# Aggregate by hour (all stations combined)
rain_by_hour <- rain_effect_hourly %>%
  group_by(hour) %>%
  summarize(
    median_change = median(rain_change_pct, na.rm = TRUE),
    .groups = "drop"
  )

p_rain <- ggplot(rain_by_hour, aes(x = hour, y = median_change)) +
  geom_line(linewidth = 1.2, color = okabe_ito$blue) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  scale_y_continuous(labels = percent_format()) +
  labs(
    title = "Rain Sensitivity by Hour",
    subtitle = "Negative values = fewer riders on rainy days",
    x = "Hour of Day", y = "% Change on Rainy Days"
  )

p_rain

Finding: The rain sensitivity by hour plot shows that subway ridership drops most sharply during the morning commute hours (7-9 AM), with declines of 20-30%. This is likely because morning commuters face strict arrival deadlines—they must get to work on time, so if it’s raining they switch to alternatives like taxis, rideshares, or working from home. In contrast, afternoon and evening ridership shows smaller drops (5-18%) because travelers have more flexibility to simply wait for the rain to stop before taking the subway.

These insights suggest that the MTA should anticipate reduced morning rush demand on rainy days and consider adjusting service accordingly, while afternoon and evening service can remain relatively stable.

Summary

Our four-story analysis reveals distinct patterns in how weather, events, and station characteristics shape Manhattan subway ridership:

Weather Effects: Ridership peaks in a “Goldilocks zone” of 60-75°F with no rain, where ridership exceeds baseline by 4-6%. Precipitation has a stronger marginal effect than temperature—a half-inch of rain causes more ridership loss than a 20°F temperature swing. Cold, wet days see the largest drops (20-35% below average).

Event Impacts: Major events create geographically distinct “blast radii.” The NYC Marathon affects stations along the route and near the finish line, while the Pride March impact is tightly concentrated in the West Village and Chelsea. These patterns suggest station-specific rather than system-wide service adjustments would be most effective.

Station Personalities: K-means clustering identifies two station archetypes based on commute flow. Residential stations (concentrated in Upper Manhattan neighborhoods) show sharp morning peaks as commuters depart for work. Commercial stations (dominating Midtown and Lower Manhattan) show evening peaks as workers head home. Major transit hubs like Penn Station behave like residential stations despite their Midtown location.

Rain Timing: Rain sensitivity peaks during morning commute hours (7-9 AM) with 20-30% ridership drops, as commuters with strict arrival deadlines switch to alternatives. Afternoon and evening drops are smaller (5-18%) because travelers can wait out the rain.

Together, these findings provide actionable insights for demand forecasting, service planning, and understanding the daily rhythm of New York City transit.

P8105 Data Science I | Fall 2025

Columbia University Mailman School of Public Health

Department of Biostatistics