The Disability Belt (R)¶
spdep LISA, CARBayes disease-mapping, spatialreg spillovers — at national county scale¶
The R counterpart to belt_python.ipynb, using the field-standard spatial stack on the real SSA/Census county data (3,099 contiguous-US counties, Dec 2024; the same SSI/SSDI county-participation measures behind Mathematica's interactive State Disability Maps). We first age-standardize the rates (indirect standardization against SSA's national disabled-workers-by-age table), then CARBayes fits the BYM disease-mapping model on that age-standardized offset for a smoothed relative-risk map and an exceedance map; spdep delineates the belt with local Moran (LISA); and spatialreg fits the spatial-lag model and its spillover impacts. Rates are per 1,000 working-age (18–64) residents; the DI figure is disabled workers only and SSI is suppressed for small counties. Maps throughout with sf/ggplot2.
options(repr.plot.width=11, repr.plot.height=6.5)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
suppressMessages({library(sf); library(spdep); library(spatialreg); library(CARBayes); library(ggplot2); library(patchwork)})
d <- read.csv("disability_belt.csv", colClasses=c(fips="character"))
Ed <- read.csv("belt_adj_edges.csv", colClasses="character") # county-pair edge list
n <- nrow(d); W <- matrix(0, n, n)
ii <- match(Ed$fips_i, d$fips); jj <- match(Ed$fips_j, d$fips)
W[cbind(ii, jj)] <- 1; W[cbind(jj, ii)] <- 1
sf0 <- st_read("belt_counties.geojson", quiet=TRUE)
sf0 <- sf0[match(d$fips, sf0$id), ] # align sf order to data/W order
sf0 <- cbind(sf0, d)
proj <- tryCatch(st_transform(sf0, 5070), error=function(e) sf0) # Albers; fall back to lat/long
lw <- mat2listw(W, style="W")
cat(n, "contiguous-US counties; national DI rate", round(1000*sum(d$DI)/sum(d$wa_pop),1), "per 1,000 working-age\n")
Warning message: "package 'spatialreg' was built under R version 4.6.1"
3099 contiguous-US counties; national DI rate 35.1 per 1,000 working-age
1. The belt¶
County SSDI rate (disabled-worker beneficiaries per 1,000 working-age (18-64) residents). The dark band through Central Appalachia and the Deep South is the disability belt.
ggplot(proj) + geom_sf(aes(fill=di_rate), color=NA) +
scale_fill_distiller(palette="OrRd", direction=1, limits=c(0, quantile(d$di_rate,0.98)), oob=scales::squish) +
labs(title="SSDI disabled-worker beneficiaries per 1,000 working-age (18-64) residents (Dec 2024)", fill="per 1,000") + theme_void()
2. Is the belt just an older workforce? — age-standardization¶
SSDI receipt rises steeply with age, and belt counties have older working-age populations. Indirect standardization removes age composition: using SSA's national age-specific disabled-worker rates $R^{\text{nat}}_a$, the expected count under the county's own age mix is $E_i=\sum_a N_{i,a}R^{\text{nat}}_a$ and the standardized participation ratio is $\text{SPR}_i=D_i/E_i$ (>1 = more than age predicts — participation, not a mortality/morbidity ratio, since this is program receipt). Both $E_i$ (E_age) and $\text{SPR}_i$ (spr) are precomputed in the data. The age-standardized map (right) is almost identical to the crude map (left) — the belt is not an age artifact.
cat(sprintf("region-wide SPR = sum(DI)/sum(E_age) = %.3f\n", sum(d$DI)/sum(d$E_age)))
cat(sprintf("corr(crude rate, age-standardized SPR) = %.3f; %%55-64 share explains only R^2 = %.3f of crude variance\n",
cor(d$di_rate, d$spr), cor(d$di_rate, d$pct_5564)^2))
pa <- ggplot(proj)+geom_sf(aes(fill=di_rate),color=NA)+scale_fill_distiller(palette="OrRd",direction=1,limits=c(0,quantile(d$di_rate,.98)),oob=scales::squish)+labs(title="Crude DI rate per 1,000 (18-64)",fill="/1,000")+theme_void()
pb <- ggplot(proj)+geom_sf(aes(fill=spr),color=NA)+scale_fill_distiller(palette="OrRd",direction=1,limits=c(0,quantile(d$spr,.98)),oob=scales::squish)+labs(title="Age-standardized participation ratio (DI / age-expected)",fill="SPR")+theme_void()
print(pa | pb)
cat("\ntop crude-rate counties: crude rate, %55-64 share, age-standardized SPR\n")
print(head(d[order(-d$di_rate), c("county","STNAME","di_rate","pct_5564","spr")], 6), row.names=FALSE)
region-wide SPR = sum(DI)/sum(E_age) = 0.986
corr(crude rate, age-standardized SPR) = 0.954; %55-64 share explains only R^2 = 0.137 of crude variance
top crude-rate counties: crude rate, %55-64 share, age-standardized SPR
county STNAME di_rate pct_5564 spr
Buchanan Virginia 176.7039 25.67616 4.169138
Dickenson Virginia 174.7700 24.49409 4.206611
Norton Virginia 149.9758 19.15820 4.349821
Mingo West Virginia 145.7536 24.89211 3.501126
Floyd Kentucky 142.9585 23.85177 3.568351
Magoffin Kentucky 141.4514 25.47663 3.370429
Covariate check. Adding the county's %(55–64) share to the spatial-lag model (via spatialreg::impacts) does not absorb the poverty effect — confirming age composition and economic disadvantage are separate channels.
sar1 <- lagsarlm(di_rate ~ poverty_pct, data=d, listw=lw, method="Matrix")
sar2 <- lagsarlm(di_rate ~ poverty_pct + pct_5564, data=d, listw=lw, method="Matrix")
i1 <- impacts(sar1, listw=lw, R=300); i2 <- impacts(sar2, listw=lw, R=300)
cat(sprintf("poverty TOTAL effect: %.2f (alone) -> %.2f (controlling for %%55-64)\n", i1$res$total[1], i2$res$total[1]))
cat(sprintf("%%55-64 own TOTAL effect: %.2f per 1,000 per pp\n", i2$res$total[2]))
poverty TOTAL effect: 3.83 (alone) -> 3.98 (controlling for %55-64)
%55-64 own TOTAL effect: 5.82 per 1,000 per pp
3. CAR disease-mapping with CARBayes — smoothed risk and exceedance¶
S.CARleroux fits a Leroux-CAR Poisson model to the DI counts with the age-standardized offset $\log E_i$, borrowing strength across neighbours (the Leroux CAR generalizes the intrinsic-CAR/BYM used from-scratch in Python — a spatial-dependence parameter $\rho$ mixes between independence and the intrinsic CAR, and it mixes far more reliably than S.CARbym's spatial/unstructured variance split at 3,000 areas). The relative risk map is the smoothed participation ratio, and the exceedance maps flag counties confidently above the threshold; we show $P(\text{RR}>1.5),\,P(\text{RR}>2),\,P(\text{RR}>3)$ (1.5 floods the East, 3 keeps only the peak, 2 is the informative middle). (CARBayes scales to thousands of counties.)
dat <- data.frame(DI=d$DI, E_age=d$E_age)
set.seed(1)
m <- S.CARleroux(DI ~ offset(log(E_age)), family="poisson", data=dat, W=W,
burnin=2000, n.sample=10000, thin=4, verbose=FALSE)
E <- d$E_age
proj$RR <- m$fitted.values / E
Efit <- matrix(E, nrow=nrow(m$samples$fitted), ncol=n, byrow=TRUE); RRs <- m$samples$fitted / Efit
cat(sprintf("CAR variance tau2 %.3f (n.eff %.0f), spatial rho %.2f; corr(RR, raw SPR) %.3f\n",
m$summary.results["tau2","Mean"], m$summary.results["tau2","n.effective"],
m$summary.results["rho","Mean"], cor(proj$RR, d$spr)))
mkexc <- function(thr){ proj$e <- colMeans(RRs > thr)
ggplot(proj)+geom_sf(aes(fill=e),color=NA)+scale_fill_distiller(palette="Reds",direction=1,limits=c(0,1))+
labs(title=sprintf("P(RR > %.1f) [%d cty >0.9]", thr, sum(colMeans(RRs>thr)>0.9)),fill="P")+theme_void() }
p1 <- ggplot(proj)+geom_sf(aes(fill=RR),color=NA)+scale_fill_distiller(palette="OrRd",direction=1,limits=c(0,quantile(proj$RR,.98)),oob=scales::squish)+labs(title="CAR smoothed relative risk (age-standardized)",fill="RR")+theme_void()
print(p1)
print(mkexc(1.5) | mkexc(2) | mkexc(3))
for(thr in c(1.5,2,3)) cat(sprintf("P(RR>%.1f)>0.9: %d counties\n", thr, sum(colMeans(RRs>thr)>0.9)))
cat("Smoothed-risk surface matches the Python BYM (corr ~0.99); the exceedance SET size differs by engine (uncertainty model).\n")
CAR variance tau2 0.803 (n.eff 264), spatial rho 0.75; corr(RR, raw SPR) 0.993
P(RR>1.5)>0.9: 205 counties P(RR>2.0)>0.9: 26 counties P(RR>3.0)>0.9: 2 counties
Smoothed-risk surface matches the Python BYM (corr ~0.99); the exceedance SET size differs by engine (uncertainty model).
4. LISA hotspots with spdep¶
moran.test gives the global autocorrelation; localmoran gives each county's local statistic, classified into High-High (belt hotspots), Low-Low, and outliers by comparing the county and its neighbours to the mean.
gI <- moran.test(d$di_rate, lw)$estimate[1]
lmc <- localmoran(d$di_rate, lw); zc <- d$di_rate - mean(d$di_rate); lag <- lag.listw(lw, d$di_rate) - mean(d$di_rate)
# one test per county: at an uncorrected 0.05 across ~3,100 counties ~155 are flagged by chance,
# so control the false-discovery rate (Bonferroni is unusable here -- see the Python notebook).
p_bh <- p.adjust(lmc[,5], method='BH'); sig <- p_bh < 0.05; sig_raw <- lmc[,5] < 0.05
clab <- ifelse(!sig, "not significant", ifelse(zc>0 & lag>0, "High-High (belt)", ifelse(zc<0 & lag<0, "Low-Low", ifelse(zc>0,"High-Low","Low-High"))))
proj$cluster <- factor(clab, levels=c("not significant","High-High (belt)","Low-Low","Low-High","High-Low"))
clab_raw <- ifelse(!sig_raw, "ns", ifelse(zc>0 & lag>0, "HH", "other"))
cat(sprintf("global Moran's I = %.2f; High-High belt counties: %d after Benjamini-Hochberg (%d uncorrected)\n",
gI, sum(clab=="High-High (belt)"), sum(clab_raw=="HH")))
ggplot(proj) + geom_sf(aes(fill=cluster), color=NA) +
scale_fill_manual(values=c("not significant"="grey85","High-High (belt)"="#c53030","Low-Low"="#2b6cb0","Low-High"="#9ecae1","High-Low"="#dd6b20")) +
labs(title="LISA clusters of SSDI receipt — the red High-High cluster is the disability belt", fill="") + theme_void()
global Moran's I = 0.67; High-High belt counties: 290 after Benjamini-Hochberg (440 uncorrected)
5. Spillovers with spatialreg¶
The spatial-lag model with impacts giving the direct / indirect / total effects — the spillover decomposition that quantifies how far poverty's effect reaches. The primary outcome is the age-standardized rate (the defensible estimand, from §2); the crude rate and crude + %(55–64) are shown as robustness — the total effect is invariant to the age treatment.
sar_adj <- lagsarlm(di_rate_adj ~ poverty_pct, data=d, listw=lw, method="Matrix")
cat(sprintf("OLS(age-adj) residual Moran p = %.3g; SAR(age-adjusted) rho = %.2f\n",
lm.morantest(lm(di_rate_adj ~ poverty_pct, data=d), lw)$p.value, sar_adj$rho))
imp <- impacts(sar_adj, listw=lw, R=500)
print(round(data.frame(direct=imp$res$direct, indirect=imp$res$indirect, total=imp$res$total, row.names="poverty"),2))
i_adj<-impacts(sar_adj,listw=lw,R=300); i_cr<-impacts(sar1,listw=lw,R=300); i_ca<-impacts(sar2,listw=lw,R=300)
cat(sprintf("\npoverty TOTAL effect is invariant to age treatment: age-adjusted %.2f | crude %.2f | crude+%%55-64 %.2f\n",
i_adj$res$total[1], i_cr$res$total[1], i_ca$res$total[1]))
cat("As in the Python notebook, poverty's effect is dominated by the INDIRECT (spillover) part: a strong spatial\n")
cat("multiplier (rho ~ 0.75) spreads it across neighbouring counties, and the age treatment barely moves it.\n")
OLS(age-adj) residual Moran p = 0; SAR(age-adjusted) rho = 0.75
direct indirect total poverty 1.14 2.74 3.89
poverty TOTAL effect is invariant to age treatment: age-adjusted 3.89 | crude 3.83 | crude+%55-64 3.98
As in the Python notebook, poverty's effect is dominated by the INDIRECT (spillover) part: a strong spatial
multiplier (rho ~ 0.75) spreads it across neighbouring counties, and the age treatment barely moves it.
6. Two programs, two logics — DI (insurance) vs SSI (means-tested)¶
DI is social insurance (needs a work history, not means-tested); SSI is means-tested welfare (low income/assets). The prediction: SSI receipt should be more tightly coupled to poverty (the means test is a poverty test), while DI, reflecting regional labour-market decline, should be more spatially contagious. We fit the same lagsarlm model to both on the counties where SSI is disclosed (SSA suppresses ~600 small counties; SSI here is the 18–64 blind/disabled caseload, crude working-age rate).
mask <- !is.na(d$SSI_1864); Wsub <- W[mask,mask]; dsub <- d[mask,]
deg <- rowSums(Wsub); keep <- deg>0; Wsub <- Wsub[keep,keep]; dsub <- dsub[keep,] # drop islands from subsetting
dsub$ssi_wa <- 1000*dsub$SSI_1864/dsub$wa_pop; lws <- mat2listw(Wsub, style="W")
cat(sprintf("SSI disclosed %d of %d counties; spatial model on %d\n", sum(mask), n, nrow(dsub)))
cat(sprintf("corr with poverty: DI %.2f SSI %.2f (SSI far more poverty-coupled -- the means test)\n",
cor(dsub$di_rate, dsub$poverty_pct), cor(dsub$ssi_wa, dsub$poverty_pct)))
sar_di <- lagsarlm(di_rate ~ poverty_pct, data=dsub, listw=lws, method="Matrix")
sar_ssi <- lagsarlm(ssi_wa ~ poverty_pct, data=dsub, listw=lws, method="Matrix")
cat(sprintf("DI : rho %.2f poverty coef %.2f\n", sar_di$rho, sar_di$coefficients["poverty_pct"]))
cat(sprintf("SSI: rho %.2f poverty coef %.2f (larger direct poverty link, weaker spillover)\n", sar_ssi$rho, sar_ssi$coefficients["poverty_pct"]))
psub <- proj[mask,]; psub <- psub[keep,]
psub$tilt <- as.numeric(scale(dsub$ssi_wa)) - as.numeric(scale(dsub$di_rate))
print(ggplot(psub)+geom_sf(aes(fill=tilt),color=NA)+scale_fill_distiller(palette="PuOr",direction=1,limits=c(-2.5,2.5),oob=scales::squish)+
labs(title="Where SSI runs relatively higher (purple) vs DI (orange) — the Delta joins Appalachia", fill="z-diff")+theme_void())
cat("The institutions show up: SSI's poverty coefficient is larger (means test), DI's spatial multiplier\n")
cat("is larger -- consistent with regional labour-market decline, though on one cross-section rho measures\n")
cat("association rather than diffusion. Same belt, different mechanisms.\n")
Warning message in mat2listw(Wsub, style = "W"): "neighbour object has 10 sub-graphs"
SSI disclosed 2493 of 3099 counties; spatial model on 2476
corr with poverty: DI 0.55 SSI 0.77 (SSI far more poverty-coupled -- the means test)
Warning message in sqrt(diag(fdHess)[-1]): "NaNs produced"
Warning message in sqrt(fdHess[1, 1]): "NaNs produced"
DI : rho 0.73 poverty coef 0.97
SSI: rho 0.52 poverty coef 1.47 (larger direct poverty link, weaker spillover)
The institutions show up: SSI's poverty coefficient is larger (means test), DI's spatial multiplier
is larger -- consistent with regional labour-market decline, though on one cross-section rho measures
association rather than diffusion. Same belt, different mechanisms.
7. Summary¶
The standard R spatial stack reproduces the from-scratch analysis on the real disability data at national county scale. Age-standardization (indirect, against SSA's national age-specific rates) shows the belt is not a composition artifact — the standardized-participation-ratio map matches the crude map (correlation ≈ 0.95) and the %55-64 share explains little of the variance. CARBayes (Leroux CAR) smooths the counts on the age-standardized offset into a relative-risk surface — matching the Python BYM (corr ≈ 0.99) — with exceedance maps at three thresholds. spdep confirms extremely strong autocorrelation and delineates the belt as a coherent High-High cluster through Appalachia and the Deep South. spatialreg quantifies the spillover ($\rho\approx0.75$) and the direct/indirect/total impacts on the age-standardized rate (robust to the age treatment), showing poverty's effect spreads across county lines. And the DI-vs-SSI contrast reads the programs' institutions off the estimates: means-tested SSI couples more tightly to poverty (larger direct coefficient), while insurance-based DI is more spatially contagious (larger $\rho$). Together — spdep, CARBayes, spatialreg, sf/ggplot2 — these are the tools an SSA analyst would use for regional disability surveillance, and the package mirror of belt_python.ipynb.
This is the capstone application of the spatial arc: the areal and econometric methods, built earlier on standard teaching datasets, applied to a genuine, high-stakes policy geography. Extensions: yearly SSA county data supports a spatiotemporal analysis of the belt's growth (the space-time CAR capstone), and a joint SSDI + SSI multivariate CAR (prototyped in the areal notebooks) would model the shared and program-specific geography directly.