Spatiotemporal Modelling — the Space-Time CAR

Space and Time Together

The method capstone of the arc, and the point where the spatial work meets the time-series work. Data indexed by both region and year — disease counts, mortality, benefit receipt — support a decomposition of the log relative risk into a spatial main effect (which regions run high, averaged over time), a temporal main effect (the trend shared by every region), and, optionally, a space-time interaction for regions whose trajectories depart from that shared wave. This is the Knorr-Held ANOVA, and $\phi$ is exactly the areal CAR machinery from the start of the arc with a random walk in time bolted alongside it.

The teaching decision here is a good one: real geography, simulated process. Georgia's 159 counties and their true populations carry a known spatial field, a known epidemic-style temporal wave that rises to a mid-series peak and recedes, and a known covariate effect — 59,114 events across 12 years. Because the truth is known, recovery is checkable rather than plausible, which is the only way to demonstrate that three components really can be separated from one noisy panel.

Does it recover the truth?

They can. All three engines recover all three components, and agree with each other: the covariate at 0.29 against a true 0.30, the spatial map at a correlation of 0.99, and the temporal wave at 1.00. That is the claim the example exists to make, and it is made against known truth rather than asserted from a plausible-looking map.

What is fitted, and what is left out

Worth being precise about what is fitted: the from-scratch sampler and the R model are both the additive form — two main effects, no interaction. R's ST.CARanova is called with interaction=FALSE deliberately, to match. The interaction is the least identified part of the decomposition, and on data generated without one it is noise to fit; interaction=TRUE and the autoregressive ST.CARar are the routes to it when the trajectories genuinely diverge.

A reparameterisation that mattered

The PyMC cross-check turns on how the temporal walk is written, and the difference is worth measuring rather than asserting. A GaussianRandomWalk without an explicit initial distribution defaults to $N(0,100)$: that free level is unidentified against the intercept, so the posterior has a long flat ridge along which the two trade off. Writing the walk instead as explicit non-centred increments constrained to sum to zero — the same form the from-scratch sampler uses — removes the ridge, and the notebook fits both versions at identical settings to show what it costs. On the same 2,000 draws the default parameterisation takes 6.4× as long (64 seconds against 10) and returns a far worse chain: minimum ESS 15 against 92, $\hat{R}$ 1.100 against 1.038. The diagnostic that names the cause is tree depth — every one of the 2,000 draws hits the maximum against none for the explicit form, meaning the sampler takes the longest trajectory it is permitted and still cannot turn around. Neither run diverges, so divergence is not the failure mode here. The published fit above, using the explicit form at four chains, reaches $\hat{R}=1.010$ and a minimum ESS of 406. Geometry, not compute — the extra time is the symptom, and the fix is a reparameterisation rather than a faster machine.

Notebooks

Downloads

Space-Time Module — Source Code

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

References

• Knorr-Held, L. (2000). Bayesian modelling of inseparable space-time variation in disease risk. Statistics in Medicine 19(17–18), 2555–2567. — the ANOVA decomposition and its interaction types
• Bernardinelli, L. et al. (1995). Bayesian analysis of space-time variation in disease risk. Statistics in Medicine 14(21–22), 2433–2443. — the space-time disease-mapping problem
• Lee, D., Rushworth, A. & Napier, G. (2018). Spatio-temporal areal unit modelling in R with conditional autoregressive priors using the CARBayesST package. Journal of Statistical Software 84(9). — ST.CARanova and ST.CARar
• Rue, H. & Held, L. (2005). Gaussian Markov Random Fields. Chapman & Hall/CRC. — the random-walk and CAR priors as GMRFs
• Betancourt, M. & Girolami, M. (2015). Hamiltonian Monte Carlo for hierarchical models. In Current Trends in Bayesian Methodology. CRC Press. — why the non-centred parameterisation fixes the geometry