"""Safer streets, except on two wheels: the analysis behind a SparkyData article. py src/build_article.py who-still-dies-on-the-streets Traffic deaths in New York are near the fewest the City has recorded. Who still dies, and how has that changed over thirty years? Three sources answer it: - NHTSA's Fatality Analysis Reporting System (card ``fars_nyc``): every death within 30 days of a crash involving a motor vehicle on a public road, by road user, 1995 to 2024 (final files to 2023, the Annual Report File for 2024), coded the same way every year except for e-bikes, which moved from mopeds to bicycles in 2022 (when crashes involving only e-bikes stopped being counted); so the long view keeps motorcycles, mopeds, motor scooters and bicycles together as "two wheels", a line that undercounts e-bike deaths from 2022, and the fall in all deaths is also given ending before 2022; stand-up scooters are personal conveyances, outside the line in every year; - the NYPD's collision files (card ``nypd_collisions``): every person recorded killed from 2013, the road user read from their own vehicle's record (e-bikes, e-scooters and mopeds apart from bicycles and motorcycles), to June 2026; every dead person's vehicle record is read (none is inferred from a missing one), the joins are checked, and the labels' history is shown by year, since e-bike and e-scooter labels are in use only from August 2019, and set against the federal record's own coding of mopeds and motorized bicycles to 2021; - the year to date: the NYPD's CompStat count of traffic fatalities to 11 October 2026, DOT's first-half release, and DOT's nine-month figures as amNY reported them on 2 October. FARS and the NYPD count differently (FARS: within 30 days, crashes with a motor vehicle, public roads; the NYPD: the deaths its crash reports hold) and DOT applies rules of its own, so the three never share a line. Counts of deaths carry no sampling margin; a change is named only when the approximate 90% interval of the ratio of two counts, as independent Poisson counts, excludes no change (``sparkylib.stats.count_change``). As a robustness check every named change is recomputed with the variance scaled by the extra year-to-year swing its series shows (``sparkylib.stats. dispersion``); a change that does not survive it is not called a rise or a fall in the story. Writes results.json, viz.json, charts.json, data.csv and build.json into output/articles/who_still_dies_on_the_streets/ and site/public/articles/who-still-dies-on-the-streets/, and copies this script there as analysis.py. Deterministic: no sampling. """ from __future__ import annotations import re import pandas as pd from article_kit import Article, lines, r from nycdata import crime, roads from sparkylib.stats import count_change, dispersion SLUG = "who-still-dies-on-the-streets" SEED = 20261014 EARLY = range(1995, 2000) # the first five years of the long view RECENT = range(2020, 2025) # the last five FIRST_NYPD, LAST_YEAR = 2013, 2025 #: The long view's groups: two wheels are FARS's motorcyclists and cyclists together, the series #: that stays whole when e-bikes moved from mopeds to bicycles in 2022. LONG_GROUPS = ["pedestrian", "vehicle occupant", "two wheels"] LONG_NAMES = {"pedestrian": "Pedestrians", "vehicle occupant": "In cars, vans, trucks and buses", "two wheels": "On two wheels (motorcycles, mopeds, motor scooters, bicycles)"} LONG_SHORT = {"pedestrian": "Pedestrians", "vehicle occupant": "In vehicles", "two wheels": "Two wheels"} TWO_WHEELS_NYPD = ["e-micromobility", "motorcycle", "bicycle"] COMPSTAT, COMPSTAT_WEEK = "nypd_compstat_citywide_2026_10_11", "2026-10-11" # the latest week before release MANUAL, MANUAL_PDF = "nhtsa_fars_manual_1975_2024", "nhtsa_fars_analytical_users_manual_1975_2024.pdf" def long_years() -> pd.DataFrame: """FARS deaths by year (rows) and road-user group (columns), with two wheels and the total.""" f = roads.fars_deaths() g = f.groupby(["year", "road_user"])["deaths"].sum().unstack(fill_value=0) g = g.reindex(columns=roads.FARS_GROUPS, fill_value=0) g["two wheels"] = g["motorcyclist"] + g["cyclist"] g["total"] = g[roads.FARS_GROUPS].sum(axis=1) return g.loc[roads.FARS_FIRST:roads.FARS_LAST] def nypd_years(k: pd.DataFrame) -> tuple[pd.DataFrame, pd.DataFrame]: """NYPD people killed by year and road-user group (own-vehicle classes), whole years 2013-2025, and January to May of every year to 2026 (the files' last complete months).""" k = k.assign(group=k["road_user"].map(roads.GROUP)) years = k[k["year"].between(FIRST_NYPD, LAST_YEAR)].groupby(["year", "group"]).size().unstack(fill_value=0) years = years.reindex(columns=roads.GROUPS_ORDER, fill_value=0) years["total"] = years.sum(axis=1) jm = k[(k["date"].dt.month <= 5) & (k["year"] >= FIRST_NYPD)] first5 = jm.groupby(["year", "group"]).size().unstack(fill_value=0).reindex(columns=roads.GROUPS_ORDER, fill_value=0) first5["total"] = first5.sum(axis=1) return years, first5 def quoted_figures(A: Article) -> dict: """Every figure the prose quotes from a report, found in its saved copy.""" Q = A.quotes() y25, h1, am = ("nyc_dot_traffic_deaths_2025_release.html", "nyc_dot_traffic_deaths_h1_2026_release.html", "amny_dot_traffic_deaths_2026-10-02.html") return { "dot_2025": Q.quote("nyc_dot_traffic_deaths_2025_release", y25, "There were 205 recorded traffic deaths in 2025"), "dot_2025_record": Q.quote("nyc_dot_traffic_deaths_2025_release", y25, "dating back to when record-keeping began in 1910"), "dot_2024": Q.quote("nyc_dot_traffic_deaths_2025_release", y25, "a 19 percent decline from 253"), "dot_peds_2025": Q.quote("nyc_dot_traffic_deaths_2025_release", y25, "from 122 pedestrian deaths in 2024 to 111 in 2025"), "dot_every_mode": Q.quote("nyc_dot_traffic_deaths_2025_release", y25, "Every major travel mode experienced a decline in traffic deaths in 2025"), "dot_motorcycles": Q.quote("nyc_dot_traffic_deaths_2025_release", y25, "and motorcycle users, by 32 percent"), "dot_h1_total": Q.quote("nyc_dot_traffic_deaths_h1_2026_release", h1, "there have been 94 traffic fatalities"), "dot_h1_total_2025": Q.quote("nyc_dot_traffic_deaths_h1_2026_release", h1, "91 and well below the average from 2021 to 2024"), "dot_h1_micro": Q.quote("nyc_dot_traffic_deaths_h1_2026_release", h1, "have increased from 16 to 18"), "dot_h1_micro_earlier": Q.quote(*roads.DOT_FIRST_HALF), "amny_peds": Q.quote("amny_dot_traffic_deaths_2026_10_02", am, "67 pedestrian traffic deaths over the first nine months"), "amny_peds_record": Q.quote("amny_dot_traffic_deaths_2026_10_02", am, "the safest for pedestrians since the city began tallying traffic deaths in 1910"), "amny_micro": Q.quote("amny_dot_traffic_deaths_2026_10_02", am, "17 deaths among e-bike riders and 11 e-scooter rider fatalities"), "amny_illegal": Q.quote("amny_dot_traffic_deaths_2026_10_02", am, "at least two of the 17 e-bike users who died this year were riding illegal devices"), "amny_illegal_scooters": Q.quote("amny_dot_traffic_deaths_2026_10_02", am, "six of the e-scooter rider deaths involved illegal two-wheelers"), "compstat_traffic": Q.quote(COMPSTAT, "nypd_compstat_citywide_2026-10-11.pdf", "Traffic Fatalities 7 6 16.7 20 19 5.3 151 155 -2.6"), # what FARS counts, from NHTSA's manual: the case criteria, and e-bikes from 2022 "fars_case": Q.quote(MANUAL, MANUAL_PDF, "the crash had to involve a motor vehicle traveling on a trafficway " "customarily open to the public"), "fars_ebikes": Q.quote(MANUAL, MANUAL_PDF, "crashes involving only motorized bicycles were no longer collected"), } def main(argv: list[str] | None = None) -> int: """Run the analysis and write every file (no options).""" A = Article(SLUG, __file__, seed=SEED) claimed: dict[str, str] = {} def claim(path: str, stated: str) -> None: """A direction the prose states (A.claim), kept here too for the robustness check.""" A.claim(path, stated) claimed[path] = stated # ---- this year: the NYPD's weekly count and DOT's releases cs = crime.compstat(week=COMPSTAT_WEEK) tf = cs["rows"]["Traffic Fatalities"] A.results["now"] = {"week_end": cs["week_end"], "traffic_fatalities": count_change(tf["ytd"], tf["ytd_prior"])} if A.results["now"]["traffic_fatalities"]["direction"] != "unclear": raise AssertionError("the prose says this year's traffic deaths are about level with last year's") A.printed("now.traffic_fatalities.now", fmt="int") A.printed("now.traffic_fatalities.before", fmt="int") A.latest(cs["week_end"], "NYPD CompStat citywide report, traffic fatalities year to date") quotes = quoted_figures(A) A.results["quotes"] = quotes def num(key: str, pattern: str) -> int: """A figure read from its saved quote, so the prose and the source cannot drift apart.""" m = re.search(pattern, quotes[key]["quote"]) if not m: raise AssertionError(f"{key}: {pattern!r} not in the quote") return int(m.group(1).replace(",", "")) dot = {"deaths_2025": num("dot_2025", r"There were (\d+) recorded traffic deaths"), "deaths_2024": num("dot_2024", r"decline from (\d+)"), "peds_2025": num("dot_peds_2025", r"to (\d+) in 2025"), "h1_2026": num("dot_h1_total", r"there have been (\d+) traffic"), "h1_2025": num("dot_h1_total_2025", r"(\d+) and well below"), "h1_micro_2026": num("dot_h1_micro", r"from \d+ to (\d+)"), "h1_micro_2025": num("dot_h1_micro", r"from (\d+) to"), "h1_micro_2024": num("dot_h1_micro_earlier", r"the (\d+) recorded at this point in 2024"), "sep_peds_2026": num("amny_peds", r"(\d+) pedestrian"), "sep_ebike_2026": num("amny_micro", r"(\d+) deaths among e-bike"), "sep_escooter_2026": num("amny_micro", r"and (\d+) e-scooter")} A.results["dot"] = dot for key in ("deaths_2025", "deaths_2024", "h1_2026", "h1_2025", "h1_micro_2026", "h1_micro_2025", "sep_peds_2026", "sep_ebike_2026", "sep_escooter_2026"): A.printed(f"dot.{key}", fmt="int") A.results["dot_changes"] = {"h1_total": count_change(dot["h1_2026"], dot["h1_2025"]), "h1_micro": count_change(dot["h1_micro_2026"], dot["h1_micro_2025"])} for key, c in A.results["dot_changes"].items(): if c["direction"] != "unclear": raise AssertionError(f"the prose says this year's {key} is about the same as last year's") # ---- thirty years: FARS g = long_years() early, recent = g.loc[EARLY.start:EARLY.stop - 1], g.loc[RECENT.start:RECENT.stop - 1] cols = roads.FARS_GROUPS + ["two wheels", "total"] A.results["fars"] = { "years": {str(y): {c: int(g.loc[y, c]) for c in cols} for y in g.index}, "early_years": [EARLY.start, EARLY.stop - 1], "recent_years": [RECENT.start, RECENT.stop - 1], "early": {c: r(early[c].mean()) for c in cols}, "recent": {c: r(recent[c].mean()) for c in cols}, "change": {c: count_change(recent[c].sum(), early[c].sum()) for c in cols}, "two_wheel_share": {str(y): r(100 * g.loc[y, "two wheels"] / g.loc[y, "total"]) for y in (roads.FARS_FIRST, roads.FARS_LAST)}, } fr = A.results["fars"] for c in ("pedestrian", "vehicle occupant", "total"): claim(f"fars.change.{c}", "down") claim("fars.change.two wheels", "up") if not 45 <= abs(fr["change"]["pedestrian"]["change_pct"]) <= 55: raise AssertionError("the prose says pedestrian deaths fell by about half") for c in ("vehicle occupant", "pedestrian", "two wheels", "total"): A.printed(f"fars.recent.{c}", f"{fr['recent'][c]:.0f}") A.printed(f"fars.early.{c}", f"{fr['early'][c]:.0f}") if not 45 <= abs(fr["change"]["pedestrian"]["change_pct"]) <= 55 or abs(fr["change"]["vehicle occupant"]["change_pct"]) <= 55: raise AssertionError("the prose says pedestrian deaths fell by about half and deaths in vehicles by more") # the shape of the two-wheel rise: a plateau to 2019, a step in 2020 plateau = g.loc[roads.FARS_FIRST:2019, "two wheels"] step = g.loc[2020:roads.FARS_LAST, "two wheels"] fr["two_wheels_plateau"] = {"first": roads.FARS_FIRST, "last": 2019, "mean": r(plateau.mean(), 0), "min": int(plateau.min()), "max": int(plateau.max())} fr["two_wheels_step"] = {"first": 2020, "last": roads.FARS_LAST, "min": int(step.min()), "max": int(step.max())} fr["two_wheels_step_vs_plateau"] = count_change(step.sum(), plateau.sum(), exposure_now=len(step), exposure_before=len(plateau)) claim("fars.two_wheels_step_vs_plateau", "up") A.printed("fars.two_wheels_plateau.mean", fmt="int") A.printed("fars.two_wheels_step.min", fmt="int") # the 2022 break (e-bikes moved, e-bike-only crashes dropped): the same comparison ending before it pre = g.loc[2020:2021, "two wheels"] fr["two_wheels_before_break"] = {"years": [2020, 2021], "mean": r(pre.mean()), "change": count_change(pre.sum(), early["two wheels"].sum(), exposure_now=len(pre), exposure_before=len(early))} claim("fars.two_wheels_before_break.change", "up") A.printed("fars.two_wheels_before_break.mean", f"{fr['two_wheels_before_break']['mean']:.0f}") # the break dropped crashes involving only e-bikes, which can only lower the later years, so the fall in # all deaths is given ending before it too (2026-10-11 outside review, R3) pre_total = g.loc[2020:2021, "total"] fr["total_before_break"] = {"years": [2020, 2021], "mean": r(pre_total.mean()), "change": count_change(pre_total.sum(), early["total"].sum(), exposure_now=len(pre_total), exposure_before=len(early))} claim("fars.total_before_break.change", "down") A.printed("fars.total_before_break.mean", f"{fr['total_before_break']['mean']:.0f}") # the headline changes' model ranges, printed in the FAQ for c in ("pedestrian", "vehicle occupant", "two wheels"): ch = fr["change"][c] for end in ("lo_pct", "hi_pct"): A.printed(f"fars.change.{c}.{end}", f"{abs(ch[end]):.1f}%") A.printed(f"fars.two_wheel_share.{roads.FARS_FIRST}", f"{fr['two_wheel_share'][str(roads.FARS_FIRST)]:.0f}%") A.printed(f"fars.two_wheel_share.{roads.FARS_LAST}", f"{fr['two_wheel_share'][str(roads.FARS_LAST)]:.0f}%") A.long_view(roads.FARS_FIRST, roads.FARS_LAST, "NHTSA FARS, New York City's five counties, deaths by road user", why_shorter=None) # the three sources' counts for one year, for the FAQ on why they differ k, _ = roads.collision_deaths() ny, first5 = nypd_years(k) A.results["sources_2024"] = {"fars": int(g.loc[2024, "total"]), "nypd": int(ny.loc[2024, "total"]), "dot": dot["deaths_2024"]} for key in ("fars", "nypd"): A.printed(f"sources_2024.{key}", fmt="int") # ---- the NYPD's files: e-bikes, e-scooters and mopeds apart # every dead person's vehicle record is read and the joins are checked (collision_deaths raises on a # repeated id, a many-to-many join or a vehicle from another crash); a check that fails or cannot run # stops the build vchecks = roads.check_own_vehicles(k) + roads.check_light_extract() A.results["vehicle_checks"] = [x.as_dict() for x in vchecks] if not all(x.passed is True for x in vchecks): raise AssertionError("a vehicle-record check failed or could not run: " + "; ".join(f"{x.name}: {x.detail}" for x in vchecks if x.passed is not True)) A.results["nypd"] = { "file_last_day": k["date"].max().date().isoformat(), "years": {str(y): {c: int(ny.loc[y, c]) for c in roads.GROUPS_ORDER + ["total"]} for y in ny.index}, "micro_before_2020_max": int(ny.loc[FIRST_NYPD:2019, "e-micromobility"].max()), "micro_peak": {"year": int(ny["e-micromobility"].idxmax()), "n": int(ny["e-micromobility"].max())}, "micro_last": int(ny.loc[LAST_YEAR, "e-micromobility"]), "bicycle_2019": int(ny.loc[2019, "bicycle"]), "bicycle_last": int(ny.loc[LAST_YEAR, "bicycle"]), "first5": {str(y): {c: int(first5.loc[y, c]) for c in roads.GROUPS_ORDER + ["total"]} for y in first5.index}, } nr = A.results["nypd"] if nr["micro_peak"]["year"] != 2023: raise AssertionError("the prose says e-bike, e-scooter and moped deaths peaked in 2023") kk = k[k["year"].between(FIRST_NYPD, LAST_YEAR)].assign(group=lambda d: d["road_user"].map(roads.GROUP)) off = kk[kk["person_type"].ne("Pedestrian") & kk["vehicle_id"].notna()] fetched, unknown = off[off["own_source"].eq("fetched by id")], off[off["own"].eq("unknown")] nr["own_vehicles"] = { "years": [FIRST_NYPD, LAST_YEAR], "off_foot_with_vehicle_id": int(len(off)), "in_light_extract": int(off["own_source"].eq("light-vehicle extract").sum()), "fetched_by_id": int(len(fetched)), "fetched_motor_vehicle": int(fetched["own"].eq("motor vehicle").sum()), "type_unknown": int(len(unknown)), "type_unknown_by_year": {str(y): int(n) for y, n in unknown.groupby("year").size().items()}, "type_unknown_2013_2016": int(unknown["year"].le(2016).sum()), "without_vehicle_id": int((kk["person_type"].ne("Pedestrian") & kk["vehicle_id"].isna()).sum()), } if nr["own_vehicles"]["fetched_motor_vehicle"] + nr["own_vehicles"]["type_unknown"] != nr["own_vehicles"]["fetched_by_id"]: raise AssertionError("the prose says every fetched record is a motor vehicle or gives no type") for key in ("fetched_by_id", "type_unknown", "type_unknown_2013_2016", "without_vehicle_id"): A.printed(f"nypd.own_vehicles.{key}", fmt="int") # when the police records began to name e-devices: the labels of every light vehicle by year (every # crash, not only fatal ones), and the month e-bike and e-scooter labels arrive in numbers lab = roads.light_labels("Y").loc[FIRST_NYPD:LAST_YEAR] e_mon = roads.light_labels("M")[["e-bike", "e-scooter"]].sum(axis=1) nr["labels"] = { "by_year": {str(y): {c: int(lab.loc[y, c]) for c in roads.LIGHT} for y in lab.index}, "e_labels_by_month_2019": {m: int(n) for m, n in e_mon.items() if m.startswith("2019")}, "e_labels_in_numbers_from": str(e_mon[e_mon.ge(50)].index.min()), "e_labels_2016_2018_max": int((lab["e-bike"] + lab["e-scooter"]).loc[2016:2018].max()), "bicycles_2016_2018_min": int(lab["bicycle"].loc[2016:2018].min()), "e_bikes_2020": int(lab.loc[2020, "e-bike"]), } if nr["labels"]["e_labels_in_numbers_from"] != "2019-08": raise AssertionError("the prose says e-bike and e-scooter labels are in use only from August 2019") A.printed("nypd.labels.e_labels_2016_2018_max", fmt="int") A.printed("nypd.labels.bicycles_2016_2018_min", fmt="int") A.printed("nypd.labels.e_bikes_2020", fmt="int") # the federal record's own coding of mopeds, motorized bicycles (e-bikes) and motor scooters, to 2021 # (from 2022 an e-bike is no longer a motor vehicle there): a check on the police labels fm = roads.fars_motored_cycles(FIRST_NYPD, 2021) ml = (fm[fm["body_typ"].isin(roads.FARS_MOPED_LIKE)].groupby("year")["deaths"].sum() .reindex(range(FIRST_NYPD, 2022), fill_value=0)) A.results["fars_moped_check"] = { "body_types": list(roads.FARS_MOPED_LIKE), "years": {str(y): {"fars_moped_like": int(ml[y]), "nypd_e_micromobility": int(ny.loc[y, "e-micromobility"])} for y in ml.index}, "fars_2013_2019_mean": r(ml.loc[:2019].mean()), "fars_2013_2019_max": int(ml.loc[:2019].max()), "fars_2021": int(ml[2021]), "nypd_2021": int(ny.loc[2021, "e-micromobility"]), } fmc = A.results["fars_moped_check"] if fmc["fars_2013_2019_max"] > nr["micro_before_2020_max"] or fmc["fars_2021"] < 4 * fmc["fars_2013_2019_max"]: raise AssertionError("the prose says the federal record also counts few such deaths before 2020 and many more by 2021") A.printed("fars_moped_check.fars_2013_2019_mean", f"{fmc['fars_2013_2019_mean']:.0f}") A.printed("fars_moped_check.fars_2021", fmt="int") # the step in 2020 itself, by vehicle: mostly motorcycles; e-devices grew most from 2021 to 2023 nr["step_2020"] = {c: int(ny.loc[2020, c] - ny.loc[2019, c]) for c in TWO_WHEELS_NYPD} em = ny["e-micromobility"] if max(nr["step_2020"], key=nr["step_2020"].get) != "motorcycle" or not em.loc[2020] < em.loc[2021] < em.loc[2022] < em.loc[2023]: raise AssertionError("the prose says the 2020 step was mostly motorcycles and e-device deaths grew from 2021 to 2023") A.printed("nypd.years.2019.motorcycle", fmt="int") A.printed("nypd.years.2020.motorcycle", fmt="int") # the pedal-bicycle base year is the series' high, and the page says so if int(ny["bicycle"].idxmax()) != 2019: raise AssertionError("the prose says 2019 was the highest year for pedal-bicycle deaths since 2013") if int(em.loc[2022]) != nr["micro_last"]: raise AssertionError("the prose says e-device deaths in the last year equal 2022's") nr["bicycle_since_2019"] = count_change(nr["bicycle_last"], nr["bicycle_2019"]) claim("nypd.bicycle_since_2019", "down") nr["micro_since_2019"] = count_change(nr["micro_last"], int(ny.loc[2019, "e-micromobility"])) claim("nypd.micro_since_2019", "up") # a fall from a peak picked out of the same series: given with its range, not called a fall in the # story, since it does not survive the robustness check below nr["micro_since_peak"] = count_change(nr["micro_last"], nr["micro_peak"]["n"]) tw = ny[TWO_WHEELS_NYPD].sum(axis=1) before = tw.loc[FIRST_NYPD:2019] nr["two_wheels"] = {"years": {str(y): int(v) for y, v in tw.items()}, "peak_year": int(tw.idxmax()), "peak": int(tw.max()), "last": int(tw.loc[LAST_YEAR]), "mean_2013_2019": r(before.mean(), 0)} if nr["two_wheels"]["peak_year"] != 2023: raise AssertionError("the prose says deaths on two wheels peaked in 2023") nr["two_wheels_since_peak"] = count_change(tw.loc[LAST_YEAR], tw.max()) claim("nypd.two_wheels_since_peak", "down") # the last year against the average of 2013-2019: a year's rate against a seven-year rate nr["two_wheels_vs_2013_2019"] = count_change(tw.loc[LAST_YEAR], before.sum(), exposure_before=len(before)) claim("nypd.two_wheels_vs_2013_2019", "up") for key in ("two_wheels.peak", "two_wheels.last", "two_wheels.mean_2013_2019"): A.printed(f"nypd.{key}", fmt="int") for key in ("micro_before_2020_max", "micro_peak.n", "micro_last", "bicycle_2019", "bicycle_last"): A.printed(f"nypd.{key}", fmt="int") jm26 = first5.loc[2026] nr["first5_2026_two_wheels"] = int(jm26[TWO_WHEELS_NYPD].sum()) nr["first5_2026_total"] = int(jm26["total"]) # this year so far in the police records: January to May on two wheels against the same months of the # two years before; given as counts, not called a change (it does not survive the extra swing) f5 = first5[TWO_WHEELS_NYPD].sum(axis=1) nr["first5_two_wheels"] = {str(y): int(f5.loc[y]) for y in (LAST_YEAR - 1, LAST_YEAR, LAST_YEAR + 1)} nr["first5_two_wheels_vs_last_year"] = count_change(f5.loc[LAST_YEAR + 1], f5.loc[LAST_YEAR]) for y in (LAST_YEAR - 1, LAST_YEAR, LAST_YEAR + 1): A.printed(f"nypd.first5_two_wheels.{y}", fmt="int") # how each death on two wheels was classed: its own vehicle's free-text type, the person type, or the # crash's vehicles (a fallback when the person has no vehicle id) on2 = kk[kk["group"].isin(TWO_WHEELS_NYPD)] route = on2["how"].map(lambda h: "own vehicle" if h == "own vehicle type" else "crash fallback" if h.startswith("no vehicle_id; crash-level") else "person type") by_route = on2.groupby([on2["year"], route]).size().unstack(fill_value=0) own_type = on2["own_vehicle_type"].fillna("").str.strip().str.lower() deaths_per_crash = k["collision_id"].value_counts() nr["classification"] = { "two_wheel_deaths": int(len(on2)), "by_route": {c: int(by_route[c].sum()) for c in by_route.columns}, "by_route_year": {str(y): {c: int(by_route.loc[y, c]) for c in by_route.columns} for y in by_route.index}, "standing_scooter": int(own_type.str.contains(r"standing s").sum()), "minibike_as_moped": int(own_type.str.contains(r"mini ?bike|minicycle").sum()), "pedicab_as_bicycle": int(own_type.str.contains(r"pedicab").sum()), # deaths that share a crash with another death, which a Poisson model treats as independent "sharing_a_crash": int(on2["collision_id"].map(deaths_per_crash).gt(1).sum()), } if nr["classification"]["by_route"].get("crash fallback", 0): raise AssertionError("the prose says no death on two wheels was classed by the crash's vehicles") A.printed("nypd.classification.sharing_a_crash", fmt="int") A.printed("nypd.classification.two_wheel_deaths", fmt="int") sens = on2[route.ne("crash fallback")].groupby(["year", "group"]).size().unstack(fill_value=0) nr["micro_without_fallback"] = {str(y): int(sens.loc[y].get("e-micromobility", 0)) for y in sens.index} if int(sens["e-micromobility"].idxmax()) != 2023: raise AssertionError("without the crash-level fallback the e-bike, e-scooter and moped peak is not 2023") # crash-level classes against own-vehicle classes, for e-bike deaths ebike = kk[kk["road_user"].eq("e-bike")] nr["ebike_crash_level_disagrees"] = {"n": int((ebike["road_user_crash"] != "e-bike").sum()), "of": int(len(ebike))} # the rise since 2020 in the police records, by vehicle: 2013-2019 against 2020-2025, a year on # average. Net accounting, as the records label the vehicles: e-device labels are in use only from # August 2019, so bicycles and e-devices together are also given, a sum the labels cannot move a7, b6 = ny.loc[FIRST_NYPD:2019], ny.loc[2020:LAST_YEAR] parts = {c: (a7[c], b6[c]) for c in TWO_WHEELS_NYPD} parts["two wheels"] = (a7[TWO_WHEELS_NYPD].sum(axis=1), b6[TWO_WHEELS_NYPD].sum(axis=1)) parts["bicycle and e-devices"] = (a7[["bicycle", "e-micromobility"]].sum(axis=1), b6[["bicycle", "e-micromobility"]].sum(axis=1)) rise = {"before": [FIRST_NYPD, 2019], "after": [2020, LAST_YEAR]} rise.update({c: r(y.mean() - x.mean()) for c, (x, y) in parts.items()}) rise["rounded"] = {c: r(rise[c], 0) for c in parts} rise["change"] = {c: count_change(y.sum(), x.sum(), exposure_now=len(y), exposure_before=len(x)) for c, (x, y) in parts.items()} nr["rise_since_2020"] = rise if not rise["e-micromobility"] > rise["two wheels"] > 0 > rise["bicycle"] or rise["motorcycle"] <= 0: raise AssertionError("the prose says the labelled e-device rise slightly exceeds the net rise, motorcycles averaged " "more and pedal bicycles fewer") for c in ("two wheels", "e-micromobility", "bicycle and e-devices", "bicycle"): claim(f"nypd.rise_since_2020.change.{c}", "up" if rise[c] > 0 else "down") for c in ("two wheels", "e-micromobility", "motorcycle", "bicycle", "bicycle and e-devices"): A.printed(f"nypd.rise_since_2020.rounded.{c}", f"{abs(rise['rounded'][c]):.0f}") for c in ("two wheels", "e-micromobility", "motorcycle", "bicycle"): A.printed(f"nypd.rise_since_2020.{c}", f"{abs(rise[c]):.1f}") A.printed("nypd.ebike_crash_level_disagrees.n", fmt="int") A.printed("nypd.ebike_crash_level_disagrees.of", fmt="int") # ---- robustness: every named change again, with the variance scaled by the extra year-to-year swing # its series shows (sparkylib.stats.dispersion, at least 1); a change named in the story must survive def phi(values) -> float: return round(max(1.0, dispersion(list(values))), 2) phis = {f"fars:{c}": phi(g[c]) for c in ("pedestrian", "vehicle occupant", "two wheels", "total")} phis.update({f"nypd:{c}": phi(ny[c]) for c in TWO_WHEELS_NYPD}) phis["nypd:two wheels"] = phi(tw) phis["nypd:bicycle and e-devices"] = phi(ny["bicycle"] + ny["e-micromobility"]) again = {f"fars.change.{c}": (fr["change"][c], 1, 1, f"fars:{c}") for c in ("pedestrian", "vehicle occupant", "two wheels", "total")} again["fars.two_wheels_step_vs_plateau"] = (fr["two_wheels_step_vs_plateau"], len(step), len(plateau), "fars:two wheels") again["fars.two_wheels_before_break.change"] = (fr["two_wheels_before_break"]["change"], len(pre), len(early), "fars:two wheels") again["fars.total_before_break.change"] = (fr["total_before_break"]["change"], len(pre_total), len(early), "fars:total") for key, c in (("bicycle_since_2019", "bicycle"), ("micro_since_2019", "e-micromobility"), ("micro_since_peak", "e-micromobility"), ("two_wheels_since_peak", "two wheels"), ("first5_two_wheels_vs_last_year", "two wheels")): again[f"nypd.{key}"] = (nr[key], 1, 1, f"nypd:{c}") again["nypd.two_wheels_vs_2013_2019"] = (nr["two_wheels_vs_2013_2019"], 1, len(before), "nypd:two wheels") for c in parts: again[f"nypd.rise_since_2020.change.{c}"] = (rise["change"][c], len(b6), len(a7), f"nypd:{c}") # keyed by the change's results path with "_" for ".", the path itself kept inside robust = {path.replace(".", "_"): {**count_change(node["now"], node["before"], exposure_now=en, exposure_before=eb, dispersion=phis[s]), "series": s, "path": path} for path, (node, en, eb, s) in again.items()} A.results["robustness"] = {"dispersion": phis, "max_dispersion": max(phis.values()), "changes": robust} unchecked = [p for p in claimed if p.replace(".", "_") not in robust] lost = [p for p, d in claimed.items() if p.replace(".", "_") in robust and robust[p.replace(".", "_")]["direction"] != d] if unchecked or lost: raise AssertionError(f"named changes not rechecked {unchecked} or not surviving the extra year-to-year swing {lost}") for key in ("nypd_micro_since_peak", "nypd_rise_since_2020_change_motorcycle", "nypd_first5_two_wheels_vs_last_year"): if robust[key]["direction"] != "unclear": raise AssertionError(f"{key}: the prose says it is not sure once the extra swing is allowed for") A.printed("robustness.max_dispersion", f"{A.results['robustness']['max_dispersion']:.1f}") # the model ranges the questions print for path in ("nypd.two_wheels_since_peak", "nypd.two_wheels_vs_2013_2019", "fars.two_wheels_before_break.change", "fars.change.total", "fars.total_before_break.change", "nypd.micro_since_peak", "nypd.bicycle_since_2019", "dot_changes.h1_micro", "nypd.rise_since_2020.change.two wheels", "nypd.rise_since_2020.change.motorcycle", "nypd.first5_two_wheels_vs_last_year", "robustness.changes.nypd_micro_since_peak", "robustness.changes.nypd_rise_since_2020_change_motorcycle", "robustness.changes.nypd_first5_two_wheels_vs_last_year"): node = A.results for part in path.split("."): node = node[part] for end in ("lo_pct", "hi_pct"): A.printed(f"{path}.{end}", f"{abs(node[end]):.1f}%") # ---- reproduction: the NYPD's own-vehicle classes against DOT's first-half counts for year, key, tol in ((2023, "h1_micro_2024", 1), (2024, "h1_micro_2024", 0), (2025, "h1_micro_2025", 0)): A.reproduce(f"E-bike, e-scooter and moped riders killed, 1 January to 29 June {year}", published=dot[key], published_moe90=None, reproduced=roads.first_half_micro(k, year), kind="count", tolerance=tol) # the federal extract against the crash files' own totals, every year A.results["fars_check"] = [x.as_dict() for x in roads.check_fars()] if not all(x["passed"] for x in A.results["fars_check"]): raise AssertionError("the federal extract does not match the crash files' totals") A.results["reproduction_source"] = "NYC DOT, first-half release of 1 July 2026 (its counts for 2023, 2024 and 2025; 23 in both 2023 and 2024)" A.results["reproduction_note"] = ("Administrative counts: the NYPD person file, each death classed by the dead person's own " "vehicle, against DOT's count for the same days; exact for 2024 and 2025, within one " "for 2023.") # ---- data.csv rows = [] for y in g.index: for c in roads.FARS_GROUPS: rows.append({"source": "NHTSA FARS", "year": int(y), "group": c, "deaths": int(g.loc[y, c])}) for y in ny.index: for c in roads.GROUPS_ORDER: rows.append({"source": "NYPD collisions, whole year", "year": int(y), "group": c, "deaths": int(ny.loc[y, c])}) for y in first5.index: for c in roads.GROUPS_ORDER: rows.append({"source": "NYPD collisions, January to May", "year": int(y), "group": c, "deaths": int(first5.loc[y, c])}) data = pd.DataFrame(rows) data["group"] = data["group"].replace({"e-micromobility": "e-bike, e-scooter or moped"}) # ---- charts fy = [str(y) for y in g.index] A.chart("thirty_years", lines( title="Deaths on foot fell by about half and in vehicles by more; deaths on two wheels rose", subtitle=f"Traffic deaths in New York City by road user, {roads.FARS_FIRST} to {roads.FARS_LAST}", x=[[y, y] for y in fy], series=[[c, LONG_NAMES[c], LONG_SHORT[c]] for c in LONG_GROUPS], cells={"": {c: {y: int(g.loc[int(y), c]) for y in fy} for c in LONG_GROUPS}}, format="num", axis="Deaths a year", x_ticks=[y for y in fy if int(y) % 5 == 0 or y == fy[-1]], note="The federal record counts a death within 30 days of a crash involving a motor vehicle on a public road. " "Two wheels are motorcycles, mopeds, motor scooters and bicycles together, e-bikes included: the federal " "record counted e-bikes with mopeds to 2021 and with bicycles from 2022, when it stopped counting crashes " "involving only e-bikes. People on stand-up scooters, skates and wheelchairs, which the record counts as " "personal conveyances, are outside the line in every year, with a few deaths of unknown type. 2024 is the " "record's provisional Annual Report File, replaced by a final file about a year later.", provenance="NHTSA Fatality Analysis Reporting System, 1995-2024")) ty = [str(y) for y in ny.index] A.chart("two_wheels", lines( title="Deaths recorded on e-bikes, e-scooters and mopeds rose from almost none to a peak in 2023", subtitle=f"People killed on two wheels, by their own vehicle, {FIRST_NYPD} to {LAST_YEAR}", x=[[y, y] for y in ty], series=[["e-micromobility", "E-bikes, e-scooters and mopeds", "E-bikes, scooters, mopeds"], ["motorcycle", "Motorcycles", "Motorcycles"], ["bicycle", "Pedal bicycles", "Bicycles"]], cells={"": {c: {y: int(ny.loc[int(y), c]) for y in ty} for c in TWO_WHEELS_NYPD}}, format="num", axis="People killed a year", x_ticks=[y for y in ty if int(y) % 2 == 1 or y == ty[-1]], note="The NYPD's person file, each death classed by the dead person's own vehicle as its record labels it; it " "reproduces DOT's first-half counts of e-bike, e-scooter and moped riders killed in 2024 and 2025 exactly. " "The police records name e-bikes and e-scooters routinely only from August 2019, so before then some " "e-bike riders may be counted with pedal bicycles. Deaths as the police records hold them; the NYPD doesn't " "publish when a later death is added.", provenance="NYPD Motor Vehicle Collisions (persons, vehicles), NYC Open Data")) A.write(data=data) print(f"{SLUG}: wrote {A.out}") return 0 if __name__ == "__main__": raise SystemExit(main())