Areal Spatial Modelling — Weights, Autocorrelation & the CAR Family

Smoothing noisy regional rates by borrowing from the neighbours

This opens a spatial-modelling arc. Areal (or lattice) data are measurements attached to regions — counties, tracts, states — that have no coordinates, only a neighbour structure. Two regions that share a border tend to be alike (spatial autocorrelation), and ignoring that both wastes information and, for rare events in small regions, leaves the raw rates dangerously noisy. This is precisely the situation for regional disability, mortality and migration incidence: an SSA analyst comparing disability rates across states or counties faces small counts, unstable rates, and spatial structure — the machinery here is built for exactly that.

We use the classic North Carolina SIDS data — sudden infant deaths across the 100 NC counties, 1974–78 — as a stand-in for any regional event-count problem (deaths per birth here; disability awards per capita there). The tools: a spatial-weights matrix, Moran's I for autocorrelation, and the conditional autoregressive (CAR) disease-mapping model that produces a smoothed relative-risk map, pulling unreliable small-county rates toward their neighbours. Maps throughout — that is the point of spatial work. This builds directly on the Besag–York–Mollié model in Bayesian Hierarchical Spatial Poisson Model (Scotland lip cancer), adding a formal autocorrelation test and an estimated spatial-dependence parameter.

import numpy as np, pandas as pd, matplotlib.pyplot as plt, geopandas as gpd
import areal as A
rng = np.random.default_rng(1)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; GREY="#718096"
d = pd.read_csv("nc_sids.csv"); W = pd.read_csv("nc_adj.csv").to_numpy(float)
gdf = gpd.read_file("nc_counties.geojson").merge(d[["fips"]].assign(order=range(len(d))),
      left_on="FIPSNO", right_on="fips").sort_values("order").reset_index(drop=True)  # align map to data order
y = d["sids"].to_numpy(float); births = d["births"].to_numpy(float)
E = A.expected_counts(y, births); SMR = y / E
def choropleth(ax, values, title, cmap="OrRd", vcenter=None):
    g = gdf.copy(); g["v"] = values
    g.plot(column="v", ax=ax, cmap=cmap, edgecolor="white", linewidth=0.3, legend=True,
           legend_kwds={"shrink":0.6}); ax.set_title(title, fontsize=10); ax.axis("off")
print(f"{len(d)} NC counties; {int(y.sum())} SIDS deaths over {int(births.sum()):,} births ({1000*y.sum()/births.sum():.2f} per 1000)")
print("Areal data: values on regions with a neighbour structure. The task is a reliable map of the underlying risk.")

100 NC counties; 667 SIDS deaths over 329,962 births (2.02 per 1000)
Areal data: values on regions with a neighbour structure. The task is a reliable map of the underlying risk.

1. Expected counts, the SMR, and why raw rates mislead

Raw death counts are not comparable across counties — a large county has more SIDS deaths simply because it has more births. To measure risk we need an expected count $E_i$: how many deaths county $i$ would have at a common reference rate. Using internal standardisation, that reference is the statewide rate $\bar r=\sum_j y_j/\sum_j\text{births}_j\approx 2.0$ per 1000, so $$E_i=\text{births}_i\times\bar r .$$ Births — not residents — are the *population at risk*** (a SIDS death can only befall a birth); picking the right denominator (births, person-years, insured workers…) is the modeller's job. By construction $\sum_i E_i=\sum_i y_i$: expected and observed totals match. (More generally, with age/sex strata $k$ one uses **indirect standardisation, $E_i=\sum_k \text{pop}_{ik}\,r_k^{\text{ref}}$, so a region is not flagged high merely for having an older population.)

The standardised mortality ratio $\text{SMR}_i=y_i/E_i$ is then the raw relative-risk estimate: $\text{SMR}=1$ is average risk, $2$ is twice expected. It is the naive risk map (below). But a county with few births swings wildly on one or two extra deaths, so the raw SMR is dominated by small-area noise, not real risk — as the scatter of SMR against county size then makes plain. That is the problem the spatial model fixes.

fig,ax=plt.subplots(figsize=(7.5,3.4))
choropleth(ax, SMR, "Raw SMR (observed / expected SIDS deaths)", cmap="OrRd")
plt.tight_layout(); plt.show()
print(f"raw SMR ranges {SMR.min():.1f} to {SMR.max():.1f}; the darkest counties are the alarming-looking extremes.")

raw SMR ranges 0.0 to 4.7; the darkest counties are the alarming-looking extremes.

fig,ax=plt.subplots(figsize=(5.6,3))
ax.scatter(births, SMR, s=22, color=BLUE, alpha=.6); ax.axhline(1,color="k",lw=.8,label="expected (SMR=1)")
ax.set_xlabel("births in county  (larger = more reliable)"); ax.set_ylabel("raw SMR"); ax.set_title("Small counties swing wildly")
ax.annotate("few births -> extreme,\nunreliable SMR", xy=(births.min(),SMR.max()), xytext=(0.30*births.max(),0.88*SMR.max()),
            fontsize=9, ha="left", arrowprops=dict(arrowstyle="->",color=GREY)); ax.legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
print("But those extremes are all SMALL counties, where one extra death moves the rate enormously -- the dark spots on")
print("the map are mostly noise, not real risk. A credible map must borrow strength across space (the neighbour graph).")

But those extremes are all SMALL counties, where one extra death moves the rate enormously -- the dark spots on
the map are mostly noise, not real risk. A credible map must borrow strength across space (the neighbour graph).

2. Is there spatial structure? — Moran's I

Before modelling, test whether neighbouring counties really are alike. Moran's I is the spatial correlation coefficient: $I=\frac{n}{S_0}\frac{\sum_{ij}w_{ij}z_iz_j}{\sum_i z_i^2}$ with $z_i=y_i-\bar y$. A permutation test — reshuffling the values over the map — gives its null distribution. The Moran scatterplot plots each county's value against the average of its neighbours; an upward slope is positive autocorrelation.

I,pval,null = A.morans_test(SMR, W, rng, nperm=999)                   # binary adjacency
Wr = W / W.sum(1, keepdims=True); lag = Wr @ (SMR - SMR.mean())        # spatially-lagged (row-standardised)
Ir,pvalr,_ = A.morans_test(SMR, Wr, rng, nperm=999)                    # row-standardised, as spdep uses
fig,ax=plt.subplots(1,2,figsize=(9,3))
ax[0].hist(null,bins=30,color=GREY,alpha=.7,density=True); ax[0].axvline(I,color=RED,lw=2,label=f"observed I={I:.2f}")
ax[0].axvline(null.mean(),color="k",ls="--",lw=1,label="null mean"); ax[0].set_xlabel("Moran's I"); ax[0].set_title(f"Permutation test (p={pval:.3f})"); ax[0].legend(frameon=False,fontsize=8)
ax[1].scatter(SMR-SMR.mean(), lag, s=18, color=BLUE, alpha=.6)
b=np.polyfit(SMR-SMR.mean(), lag, 1)[0]; xx=np.linspace((SMR-SMR.mean()).min(),(SMR-SMR.mean()).max(),50)
ax[1].plot(xx,b*xx,color=RED,lw=2,label=f"slope={b:.2f}"); ax[1].axhline(0,color="k",lw=.5); ax[1].axvline(0,color="k",lw=.5)
ax[1].set_xlabel("county SMR (centred)"); ax[1].set_ylabel("mean of neighbours' SMR"); ax[1].set_title("Moran scatterplot"); ax[1].legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
print(f"Moran's I = {I:.2f} (binary adjacency, p = {pval:.3f})   {Ir:.2f} (row-standardised, p = {pvalr:.3f})")
print(f"Both are Moran's I; they differ only in how W is normalised, and the scatterplot slope ({b:.2f}) is the")
print("row-standardised version -- which is what spdep reports in the R notebook. SIDS rates are significantly")
print("autocorrelated either way: neighbouring counties resemble each other far more than a random map would.")

Moran's I = 0.21 (binary adjacency, p = 0.002)   0.23 (row-standardised, p = 0.001)
Both are Moran's I; they differ only in how W is normalised, and the scatterplot slope (0.23) is the
row-standardised version -- which is what spdep reports in the R notebook. SIDS rates are significantly
autocorrelated either way: neighbouring counties resemble each other far more than a random map would.

3. The BYM disease-mapping model — a smoothed risk map

The Besag–York–Mollié model treats the counts as $y_i\sim\text{Poisson}(E_i\,\theta_i)$, with the log relative risk $\log\theta_i$ built from a covariate, a spatially structured CAR effect $\phi$ (borrowing from neighbours) and an unstructured effect $\theta^{\text u}$. The expected count enters as an offset: $$\log\mu_i=\underbrace{\log E_i}_{\text{coefficient fixed at }1}+\ \alpha+x_i\beta+\phi_i+\theta^{\text u}_i .$$ An offset is a known exposure whose coefficient is pinned at 1 rather than estimated — and pinning it at 1 is exactly what makes $\text{RR}_i=\exp(\alpha+x_i\beta+\phi_i+\theta^{\text u}_i)$ a relative risk, a multiplier on the expected count, so the whole map reads relative to baseline. We include the county's proportion of nonwhite births as the covariate; the fitted RR is the smoothed map — small-county noise pulled toward the local level.

The covariate, and what each parameter means

Why nonwhite births? SIDS incidence in the US has historically run about twice as high among nonwhite infants, so a county's proportion of nonwhite births ($\text{NWBIR74}/\text{BIR74}$, standardised) is the strongest available ecological predictor of its SIDS rate — and the only real covariate this dataset carries. Including it lets $x_i\beta$ absorb the compositional part of the variation (nonwhite-birth share is itself spatially clustered), so the spatial effect $\phi$ is left to capture residual clustering beyond demographics. Two caveats: it is an ecological covariate — county-level, not individual, so a positive $\beta$ reflects correlated socioeconomic factors (healthcare access, poverty, sleeping environment), not a causal race effect, and reading it at the individual level would be the ecological fallacy — and it is associational, not causal. For regional disability incidence the analogues would be age structure, industry mix and poverty rate.

The parameters.

symbol	role	how to read it
$\alpha$	intercept	baseline log-risk; $\exp(\alpha)\approx$ overall relative risk $\approx 1$
$\beta$	covariate effect	per 1 SD of nonwhite-birth share; $\exp(\beta)$ is the relative-risk multiplier
$\phi_i$	structured (CAR) effect	each county's risk deviation borrowed from neighbours; sum-to-zero
$\theta^{\text u}_i$	unstructured effect	independent county heterogeneity, $\sim N(0,\sigma^2)$
$\tau$	spatial SD	how large the spatial variation is
$\sigma$	unstructured SD	how large the non-spatial noise is
$\rho$	dependence (proper-CAR fit)	strength of neighbour borrowing, 0 to 1

Two readings carry the interpretation: $\exp(\beta)$ turns the coefficient into a risk multiplier (e.g. $\beta=0.4\Rightarrow$ a 1-SD-higher county has $\approx 49\%$ higher risk), and the ratio of $\tau$ to $\sigma$ is a variance partition — $\tau\gg\sigma$ means the leftover variation is overwhelmingly spatial, which is exactly what justifies the CAR smoothing.

x = ((d["nonwhite"]/d["births"]).to_numpy()); x=(x-x.mean())/x.std()
res = A.bym_gibbs(y, E, x, W, rng, draws=2500, burn=2500, rho=1.0)
RR = res["RR"].mean(0); be=res["beta"].mean(0)[0]; bl,bh=np.percentile(res["beta"],[2.5,97.5])
print(f"covariate (nonwhite births): beta = {be:+.2f}  95% CrI [{bl:+.2f}, {bh:+.2f}]  (higher nonwhite share -> higher SIDS risk)")
print(f"variance split: spatial SD {res['spatial_sd'].mean():.2f}  vs  unstructured SD {res['unstruct_sd'].mean():.2f}  -> structure dominates")
fig,ax=plt.subplots(1,2,figsize=(9.8,3))
vmax=max(SMR.max(),RR.max())
choropleth(ax[0], SMR, "Raw SMR (noisy)", cmap="OrRd")
choropleth(ax[1], RR, "BYM smoothed relative risk", cmap="OrRd")
plt.tight_layout(); plt.show()
fig,ax=plt.subplots(figsize=(5.2,2.7))
ax.scatter(SMR, RR, s=22, color=BLUE, alpha=.6); lim=[0,SMR.max()*1.05]; ax.plot(lim,lim,"k--",lw=1)
ax.axhline(1,color=GREY,lw=.6); ax.set_xlabel("raw SMR"); ax.set_ylabel("BYM relative risk")
ax.set_title(f"Shrinkage: raw sd {SMR.std():.2f} -> smoothed sd {RR.std():.2f}"); plt.tight_layout(); plt.show()
print("The smoothed map keeps the broad high-risk band but tames the small-county spikes: points sit below the")
print("45-degree line at the extremes (pulled toward 1) and the map is now a credible picture of spatial risk --")
print("the same shrinkage logic as the Scotland lip-cancer BYM, here with an estimated covariate effect.")

covariate (nonwhite births): beta = +0.43  95% CrI [+0.29, +0.55]  (higher nonwhite share -> higher SIDS risk)
variance split: spatial SD 0.34  vs  unstructured SD 0.14  -> structure dominates

The smoothed map keeps the broad high-risk band but tames the small-county spikes: points sit below the
45-degree line at the extremes (pulled toward 1) and the map is now a credible picture of spatial risk --
the same shrinkage logic as the Scotland lip-cancer BYM, here with an estimated covariate effect.

4. How strong is the spatial dependence? — proper CAR

The BYM above uses the intrinsic CAR ($\rho=1$): every county borrows the full neighbourhood mean. A proper CAR instead estimates $\rho\in(0,1)$, the strength of spatial dependence — $\rho\to0$ means the map is essentially unstructured, $\rho\to1$ the strongly-smoothed intrinsic limit. Estimating it puts a number on how spatial the process is.

resp = A.bym_gibbs(y, E, x, W, rng, draws=2000, burn=2500, sample_rho=True)
rho=resp["rho"]; rl,rh=np.percentile(rho,[2.5,97.5])
fig,ax=plt.subplots(figsize=(5.2,2.5)); ax.hist(rho,bins=40,color=GREEN,alpha=.8,density=True)
ax.axvline(rho.mean(),color=RED,lw=2,label=f"mean {rho.mean():.2f}"); ax.set_xlabel(r"spatial-dependence parameter $\rho$")
ax.set_title(f"Proper-CAR spatial dependence: rho = {rho.mean():.2f}  95% CrI [{rl:.2f}, {rh:.2f}]"); ax.legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
print(f"rho = {rho.mean():.2f} [{rl:.2f}, {rh:.2f}]: moderate spatial dependence, though with only 100 counties it is")
print("imprecisely estimated (a common feature -- the intrinsic rho=1 default is often used precisely because rho is")
print("hard to pin down). Either way the covariate effect and the broad risk pattern are stable.")

rho = 0.47 [0.00, 0.92]: moderate spatial dependence, though with only 100 counties it is
imprecisely estimated (a common feature -- the intrinsic rho=1 default is often used precisely because rho is
hard to pin down). Either way the covariate effect and the broad risk pattern are stable.

5. A policy map — where is risk elevated with confidence?

A point estimate of relative risk hides its uncertainty. The exceedance probability $P(\text{RR}_i>1.5\mid\text{data})$ — computed straight from the posterior draws — flags the counties the model is confident are high-risk, not merely those with a high point estimate. This is the map an analyst actually acts on: which regions warrant attention.

exc = (res["RR"] > 1.5).mean(0)
fig,ax=plt.subplots(1,2,figsize=(9.8,3))
choropleth(ax[0], RR, "Relative risk (posterior mean)", cmap="OrRd")
choropleth(ax[1], exc, r"P(relative risk > 1.5)  — exceedance", cmap="Reds")
plt.tight_layout(); plt.show()
nhi=(exc>0.9).sum()
print(f"{nhi} counties have P(RR>1.5) above 0.9 -- flagged as elevated with high posterior confidence, distinct from")
print("counties whose point estimate is high only because of small-sample noise. Exceedance probabilities turn the")
print("smoothed map into a decision tool -- exactly what regional disability/mortality surveillance needs.")

6 counties have P(RR>1.5) above 0.9 -- flagged as elevated with high posterior confidence, distinct from
counties whose point estimate is high only because of small-sample noise. Exceedance probabilities turn the
smoothed map into a decision tool -- exactly what regional disability/mortality surveillance needs.

6. Cross-check in PyMC

PyMC has a built-in intrinsic-CAR distribution (pm.ICAR). We fit the same BYM model and confirm the covariate effect and the smoothed risks match the from-scratch sampler.

import pymc as pm
adj=[np.where(W[i]>0)[0].tolist() for i in range(len(y))]
with pm.Model() as mod:
    a0=pm.Normal("a0",0,5); b=pm.Normal("b",0,5); tau=pm.HalfNormal("tau",2); sig=pm.HalfNormal("sig",2)
    phi=pm.ICAR("phi", W=W.astype(int)); th=pm.Normal("th",0,1,shape=len(y))
    eta=np.log(E)+a0+b*x+tau*phi+sig*th
    pm.Poisson("y", mu=pm.math.exp(eta), observed=y)
    idata=pm.sample(3000,tune=3000,chains=4,target_accept=0.97,random_seed=3,progressbar=False)
b_pm=float(idata.posterior["b"].mean())
RR_pm=np.exp(idata.posterior["a0"].values.ravel()[:,None]+idata.posterior["b"].values.ravel()[:,None]*x[None,:]
             +(idata.posterior["tau"].values.ravel()[:,None]*idata.posterior["phi"].stack(s=("chain","draw")).values.T)
             +(idata.posterior["sig"].values.ravel()[:,None]*idata.posterior["th"].stack(s=("chain","draw")).values.T)).mean(0)
import arviz as az
V=["a0","b","tau","sig","phi","th"]
maxrhat=max(float(v.max()) for v in az.rhat(idata,var_names=V).data_vars.values())
miness=min(float(v.min()) for v in az.ess(idata,var_names=V).data_vars.values())
print(f"covariate beta:  from-scratch {be:+.2f}   PyMC {b_pm:+.2f}")
print(f"relative-risk agreement: corr {np.corrcoef(RR,RR_pm)[0,1]:.3f}")
print(f"PyMC convergence over all {sum(int(np.prod(v.shape)) for v in az.rhat(idata,var_names=V).data_vars.values())} parameters: max r-hat {maxrhat:.3f}, min ESS {miness:.0f}")
print("PyMC's ICAR reproduces the from-scratch BYM: the smoothed map and the covariate effect agree.")

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

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [a0, b, tau, sig, phi, th]

Sampling 4 chains for 3_000 tune and 3_000 draw iterations (12_000 + 12_000 draws total) took 55 seconds.

covariate beta:  from-scratch +0.43   PyMC +0.41
relative-risk agreement: corr 0.995

PyMC convergence over all 204 parameters: max r-hat 1.005, min ESS 844
PyMC's ICAR reproduces the from-scratch BYM: the smoothed map and the covariate effect agree.

7. Summary

Areal data live on a graph of regions, and their raw rates are noisy where regions are small. The workflow here is the backbone of spatial epidemiology and regional analysis: build the spatial-weights matrix, confirm structure with Moran's I (the NC SIDS rates were significantly autocorrelated), and fit a CAR/BYM model that borrows strength across neighbours to produce a smoothed relative-risk map — with a covariate effect (higher nonwhite-birth share, higher SIDS risk), an estimated spatial-dependence $\rho$, and exceedance probabilities that flag genuinely elevated regions. Every step ended in a map, because the map is the deliverable. A PyMC ICAR fit reproduced it.

This is the areal foundation for the arc, and the direct tool for regional disability, mortality and migration incidence in SSA-style work — small-area counts, spatial smoothing, and confident flagging of high-risk regions. It extends the Besag–York–Mollié model of Bayesian Hierarchical Spatial Poisson Model (Scotland) with autocorrelation testing and a proper-CAR $\rho$. One deliberate distinction sets up the next subsection: the CAR here is a hierarchical random effect (smoothing), whereas the spatial-econometrics models to come use the same weights matrix in a simultaneous autoregressive lag of the outcome itself (spillovers) — same neighbour graph, different philosophy. The remaining subsections — geostatistics/kriging, point processes, and a spatiotemporal capstone — complete the spatial toolkit.