BYM Decomposition of Spatial Effects

Why Decompose an Areal Effect?

In an areal GLMM, observations belong to geographic units such as districts, counties, or service regions. After accounting for measured predictors, outcomes may still vary across areas. Some of this variation may be spatially patterned: neighboring areas tend to have similar latent effects. Other variation may be specific to an area and unrelated to its neighbors. The BYM decomposition represents both possibilities in one area-level effect.

The decomposition is useful because a single spatial effect can blur two distinct patterns. A strongly varying effect need not mean that the map has strong spatial structure: it could instead reflect substantial area-specific noise. Conversely, smooth spatial variation can be important even when the area-specific noise is modest. The model separates these components, but their variance parameters require careful interpretation.

Model Definition and Assumptions

Let $i=1,\ldots,m$ index areas and $j$ index observations within area $i$. Write the conditional GLMM linear predictor as

g\{\operatorname{E}(Y_{ij}\mid u_i,v_i)\} = x_{ij}^{\mathsf T}\beta + u_i+v_i .

Here $x_{ij}$ is the covariate vector, $\beta$ is the fixed-effect coefficient vector, $u_i$ is the structured spatial effect, and $v_i$ is the unstructured area effect. The response distribution and link $g$ are specified as part of the GLMM. Conditional on all random effects and covariates, observations are assumed to follow the chosen response model; the area effects enter through the linear predictor.

Definition: In the BYM decomposition, $u=(u_1,\ldots,u_m)^{\mathsf T}$ is an intrinsic conditional autoregressive (ICAR) effect on the area adjacency graph, and $v_i$ are independent Gaussian effects: $$ u\sim\text{ICAR}(\tau_u),\qquad v_i\overset{\text{ind}}{\sim}N(0,\tau_v^{-1}),\qquad u\ \perp\ v. $$ The BYM area effect is $b_i=u_i+v_i$. The precision parameters $\tau_u,\tau_v$ are positive; their reciprocals are variance parameters.

For a symmetric adjacency matrix $W$, let $w_{ik}=1$ when areas $i$ and $k$ are neighbors and $0$ otherwise, and let $d_i=\sum_k w_{ik}$. The ICAR precision matrix is $Q=D-W$, where $D$ is diagonal with entries $d_i$. For a connected graph without isolated areas, its conditional distributions are

u_i\mid u_{-i}\sim N\left(\frac{1}{d_i}\sum_{k:w_{ik}=1}u_k,\, \frac{1}{\tau_u d_i}\right).

The corresponding joint density is proportional to $\exp(-\tau_u u^{\mathsf T}Qu/2)$. Since $Q$ has a constant-vector null space, this is not a proper density over all of $\mathbb{R}^m$. A common identification constraint is $\sum_i u_i=0$. For a graph with multiple connected components, a separate sum-to-zero constraint is needed within each component. Isolated areas require separate treatment because their ICAR conditional distributions are not defined.

The usual BYM assumptions are that the adjacency graph is specified meaningfully, the structured and unstructured components are independent a priori, and both have mean zero under the chosen constraint. Conditional outcome distributions must also be suitable for the data. The random effects describe residual area variation; they do not remove the need to check the fixed-effects specification or the response model.

How the Two Components Contribute

Under the sum-to-zero constraint, the ICAR covariance on the constrained subspace is $\tau_u^{-1}Q^{+}$, where $Q^{+}$ is the Moore–Penrose pseudoinverse. Because $u$ and $v$ are independent, the covariance of $b=u+v$ is the sum of their covariance matrices:

\operatorname{Cov}(b) =\tau_u^{-1}Q^{+}+\tau_v^{-1}I.

Thus the unstructured component contributes the same variance $\tau_v^{-1}$ to every area and no covariance between distinct areas. The structured component contributes area-specific marginal variances $\tau_u^{-1}(Q^{+})_{ii}$ and spatial covariances $\tau_u^{-1}(Q^{+})_{ik}$. Those structured marginal variances generally differ by area and depend on the graph. In particular, $\tau_u^{-1}$ is a precision-scaled ICAR parameter, not the marginal variance of every $u_i$.

This distinction is central to interpretation. A comparison of $\tau_u^{-1}$ and $\tau_v^{-1}$ alone does not, in general, give a clean proportion of total area-level variation attributable to spatial structure. The graph affects the marginal variances of $u$, and the ICAR constraint also induces dependence. Interpret estimated component parameters in relation to the specified graph and any scaling or reparameterization used.

Worked Example: Total Variance on a Three-Area Chain

Consider three hypothetical areas arranged in a chain: area 1 neighbors area 2, and area 2 neighbors areas 1 and 3. Suppose $\tau_u=4$ and $\tau_v=1$. The graph is connected, so impose $\sum_i u_i=0$. For this graph,

Q^{+}=\frac{1}{9} \begin{pmatrix} 5&-1&-4\\ -1&2&-1\\ -4&-1&5 \end{pmatrix}.

The structured covariance is $\tau_u^{-1}Q^{+}=Q^{+}/4$. Its diagonal entries are $5/36,\ 2/36,\ 5/36$, or approximately $0.1389,\ 0.05556,\ 0.1389$. The unstructured variance is $\tau_v^{-1}=1$ in every area. Therefore the total marginal variances are

\begin{aligned} \operatorname{Var}(b_1)&=\frac{5}{36}+1=1.1389,\\ \operatorname{Var}(b_2)&=\frac{2}{36}+1=1.0556,\\ \operatorname{Var}(b_3)&=\frac{5}{36}+1=1.1389. \end{aligned}

The structured covariance between areas 1 and 2 is $-1/36$, and between areas 1 and 3 it is $-1/9$; the unstructured component adds no off-diagonal covariance. These negative values are compatible with the sum-to-zero constraint, not evidence that the two components are negatively associated. This calculation also shows why the structured precision parameter alone is not the marginal variance of an area effect.

A Short Derivation: Why Covariances Add

The covariance calculation follows directly from the model definition. For any two areas $i$ and $k$, $b_i=u_i+v_i$ and $b_k=u_k+v_k$. Expanding their covariance gives

\operatorname{Cov}(b_i,b_k) =\operatorname{Cov}(u_i,u_k) +\operatorname{Cov}(u_i,v_k) +\operatorname{Cov}(v_i,u_k) +\operatorname{Cov}(v_i,v_k).

Independence of the two components makes the middle terms zero. The independent $v_i$ effects have covariance zero for $i\ne k$, and variance $\tau_v^{-1}$ for $i=k$. The remaining term is the ICAR covariance. In matrix form this yields $\tau_u^{-1}Q^++\tau_v^{-1}I$. This derivation depends on the independence assumption; a model that correlates its two components would have additional cross-covariance terms.

Worked Interpretation in a Count GLMM

For a Poisson model with a log link and an exposure offset $E_i$, one possible specification is

Y_i\mid u_i,v_i\sim\operatorname{Poisson}(\mu_i),\qquad \log(\mu_i)=\log(E_i)+\beta_0+\beta_1x_i+u_i+v_i.

Conditional on the area effects, the expected count is $\mu_i$; the expected rate per exposure unit is $\exp(\beta_0+\beta_1x_i+u_i+v_i)$. The next example uses hypothetical parameter values to illustrate the decomposition, not to describe an observed study.

Worked Example: Combining Fixed and Area Effects

Suppose a hypothetical area has exposure $E_i=100$, covariate $x_i=1$, $\beta_0=-2$, $\beta_1=0.4$, structured effect $u_i=0.2$, and unstructured effect $v_i=0.1$. Its log expected count is

\log(\mu_i)=\log(100)-2+0.4(1)+0.2+0.1 =4.6052-1.3=3.3052.

Thus $\mu_i=\exp(3.3052)\approx27.25$ expected events, or an expected rate of $27.25/100=0.2725$ events per exposure unit. The combined latent area effect is $b_i=0.2+0.1=0.3$, which multiplies the conditional rate by $\exp(0.3)\approx1.350$ relative to the same covariate and fixed effects with $b_i=0$. This is a conditional model interpretation; it is not a causal effect and does not describe a marginal rate averaged over the random-effect distribution.

Conditional Behavior, Diagnostics, and Limitations

The conditional ICAR mean is a neighbor average, while its conditional variance decreases as an area has more neighbors: $\operatorname{Var}(u_i\mid u_{-i})=1/(\tau_u d_i)$. That formula describes the structured component conditional on the other structured effects. It is not the variance of $b_i$, which also includes $v_i$, and it does not say that neighboring areas must have identical effects.

Worked Example: Reading a Local Conditional Distribution

Suppose a hypothetical area has two neighbors with structured effects $-0.4$ and $0.2$, and $\tau_u=4$. Its conditional structured-effect mean and variance are

\operatorname{E}(u_i\mid u_{-i})=\frac{-0.4+0.2}{2}=-0.1,\qquad \operatorname{Var}(u_i\mid u_{-i})=\frac{1}{4(2)}=0.125.

The conditional standard deviation is $\sqrt{0.125}\approx0.3536$. If this area also has $v_i=0.15$, its combined effect at this particular value of $u_i$ is $u_i+v_i$; the value $v_i$ does not change the conditional mean or variance formula for $u_i$. The example concerns the prior conditional structure, not a posterior estimate based on outcomes.

In an applied analysis, inspect whether residual spatial pattern remains after fitting, whether a small number of areas dominate the estimated spatial component, and whether estimates change substantially under plausible graph or prior choices. Check for isolated areas, disconnected graph components, and sensitivity to constraints. Poor fit of the response distribution, omitted covariates, or an inappropriate adjacency definition can all be mistaken for a need for stronger spatial structure.

The two components can be difficult to distinguish from the data: a spatially smooth realization of independent effects can resemble a structured effect, and a structured field can include local variation. The decomposition does not guarantee precise estimation of both variance parameters. In addition, with a constrained ICAR field, component variances can be sensitive to graph size and shape. Report the adjacency definition, constraints, parameterization, and uncertainty rather than presenting the two precision estimates as universally comparable shares of variance.

How BYM Differs From Common Alternatives

Model for area variation	What it represents	Key distinction from BYM
Independent Gaussian area effects only	Area-specific variation with no adjacency-based dependence	Does not represent spatial smoothing across neighbors.
ICAR effect only	Structured variation linked to the adjacency graph	Has no separate independent Gaussian component for local noise.
BYM decomposition	Structured ICAR variation plus independent Gaussian variation	Allows both patterns, but their contributions can be confounded and graph-dependent.
Proper CAR effect	A spatially dependent effect with a proper joint distribution under its parameter conditions	Uses a different covariance specification; it is not simply the same as the BYM sum.

The right choice depends on the scientific meaning of the graph, the desired spatial covariance structure, and whether the data can distinguish the components. BYM does not make those choices automatic. It is also not a substitute for an appropriate likelihood or for checking the fixed-effects part of the GLMM.

Common Mistake

Wrong reasoning: “The structured variance is $1/\tau_u$, the unstructured variance is $1/\tau_v$, so their ratio directly gives the fraction of area-level variation that is spatial.”
Correction: $1/\tau_v$ is the marginal variance of each $v_i$, but the ICAR marginal variance in area $i$ is $(Q^+)_{ii}/\tau_u$, which depends on the graph. Compare appropriate marginal variance summaries for the stated graph and parameterization, and include uncertainty.
Wrong reasoning: “The BYM sum-to-zero constraint means the observed area effects must add to zero.”
Correction: The constraint applies to the structured component $u$. The independent component $v$, and therefore the sum $b=u+v$, is not constrained to sum to zero.

Key takeaway: BYM writes each area effect as a graph-structured ICAR term plus independent Gaussian noise. The components are conceptually distinct, but their marginal variance contributions depend on the graph, constraints, and parameterization.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

If $\tau_v=2$, what is the marginal variance of each unstructured effect $v_i$?
Why does the BYM model require a constraint or another treatment of the ICAR null space?
In the three-area chain example, what is the total covariance between areas 1 and 2?
What does a positive $v_i$ represent in a GLMM linear predictor?
Name one limitation of interpreting the BYM precision parameters as direct proportions of total spatial variation.

Solutions. (1) $\operatorname{Var}(v_i)=1/\tau_v=1/2=0.5$. (2) The ICAR precision matrix is singular because a constant shift in all $u_i$ has zero ICAR penalty; a sum-to-zero constraint on a connected graph identifies the structured field. (3) It is the structured covariance $-1/36$ plus zero unstructured covariance, hence $-1/36\approx-0.02778$. (4) Holding the other terms fixed, it adds $v_i$ on the link scale; for a log link, a value $v_i$ multiplies the conditional mean by $\exp(v_i)$. (5) The structured marginal variances depend on the graph through the diagonal of $Q^+$, so $1/\tau_u$ is not a common area-level marginal variance directly comparable to $1/\tau_v$.