Sampling from Gaussian CAR and ICAR Models

Purpose

This vignette demonstrates how to sample Gaussian random effects defined on a graph using trafficCAR. The emphasis is on models specified through a sparse precision matrix:

\[ x \sim \mathrm{N}(Q^{-1} b,\; Q^{-1}), \]

where \(Q\) is a sparse precision matrix and \(b\) is an optional linear term. Both proper CAR and intrinsic CAR (ICAR) specifications are illustrated.

library(trafficCAR)
library(Matrix)

CAR and ICAR precision matrices

Let \(A\) be a symmetric adjacency matrix with zero diagonal, and let \(D = \mathrm{diag}(A\mathbf{1})\). A common CAR family is

\[ Q = \tau (D - \rho A), \]

with \(\tau > 0\). For a proper CAR, choose \(\rho\) so that \(D - \rho A\) is positive definite. For an ICAR, the conventional choice is \(\rho = 1\), which yields a singular precision matrix.

The function car_precision() constructs \(Q\) from \(A\).

A small graph example

Construct a simple path graph with 3 nodes.

A_small <- matrix(
  c(0, 1, 0,
    1, 0, 1,
    0, 1, 0),
  nrow = 3,
  byrow = TRUE
)
A_small
#>      [,1] [,2] [,3]
#> [1,]    0    1    0
#> [2,]    1    0    1
#> [3,]    0    1    0

Proper CAR

Q_car <- car_precision(A_small, type = "proper", rho = 0.5, tau = 1)
Q_car
#> 3 x 3 sparse Matrix of class "dsCMatrix"
#>                    
#> [1,]  1.0 -0.5  .  
#> [2,] -0.5  2.0 -0.5
#> [3,]  .   -0.5  1.0

ICAR

Q_icar <- car_precision(A_small, type = "icar", tau = 1)
Q_icar
#> 3 x 3 sparse Matrix of class "dsCMatrix"
#>              
#> [1,]  1 -1  .
#> [2,] -1  2 -1
#> [3,]  . -1  1

Sampling from a sparse precision matrix

Sampling is performed by rmvnorm_prec(), which draws from \(\mathcal{N}(Q^{-1} b, Q^{-1})\) using sparse factorizations. A dense covariance matrix is not formed.

Zero-mean sampling

x0 <- trafficCAR:::rmvnorm_prec(Q_car)
x0
#> [1] -1.90422144 -0.06060391  0.06454835

Nonzero mean (linear term)

b <- c(1, 0, -1)
x1 <- trafficCAR:::rmvnorm_prec(Q_car, b = b)
x1
#> [1]  2.6066965  0.3768962 -0.5363042

Notes on ICAR sampling

For ICAR models, \(Q\) is singular and the implied Gaussian distribution is improper. Samples are identifiable only up to an additive constant. A common post-processing step is to center a draw.

# ICAR precision matrices are singular, so direct Cholesky-based sampling may fail.
# A practical approach for simulation is to add a small ridge and then center.
# This yields an approximate draw on the ICAR subspace.

n_icar <- nrow(Q_icar)
eps <- 1e-6
Q_icar_ridge <- Q_icar + Matrix::Diagonal(n_icar, x = eps)

x_icar <- trafficCAR:::rmvnorm_prec(Q_icar_ridge)
x_icar_centered <- x_icar - mean(x_icar)
x_icar_centered
#> [1] -0.6858918  0.3233726  0.3625193

A road-like graph example

This section illustrates CAR/ICAR sampling on a graph resembling a small road network. If sf is available, we build an adjacency matrix from the bundled roads_small dataset. Otherwise, we fall back to a small synthetic graph.

A_road <- NULL

has_sf <- requireNamespace("sf", quietly = TRUE)

if (has_sf) {
  data("roads_small", package = "trafficCAR")
  roads <- roads_small

  segments <- roads_to_segments(
    roads,
    crs_m = 3857,
    split_at_intersections = TRUE
  )

  if (nrow(segments) > 200) {
    segments <- segments[seq_len(200), ]
  }

  adjacency <- build_adjacency(segments, crs_m = 3857)

  # Drop isolated segments to keep CAR models well-defined.
  if (any(adjacency$isolates)) {
    segments <- segments[!adjacency$isolates, ]
    adjacency <- build_adjacency(segments, crs_m = 3857)
  }

  A_road <- adjacency$A
}

Option 2: Synthetic road-like graph (fallback)

The following graph mimics a small road network with a T-junction.

if (is.null(A_road)) {
  # nodes: 1--2--3 and 2--4 (T-junction at node 2)
  edges <- matrix(c(
    1, 2,
    2, 3,
    2, 4
  ), ncol = 2, byrow = TRUE)

  if (requireNamespace("igraph", quietly = TRUE)) {
    g <- igraph::graph_from_edgelist(edges, directed = FALSE)
    A_road <- igraph::as_adjacency_matrix(g, sparse = TRUE, attr = NULL)
  } else {
    # minimal fallback without igraph (dense, small only)
    n <- max(edges)
    A_road <- matrix(0, n, n)
    for (k in seq_len(nrow(edges))) {
      i <- edges[k, 1]; j <- edges[k, 2]
      A_road[i, j] <- 1
      A_road[j, i] <- 1
    }
    A_road <- Matrix(A_road, sparse = TRUE)
  }
}

c(n = nrow(A_road), nnz = Matrix::nnzero(A_road))
#>   n nnz 
#> 183 364

CAR/ICAR sampling on the road graph

# proper CAR (users should choose rho consistent with their graph and model)
Q_car_road <- car_precision(A_road, type = "proper", rho = 0.5, tau = 1)

# ICAR
Q_icar_road <- car_precision(A_road, type = "icar", tau = 1)

x_car_road <- trafficCAR:::rmvnorm_prec(Q_car_road)

n_road <- nrow(Q_icar_road)
eps <- 1e-6
Q_icar_road_ridge <- Q_icar_road + Matrix::Diagonal(n_road, x = eps)
x_icar_road <- trafficCAR:::rmvnorm_prec(Q_icar_road_ridge)
x_icar_road <- x_icar_road - mean(x_icar_road)

summary(x_car_road)
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -2.22146 -0.48870  0.03991  0.03359  0.57647  2.30052
summary(x_icar_road)
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -1381.57   -60.11   100.95     0.00   247.66   940.77

Summary