The Disability Belt Over Time (R) — CARBayesST::ST.CARlinear¶
Spatially-varying time trend on the real SSDI panel — five snapshots, 2004–2024¶
The panel is five county-level snapshots — 2004, 2009, 2014, 2019 and 2024 — not twenty-one annual observations. That is ample for a linear trend per county, which is what is fitted, but it is worth being exact about: what follows estimates the direction and slope of each county's drift across two decades, not its year-to-year path.
The R counterpart to belt_spacetime_python.ipynb. ST.CARlinear fits exactly the from-scratch model — a persistent spatial intercept $\phi_i$ (the time-averaged belt) and a spatially-varying slope $\delta_i$ (each county's trend relative to the national trajectory) — on the identical five-snapshot panel, with the year-specific age-standardized offset $\log E_{it}$. Mapping the slope $\delta_i$ answers where the belt grew, held, or receded.
options(repr.plot.width=13, repr.plot.height=5, warn=-1)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
suppressMessages({library(CARBayesST); library(sf); library(ggplot2); library(patchwork)})
P <- read.csv("belt_panel.csv", colClasses=c(fips="character"))
years <- c(2004,2009,2014,2019,2024); fips <- sort(unique(P$fips))
Ed <- read.csv("belt_panel_adj_edges.csv", colClasses="character") # county-pair edge list
nW <- length(fips); W <- matrix(0, nW, nW)
ii <- match(Ed$fips_i, fips); jj <- match(Ed$fips_j, fips)
W[cbind(ii, jj)] <- 1; W[cbind(jj, ii)] <- 1
# order data SPACE-within-TIME (all counties for t1, then t2, ...), counties in the W (sorted-fips) order
dat <- do.call(rbind, lapply(years, function(yr){ s<-P[P$year==yr,]; s<-s[match(fips,s$fips),]
data.frame(fips=fips, year=yr, DI=s$DI, E_age=s$E_age, spr=s$spr) }))
cat(sprintf("panel %d counties x %d years = %d rows; W %d x %d\n", length(fips), length(years), nrow(dat), nrow(W), ncol(W)))
geo <- st_read("belt_counties.geojson", quiet=TRUE); geo <- geo[geo$id %in% fips, ]
geo <- geo[match(fips, geo$id), ]
natrate <- sapply(years, function(y) 1000*sum(P$DI[P$year==y])/sum(P$wa_pop[P$year==y]))
cat("national crude DI rate /1,000:", paste(round(natrate,1), collapse=", "), "\n")
panel 3035 counties x 5 years = 15175 rows; W 3035 x 3035
national crude DI rate /1,000: 33.2, 39.7, 44.6, 41.2, 35.1
1. The panel and the national arc¶
National SSDI rose to a 2014 peak and fell back. The age-standardized participation ratio (SPR) mapped at each snapshot shows the belt holding its ground through the whole cycle.
maps <- lapply(years, function(yr){ geo$spr <- dat$spr[dat$year==yr]
ggplot(geo)+geom_sf(aes(fill=spr),color=NA)+scale_fill_distiller(palette="OrRd",direction=1,limits=c(0,3),oob=scales::squish)+
labs(title=yr,fill="SPR")+theme_void()+theme(plot.title=element_text(size=10),legend.position="none") })
print(wrap_plots(maps, nrow=1) + plot_annotation(title="Age-standardized participation ratio (SPR) — the belt at each snapshot, 2004–2024"))
options(repr.plot.width=7, repr.plot.height=4)
di <- sapply(years, function(y) sum(P$DI[P$year==y])/1e6)
print(ggplot(data.frame(year=years,di=di,rate=natrate))+geom_line(aes(year,di),color="#2b6cb0",linewidth=1)+geom_point(aes(year,di),color="#2b6cb0",size=2)+
geom_vline(xintercept=2014,linetype=3)+labs(y="disabled workers (millions)",title="National SSDI: rise to 2014 peak, then decline")+theme_minimal())
2. Fit ST.CARlinear¶
The spatially-varying linear-trend model: DI ~ offset(log(E_age)), Poisson, with the $3{,}035\times3{,}035$ spatial weight matrix. It returns a spatial intercept $\phi_i$ (persistent belt), a global time slope, and spatial slope deviations $\delta_i$ (Leroux CAR on both).
set.seed(1)
m <- ST.CARlinear(DI ~ offset(log(E_age)), family="poisson", data=dat, W=W,
burnin=2000, n.sample=10000, thin=4, verbose=FALSE)
cat("samples:", paste(names(m$samples), collapse=", "), "\n")
phi <- colMeans(m$samples$phi) # persistent spatial level
delta <- colMeans(m$samples$delta) # spatial slope deviation (trend vs national)
cat(sprintf("global time slope (alpha): %.3f; level SD %.2f >> slope SD %.3f (belt entrenched, not shifting)\n",
mean(m$samples$alpha), sd(phi), sd(delta)))
samples: beta, alpha, phi, delta, tau2, rho, fitted, Y
global time slope (alpha): 0.145; level SD 0.45 >> slope SD 0.209 (belt entrenched, not shifting)
3. The persistent belt and the trend map¶
$\phi_i$ is the time-averaged belt (left); $\delta_i$ the county trend relative to the nation (right) — red grew, blue receded. As in the Python fit, the coalfield core is a high, stable plateau, while the movement is idiosyncratic and local.
geo$phi <- phi; geo$delta <- delta; lim <- quantile(abs(delta),0.98)
p1 <- ggplot(geo)+geom_sf(aes(fill=exp(phi)),color=NA)+scale_fill_distiller(palette="OrRd",direction=1,limits=c(0.4,quantile(exp(phi),.98)),oob=scales::squish)+labs(title="Persistent belt exp(phi)",fill="RR")+theme_void()
p2 <- ggplot(geo)+geom_sf(aes(fill=delta),color=NA)+scale_fill_gradient2(low="#2b6cb0",mid="grey90",high="#c53030",limits=c(-lim,lim),oob=scales::squish)+labs(title="Trend deviation delta_i (red grew / blue receded)",fill="slope")+theme_void()
options(repr.plot.width=13, repr.plot.height=5); print(p1 | p2)
db <- read.csv("disability_belt.csv", colClasses=c(fips="character"))
info <- merge(data.frame(fips=fips,phi=phi,delta=delta), db[,c('fips','county','STNAME')], by='fips')
cat('persistent belt core (highest phi):\n'); print(head(info[order(-info$phi),c('county','STNAME','phi','delta')],6), row.names=FALSE)
core <- info[info$phi>quantile(info$phi,0.9),]
cat(sprintf('\nbelt core mean trend delta = %+.3f (stable-to-receding; entrenched, not expanding)\n', mean(core$delta)))
persistent belt core (highest phi):
county STNAME phi delta
Dickenson Virginia 1.324759 -0.3021643
Buchanan Virginia 1.302006 -0.6373676
Norton Virginia 1.238376 -0.1224567
Floyd Kentucky 1.147661 -0.2075540
Leslie Kentucky 1.095059 -0.3223878
Martin Kentucky 1.087653 -0.2789697
belt core mean trend delta = -0.078 (stable-to-receding; entrenched, not expanding)
4. SSI over time — a different program, a different signature¶
DI is social insurance; SSI is means-tested welfare. Both national rolls peak around 2014, but DI booms and busts while SSI is flatter and gently declining (2024 below 2004). We fit ST.CARlinear to SSI (18–64 blind/disabled) on the counties SSA discloses every year. (SSI is population-standardized, not age-standardized — SSA publishes no SSI-by-age, and its age gradient is far flatter than DI's.)
disc <- aggregate(!is.na(P$SSI_1864), list(fips=P$fips), sum); sf <- sort(disc$fips[disc$x==5])
keep <- fips %in% sf; Ws <- W[keep,keep]; sf <- fips[keep]; deg <- rowSums(Ws)
Ws <- Ws[deg>0,deg>0]; sf <- sf[deg>0]; ns <- length(sf)
Ps <- P[P$fips %in% sf,]
natssi <- sapply(years,function(y) 1000*sum(Ps$SSI_1864[Ps$year==y])/sum(Ps$wa_pop[Ps$year==y]))
natdi2 <- sapply(years,function(y) 1000*sum(Ps$DI[Ps$year==y])/sum(Ps$wa_pop[Ps$year==y]))
options(repr.plot.width=7.5, repr.plot.height=4)
print(ggplot(data.frame(year=rep(years,2),rate=c(natdi2,natssi),prog=rep(c("DI (insurance)","SSI (means-tested)"),each=length(years))))+
geom_line(aes(year,rate,color=prog),linewidth=1)+geom_point(aes(year,rate,color=prog),size=2)+geom_vline(xintercept=2014,linetype=3)+
scale_color_manual(values=c("DI (insurance)"="#2b6cb0","SSI (means-tested)"="#c53030"),name="")+
labs(y="rate per 1,000 working-age",title="Two programs, two time signatures")+theme_minimal()+theme(legend.position="top"))
dsi <- do.call(rbind, lapply(years, function(yr){ s<-Ps[Ps$year==yr,]; s<-s[match(sf,s$fips),]
data.frame(fips=sf, SSI=s$SSI_1864, E=s$wa_pop*natssi[match(yr,years)]/1000) }))
set.seed(1)
ms <- ST.CARlinear(SSI ~ offset(log(E)), family="poisson", data=dsi, W=Ws, burnin=2000, n.sample=8000, thin=4, verbose=FALSE)
phis <- colMeans(ms$samples$phi); deltas <- colMeans(ms$samples$delta)
cat(sprintf("SSI: level SD %.2f >> slope SD %.3f (also entrenched); n=%d counties\n", sd(phis), sd(deltas), ns))
SSI: level SD 0.61 >> slope SD 0.301 (also entrenched); n=2842 counties
gs <- geo[match(sf, geo$id),]; gs$phis <- phis; gs$deltas <- deltas; lim <- quantile(abs(deltas),0.98)
q1 <- ggplot(gs)+geom_sf(aes(fill=exp(phis)),color=NA)+scale_fill_distiller(palette="BuPu",direction=1,limits=c(0.4,quantile(exp(phis),.98)),oob=scales::squish)+labs(title="Persistent SSI belt exp(phi)",fill="RR")+theme_void()
q2 <- ggplot(gs)+geom_sf(aes(fill=deltas),color=NA)+scale_fill_gradient2(low="#2b6cb0",mid="grey90",high="#c53030",limits=c(-lim,lim),oob=scales::squish)+labs(title="SSI trend deviation (red grew / blue receded)",fill="slope")+theme_void()
options(repr.plot.width=13, repr.plot.height=5); print(q1 | q2)
phi_di <- phi[match(sf,fips)]; delta_di <- delta[match(sf,fips)]
cat(sprintf("DI vs SSI share the PERSISTENT belt: corr(phi_DI, phi_SSI) = %.2f\n", cor(phi_di,phis)))
cat(sprintf("DI vs SSI dynamics agree less: corr(delta_DI, delta_SSI) = %.2f\n", cor(delta_di,deltas)))
DI vs SSI share the PERSISTENT belt: corr(phi_DI, phi_SSI) = 0.80
DI vs SSI dynamics agree less: corr(delta_DI, delta_SSI) = 0.57
5. Does the belt spread? — testing the spillover reading of $\rho$¶
The cross-sectional notebook found a large spatial multiplier and was careful not to read it as diffusion: on one snapshot, an omitted regional factor and genuine county-to-county spread look identical. A panel separates them. If receipt propagated across county lines the belt would be expanding at its edges, so the counties ringing the core — exposed to it but not part of it — should carry the largest positive slope deviations. ST.CARlinear supplies posterior draws of $\delta_i$, so the test runs here independently of the Python engine and on a different spatial prior (Leroux rather than intrinsic CAR).
core_m <- phi > quantile(phi, 0.90) # the persistent belt core
ring_m <- (as.vector(W %*% as.numeric(core_m)) > 0) & !core_m # touches the core, not in it
rest_m <- !core_m & !ring_m
D <- m$samples$delta # posterior draws of the slope deviation
grp <- list(`belt core`=core_m, `ring adjacent to the core`=ring_m, `rest of the country`=rest_m)
cat(sprintf("%-28s%9s%14s%13s\n", "", "counties", "mean slope", "P(slope>0)"))
for (nm in names(grp)) {
mk <- grp[[nm]]; dd <- rowMeans(D[, mk, drop=FALSE])
cat(sprintf("%-28s%9d%+14.3f%13.2f\n", nm, sum(mk), mean(dd), mean(dd > 0)))
}
contrast <- rowMeans(D[, ring_m, drop=FALSE]) - rowMeans(D[, rest_m, drop=FALSE])
cat(sprintf("\nring minus rest: %+.3f P(ring grew faster than the rest) = %.2f\n",
mean(contrast), mean(contrast > 0)))
cat("\nThe ring is not the fastest-growing group in the country -- it is slower than the rest, which is\n")
cat("the opposite of what diffusion predicts. Two decades is ample time for a contagious process to show\n")
cat("itself at the boundary; what shows up there is retreat. This reproduces the Python notebook on a\n")
cat("different spatial prior, so the conclusion is not an artefact of either model: read the spatial\n")
cat("multiplier as how the association is DISTRIBUTED across space, not as a mechanism by which receipt\n")
cat("travels. The belt is a place where disadvantage persists, not a front that is advancing.\n")
counties mean slope P(slope>0)
belt core 304 -0.078 0.00 ring adjacent to the core 395 -0.022 0.00 rest of the country 2336 +0.014 1.00
ring minus rest: -0.036 P(ring grew faster than the rest) = 0.00
The ring is not the fastest-growing group in the country -- it is slower than the rest, which is
the opposite of what diffusion predicts. Two decades is ample time for a contagious process to show
itself at the boundary; what shows up there is retreat. This reproduces the Python notebook on a
different spatial prior, so the conclusion is not an artefact of either model: read the spatial
multiplier as how the association is DISTRIBUTED across space, not as a mechanism by which receipt
travels. The belt is a place where disadvantage persists, not a front that is advancing.
6. Summary¶
CARBayesST::ST.CARlinear reproduces the from-scratch result on the identical panel: the disability belt is entrenched, not expanding. The persistent spatial level $\phi_i$ (the coalfield core — Buchanan, Dickenson, Norton VA; Floyd, Pike KY; Mingo WV) dominates the county-specific slope $\delta_i$, the core's trend is flat-to-slightly-negative, and the national rolls' rise-to-2014-peak-and-fall lifted the whole country without redrawing the map. The year-specific age-standardization (three Census vintages × the yearly SSA age distribution) removes the baby-boomer ageing wave that would otherwise fake a growing belt. SSI (§4) shares the same persistent geography but a different time signature — flatter than DI, gently declining rather than boom-bust — and is likewise entrenched. This is the package mirror of belt_spacetime_python.ipynb and the temporal completion of the disability-belt application — both programs mapped across space and time. Extension: ST.CARanova for non-linear county trajectories.