Areal Spatial Modelling — Weights, Autocorrelation & the CAR Family

Regions, Not Points

Areal data attach a measurement to a region — a county, tract or state — rather than to a point. Regions have no natural coordinates, only a neighbour structure, encoded in a weights matrix $W$ with $w_{ij}=1$ when two regions share a border. Regions that touch tend to be alike, and ignoring that both wastes information and, for rare events in small areas, leaves the raw rates dangerously noisy: North Carolina's SIDS rates run from a standardised ratio of 0 to 4.7 across 100 counties and 667 deaths, most of that spread being sample size rather than risk.

Is there spatial structure?

Moran's I is the spatial analogue of a correlation coefficient, and a permutation test says whether it beats chance. One wrinkle is worth naming because it looks like disagreement and is not: the statistic depends on how $W$ is normalised. On the binary adjacency it is 0.21; row-standardised — the spdep default, and the slope of the Moran scatterplot — it is 0.23. Both are Moran's I, both give $p\le0.002$ against 999 reshuffles, and the row-standardised value is what the R notebook reports. The structure is real, and it is what the CAR model exploits.

The BYM disease-mapping model

The workhorse is the Besag–York–Mollié convolution: a Poisson count with an offset for the expected count, a spatially structured effect $\phi$ carrying a conditional autoregressive prior, and an unstructured effect $\theta$ for everything local. Each region's $\phi_i$ is pulled toward its neighbours' mean with a precision proportional to its neighbour count, so evidence is borrowed exactly along the adjacency graph. The sampler is Metropolis-within-Gibbs — conjugate inverse-gamma draws for the variances, Metropolis steps for the rest, with every random-walk scale tuned during burn-in.

On the NC-SIDS data the covariate effect is clear — a higher nonwhite-births share goes with higher SIDS risk — and three independent engines agree on it: +0.42 from the from-scratch sampler, +0.41 from PyMC's ICAR, and +0.40 [0.28, 0.54] from R's CARBayes. The relative-risk maps agree to a correlation of 0.996. The smoothing is the point: the raw standardised ratio has a spread of 0.78, the smoothed relative risk 0.44, with the extremes pulled toward 1 — the small-county spikes tamed while the broad high-risk band survives.

How strong is the dependence, and where is risk elevated?

Fixing $\rho=1$ gives the intrinsic CAR, which borrows the whole neighbourhood mean and is the BYM default. Letting it go free gives a proper CAR that also estimates how strong the dependence is — a quantity the intrinsic version simply assumes. And because a smoothed map is only useful if it supports a decision, the last step converts the posterior into an exceedance probability: R flags the 7 counties with $P(\text{RR}>1.5)>0.9$ — elevated with confidence, rather than merely elevated in a point estimate.

Notebooks

Downloads

Areal Module — Source Code

"""
areal.py -- AREAL / LATTICE spatial models: weights, autocorrelation, and the CAR family (scratch).

Backs the notebooks in  "Areal Spatial Modelling -- Weights, Autocorrelation & the CAR Family".

Areal data are measurements attached to REGIONS -- counts or rates for counties, tracts, states.
The regions have no natural coordinates, only a NEIGHBOUR structure, encoded in a spatial-weights
matrix W (w_ij = 1 if regions i and j share a border). Two regions that touch tend to be alike:
SPATIAL AUTOCORRELATION. Ignoring it wastes information and, for rare events in small regions,
leaves the raw rates dangerously noisy.

MORAN'S I measures the autocorrelation -- the spatial analogue of a correlation coefficient -- and a
permutation test says whether it is more than chance.

The workhorse model is the CONDITIONAL AUTOREGRESSIVE (CAR) prior on region-level random effects.
For disease/event mapping the counts are Poisson,

    y_i ~ Poisson(E_i * exp(alpha + x_i' beta + phi_i + theta_i)),

with E_i an expected count (offset), phi a SPATIALLY STRUCTURED effect and theta an unstructured one
-- the Besag-York-Mollie (BYM) convolution. phi is given a CAR prior whose conditional is

    phi_i | phi_{-i} ~ N( rho * (sum_j w_ij phi_j) / n_i ,  tau^2 / n_i ),

where n_i is the number of neighbours. rho = 1 is the INTRINSIC CAR (the BYM default -- borrows the
whole neighbourhood mean); rho < 1 is a PROPER CAR that also estimates the STRENGTH of spatial
dependence. The result is a smoothed map of RELATIVE RISK that pulls unreliable small-area rates
toward their neighbours -- exactly what is wanted for regional disability/mortality incidence.

The sampler is Metropolis-within-Gibbs: the variance components are conjugate (inverse-gamma), while
the non-conjugate Poisson forces Metropolis steps for alpha, beta and the random effects (with the
CAR conditional as an efficient proposal). Every random-walk scale is tuned during burn-in toward a
0.44 acceptance rate with a diminishing step -- the informative Poisson likelihood puts the posterior
sd of a coefficient two orders of magnitude below a plausible-looking fixed scale, and an untuned
chain simply sits still. This is the areal counterpart of the disease-mapping work
in "Bayesian Hierarchical Spatial Poisson Model" (Scotland lip cancer), extended here with a
formal autocorrelation test and an estimated spatial-dependence parameter.
"""

import numpy as np


# --------------------------------------------------------------------------- #
#  spatial autocorrelation                                                      #
# --------------------------------------------------------------------------- #

def morans_i(y, W):
    """Moran's I for values y on a binary weights matrix W."""
    y = np.asarray(y, float); z = y - y.mean(); n = len(y); S0 = W.sum()
    return (n / S0) * (z @ W @ z) / (z @ z)


def morans_test(y, W, rng, nperm=999):
    """permutation test for Moran's I: returns (I, p-value, permutation-null draws)."""
    I = morans_i(y, W); n = len(y)
    null = np.array([morans_i(y[rng.permutation(n)], W) for _ in range(nperm)])
    p = (1 + (null >= I).sum()) / (nperm + 1)
    return I, p, null


# --------------------------------------------------------------------------- #
#  BYM / CAR Poisson disease-mapping sampler                                     #
# --------------------------------------------------------------------------- #

def bym_gibbs(y, E, X, W, rng, draws=4000, burn=4000, rho=1.0, sample_rho=False,
              a=1.0, b=0.01, beta_sd=5.0):
    """Poisson BYM: y_i ~ Poisson(E_i exp(alpha + x_i'beta + phi_i + theta_i)).
      phi : CAR spatial effect (rho=1 intrinsic; sample_rho=True estimates a proper-CAR rho in (0,1))
      theta : unstructured N(0, sig2) effect
    Returns posterior draws of alpha, beta, the spatial SD and unstructured SD, rho, and the
    relative risks RR_i = exp(alpha + x_i'beta + phi_i + theta_i)."""
    y = np.asarray(y, float); E = np.asarray(E, float); X = np.atleast_2d(np.asarray(X, float))
    if X.shape[0] != len(y):
        X = X.T
    n, p = X.shape; nnb = W.sum(1)                                    # neighbour counts
    alpha = np.log((y.sum() + 1) / (E.sum() + 1)); beta = np.zeros(p)
    phi = np.zeros(n); theta = np.zeros(n); tau2 = 0.5; sig2 = 0.1
    sphi = 0.4 * np.ones(n); sth = 0.4 * np.ones(n); sal = 0.1; sbe = 0.1 * np.ones(p)
    A = np.empty(draws); B = np.empty((draws, p)); TAU = np.empty(draws); SIG = np.empty(draws)
    RHO = np.empty(draws); RR = np.empty((draws, n))
    def loglam(al, be, ph, th):
        return np.log(E) + al + X @ be + ph + th
    acc = np.zeros(n); acc_th = np.zeros(n); acc_al = 0.0; acc_be = np.zeros(p); nadapt = 0
    for it in range(draws + burn):
        eta = loglam(alpha, beta, phi, theta); lam = np.exp(eta)
        # ---- alpha (RW-Metropolis) ----
        al_p = alpha + sal * rng.standard_normal(); lam_p = np.exp(loglam(al_p, beta, phi, theta))
        if np.log(rng.random()) < (y @ (np.log(lam_p) - np.log(lam)) - (lam_p - lam).sum() - (al_p**2 - alpha**2) / (2 * 100)):
            alpha = al_p; lam = lam_p; acc_al += 1
        # ---- beta (componentwise RW-Metropolis) ----
        for j in range(p):
            be_p = beta.copy(); be_p[j] += sbe[j] * rng.standard_normal()
            lam_p = np.exp(loglam(alpha, be_p, phi, theta))
            if np.log(rng.random()) < (y @ (np.log(lam_p) - np.log(lam)) - (lam_p - lam).sum() - (be_p[j]**2 - beta[j]**2) / (2 * beta_sd**2)):
                beta = be_p; lam = lam_p; acc_be[j] += 1
        # ---- phi (Metropolis with CAR-conditional prior mean) ----
        for i in range(n):
            nbmean = (W[i] @ phi) / max(nnb[i], 1.0)                   # neighbours as they stand NOW
            phi_p = phi[i] + sphi[i] * rng.standard_normal()
            dloglik = y[i] * (phi_p - phi[i]) - (np.exp(eta[i] + phi_p - phi[i]) - lam[i])
            dprior = -(nnb[i] / (2 * tau2)) * ((phi_p - rho * nbmean)**2 - (phi[i] - rho * nbmean)**2)
            if np.log(rng.random()) < dloglik + dprior:
                eta[i] += phi_p - phi[i]; phi[i] = phi_p; lam[i] = np.exp(eta[i]); acc[i] += 1
        phi -= phi.mean()                                             # sum-to-zero identification
        eta = loglam(alpha, beta, phi, theta); lam = np.exp(eta)
        # ---- theta (Metropolis, N(0,sig2) prior) ----
        for i in range(n):
            th_p = theta[i] + sth[i] * rng.standard_normal()
            dloglik = y[i] * (th_p - theta[i]) - (np.exp(eta[i] + th_p - theta[i]) - lam[i])
            dprior = -(th_p**2 - theta[i]**2) / (2 * sig2)
            if np.log(rng.random()) < dloglik + dprior:
                eta[i] += th_p - theta[i]; theta[i] = th_p; lam[i] = np.exp(eta[i]); acc_th[i] += 1
        # ---- variance components (Gibbs) ----
        quad = 0.5 * (nnb * phi * phi).sum() - 0.5 * rho * (phi * (W @ phi)).sum()
        tau2 = 1.0 / rng.gamma(a + (n - 1) / 2, 1.0 / (b + max(quad, 1e-6)))
        sig2 = 1.0 / rng.gamma(a + n / 2, 1.0 / (b + 0.5 * (theta @ theta)))
        # ---- rho (proper CAR, Metropolis on logit) ----
        if sample_rho:
            rho_p = np.clip(rho + 0.05 * rng.standard_normal(), 1e-3, 0.999)
            def carlp(r):
                Q = np.diag(nnb) - r * W
                sgn, ld = np.linalg.slogdet(Q)
                return 0.5 * ld - (0.5 / tau2) * (phi @ Q @ phi)
            if np.log(rng.random()) < carlp(rho_p) - carlp(rho):
                rho = rho_p
        if it < burn and it % 100 == 99:                              # adapt every scale during burn-in
            nadapt += 1; g = 1.0 / np.sqrt(nadapt)                     # diminishing step, so it settles
            sphi *= np.exp((acc / 100 - 0.44) * g); sth *= np.exp((acc_th / 100 - 0.44) * g)
            sal *= np.exp((acc_al / 100 - 0.44) * g); sbe *= np.exp((acc_be / 100 - 0.44) * g)
            acc[:] = 0; acc_th[:] = 0; acc_al = 0.0; acc_be[:] = 0
        if it >= burn:
            k = it - burn; A[k] = alpha; B[k] = beta; TAU[k] = np.sqrt(tau2); SIG[k] = np.sqrt(sig2)
            RHO[k] = rho; RR[k] = np.exp(alpha + X @ beta + phi + theta)
    return dict(alpha=A, beta=B, spatial_sd=TAU, unstruct_sd=SIG, rho=RHO, RR=RR)


def expected_counts(y, pop):
    """internally-standardised expected counts: E_i = pop_i * (sum y / sum pop)."""
    y = np.asarray(y, float); pop = np.asarray(pop, float)
    return pop * (y.sum() / pop.sum())

References

• Besag, J., York, J. & Mollié, A. (1991). Bayesian image restoration, with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics 43(1), 1–20. — the BYM convolution model
• Besag, J. (1974). Spatial interaction and the statistical analysis of lattice systems. JRSS-B 36(2), 192–236. — conditional autoregression itself
• Moran, P. A. P. (1950). Notes on continuous stochastic phenomena. Biometrika 37(1/2), 17–23. — the autocorrelation statistic
• Banerjee, S., Carlin, B. P. & Gelfand, A. E. (2014). Hierarchical Modeling and Analysis for Spatial Data, 2nd ed. CRC Press. — the standard reference for the whole areal chapter
• Cressie, N. & Chan, N. H. (1989). Spatial modeling of regional variables. JASA 84(406), 393–401. — the North Carolina SIDS analysis
• Lee, D. (2013). CARBayes: an R package for Bayesian spatial modeling with conditional autoregressive priors. Journal of Statistical Software 55(13). — the package used for the R cross-check