Mapping Permutation Importance Across Space

Compute permutation importance per region instead of globally, map the result, and reveal features that dominate in one part of a study area and do nothing in another.

Compute permutation importance separately inside each region, then map the result. A single global table says “elevation matters 0.29”; a map says “elevation carries the model in the uplands and does nothing on the coastal plain, where impervious surface takes over”. The second statement is actionable — it tells you where the model is relying on something that may not transfer, and where a regional sub-model or an extra feature would help.

This page is the per-region method. For the global machinery it builds on — block permutation, SHAP, collinearity — see Model Explainability for Spatial Predictions.

Permute inside one region, score that region, repeat The study area is divided into six regions. For each region in turn, the values of a single feature are shuffled among the rows belonging to that region only, and the score drop measured on that region's rows becomes the region's importance value for that feature. The other regions are untouched during each measurement. One model, one feature, six independent measurements region 2 shuffle elevation here only score region 2 → drop 0.31 repeat 0.340.31 0.020.11 an importance surface, not a number high low The global mean of these four numbers is 0.20 — a value none of the regions has.

Why This Fails in Geospatial ML Pipelines

Global importance averages away the thing you most want to see. Geospatial processes are rarely stationary: the drivers of flooding in a mountain catchment are not the drivers on a floodplain, and a model fitted across both learns both rules and applies whichever the features indicate. Averaged, the two sets of drivers each look moderately important, and the resulting bar chart describes no place in particular. Teams then make feature decisions — dropping a “weak” variable — that are wrong for half the study area.

The second failure is refitting per region. It seems natural: subset the data, fit, inspect. But now each region has its own model, and comparing their importances compares models, not places. The question worth answering is how the deployed model behaves in each region, which requires one fixed model and per-region permutation.

Third, small regions produce noise that looks like signal. A region with forty samples has a score estimate with wide error bars, so its importance value swings between repeats. Mapped without a sample-count layer, those swings render as bright patches that invite over-interpretation — the same trap as an unsmoothed rate map in epidemiology, and the same fix: report the denominator.

Core Principles

  • One model, many measurements. Never refit per region.
  • Permute within the region only. Rows outside it stay untouched.
  • Score on that region’s rows. A global score would dilute the effect.
  • Repeat and report variability. Ten repeats minimum, with a standard deviation.
  • Carry the sample count into the map. Thin regions must be visually discountable.
  • Use the same regions as your evaluation folds where possible, so the two analyses speak to each other.

Production-Ready Code

from __future__ import annotations

import logging
import geopandas as gpd
import numpy as np
import pandas as pd
from sklearn.metrics import roc_auc_score

logger = logging.getLogger(__name__)


def regional_permutation_importance(
    model, x: pd.DataFrame, y: np.ndarray, regions: np.ndarray,
    features: list[str] | None = None, n_repeats: int = 10,
    min_samples: int = 100, seed: int = 0,
) -> pd.DataFrame:
    """Permutation importance computed inside each region, for one fixed model.

    Args:
        model: A FITTED model. It is never refitted.
        regions: Region label per row of `x`.
        min_samples: Regions below this are reported with importance NaN.

    Returns:
        Long DataFrame: region, feature, auc_drop, auc_drop_std, n_samples.
    """
    rng = np.random.default_rng(seed)
    features = features or list(x.columns)
    rows = []

    for region in pd.unique(regions):
        mask = regions == region
        n = int(mask.sum())
        xr, yr = x.loc[mask], y[mask]

        if n < min_samples or len(np.unique(yr)) < 2:
            logger.info("region %s: %d samples, skipped", region, n)
            for f in features:
                rows.append({"region": region, "feature": f, "auc_drop": np.nan,
                             "auc_drop_std": np.nan, "n_samples": n})
            continue

        base = roc_auc_score(yr, model.predict_proba(xr)[:, 1])
        for f in features:
            drops = []
            for _ in range(n_repeats):
                xp = xr.copy()
                xp[f] = rng.permutation(xp[f].to_numpy())
                drops.append(base - roc_auc_score(yr, model.predict_proba(xp)[:, 1]))
            rows.append({"region": region, "feature": f,
                         "auc_drop": float(np.mean(drops)),
                         "auc_drop_std": float(np.std(drops)),
                         "n_samples": n})
        logger.info("region %s: %d samples, base AUC %.3f", region, n, base)

    return pd.DataFrame(rows)


def to_importance_layer(result: pd.DataFrame, region_geoms: gpd.GeoDataFrame,
                        feature: str, region_col: str = "region") -> gpd.GeoDataFrame:
    """Join one feature's per-region importance onto the region polygons for mapping."""
    sub = result.loc[result["feature"] == feature,
                     ["region", "auc_drop", "auc_drop_std", "n_samples"]]
    out = region_geoms.merge(sub, left_on=region_col, right_on="region", how="left")
    out["reliable"] = out["n_samples"] >= 100
    if out["auc_drop"].isna().all():
        raise ValueError(f"no region produced an importance for {feature!r}")
    return out


def dominant_feature_map(result: pd.DataFrame, region_geoms: gpd.GeoDataFrame,
                         region_col: str = "region") -> gpd.GeoDataFrame:
    """Which feature dominates in each region — the single most useful map here."""
    best = (result.dropna(subset=["auc_drop"])
                  .sort_values("auc_drop", ascending=False)
                  .groupby("region", as_index=False)
                  .first()[["region", "feature", "auc_drop", "n_samples"]])
    return region_geoms.merge(best, left_on=region_col, right_on="region", how="left")

Step-by-Step Walkthrough

Step 1 — Define regions that mean something

import geopandas as gpd

samples = gpd.read_file("samples.gpkg")
ecoregions = gpd.read_file("ecoregions.gpkg").to_crs(samples.crs)

joined = gpd.sjoin(samples, ecoregions[["eco_name", "geometry"]],
                   how="left", predicate="within")
regions = joined["eco_name"].fillna("unassigned").to_numpy()
print(pd.Series(regions).value_counts().head())

Ecoregions, administrative units and coarse H3 cells all work. What matters is that the boundary corresponds to something that could plausibly change the process — using arbitrary squares produces a map that is harder to interpret, though still valid.

Step 2 — Run the per-region measurement

result = regional_permutation_importance(
    model, X_test, y_test, regions[test_mask],
    features=["elevation", "slope_deg", "tpi_11", "ndvi_mean",
              "impervious_pct", "rain_24h", "upslope_area"],
    n_repeats=10, min_samples=150, seed=0,
)
print(result.pivot_table(index="region", columns="feature", values="auc_drop").round(3))

Step 3 — Map the surface for one feature

layer = to_importance_layer(result, ecoregions, feature="upslope_area",
                            region_col="eco_name")
layer.to_file("artifacts/importance_upslope_area.gpkg", driver="GPKG")

print(layer.loc[layer["reliable"], ["eco_name", "auc_drop", "n_samples"]]
           .sort_values("auc_drop", ascending=False).head())

Render with the unreliable regions hatched or greyed rather than dropped. A blank region is information — it says the model could not be evaluated there.

Step 4 — Map the dominant feature

dominant = dominant_feature_map(result, ecoregions, region_col="eco_name")
print(dominant[["eco_name", "feature", "auc_drop"]].to_string(index=False))

This single map usually changes the conversation. When one half of the study area is driven by terrain and the other by land use, the options — a regional interaction term, a stratified model, an extra feature for the weak half — become obvious in a way no global chart makes them.

The dominant-feature map is the output people act on On the left a global bar chart shows three features with similar importance. On the right a map of the same study area shows the upland region dominated by upslope area, the lowland region dominated by impervious surface, and the coastal region dominated by rainfall, information the bar chart cannot express. global ranking dominant feature by region upslope area impervious % rainfall 24h “all three matter, roughly equally” no action follows from this uplands · upslope area (0.34) lowlands · impervious % (0.29) coast · rainfall 24h (0.22) “three processes, one model” now the next experiment is obvious

Verification

import numpy as np
import pandas as pd


def test_permutation_stays_inside_its_region():
    """Rows outside the region under test must be untouched."""
    x = pd.DataFrame({"a": np.arange(20.0), "b": np.arange(20.0)})
    regions = np.array(["north"] * 10 + ["south"] * 10)
    rng = np.random.default_rng(0)

    xp = x.copy()
    mask = regions == "north"
    xp.loc[mask, "a"] = rng.permutation(xp.loc[mask, "a"].to_numpy())
    assert np.array_equal(xp.loc[~mask, "a"], x.loc[~mask, "a"])


def test_thin_regions_report_nan_not_a_number():
    res = regional_permutation_importance(model, X_small, y_small, regions_small,
                                          min_samples=1_000)
    assert res["auc_drop"].isna().all()
    assert (res["n_samples"] < 1_000).all()


def test_model_is_not_refitted():
    """The importance routine must not mutate the model."""
    before = model.get_booster().save_raw() if hasattr(model, "get_booster") else None
    regional_permutation_importance(model, X_test.head(500), y_test[:500],
                                    regions[:500], n_repeats=2)
    after = model.get_booster().save_raw() if hasattr(model, "get_booster") else None
    assert before == after


def test_dominant_map_has_one_row_per_region():
    dom = dominant_feature_map(result, ecoregions, region_col="eco_name")
    assert len(dom) == len(ecoregions)

Beyond tests, compare the regional map against the residual map from the same held-out set. Regions with large residuals and uniformly low importance are the model’s blind spots: it is neither accurate there nor relying on anything in particular, which usually means a driver is missing from the feature set entirely.

Reading the importance map against the residual map A two by two matrix crosses low and high regional residuals with low and high regional importance. Low residual and high importance is healthy. High residual with low importance indicates a missing driver. High residual with high importance indicates over-reliance on a feature that does not transfer. Low residual with low importance indicates the region is easy for reasons the features do not capture. Cross the two maps before drawing conclusions low importance high importance low residual high residual easy region accurate, but not via these features healthy accurate and clearly driven missing driver nothing in the feature set explains it over-reliant leaning hard on a feature that misleads here

FAQ

How large should the regions be?

Large enough for a few hundred samples and both classes, small enough that the process could plausibly differ between them. Ecoregions, administrative units or coarse hexagons all work; report the sample count so thin regions can be discounted.

Should the model be refitted per region?

No — refitting compares different models rather than different places. Keep one fixed, deployed model and permute inside each region, which measures how that model behaves locally.

What does a region with near-zero importance for every feature mean?

Usually that the model predicts a near-constant value there, so shuffling changes nothing. That is a real finding: the model has no local discrimination in that region, which the global score hides completely.

Part of: Model Explainability for Spatial Predictions Part of: Training Geospatial Predictive Models in Python