Spatiotemporal Modelling — the Space-Time CAR

Space + time, the capstone of the spatial arc

The final piece: data indexed by both region and time. Disease counts, mortality, disability receipt, unemployment — all vary across areas and years, and modelling them jointly unifies the areal CAR work (space) with the time-series work (time). The workhorse is the space-time ANOVA (Knorr-Held) decomposition of the log relative risk: $$y_{it}\sim\text{Poisson}\!\big(E_{it}\,e^{\alpha+x_i'\beta+\phi_i+\gamma_t}\big),$$ with a spatial main effect $\phi_i$ (a CAR/ICAR field — which regions run high), a temporal main effect $\gamma_t$ (a random walk — the shared trend across regions), and optionally a space-time interaction $\delta_{it}$. It is the areal-CAR machinery in space and a random-walk smoother in time.

Data. We use real Georgia county geography (159 counties, real adjacency and population) and simulate a 12-year disease process on it — a smooth spatial field, an epidemic-style temporal rise-and-fall, and a covariate — so we can validate that the model recovers the truth. The marquee real application is the disability belt over time (whether the belt is expanding or stable), which awaits the multi-year SSA county panel; the method here is exactly what that analysis needs. We fit from scratch, map the space-time structure, and cross-check in PyMC.

import numpy as np, pandas as pd, matplotlib.pyplot as plt, geopandas as gpd, contextily as cx
import spacetime as ST
rng = np.random.default_rng(0)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; GREY="#cccccc"
d = pd.read_csv("ga_pop.csv", dtype={"fips":str}); W = pd.read_csv("ga_adj.csv").to_numpy(float); n=len(d)
gdf = gpd.read_file("ga_counties.geojson").merge(d.assign(order=range(n)), left_on="id", right_on="fips").sort_values("order").reset_index(drop=True)
pop = d["pop"].to_numpy()
y, E, x, truth = ST.simulate_st(pop, W, rng, T=12)               # 12-year simulated process on real GA geography
pd.DataFrame(y, columns=[f"y{t}" for t in range(12)]).assign(fips=d["fips"], x=x, pop=pop).to_csv("ga_st_counts.csv", index=False)
def choro(ax, vals, title, cmap="OrRd", vmin=None, vmax=None):
    g=gdf.copy(); g["v"]=vals; g.plot(column="v",ax=ax,cmap=cmap,edgecolor="white",linewidth=0.2,legend=True,legend_kwds={"shrink":0.6},vmin=vmin,vmax=vmax)
    ax.set_title(title,fontsize=10); ax.axis("off")
print(f"Georgia: {n} counties x 12 years; {int(y.sum()):,} simulated events; real population {pop.sum():,}")
print("Real geography + population, simulated space-time disease process (known truth for validation).")

Georgia: 159 counties x 12 years; 59,114 simulated events; real population 11,029,227
Real geography + population, simulated space-time disease process (known truth for validation).

1. The setting

Georgia's 159 counties (located on the map for context) carry a simulated 12-year disease process. The right panel shows the observed rate over the whole period — a spatial pattern is visible, but it is blurred by year-to-year Poisson noise and the shared temporal wave. The model will separate space, time, and covariate.

loc = gdf.to_crs(3857)
fig,ax=plt.subplots(1,2,figsize=(9.8,3.8))
loc.plot(ax=ax[0], color="none", edgecolor=BLUE, linewidth=0.5)
try: cx.add_basemap(ax[0], source=cx.providers.OpenStreetMap.Mapnik, crs=3857, attribution_size=5)
except Exception as e: print("basemap unavailable:", str(e)[:60])
ax[0].set_title("Georgia counties (location)"); ax[0].axis("off")
overall = y.sum(1)/(E*12)                                        # observed events per 1000 per year, all years
choro(ax[1], overall, "Observed rate per 1,000/yr (all 12 years)", cmap="OrRd")
plt.tight_layout(); plt.show()
print("The overall map hints at a spatial pattern, but it mixes the true spatial signal with Poisson noise and the")
print("temporal wave. Disentangling the three is the job of the space-time model.")

The overall map hints at a spatial pattern, but it mixes the true spatial signal with Poisson noise and the
temporal wave. Disentangling the three is the job of the space-time model.

2. Fit and recovery

We fit the space-time model (spatial ICAR $\phi$, temporal random-walk $\gamma$, covariate $\beta$) from scratch. Because the data are simulated, we can check it against the truth: the estimated spatial field, temporal wave, and covariate should match what generated the data.

res = ST.st_gibbs(y, E, x[:,None], W, rng, draws=2500, burn=2500)
phi_e=res["phi"].mean(0); gam_e=res["gamma"].mean(0); be=res["beta"].mean(0)[0]
print(f"covariate beta: true {truth['beta']:.2f}, estimated {be:.2f} [{np.percentile(res['beta'],2.5):.2f}, {np.percentile(res['beta'],97.5):.2f}]")
print(f"spatial field recovery: corr {np.corrcoef(truth['phi'],phi_e)[0,1]:.2f}   temporal wave recovery: corr {np.corrcoef(truth['gamma'],gam_e)[0,1]:.2f}")
fig,ax=plt.subplots(1,3,figsize=(12,3))
vmax=max(abs(truth['phi']).max(),abs(phi_e).max())
choro(ax[0], truth['phi'], "TRUE spatial field", cmap="RdBu_r", vmin=-vmax, vmax=vmax)
choro(ax[1], phi_e, "ESTIMATED spatial field", cmap="RdBu_r", vmin=-vmax, vmax=vmax)
ax[2].plot(range(12), truth['gamma'], "o-", color=GREEN, label="true temporal wave")
ax[2].plot(range(12), gam_e, "s--", color=RED, label="estimated"); ax[2].axhline(0,color="k",lw=.5)
ax[2].set_xlabel("year"); ax[2].set_ylabel(r"temporal effect $\gamma_t$"); ax[2].set_title("Temporal trend recovered"); ax[2].legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
print(f"The model recovers all three components -- spatial map (corr {np.corrcoef(truth['phi'],phi_e)[0,1]:.2f}), the epidemic-style")
print(f"temporal wave (corr ~1.0), and the covariate ({be:.2f} vs {truth['beta']:.2f}) -- separating space, time, and covariate")
print("from a single noisy space-time count panel.")

covariate beta: true 0.30, estimated 0.29 [0.25, 0.32]
spatial field recovery: corr 0.99   temporal wave recovery: corr 1.00

The model recovers all three components -- spatial map (corr 0.99), the epidemic-style
temporal wave (corr ~1.0), and the covariate (0.29 vs 0.30) -- separating space, time, and covariate
from a single noisy space-time count panel.

3. The decomposition, mapped through time

The fitted relative risk $\text{RR}_{it}=e^{\alpha+x_i\beta+\phi_i+\gamma_t}$ is a surface over space and time. With a shared temporal main effect, the spatial pattern holds its shape while the whole map rises and falls with $\gamma_t$ — the epidemic sweeping the state. A few years show it (a space-time interaction term, an extension, would let individual counties depart from the shared wave).

al=res["alpha"].mean(); years=[0,4,6,11]
RR={t:np.exp(al+be*x+phi_e+gam_e[t]) for t in years}
vmax=max(v.max() for v in RR.values())
fig,ax=plt.subplots(1,4,figsize=(12,2.7))
for k,t in enumerate(years):
    choro(ax[k], RR[t], f"year {t}  (RR)", cmap="OrRd", vmin=0, vmax=vmax)
plt.tight_layout(); plt.show()
print("The high-risk geography is stable while intensity rises to the mid-series peak (year ~6) and recedes -- the")
print("space-time surface. Spatial main effect = WHERE, temporal main effect = WHEN, and together the risk everywhere.")

The high-risk geography is stable while intensity rises to the mid-series peak (year ~6) and recedes -- the
space-time surface. Spatial main effect = WHERE, temporal main effect = WHEN, and together the risk everywhere.

4. Cross-check in PyMC

The same model in PyMC: a Poisson likelihood on the count panel with an intercept, covariate, an intrinsic-CAR spatial field, and a Gaussian random-walk temporal field. We confirm the covariate and the spatial/temporal magnitudes.

import pymc as pm
adj=W.astype(int)
with pm.Model() as mod:
    a0=pm.Normal("a0",0,5); b=pm.Normal("b",0,5); tp=pm.HalfNormal("tp",2); tg=pm.HalfNormal("tg",2)
    phi=pm.ICAR("phi",W=adj)
    import pytensor.tensor as pt
    # explicit non-centred random walk, summing to zero -- the same parameterisation the
    # from-scratch sampler uses. GaussianRandomWalk without an init_dist defaults to
    # Normal(0, 100), whose free level fights the intercept and wrecks the geometry.
    dg=pm.Normal("dg",0,1,shape=y.shape[1]-1)
    g=pt.concatenate([pt.zeros(1), pt.cumsum(dg)]); g=g-g.mean()
    eta=np.log(E)[:,None]+a0+b*x[:,None]+tp*phi[:,None]+tg*g[None,:]
    pm.Poisson("y", mu=pm.math.exp(eta), observed=y)
    idata=pm.sample(2000,tune=2000,chains=4,target_accept=0.95,random_seed=3,progressbar=False)
b_pm=float(idata.posterior["b"].mean())
phi_pm=(idata.posterior["tp"].values.reshape(-1,1)*idata.posterior["phi"].stack(s=("chain","draw")).values.T).mean(0)
import arviz as az
V=["a0","b","tp","tg","phi","dg"]
print(f"PyMC convergence: max r-hat {max(float(v.max()) for v in az.rhat(idata,var_names=V).data_vars.values()):.3f}, "
      f"min ESS {min(float(v.min()) for v in az.ess(idata,var_names=V).data_vars.values()):.0f}")
print(f"covariate beta:  from-scratch {be:.2f}   PyMC {b_pm:.2f}   (true {truth['beta']:.2f})")
print(f"spatial field agreement (from-scratch vs PyMC): corr {np.corrcoef(phi_e,phi_pm)[0,1]:.2f}")
print("PyMC's ICAR-spatial + random-walk-temporal Poisson reproduces the from-scratch space-time decomposition.")

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, tp, tg, phi, dg]

Sampling 4 chains for 2_000 tune and 2_000 draw iterations (8_000 + 8_000 draws total) took 19 seconds.

PyMC convergence: max r-hat 1.010, min ESS 406
covariate beta:  from-scratch 0.29   PyMC 0.29   (true 0.30)
spatial field agreement (from-scratch vs PyMC): corr 1.00
PyMC's ICAR-spatial + random-walk-temporal Poisson reproduces the from-scratch space-time decomposition.

Geometry, not compute — the two parameterisations timed side by side

The cell above uses an explicit non-centred random walk constrained to sum to zero. The obvious alternative is PyMC's GaussianRandomWalk, which without an init_dist defaults to $\mathcal N(0,100)$ — and that free level is unidentified against the intercept, so the sampler has to explore a long flat ridge instead of a well-conditioned posterior.

The claim worth checking is that this is a geometry problem rather than a compute one. Both versions are fitted below at identical settings — same draws, same chains, same target_accept — and timed.

import time, arviz as az
import pytensor.tensor as pt

def fit(kind, draws=1000, tune=1000, chains=2, seed=3):
    with pm.Model():
        a0 = pm.Normal("a0", 0, 5); b = pm.Normal("b", 0, 5)
        tp = pm.HalfNormal("tp", 2); tg = pm.HalfNormal("tg", 2)
        phi = pm.ICAR("phi", W=adj)
        if kind == "explicit":                       # non-centred increments, summing to zero
            dg = pm.Normal("dg", 0, 1, shape=y.shape[1]-1)
            g = pt.concatenate([pt.zeros(1), pt.cumsum(dg)]); g = g - g.mean()
        else:                                        # library default: free level, N(0,100)
            g = pm.GaussianRandomWalk("g", sigma=1.0, shape=y.shape[1])
        eta = np.log(E)[:,None] + a0 + b*x[:,None] + tp*phi[:,None] + tg*g[None,:]
        pm.Poisson("y", mu=pm.math.exp(eta), observed=y)
        t0 = time.perf_counter()
        idt = pm.sample(draws, tune=tune, chains=chains, target_accept=0.95,
                        random_seed=seed, progressbar=False)
        el = time.perf_counter() - t0
    V = ["a0","b","tp","tg"]
    rh = max(float(v.max()) for v in az.rhat(idt, var_names=V).data_vars.values())
    es = min(float(v.min()) for v in az.ess(idt, var_names=V).data_vars.values())
    dv = int(idt.sample_stats["diverging"].values.sum())
    td = int((idt.sample_stats["tree_depth"].values >= 10).sum()) if "tree_depth" in idt.sample_stats else -1
    return el, rh, es, dv, td, chains*draws

rows = [("explicit sum-to-zero", *fit("explicit")), ("GaussianRandomWalk default", *fit("default"))]
print(f"{'parameterisation':28s}{'draws':>7s}{'seconds':>9s}{'max r-hat':>11s}{'min ESS':>9s}{'diverg':>8s}{'maxtree':>9s}")
for nm, el, rh, es, dv, td, nd in rows:
    print(f"{nm:28s}{nd:>7d}{el:>9.0f}{rh:>11.3f}{es:>9.0f}{dv:>8d}{td:>9d}")

(e1,r1,s1,d1,t1,n1),(e2,r2,s2,d2,t2,n2) = [r[1:] for r in rows]
print(f"\nSame model, same data, same {n1} draws. The default parameterisation takes {e2/e1:.1f} times as long")
print(f"and returns a far worse chain: minimum ESS {s2:.0f} against {s1:.0f} and r-hat {r2:.3f} against {r1:.3f}.")
print(f"The diagnostic that names the cause is tree depth -- {t2} of {n2} draws hit the maximum against {t1}")
print(f"for the explicit version. (Neither run diverges: {d2} and {d1}. Divergence is not the failure mode here.)")
print("Hitting the cap every draw means the sampler is taking the longest trajectory it is allowed and still")
print("not turning around -- the signature of an unidentified direction, here the random walk's free level")
print("fighting the intercept. The extra time is the symptom rather than the cause, which is why the fix is a")
print("reparameterisation and not a faster machine. For reference the published 4-chain fit above, using the")
print("explicit form, reaches r-hat 1.010 and minimum ESS 406.")

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (2 chains in 2 jobs)

NUTS: [a0, b, tp, tg, phi, dg]

Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 9 seconds.

We recommend running at least 4 chains for robust computation of convergence diagnostics

The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details

The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\pymc\distributions\timeseries.py:290: UserWarning: Initial distribution not specified, defaulting to `Normal.dist(0, 100)`.You can specify an init_dist manually to suppress this warning.
  warnings.warn(

Initializing NUTS using jitter+adapt_diag...

Rewrite failure due to: local_subtensor_merge_slice

node: Subtensor{i}(Subtensor{start:}.0, 0)

TRACEBACK:

Traceback (most recent call last):
  File "C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\pytensor\graph\rewriting\basic.py", line 1920, in process_node
    replacements = node_rewriter.transform(
                   ^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\pytensor\graph\rewriting\basic.py", line 993, in transform
    return self.fn(fgraph, node)
           ^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\pytensor\tensor\rewriting\subtensor.py", line 973, in local_subtensor_merge_slice
    return _local_subtensor_merge_rewrite(fgraph, node, merge_integer_index=False)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\pytensor\tensor\rewriting\subtensor.py", line 930, in _local_subtensor_merge_rewrite
    merged = merge_two_slices(
             ^^^^^^^^^^^^^^^^^
  File "C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\pytensor\tensor\rewriting\subtensor.py", line 1612, in merge_two_slices
    p_val = sl1.start + sl2 * sl1.step
            ~~~~~~~~~~^~~~~~~~~~~~~~~~
OverflowError: Python integer 163 out of bounds for int8

Multiprocess sampling (2 chains in 2 jobs)

NUTS: [a0, b, tp, tg, phi, g]

Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 63 seconds.

Chain 0 reached the maximum tree depth. Increase `max_treedepth`, increase `target_accept` or reparameterize.

Chain 1 reached the maximum tree depth. Increase `max_treedepth`, increase `target_accept` or reparameterize.

We recommend running at least 4 chains for robust computation of convergence diagnostics

The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details

The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details

parameterisation              draws  seconds  max r-hat  min ESS  diverg  maxtree
explicit sum-to-zero           2000       10      1.038       92       0        0
GaussianRandomWalk default     2000       64      1.100       15       0     2000

Same model, same data, same 2000 draws. The default parameterisation takes 6.4 times as long
and returns a far worse chain: minimum ESS 15 against 92 and r-hat 1.100 against 1.038.
The diagnostic that names the cause is tree depth -- 2000 of 2000 draws hit the maximum against 0
for the explicit version. (Neither run diverges: 0 and 0. Divergence is not the failure mode here.)
Hitting the cap every draw means the sampler is taking the longest trajectory it is allowed and still
not turning around -- the signature of an unidentified direction, here the random walk's free level
fighting the intercept. The extra time is the symptom rather than the cause, which is why the fix is a
reparameterisation and not a faster machine. For reference the published 4-chain fit above, using the
explicit form, reaches r-hat 1.010 and minimum ESS 406.

5. Summary

The space-time CAR closes the spatial arc by adding time to the areal model. Its ANOVA decomposition splits a region-by-year count panel into a spatial main effect (a CAR field — where risk is high), a temporal main effect (a random walk — when it is high, here an epidemic-style rise and fall), a covariate effect, and optionally a space-time interaction. On real Georgia geography with a simulated process, the from-scratch sampler recovered all three cleanly (spatial-field correlation ≈ 0.99, temporal wave ≈ 1.0, covariate on target), and PyMC reproduced it with an intrinsic-CAR spatial field and a Gaussian-random-walk temporal field.

This is the synthesis of the whole arc: the CAR of the areal notebooks in space, and a random walk — the temporal cousin of the CAR, and kin to the state-space/BVAR time-series work — in time. Two caveats worth carrying: the spatial main effect and a spatially-structured covariate can be confounded (the field absorbs the covariate's effect — the reason the simulated covariate here is non-spatial), a live issue in spatial regression; and the additive model assumes every region follows the same temporal wave, which a space-time interaction relaxes. The obvious real application is the disability belt over time — is it expanding, stable, or shifting? — using SSA's yearly county panel, the natural continuation of The Disability Belt — a Spatial Application. With that, the spatial arc spans areal, econometric, geostatistical, point-process, and spatiotemporal models, each from scratch, cross-checked, and mapped.