""" spacetime.py -- SPATIOTEMPORAL areal models: the space-time CAR (from scratch). Backs the notebooks in "Spatiotemporal Modelling -- the Space-Time CAR". The capstone of the spatial arc: data indexed by BOTH region and time -- disease counts, mortality, disability receipt, unemployment across areas over years. It unifies the areal CAR models (space) with the time-series work (time). The workhorse is the space-time ANOVA / Knorr-Held decomposition of the log relative risk: y_it ~ Poisson( E_it * exp(alpha + x_i' beta + phi_i + gamma_t) ) phi_i : SPATIAL main effect -- a CAR/ICAR field (which regions run high, averaged over time) gamma_t : TEMPORAL main effect -- a random walk (the shared trend across all regions) st_gibbs fits this ADDITIVE form -- the two main effects, no interaction. The full Knorr-Held ANOVA adds a third component delta_it, a space-time INTERACTION letting individual regions depart from the shared temporal wave; CARBayesST reaches it with interaction=TRUE, and ST.CARar is the autoregressive alternative. The additive model is the right starting point because the interaction is the least identified part of the decomposition, and on data generated without one it is simply noise to fit. phi is the areal-CAR machinery of the earlier notebooks; gamma is a first-order random walk in time (a Bayesian smooth trend, the temporal analogue of the CAR); together they are a hierarchical model in space AND time. Fitting is Metropolis-within-Gibbs: phi_i is updated against the Poisson likelihood summed over all TIMES for region i (with the ICAR neighbour prior), gamma_t against the likelihood summed over all REGIONS at time t (with the random-walk prior), and the two smoothing variances are conjugate inverse-gamma. This is the model CARBayesST fits (ST.CARanova / ST.CARar) and the space-time extension of everything in the arc. """ import numpy as np def smooth_field(W, rng, k=8, scale=0.6): """a spatially smooth field: Gaussian noise repeatedly averaged with neighbours over the graph.""" n = W.shape[0]; f = rng.standard_normal(n); nnb = np.maximum(W.sum(1), 1) for _ in range(k): f = 0.5 * f + 0.5 * (W @ f) / nnb f = (f - f.mean()) / f.std() return scale * f def simulate_st(pop, W, rng, T=12, alpha=-1.0, beta=0.3, spatial_scale=0.6, temporal_amp=0.7): """simulate a space-time Poisson process on real geography. Temporal effect is an epidemic-style rise-and-fall; spatial effect a smooth field; covariate = standardised log-population. Returns y (N,T), offset E (N,), covariate x (N,), and the true (phi, gamma, beta).""" n = len(pop); E = np.asarray(pop, float) / 1000.0 x = rng.standard_normal(n) # a county covariate, NOT spatially structured (avoids confounding with phi) phi = smooth_field(W, rng, scale=spatial_scale) tt = np.arange(T) gamma = temporal_amp * np.exp(-((tt - T * 0.45) / (T * 0.22)) ** 2) # rise-peak-fall gamma = gamma - gamma.mean() eta = alpha + beta * x[:, None] + phi[:, None] + gamma[None, :] y = rng.poisson(E[:, None] * np.exp(eta)) return y, E, x, dict(phi=phi, gamma=gamma, beta=beta) def st_gibbs(y, E, X, W, rng, draws=3000, burn=3000, a=1.0, b=0.01, beta_sd=5.0): """space-time ANOVA Poisson CAR: y_it ~ Poisson(E_i exp(alpha + x_i'beta + phi_i + gamma_t)). phi ~ ICAR (spatial), gamma ~ RW1 (temporal). Returns draws of alpha, beta, phi, gamma, and the spatial and temporal smoothing SDs.""" y = np.asarray(y, float); E = np.asarray(E, float); X = np.atleast_2d(np.asarray(X, float)) if X.shape[0] != y.shape[0]: X = X.T n, T = y.shape; p = X.shape[1]; nnb = W.sum(1) alpha = np.log(y.sum() / (E.sum() * T) + 1e-9); beta = np.zeros(p) phi = np.zeros(n); gamma = np.zeros(T); t2p = 0.3; t2g = 0.3 sal = 0.05; sbe = 0.05 * np.ones(p); sp = 0.3 * np.ones(n); sg = 0.2 * np.ones(T); accp = np.zeros(n) acca = 0.0; accb = np.zeros(p); accg = np.zeros(T); nadapt = 0 A = np.empty(draws); B = np.empty((draws, p)); PHI = np.empty((draws, n)); GAM = np.empty((draws, T)) SDP = np.empty(draws); SDG = np.empty(draws) for it in range(draws + burn): eta = alpha + (X @ beta)[:, None] + phi[:, None] + gamma[None, :]; lam = E[:, None] * 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)/(2*100)): alpha += d; lam = lp; acca += 1 # ---- beta ---- for j in range(p): d = sbe[j] * rng.standard_normal(); lp = lam * np.exp(d * X[:, j][:, None]) dll = d * (y * X[:, j][:, None]).sum() - (lp - lam).sum() - (( (beta[j]+d)**2 - beta[j]**2)/(2*beta_sd**2)) if np.log(rng.random()) < dll: beta[j] += d; lam = lp; accb[j] += 1 # ---- phi_i (ICAR; likelihood summed over time) ---- yr = y.sum(1) for i in range(n): nbm = (W[i] @ phi) / max(nnb[i], 1.0) # neighbours as they stand NOW d = sp[i] * rng.standard_normal() 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()) < dll + dpr: phi[i] += d; lam[i] *= np.exp(d); accp[i] += 1 phi -= phi.mean() lam = E[:, None] * np.exp(alpha + (X @ beta)[:, None] + phi[:, None] + gamma[None, :]) # ---- gamma_t (RW1; likelihood summed over regions) ---- yc = y.sum(0) for t in range(T): d = sg[t] * rng.standard_normal() dll = yc[t] * d - lam[:, t].sum() * (np.exp(d) - 1) q = 0.0 if t > 0: q += (gamma[t] + d - gamma[t-1])**2 - (gamma[t] - gamma[t-1])**2 if t < T-1: q += (gamma[t+1] - gamma[t] - d)**2 - (gamma[t+1] - gamma[t])**2 dpr = -q / (2 * t2g) if np.log(rng.random()) < dll + dpr: gamma[t] += d; lam[:, t] *= np.exp(d); accg[t] += 1 gamma -= gamma.mean() # ---- smoothing variances ---- 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))) qg = 0.5 * np.sum(np.diff(gamma) ** 2) t2g = 1.0 / rng.gamma(a + (T - 1) / 2, 1.0 / (b + max(qg, 1e-6))) if it < burn and it % 100 == 99: # adapt every scale during burn-in nadapt += 1; g = 1.0 / np.sqrt(nadapt) sp *= np.exp((accp / 100 - 0.44) * g); sg *= np.exp((accg / 100 - 0.44) * g) sal *= np.exp((acca / 100 - 0.44) * g); sbe *= np.exp((accb / 100 - 0.44) * g) accp[:] = 0; accg[:] = 0; acca = 0.0; accb[:] = 0 if it >= burn: k = it - burn; A[k] = alpha; B[k] = beta; PHI[k] = phi; GAM[k] = gamma SDP[k] = np.sqrt(t2p); SDG[k] = np.sqrt(t2g) return dict(alpha=A, beta=B, phi=PHI, gamma=GAM, spatial_sd=SDP, temporal_sd=SDG)