Riebler 2016 - Classical Disease Mapping Models Besag BYM Leroux Dean

Summary

Disease mapping uses $y_i\mid {\theta }_i\sim \text{Poisson}(E_i{\theta }_i)$ with $\log {\theta }_i=\mu +z_i^{\top }\beta +b_i$. Four latent models for $b$: Besag (intrinsic GMRF), BYM (structured + unstructured), Leroux (mixed precision) and Dean (mixed variance). None has a scaled structured part; Leroux cannot be scaled.

Overview

Raw SMRs ${\hat{\theta }}_i=y_i/E_i$ have $\mathrm{Var}({\hat{\theta }}_i)={\theta }_i/E_i$, so small-$E_i$ areas are noisy; borrow strength over neighbours. Weakly informative $N(0,100)$ priors on $\mu ,\beta$.

Main Content

Poisson disease-mapping likelihood

$y_i\mid {\theta }_i\sim \text{Poisson}(E_i{\theta }_i)$, ${\eta }_i=\log {\theta }_i=\mu +z_i^{\top }\beta +b_i$; $E_i$ expected cases from population at risk; $b$ spatial random effect for unexplained heterogeneity.

Besag (ICAR) model, eq. 1-2

$b_i\mid b_{-i},{\tau }_b\sim N(\frac{1}{n_{{\delta }_i}}{\sum }_{j\in {\delta }_i}{b}_j,\frac{1}{n_{{\delta }_i}{\tau }_b})$; joint $\pi (b\mid {\tau }_b)\propto \exp (-\frac{{\tau }_b}{2}{\sum }_{i\sim j}(b_i-b_j{)}^2)=\exp (-\frac{{\tau }_b}{2}{b}^{\top }Qb)$ with $Q_{ii}=n_{{\delta }_i}$, $Q_{ij}=-1$ if $i\sim j$, else 0. $Q$ is singular: rank deficiency equals the number of connected components $I$ (1 for a connected map), $1$ is a null vector, density $\pi (b\mid {\tau }_b)=K{\tau }_b^{(n-I)/2}\exp (-\frac{{\tau }_b}{2}{b}^{\top }Qb)$; impose sum-to-zero on each component. Cannot represent pure overdispersion.

BYM model

$b=v+u$, $v\sim N(0,{\tau }_v^{-1}I)$ (overdispersion), $u\sim N(0,{\tau }_u^{-1}{Q}^{-})$ (Besag); $\mathrm{Var}(b)={\tau }_v^{-1}I+{\tau }_u^{-1}{Q}^{-}$. Not identifiable as a split; ${\tau }_v^{-1}$ is a marginal variance whereas ${\tau }_u^{-1}$ is a conditional variance given neighbours, so priors on the two are not on the same scale.

Leroux model (eq. 3)

$\mathrm{Var}(b\mid {\tau }_b,\varphi )={\tau }_b^{-1}((1-\varphi )I+\varphi Q{)}^{-1}$, $\varphi \in [0,1]$: pure overdispersion at $\varphi =0$, Besag at $\varphi =1$; the precision matrix is a weighted average of $I$ and $Q$.

Dean model (eq. 4-5)

$b={\tau }_b^{-1/2}(\sqrt{1-\varphi }v+\sqrt{\varphi }u)$, $\mathrm{Var}(b)={\tau }_b^{-1}((1-\varphi )I+\varphi Q^{-})$: a BYM reparameterisation with ${\tau }_u^{-1}={\tau }_b^{-1}\varphi$, ${\tau }_v^{-1}={\tau }_b^{-1}(1-\varphi )$; the covariance is a weighted average.

Examples

Neighbour-list format for INLA (Appendix B): first line $n$; each following line = region index, number of neighbours, neighbour indices. In PyMC, pm.ICAR(W=adjacency, sigma=...) implements the Besag density (E18).

Connections

Riebler 2016 - Scaling and the BYM2 Construction fixes the unscaled-precision problem; Spatial Models - BYM gives the basic BYM tutorial.

See Also

Riebler 2016 - BYM2 Intuitive Spatial Model - Overview