Stage 2: Collision Risk Model

Analysis of Stage 2 collision risk model outputs, feature behaviour, metrics, ranking patterns, and diagnostic comparisons.

Last updated: 2026-07-04 · Full GB output rebuild: 2026-07-04.

1 Overview

Stage 2 fits a Poisson collision model on all OS Open Roads network data links × AADF years using estimated AADT as the exposure offset. This is the main exposure-adjusted collision risk model. Two complementary models are used:

• Poisson GLM (statsmodels) — interpretable coefficients, fast inference
• XGBoost Poisson — captures non-linear interactions, higher predictive power

Both models are trained on link × year data then pooled to a single stable risk score per link. The pooled output (risk_scores.parquet) has one row per link with no year dimension.

1.1 Current full-GB collision retrain

The current full-GB Stage 2 collision retrain uses the refreshed GB Open Roads and context feature tables. The run built 39,412,990 link-year rows for 2015-2024 and scored all 3,941,299 road links. The model included GB population, deprivation, and rural/urban context fields, and excluded the retired English-only IMD fields.

The XGBoost Poisson model achieved pseudo-R² = 0.360 with full-zero training over the full 39,412,990 link-year population: 945,373 positive road-link × year rows and 38,467,617 zero rows. The production XGBoost fit uses tree_method="hist" and max_bin=128; sampled-zero XGBoost is available only as an explicit memory fallback. The GLM baseline converged with pseudo-R² = 0.505 on its 3,781,492-row 1:3 zero-collision sampled surface. These values are not directly comparable because the GLM and XGBoost evaluation surfaces differ.

Outputs were written to data/models/risk_scores.parquet, data/models/collision_glm.pkl, data/models/collision_xgb.json, and data/models/collision_metrics.json.

from pathlib import Path
import json
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
from IPython.display import HTML, display

from road_risk.config import _ROOT as ROOT, cfg
RANDOM_STATE = 42
GB_BBOX = cfg["study_area"]["bbox_wgs84"]
GB_MAP_EXTENT = {
    "min_lon": min(GB_BBOX["min_lon"], -8.7),
    "max_lon": max(GB_BBOX["max_lon"], 2.1),
    "min_lat": min(GB_BBOX["min_lat"], 49.8),
    "max_lat": max(GB_BBOX["max_lat"], 60.9),
}

try:
    import geopandas as gpd
    HAS_GPD = True
except Exception:
    HAS_GPD = False

risk = pd.read_parquet(ROOT / "data/models/risk_scores.parquet")

with open(ROOT / "data/models/collision_metrics.json") as f:
    metrics = json.load(f)

XGB_CURRENT_PR2 = metrics["xgb"]["pseudo_r2"]

if HAS_GPD:
    or_path = ROOT / "data/processed/shapefiles/openroads.parquet"
    openroads = gpd.read_parquet(or_path) if or_path.exists() else None
else:
    openroads = None


def set_gb_wgs84_extent(ax):
    """Use a stable GB frame and latitude-aware WGS84 aspect."""
    ax.set_xlim(GB_MAP_EXTENT["min_lon"], GB_MAP_EXTENT["max_lon"])
    ax.set_ylim(GB_MAP_EXTENT["min_lat"], GB_MAP_EXTENT["max_lat"])
    mid_lat = (GB_MAP_EXTENT["min_lat"] + GB_MAP_EXTENT["max_lat"]) / 2
    ax.set_aspect(1 / np.cos(np.deg2rad(mid_lat)), adjustable="box")
    ax.set_anchor("C")


MAP_FIGSIZE = (7.4, 9.6)
MAP_DPI = 170
MAP_LEGEND_KWDS = {
    "shrink": 0.58,
    "fraction": 0.035,
    "pad": 0.015,
}
RISK_CMAP = mcolors.LinearSegmentedColormap.from_list(
    "orr_risk_percentile",
    ["#2166ac", "#67a9cf", "#fddbc7", "#ef8a62", "#b2182b"],
)


def clipped_upper(series, q=0.99, fallback=None):
    value = series.dropna().quantile(q)
    if pd.isna(value) or value <= 0:
        return fallback
    return float(value)


def show_table(df, caption=None, formats=None, index=True):
    """Render compact reader-facing tables instead of plain console strings."""
    table = df.copy()
    if formats:
        for col, fmt in formats.items():
            if col in table.columns:
                table[col] = table[col].map(
                    lambda x, fmt=fmt: "" if pd.isna(x) else fmt.format(x)
                )
    html = table.to_html(
        classes="orr-table",
        border=0,
        escape=False,
        index=index,
    )
    if caption:
        html = html.replace(">", f"><caption>{caption}</caption>", 1)
    display(HTML(html))


def show_metric_table(rows, caption=None):
    show_table(pd.DataFrame(rows, columns=["Metric", "Value"]), caption=caption, index=False)

cls_order = ["Motorway", "A Road", "B Road", "Classified Unnumbered",
             "Not Classified", "Unclassified", "Unknown"]
cls_in_data = [c for c in cls_order if c in risk["road_classification"].unique()]

show_metric_table(
    [
        ("Links scored", f"{len(risk):,}"),
        ("Links with collisions", f"{(risk['collision_count'] > 0).sum():,} ({(risk['collision_count'] > 0).mean():.1%})"),
        ("Total collisions", f"{risk['collision_count'].sum():,}"),
        ("Total fatals", f"{risk['fatal_count'].sum():,}"),
        ("GLM pseudo-R²", f"{metrics['glm']['pseudo_r2']:.3f}"),
        ("XGB pseudo-R²", f"{XGB_CURRENT_PR2:.3f} (current clean GB run)"),
        ("GLM converged", metrics["glm"]["converged"]),
    ],
    caption="Current Stage 2 model outputs",
)

---

2 Model performance

2.1 GLM summary

g = metrics["glm"]
show_metric_table(
    [
        ("Training rows (downsampled for GLM)", f"{g['n_obs']:,}"),
        ("Full dataset rows", f"{g['n_full']:,}"),
        ("Positive link-years", f"{g['n_pos']:,} ({g['n_pos']/g['n_full']:.2%})"),
        ("Pseudo-R² (1 - D/D₀)", f"{g['pseudo_r2']:.3f}"),
        ("AIC", f"{g['aic']:,.0f}"),
        ("Features", f"{len(g['features'])}"),
    ],
    caption="GLM training surface",
)
GLM_FEATURE_DETAILS = {
    "road_class_ord": (
        "Ordinal road-class code, with motorways highest and unknown lowest.",
        "OS Open Roads `road_classification`; encoded in `road_risk.model.constants`.",
    ),
    "form_of_way_ord": (
        "Ordinal carriageway/form code distinguishing duals, slips, roundabouts, and single carriageways.",
        "OS Open Roads `form_of_way`; encoded in `road_risk.model.constants`.",
    ),
    "is_motorway": (
        "Binary flag for motorway links.",
        "Derived from OS Open Roads `road_classification`.",
    ),
    "is_a_road": (
        "Binary flag for A-road links.",
        "Derived from OS Open Roads `road_classification`.",
    ),
    "is_slip_road": (
        "Binary flag for slip roads.",
        "Derived from OS Open Roads `form_of_way`.",
    ),
    "is_roundabout": (
        "Binary flag for roundabout links.",
        "Derived from OS Open Roads `form_of_way`.",
    ),
    "is_dual": (
        "Binary flag for dual or collapsed-dual carriageway links.",
        "Derived from OS Open Roads `form_of_way`.",
    ),
    "is_trunk": (
        "Binary flag for trunk-road links.",
        "OS Open Roads `is_trunk`.",
    ),
    "is_primary": (
        "Binary flag for primary-route links.",
        "OS Open Roads `is_primary`.",
    ),
    "log_link_length": (
        "Natural log of link length in kilometres; separate from the exposure offset so length can depart from a strict linear effect.",
        "OS Open Roads `link_length_km`.",
    ),
    "is_covid": (
        "Binary flag for Covid-period years.",
        "AADF/model year; 2020 and 2021 in `road_risk.model.constants`.",
    ),
    "year_norm": (
        "Training year scaled from 0 to 1 across the model period.",
        "AADF/model year in the Stage 2 link-year table.",
    ),
    "hgv_proportion": (
        "Heavy-goods-vehicle share of traffic.",
        "`data/features/road_traffic_features.parquet`, derived from AADF traffic composition.",
    ),
    "degree_mean": (
        "Mean graph degree of the link's start and end nodes; a junction-complexity proxy.",
        "`data/features/network_features.parquet`, built from the OS Open Roads graph.",
    ),
    "betweenness": (
        "Approximate betweenness centrality averaged across the link endpoints.",
        "`data/features/network_features.parquet`, built with NetworkX on OS Open Roads.",
    ),
    "betweenness_relative": (
        "Log centrality relative to the mean for the same road class.",
        "`data/features/network_features.parquet`, derived from `betweenness` and road class.",
    ),
    "dist_to_major_km": (
        "Graph distance to the nearest motorway or A-road node.",
        "`data/features/network_features.parquet`, built from the OS Open Roads graph.",
    ),
    "pop_density_per_km2": (
        "Population density at the road-link centroid.",
        "GB OA population-density context joined to road-link centroids.",
    ),
    "overall_decile_within_country": (
        "Overall deprivation decile within each nation; 1 is most deprived, 10 least deprived.",
        "England IoD 2025, Wales WIMD 2019, Scotland SIMD 2020v2; joined by link centroid.",
    ),
    "income_decile_within_country": (
        "Income deprivation decile within each nation; 1 is most deprived, 10 least deprived.",
        "England/Wales/Scotland deprivation context where available; joined by link centroid.",
    ),
    "employment_decile_within_country": (
        "Employment deprivation decile within each nation; 1 is most deprived, 10 least deprived.",
        "England/Wales/Scotland deprivation context where available; joined by link centroid.",
    ),
    "deprivation_country_england": (
        "Country indicator for England deprivation assignments.",
        "GB deprivation assignment provenance.",
    ),
    "deprivation_country_wales": (
        "Country indicator for Wales deprivation assignments.",
        "GB deprivation assignment provenance.",
    ),
    "deprivation_country_scotland": (
        "Country indicator for Scotland deprivation assignments.",
        "GB deprivation assignment provenance.",
    ),
}


def describe_glm_feature(feature):
    if feature.endswith("_missing"):
        base = feature.removesuffix("_missing")
        treatment = "Missingness indicator: 1 where the raw feature was unavailable."
    elif feature.endswith("_imputed"):
        base = feature.removesuffix("_imputed")
        treatment = "Median-imputed numeric value used by the GLM."
    else:
        base = feature
        treatment = "Core model feature."
    description, source = GLM_FEATURE_DETAILS.get(
        base,
        ("Model input from the Stage 2 training table.", "Stage 2 modelling dataset."),
    )
    return {
        "Feature": feature,
        "Description": description,
        "Source / construction": source,
        "GLM treatment": treatment,
    }


show_table(
    pd.DataFrame([describe_glm_feature(feature) for feature in g["features"]]),
    caption="GLM feature list with source and treatment",
    index=False,
)

2.2 GLM coefficients

glm_path = ROOT / "data/models/collision_glm.pkl"
if glm_path.exists():
    try:
        import statsmodels.api as sm
        glm_result = sm.load(str(glm_path))

        coef = pd.DataFrame({
            "coef":    glm_result.params,
            "ci_low":  glm_result.conf_int()[0],
            "ci_high": glm_result.conf_int()[1],
            "pvalue":  glm_result.pvalues,
        }).drop(index="const", errors="ignore").sort_values("coef")

        bar_colors = ["#e63946" if p < 0.05 else "#adb5bd" for p in coef["pvalue"]]

        ci_extent = np.nanmax(np.abs(coef[["ci_low", "ci_high"]].to_numpy()))
        coef_extent = np.nanmax(np.abs(coef["coef"].to_numpy()))
        x_extent = min(max(ci_extent, coef_extent) * 1.12, 3.0)

        fig, ax = plt.subplots(figsize=(11, max(6, len(coef) * 0.48)))
        ax.barh(coef.index, coef["coef"], color=bar_colors, alpha=0.85)
        xerr_low = (coef["coef"] - coef["ci_low"]).clip(lower=0)
        xerr_high = (coef["ci_high"] - coef["coef"]).clip(lower=0)
        ax.errorbar(
            coef["coef"], range(len(coef)),
            xerr=[xerr_low, xerr_high],
            fmt="none", color="black", linewidth=0.8, capsize=3,
        )
        ax.axvline(0, color="black", linewidth=0.8, linestyle="--")
        ax.set_xlim(-x_extent, x_extent)
        ax.grid(axis="x", color="#e5e7eb", linewidth=0.8)
        ax.set_xlabel("Coefficient (log-rate scale)\nRed = significant p<0.05")
        ax.set_title(
            "Poisson GLM coefficients\n"
            "(positive = higher collision rate per vehicle-km)"
        )
        plt.tight_layout()
        plt.show()

        sig = coef[coef["pvalue"] < 0.05].copy()
        sig["IRR"] = np.exp(sig["coef"])
        show_table(
            sig[["coef", "IRR", "pvalue"]],
            caption="Significant coefficients (p<0.05) with incidence rate ratios",
            formats={"coef": "{:.4f}", "IRR": "{:.3f}", "pvalue": "{:.2e}"},
        )

    except Exception as e:
        print(f"Could not load GLM: {e}")
else:
    print("collision_glm.pkl not found — run: python -m road_risk.model --stage collision")

2.3 Deprivation signal

The GB deprivation features are within-country deciles, running 1 (most deprived) to 10 (least deprived). They are contextual area features, not GB-wide absolute deprivation ranks and not road attributes. See deprivation.

deprivation_terms = [
    "overall_decile_within_country",
    "income_decile_within_country",
    "employment_decile_within_country",
    "deprivation_country_england",
    "deprivation_country_wales",
    "deprivation_country_scotland",
]
try:
    rows = []
    for base in deprivation_terms:
        for t in [base, f"{base}_imputed"]:
            if t in glm_result.params.index:
                rows.append((t, glm_result.params[t],
                             np.exp(glm_result.params[t]), glm_result.pvalues[t]))
    if rows:
        deprivation_df = pd.DataFrame(rows, columns=["term", "coef", "IRR", "pvalue"])
        show_table(
            deprivation_df,
            caption="GB deprivation terms in the fitted GLM",
            formats={"coef": "{:.4f}", "IRR": "{:.3f}", "pvalue": "{:.2e}"},
            index=False,
        )
    else:
        print("No GB deprivation terms in the fitted GLM.")
except NameError:
    print("Run the GLM coefficients cell first.")

With ~39.4M link-years, statistical significance is near-guaranteed; the practical effect is modest (each decile step shifts the rate a few percent) and is associational, not causal — LSOA-level context, not a road attribute.

2.4 XGBoost feature importance

xgb_path = ROOT / "data/models/collision_xgb.json"
if xgb_path.exists():
    try:
        from xgboost import XGBRegressor
        xgb = XGBRegressor()
        xgb.load_model(str(xgb_path))

        importance = pd.Series(
            xgb.feature_importances_,
            index=metrics["xgb"]["features"],
        ).sort_values()

        fig, ax = plt.subplots(figsize=(8, max(5, len(importance) * 0.45)))
        ax.barh(importance.index, importance.values, color="#457b9d", alpha=0.85)
        ax.set_xlabel("Feature importance (gain, normalised)")
        ax.set_title("XGBoost Poisson — feature importance")
        plt.tight_layout()
        plt.show()

    except Exception as e:
        print(f"Could not load XGBoost: {e}")
else:
    print("collision_xgb.json not found.")

2.5 GLM vs XGBoost agreement

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

has_both = risk[["predicted_glm", "predicted_xgb"]].gt(0).all(axis=1)
sample = risk[has_both].sample(min(50_000, has_both.sum()), random_state=RANDOM_STATE)

axes[0].scatter(np.log(sample["predicted_glm"]), np.log(sample["predicted_xgb"]),
                s=2, alpha=0.15, color="#1d3557")
mn = min(np.log(sample["predicted_glm"]).min(), np.log(sample["predicted_xgb"]).min())
mx = max(np.log(sample["predicted_glm"]).max(), np.log(sample["predicted_xgb"]).max())
axes[0].plot([mn, mx], [mn, mx], "r--", linewidth=1)
corr = np.corrcoef(np.log(sample["predicted_glm"]), np.log(sample["predicted_xgb"]))[0, 1]
axes[0].text(0.05, 0.92, f"r = {corr:.3f}", transform=axes[0].transAxes, fontsize=10)
axes[0].set_xlabel("log(GLM predicted rate)")
axes[0].set_ylabel("log(XGBoost predicted rate)")
axes[0].set_title("GLM vs XGBoost predictions (50k sample, log scale)")

pct_data = [risk[risk["road_classification"] == c]["risk_percentile"].values
            for c in cls_in_data]
axes[1].boxplot(pct_data, labels=cls_in_data, patch_artist=True,
                boxprops=dict(facecolor="#a8dadc", color="#1d3557"),
                medianprops=dict(color="#e63946", linewidth=2))
axes[1].set_ylabel("Risk percentile")
axes[1].set_title("Risk percentile distribution by road class")
axes[1].tick_params(axis="x", rotation=30)

plt.tight_layout()
plt.show()

---

3 Risk score distribution

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

pos = risk[risk["predicted_glm"] > 0]
axes[0].hist(np.log(pos["predicted_glm"]), bins=60, color="#1d3557", edgecolor="none")
axes[0].set_xlabel("log(predicted collisions/year)")
axes[0].set_ylabel("Link count")
axes[0].set_title("Predicted collision rate — all links (log scale)")

axes[1].hist(risk["collision_count"].clip(0, 20), bins=21,
             range=(-0.5, 20.5), color="#457b9d", edgecolor="white", linewidth=0.3)
axes[1].set_yscale("log")
axes[1].set_xlabel("Observed collisions (pooled years, capped at 20)")
axes[1].set_ylabel("Link count (log scale)")
axes[1].set_title(
    f"Observed collision counts\n"
    f"({(risk['collision_count'] == 0).mean():.1%} of links: zero collisions)"
)

plt.tight_layout()
plt.show()

3.1 Risk percentile by road class

cls_med = (
    risk.groupby("road_classification")["risk_percentile"]
    .median()
    .reindex(cls_in_data)
    .dropna()
    .sort_values(ascending=False)
)
colors = ["#e63946" if rc in ["Motorway", "A Road"] else "#457b9d"
          for rc in cls_med.index]

fig, ax = plt.subplots(figsize=(9, 4))
ax.bar(cls_med.index, cls_med.values, color=colors)
ax.axhline(50, color="black", linestyle="--", linewidth=0.8, alpha=0.5)
ax.set_ylabel("Median risk percentile")
ax.set_title("Median risk percentile by road class  (red = major roads)")
ax.tick_params(axis="x", rotation=30)
plt.tight_layout()
plt.show()

show_table(
    risk.groupby("road_classification")["risk_percentile"]
    .agg(["mean", "median", "count"])
    .reindex(cls_in_data),
    caption="Risk percentile by road class",
    formats={"mean": "{:.1f}", "median": "{:.1f}", "count": "{:,.0f}"},
)

---

4 Residual analysis

Residuals (observed − predicted × n_years) identify where the model systematically under- or over-predicts. Large positive residuals flag links that had more collisions than the model expected given their traffic and road characteristics.

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

axes[0].hist(risk["residual_glm"].clip(-5, 10), bins=80,
             color="#457b9d", edgecolor="none")
axes[0].axvline(0, color="#e63946", linewidth=1.5, linestyle="--")
axes[0].set_yscale("log")
axes[0].set_xlabel("Residual (observed − predicted, pooled years)")
axes[0].set_ylabel("Link count (log scale)")
axes[0].set_title(
    f"GLM residuals\n"
    f"median={risk['residual_glm'].median():.3f}  std={risk['residual_glm'].std():.2f}"
)

pos = risk[risk["predicted_glm"] > 0].sample(50_000, random_state=RANDOM_STATE)
axes[1].scatter(np.log(pos["predicted_glm"]), pos["residual_glm"].clip(-5, 10),
                s=2, alpha=0.15, color="#1d3557")
axes[1].axhline(0, color="#e63946", linewidth=1, linestyle="--")
axes[1].set_xlabel("log(predicted collision rate)")
axes[1].set_ylabel("Residual (clipped −5 to +10)")
axes[1].set_title("Residuals vs predicted rate (overdispersion check)")

plt.tight_layout()
plt.show()

4.1 Residuals by road class

res_by_class = (
    risk.groupby("road_classification")["residual_glm"]
    .agg(["mean", "median", "std", "count"])
    .round(3)
    .reindex(cls_in_data)
    .dropna()
)
show_table(
    res_by_class,
    caption="GLM residuals by road class",
    formats={"mean": "{:.3f}", "median": "{:.3f}", "std": "{:.3f}", "count": "{:,.0f}"},
)

fig, ax = plt.subplots(figsize=(9, 4))
ax.bar(res_by_class.index, res_by_class["median"], color="#457b9d")
ax.axhline(0, color="black", linewidth=0.8, linestyle="--")
ax.set_ylabel("Median residual (observed − predicted)")
ax.set_title("Median GLM residual by road class\n(positive = model under-predicts risk)")
ax.tick_params(axis="x", rotation=30)
plt.tight_layout()
plt.show()

---

5 Geographic risk maps

5.1 Full-network risk percentile

if HAS_GPD and openroads is not None:
    or_risk = openroads[["link_id", "road_classification", "geometry"]].merge(
        risk[["link_id", "risk_percentile", "residual_glm", "collision_count"]],
        on="link_id", how="inner",
    )
    or_risk = or_risk[or_risk.geometry.notna()].copy()
    if or_risk.crs is not None and or_risk.crs.to_epsg() != 4326:
        or_risk = or_risk.to_crs(4326)

    fig, ax = plt.subplots(figsize=MAP_FIGSIZE, dpi=MAP_DPI)
    or_risk.plot(
        column="risk_percentile", ax=ax,
        cmap=RISK_CMAP, vmin=0, vmax=100,
        linewidth=0.12, alpha=0.88, legend=True, rasterized=True,
        legend_kwds={**MAP_LEGEND_KWDS, "label": "Risk percentile"},
    )
    ax.set_title("Risk percentile — all scored GB links")
    set_gb_wgs84_extent(ax)
    ax.set_axis_off()
    plt.tight_layout()
    plt.show()

    major = or_risk[or_risk["road_classification"].isin(["Motorway", "A Road"])].copy()
    major["risk_percentile_within_major_roads"] = major["risk_percentile"].rank(pct=True) * 100
    fig, ax = plt.subplots(figsize=MAP_FIGSIZE, dpi=MAP_DPI)
    major.plot(
        column="risk_percentile_within_major_roads", ax=ax,
        cmap=RISK_CMAP, vmin=0, vmax=100,
        linewidth=0.55, alpha=0.95, legend=True, rasterized=True,
        legend_kwds={**MAP_LEGEND_KWDS, "label": "Percentile within motorways and A roads"},
    )
    ax.set_title("Relative risk percentile — Motorways and A Roads")
    set_gb_wgs84_extent(ax)
    ax.set_axis_off()
    plt.tight_layout()
    plt.show()

5.2 Top 1% risk links

if HAS_GPD and openroads is not None:
    top1 = or_risk[or_risk["risk_percentile"] >= 99].copy()
    rng  = np.random.default_rng(RANDOM_STATE)
    bg   = or_risk.iloc[rng.choice(len(or_risk), size=min(150_000, len(or_risk)), replace=False)]
    collision_cap = clipped_upper(top1["collision_count"], q=0.99, fallback=30)
    residual_low = float(top1["residual_glm"].quantile(0.02))
    residual_high = float(top1["residual_glm"].quantile(0.98))
    residual_limit = max(abs(residual_low), abs(residual_high), 5)
    top1["collision_count_capped"] = top1["collision_count"].clip(upper=collision_cap)
    top1["residual_glm_capped"] = top1["residual_glm"].clip(
        lower=-residual_limit, upper=residual_limit
    )

    fig, ax = plt.subplots(figsize=MAP_FIGSIZE, dpi=MAP_DPI)
    bg.plot(ax=ax, color="#94a3b8", linewidth=0.2, alpha=0.55, zorder=1)
    top1.plot(
        column="collision_count_capped", ax=ax, cmap="OrRd",
        vmin=0, vmax=collision_cap, linewidth=1.15, alpha=1.0, legend=True, zorder=2,
        legend_kwds={
            **MAP_LEGEND_KWDS,
            "label": f"Observed collisions (pooled, capped at {collision_cap:.0f})",
        },
    )
    ax.set_title(f"Top 1% risk links ({len(top1):,}) — observed collisions")
    set_gb_wgs84_extent(ax)
    ax.set_axis_off()
    plt.tight_layout()
    plt.show()

    fig, ax = plt.subplots(figsize=MAP_FIGSIZE, dpi=MAP_DPI)
    bg.plot(ax=ax, color="#94a3b8", linewidth=0.2, alpha=0.55, zorder=1)
    top1.plot(
        column="residual_glm_capped", ax=ax, cmap="RdBu_r",
        norm=mcolors.TwoSlopeNorm(vmin=-residual_limit, vcenter=0, vmax=residual_limit),
        linewidth=1.15, alpha=1.0, legend=True, zorder=2,
        legend_kwds={
            **MAP_LEGEND_KWDS,
            "label": f"GLM residual (capped +/-{residual_limit:.0f})",
        },
    )
    ax.set_title("Top 1% risk links — GLM residual")
    set_gb_wgs84_extent(ax)
    ax.set_axis_off()
    plt.tight_layout()
    plt.show()

    show_table(
        top1["road_classification"].value_counts().rename_axis("road_classification").reset_index(name="links"),
        caption=f"Top 1% risk links by road class ({len(top1):,} total)",
        formats={"links": "{:,.0f}"},
        index=False,
    )

5.3 Excess-risk map (residual > 2)

Links where observed collisions substantially exceeded the model’s prediction given traffic volume and road type — potential candidates for engineering review.

if HAS_GPD and openroads is not None:
    excess = or_risk[or_risk["residual_glm"] > 2].copy()
    excess_cap = clipped_upper(excess["residual_glm"], q=0.99, fallback=40)
    excess["residual_glm_capped"] = excess["residual_glm"].clip(upper=excess_cap)
    if "bg" not in globals():
        rng = np.random.default_rng(RANDOM_STATE)
        bg = or_risk.iloc[rng.choice(len(or_risk), size=min(150_000, len(or_risk)), replace=False)]

    fig, ax = plt.subplots(figsize=MAP_FIGSIZE, dpi=MAP_DPI)
    bg.plot(ax=ax, color="#94a3b8", linewidth=0.2, alpha=0.55, zorder=1)
    excess.plot(column="residual_glm_capped", ax=ax, cmap="PuRd",
                vmin=2, vmax=excess_cap, linewidth=1.25, alpha=1.0, legend=True,
                zorder=2, legend_kwds={
                    **MAP_LEGEND_KWDS,
                    "label": f"Residual (observed - predicted, capped at {excess_cap:.0f})",
                })
    ax.set_title(
        f"Excess-risk links — residual > 2  ({len(excess):,} links)\n"
        "Had substantially more collisions than road type and traffic predict"
    )
    set_gb_wgs84_extent(ax)
    ax.set_axis_off()
    plt.tight_layout()
    plt.show()

---

6 Summary

show_metric_table(
    [
        ("Links scored", f"{len(risk):,}"),
        ("Links with collisions", f"{(risk['collision_count'] > 0).sum():,} ({(risk['collision_count'] > 0).mean():.1%})"),
        ("Poisson GLM pseudo-R²", f"{metrics['glm']['pseudo_r2']:.3f}"),
        ("XGBoost pseudo-R²", f"{XGB_CURRENT_PR2:.3f} (current clean GB run)"),
        ("Top 1% risk links", f"{(risk['risk_percentile'] >= 99).sum():,}"),
        ("Excess-risk links (residual > 2)", f"{(risk['residual_glm'] > 2).sum():,}"),
    ],
    caption="Stage 2 collision model summary",
)
show_table(
    risk[risk["residual_glm"] > 2]["road_classification"]
    .value_counts()
    .rename_axis("road_classification")
    .reset_index(name="links"),
    caption="Excess-risk links by road class",
    formats={"links": "{:,.0f}"},
    index=False,
)

Key observations:

• XGBoost’s current full-GB run has pseudo-R² 0.360, with temporal features included and full-zero training, while the GLM sits at pseudo-R² 0.505 on its own 1:3 zero-collision sampled surface. Earlier repo docs cited an XGBoost pseudo-R² around 0.86, but that figure came from a pre-fix evaluation surface that was later superseded after a feature-table leakage diagnosis. The two current pseudo-R² values are also still computed on different row subsets and with different null models (GLM: in-sample on a downsampled training set; XGBoost: out-of-sample on the full zero-heavy test set), so the gap should not be read as pure modelling lift. These metrics are from the current full-GB run after the counted-only Stage 1a filter, the Stage 2 post-event provenance guard, OSM feature enrichment, deprivation features, terrain grade, and the GLM optional-feature imputation refactor. See feature engineering for the full caveat.
• The GLM provides interpretable incidence rate ratios per feature; XGBoost predictions complement these with higher predictive power for the risk map.
• Residuals show where traffic volume and road geometry alone under-predict observed collisions — likely reflecting local factors not in the feature set (junction design, sight lines, speed compliance).
• The regenerated top 1% output contains 39,413 links: 28,435 A roads, 4,234 B roads, 4,168 classified unnumbered roads, 2,001 motorways, 568 unclassified roads, 6 unknown-class links, and 1 not-classified link.
• Excess-risk links (residual > 2) are candidates for targeted intervention where infrastructure or behaviour factors drive risk beyond what the model can attribute to road class and traffic alone.