Riebler 2016 - Penalised Complexity Priors for BYM2

Summary

PC priors place an exponential prior on the Kullback-Leibler “distance” from a flexible model to a base model. For BYM2: $\pi ({\tau }_b)=\frac{\theta }{2}{\tau }_b^{-3/2}\exp (-\theta {\tau }_b^{-1/2})$ (type-2 Gumbel), i.e. an exponential prior on $\sigma ={\tau }_b^{-1/2}$, calibrated by $\mathrm{Pr}(\sigma >U)=\alpha$ giving $\theta =-\ln \alpha /U$; and a tabulated prior on $\varphi$ (via $\text{logit}\varphi$) calibrated by $\mathrm{Pr}(\varphi <U)=\alpha$, shrinking to $\varphi =0$.

Overview

Four principles (Simpson et al.): (1) Occam’s razor, (2) KLD measures complexity, (3) constant-rate penalisation, (4) user-defined scaling. BYM2 is a flexible extension of two nested base models: level 1 = constant risk (infinite precision); level 2 (given variance) = unstructured noise only ($\varphi =0$).

Main Content

PC prior distance

$d(\cdot )=\sqrt{2\text{KLD}(\text{flexible}\Vert \text{base})}$; prior: $d\sim \text{Exp}(\lambda )$, then transform by change of variables.

PC prior for the precision ${\tau }_b$ (Sec. 4.1, App. C.1)

Base: large fixed precision ${\tau }_0\gg {\tau }_b$. The $(n-1)$-dimensional stochastic part gives $\text{KLD}({\tau }_b)=\frac{n-1}{2}(\frac{{\tau }_0}{{\tau }_b}-1-\log \frac{{\tau }_0}{{\tau }_b})\to \frac{n-1}{2}\frac{{\tau }_0}{{\tau }_b}$, so $d({\tau }_b)\to \sqrt{(n-1){\tau }_0/{\tau }_b}$. With $\theta =\lambda \sqrt{(n-1){\tau }_0}$ held fixed as ${\tau }_0\to \infty$:

$$\pi ({\tau }_b)=\frac{\theta }{2}{\tau }_b^{-3/2}\exp (-\theta {\tau }_b^{-1/2}),\theta =-\frac{\ln \alpha }{U}\text{from }\mathrm{Pr}({\tau }_b^{-1/2}>U)=\alpha .$$

Interpretation: $\sigma ={\tau }_b^{-1/2}\sim \text{Exp}(\theta )$. Example: $\mathrm{Pr}(\sigma >1)=0.01$ is about “99% probability that residual relative risk is below 2”. The paper’s Fig. 2 uses $U=1,\alpha =0.01$ (so $\theta \approx 4.605$).

PC prior for the mixing parameter $\varphi$ (Sec. 4.2, App. C.2)

With ${\Sigma }_{base}=I$, ${\Sigma }_{flex}(\varphi )=(1-\varphi )I+\varphi Q^{\star -}$:

$$\text{KLD}(\varphi )=\frac{1}{2}(\mathrm{tr}{\Sigma }_{flex}-n-\log |{\Sigma }_{flex}|)=\frac{1}{2}(\varphi (\mathrm{tr}{Q}^{\star -}-n)-\sum _{i=1}^n\log (1-\varphi +\varphi {\tilde{\gamma }}_i)),$$

${\tilde{\gamma }}_i=1/{\gamma }_i$ for eigenvalues ${\gamma }_i>0$ of $Q^{\star }$ and $0$ for the zero eigenvalue. Put $d(\varphi )=\sqrt{2\text{KLD}(\varphi )}\sim \text{Exp}(\lambda )$ and transform to the logit scale; no closed form. Calibrate $\lambda$ by $\mathrm{Pr}(\varphi <U)=\alpha$; INLA: INLA:::inla.pc.bym.phi(Q, rankdef=1, u=0.5, alpha=2/3). Reasonable default: $\mathrm{Pr}(\varphi <0.5)=2/3$ (conservative: unstructured part dominates).

(The printed formula in the paper reads $\frac{1}{2}(n\varphi \mathrm{tr}(Q_{\star }^{-1})/n\dots )$ with an extraction ambiguity; the expression above is the algebraic expansion of the trace form and is what the code below uses.)

Why not gamma priors: gamma$(1,0.01)$ has no density at the base model, so the base model can never be reached (overfitting); uniform priors on $\varphi$ ignore the shrinkage principle. Figure 2: PC(${\tau }_b$; $U=1,\alpha =0.01$) vs Gamma(1, rate .01/.02); PC($\varphi$; $U=.5,\alpha =2/3$) vs Uniform(0,1).

Examples

Distances and PC density on a 16-node grid

Using the scaled grid from Riebler 2016 - Scaling and the BYM2 Construction: $d(\varphi )\approx$ 0.029, 0.280, 1.244, 2.298, 2.876 at $\varphi =0.01,0.1,0.5,0.9,0.99$ (verified numerically); $\pi (\tau )$ at $\tau =1,10,100$ with $\theta =4.605$ equals 0.0230, 0.0170, 0.00145 (units: per unit precision).

import numpy as np
def pc_tau_pdf(tau, U=1., alpha=.01):
    th = -np.log(alpha)/U; return th/2*tau**-1.5*np.exp(-th*tau**-.5)
def d_phi(phi, Qs):
    n = Qs.shape[0]; ev = np.linalg.eigvalsh(Qs); g = np.where(ev > 1e-9, 1/np.maximum(ev, 1e-9), 0.)
    return np.sqrt(2*0.5*(phi*(g.sum()-n) - np.sum(np.log(1-phi+phi*g))))

PyMC: sigma = pm.Exponential("sigma", rate=-np.log(alpha)/U) implements the PC prior on ${\tau }_b$ in standard-deviation form; the $\varphi$ prior needs a tabulated log-density via pm.CustomDist or a Beta approximation.

Connections

Riebler 2016 - Simulation Study and Sardinia IDDM Application tests these priors; compare to default priors in Prior Predictive Checking.

See Also

Riebler 2016 - BYM2 Intuitive Spatial Model - Overview