Analysis

59 - National Service Elasticity of Ridership

Coverage: Coverage window unavailable for this page.

Built 2026-06-15 11:52 UTC · Commit e5cf673

Page Navigation

Analysis Navigation

Data Provenance

flowchart LR
  59_national_service_elasticity(["59 - National Service Elasticity of Ridership"])
  t_ntd_annual_service[("ntd_annual_service")] --> 59_national_service_elasticity
  06_ntd_service[["NTD Annual Service ETL"]] --> t_ntd_annual_service
  u1_06_ntd_service[/"data/ntd-annual-service/2023_TS2.2_Service_Data.xlsx"/] --> 06_ntd_service
  d1_59_national_service_elasticity(("numpy (lib)")) --> 59_national_service_elasticity
  d2_59_national_service_elasticity(("polars (lib)")) --> 59_national_service_elasticity
  d3_59_national_service_elasticity(("statsmodels (lib)")) --> 59_national_service_elasticity
  d4_59_national_service_elasticity(("scipy (lib)")) --> 59_national_service_elasticity
  classDef page fill:#dbeafe,stroke:#1d4ed8,color:#1e3a8a,stroke-width:2px;
  classDef table fill:#ecfeff,stroke:#0e7490,color:#164e63;
  classDef dep fill:#fff7ed,stroke:#c2410c,color:#7c2d12,stroke-dasharray: 4 2;
  classDef file fill:#eef2ff,stroke:#6366f1,color:#3730a3;
  classDef api fill:#f0fdf4,stroke:#16a34a,color:#14532d;
  classDef pipeline fill:#f5f3ff,stroke:#7c3aed,color:#4c1d95;
  class 59_national_service_elasticity page;
  class t_ntd_annual_service table;
  class d1_59_national_service_elasticity,d2_59_national_service_elasticity,d3_59_national_service_elasticity,d4_59_national_service_elasticity dep;
  class u1_06_ntd_service file;
  class 06_ntd_service pipeline;

Findings

Findings: National Service Elasticity of Ridership

Summary

Across 8,323 year-over-year changes at 497 sizeable US transit agencies (1991–2024, pandemic years set aside), a 1% change in service hours is associated with about a 0.48% change in ridership in the same direction, the same year. Cutting service reliably loses riders — agencies that cut service by 10% or more lost a median of 6.3% of ridership — but ridership moves much less than one-for-one with service, and service changes explain only a small share (R² ≈ 0.12) of year-to-year ridership swings. Most of what drives ridership lies outside the agency's service-level dial. Pittsburgh (PRT) sits squarely on this national pattern.

Key Numbers

Headline elasticity: 0.48 (95% CI [0.46, 0.52]) — % ridership change per 1% service change, non-pandemic years, sizeable agencies (n = 8,323).
Cuts only: 0.36 (95% CI [0.30, 0.42], n = 2,945) — among service-cut episodes alone the slope is a bit shallower.
Sensitivity (base VRH ≥ 250,000): 0.52 — essentially unchanged, so the result is not an artifact of the size threshold.
R² ≈ 0.12 — service change explains ~12% of the variance in same-year ridership change. Pearson r = 0.35 (p < 1e-200).
Dose–response gradient (median ridership change): deep cut (≤ −10%) → −6.3%; cut (−10 to −5%) → −2.6%; mild cut (−5 to 0%) → −0.9%; flat (0 to +5%) → +0.9%; growth (> +5%) → +5.5%. Monotonic across all five buckets.
Pandemic transitions (2020–21): slope ≈ 1.0 — service and ridership co-moved far more tightly, confirming the pandemic was a different regime and correctly excluded.

Observations

Service cuts cost riders, but less than proportionally. The 0.48 slope means an agency cutting 10% of service hours typically loses about 5% of ridership in the same year. This matches the published short-run transit service elasticity range (~0.3–0.5), which is reassuring for the data.
The relationship is stable across three decades. The yearly-fit elasticity hovers around 0.3–0.6 from 1992 through 2019 with no trend — the service–ridership link is a structural feature, not a recent phenomenon. Only 2020–2021 break out (slope ≈ 1.0+).
Service is a minor driver of year-to-year ridership. With R² ≈ 0.12, roughly seven-eighths of the variation in annual ridership change is unrelated to that year's service change — reflecting fuel prices, the economy, fares, demographics, and competing travel modes. Restoring service helps, but it is not the main lever on ridership.
The dose–response curve is clean and monotonic. Bigger cuts → bigger ridership losses, bigger expansions → bigger gains, with no reversal. This is the clearest single picture of the relationship.
PRT is a typical agency on this measure. Across its history PRT's service-vs- ridership changes scatter around the national line. Its pre-pandemic cut years (e.g. 1992: −9% VRH / −11% UPT; 2003: −5.5% / −7.5%) track the national slope closely. Its 2022–2023 ridership rebounds (+44%, +17% on roughly flat service) are pandemic- recovery effects, not a service response.

Discussion

The dominant message echoes Analysis 39 from a much larger evidence base: service level is a real but modest lever on ridership. An agency that cuts service will lose riders along a predictable ~0.5 slope, and a deep cut compounds into a feedback risk (less service → longer waits → riders leave → cuts look justified). But the low R² means service restoration alone cannot recover ridership when the larger demand drivers have shifted — which is exactly PRT's post-2019 situation, where ridership fell far more than service. The elasticity quantifies the portion of ridership that is in the agency's hands: cutting 10% of service is worth roughly 5% of riders, no small thing, but the other half of any large ridership swing comes from forces the service budget does not touch.

Caveats

Association, not causation. Service and ridership are jointly determined. Agencies often cut service because ridership is already falling (reverse causality), and both respond to shared shocks. The 0.48 slope is descriptive co-movement, not a controlled "cut X% → lose Y%" causal estimate. The shallower cuts-only slope (0.36) is consistent with this — cuts are often reactive to demand already softening.
Ecological / agency-level. Results describe agencies, not individual riders. No inference about why any particular person stopped riding is supported.
System-level, all modes. VRH and UPT aggregate bus, rail, and demand-response into one figure per agency. Mode-specific elasticities (a rail cut vs a bus cut) cannot be separated here.
Recovery years remain in the sample. Only 2020 and 2021 are flagged pandemic; 2022–2023 ridership rebounds (large positive intercepts in those years' fits) stay in the analytic set and, if anything, attenuate the slope. Excluding them would raise the estimate slightly, so 0.48 is conservative.
Contemporaneous only. This measures same-year co-movement; lagged effects (a cut this year shedding riders over several years) are not captured and would tend to make the true multi-year elasticity larger.
VRH measures scheduled service, not quality. Reliability, frequency, and coverage changes that don't move total VRH are invisible to this metric.

Validation

Data source verified. VRH and UPT from ntd_annual_service (NTD TS2.2 workbooks, pipeline 06). Columns confirmed against data/DATA_DICTIONARY.md. PRT (ntd_id 30022) annual values spot-checked against the raw table.
Geographic/temporal scope. Both metrics come from the same table, same agency-year grain; changes use strictly consecutive years only (gap years excluded).
Null/zero handling. Only agency-years with VRH > 0 and UPT > 0 in both endpoints enter a change pair; nulls and zeros are dropped, not treated as observations.
Aggregates sanity-checked. Headline elasticity (0.48) lands inside the published short-run transit service-elasticity range (~0.3–0.5). PRT's annual VRH/UPT match Analysis 39's figures.
Direction of effects checked. Slope is positive (cut service → lose riders, grow service → gain riders) — the expected direction. The dose–response gradient is monotonic with no reversal.
Surprising results investigated. The pandemic subset's near-1.0 slope was investigated and is expected (joint exogenous collapse); it is reported separately, not folded into the headline.
Small-sample handling. Agencies below 100,000 base VRH are excluded; the per-year elasticity chart further drops years with n < 10. The 250,000-VRH sensitivity fit confirms the estimate is not threshold-driven.
Ecological framing. All statements are framed as agency-level associations.

Output

image service_vs_ridership_scatter.png
Each point is one agency-year; service change (x) vs ridership change (y) with the fitted 0.48 elasticity line and a 1:1 reference. Pandemic years excluded.

image dose_response_bins.png
Median ridership change (with interquartile range) grouped by size of service change, from deep cuts to growth — a clean monotonic gradient.

image elasticity_over_time.png
The fitted service-to-ridership elasticity estimated separately each year (1991-2024), showing it holds steady near 0.5 except during the pandemic.

image prt_cut_episodes.png
Pittsburgh (PRT) year-over-year service and ridership changes plotted against the national pattern, with cut, growth, and pandemic years marked.

No interactive outputs declared.

data agency_year_changes.csv

The per-transition panel — agency, year, prior and current VRH/UPT, and percent changes — used for all fits.

Preview CSV

Expand to load preview.

data elasticity_by_year.csv

Per-year fitted elasticity slope, 95% confidence interval, R-squared, and sample size.

Preview CSV

Expand to load preview.

data summary.csv

Headline elasticity estimates for the all-years, cuts-only, size-sensitivity, and pandemic subsets.

Preview CSV

Expand to load preview.

Methods

Methods: National Service Elasticity of Ridership

Question

Nationally, across all US transit agencies and the full NTD history (1991–2024), when an agency cuts (or expands) service, how much does its ridership move in the same year? Put differently: what is the agency-level service elasticity of ridership — the percent change in unlinked passenger trips (UPT) associated with a 1% change in vehicle revenue hours (VRH) — and how stable is it over time?

This generalizes Analysis 39, which compared a single 2019→2024 snapshot. Here we use every consecutive-year change across three decades and many cut episodes, not just the COVID window.

Approach

Build a year-over-year change panel. For each agency, take every pair of consecutive reporting years where VRH and UPT are both present and positive in both years. For each pair compute vrh_pct = ΔVRH/VRH_prev × 100 and upt_pct = ΔUPT/UPT_prev × 100. One row per agency-year-transition.
Restrict to sizeable agencies. Keep transitions where the base-year VRH ≥ 100,000 hours (~the largest ~340 agencies in 2019). Tiny rural systems produce wild percent swings on small denominators that would dominate the regression. The threshold is reported and a sensitivity check at 250,000 is run.
Flag and separate the pandemic. Transitions ending in 2020 or 2021 are marked pandemic. The pandemic collapsed service and ridership jointly for exogenous public-health reasons, not as a service→ridership response, so these transitions are excluded from the headline elasticity and reported separately.
Trim reporting discontinuities. Drop transitions with |vrh_pct| > 50% or |upt_pct| > 75% as likely mergers, definition changes, or first-year reporting artifacts rather than genuine service decisions. The count dropped is reported.
Estimate elasticity (OLS). Regress upt_pct ~ vrh_pct on the non-pandemic, trimmed panel. The slope is the elasticity; report slope, 95% CI, R², n, and the Pearson correlation. Refit on the cuts-only subset (vrh_pct < 0), since the question is specifically about service reductions.
Elasticity over time. Fit the same slope separately within each year (1991–2024) to check whether the service–ridership relationship is stable, strengthening, or weakening; pandemic years are flagged on the chart.
Dose–response bins. Bucket transitions by service-change size (deep cut ≤ −10%, cut −10 to −5%, mild cut −5 to 0%, roughly flat 0 to +5%, growth > +5%) and show the distribution (median, IQR) of ridership change in each bucket.
PRT overlay. Plot PRT's own consecutive-year changes against the national regression line, highlighting its service-cut years.

Framing. This is an agency-level association, not an individual rider behavior model and not a causal demand elasticity. Service and ridership are jointly determined — agencies frequently cut service because ridership is already falling (reverse causality), and both respond to shared shocks (fuel prices, recessions, local economies). The slope describes how service and ridership co-move across agencies and years; it is not a controlled estimate of "cut service by 1% → lose X% of riders."

Data

Name	Description	Source
`ntd_annual_service`	Annual system-level VRH, UPT per agency (1991–2024)	`prt.db` table (pipeline 06, NTD TS2.2 workbooks)

Columns used: ntd_id, agency_name, year, vrh, upt. System-level (all modes aggregated). Inclusion: VRH & UPT present and > 0 in both years of a transition; base-year VRH ≥ 100,000.

Output

output/service_vs_ridership_scatter.png — Δ%VRH vs Δ%UPT for all non-pandemic sizeable transitions, with the fitted elasticity line, y=x reference, and the service-cut region shaded.
output/dose_response_bins.png — distribution of ridership change within each service-change bucket (median, IQR, n per bin).
output/elasticity_over_time.png — per-year fitted elasticity slope (1991–2024) with 95% CI band; pandemic years flagged.
output/prt_cut_episodes.png — PRT's consecutive-year service vs ridership changes over its history, cut years highlighted, against the national line.
output/agency_year_changes.csv — the per-transition change panel used.
output/elasticity_by_year.csv — per-year slope, CI, n, R².
output/summary.csv — headline elasticity estimates (all, cuts-only, sensitivity, pandemic-subset).

  Source Code
    
      
      """Analysis 59: estimate the agency-level service elasticity of ridership across the full NTD panel (1991-2024)."""

import numpy as np
import polars as pl
import statsmodels.api as sm

from prt_otp_analysis.common import (
    COLOR_GRAY,
    COLOR_NEGATIVE,
    COLOR_PRIMARY,
    analysis_dir,
    correlate,
    get_db,
    phase,
    run_analysis,
    save_chart,
    save_csv,
    setup_plotting,
    validate,
)
from prt_otp_analysis.common.schemas import NTD_ANNUAL_SERVICE

OUT = analysis_dir(__file__)

PRT_NTD_ID = 30022
MIN_BASE_VRH = 100_000        # exclude tiny agencies whose % swings dominate
SENSITIVITY_VRH = 250_000     # robustness threshold
VRH_TRIM_PCT = 50.0           # drop |vrh_pct| above this (reporting discontinuities)
UPT_TRIM_PCT = 75.0           # drop |upt_pct| above this
PANDEMIC_END_YEARS = (2020, 2021)


def load_changes(conn) -> pl.DataFrame:
    """Build the consecutive-year change panel: one row per agency-year transition.

    Returns columns: ntd_id, agency_name, year (end year), vrh_prev, vrh,
    upt_prev, upt, vrh_pct, upt_pct, base_vrh, pandemic.
    """
    rows = conn.execute("""
        SELECT ntd_id, agency_name, year, vrh, upt
        FROM ntd_annual_service
        WHERE vrh > 0 AND upt > 0
        ORDER BY ntd_id, year
    """).fetchall()
    df = pl.DataFrame([dict(r) for r in rows])
    validate(df, NTD_ANNUAL_SERVICE, subset=True)

    # Lag within agency; only keep strictly consecutive years (year - prev == 1).
    df = df.sort("ntd_id", "year").with_columns(
        vrh_prev=pl.col("vrh").shift(1).over("ntd_id"),
        upt_prev=pl.col("upt").shift(1).over("ntd_id"),
        year_prev=pl.col("year").shift(1).over("ntd_id"),
    )
    df = df.filter(
        pl.col("year_prev").is_not_null() & (pl.col("year") - pl.col("year_prev") == 1)
    )

    df = df.with_columns(
        base_vrh=pl.col("vrh_prev"),
        vrh_pct=(pl.col("vrh") - pl.col("vrh_prev")) / pl.col("vrh_prev") * 100,
        upt_pct=(pl.col("upt") - pl.col("upt_prev")) / pl.col("upt_prev") * 100,
        pandemic=pl.col("year").is_in(PANDEMIC_END_YEARS),
    )
    return df


def fit_elasticity(change_df: pl.DataFrame) -> dict[str, float]:
    """OLS of upt_pct ~ vrh_pct. Returns slope (elasticity), intercept, CI, R², n."""
    x = change_df["vrh_pct"].to_numpy()
    y = change_df["upt_pct"].to_numpy()
    n = len(x)
    if n < 3:
        nan = float("nan")
        return {"slope": nan, "intercept": nan, "ci_lo": nan, "ci_hi": nan, "r2": nan, "n": n}
    model = sm.OLS(y, sm.add_constant(x)).fit()
    ci = model.conf_int(alpha=0.05)
    return {
        "slope": float(model.params[1]),
        "intercept": float(model.params[0]),
        "ci_lo": float(ci[1][0]),
        "ci_hi": float(ci[1][1]),
        "r2": float(model.rsquared),
        "n": n,
    }


@run_analysis(59, "National Service Elasticity of Ridership")
def main():
    plt = setup_plotting()
    conn = get_db()

    with phase("Building year-over-year change panel"):
        raw_df = load_changes(conn)
        conn.close()
        print(f"   {len(raw_df):,} consecutive-year transitions (VRH & UPT > 0, all agencies)")

        sized_df = raw_df.filter(pl.col("base_vrh") >= MIN_BASE_VRH)
        print(f"   {len(sized_df):,} transitions from agencies with base VRH >= {MIN_BASE_VRH:,}")

        # Trim reporting discontinuities.
        trimmed_df = sized_df.filter(
            (pl.col("vrh_pct").abs() <= VRH_TRIM_PCT) & (pl.col("upt_pct").abs() <= UPT_TRIM_PCT)
        )
        dropped = len(sized_df) - len(trimmed_df)
        print(f"   Dropped {dropped} transitions exceeding trim bounds "
              f"(|VRH|>{VRH_TRIM_PCT:.0f}% or |UPT|>{UPT_TRIM_PCT:.0f}%)")

        analytic_df = trimmed_df.filter(~pl.col("pandemic"))
        pandemic_df = trimmed_df.filter(pl.col("pandemic"))
        cuts_df = analytic_df.filter(pl.col("vrh_pct") < 0)
        print(f"   Analytic (non-pandemic) sample: {len(analytic_df):,} transitions "
              f"({analytic_df['ntd_id'].n_unique()} agencies, "
              f"{analytic_df['year'].min()}-{analytic_df['year'].max()})")
        print(f"   Of those, {len(cuts_df):,} are service cuts (VRH change < 0)")

    # ---- Headline elasticity fits ----
    with phase("Estimating elasticity"):
        all_fit = fit_elasticity(analytic_df)
        cuts_fit = fit_elasticity(cuts_df)
        sens_fit = fit_elasticity(
            trimmed_df.filter((pl.col("base_vrh") >= SENSITIVITY_VRH) & ~pl.col("pandemic"))
        )
        pandemic_fit = fit_elasticity(pandemic_df)
        corr = correlate(analytic_df, "vrh_pct", "upt_pct")

        print(f"\n   All non-pandemic transitions  : elasticity = {all_fit['slope']:.3f} "
              f"[{all_fit['ci_lo']:.3f}, {all_fit['ci_hi']:.3f}]  "
              f"R²={all_fit['r2']:.3f}  n={all_fit['n']:,}")
        print(f"   Service cuts only             : elasticity = {cuts_fit['slope']:.3f} "
              f"[{cuts_fit['ci_lo']:.3f}, {cuts_fit['ci_hi']:.3f}]  "
              f"R²={cuts_fit['r2']:.3f}  n={cuts_fit['n']:,}")
        print(f"   Sensitivity (base VRH>={SENSITIVITY_VRH:,}) : elasticity = {sens_fit['slope']:.3f} "
              f"[{sens_fit['ci_lo']:.3f}, {sens_fit['ci_hi']:.3f}]  n={sens_fit['n']:,}")
        print(f"   Pandemic transitions (2020-21): elasticity = {pandemic_fit['slope']:.3f} "
              f"[{pandemic_fit['ci_lo']:.3f}, {pandemic_fit['ci_hi']:.3f}]  n={pandemic_fit['n']:,}")
        print(f"   Pearson r (non-pandemic): {corr['pearson_r']:.3f} (p={corr['pearson_p']:.1e})")

        summary_df = pl.DataFrame([
            {"sample": "all_non_pandemic", **all_fit},
            {"sample": "cuts_only", **cuts_fit},
            {"sample": f"sensitivity_vrh_{SENSITIVITY_VRH}", **sens_fit},
            {"sample": "pandemic_2020_2021", **pandemic_fit},
        ])

    # ---- Per-year elasticity ----
    with phase("Fitting per-year elasticity"):
        year_rows = []
        for yr in sorted(trimmed_df["year"].unique().to_list()):
            yr_df = trimmed_df.filter(pl.col("year") == yr)
            fit = fit_elasticity(yr_df)
            year_rows.append({"year": yr, "pandemic": yr in PANDEMIC_END_YEARS, **fit})
        year_df = pl.DataFrame(year_rows)

    # ---- Save CSVs ----
    print()
    save_csv(
        analytic_df.select(
            "ntd_id", "agency_name", "year", "vrh_prev", "vrh", "upt_prev", "upt",
            "vrh_pct", "upt_pct", "base_vrh",
        ).sort("year", "ntd_id"),
        OUT / "agency_year_changes.csv",
    )
    save_csv(year_df, OUT / "elasticity_by_year.csv")
    save_csv(summary_df, OUT / "summary.csv")

    # ---- Chart 1: scatter with elasticity line ----
    with phase("Generating service-vs-ridership scatter"):
        x = analytic_df["vrh_pct"].to_numpy()
        y = analytic_df["upt_pct"].to_numpy()

        fig, ax = plt.subplots(figsize=(11, 9))
        ax.axvspan(-VRH_TRIM_PCT, 0, color=COLOR_NEGATIVE, alpha=0.05, zorder=0)
        ax.text(-VRH_TRIM_PCT + 1, UPT_TRIM_PCT - 4, "service cut →", color=COLOR_NEGATIVE,
                fontsize=9, alpha=0.8, va="top")

        ax.scatter(x, y, s=14, alpha=0.30, color=COLOR_PRIMARY, edgecolors="none", zorder=2)

        line_x = np.array([-VRH_TRIM_PCT, VRH_TRIM_PCT])
        ax.plot(line_x, all_fit["intercept"] + all_fit["slope"] * line_x,
                color=COLOR_NEGATIVE, linewidth=2.5, zorder=4,
                label=f"Elasticity = {all_fit['slope']:.2f}  (R²={all_fit['r2']:.2f}, n={all_fit['n']:,})")
        ax.plot(line_x, line_x, color=COLOR_GRAY, linestyle="--", linewidth=1, zorder=3,
                label="1:1 (ridership tracks service exactly)")

        ax.axhline(0, color="#CCCCCC", linewidth=0.8, zorder=1)
        ax.axvline(0, color="#CCCCCC", linewidth=0.8, zorder=1)
        ax.set_xlim(-VRH_TRIM_PCT, VRH_TRIM_PCT)
        ax.set_ylim(-UPT_TRIM_PCT, UPT_TRIM_PCT)
        ax.set_xlabel("Service change, year over year — Vehicle Revenue Hours (%)")
        ax.set_ylabel("Ridership change, year over year — Unlinked Passenger Trips (%)")
        ax.set_title("When agencies change service, how does ridership move?\n"
                     f"US transit agencies, base VRH ≥ {MIN_BASE_VRH:,}, 1991–2024, pandemic years excluded")
        ax.legend(loc="lower right", fontsize=9)
        save_chart(fig, OUT / "service_vs_ridership_scatter.png")

    # ---- Chart 2: dose-response bins ----
    with phase("Generating dose-response bins"):
        bins = [
            ("Deep cut\n(≤ −10%)", pl.col("vrh_pct") <= -10),
            ("Cut\n(−10 to −5%)", (pl.col("vrh_pct") > -10) & (pl.col("vrh_pct") <= -5)),
            ("Mild cut\n(−5 to 0%)", (pl.col("vrh_pct") > -5) & (pl.col("vrh_pct") < 0)),
            ("Flat / small grow\n(0 to +5%)", (pl.col("vrh_pct") >= 0) & (pl.col("vrh_pct") <= 5)),
            ("Growth\n(> +5%)", pl.col("vrh_pct") > 5),
        ]
        labels, medians, q1s, q3s, ns = [], [], [], [], []
        for label, cond in bins:
            sub = analytic_df.filter(cond)["upt_pct"]
            if len(sub) == 0:
                continue
            labels.append(label)
            medians.append(sub.median())
            q1s.append(sub.quantile(0.25))
            q3s.append(sub.quantile(0.75))
            ns.append(len(sub))

        fig, ax = plt.subplots(figsize=(11, 7))
        xpos = range(len(labels))
        lower = [m - q for m, q in zip(medians, q1s)]
        upper = [q - m for q, m in zip(q3s, medians)]
        colors = [COLOR_NEGATIVE if m < 0 else COLOR_PRIMARY for m in medians]
        ax.bar(xpos, medians, color=colors, alpha=0.85,
               yerr=[lower, upper], capsize=6, error_kw={"alpha": 0.5, "lw": 1.2})
        ax.axhline(0, color="#999999", linewidth=0.9)
        for i, (m, n) in enumerate(zip(medians, ns)):
            ax.text(i, m + (1 if m >= 0 else -1), f"{m:+.1f}%\n(n={n})",
                    ha="center", va="bottom" if m >= 0 else "top", fontsize=8)
        ax.set_xticks(list(xpos))
        ax.set_xticklabels(labels, fontsize=9)
        ax.set_ylabel("Ridership change — median (bars), IQR (whiskers) (%)")
        ax.set_title("Ridership change by size of service change\n"
                     "US transit agencies, 1991–2024, pandemic years excluded")
        save_chart(fig, OUT / "dose_response_bins.png")

    # ---- Chart 3: elasticity over time ----
    with phase("Generating elasticity-over-time chart"):
        plot_df = year_df.filter(pl.col("n") >= 10).sort("year")
        yrs = plot_df["year"].to_list()
        slopes = plot_df["slope"].to_list()
        lo = plot_df["ci_lo"].to_list()
        hi = plot_df["ci_hi"].to_list()
        is_pan = plot_df["pandemic"].to_list()

        fig, ax = plt.subplots(figsize=(13, 6.5))
        ax.fill_between(yrs, lo, hi, color=COLOR_PRIMARY, alpha=0.15, zorder=1)
        ax.plot(yrs, slopes, color=COLOR_PRIMARY, linewidth=1.8, marker="o", markersize=4, zorder=2)
        for yr, s, pan in zip(yrs, slopes, is_pan):
            if pan:
                ax.scatter([yr], [s], color=COLOR_NEGATIVE, s=70, zorder=3, edgecolors="white")
        ax.axhline(all_fit["slope"], color=COLOR_GRAY, linestyle="--", linewidth=1.2,
                   label=f"Pooled non-pandemic elasticity = {all_fit['slope']:.2f}")
        ax.axhline(0, color="#CCCCCC", linewidth=0.8)
        ax.scatter([], [], color=COLOR_NEGATIVE, s=70, edgecolors="white",
                   label="Pandemic year (2020, 2021)")
        ax.set_xlabel("Year")
        ax.set_ylabel("Fitted elasticity (Δ%UPT per 1% Δ%VRH)")
        ax.set_title("Service elasticity of ridership over time\n"
                     "Slope of yearly Δridership-vs-Δservice fit, sizeable US agencies")
        ax.legend(loc="upper left", fontsize=9)
        save_chart(fig, OUT / "elasticity_over_time.png")

    # ---- Chart 4: PRT overlay ----
    with phase("Generating PRT cut-episodes chart"):
        prt_df = trimmed_df.filter(pl.col("ntd_id") == PRT_NTD_ID).sort("year")
        fig, ax = plt.subplots(figsize=(11, 9))

        bg = analytic_df
        ax.scatter(bg["vrh_pct"].to_numpy(), bg["upt_pct"].to_numpy(),
                   s=10, alpha=0.12, color=COLOR_GRAY, edgecolors="none", zorder=1)
        line_x = np.array([-30, 30])
        ax.plot(line_x, all_fit["intercept"] + all_fit["slope"] * line_x,
                color="#888888", linewidth=1.8, linestyle="--", zorder=2,
                label=f"National elasticity = {all_fit['slope']:.2f}")

        for row in prt_df.iter_rows(named=True):
            is_cut = row["vrh_pct"] < 0
            is_pan = row["pandemic"]
            color = COLOR_NEGATIVE if is_pan else (COLOR_PRIMARY if is_cut else "#6ACC65")
            ax.scatter(row["vrh_pct"], row["upt_pct"], s=90, color=color,
                       edgecolors="black", linewidths=0.8, zorder=4)
            ax.annotate(str(row["year"]), (row["vrh_pct"], row["upt_pct"]),
                        textcoords="offset points", xytext=(6, 5), fontsize=8)

        ax.scatter([], [], s=90, color=COLOR_PRIMARY, edgecolors="black", label="PRT service-cut year")
        ax.scatter([], [], s=90, color="#6ACC65", edgecolors="black", label="PRT service-growth year")
        ax.scatter([], [], s=90, color=COLOR_NEGATIVE, edgecolors="black", label="PRT pandemic year")
        ax.axhline(0, color="#CCCCCC", linewidth=0.8)
        ax.axvline(0, color="#CCCCCC", linewidth=0.8)
        ax.set_xlim(-30, 30)
        ax.set_ylim(-60, 30)
        ax.set_xlabel("PRT service change, year over year — VRH (%)")
        ax.set_ylabel("PRT ridership change, year over year — UPT (%)")
        ax.set_title("Pittsburgh (PRT) service vs ridership changes against the national pattern")
        ax.legend(loc="lower right", fontsize=9)
        save_chart(fig, OUT / "prt_cut_episodes.png")

    # ---- Console: PRT episode detail ----
    print("\n   PRT consecutive-year changes:")
    for row in prt_df.iter_rows(named=True):
        tag = " [pandemic]" if row["pandemic"] else (" [cut]" if row["vrh_pct"] < 0 else "")
        print(f"     {row['year']}: VRH {row['vrh_pct']:>+6.1f}%   UPT {row['upt_pct']:>+6.1f}%{tag}")


if __name__ == "__main__":
    main()

    

    

Sources

Name	Type	Why It Matters	Owner	Freshness	Caveat
ntd_annual_service	table	Primary analytical table used in this page's computations.	Produced by NTD Annual Service ETL.	Updated when the producing pipeline step is rerun.	Coverage depends on upstream source availability and ETL assumptions.
Upstream sources (1) file data/ntd-annual-service/2023_TS2.2_Service_Data.xlsx — NTD TS2.2 workbook with annual service data by system.
numpy	dependency	Runtime dependency required for this page's pipeline or analysis code.	Open-source Python ecosystem maintainers.	Version pinned by project environment until dependency updates are applied.	Library updates may change behavior or defaults.
polars	dependency	Runtime dependency required for this page's pipeline or analysis code.	Open-source Python ecosystem maintainers.	Version pinned by project environment until dependency updates are applied.	Library updates may change behavior or defaults.
statsmodels	dependency	Runtime dependency required for this page's pipeline or analysis code.	Open-source Python ecosystem maintainers.	Version pinned by project environment until dependency updates are applied.	Library updates may change behavior or defaults.
scipy	dependency	Runtime dependency required for this page's pipeline or analysis code.	Open-source Python ecosystem maintainers.	Version pinned by project environment until dependency updates are applied.	Library updates may change behavior or defaults.