Spatial Epi Public Health GIS

Applying FDR Correction to Local Moran’s I Results

This guide, part of Global & Local Moran’s I Implementation, deals with the multiple-testing problem that every Local Moran’s I map carries and most published maps ignore.

Problem Context & Constraints

Local Moran’s I produces one statistic and one pseudo p-value per area. A county map has 60 to 250 areas, a tract map several thousand. Testing each at 0.05 means that, with no clustering at all, about 5 percent of areas will be labelled significant anyway — 50 spurious high-high or low-low areas on a map of a thousand tracts. Because neighbouring statistics share data, those chance flags tend to appear in small contiguous groups, which look exactly like real clusters.

Three corrections are in common use. Bonferroni divides the significance level by the number of tests; it controls the chance of any false flag but is so strict on large maps that real clusters vanish. False discovery rate (FDR) control, usually by the Benjamini–Hochberg procedure, bounds the expected proportion of flagged areas that are false; it is much less conservative and fits the exploratory purpose of a LISA map. Neighbourhood-based adjustments, such as dividing by the typical neighbour count, are heuristics without a formal guarantee.

FDR is the sensible default for surveillance mapping, but applying it to Local Moran’s I has two practical complications. Permutation p-values are discrete: with 999 permutations, the smallest possible value is 0.001, and FDR thresholds on a large map can fall below that, making it impossible for any area to pass. And local statistics are correlated, which the standard Benjamini–Hochberg procedure tolerates for positive dependence but which should be acknowledged.

Chance Clusters on a Map With No Clustering Two maps of the same randomly generated data with no true clustering. On the left, uncorrected Local Moran's I at the five percent level flags several small groups of adjacent areas as high-high or low-low, which look like clusters. On the right, after false discovery rate correction, no areas are flagged. Because local statistics share neighbours, chance flags cluster together and are easily mistaken for real patterns. Random data, no clusters — two different maps Uncorrected, p < 0.05 7 areas flagged, all by chance FDR-corrected, q < 0.05 no areas flagged the correct answer for this data chance flags form contiguous groups because neighbouring statistics share data

Prerequisites

Step-by-Step Solution

Choose the permutation count from the map size so that FDR thresholds are reachable, run Local Moran’s I, apply Benjamini–Hochberg to the pseudo p-values, and keep both the raw and adjusted labels.

# Local Moran's I with Benjamini-Hochberg false discovery rate control.
# Pinned: numpy==1.26.4, geopandas==0.14.4, libpysal==4.10, esda==2.5.1, statsmodels==0.14.2
import logging
import math

import numpy as np
import geopandas as gpd
from libpysal.weights import Queen
from esda.moran import Moran_Local
from statsmodels.stats.multitest import multipletests

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("lisa_fdr")

QUAD = {1: "HH", 2: "LH", 3: "LL", 4: "HL"}


def permutations_for(n_areas: int, q: float = 0.05, target_hits: int = 5) -> int:
    """Permutations needed so that target_hits areas could pass BH at level q."""
    smallest_threshold = q * target_hits / n_areas
    perms = math.ceil(1 / smallest_threshold) * 2 - 1
    return max(999, min(perms, 99_999))


def lisa_fdr(gdf: gpd.GeoDataFrame, col: str, q: float = 0.05, seed: int = 20260918) -> gpd.GeoDataFrame:
    w = Queen.from_dataframe(gdf, use_index=False)
    w.transform = "r"
    islands = w.islands
    if islands:
        log.warning("%d areas have no neighbours; excluded from testing", len(islands))

    perms = permutations_for(len(gdf), q)
    lm = Moran_Local(gdf[col].to_numpy(dtype=float), w, permutations=perms, seed=seed)

    testable = np.ones(len(gdf), dtype=bool)
    testable[islands] = False
    p = lm.p_sim.astype(float)
    reject = np.zeros(len(gdf), dtype=bool)
    q_vals = np.full(len(gdf), np.nan)
    reject[testable], q_vals[testable], _, _ = multipletests(p[testable], alpha=q, method="fdr_bh")

    out = gdf.copy()
    out["lisa_i"] = lm.Is
    out["p_sim"] = p
    out["q_bh"] = q_vals
    out["quad"] = [QUAD[k] for k in lm.q]
    out["label_raw"] = np.where((p <= q) & testable, out["quad"], "ns")
    out["label_fdr"] = np.where(reject, out["quad"], "ns")
    log.info("n=%d perms=%d raw-significant=%d fdr-significant=%d",
             len(gdf), perms, int(((p <= q) & testable).sum()), int(reject.sum()))
    return out

The permutation count matters more than it seems. Benjamini–Hochberg compares the k-th smallest p-value to q·k/n. On a map of 2,000 tracts at q = 0.05, the first comparison is against 0.000025 — forty times smaller than the smallest p-value 999 permutations can produce. With too few permutations, FDR rejects nothing, and the analyst concludes there is no clustering when the test never had the resolution to say so. permutations_for sizes the count so that a handful of areas could, in principle, pass.

Keep label_raw alongside label_fdr. The raw labels are useful for exploration and hypothesis generation; the corrected labels are what goes in a report that says “significant”. Keeping both also shows how much the correction changed: a map where 164 raw flags shrink to 31 corrected ones tells a reader that most of the uncorrected pattern was noise, which is itself worth knowing when an older uncorrected map is still circulating among partners.

How Benjamini–Hochberg Draws the Line Sorted p-values plotted by rank rise slowly from near zero. A straight threshold line rises from zero with slope q divided by n. The largest rank where a p-value sits below the line sets the cut-off, and every area at or below that rank is flagged. A horizontal floor marks the smallest p-value that the chosen number of permutations can produce; when the threshold line starts below that floor, no area can pass. The cut-off is the last rank under the line rank of p-value (smallest first) q · k / n cut-off permutation floor 1/(perms+1) if the dashed line starts under the floor, add permutations

Validation & Edge Cases

Nothing survives correction. An empty corrected map is a legitimate result, but check the permutation count first. If permutations_for hit its cap of 99,999 and the map is very large, consider testing at a coarser geography or focusing on a pre-specified region rather than concluding there is no clustering.

Islands. Areas with no neighbours have undefined local statistics. Exclude them from the set of tests, as the code does, so they neither inflate n nor contribute arbitrary p-values, and list them in the output.

Correlated tests. Benjamini–Hochberg controls the FDR under positive regression dependence, which local Moran statistics on a contiguity graph roughly satisfy. Under arbitrary dependence the Benjamini–Yekutieli variant (method="fdr_by") is valid but much more conservative. Use Benjamini–Hochberg by default and state the assumption; switch to Benjamini–Yekutieli only when a reviewer requires a guarantee under any dependence.

Two-sided versus one-sided interest. The pseudo p-value from esda is one-sided in the direction of the observed statistic, folded so it reflects extremeness either way. If only high-high clusters matter for the program, test only those areas — but decide that before looking at the map, and correct over that pre-specified subset.

Seeds and reproducibility. Permutation p-values change slightly between runs; near the FDR cut-off, an area can flip between significant and not. Fix the seed for reproducibility, and for areas whose q-values sit within 10 percent of the threshold, report them as borderline rather than presenting a hard boundary that a rerun could move.

Correcting over the wrong family. The set of tests being corrected should match the question. A statewide map corrected over all tracts, then filtered to one county for a local report, keeps the statewide correction; re-running the correction over just that county’s tracts gives a less strict threshold and more flags. Neither is wrong, but they answer different questions, and switching between them after seeing the results is a quiet way to manufacture significance. Fix the family — statewide, region or county — in the analysis protocol.

Reporting q-values rather than labels. A binary flagged/not-flagged map hides how close each area came. Publishing the q-value for each area, or shading by q-value bands, lets readers see that an area at q = 0.052 is nearly as strong as one at 0.048. This is especially useful for program staff who must decide whether a borderline area warrants a phone call.

Comparing maps across time. Each period’s map is corrected within itself. An area that passes in March and not in April may simply have moved across the threshold by noise. For trend statements, use a space-time method rather than comparing corrected LISA maps month by month.

INFO lisa_fdr n=1987 perms=7951 raw-significant=164 fdr-significant=31
WARNING lisa_fdr 3 areas have no neighbours; excluded from testing
Three Corrections, Three Answers Horizontal bars showing flagged areas on a map of about two thousand tracts with one real cluster. No correction flags one hundred and sixty-four areas, most of them chance. False discovery rate control flags thirty-one, concentrated on the real cluster. Bonferroni flags six, only the core of the real cluster. FDR keeps the cluster's extent while removing most chance flags. Flagged areas, one real cluster on ~2,000 tracts none 164 FDR (BH) 31 Bonferroni 6 FDR keeps the cluster's extent while removing most chance flags

Compliance Notes

Record the correction method and level, the permutation count and seed, the number of tests after excluding islands, and both raw and corrected counts. Any published map that uses the word “significant” should state which labels it shows.

Presenting an uncorrected LISA map as a map of significant clusters overstates the evidence, and decisions such as targeted screening or resource allocation made from it will be partly driven by chance. Where uncorrected maps are shared for exploration, label them as such in the legend, not only in the methods text.