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