import numpy as np, pandas as pd, matplotlib.pyplot as plt, geopandas as gpd
import belt as B
rng = np.random.default_rng(1)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#cfcfcf"
d = pd.read_csv("disability_belt.csv", dtype={"fips":str})
_e = pd.read_csv("belt_adj_edges.csv", dtype=str)                  # county-pair edge list
W = B.adjacency_from_edges(_e["fips_i"], _e["fips_j"], d["fips"]); Wrs = B.row_standardise(W); n=len(d)
gdf = gpd.read_file("belt_counties.geojson").set_crs(4326).to_crs(5070)          # Albers equal-area CONUS
gdf = gdf.merge(d.assign(order=range(n)), left_on="id", right_on="fips").sort_values("order").reset_index(drop=True)
def choropleth(ax, values, title, cmap="OrRd", vmax=None, categorical=False, cats=None):
    g=gdf.copy(); g["v"]=values
    if categorical:
        g.plot(column="v", ax=ax, categorical=True, legend=True, edgecolor="none",
               legend_kwds={"loc":"lower left","fontsize":7}, cmap=cmap)
    else:
        g.plot(column="v", ax=ax, cmap=cmap, edgecolor="none", legend=True,
               legend_kwds={"shrink":0.5}, vmax=vmax)
    ax.set_title(title, fontsize=11); ax.axis("off")
print(f"{n} contiguous-US counties;  {int(d['DI'].sum()):,} SSDI disabled-worker beneficiaries")
print(f"DI rate per 1,000 working-age (18-64): mean {d['di_rate'].mean():.1f}, max {d['di_rate'].max():.1f} (Buchanan Co., VA)")

3099 contiguous-US counties;  6,993,541 SSDI disabled-worker beneficiaries
DI rate per 1,000 working-age (18-64): mean 52.5, max 176.7 (Buchanan Co., VA)

fig,ax=plt.subplots(figsize=(12,6.5))
choropleth(ax, d["di_rate"].values, "SSDI disabled-worker beneficiaries per 1,000 working-age (18-64) residents (Dec 2024)", cmap="OrRd", vmax=np.percentile(d["di_rate"],98))
plt.tight_layout(); plt.show()
fig,ax=plt.subplots(figsize=(12,6.5))
ssi=d["ssi_rate"].values
choropleth(ax, np.where(np.isnan(ssi),0,ssi), "SSI blind/disabled recipients per 1,000 working-age (18-64) residents", cmap="BuPu", vmax=np.nanpercentile(ssi,98))
plt.tight_layout(); plt.show()
print(f"national DI rate {1000*d['DI'].sum()/d['wa_pop'].sum():.1f}/1,000 working-age;  the belt counties run 3-5x higher.")
print("SSDI and SSI trace the same belt -- Central Appalachia and the Deep South -- the geography this project explains.")

national DI rate 35.1/1,000 working-age;  the belt counties run 3-5x higher.
SSDI and SSI trace the same belt -- Central Appalachia and the Deep South -- the geography this project explains.

# national age-specific DI rate per 1,000 (SSA SSDI Annual Statistical Report 2024, Table 19; disabled
# workers by age, Dec 2024) -- steeply increasing:
bands   = ["18-29","30-34","35-39","40-44","45-49","50-54","55-59","60-64"]
popcols = ["pop_1829","pop_3034","pop_3539","pop_4044","pop_4549","pop_5054","pop_5559","pop_6064"]
Rnat    = np.array([2.06,6.15,10.92,18.50,28.46,43.02,72.64,158.59]) / 1000.0        # per person
E_age, spr = B.indirect_standardize(d["DI"].values, d[popcols].values, Rnat)         # E_i and SPR (participation ratio)
print(f"national age gradient: {1000*Rnat[0]:.1f}/1,000 (18-29)  ->  {1000*Rnat[-1]:.1f}/1,000 (60-64)  = {Rnat[-1]/Rnat[0]:.0f}x")
print(f"region-wide SPR = sum(DI)/sum(E) = {d['DI'].sum()/E_age.sum():.3f}  (national rates well-calibrated to CONUS)")
print(f"corr(crude rate, age-standardized SPR) = {np.corrcoef(d['di_rate'],spr)[0,1]:.3f}  -> the geography barely moves")
print(f"county %55-64 share explains only R^2 = {np.corrcoef(d['di_rate'],d['pct_5564'])[0,1]**2:.3f} of crude-rate variance")
fig,ax=plt.subplots(1,2,figsize=(14,4.6))
choropleth(ax[0], d["di_rate"].values, "Crude DI rate per 1,000 (18-64)", cmap="OrRd", vmax=np.percentile(d["di_rate"],98))
choropleth(ax[1], spr, "Age-standardized participation ratio (DI ÷ age-expected)", cmap="OrRd", vmax=np.percentile(spr,98))
plt.tight_layout(); plt.show()
print("\ntop crude-rate counties: crude rate, %55-64 share, and age-standardized SPR")
print(d.assign(spr=spr).nlargest(6,"di_rate")[["county","STNAME","di_rate","pct_5564","spr"]].round(1).to_string(index=False))
print("\nBuchanan Co. VA is 5x the national CRUDE rate; age-standardized it is still SPR ~4.2 -- four times the disabled")
print("workers its age structure predicts. Age composition explains only a modest slice; the belt is genuine excess.")

national age gradient: 2.1/1,000 (18-29)  ->  158.6/1,000 (60-64)  = 77x
region-wide SPR = sum(DI)/sum(E) = 0.986  (national rates well-calibrated to CONUS)
corr(crude rate, age-standardized SPR) = 0.954  -> the geography barely moves
county %55-64 share explains only R^2 = 0.137 of crude-rate variance

top crude-rate counties: crude rate, %55-64 share, and age-standardized SPR
   county        STNAME  di_rate  pct_5564  spr
 Buchanan      Virginia    176.7      25.7  4.2
Dickenson      Virginia    174.8      24.5  4.2
   Norton      Virginia    150.0      19.2  4.3
    Mingo West Virginia    145.8      24.9  3.5
    Floyd      Kentucky    143.0      23.9  3.6
 Magoffin      Kentucky    141.5      25.5  3.4

Buchanan Co. VA is 5x the national CRUDE rate; age-standardized it is still SPR ~4.2 -- four times the disabled
workers its age structure predicts. Age composition explains only a modest slice; the belt is genuine excess.

yb=d["di_rate"].values
Xp =np.column_stack([np.ones(n), d["poverty_pct"].values])
Xpa=np.column_stack([np.ones(n), d["poverty_pct"].values, d["pct_5564"].values])
sp_ =B.sar_gibbs(yb,Xp ,W,rng,draws=1200,burn=700); _,_,Tp =B.sar_effects(sp_["beta"],sp_["rho"],sp_["evals"],1)
spa =B.sar_gibbs(yb,Xpa,W,rng,draws=1200,burn=700); _,_,Tpa=B.sar_effects(spa["beta"],spa["rho"],spa["evals"],1)
_,_,Ta=B.sar_effects(spa["beta"],spa["rho"],spa["evals"],2)
print(f"poverty TOTAL effect on DI rate:  {Tp.mean():+.2f} (poverty alone)  ->  {Tpa.mean():+.2f} (controlling for %55-64)")
print(f"%55-64 share own TOTAL effect: {Ta.mean():+.2f} per 1,000 per pp  (older workforce does raise receipt -- but does not explain away poverty)")
print("Poverty's effect survives age controls intact: age composition and economic disadvantage are separate channels.")

poverty TOTAL effect on DI rate:  +3.85 (poverty alone)  ->  +3.97 (controlling for %55-64)
%55-64 share own TOTAL effect: +5.81 per 1,000 per pp  (older workforce does raise receipt -- but does not explain away poverty)
Poverty's effect survives age controls intact: age composition and economic disadvantage are separate channels.

bym=B.bym_gibbs(d["DI"].values, d["E_age"].values, W, rng, draws=1200, burn=800)
RR=np.exp(bym["alpha"][:,None]+bym["phi"]+bym["theta"]); rr=RR.mean(0)
print(f"spatial SD {bym['spatial_sd'].mean():.2f} dominates unstructured SD {bym['hetero_sd'].mean():.2f}  -> risk is spatially structured, not scattered noise")
print(f"corr(BYM smoothed RR, raw SPR) = {np.corrcoef(rr,d['spr'])[0,1]:.3f}  (smoothing stabilises small counties, preserves the belt)")
fig,ax=plt.subplots(figsize=(12,6))
choropleth(ax, rr, "BYM smoothed relative risk (age-standardized)", cmap="OrRd", vmax=np.percentile(rr,98))
plt.tight_layout(); plt.show()
# exceedance at three thresholds -- the choice tunes stringency
fig,ax=plt.subplots(1,3,figsize=(18,4.4))
for k,c in enumerate([1.5,2,3]):
    exc=(RR>c).mean(0)
    choropleth(ax[k], exc, f"P(RR > {c})   [{(exc>0.9).sum()} counties >0.9]", cmap="Reds", vmax=1.0)
    print(f"P(RR>{c})>0.9 : {(exc>0.9).sum()} counties")
plt.tight_layout(); plt.show()
print("\nhighest smoothed-RR counties (the age-standardized belt core):")
print(d.assign(rr=rr,P=(RR>2).mean(0)).nlargest(6,"rr")[["county","STNAME","rr","P"]].round(2).to_string(index=False))
print("\nThe exceedance maps are the policy deliverable: at RR>2, the counties where SSDI receipt is, with high posterior")
print("probability, at least double what age structure alone predicts. The smoothed-risk SURFACE is cross-checked by")
print("CARBayes' Leroux CAR (R notebook, corr ~0.99) and by PyMC's pm.ICAR at county scale (NC-SIDS areal notebook).")
print("(The confident-set SIZE is model-dependent -- it reflects how much variance each model assigns the random")
print(" effects -- so exact exceedance counts differ across engines even when the risk map agrees.)")

spatial SD 0.47 dominates unstructured SD 0.20  -> risk is spatially structured, not scattered noise
corr(BYM smoothed RR, raw SPR) = 0.999  (smoothing stabilises small counties, preserves the belt)

P(RR>1.5)>0.9 : 890 counties

P(RR>2)>0.9 : 284 counties

P(RR>3)>0.9 : 20 counties

highest smoothed-RR counties (the age-standardized belt core):
   county   STNAME   rr   P
   Norton Virginia 4.33 1.0
Dickenson Virginia 4.19 1.0
 Buchanan Virginia 4.15 1.0
    Perry  Alabama 3.62 1.0
     Wise Virginia 3.59 1.0
    Floyd Kentucky 3.56 1.0

The exceedance maps are the policy deliverable: at RR>2, the counties where SSDI receipt is, with high posterior
probability, at least double what age structure alone predicts. The smoothed-risk SURFACE is cross-checked by
CARBayes' Leroux CAR (R notebook, corr ~0.99) and by PyMC's pm.ICAR at county scale (NC-SIDS areal notebook).
(The confident-set SIZE is model-dependent -- it reflects how much variance each model assigns the random
 effects -- so exact exceedance counts differ across engines even when the risk map agrees.)

I = B.morans_i(d["di_rate"].values, Wrs)
NPERM = 9999                                                       # p-value floor 1/(NPERM+1)
Ii, p, quad = B.local_moran(d["di_rate"].values, Wrs, rng, nperm=NPERM)
labels = np.array(["not significant","High-High (belt)","Low-Low","Low-High","High-Low"])
# one test per county, so correct for multiplicity: at an uncorrected 0.05 across 3,099
# counties ~155 would be flagged by chance alone.
sig = B.fdr_bh(p, 0.05)
cl_raw = np.where(p<0.05, quad, 0)
cl = np.where(sig, quad, 0)                                        # 0 = not significant
d["cluster"] = labels[cl]
print(f"global Moran's I = {I:.2f}  -> extremely strong spatial clustering of DI receipt")
for k in [1,2,3,4]: print(f"  {labels[k]:20s}: {(cl==k).sum()} counties")
from matplotlib.colors import ListedColormap
order=["not significant","High-High (belt)","Low-Low","Low-High","High-Low"]
cmap=ListedColormap([GREY,RED,BLUE,"#9ecae1",ORANGE])
gdf2=gdf.copy(); gdf2["cluster"]=pd.Categorical(d["cluster"],categories=order)
fig,ax=plt.subplots(figsize=(12,6.5))
gdf2.plot(column="cluster",ax=ax,cmap=cmap,edgecolor="none",legend=True,legend_kwds={"loc":"lower left","fontsize":8})
ax.set_title("LISA clusters of SSDI receipt — the red High-High cluster IS the disability belt",fontsize=11); ax.axis("off")
plt.tight_layout(); plt.show()
print(f"\nMultiplicity matters at this scale. An uncorrected p<0.05 flags {(cl_raw==1).sum()} High-High counties, but")
print(f"~{0.05*len(p):.0f} of 3,099 would be flagged by chance alone; after Benjamini-Hochberg control at q=0.05,")
print(f"{(cl==1).sum()} survive -- close to the uncorrected count minus that expected false-positive share.")
print("Bonferroni is not an option here: 0.05/3099 is below the 1/(NPERM+1) floor a permutation p-value can reach.")
print(f"So the belt is real and large -- {(cl==1).sum()} counties through Appalachia and the Deep South, delineated")
print("statistically rather than by eye -- but its headline size shrinks by a third once multiplicity is handled.")

global Moran's I = 0.67  -> extremely strong spatial clustering of DI receipt
  High-High (belt)    : 239 counties
  Low-Low             : 287 counties
  Low-High            : 8 counties
  High-Low            : 2 counties

Multiplicity matters at this scale. An uncorrected p<0.05 flags 401 High-High counties, but
~155 of 3,099 would be flagged by chance alone; after Benjamini-Hochberg control at q=0.05,
239 survive -- close to the uncorrected count minus that expected false-positive share.
Bonferroni is not an option here: 0.05/3099 is below the 1/(NPERM+1) floor a permutation p-value can reach.
So the belt is real and large -- 239 counties through Appalachia and the Deep South, delineated
statistically rather than by eye -- but its headline size shrinks by a third once multiplicity is handled.

specs=[("age-adjusted rate ~ poverty (primary)", d["di_rate_adj"].values, np.column_stack([np.ones(n), d["poverty_pct"].values])),
       ("crude rate ~ poverty",                   d["di_rate"].values,     np.column_stack([np.ones(n), d["poverty_pct"].values])),
       ("crude rate ~ poverty + %55-64",          d["di_rate"].values,     np.column_stack([np.ones(n), d["poverty_pct"].values, d["pct_5564"].values]))]
print("poverty's TOTAL effect (per 1,000 per +1pp poverty) is invariant to how age is handled:")
res=[]
for tag,yy,XX in specs:
    s=B.sar_gibbs(yy,XX,W,rng,draws=1500,burn=800); Dd,Ii,Tt=B.sar_effects(s["beta"],s["rho"],s["evals"],1)
    res.append((s,Dd,Ii,Tt)); print(f"  {tag:38s}: rho {s['rho'].mean():.2f},  poverty total {Tt.mean():+.2f}")
sar,D,Ind,T=res[0]; rho=sar["rho"]; y=specs[0][1]; X=specs[0][2]      # primary = age-adjusted (carried to PyMC)
b_ols,se_ols,resid=B.ols(y,X); print(f"\n(OLS on the age-adjusted rate leaves residual Moran's I {B.morans_i(resid,Wrs):.2f} -> a spatial model is needed)")
print(f"PRIMARY (age-adjusted): rho {rho.mean():.2f} [{np.percentile(rho,2.5):.2f},{np.percentile(rho,97.5):.2f}];  poverty direct {D.mean():+.2f}, indirect {Ind.mean():+.2f}, total {T.mean():+.2f}")
fig,ax=plt.subplots(1,2,figsize=(12,3.8))
ax[0].hist(rho,bins=40,color=BLUE,alpha=.8,density=True); ax[0].axvline(rho.mean(),color=RED,lw=2)
ax[0].set_xlabel(r"spatial spillover $\rho$"); ax[0].set_title(f"Strong spillover: rho = {rho.mean():.2f}")
ax[1].bar(["direct","indirect\n(spillover)","total"],[D.mean(),Ind.mean(),T.mean()],color=[BLUE,ORANGE,GREEN])
ax[1].set_ylabel("age-adj DI/1,000 per +1pp poverty"); ax[1].set_title("Poverty's effect is mostly SPILLOVER")
plt.tight_layout(); plt.show()
print(f"rho = {rho.mean():.2f} is large, so poverty's TOTAL effect")
print(f"({T.mean():+.2f} per 1,000 per point of poverty) is ~{T.mean()/D.mean():.1f}x its direct effect -- most of the impact SPILLS")
print("onto neighbours -- the policy-relevant quantity an ordinary regression understates, and the age treatment")
print("barely moves it. Read rho as association, not diffusion: this is one observational cross-section with a")
print("single covariate, and poverty, labour-market decline and receipt are all spatially smooth, so an omitted")
print("regional factor produces the same signature as genuine spillover. The lag specification fits better than")
print("the alternatives; that is a statement about error structure, not evidence that receipt spreads county to")
print("county. Establishing the latter needs variation the map alone does not contain -- and the panel in the")
print("spatio-temporal notebook supplies exactly that: if receipt diffused, the belt would be expanding at its")
print("edges over two decades, and it is not.")

poverty's TOTAL effect (per 1,000 per +1pp poverty) is invariant to how age is handled:

  age-adjusted rate ~ poverty (primary) : rho 0.75,  poverty total +3.90

  crude rate ~ poverty                  : rho 0.76,  poverty total +3.84

  crude rate ~ poverty + %55-64         : rho 0.71,  poverty total +4.03

(OLS on the age-adjusted rate leaves residual Moran's I 0.62 -> a spatial model is needed)
PRIMARY (age-adjusted): rho 0.75 [0.73,0.77];  poverty direct +1.14, indirect +2.76, total +3.90

rho = 0.75 is large, so poverty's TOTAL effect
(+3.90 per 1,000 per point of poverty) is ~3.4x its direct effect -- most of the impact SPILLS
onto neighbours -- the policy-relevant quantity an ordinary regression understates, and the age treatment
barely moves it. Read rho as association, not diffusion: this is one observational cross-section with a
single covariate, and poverty, labour-market decline and receipt are all spatially smooth, so an omitted
regional factor produces the same signature as genuine spillover. The lag specification fits better than
the alternatives; that is a statement about error structure, not evidence that receipt spreads county to
county. Establishing the latter needs variation the map alone does not contain.

import pymc as pm, pytensor.tensor as pt
evals=sar["evals"]; lo,hi=1/evals.min()+1e-4,1/evals.max()-1e-4; Wy=Wrs@y
with pm.Model() as mod:
    rho_=pm.Uniform("rho",lo,hi); beta=pm.Normal("beta",0,100,shape=2); sig=pm.HalfNormal("sig",30)
    resid=(y-rho_*Wy)-(beta[0]+beta[1]*X[:,1])
    ll=pt.sum(pt.log(1-rho_*evals)) - 0.5*n*pt.log(2*np.pi*sig**2) - pt.sum(resid**2)/(2*sig**2)
    pm.Potential("lik",ll)
    idata=pm.sample(800,tune=1200,chains=4,target_accept=0.9,random_seed=3,progressbar=False)
print(f"rho:          from-scratch {rho.mean():.2f}   PyMC {float(idata.posterior['rho'].mean()):.2f}")
print(f"poverty coef: from-scratch {sar['beta'].mean(0)[1]:+.2f}   PyMC {float(idata.posterior['beta'].mean(('chain','draw'))[1]):+.2f}")
print("PyMC agrees: the strong spatial spillover in disability receipt is robust to the sampler.")

g++ not available, if using conda: `conda install gxx`

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [rho, beta, sig]

Sampling 4 chains for 1_200 tune and 800 draw iterations (4_800 + 3_200 draws total) took 6 seconds.

rho:          from-scratch 0.75   PyMC 0.75
poverty coef: from-scratch +0.97   PyMC +0.97
PyMC agrees: the strong spatial spillover in disability receipt is robust to the sampler.

mask=d["SSI_1864"].notna().to_numpy()
Ws=W[np.ix_(mask,mask)]; ds=d[mask].reset_index(drop=True)
keep=Ws.sum(1)>0; Ws=Ws[np.ix_(keep,keep)]; ds=ds[keep].reset_index(drop=True); ns=len(ds)   # drop islands created by subsetting
ds["ssi_wa"]=1000*ds["SSI_1864"]/ds["wa_pop"]; Wrs_s=B.row_standardise(Ws)
print(f"SSI disclosed for {mask.sum()} of {n} counties; SSI spatial model on {ns} (islands dropped)")
print(f"corr with poverty:   DI {ds['di_rate'].corr(ds['poverty_pct']):.2f}   SSI {ds['ssi_wa'].corr(ds['poverty_pct']):.2f}   (SSI far more poverty-coupled -- the means test)")
Xs=np.column_stack([np.ones(ns), ds["poverty_pct"].values]); out={}
for nm,yy in [("DI",ds["di_rate"].values), ("SSI",ds["ssi_wa"].values)]:
    s=B.sar_gibbs(yy,Xs,Ws,rng,draws=1500,burn=800); Dd,Ii,Tt=B.sar_effects(s["beta"],s["rho"],s["evals"],1); out[nm]=(s,Dd,Ii,Tt)
    print(f"  {nm:3s}: rho {s['rho'].mean():.2f}   poverty coef {s['beta'].mean(0)[1]:+.2f}   (direct {Dd.mean():+.2f}, indirect {Ii.mean():+.2f}, total {Tt.mean():+.2f})")
# where does SSI run relatively higher than DI? (z-score tilt) -- the Deep-South/Delta shift
ds["tilt"]=((ds["ssi_wa"]-ds["ssi_wa"].mean())/ds["ssi_wa"].std()) - ((ds["di_rate"]-ds["di_rate"].mean())/ds["di_rate"].std())
gsub=gdf[mask].reset_index(drop=True); gsub=gsub[keep].reset_index(drop=True); gsub["tilt"]=ds["tilt"].values
fig,ax=plt.subplots(1,2,figsize=(14,4.6))
gsub.plot(column="tilt",ax=ax[0],cmap="PuOr_r",edgecolor="none",legend=True,legend_kwds={"shrink":0.5},vmin=-2.5,vmax=2.5)
ax[0].set_title("Where SSI runs relatively higher (purple) vs DI (orange)",fontsize=10); ax[0].axis("off")
ax[1].bar(["DI\ncoef","SSI\ncoef"],[out["DI"][0]["beta"].mean(0)[1],out["SSI"][0]["beta"].mean(0)[1]],color=[BLUE,RED])
ax[1].set_ylabel("poverty coefficient"); ax[1].set_title("SSI's direct poverty link is larger")
plt.tight_layout(); plt.show()
print("\nThe programs' institutions show up in the estimates. SSI's poverty COEFFICIENT is larger (the means test binds")
print("receipt to individual poverty), but DI's spatial SPILLOVER rho is larger -- DI reflects regional labour-market")
print("decline that spreads across county lines, while SSI is set locally by household income. Same belt, different")
print("mechanisms: SSI tilts toward the deepest-poverty areas (the Mississippi Delta joins Appalachia).")

SSI disclosed for 2493 of 3099 counties; SSI spatial model on 2476 (islands dropped)
corr with poverty:   DI 0.55   SSI 0.77   (SSI far more poverty-coupled -- the means test)

  DI : rho 0.72   poverty coef +0.98   (direct +1.16, indirect +2.37, total +3.53)

  SSI: rho 0.52   poverty coef +1.47   (direct +1.58, indirect +1.47, total +3.04)

The programs' institutions show up in the estimates. SSI's poverty COEFFICIENT is larger (the means test binds
receipt to individual poverty), but DI's spatial SPILLOVER rho is larger -- DI reflects regional labour-market
decline that spreads across county lines, while SSI is set locally by household income. Same belt, different
mechanisms: SSI tilts toward the deepest-poverty areas (the Mississippi Delta joins Appalachia).

The Disability Belt — a Spatial Application¶

Mapping and explaining SSDI/SSI receipt across US counties¶

1. The belt¶

2. Is the belt just an older workforce? — age-standardization¶

3. Smoothing the map — Bayesian disease mapping (BYM) and exceedance¶

4. Is the belt a statistical cluster, and where? — Moran's I and LISA hotspots¶

5. What drives it? — a spatial-lag model of poverty and spillovers¶

6. Cross-check in PyMC¶

7. Two programs, two logics — DI (insurance) vs SSI (means-tested)¶

8. Summary¶