The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.