The Disability Belt — a Spatial Application

A Real Question, Four Alternative Explanations

The arc's tools, turned on a real question. Across 3,099 contiguous-US counties, 6,993,541 disabled-worker beneficiaries receive SSDI at a national rate of 35.1 per 1,000 working-age adults — but Buchanan County, Virginia runs at 176.7. That band through Central Appalachia and the Deep South is the "disability belt". Everything here is an attempt to say something defensible about it rather than to point at a map.

Is it just an older workforce?

The first alternative explanation is age: disability receipt rises steeply with age (the national gradient runs 2.1 per 1,000 at 18–29 to 158.6 at 60–64, a factor of 77), so a belt of older workforces would produce a belt of receipt. Indirect age-standardisation settles it. The geography barely moves — crude rates and age-standardised ratios correlate 0.954, and a county's 55–64 share explains only $R^2=0.137$ of crude-rate variance. Buchanan is still at a standardised participation ratio of 4.2: four times the disabled workers its own age structure predicts. The belt is genuine excess, not demography.

Is it small-county noise?

The second is small-county noise. Bayesian disease mapping (BYM) smooths the age-standardised rate with a spatial field, and the variance split answers the question directly: a spatial SD of 0.47 against an unstructured SD of 0.20 means risk is spatially structured rather than scattered. The deliverable is the exceedance map — the counties where receipt is, with posterior probability above 0.9, at least double what age predicts. Both notebooks are careful that the size of that confident set is model-dependent (284 counties from scratch, 26 under CARBayes' Leroux CAR) even though the risk surface agrees at a correlation of ~0.99 — a difference in how much variance each model assigns the random effects, not a disagreement about where the belt is.

Is the cluster an artefact of testing 3,099 counties?

The third is multiplicity, and it bites hardest. A LISA map runs one hypothesis test per county, so across 3,099 counties an uncorrected $p<0.05$ flags roughly 155 by chance alone. Uncorrected, the belt comes to 401 High-High counties; under Benjamini–Hochberg control at $q=0.05$ it is 239 — close to the uncorrected count minus that expected false-positive share. Bonferroni is not even available: $0.05/3099$ sits below the $1/(B+1)$ floor any permutation p-value can reach. So the belt is real, coherent and large — and a third smaller than the naive count. Global Moran's I is 0.67 either way.

What drives it, and how far does it reach?

What drives it? A spatial-lag model on the age-standardised rate gives $\rho=0.75$ [0.73, 0.77] and decomposes poverty's effect into direct +1.14, indirect +2.76, total +3.90 per 1,000 per point of poverty — so roughly 70% of the association is spillover, the part an ordinary regression discards. R reproduces it at 1.14 / 2.74 / 3.89 and PyMC agrees on $\rho$ to two decimals. The total is also invariant to how age is handled (3.83 crude, 3.90 age-adjusted, 3.98 controlling for the 55–64 share). Read $\rho$ as association, not diffusion: this is one observational cross-section in which poverty, labour-market decline and receipt are all spatially smooth, so an omitted regional factor leaves the same signature as genuine spillover — a question one snapshot cannot settle, but the panel below can.

Two programmes, two logics

The two programmes then separate in a way their statutes predict. SSI — means-tested — is far more tightly coupled to poverty (correlation 0.77 against DI's 0.55) and carries the larger poverty coefficient (+1.47 against +0.97), because the means test binds receipt to household income. DI — social insurance, earned through work history — carries the larger spatial multiplier ($\rho=0.73$ against 0.52), consistent with regional labour-market decline rather than individual circumstance. Same belt, different mechanisms, with SSI tilting toward the deepest-poverty areas so the Mississippi Delta joins Appalachia.

Is the belt growing?

Finally, is the belt growing? A spatially-varying linear trend over five snapshots — 2004, 2009, 2014, 2019, 2024 — with year-specific age-standardisation so the baby-boomer ageing wave cannot masquerade as belt dynamics. The answer is unambiguous in both engines: the level SD dwarfs the trend SD (0.53 vs 0.126 from scratch; 0.45 vs 0.209 in CARBayesST), and the belt core's mean trend deviation is negative (−0.020 and −0.078). National rolls rose to a peak around 2014 and receded. The belt is entrenched, not expanding — it holds its ground while the national cycle moves around it. DI and SSI share the persistent geography (correlation 0.80) far more than they share its dynamics (0.57).

Does the belt spread?

That entrenchment does more than describe the belt — it tests the reading of $\rho$ that the cross-section could not. If receipt genuinely propagated across county lines, the belt would be expanding at its edges, so counties ringing the core should show the largest positive trend deviations. They show the opposite. The 394 counties adjacent to the core sit at a mean trend of −0.009 against +0.004 for the rest of the country, and the posterior probability that the ring grew faster than the rest is 0.00. The belt's edge is receding relative to the nation, not advancing into it. Two decades is ample time for a contagious process to show itself at the boundary; what shows up there is retreat. So the large spatial multiplier is better read as a statement about how the association is distributed across space — poverty in one county tracks receipt in its neighbours — than as a mechanism by which receipt travels. The belt is a place where disadvantage persists, not a front that is advancing.

What the analysis does not support

Stated plainly, because they bound what the analysis supports: rates use working-age (18–64) denominators and are age-standardised within that band, so the DI numerator is mildly inflated by disabled workers past 64 — a level effect, not a geographic one; SSI is left crude and is suppressed for roughly 600 small counties; the associations are ecological, about counties rather than people; and receipt reflects both underlying need and the geography of programme administration. The panel is five snapshots, not twenty-one annual observations — ample for a trend, but it estimates each county's direction of drift, not its year-to-year path.

Notebooks

Downloads

Belt Module — Source Code

"""
belt.py -- spatial analysis of the DISABILITY BELT (application; scales to ~3000 counties).

Backs the notebooks in  "The Disability Belt -- a Spatial Application".

Applies the arc's areal + econometric tools to a real, policy-relevant problem: the geography of
SSDI ("DI", disabled-worker beneficiaries) and SSI (blind/disabled recipients) receipt across the
~3,100 contiguous-US counties (SSA county data, Dec 2024; Census population + SAIPE poverty/income).
The "disability belt" -- Central Appalachia and the Deep South -- is the motivating pattern.

Two tools, both built to scale to thousands of counties:

  LOCAL MORAN (LISA):  formally delineates the belt. The local statistic I_i = z_i * sum_j w_ij z_j
  flags each county as part of a High-High cluster (a hotspot -- the belt), a Low-Low cluster, or a
  spatial outlier, with a permutation p-value. This turns "the belt" from an eyeballed pattern into
  a set of statistically-significant hotspot counties.

  SPATIAL LAG (SAR):  y = rho W y + X beta + eps, to ask what DRIVES the rate (poverty, income) and
  whether it SPILLS across county lines. At n ~ 3000 the exact log-determinant log|I - rho W| =
  sum_i log(1 - rho lambda_i) is computed once from the eigenvalues of the row-standardised W, so the
  Metropolis step is O(n). The LeSage-Pace effects are also O(n) per draw: with a row-standardised W,
  W 1 = 1 so the average TOTAL effect has the closed form beta / (1 - rho), and the average DIRECT
  effect is beta * (1/n) sum_i 1/(1 - rho lambda_i) via the same eigenvalues -- no n-by-n inverse.
"""

import numpy as np


def adjacency_from_edges(fips_i, fips_j, fips):
    """expand a county-pair edge list into the dense symmetric 0/1 adjacency, aligned to `fips`.
    A contiguity graph on ~3,000 counties has only ~9,200 edges, so the pair list is ~160x smaller
    than the matrix it encodes -- the dense form is only ever needed in memory, never on disk."""
    pos = {str(f): k for k, f in enumerate(fips)}
    n = len(pos); W = np.zeros((n, n))
    for a, b in zip(fips_i, fips_j):
        i = pos.get(str(a)); j = pos.get(str(b))
        if i is not None and j is not None:
            W[i, j] = 1.0; W[j, i] = 1.0
    return W


def row_standardise(W):
    W = np.asarray(W, float); rs = W.sum(1, keepdims=True); rs[rs == 0] = 1
    return W / rs


def morans_i(y, Wrs):
    z = np.asarray(y, float) - np.mean(y); return (len(y) / Wrs.sum()) * (z @ Wrs @ z) / (z @ z)


def local_moran(y, Wrs, rng, nperm=499):
    """local Moran's I per county with a permutation p-value and cluster quadrant.
    quadrant: 1=High-High (hotspot), 2=Low-Low, 3=Low-High, 4=High-Low. Returns (Ii, p, quadrant)."""
    y = np.asarray(y, float); z = (y - y.mean()); s2 = (z @ z) / len(y); zz = z / s2
    lag = Wrs @ z; Ii = zz * lag
    n = len(y); nb = (Wrs > 0)
    perm = np.empty((nperm, n))
    for b in range(nperm):
        zp = z[rng.permutation(n)]
        perm[b] = zz * (Wrs @ zp)
    p = (1 + (perm >= Ii).sum(0)) / (nperm + 1)
    p = np.minimum(p, 1 - p) * 2                                    # two-sided-ish
    quad = np.where(z > 0, np.where(lag > 0, 1, 4), np.where(lag > 0, 3, 2))
    return Ii, p, quad


def fdr_bh(p, q=0.05):
    """Benjamini-Hochberg mask at false-discovery rate q. A LISA map runs one test per county, so at
    an uncorrected 0.05 across ~3,000 counties roughly 150 are flagged by chance; FDR is the standard
    correction here (Bonferroni is both far too strict and, with a permutation p-value floor of
    1/(nperm+1), usually unreachable)."""
    p = np.asarray(p, float); m = len(p)
    ps = np.sort(p); ok = np.where(ps <= np.arange(1, m + 1) * q / m)[0]
    return p <= ps[ok.max()] if len(ok) else np.zeros(m, bool)


def sar_gibbs(y, X, W, rng, draws=3000, burn=1500, b_sd=1e3):
    """Bayesian spatial-lag (SAR) model at scale. Conditional on rho the model is OLS on (I-rhoW)y,
    so beta,sigma^2 are conjugate; rho is Metropolis with the eigenvalue log-determinant. Returns
    draws of beta, rho, sigma, plus the eigenvalues of the row-standardised W (for the effects)."""
    y = np.asarray(y, float); X = np.asarray(X, float); n, k = X.shape
    Wrs = row_standardise(W); Wy = Wrs @ y
    evals = np.linalg.eigvals(Wrs).real
    lo, hi = 1.0 / evals.min() + 1e-4, 1.0 / evals.max() - 1e-4
    rho = 0.0; sig2 = np.var(y); step = 0.05; P0 = np.eye(k) / b_sd ** 2; acc = 0
    B = np.empty((draws, k)); R = np.empty(draws); S = np.empty(draws)
    def ld(r): return np.sum(np.log(1 - r * evals))
    for it in range(draws + burn):
        ystar = y - rho * Wy                                        # (I-rhoW)y
        V = np.linalg.inv(X.T @ X / sig2 + P0); m = V @ (X.T @ ystar / sig2)
        beta = m + np.linalg.cholesky(V) @ rng.standard_normal(k)
        resid = ystar - X @ beta; sig2 = (resid @ resid) / rng.chisquare(n)
        rp = rho + step * rng.standard_normal()
        if lo < rp < hi:
            rc = ystar - X @ beta; rn = (y - rp * Wy) - X @ beta
            if np.log(rng.random()) < (ld(rp) - ld(rho)) - (rn @ rn - rc @ rc) / (2 * sig2):
                rho = rp; acc += 1
        if it < burn and it % 200 == 199:
            step *= np.exp((acc / 200 - 0.4) * 0.5); acc = 0
        if it >= burn:
            j = it - burn; B[j] = beta; R[j] = rho; S[j] = np.sqrt(sig2)
    return dict(beta=B, rho=R, sigma=S, evals=evals)


def sar_effects(beta_draws, rho_draws, evals, col):
    """LeSage-Pace direct/indirect/total effects for covariate `col`, O(n) per draw via eigenvalues
    (row-standardised W): total = beta/(1-rho), direct = beta*(1/n)sum 1/(1-rho*lambda_i)."""
    n = len(evals); b = beta_draws[:, col]
    direct = b * np.array([np.mean(1.0 / (1 - r * evals)) for r in rho_draws])
    total = b / (1 - rho_draws)
    return direct, total - direct, total


def ols(y, X):
    XtXi = np.linalg.inv(X.T @ X); b = XtXi @ (X.T @ y)
    r = y - X @ b; s2 = (r @ r) / (len(y) - X.shape[1])
    return b, np.sqrt(np.diag(s2 * XtXi)), r


def bym_gibbs(y, E, W, rng, draws=1500, burn=1000, a=1.0, b=0.01):
    """Besag-York-Mollie disease-mapping model, scaled to thousands of counties. With an expected-count
    offset E_i (here the AGE-STANDARDIZED expected disabled workers):

        y_i ~ Poisson( E_i * exp(alpha + phi_i + theta_i) )

    phi ~ ICAR spatial field (smooths toward neighbours), theta ~ iid N(0, tau_theta^2) unstructured
    heterogeneity. The fitted relative risk RR_i = exp(alpha + phi_i + theta_i) is a SMOOTHED
    standardized participation ratio -- it borrows strength from neighbours, stabilising small-county
    rates -- and its posterior gives EXCEEDANCE probabilities P(RR_i > c). Metropolis-within-Gibbs: alpha/phi_i/theta_i Metropolis on
    the Poisson likelihood (phi_i with the ICAR neighbour-mean conditional; the neighbour sum is
    updated incrementally on acceptance, so a sweep is O(edges) not O(n^2)), tau^2 conjugate InvGamma.
    Returns draws of alpha, phi, theta and the spatial/unstructured SDs."""
    y = np.asarray(y, float); E = np.asarray(E, float); W = np.asarray(W, float)
    n = len(y); nnb = W.sum(1); nnb[nnb == 0] = 1
    nbr = [np.where(W[i] > 0)[0] for i in range(n)]
    alpha = np.log(max(y.sum(), 1) / E.sum()); phi = np.zeros(n); theta = np.zeros(n)
    t2p = 0.3; t2t = 0.1; sal = 0.02; sp = 0.4 * np.ones(n); st = 0.4 * np.ones(n)
    accp = np.zeros(n); acct = np.zeros(n)
    A = np.empty(draws); PHI = np.empty((draws, n)); TH = np.empty((draws, n))
    SDP = np.empty(draws); SDT = np.empty(draws)
    nbsum = W @ phi
    for it in range(draws + burn):
        lam = E * np.exp(alpha + phi + theta)
        # ---- alpha (Poisson intercept, vague N(0,100) prior) ----
        d = sal * rng.standard_normal(); lp = lam * np.exp(d)
        if np.log(rng.random()) < (y.sum() * d - (lp - lam).sum() - ((alpha + d) ** 2 - alpha ** 2) / 200):
            alpha += d; lam = lp
        # ---- phi_i  (ICAR; incremental neighbour sum) ----
        for i in range(n):
            d = sp[i] * rng.standard_normal(); nbmean = nbsum[i] / nnb[i]
            dll = y[i] * d - lam[i] * (np.exp(d) - 1)
            dpr = -(nnb[i] / (2 * t2p)) * ((phi[i] + d - nbmean) ** 2 - (phi[i] - nbmean) ** 2)
            if np.log(rng.random()) < dll + dpr:
                phi[i] += d; lam[i] *= np.exp(d); nbsum[nbr[i]] += d; accp[i] += 1
        phi -= phi.mean(); nbsum = W @ phi
        # ---- theta_i  (iid heterogeneity) ----
        for i in range(n):
            d = st[i] * rng.standard_normal()
            dll = y[i] * d - lam[i] * (np.exp(d) - 1)
            dpr = -((theta[i] + d) ** 2 - theta[i] ** 2) / (2 * t2t)
            if np.log(rng.random()) < dll + dpr:
                theta[i] += d; lam[i] *= np.exp(d); acct[i] += 1
        # ---- variances (conjugate InvGamma) ----
        qp = 0.5 * (nnb * phi * phi).sum() - 0.5 * (phi * (W @ phi)).sum()
        t2p = 1.0 / rng.gamma(a + (n - 1) / 2, 1.0 / (b + max(qp, 1e-6)))
        t2t = 1.0 / rng.gamma(a + n / 2, 1.0 / (b + 0.5 * (theta @ theta)))
        if it < burn and it % 100 == 99:
            sp *= np.exp((accp / 100 - 0.44) * 0.4); st *= np.exp((acct / 100 - 0.44) * 0.4)
            accp[:] = 0; acct[:] = 0
        if it >= burn:
            j = it - burn; A[j] = alpha; PHI[j] = phi; TH[j] = theta
            SDP[j] = np.sqrt(t2p); SDT[j] = np.sqrt(t2t)
    return dict(alpha=A, phi=PHI, theta=TH, spatial_sd=SDP, hetero_sd=SDT)


def st_linear_gibbs(y, E, W, cvec, rng, draws=1500, burn=1000, a=1.0, b=0.01):
    """SPATIO-TEMPORAL model of the belt over time: a spatially-varying linear time trend (the model
    CARBayesST fits with ST.CARlinear). For county i, year t (centred time c_t):

        y_it ~ Poisson( E_it * exp( alpha + phi_i + gamma_t + b_i * c_t ) )

    phi_i ~ ICAR persistent spatial level (the time-averaged belt);
    gamma_t (sum-to-zero) the NATIONAL trajectory over time (soaks up the shared trend, incl. its
        linear part, so b is a pure departure);
    b_i ~ ICAR (sum-to-zero) the county-specific TREND DEVIATION -- b_i>0 means county i's relative
        receipt GREW faster than the nation, b_i<0 that it receded. Mapping b_i answers "where has the
        belt expanded / held / contracted." E_it is the year-specific AGE-STANDARDIZED expected count,
        so the national ageing wave is already removed. Metropolis-within-Gibbs with incremental
        neighbour sums (a sweep is O(edges)); tau^2 conjugate InvGamma. Returns draws of alpha, phi,
        gamma, b and the level/trend SDs."""
    y=np.asarray(y,float); E=np.asarray(E,float); W=np.asarray(W,float); c=np.asarray(cvec,float)
    n,T=y.shape; nnb=W.sum(1); nnb[nnb==0]=1; nbr=[np.where(W[i]>0)[0] for i in range(n)]
    logE=np.log(np.maximum(E,1e-9))
    alpha=np.log(max(y.sum(),1)/E.sum()); phi=np.zeros(n); gamma=np.zeros(T); bb=np.zeros(n)
    t2p=0.3; t2b=0.05; sal=0.02; sp=0.3*np.ones(n); sg=0.1*np.ones(T); sb=0.1*np.ones(n)
    accp=np.zeros(n); accb=np.zeros(n)
    A=np.empty(draws); PHI=np.empty((draws,n)); GAM=np.empty((draws,T)); B=np.empty((draws,n))
    SDP=np.empty(draws); SDB=np.empty(draws)
    nbsum_p=W@phi; nbsum_b=W@bb
    def eta(): return logE+alpha+phi[:,None]+gamma[None,:]+np.outer(bb,c)
    for it in range(draws+burn):
        lam=np.exp(eta())
        # ---- alpha ----
        d=sal*rng.standard_normal(); lp=lam*np.exp(d)
        if np.log(rng.random())<(y.sum()*d-(lp-lam).sum()-((alpha+d)**2-alpha**2)/200): alpha+=d; lam=lp
        # ---- phi_i (ICAR level; likelihood summed over time) ----
        yr=y.sum(1)
        for i in range(n):
            d=sp[i]*rng.standard_normal(); nbm=nbsum_p[i]/nnb[i]
            dll=yr[i]*d-lam[i].sum()*(np.exp(d)-1)
            dpr=-(nnb[i]/(2*t2p))*((phi[i]+d-nbm)**2-(phi[i]-nbm)**2)
            if np.log(rng.random())<dll+dpr: phi[i]+=d; lam[i]*=np.exp(d); nbsum_p[nbr[i]]+=d; accp[i]+=1
        phi-=phi.mean(); nbsum_p=W@phi
        # ---- gamma_t (national trajectory, sum-to-zero, vague) ----
        lam=np.exp(eta()); yc=y.sum(0)
        for t in range(T):
            d=sg[t]*rng.standard_normal()
            dll=yc[t]*d-lam[:,t].sum()*(np.exp(d)-1)-((gamma[t]+d)**2-gamma[t]**2)/(2*10.0)
            if np.log(rng.random())<dll: gamma[t]+=d; lam[:,t]*=np.exp(d)
        gamma-=gamma.mean()
        # ---- b_i (ICAR trend deviation; likelihood weighted by c_t) ----
        lam=np.exp(eta()); yb=(y*c[None,:]).sum(1)
        for i in range(n):
            d=sb[i]*rng.standard_normal(); nbm=nbsum_b[i]/nnb[i]
            dll=yb[i]*d-(lam[i]*(np.exp(d*c)-1)).sum()
            dpr=-(nnb[i]/(2*t2b))*((bb[i]+d-nbm)**2-(bb[i]-nbm)**2)
            if np.log(rng.random())<dll+dpr: bb[i]+=d; lam[i]*=np.exp(d*c); nbsum_b[nbr[i]]+=d; accb[i]+=1
        bb-=bb.mean(); nbsum_b=W@bb
        # ---- smoothing variances ----
        qp=0.5*(nnb*phi*phi).sum()-0.5*(phi*(W@phi)).sum(); t2p=1.0/rng.gamma(a+(n-1)/2,1.0/(b+max(qp,1e-6)))
        qb=0.5*(nnb*bb*bb).sum()-0.5*(bb*(W@bb)).sum();     t2b=1.0/rng.gamma(a+(n-1)/2,1.0/(b+max(qb,1e-6)))
        if it<burn and it%100==99:
            sp*=np.exp((accp/100-0.44)*0.4); sb*=np.exp((accb/100-0.44)*0.4); accp[:]=0; accb[:]=0
        if it>=burn:
            j=it-burn; A[j]=alpha; PHI[j]=phi; GAM[j]=gamma; B[j]=bb; SDP[j]=np.sqrt(t2p); SDB[j]=np.sqrt(t2b)
    return dict(alpha=A, phi=PHI, gamma=GAM, b=B, level_sd=SDP, trend_sd=SDB)


def indirect_standardize(counts, band_pops, nat_rates):
    """INDIRECT age-standardization (the disease-mapping standard, and the E_i offset used throughout
    the areal arc). band_pops is (n, A) county population by age band; nat_rates is (A,) the NATIONAL
    age-specific rate per band. The expected count is what each county WOULD have under national
    age-specific rates applied to ITS OWN age composition:

        E_i = sum_a  band_pops[i, a] * nat_rates[a]

    and the standardized ratio SPR_i = counts_i / E_i (>1 = more than age structure alone predicts).
    SPR = "standardized participation ratio" (this is program RECEIPT, not deaths/disease, so not an
    SMR/SIR). Indirect (not direct) because the county's own age-specific counts are unknown -- only the total
    counts_i and the age composition band_pops are observed; that is exactly the SSA county situation.
    Returns (E, SMR)."""
    band_pops = np.asarray(band_pops, float); nat_rates = np.asarray(nat_rates, float)
    E = band_pops @ nat_rates
    return E, np.asarray(counts, float) / E

References

• Anselin, L. (1995). Local indicators of spatial association — LISA. Geographical Analysis 27(2), 93–115. — the local statistic used to delineate the belt
• Caldas de Castro, M. & Singer, B. H. (2006). Controlling the false discovery rate: a new application to account for multiple and dependent tests in local statistics of spatial association. Geographical Analysis 38(2), 180–208. — why FDR rather than Bonferroni on a LISA map
• Besag, J., York, J. & Mollié, A. (1991). Bayesian image restoration, with two applications in spatial statistics. AISM 43(1), 1–20. — the smoothing model behind the exceedance maps
• LeSage, J. & Pace, R. K. (2009). Introduction to Spatial Econometrics. CRC Press. — the direct/indirect decomposition, and the eigenvalue methods that make it feasible at this scale
• Autor, D. H. & Duggan, M. G. (2003). The rise in the disability rolls and the decline in unemployment. Quarterly Journal of Economics 118(1), 157–206. — the labour-market reading of regional variation in DI receipt
• Social Security Administration (2024). SSI Recipients by State and County and Annual Statistical Report on the Social Security Disability Insurance Program. — the county counts and the national age-specific rates