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