Joint Disease Mapping — the Shared-Component CAR

Python · PyMC · R (CARBayes)  ·  Download shared-component module

One Geography, Two Outcomes

The foundations example mapped one outcome. Related outcomes often share a geography — two diseases, incidence and mortality, or the same disease in two periods — and modelling them together both borrows strength and answers a question a single map cannot: is the geography of risk the same? Here the outcomes are North Carolina SIDS deaths in 1974–78 (667) and 1979–84 (836) across the same 100 counties.

The shared-component model splits each outcome's spatially structured log-risk into a common field ss and an outcome-specific field ϕk\phi_k, each an intrinsic CAR. Because the three fields are independent and share the same neighbourhood precision, the correlation between the two periods' spatial effects has a clean closed form, and the shared field is the stable latent risk surface both periods draw on. The sampler is Metropolis-within-Gibbs, with the shared field updated against both Poisson likelihoods at once.

yki∼Poisson(Eki eαk+si+ϕki),corr=τs2(τs2+τ12)(τs2+τ22)y_{ki} \sim \text{Poisson}\big(E_{ki}\,e^{\alpha_k + s_i + \phi_{ki}}\big), \qquad \text{corr} = \frac{\tau_s^2}{\sqrt{(\tau_s^2+\tau_1^2)(\tau_s^2+\tau_2^2)}}

The raw maps correlate only 0.21, but that number is mostly small-area noise. Fitting the model, most of each period's spatial signal turns out to be shared — 73% and 87% — for a cross-period correlation of 0.78, and the smoothed maps correlate 0.84. The high-risk core is largely stable; what changes between periods is a smaller residue. A single-period map could not have separated the two.

NC SIDS, 1974–78 vs 1979–84correlation between periods
raw standardised ratios0.21 — mostly small-area noise
smoothed relative risks0.84
spatial effects, shared-component (from scratch)0.78 [0.48, 0.95]
spatial effects, shared-component (PyMC, half-normal priors)0.65
between-outcome Σ\Sigma, Leroux CAR (R, different model)0.53 [0.11, 0.86]

How much of this is actually identified?

That said, the shared fraction should be read as an indication, not an estimate. Simulating eight datasets from a known shared fraction of 0.50 and refitting, the estimate comes back anywhere in 0.14 to 0.86 (sd 0.25) — the split between a "shared" and a "specific" field is close to unidentified, because the data see only their sum. The cross-period correlation is roughly twice as stable (0.25–0.63, sd 0.13), since it depends on the same variances only through a ratio. Sweeping the inverse-gamma prior on the spatial variances moves the shared fraction from 0.55 to 0.70 on the real data, which is the honest measure of how much the prior is doing.

Why three engines give three numbers

This is why the three engines land at 0.78 (from scratch), 0.65 (PyMC) and 0.53 (R) rather than on one number — and why that is not a discrepancy to explain away. The from-scratch sampler and PyMC fit the same model under different priors (inverse-gamma on the variances against half-normal on the standard deviations), and PyMC's value sits beside what the from-scratch sampler gives under a comparable prior (0.62). R fits a genuinely different parametrisation — MVS.CARleroux, a joint Leroux CAR reading the correlation off a between-outcome covariance Σ\Sigma — and reports 0.53 with a 95% interval of [0.11, 0.86]. All three intervals overlap heavily, and the shared verdict is the robust one: related but not identical.

Notebooks

Python The shared-component sampler from scratch, the variance decomposition into shared and period-specific fields, an identification check that sweeps the variance prior and simulates from a known shared fraction, and a PyMC cross-check (Python) Open → R Independent R implementation of the alternative parametrisation — CARBayes MVS.CARleroux fits a joint Leroux CAR and reads the cross-period correlation off the between-outcome covariance (R) Open →

Downloads

sp_mcar.py Shared-component CAR for two Poisson outcomes — a common field updated against both likelihoods, outcome-specific fields, adaptive proposals, and the derived shared fractions and cross-outcome correlation (NumPy) nc_sids2.csv North Carolina SIDS in two periods — births and deaths for 1974–78 and 1979–84 across 100 counties nc_adj.csv The 100 × 100 binary adjacency matrix nc_counties.geojson County boundaries for the choropleth maps

Shared-Component Module — Source Code

"""
mcar.py -- MULTIVARIATE / JOINT areal disease mapping: the shared-component CAR (from scratch).

Backs the notebooks in  "Joint Disease Mapping -- the Shared-Component CAR".

The foundations project mapped ONE outcome. Often several related outcomes live on the same
regions -- two diseases, two time periods, incidence and mortality -- and modelling them JOINTLY
both borrows strength and answers a new question: do they share a spatial pattern? Here the two
outcomes are SIDS counts in North Carolina over 1974-78 and 1979-84 across the 100 counties, and the
question is whether the geography of risk is STABLE across the two periods.

The SHARED-COMPONENT model (Knorr-Held & Best) splits each outcome's spatially structured log-risk
into a COMMON field and an outcome-SPECIFIC field:

    y1_i ~ Poisson(E1_i * exp(alpha1 + s_i + phi1_i))
    y2_i ~ Poisson(E2_i * exp(alpha2 + s_i + phi2_i))

with s, phi1, phi2 each an intrinsic-CAR field (borrowing from neighbours). The SHARED field s is
common to both outcomes -- the stable, latent spatial risk surface -- while phi1, phi2 capture what
is idiosyncratic to each. Because the three fields are independent, the total spatial effect of
period k, s + phi_k, has variance tau_s^2 + tau_k^2, and the two periods' spatial effects covary
only through s, so their CORRELATION is

    corr = tau_s^2 / sqrt( (tau_s^2 + tau_1^2)(tau_s^2 + tau_2^2) ),

and the SHARED FRACTION for period k is tau_s^2 / (tau_s^2 + tau_k^2). A high shared fraction means
the spatial pattern is common to both periods (stable geography); a low one means each period has its
own map. The sampler is Metropolis-within-Gibbs, as in the foundations project, except that the
shared field s is updated against BOTH Poisson likelihoods at once. This is the shared-component
counterpart to the multivariate-CAR (MVS.CARleroux) model fitted by CARBayes in the R notebook.
"""

import numpy as np


def sharedcar_gibbs(y1, E1, y2, E2, W, rng, draws=4000, burn=4000, a=1.0, b=0.01):
    """Shared-component CAR for two Poisson outcomes on the same regions. Returns posterior draws of
    the intercepts, the three spatial SDs (shared, and period-specific 1 & 2), the derived shared
    fractions and cross-period correlation, and the two relative-risk surfaces."""
    y1 = np.asarray(y1, float); E1 = np.asarray(E1, float); y2 = np.asarray(y2, float); E2 = np.asarray(E2, float)
    n = len(y1); nnb = W.sum(1)
    a1 = np.log((y1.sum() + 1) / (E1.sum() + 1)); a2 = np.log((y2.sum() + 1) / (E2.sum() + 1))
    s = np.zeros(n); p1 = np.zeros(n); p2 = np.zeros(n)
    ts2 = 0.3; t12 = 0.3; t22 = 0.3
    ss = 0.4 * np.ones(n); s1 = 0.4 * np.ones(n); s2 = 0.4 * np.ones(n)
    accs = np.zeros(n); acc1 = np.zeros(n); acc2 = np.zeros(n)
    acca = np.zeros(2); sa = np.array([0.05, 0.05]); nadapt = 0
    A1 = np.empty(draws); A2 = np.empty(draws)
    TS = np.empty(draws); T1 = np.empty(draws); T2 = np.empty(draws)
    SH1 = np.empty(draws); SH2 = np.empty(draws); COR = np.empty(draws)
    RR1 = np.empty((draws, n)); RR2 = np.empty((draws, n)); S = np.empty((draws, n))
    for it in range(draws + burn):
        eta1 = a1 + s + p1; lam1 = E1 * np.exp(eta1)
        eta2 = a2 + s + p2; lam2 = E2 * np.exp(eta2)
        # ---- intercepts (RW-Metropolis) ----
        ap = a1 + sa[0] * rng.standard_normal(); lp = E1 * np.exp(ap + s + p1)
        if np.log(rng.random()) < (y1 @ (np.log(lp) - np.log(lam1)) - (lp - lam1).sum() - (ap**2 - a1**2)/(2*100)):
            a1 = ap; lam1 = lp; acca[0] += 1
        ap = a2 + sa[1] * rng.standard_normal(); lp = E2 * np.exp(ap + s + p2)
        if np.log(rng.random()) < (y2 @ (np.log(lp) - np.log(lam2)) - (lp - lam2).sum() - (ap**2 - a2**2)/(2*100)):
            a2 = ap; lam2 = lp; acca[1] += 1
        # ---- shared field s (updated against BOTH likelihoods) ----
        for i in range(n):
            nbm_i = (W[i] @ s) / max(nnb[i], 1.0)
            d = ss[i] * rng.standard_normal()
            dll = y1[i]*d - lam1[i]*(np.exp(d)-1) + y2[i]*d - lam2[i]*(np.exp(d)-1)
            dpr = -(nnb[i]/(2*ts2)) * ((s[i]+d - nbm_i)**2 - (s[i] - nbm_i)**2)
            if np.log(rng.random()) < dll + dpr:
                s[i] += d; lam1[i] *= np.exp(d); lam2[i] *= np.exp(d); accs[i] += 1
        s -= s.mean()
        lam1 = E1 * np.exp(a1 + s + p1); lam2 = E2 * np.exp(a2 + s + p2)
        # ---- period-specific fields ----
        for i in range(n):
            nbm_i = (W[i] @ p1) / max(nnb[i], 1.0)
            d = s1[i] * rng.standard_normal()
            dll = y1[i]*d - lam1[i]*(np.exp(d)-1)
            dpr = -(nnb[i]/(2*t12)) * ((p1[i]+d - nbm_i)**2 - (p1[i] - nbm_i)**2)
            if np.log(rng.random()) < dll + dpr:
                p1[i] += d; lam1[i] *= np.exp(d); acc1[i] += 1
        p1 -= p1.mean()
        for i in range(n):
            nbm_i = (W[i] @ p2) / max(nnb[i], 1.0)
            d = s2[i] * rng.standard_normal()
            dll = y2[i]*d - lam2[i]*(np.exp(d)-1)
            dpr = -(nnb[i]/(2*t22)) * ((p2[i]+d - nbm_i)**2 - (p2[i] - nbm_i)**2)
            if np.log(rng.random()) < dll + dpr:
                p2[i] += d; lam2[i] *= np.exp(d); acc2[i] += 1
        p2 -= p2.mean()
        # ---- variance components (Gibbs, ICAR quadratic form) ----
        qs = 0.5 * (nnb * s * s).sum() - 0.5 * (s * (W @ s)).sum()
        q1 = 0.5 * (nnb * p1 * p1).sum() - 0.5 * (p1 * (W @ p1)).sum()
        q2 = 0.5 * (nnb * p2 * p2).sum() - 0.5 * (p2 * (W @ p2)).sum()
        ts2 = 1.0 / rng.gamma(a + (n-1)/2, 1.0/(b + max(qs, 1e-6)))
        t12 = 1.0 / rng.gamma(a + (n-1)/2, 1.0/(b + max(q1, 1e-6)))
        t22 = 1.0 / rng.gamma(a + (n-1)/2, 1.0/(b + max(q2, 1e-6)))
        if it < burn and it % 100 == 99:                               # adapt every scale during burn-in
            nadapt += 1; g = 1.0 / np.sqrt(nadapt)
            ss *= np.exp((accs/100 - 0.44) * g); s1 *= np.exp((acc1/100 - 0.44) * g)
            s2 *= np.exp((acc2/100 - 0.44) * g); sa *= np.exp((acca/100 - 0.44) * g)
            accs[:] = 0; acc1[:] = 0; acc2[:] = 0; acca[:] = 0
        if it >= burn:
            k = it - burn; A1[k] = a1; A2[k] = a2
            TS[k] = np.sqrt(ts2); T1[k] = np.sqrt(t12); T2[k] = np.sqrt(t22)
            SH1[k] = ts2/(ts2+t12); SH2[k] = ts2/(ts2+t22)
            COR[k] = ts2/np.sqrt((ts2+t12)*(ts2+t22))
            RR1[k] = np.exp(a1 + s + p1); RR2[k] = np.exp(a2 + s + p2); S[k] = s
    return dict(alpha1=A1, alpha2=A2, sd_shared=TS, sd_spec1=T1, sd_spec2=T2,
                shared_frac1=SH1, shared_frac2=SH2, corr=COR, RR1=RR1, RR2=RR2, shared=S)


def expected_counts(y, pop):
    y = np.asarray(y, float); pop = np.asarray(pop, float); return pop * (y.sum() / pop.sum())

References