Tutorials › Generalized Linear Mixed Models (GLMM) › Proper CAR Models and Spatial Dependence

Spatial and Temporal Correlation · Tutorial 584 of 698

Proper CAR Models and Spatial Dependence

Learn how proper CAR conditionals combine into a nonsingular joint Gaussian model and how the adjacency matrix determines the allowable dependence parameter.

Advanced 12 min read

What You'll Learn

  • Specify proper CAR conditional means and variances using a symmetric adjacency matrix.
  • Derive the joint precision matrix from the conditional distributions.
  • Use normalized-adjacency eigenvalues to find the admissible range of the dependence parameter.
  • Distinguish a proper CAR model from its intrinsic boundary case.
  • Check numerical examples, including conditional and marginal variances.
  • Diagnose invalid adjacency structures and nearly singular precision matrices.

From Local Neighbourhoods to a Joint Spatial Model

A conditional autoregressive (CAR) model describes each area's latent effect through its neighbours. The conditional statements are local, but a spatial GLMM requires a coherent joint distribution for the entire vector of effects. In a proper CAR model, those conditionals imply a multivariate normal distribution with a positive-definite precision matrix. The adjacency structure determines which effects interact; the dependence parameter determines how strongly they interact and whether the joint distribution exists.

Let \(u=(u_1,\ldots,u_m)^\mathsf{T}\) be the latent effects for \(m\) areas. Let \(W=(w_{ij})\) be a fixed, symmetric, nonnegative adjacency matrix, with \(w_{ii}=0\). Usually \(w_{ij}=1\) when areas \(i\) and \(j\) are neighbours and \(0\) otherwise, though nonnegative weights can be used. Define \(d_i=\sum_j w_{ij}\) and \(D=\operatorname{diag}(d_1,\ldots,d_m)\). For the standard specification below, every area must have \(d_i>0\); an area with no neighbours needs separate treatment rather than division by zero.

Definition: A proper CAR model specifies, for each area \(i\), $$ u_i\mid u_{-i}\sim N\left(\rho\,\frac{\sum_{j\ne i}w_{ij}u_j}{d_i},\frac{\tau^2}{d_i}\right), $$ where \(u_{-i}\) denotes all effects except \(u_i\), \(\tau^2>0\) is a scale parameter, and \(\rho\) is the spatial-dependence parameter. The conditionals are proper when they are compatible with a nonsingular joint distribution; for this specification that requires the joint precision matrix to be positive definite.

The conditional mean is a weighted neighbour average multiplied by \(\rho\). At \(\rho=0\), the effects are independent under the joint model, with \(\operatorname{Var}(u_i)=\tau^2/d_i\). Thus, even at zero spatial dependence, areas with different numbers or total weights of neighbours have different marginal variances in this parameterization. The conditional variance \(\tau^2/d_i\) is not generally the marginal variance once \(\rho\ne0\).

Deriving the Joint Precision and Its Validity Condition

For a multivariate normal vector with precision matrix \(Q\), the conditional variance of component \(i\), given the other components, is \(1/Q_{ii}\), and its conditional mean is \(-Q_{ii}^{-1}\sum_{j\ne i}Q_{ij}u_j\). Matching these expressions to the CAR conditionals gives the only compatible symmetric precision matrix:

$$ Q=\frac{1}{\tau^2}(D-\rho W),\qquad u\sim N\left(0,Q^{-1}\right). $$

Indeed, \(Q_{ii}=d_i/\tau^2\) and \(Q_{ij}=-\rho w_{ij}/\tau^2\) for \(i\ne j\). Therefore the conditional variance is \(\tau^2/d_i\), and the conditional mean is \(\rho\sum_j w_{ij}u_j/d_i\), as specified. Symmetry of \(W\) matters: it ensures the resulting precision is symmetric and the local conditional specifications are compatible with one joint Gaussian distribution.

To find when this precision is positive definite, form the symmetric normalized adjacency matrix \(S=D^{-1/2}WD^{-1/2}\). Since \(D-\rho W=D^{1/2}(I-\rho S)D^{1/2}\), and \(D\) is positive definite when every \(d_i>0\), positive definiteness is equivalent to positive definiteness of \(I-\rho S\). Let the eigenvalues of \(S\) be \(\lambda_1,\ldots,\lambda_m\), with smallest and largest eigenvalues \(\lambda_{\min}\) and \(\lambda_{\max}\). The eigenvalues of \(I-\rho S\) are \(1-\rho\lambda_k\). Every one must be strictly positive.

$$ Q\text{ is positive definite} \quad\Longleftrightarrow\quad 1-\rho\lambda_k>0\ \text{for every }k \quad\Longleftrightarrow\quad \frac{1}{\lambda_{\min}}<\rho<\frac{1}{\lambda_{\max}}, $$

where the last interval applies when the graph has at least one edge, so that the eigenvalues include both a positive and a negative value. More generally, the unambiguous condition is \(1-\rho\lambda_k>0\) for every eigenvalue. For a symmetric, nonnegative adjacency matrix with no isolated areas, the largest eigenvalue of \(S\) is \(1\); hence a common restriction for positive spatial dependence is \(0\le\rho<1\). The complete admissible interval can extend below zero, and its lower endpoint depends on the most negative eigenvalue. The upper endpoint is excluded: at an endpoint, at least one eigenvalue of \(Q\) is zero, so the joint distribution is singular rather than proper.

Derivation: A symmetric matrix is positive definite exactly when all its eigenvalues are positive. The factorization \(D-\rho W=D^{1/2}(I-\rho S)D^{1/2}\) preserves positive definiteness because \(D^{1/2}\) is nonsingular. Diagonalizing \(S\) shows that \(I-\rho S\) has eigenvalues \(1-\rho\lambda_k\). Requiring all of them to be positive yields the parameter restriction. This eigenvalue check is the precise reason a dependence parameter cannot be chosen independently of the adjacency structure.

Worked Example: Conditional Distributions on a Three-Area Path

Hypothetical setup. Three areas form a path: area 1 neighbours area 2, and area 2 neighbours areas 1 and 3. Let \(\tau^2=4\) and \(\rho=0.5\). Then \(W=\begin{pmatrix}0&1&0\\1&0&1\\0&1&0\end{pmatrix}\) and \(D=\operatorname{diag}(1,2,1)\).

The conditional distributions are \(u_1\mid u_{-1}\sim N(0.5u_2,4)\), \(u_2\mid u_{-2}\sim N(0.25(u_1+u_3),2)\), and \(u_3\mid u_{-3}\sim N(0.5u_2,4)\). For example, if \(u_2=2\), the conditional mean of \(u_1\) is \(0.5(2)=1\); the conditional variance remains \(4\).

The joint precision is $$ Q=\frac14 \begin{pmatrix} 1&-0.5&0\\ -0.5&2&-0.5\\ 0&-0.5&1 \end{pmatrix}. $$ The normalized adjacency matrix has eigenvalues \(1,0,-1\), so the admissible interval is \(-1<\rho<1\). The chosen value \(0.5\) is valid. The determinant of \(D-\rho W\) is \(2(1-\rho^2)\), which at \(\rho=0.5\) is \(2(0.75)=1.5>0\); the eigenvalue check establishes positive definiteness, not just a positive determinant.

Worked Example: How the Graph Restricts the Parameter

Hypothetical setup. Three areas each neighbour the other two, so \(W\) has zeros on the diagonal and ones off the diagonal, and \(D=2I\). The normalized adjacency matrix is \(S=W/2\). Its eigenvalues are \(1,-0.5,-0.5\).

The eigenvalue conditions are \(1-\rho>0\) and \(1+0.5\rho>0\). Together these give \(-2<\rho<1\). If positive dependence is required, the model may instead impose \(0\le\rho<1\). For \(\rho=1.2\), the condition associated with eigenvalue \(1\) is \(1-1.2=-0.2<0\); the precision is not positive definite, so these conditionals do not define a proper joint Gaussian model. At \(\rho=1\), that eigenvalue is zero, making the precision singular.

This example also shows why a blanket lower bound such as \(\rho>-1\) is not generally the full admissible interval. The graph's eigenvalues determine the lower bound; the upper bound is \(1\) here because \(\lambda_{\max}=1\).

Marginal Dependence Is Not the Conditional Dependence

The off-diagonal entries of \(Q\) encode conditional relationships, but marginal covariances are entries of \(Q^{-1}\). These are related but are not interchangeable. In particular, a conditional variance is not a marginal variance, and a conditional mean that uses only immediate neighbours does not imply that marginal dependence stops at those neighbours.

Worked Example: Recovering Marginal Covariances

Hypothetical setup. Use the three-area path from the first example, but take \(\tau^2=1\) and \(\rho=0.5\). Its precision matrix is $$ Q= \begin{pmatrix} 1&-0.5&0\\ -0.5&2&-0.5\\ 0&-0.5&1 \end{pmatrix}. $$ The determinant is \(1.5\), and direct matrix inversion gives $$ Q^{-1}=\frac{1}{1.5} \begin{pmatrix} 1.75&0.5&0.25\\ 0.5&1&0.5\\ 0.25&0.5&1.75 \end{pmatrix}. $$

Thus \(\operatorname{Var}(u_1)=1.75/1.5=1.167\), while the conditional variance of \(u_1\) given its neighbours is \(1\). Also, \(\operatorname{Cov}(u_1,u_3)=0.25/1.5=0.1667\), although areas 1 and 3 are not adjacent. Their nonzero marginal covariance arises through their shared connection to area 2. The precision matrix has a zero in position \((1,3)\), indicating conditional independence of \(u_1\) and \(u_3\) given \(u_2\) in this Gaussian model, not marginal independence.

Assumptions, Diagnostics, and Alternatives

The standard derivation rests on a fixed, symmetric adjacency matrix with nonnegative weights, positive row sums, a common \(\rho\), and a common scale \(\tau^2>0\) in the stated conditional-variance parameterization. Before fitting a model, check that the area ordering in \(W\) matches the ordering of \(u\), that the weights are symmetric, and that no row has sum zero. If weights are asymmetric, the conditionals as written do not produce the symmetric precision matrix derived here; symmetrizing weights changes the spatial structure and should be a modeling decision, not an automatic repair.

A practical validity check is to compute the eigenvalues of \(S\) and verify \(1-\rho\lambda_k>0\) for every \(k\). Also inspect the smallest eigenvalue of \(Q\), or its condition number: values close to zero indicate a nearly singular precision and can cause unstable computation and large marginal variances. A fitted \(\rho\) near its permitted boundary is a reason to check sensitivity to the adjacency specification and to examine whether the data support a boundary-near dependence estimate. An isolated area requires an explicit alternative, such as an independent effect for that area; silently assigning it a zero conditional variance is not valid.

A proper CAR model is useful when adjacency is substantively meaningful, as for areas sharing borders. A spatial covariance model based on distances instead encodes dependence through locations and a covariance function; it is a natural alternative when continuous separation is central. The intrinsic CAR model is another adjacency-based alternative: it lies at a singular boundary and needs special treatment for its null space, rather than being a proper joint Gaussian distribution. Independent area effects are simpler, but omit spatial dependence. These alternatives answer different structural assumptions; none is a universal correction for misspecified adjacency.

A proper CAR model also has limitations. Its local conditional specification does not guarantee that adjacency captures the scientifically relevant spatial relation, and the single parameter \(\rho\) may be too restrictive when dependence varies across space or direction. The marginal variances depend on the whole graph and on \(\rho\), so interpreting \(\tau^2\) as a marginal variance is generally wrong. Finally, a proper prior distribution for the spatial effects does not by itself ensure that the GLMM likelihood contains enough information to estimate all fixed effects and covariance parameters reliably.

Key takeaway: For the standard proper CAR conditionals, the joint precision is \(Q=\tau^{-2}(D-\rho W)\). The model is proper exactly when this matrix is positive definite, equivalently when \(1-\rho\lambda_k>0\) for every eigenvalue of \(D^{-1/2}WD^{-1/2}\). Check the graph and the eigenvalue restriction before interpreting or fitting the model.

Check Your Understanding

For each question, work through the calculation before reading the solution.

  1. For a three-area path with \(\tau^2=9\) and \(\rho=0.4\), what is the conditional distribution of \(u_2\) given \(u_1=2\) and \(u_3=-1\)?
    Solution: The centre area has degree \(2\), so its conditional mean is \(0.4(2+(-1))/2=0.2\), and its conditional variance is \(9/2=4.5\). Hence \(u_2\mid u_1,u_3\sim N(0.2,4.5)\).
  2. For the equal-weight three-area complete graph, is \(\rho=-2\) proper?
    Solution: The normalized adjacency eigenvalues are \(1,-0.5,-0.5\). At \(\rho=-2\), \(1-\rho(-0.5)=1-1=0\). The precision is singular, so \(\rho=-2\) is not proper; the strict admissible interval is \((-2,1)\).
  3. In the three-area path with \(\rho=0.5\), are areas 1 and 3 marginally independent?
    Solution: No. The inverse precision matrix has \(\operatorname{Cov}(u_1,u_3)=0.25/1.5=0.1667\) when \(\tau^2=1\), so they are marginally dependent. They are conditionally independent given \(u_2\), because their off-diagonal precision entry is zero.
  4. What should be checked if the adjacency matrix is asymmetric or contains a row whose sum is zero?
    Solution: The stated derivation assumes symmetric weights and positive row sums. Asymmetry breaks the symmetric-precision argument, while a zero row sum makes the conditional mean and variance formulas undefined. Revisit the adjacency specification or model the affected area separately; do not apply the standard proper CAR formula unchanged.
Clinical Biostats