Riebler 2016 - Scaling and the BYM2 Construction

Summary

The marginal variances of an intrinsic GMRF depend on the graph, so a fixed prior on ${\tau }_b$ gives different smoothing on different maps. Scale $Q$ so the geometric mean of the marginal variances (the generalised variance ${\sigma }_{GV}^2$) is 1; then BYM2 $b={\tau }_b^{-1/2}(\sqrt{1-\varphi }v+\sqrt{\varphi }{u}^{\star })$ has marginal variance about $1/{\tau }_b$ and $\varphi$ is the spatial fraction.

Overview

Graph-dependence example (Fig. 1): with ${\tau }_b=1$, ${\sigma }_{GV}^2\approx 0.29$ for Mittelfranken (12 districts), $0.40$ for Arnsberg (12 districts), $0.56$ for all 544 German districts, $4.78$ for the 3,111 continental US counties. Equal $n$ does not mean equal scale.

Main Content

Generalised variance and scaling (eq. 6)

$${\sigma }_{GV}^2(u)=\exp \left(\frac{1}{n}\sum _{i=1}^n\log \right({\tau }_b^{-1}[Q^{-}{]}_{ii})).$$

Scaled structure matrix $Q^{\star }={\sigma }_{GV}^2({\tau }_b=1)\cdot Q$ so that the geometric mean of $[Q^{\star -}{]}_{ii}$ equals 1.

Computation: avoid dense generalised inverse; for null space $V$ (e.g. $V=1$) add small $ϵ$ to the diagonal and use the conditional-on-constraints marginal variances $\mathrm{Diag}(\mathrm{Var}(b\mid V^{\top }b=0))$ with sparse algebra (Rue & Held 2005, Sec. 3.2). INLA: inla.scale.model(Q, constr=list(A=matrix(1,1,n), e=0)). Scaling also matters for hypothesis testing on IGMRFs. The Leroux model cannot be scaled because the scale depends on $\varphi$.

BYM2 model (eq. 7-8)

$$b=\frac{1}{\sqrt{{\tau }_b}}(\sqrt{1-\varphi }v+\sqrt{\varphi }{u}^{\star }),v\sim N(0,I),u^{\star }\sim N(0,Q^{\star -}),\mathrm{Var}(b)={\tau }_b^{-1}((1-\varphi )I+\varphi Q^{\star -}).$$

$\varphi \in [0,1]$ is the proportion of marginal variance explained by the structured effect; $\varphi =0$ pure overdispersion, $\varphi =1$ Besag. Marginal variances are approximately $(1-\varphi )/{\tau }_b+\varphi /{\tau }_b=1/{\tau }_b$ (approximate because $Q^{\star -}$ equalises variances only in geometric mean).

Sparsity-preserving augmented parameterisation (Sec. 3.4)

The covariance in (8) contains $Q^{\star -}$, so its inverse is dense. Let $w=(w_1,w_2)$ with $w_2=u^{\star }$, $w_1=b$: $w_1\mid w_2\sim N(\sqrt{\varphi /{\tau }_b}{w}_2,\frac{1-\varphi }{{\tau }_b}I)$. Then $w$ has precision

$$\left(\begin{array}{cc}\frac{{\tau }_b}{1-\varphi }I& -\frac{\sqrt{\varphi {\tau }_b}}{1-\varphi }I\\ -\frac{\sqrt{\varphi {\tau }_b}}{1-\varphi }I& Q^{\star }+\frac{\varphi }{1-\varphi }I\end{array}\right),$$

sparse; $w_1$ has the correct marginal and $w_2$ is the structured component.

Examples

Scaling factor on a 4x4 grid

Compute the geometric-mean marginal variance of the Besag model on a 16-node grid graph and rescale: ${\sigma }_{GV}^2\approx 0.4716$; after scaling the geometric mean is 1 (verified).

import numpy as np
k = 4; n = k*k; A = np.zeros((n, n))
for i in range(k):
    for j in range(k):
        a = i*k + j
        if i+1 < k: A[a, a+k] = A[a+k, a] = 1
        if j+1 < k: A[a, a+1] = A[a+1, a] = 1
Q = np.diag(A.sum(1)) - A                       # Besag precision (connected graph, null space = 1)
c = np.exp(np.mean(np.log(np.diag(np.linalg.pinv(Q)))))   # sigma^2_GV for tau=1
Qs = Q * c                                      # scaled
print(c, np.exp(np.mean(np.log(np.diag(np.linalg.pinv(Qs))))))   # ~0.4716, 1.0

PyMC 5 BYM2 sketch (E18; untested; pm.ICAR requires the graph and a sum-to-zero soft constraint):

import pymc as pm
with pm.Model():
    mu = pm.Normal("mu", 0, 1)
    tau_sd = pm.Exponential("sigma", 1/0.31)          # PC prior on sd: P(sigma > U)=alpha  -> rate = -log(alpha)/U
    phi = pm.Beta("phi", 0.5, 0.5)                    # placeholder; PC prior for phi is tabulated (see PC note)
    v = pm.Normal("v", 0, 1, shape=n)
    u = pm.ICAR("u", W=A, sigma=1.0, zero_sum_stdev=0.001) / np.sqrt(c)   # scale: Q* = c Q -> sd(u*) = sd(u)/sqrt(c)
    b = tau_sd * (pm.math.sqrt(1 - phi) * v + pm.math.sqrt(phi) * u)
    pm.Poisson("y", mu=E * pm.math.exp(mu + b), observed=y)

Check: with $Q^{\star }=cQ$ and $c={\sigma }_{GV}^2$ from the code above, $\mathrm{Var}(u^{\star })=Q^{-}/c$, so a unit-scale ICAR draw is divided by $\sqrt{c}$, as in the sketch.

Marketing note: with geo-level lift estimates $b_g$ on an adjacency graph, BYM2 yields a single interpretable $\varphi$ (“how much of between-DMA heterogeneity is geographic”) portable across country maps.

Connections

Riebler 2016 - Penalised Complexity Priors for BYM2 puts priors on $({\tau }_b,\varphi )$; Spatial Models - BYM explains the unscaled baseline.

See Also

Riebler 2016 - BYM2 Intuitive Spatial Model - Overview