""" 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())=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