Conditional autoregressive (CAR) models

spatial

regression

autocorrelation

areal data

Build areal spatial dependence from local conditional distributions, fit a proper CAR by hand, contrast it with SAR, and use the intrinsic CAR as a smoother.

Author

Tidy Ecology

Published

2026-08-13

The simultaneous autoregressive models build spatial dependence from the top down: one global equation, y = rho W y + X beta + eps, solved all at once. Conditional autoregressive (CAR) models take the opposite route. They specify, for each cell, only its distribution given its neighbours, and let a theorem stitch those local pieces into a coherent joint distribution. This bottom-up construction, due to Besag (1974), is the engine behind most Bayesian disease mapping and behind the spatial random effects in ecological hierarchical models. It is worth understanding directly.

This tutorial writes the conditional model down, assembles its precision matrix, fits a proper CAR by hand in base R, and then does three things the textbooks often skip: it shows why the CAR autocorrelation parameter is not comparable across software packages, it contrasts the covariance a CAR implies with the one a SAR implies, and it uses the intrinsic CAR as a practical smoother.

The conditional specification

A CAR model says that, conditional on every other cell, cell i is normal with a mean pulled towards its neighbours:

\[ y_i \mid y_{-i} \sim \mathcal{N}\!\left(\mu_i + \rho \sum_j \frac{a_{ij}}{d_i}(y_j - \mu_j),\ \frac{\tau^2}{d_i}\right), \]

where a_ij marks neighbours, d_i is the number of neighbours of cell i, and mu = X beta. The Hammersley-Clifford theorem guarantees these conditionals correspond to a single joint Gaussian, and its precision matrix (the inverse covariance) is wonderfully sparse:

\[ Q = \frac{1}{\tau^2}(D - \rho A), \]

with D the diagonal matrix of neighbour counts and A the adjacency matrix. The model is proper (a valid distribution) when D - rho A is positive definite, which pins rho inside a range set by the eigenvalues.

Setting up the lattice

side <- 20; n <- side * side
coords <- expand.grid(col = 1:side, row = 1:side)
A <- matrix(0, n, n)
idx <- function(r, c) (c - 1) * side + r
for (r in 1:side) for (c in 1:side) {
  i <- idx(r, c)
  if (r > 1)    A[i, idx(r - 1, c)] <- 1
  if (r < side) A[i, idx(r + 1, c)] <- 1
  if (c > 1)    A[i, idx(r, c - 1)] <- 1
  if (c < side) A[i, idx(r, c + 1)] <- 1
}
deg <- rowSums(A); D <- diag(deg)
# eigenvalues of the row-standardised operator fix the admissible rho range
ev <- eigen(diag(1 / sqrt(deg)) %*% A %*% diag(1 / sqrt(deg)), symmetric = TRUE)$values
rho_lo <- 1 / min(ev); rho_hi <- 1 / max(ev)
c(lower = rho_lo, upper = rho_hi)

lower upper 
   -1     1

Simulating a proper CAR field

We draw one realisation with an intercept, a covariate, a strong conditional dependence of rho = 0.9, and a conditional scale of tau = 0.4. Sampling is direct from the joint covariance tau^2 (D - rho A)^{-1}.

set.seed(4231)
cc <- scale(coords$col)[, 1]
x  <- as.numeric(scale(0.8 * cos(pi * coords$row / 7) + 0.6 * cc + rnorm(n, 0, 0.6)))
X  <- cbind(1, x)
b0 <- 0.5; b1 <- 0.7; rho <- 0.9; tau <- 0.4
Q  <- (D - rho * A) / tau^2
set.seed(11)
y  <- as.numeric(X %*% c(b0, b1) + t(chol(solve(Q))) %*% rnorm(n))

Fitting a CAR by hand

The log-likelihood needs the log-determinant of the precision. Factor it as |D - rho A| = |D| * |I - rho D^{-1}A|, so the determinant reduces to a sum over the same eigenvalues we already have. For each candidate rho we solve the generalised least squares problem for beta, read off tau from the quadratic form, and search over the single remaining parameter.

nll_car <- function(r) {
  if (r <= rho_lo || r >= rho_hi) return(1e10)
  P  <- D - r * A                               # precision up to 1/tau^2
  be <- solve(t(X) %*% P %*% X, t(X) %*% P %*% y)
  e  <- y - X %*% be; t2 <- as.numeric(t(e) %*% P %*% e) / n
  ldetP <- sum(log(deg)) + sum(log(1 - r * ev)) # log|D - rho A|
  -(-n / 2 * log(2 * pi) - n / 2 - n / 2 * log(t2) + 0.5 * ldetP)
}
opt <- optimise(nll_car, c(rho_lo + 1e-4, rho_hi - 1e-4))
car_rho <- opt$minimum
P  <- D - car_rho * A
be <- solve(t(X) %*% P %*% X, t(X) %*% P %*% y)
car_b0 <- be[1]; car_b1 <- be[2]
car_tau <- sqrt(as.numeric(t(y - X %*% be) %*% P %*% (y - X %*% be)) / n)
car_ll  <- -opt$objective
round(c(rho = car_rho, tau = car_tau, b0 = car_b0, b1 = car_b1), 4)

   rho    tau     b0     b1 
0.8184 0.3929 0.4774 0.6893

The fit recovers the generating values well: a conditional dependence of 0.82 against the true 0.9, a scale of 0.39 against 0.4, and a slope of 0.69 against 0.7. The likelihood machinery is the same log-determinant plus generalised least squares we used for the error model in the SAR tutorial, only the precision differs.

A warning about conventions

Reach for a package and you may be surprised. The spautolm function in spdep, with family = "CAR", does not use the neighbour-count precision D - rho A. It uses a symmetric-weights convention with precision I - lambda A. That is a different model, so its lambda is not on the same scale as our rho and the two numbers will never agree. To show this is a convention and not an error, we can fit the symmetric version by hand as well.

evA <- eigen(A, symmetric = TRUE)$values
loA <- 1 / min(evA); hiA <- 1 / max(evA)
nll_sym <- function(lam) {
  if (lam <= loA || lam >= hiA) return(1e10)
  P  <- diag(n) - lam * A
  be <- solve(t(X) %*% P %*% X, t(X) %*% P %*% y); e <- y - X %*% be
  s2 <- as.numeric(t(e) %*% P %*% e) / n
  -(-n / 2 * log(2 * pi) - n / 2 - n / 2 * log(s2) + 0.5 * sum(log(1 - lam * evA)))
}
car_sym_lambda <- optimise(nll_sym, c(loA + 1e-6, hiA - 1e-6))$minimum
round(car_sym_lambda, 4)

[1] 0.2122

The symmetric convention reports 0.212, and this by-hand value matches spdep’s spautolm to numerical precision. It describes the same data as our 0.82, just on a different scale. The lesson is practical: never compare a CAR autocorrelation parameter across papers or packages without first checking which weights convention each one uses.

CAR is not SAR in disguise

Both models encode neighbourhood dependence, but they imply different covariances, and neither is a reparametrisation of the other except in special cases (Ver Hoef et al. 2018). A useful way to see the difference is to trace how the correlation between a central cell and the rest decays with graph distance, the number of steps through the lattice.

centre <- idx(10, 10)
gdist  <- rep(NA_integer_, n); gdist[centre] <- 0L; step <- 0L; front <- centre
repeat {
  nb <- setdiff(unique(unlist(lapply(front, function(i) which(A[i, ] > 0)))),
                which(!is.na(gdist)))
  if (!length(nb)) break
  step <- step + 1L; gdist[nb] <- step; front <- nb
}
Ccar <- cov2cor(solve(D - 0.9 * A))[centre, ]                       # CAR reach
Wrs  <- A / deg
Csar <- cov2cor(solve(crossprod(diag(n) - 0.6 * Wrs)))[centre, ]    # SAR reach
reach <- aggregate(data.frame(CAR = Ccar, SAR = Csar), list(dist = gdist), mean)
round(head(reach, 5), 3)

  dist   CAR   SAR
1    0 1.000 1.000
2    1 0.346 0.350
3    2 0.168 0.132
4    3 0.086 0.046
5    4 0.045 0.016

At the first step the two are tuned to match, near 0.35. Further out they part company: at three steps the CAR retains a correlation of 0.086 against the SAR’s 0.046. The conditional model here carries dependence further across the map than the simultaneous one at matched short range. Which shape suits your system is an empirical question, and one more reason the two model families are worth keeping distinct.

Two panels. The left panel is a line chart of correlation against graph distance for CAR and SAR, both starting at one and decaying, with the CAR line staying above the SAR line beyond the first step. The right panel is a grid map of a single simulated proper-CAR field showing smooth spatial patches.

The intrinsic CAR as a smoother

Push rho all the way to its boundary and the precision becomes D - A, the graph Laplacian. This is the intrinsic CAR (ICAR): an improper distribution (it has no overall level, hence the usual sum-to-zero constraint) that penalises differences between neighbours. It is the most widely used spatial random effect in ecology and epidemiology, because as a smoother it borrows strength across the map. To see it work, we take a noisy signal, one observation per cell, and recover the underlying smooth field.

set.seed(4232)
truef <- as.numeric(scale(1.4 * sin(pi * coords$col / 8) * cos(pi * coords$row / 9)))
yn    <- truef + rnorm(n, 0, 1.0)          # noisy observation per cell
R_icar <- D - A                            # intrinsic precision (graph Laplacian)
# GMRF smoother S = (I + kappa R)^{-1}; choose kappa by generalised cross-validation
gcv <- function(kappa) {
  S <- solve(diag(n) + kappa * R_icar); fit <- S %*% yn
  n * sum((yn - fit)^2) / (n - sum(diag(S)))^2
}
ks <- 10^seq(-2, 2, length = 40); k_opt <- ks[which.min(sapply(ks, gcv))]
Sopt   <- solve(diag(n) + k_opt * R_icar)
fit_ic <- as.numeric(Sopt %*% yn)
icar_kappa    <- k_opt
icar_edf      <- sum(diag(Sopt))
icar_rmse_raw <- sqrt(mean((yn - truef)^2))
icar_rmse_sm  <- sqrt(mean((fit_ic - truef)^2))
round(c(kappa = icar_kappa, edf = icar_edf,
        rmse_raw = icar_rmse_raw, rmse_smooth = icar_rmse_sm), 3)

      kappa         edf    rmse_raw rmse_smooth 
      0.438     171.406       0.894       0.422

Cross-validation settles on a smoother with about 171 effective degrees of freedom out of 400 cells. It cuts the root mean squared error against the truth from 0.89 for the raw observations to 0.42, roughly halving it, simply by letting each cell learn from its neighbours. This borrowing of strength is exactly what an ICAR random effect does inside a larger model.

Two grid maps side by side. The left panel shows noisy per-cell observations with a speckled appearance. The right panel shows the ICAR-smoothed field, which reveals a clean smooth spatial pattern with the noise removed.

What to carry forward

A CAR model is built from local conditional rules and summarised by a sparse precision matrix D - rho A. It is the natural spatial term for hierarchical and generalised linear mixed models, where the conditional form slots neatly into the wider structure. Three points are worth remembering. The autocorrelation parameter depends on the weights convention, so it is not portable between software. A CAR and a SAR on the same neighbourhood imply genuinely different covariances, and the choice between them is substantive, not cosmetic. And the intrinsic CAR, the boundary case, is less a model of correlation than a smoother, penalising neighbour differences to pool information.

That smoothing power comes with a hidden cost. When the spatial random effect is flexible enough to mimic a spatially structured covariate, it can quietly absorb that covariate’s effect. This is spatial confounding, and checking a spatial regression takes it apart in detail.

References

Besag, J. (1974). Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society: Series B, 36(2), 192-236. https://doi.org/10.1111/j.2517-6161.1974.tb00999.x

Dormann, C. F., McPherson, J. M., Araujo, M. B., Bivand, R., Bolliger, J., Carl, G., et al. (2007). Methods to account for spatial autocorrelation in the analysis of species distributional data: a review. Ecography, 30(5), 609-628. https://doi.org/10.1111/j.2007.0906-7590.05171.x

Ver Hoef, J. M., Peterson, E. E., Hooten, M. B., Hanks, E. M., & Fortin, M.-J. (2018). Spatial autoregressive models for statistical inference from ecological data. Ecological Monographs, 88(1), 36-59. https://doi.org/10.1002/ecm.1283

--- title: "Conditional autoregressive (CAR) models" description: "Build areal spatial dependence from local conditional distributions, fit a proper CAR by hand, contrast it with SAR, and use the intrinsic CAR as a smoother." date: 2026-08-13 categories: [spatial, regression, autocorrelation, areal data] image: thumbnail.png --- The [simultaneous autoregressive models](../spatial-autoregressive-models-sar/) build spatial dependence from the top down: one global equation, `y = rho W y + X beta + eps`, solved all at once. Conditional autoregressive (CAR) models take the opposite route. They specify, for each cell, only its distribution *given its neighbours*, and let a theorem stitch those local pieces into a coherent joint distribution. This bottom-up construction, due to Besag (1974), is the engine behind most Bayesian disease mapping and behind the spatial random effects in ecological hierarchical models. It is worth understanding directly. This tutorial writes the conditional model down, assembles its precision matrix, fits a proper CAR by hand in base R, and then does three things the textbooks often skip: it shows why the CAR autocorrelation parameter is not comparable across software packages, it contrasts the covariance a CAR implies with the one a SAR implies, and it uses the intrinsic CAR as a practical smoother. ## The conditional specification A CAR model says that, conditional on every other cell, cell `i` is normal with a mean pulled towards its neighbours: $$ y_i \mid y_{-i} \sim \mathcal{N}\!\left(\mu_i + \rho \sum_j \frac{a_{ij}}{d_i}(y_j - \mu_j),\ \frac{\tau^2}{d_i}\right), $$ where `a_ij` marks neighbours, `d_i` is the number of neighbours of cell `i`, and `mu = X beta`. The Hammersley-Clifford theorem guarantees these conditionals correspond to a single joint Gaussian, and its precision matrix (the inverse covariance) is wonderfully sparse: $$ Q = \frac{1}{\tau^2}(D - \rho A), $$ with `D` the diagonal matrix of neighbour counts and `A` the adjacency matrix. The model is *proper* (a valid distribution) when `D - rho A` is positive definite, which pins `rho` inside a range set by the eigenvalues. ## Setting up the lattice ```{r setup, include=FALSE} library(ggplot2) library(grid) knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, dev = "png", dev.args = list(bg = "white"), dpi = 150, fig.width = 8, fig.height = 3.5) theme_te <- function() { theme_minimal(base_size = 11) + theme(plot.background = element_rect(fill = "white", colour = NA), panel.background = element_rect(fill = "white", colour = NA), panel.grid = element_blank(), plot.title = element_text(face = "bold", size = 11), plot.subtitle = element_text(size = 9, colour = "grey30"), axis.text = element_blank(), axis.ticks = element_blank(), legend.key.width = unit(0.35, "cm")) } tilemap <- function(df, value, title, subtitle = NULL) { ggplot(df, aes(col, row, fill = value)) + geom_tile() + coord_equal() + scale_fill_viridis_c(option = "mako", direction = -1) + labs(title = title, subtitle = subtitle, x = NULL, y = NULL, fill = NULL) + theme_te() } two_panel <- function(g1, g2) { grid.newpage(); pushViewport(viewport(layout = grid.layout(1, 2))) print(g1, vp = viewport(layout.pos.row = 1, layout.pos.col = 1)) print(g2, vp = viewport(layout.pos.row = 1, layout.pos.col = 2)) } ``` ```{r lattice} side <- 20; n <- side * side coords <- expand.grid(col = 1:side, row = 1:side) A <- matrix(0, n, n) idx <- function(r, c) (c - 1) * side + r for (r in 1:side) for (c in 1:side) { i <- idx(r, c) if (r > 1) A[i, idx(r - 1, c)] <- 1 if (r < side) A[i, idx(r + 1, c)] <- 1 if (c > 1) A[i, idx(r, c - 1)] <- 1 if (c < side) A[i, idx(r, c + 1)] <- 1 } deg <- rowSums(A); D <- diag(deg) # eigenvalues of the row-standardised operator fix the admissible rho range ev <- eigen(diag(1 / sqrt(deg)) %*% A %*% diag(1 / sqrt(deg)), symmetric = TRUE)$values rho_lo <- 1 / min(ev); rho_hi <- 1 / max(ev) c(lower = rho_lo, upper = rho_hi) ``` ## Simulating a proper CAR field We draw one realisation with an intercept, a covariate, a strong conditional dependence of `rho = 0.9`, and a conditional scale of `tau = 0.4`. Sampling is direct from the joint covariance `tau^2 (D - rho A)^{-1}`. ```{r simulate} set.seed(4231) cc <- scale(coords$col)[, 1] x <- as.numeric(scale(0.8 * cos(pi * coords$row / 7) + 0.6 * cc + rnorm(n, 0, 0.6))) X <- cbind(1, x) b0 <- 0.5; b1 <- 0.7; rho <- 0.9; tau <- 0.4 Q <- (D - rho * A) / tau^2 set.seed(11) y <- as.numeric(X %*% c(b0, b1) + t(chol(solve(Q))) %*% rnorm(n)) ``` ## Fitting a CAR by hand The log-likelihood needs the log-determinant of the precision. Factor it as `|D - rho A| = |D| * |I - rho D^{-1}A|`, so the determinant reduces to a sum over the same eigenvalues we already have. For each candidate `rho` we solve the generalised least squares problem for `beta`, read off `tau` from the quadratic form, and search over the single remaining parameter. ```{r car-fit} nll_car <- function(r) { if (r <= rho_lo || r >= rho_hi) return(1e10) P <- D - r * A # precision up to 1/tau^2 be <- solve(t(X) %*% P %*% X, t(X) %*% P %*% y) e <- y - X %*% be; t2 <- as.numeric(t(e) %*% P %*% e) / n ldetP <- sum(log(deg)) + sum(log(1 - r * ev)) # log|D - rho A| -(-n / 2 * log(2 * pi) - n / 2 - n / 2 * log(t2) + 0.5 * ldetP) } opt <- optimise(nll_car, c(rho_lo + 1e-4, rho_hi - 1e-4)) car_rho <- opt$minimum P <- D - car_rho * A be <- solve(t(X) %*% P %*% X, t(X) %*% P %*% y) car_b0 <- be[1]; car_b1 <- be[2] car_tau <- sqrt(as.numeric(t(y - X %*% be) %*% P %*% (y - X %*% be)) / n) car_ll <- -opt$objective round(c(rho = car_rho, tau = car_tau, b0 = car_b0, b1 = car_b1), 4) ``` The fit recovers the generating values well: a conditional dependence of `r round(car_rho, 2)` against the true `r rho`, a scale of `r round(car_tau, 2)` against `r tau`, and a slope of `r round(car_b1, 2)` against `r b1`. The likelihood machinery is the same log-determinant plus generalised least squares we used for the error model in the SAR tutorial, only the precision differs. ## A warning about conventions Reach for a package and you may be surprised. The `spautolm` function in spdep, with `family = "CAR"`, does not use the neighbour-count precision `D - rho A`. It uses a *symmetric-weights* convention with precision `I - lambda A`. That is a different model, so its `lambda` is not on the same scale as our `rho` and the two numbers will never agree. To show this is a convention and not an error, we can fit the symmetric version by hand as well. ```{r car-symmetric} evA <- eigen(A, symmetric = TRUE)$values loA <- 1 / min(evA); hiA <- 1 / max(evA) nll_sym <- function(lam) { if (lam <= loA || lam >= hiA) return(1e10) P <- diag(n) - lam * A be <- solve(t(X) %*% P %*% X, t(X) %*% P %*% y); e <- y - X %*% be s2 <- as.numeric(t(e) %*% P %*% e) / n -(-n / 2 * log(2 * pi) - n / 2 - n / 2 * log(s2) + 0.5 * sum(log(1 - lam * evA))) } car_sym_lambda <- optimise(nll_sym, c(loA + 1e-6, hiA - 1e-6))$minimum round(car_sym_lambda, 4) ``` The symmetric convention reports `r round(car_sym_lambda, 3)`, and this by-hand value matches spdep's `spautolm` to numerical precision. It describes the same data as our `r round(car_rho, 2)`, just on a different scale. The lesson is practical: never compare a CAR autocorrelation parameter across papers or packages without first checking which weights convention each one uses. ## CAR is not SAR in disguise Both models encode neighbourhood dependence, but they imply different covariances, and neither is a reparametrisation of the other except in special cases (Ver Hoef et al. 2018). A useful way to see the difference is to trace how the correlation between a central cell and the rest decays with graph distance, the number of steps through the lattice. ```{r reach} centre <- idx(10, 10) gdist <- rep(NA_integer_, n); gdist[centre] <- 0L; step <- 0L; front <- centre repeat { nb <- setdiff(unique(unlist(lapply(front, function(i) which(A[i, ] > 0)))), which(!is.na(gdist))) if (!length(nb)) break step <- step + 1L; gdist[nb] <- step; front <- nb } Ccar <- cov2cor(solve(D - 0.9 * A))[centre, ] # CAR reach Wrs <- A / deg Csar <- cov2cor(solve(crossprod(diag(n) - 0.6 * Wrs)))[centre, ] # SAR reach reach <- aggregate(data.frame(CAR = Ccar, SAR = Csar), list(dist = gdist), mean) round(head(reach, 5), 3) ``` At the first step the two are tuned to match, near `r round(reach$CAR[reach$dist == 1], 2)`. Further out they part company: at three steps the CAR retains a correlation of `r round(reach$CAR[reach$dist == 3], 3)` against the SAR's `r round(reach$SAR[reach$dist == 3], 3)`. The conditional model here carries dependence further across the map than the simultaneous one at matched short range. Which shape suits your system is an empirical question, and one more reason the two model families are worth keeping distinct. ```{r fig1-reach} #| echo: false #| fig-alt: "Two panels. The left panel is a line chart of correlation against graph distance for CAR and SAR, both starting at one and decaying, with the CAR line staying above the SAR line beyond the first step. The right panel is a grid map of a single simulated proper-CAR field showing smooth spatial patches." long <- data.frame(dist = rep(reach$dist, 2), corr = c(reach$CAR, reach$SAR), model = rep(c("CAR", "SAR"), each = nrow(reach))) long <- long[long$dist <= 8, ] g1 <- ggplot(long, aes(dist, corr, colour = model)) + geom_line(linewidth = 0.9) + geom_point(size = 1.6) + scale_colour_manual(values = c(CAR = "#275139", SAR = "#c98a1a")) + labs(title = "Correlation with the centre cell", subtitle = "by number of steps through the lattice", x = "graph distance", y = "correlation", colour = NULL) + theme_te() + theme(axis.text = element_text(size = 8), panel.grid.major.y = element_line(colour = "grey92"), legend.position = c(0.8, 0.8)) g2 <- tilemap(coords, y, "A proper-CAR field", "one realisation, rho = 0.9") two_panel(g1, g2) ``` ## The intrinsic CAR as a smoother Push `rho` all the way to its boundary and the precision becomes `D - A`, the graph Laplacian. This is the *intrinsic* CAR (ICAR): an improper distribution (it has no overall level, hence the usual sum-to-zero constraint) that penalises differences between neighbours. It is the most widely used spatial random effect in ecology and epidemiology, because as a smoother it borrows strength across the map. To see it work, we take a noisy signal, one observation per cell, and recover the underlying smooth field. ```{r icar} set.seed(4232) truef <- as.numeric(scale(1.4 * sin(pi * coords$col / 8) * cos(pi * coords$row / 9))) yn <- truef + rnorm(n, 0, 1.0) # noisy observation per cell R_icar <- D - A # intrinsic precision (graph Laplacian) # GMRF smoother S = (I + kappa R)^{-1}; choose kappa by generalised cross-validation gcv <- function(kappa) { S <- solve(diag(n) + kappa * R_icar); fit <- S %*% yn n * sum((yn - fit)^2) / (n - sum(diag(S)))^2 } ks <- 10^seq(-2, 2, length = 40); k_opt <- ks[which.min(sapply(ks, gcv))] Sopt <- solve(diag(n) + k_opt * R_icar) fit_ic <- as.numeric(Sopt %*% yn) icar_kappa <- k_opt icar_edf <- sum(diag(Sopt)) icar_rmse_raw <- sqrt(mean((yn - truef)^2)) icar_rmse_sm <- sqrt(mean((fit_ic - truef)^2)) round(c(kappa = icar_kappa, edf = icar_edf, rmse_raw = icar_rmse_raw, rmse_smooth = icar_rmse_sm), 3) ``` Cross-validation settles on a smoother with about `r round(icar_edf, 0)` effective degrees of freedom out of `r n` cells. It cuts the root mean squared error against the truth from `r round(icar_rmse_raw, 2)` for the raw observations to `r round(icar_rmse_sm, 2)`, roughly halving it, simply by letting each cell learn from its neighbours. This borrowing of strength is exactly what an ICAR random effect does inside a larger model. ```{r fig2-smooth} #| echo: false #| fig-alt: "Two grid maps side by side. The left panel shows noisy per-cell observations with a speckled appearance. The right panel shows the ICAR-smoothed field, which reveals a clean smooth spatial pattern with the noise removed." lim <- range(c(yn, fit_ic)) p1 <- tilemap(coords, yn, "Noisy observations", "one draw per cell") + scale_fill_viridis_c(option = "mako", direction = -1, limits = lim) p2 <- tilemap(coords, fit_ic, "ICAR smoother", paste0("edf = ", round(icar_edf, 0))) + scale_fill_viridis_c(option = "mako", direction = -1, limits = lim) two_panel(p1, p2) ``` ## What to carry forward A CAR model is built from local conditional rules and summarised by a sparse precision matrix `D - rho A`. It is the natural spatial term for hierarchical and generalised linear mixed models, where the conditional form slots neatly into the wider structure. Three points are worth remembering. The autocorrelation parameter depends on the weights convention, so it is not portable between software. A CAR and a SAR on the same neighbourhood imply genuinely different covariances, and the choice between them is substantive, not cosmetic. And the intrinsic CAR, the boundary case, is less a model of correlation than a smoother, penalising neighbour differences to pool information. That smoothing power comes with a hidden cost. When the spatial random effect is flexible enough to mimic a spatially structured covariate, it can quietly absorb that covariate's effect. This is spatial confounding, and [checking a spatial regression](../checking-a-spatial-regression/) takes it apart in detail. ## Related tutorials - [Spatial lag and spatial error models (SAR)](../spatial-autoregressive-models-sar/) build areal dependence from a single simultaneous equation. - [Moran eigenvector spatial filtering](../moran-eigenvector-spatial-filtering/) offers a third route that fits inside a generalised linear model. - [Checking a spatial regression](../checking-a-spatial-regression/) compares CAR, SAR and filtering and exposes spatial confounding. - [Generalised least squares for spatial correlation](../gls-spatial-correlation/) covers the continuous-distance analogue. ## References Besag, J. (1974). Spatial interaction and the statistical analysis of lattice systems. *Journal of the Royal Statistical Society: Series B*, 36(2), 192-236. https://doi.org/10.1111/j.2517-6161.1974.tb00999.x Dormann, C. F., McPherson, J. M., Araujo, M. B., Bivand, R., Bolliger, J., Carl, G., et al. (2007). Methods to account for spatial autocorrelation in the analysis of species distributional data: a review. *Ecography*, 30(5), 609-628. https://doi.org/10.1111/j.2007.0906-7590.05171.x Ver Hoef, J. M., Peterson, E. E., Hooten, M. B., Hanks, E. M., & Fortin, M.-J. (2018). Spatial autoregressive models for statistical inference from ecological data. *Ecological Monographs*, 88(1), 36-59. https://doi.org/10.1002/ecm.1283