Spatial lag and spatial error models (SAR)

spatial

regression

autocorrelation

areal data

Fit simultaneous autoregressive models to areal ecological data by hand, separate spillover from nuisance autocorrelation, and read the spatial multiplier.

Author

Tidy Ecology

Published

2026-08-13

Count a species across a grid of survey cells, or average a trait across administrative regions, and you have areal data: one value per polygon, with a neighbourhood structure linking touching units. Ordinary least squares treats those cells as independent. They rarely are. Neighbouring cells share climate, share dispersers, share history, and so their residuals travel together. We met this problem for point locations in the Moran’s I tutorial and handled continuous-distance correlation with generalised least squares. Here the geometry is a lattice, and the natural tools are the simultaneous autoregressive (SAR) models: the spatial lag model and the spatial error model.

The two models look similar and are often confused, but they say different things about the world. One describes a process that spills across boundaries; the other describes a nuisance that merely clusters. Getting the distinction right changes both your coefficients and their standard errors, so this tutorial builds both from the likelihood up, in base R, and shows what ordinary least squares gets wrong.

The lattice and its weights

Everything starts from a neighbourhood. We use a regular 20 by 20 grid and rook contiguity: two cells are neighbours if they share an edge. The adjacency matrix A holds a 1 for every neighbour pair. Dividing each row by its number of neighbours gives the row-standardised weights matrix W, so that W %*% y returns, for every cell, the average of its neighbours.

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)          # number of neighbours per cell
W   <- A / deg             # row-standardised weights

The likelihood needs the eigenvalues of W. Because W is similar to a symmetric matrix, its eigenvalues are real and fall in the interval from minus one to one, which also fixes the admissible range of the autocorrelation parameter.

ev <- eigen(diag(1 / sqrt(deg)) %*% A %*% diag(1 / sqrt(deg)), symmetric = TRUE)$values
rho_lo <- 1 / min(ev); rho_hi <- 1 / max(ev)
range(ev)

[1] -1  1

A reusable Moran’s I

We lean on Moran’s I throughout, so here is a compact base R version returning the statistic and its analytic p-value under the randomisation assumption. It matches the value from the spdep package exactly.

moran_test <- function(z, W) {
  n <- length(z); z <- z - mean(z)
  S0 <- sum(W); S1 <- 0.5 * sum((W + t(W))^2)
  S2 <- sum((rowSums(W) + colSums(W))^2)
  I  <- (n / S0) * as.numeric(t(z) %*% W %*% z) / sum(z^2)
  EI <- -1 / (n - 1); b2 <- n * sum(z^4) / (sum(z^2)^2)
  VI <- (n * ((n^2 - 3 * n + 3) * S1 - n * S2 + 3 * S0^2)
         - b2 * ((n^2 - n) * S1 - 2 * n * S2 + 6 * S0^2)) /
        ((n - 1) * (n - 2) * (n - 3) * S0^2) - EI^2
  list(I = I, p = pnorm((I - EI) / sqrt(VI), lower.tail = FALSE))
}

Two models, two stories

Write X for the design matrix holding an intercept and one covariate, and beta for the coefficients. The spatial lag model places the neighbour average of the response on the right-hand side:

\[ y = \rho W y + X\beta + \varepsilon. \]

Solving for y gives the reduced form y = (I - rho W)^{-1}(X beta + eps). The inverse (I - rho W)^{-1} is the spatial multiplier: a nudge to one cell propagates to its neighbours, then their neighbours, and so on. This is a model of a genuine process, dispersal, contagion, competition, that crosses cell boundaries.

The spatial error model puts the autocorrelation in the disturbance instead:

\[ y = X\beta + u, \qquad u = \lambda W u + \varepsilon. \]

Here the neighbourhood structure is a nuisance. The mean model X beta is what you care about; the clustered errors only spoil the standard errors. Both models are estimated by maximum likelihood, and both reduce to a one-dimensional search once you concentrate out the regression coefficients.

Simulating a spillover process

We generate data from the lag process with a modest spatial covariate, an autocorrelation of rho = 0.6, and a real slope of b1 = 0.8.

set.seed(4230)
cc <- scale(coords$col)[, 1]; rr <- scale(coords$row)[, 1]
x  <- as.numeric(scale(0.9 * sin(pi * coords$col / 6) + 0.7 * cc - 0.5 * rr +
                       rnorm(n, 0, 0.6)))
X  <- cbind(1, x)
b0 <- 1.0; b1 <- 0.8; rho <- 0.6; sig <- 0.5
Ainv <- solve(diag(n) - rho * W)
y <- as.numeric(Ainv %*% (X %*% c(b0, b1) + rnorm(n, 0, sig)))

Where ordinary least squares fails

Fit the naive regression and look at its residuals.

ols <- lm(y ~ x)
sar_ols_int   <- coef(ols)[1]
sar_ols_slope <- coef(ols)["x"]
sar_mi_ols <- moran_test(resid(ols), W)$I
sar_mp_ols <- moran_test(resid(ols), W)$p
round(c(intercept = sar_ols_int, slope = sar_ols_slope), 3)

intercept.(Intercept)               slope.x 
                2.502                 1.671

The true slope is 0.8, but ordinary least squares returns 1.67: more than double, because the omitted term rho W y is correlated with x through the multiplier and its pull is absorbed into the slope. The residuals carry a Moran’s I of 0.245 (p around 4.9^{-12}), so the independence assumption is plainly violated. Neither the coefficient nor its standard error can be trusted.

The spatial lag model by hand

The concentrated log-likelihood for the lag model, for a candidate rho, filters the response with A(rho) = I - rho W, runs a plain regression of the filtered response on X, and adds the log-determinant log|I - rho W|, which the eigenvalues make cheap.

nll_lag <- function(par) {
  be <- par[1:2]; r <- par[3]; s2 <- exp(2 * par[4])
  if (r <= rho_lo || r >= rho_hi) return(1e10)
  Ay <- y - r * (W %*% y)
  e  <- Ay - X %*% be
  -(-n / 2 * log(2 * pi) - n / 2 * log(s2) +
      sum(log(1 - r * ev)) - sum(e^2) / (2 * s2))
}
fit_lag <- optim(c(coef(ols), 0.1, log(sd(resid(ols)))), nll_lag,
                 method = "BFGS", hessian = TRUE)
se_lag <- sqrt(diag(solve(fit_lag$hessian)))
sar_slm_b0    <- fit_lag$par[1]; sar_slm_b1  <- fit_lag$par[2]
sar_slm_rho   <- fit_lag$par[3]; sar_slm_serho <- se_lag[3]
sar_slm_sigma <- exp(fit_lag$par[4])
sar_slm_ll    <- -fit_lag$value
sar_slm_aic   <- 2 * 4 + 2 * fit_lag$value
sar_mi_slm    <- moran_test(as.numeric((diag(n) - sar_slm_rho * W) %*% y -
                                       X %*% fit_lag$par[1:2]), W)$I
data.frame(estimate = round(c(b0 = sar_slm_b0, b1 = sar_slm_b1,
                              rho = sar_slm_rho, sigma = sar_slm_sigma), 4))

               estimate
b0.(Intercept)   1.0141
b1.x             0.8168
rho              0.5946
sigma            0.4949

The estimates land on the truth: rho of 0.595 against 0.6, a slope of 0.817 against 0.8, and an intercept of 1.01 against 1. These match the spdep package’s lagsarlm to within numerical precision. The filtered residuals now carry a Moran’s I of only 0.011, so the spatial structure has been accounted for rather than left in the errors.

Two grid maps side by side. The left panel shows the covariate as a smooth spatial gradient; the right panel shows the response, which follows the covariate but with the added smoothing of the spatial multiplier.

Reading the spatial multiplier

In a plain regression, a one-unit change in x shifts the response by the slope, full stop. In a lag model the effect ripples outward, so the total effect is larger than the coefficient. The impact matrix (I - rho W)^{-1} b1 splits into a direct effect (the average change in a cell from its own covariate, including feedback loops), an indirect effect (the spillover onto other cells), and their sum.

M <- solve(diag(n) - sar_slm_rho * W) * sar_slm_b1
sar_direct   <- mean(diag(M))
sar_total    <- mean(rowSums(M))
sar_indirect <- sar_total - sar_direct
round(c(direct = sar_direct, indirect = sar_indirect, total = sar_total), 3)

  direct indirect    total 
   0.915    1.100    2.015

The direct effect, 0.92, sits just above the coefficient because of feedback through the neighbourhood. The indirect effect, 1.1, is the spillover, and it is larger than the direct part here. The total, 2.01, equals the tidy closed form b1 / (1 - rho) exactly, since each row of the multiplier sums to 1 / (1 - rho) for row-standardised weights. Reporting only the coefficient would understate the ecological effect by more than half. This decomposition is the single most important reason to prefer a lag model when the process really does cross boundaries.

The spatial error model, and why it is the wrong fit here

Fit the error model to the same data. Its concentrated likelihood filters both sides by B(lambda) = I - lambda W and runs generalised least squares.

nll_err <- function(par) {
  lam <- par[1]; s2 <- exp(2 * par[2])
  if (lam <= rho_lo || lam >= rho_hi) return(1e10)
  B <- diag(n) - lam * W; BtB <- crossprod(B)
  be <- solve(t(X) %*% BtB %*% X, t(X) %*% BtB %*% y)
  e  <- B %*% (y - X %*% be)
  -(-n / 2 * log(2 * pi) - n / 2 * log(s2) +
      sum(log(1 - lam * ev)) - sum(e^2) / (2 * s2))
}
fit_err <- optim(c(0.1, log(sd(resid(ols)))), nll_err, method = "BFGS")
lam_hat <- fit_err$par[1]
Bh <- diag(n) - lam_hat * W; BtB <- crossprod(Bh)
be_err <- solve(t(X) %*% BtB %*% X, t(X) %*% BtB %*% y)
sar_sem_lambda <- lam_hat; sar_sem_b1 <- be_err[2]
sar_sem_aic <- 2 * 4 + 2 * fit_err$value
round(c(lambda = sar_sem_lambda, b1 = sar_sem_b1, AIC = sar_sem_aic), 3)

 lambda      b1     AIC 
  0.894   0.748 722.758

The error model reaches for a large lambda of 0.89 to mop up the leftover structure, but it recovers the wrong slope, 0.75, and its AIC of 723 is far worse than the lag model’s 623. That gap is the model telling you the autocorrelation here is a process, not a nuisance. On data actually generated by an error process the verdict would flip. Deciding between the two from data, rather than assuming, is the job of the formal tests in checking a spatial regression.

Two residual grid maps side by side. The left panel shows ordinary least squares residuals with visible patches of like colour, indicating spatial clustering. The right panel shows spatial lag model residuals that look like random noise with no patches.

How much does it matter? A short simulation

One dataset can mislead. We repeat the whole exercise 400 times, comparing the ordinary least squares slope with the lag estimate, and asking how often each ninety-five per cent interval actually covers the true slope.

set.seed(99); B <- 400
ols_b <- slm_b <- ols_cov <- slm_cov <- numeric(B)
for (b in 1:B) {
  yb <- as.numeric(Ainv %*% (X %*% c(b0, b1) + rnorm(n, 0, sig)))
  o  <- lm(yb ~ x); ob <- coef(o)["x"]
  ose <- summary(o)$coefficients["x", "Std. Error"]
  ols_b[b] <- ob; ols_cov[b] <- (b1 >= ob - 1.96 * ose & b1 <= ob + 1.96 * ose)
  fl <- optim(c(coef(o), 0.1, log(sd(resid(o)))), function(par) {
    be <- par[1:2]; r <- par[3]; s2 <- exp(2 * par[4])
    if (r <= rho_lo || r >= rho_hi) return(1e10)
    Ay <- yb - r * (W %*% yb); e <- Ay - X %*% be
    -(-n / 2 * log(2 * pi) - n / 2 * log(s2) + sum(log(1 - r * ev)) - sum(e^2) / (2 * s2))
  }, method = "BFGS", hessian = TRUE)
  sb <- fl$par[2]; sse <- sqrt(diag(solve(fl$hessian)))[2]
  slm_b[b] <- sb; slm_cov[b] <- (b1 >= sb - 1.96 * sse & b1 <= sb + 1.96 * sse)
}
sar_mc_ols_bias <- 100 * (mean(ols_b) - b1) / b1
sar_mc_ols_cov  <- mean(ols_cov)
sar_mc_slm_bias <- 100 * (mean(slm_b) - b1) / b1
sar_mc_slm_cov  <- mean(slm_cov)
round(c(ols_bias_pct = sar_mc_ols_bias, ols_cover = sar_mc_ols_cov,
        slm_bias_pct = sar_mc_slm_bias, slm_cover = sar_mc_slm_cov), 3)

ols_bias_pct    ols_cover slm_bias_pct    slm_cover 
     108.085        0.000        0.716        0.940

Ordinary least squares overstates the slope by 108 per cent on average, and its confidence interval covers the truth in 0 per cent of runs: a total failure of inference. The lag model is essentially unbiased, off by 0.7 per cent, with 94 per cent coverage, close to the nominal ninety-five.

What to carry forward

For areal data on a lattice, the choice is not whether to model space but how. A spatial lag model treats neighbourhood dependence as a process and rewards you with an interpretable multiplier that separates direct effects from spillover. A spatial error model treats it as a nuisance and simply repairs the standard errors. Fitting the wrong one, or fitting neither, distorts both estimates and inference, as the simulation makes plain.

Two honest cautions. Every result here assumes the neighbourhood matrix W is correct; a different contiguity rule gives different numbers, and the choice is rarely obvious. And these models assume the spatial dependence is stationary across the map. The conditional autoregressive models tutorial offers a different route to the same goal, built on local conditional distributions rather than a simultaneous system, and checking a spatial regression shows how to test which specification your data actually support and where even a good spatial model can quietly mislead you.

References

Anselin, L. (1988). Spatial Econometrics: Methods and Models. Kluwer. https://doi.org/10.1007/978-94-015-7799-1

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

Kissling, W. D., & Carl, G. (2008). Spatial autocorrelation and the selection of simultaneous autoregressive models. Global Ecology and Biogeography, 17(1), 59-71. https://doi.org/10.1111/j.1466-8238.2007.00334.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: "Spatial lag and spatial error models (SAR)" description: "Fit simultaneous autoregressive models to areal ecological data by hand, separate spillover from nuisance autocorrelation, and read the spatial multiplier." date: 2026-08-13 categories: [spatial, regression, autocorrelation, areal data] image: thumbnail.png --- Count a species across a grid of survey cells, or average a trait across administrative regions, and you have *areal* data: one value per polygon, with a neighbourhood structure linking touching units. Ordinary least squares treats those cells as independent. They rarely are. Neighbouring cells share climate, share dispersers, share history, and so their residuals travel together. We met this problem for point locations in the [Moran's I tutorial](../spatial-autocorrelation-morans-i/) and handled continuous-distance correlation with [generalised least squares](../gls-spatial-correlation/). Here the geometry is a lattice, and the natural tools are the *simultaneous autoregressive* (SAR) models: the spatial lag model and the spatial error model. The two models look similar and are often confused, but they say different things about the world. One describes a process that spills across boundaries; the other describes a nuisance that merely clusters. Getting the distinction right changes both your coefficients and their standard errors, so this tutorial builds both from the likelihood up, in base R, and shows what ordinary least squares gets wrong. ## The lattice and its weights Everything starts from a neighbourhood. We use a regular 20 by 20 grid and *rook* contiguity: two cells are neighbours if they share an edge. The adjacency matrix `A` holds a 1 for every neighbour pair. Dividing each row by its number of neighbours gives the *row-standardised* weights matrix `W`, so that `W %*% y` returns, for every cell, the average of its neighbours. ```{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, diverging = FALSE) { g <- ggplot(df, aes(col, row, fill = value)) + geom_tile() + coord_equal() + labs(title = title, subtitle = subtitle, x = NULL, y = NULL, fill = NULL) + theme_te() if (diverging) g + scale_fill_gradient2(low = "#3b4cc0", mid = "white", high = "#b40426", midpoint = 0) else g + scale_fill_viridis_c(option = "mako", direction = -1) } 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) # number of neighbours per cell W <- A / deg # row-standardised weights ``` The likelihood needs the eigenvalues of `W`. Because `W` is similar to a symmetric matrix, its eigenvalues are real and fall in the interval from minus one to one, which also fixes the admissible range of the autocorrelation parameter. ```{r eigen} ev <- eigen(diag(1 / sqrt(deg)) %*% A %*% diag(1 / sqrt(deg)), symmetric = TRUE)$values rho_lo <- 1 / min(ev); rho_hi <- 1 / max(ev) range(ev) ``` ## A reusable Moran's I We lean on Moran's I throughout, so here is a compact base R version returning the statistic and its analytic p-value under the randomisation assumption. It matches the value from the spdep package exactly. ```{r moran} moran_test <- function(z, W) { n <- length(z); z <- z - mean(z) S0 <- sum(W); S1 <- 0.5 * sum((W + t(W))^2) S2 <- sum((rowSums(W) + colSums(W))^2) I <- (n / S0) * as.numeric(t(z) %*% W %*% z) / sum(z^2) EI <- -1 / (n - 1); b2 <- n * sum(z^4) / (sum(z^2)^2) VI <- (n * ((n^2 - 3 * n + 3) * S1 - n * S2 + 3 * S0^2) - b2 * ((n^2 - n) * S1 - 2 * n * S2 + 6 * S0^2)) / ((n - 1) * (n - 2) * (n - 3) * S0^2) - EI^2 list(I = I, p = pnorm((I - EI) / sqrt(VI), lower.tail = FALSE)) } ``` ## Two models, two stories Write `X` for the design matrix holding an intercept and one covariate, and `beta` for the coefficients. The **spatial lag model** places the neighbour average of the *response* on the right-hand side: $$ y = \rho W y + X\beta + \varepsilon. $$ Solving for `y` gives the reduced form `y = (I - rho W)^{-1}(X beta + eps)`. The inverse `(I - rho W)^{-1}` is the *spatial multiplier*: a nudge to one cell propagates to its neighbours, then their neighbours, and so on. This is a model of a genuine process, dispersal, contagion, competition, that crosses cell boundaries. The **spatial error model** puts the autocorrelation in the disturbance instead: $$ y = X\beta + u, \qquad u = \lambda W u + \varepsilon. $$ Here the neighbourhood structure is a nuisance. The mean model `X beta` is what you care about; the clustered errors only spoil the standard errors. Both models are estimated by maximum likelihood, and both reduce to a one-dimensional search once you concentrate out the regression coefficients. ## Simulating a spillover process We generate data from the lag process with a modest spatial covariate, an autocorrelation of `rho = 0.6`, and a real slope of `b1 = 0.8`. ```{r simulate} set.seed(4230) cc <- scale(coords$col)[, 1]; rr <- scale(coords$row)[, 1] x <- as.numeric(scale(0.9 * sin(pi * coords$col / 6) + 0.7 * cc - 0.5 * rr + rnorm(n, 0, 0.6))) X <- cbind(1, x) b0 <- 1.0; b1 <- 0.8; rho <- 0.6; sig <- 0.5 Ainv <- solve(diag(n) - rho * W) y <- as.numeric(Ainv %*% (X %*% c(b0, b1) + rnorm(n, 0, sig))) ``` ## Where ordinary least squares fails Fit the naive regression and look at its residuals. ```{r ols} ols <- lm(y ~ x) sar_ols_int <- coef(ols)[1] sar_ols_slope <- coef(ols)["x"] sar_mi_ols <- moran_test(resid(ols), W)$I sar_mp_ols <- moran_test(resid(ols), W)$p round(c(intercept = sar_ols_int, slope = sar_ols_slope), 3) ``` The true slope is `r b1`, but ordinary least squares returns `r round(sar_ols_slope, 2)`: more than double, because the omitted term `rho W y` is correlated with `x` through the multiplier and its pull is absorbed into the slope. The residuals carry a Moran's I of `r round(sar_mi_ols, 3)` (p around `r signif(sar_mp_ols, 2)`), so the independence assumption is plainly violated. Neither the coefficient nor its standard error can be trusted. ## The spatial lag model by hand The concentrated log-likelihood for the lag model, for a candidate `rho`, filters the response with `A(rho) = I - rho W`, runs a plain regression of the filtered response on `X`, and adds the log-determinant `log|I - rho W|`, which the eigenvalues make cheap. ```{r slm} nll_lag <- function(par) { be <- par[1:2]; r <- par[3]; s2 <- exp(2 * par[4]) if (r <= rho_lo || r >= rho_hi) return(1e10) Ay <- y - r * (W %*% y) e <- Ay - X %*% be -(-n / 2 * log(2 * pi) - n / 2 * log(s2) + sum(log(1 - r * ev)) - sum(e^2) / (2 * s2)) } fit_lag <- optim(c(coef(ols), 0.1, log(sd(resid(ols)))), nll_lag, method = "BFGS", hessian = TRUE) se_lag <- sqrt(diag(solve(fit_lag$hessian))) sar_slm_b0 <- fit_lag$par[1]; sar_slm_b1 <- fit_lag$par[2] sar_slm_rho <- fit_lag$par[3]; sar_slm_serho <- se_lag[3] sar_slm_sigma <- exp(fit_lag$par[4]) sar_slm_ll <- -fit_lag$value sar_slm_aic <- 2 * 4 + 2 * fit_lag$value sar_mi_slm <- moran_test(as.numeric((diag(n) - sar_slm_rho * W) %*% y - X %*% fit_lag$par[1:2]), W)$I data.frame(estimate = round(c(b0 = sar_slm_b0, b1 = sar_slm_b1, rho = sar_slm_rho, sigma = sar_slm_sigma), 4)) ``` The estimates land on the truth: `rho` of `r round(sar_slm_rho, 3)` against `r rho`, a slope of `r round(sar_slm_b1, 3)` against `r b1`, and an intercept of `r round(sar_slm_b0, 2)` against `r b0`. These match the spdep package's `lagsarlm` to within numerical precision. The filtered residuals now carry a Moran's I of only `r round(sar_mi_slm, 3)`, so the spatial structure has been accounted for rather than left in the errors. ```{r fig1-maps} #| echo: false #| fig-alt: "Two grid maps side by side. The left panel shows the covariate as a smooth spatial gradient; the right panel shows the response, which follows the covariate but with the added smoothing of the spatial multiplier." df <- data.frame(coords, x = x, y = y) two_panel( tilemap(df, df$x, "Covariate x", "smooth spatial gradient"), tilemap(df, df$y, "Response y", "lag process, rho = 0.6") ) ``` ## Reading the spatial multiplier In a plain regression, a one-unit change in `x` shifts the response by the slope, full stop. In a lag model the effect ripples outward, so the total effect is larger than the coefficient. The impact matrix `(I - rho W)^{-1} b1` splits into a *direct* effect (the average change in a cell from its own covariate, including feedback loops), an *indirect* effect (the spillover onto other cells), and their sum. ```{r impacts} M <- solve(diag(n) - sar_slm_rho * W) * sar_slm_b1 sar_direct <- mean(diag(M)) sar_total <- mean(rowSums(M)) sar_indirect <- sar_total - sar_direct round(c(direct = sar_direct, indirect = sar_indirect, total = sar_total), 3) ``` The direct effect, `r round(sar_direct, 2)`, sits just above the coefficient because of feedback through the neighbourhood. The indirect effect, `r round(sar_indirect, 2)`, is the spillover, and it is larger than the direct part here. The total, `r round(sar_total, 2)`, equals the tidy closed form `b1 / (1 - rho)` exactly, since each row of the multiplier sums to `1 / (1 - rho)` for row-standardised weights. Reporting only the coefficient would understate the ecological effect by more than half. This decomposition is the single most important reason to prefer a lag model when the process really does cross boundaries. ## The spatial error model, and why it is the wrong fit here Fit the error model to the same data. Its concentrated likelihood filters *both* sides by `B(lambda) = I - lambda W` and runs generalised least squares. ```{r sem} nll_err <- function(par) { lam <- par[1]; s2 <- exp(2 * par[2]) if (lam <= rho_lo || lam >= rho_hi) return(1e10) B <- diag(n) - lam * W; BtB <- crossprod(B) be <- solve(t(X) %*% BtB %*% X, t(X) %*% BtB %*% y) e <- B %*% (y - X %*% be) -(-n / 2 * log(2 * pi) - n / 2 * log(s2) + sum(log(1 - lam * ev)) - sum(e^2) / (2 * s2)) } fit_err <- optim(c(0.1, log(sd(resid(ols)))), nll_err, method = "BFGS") lam_hat <- fit_err$par[1] Bh <- diag(n) - lam_hat * W; BtB <- crossprod(Bh) be_err <- solve(t(X) %*% BtB %*% X, t(X) %*% BtB %*% y) sar_sem_lambda <- lam_hat; sar_sem_b1 <- be_err[2] sar_sem_aic <- 2 * 4 + 2 * fit_err$value round(c(lambda = sar_sem_lambda, b1 = sar_sem_b1, AIC = sar_sem_aic), 3) ``` The error model reaches for a large `lambda` of `r round(sar_sem_lambda, 2)` to mop up the leftover structure, but it recovers the wrong slope, `r round(sar_sem_b1, 2)`, and its AIC of `r round(sar_sem_aic, 0)` is far worse than the lag model's `r round(sar_slm_aic, 0)`. That gap is the model telling you the autocorrelation here is a *process*, not a nuisance. On data actually generated by an error process the verdict would flip. Deciding between the two from data, rather than assuming, is the job of the formal tests in [checking a spatial regression](../checking-a-spatial-regression/). ```{r fig2-resid} #| echo: false #| fig-alt: "Two residual grid maps side by side. The left panel shows ordinary least squares residuals with visible patches of like colour, indicating spatial clustering. The right panel shows spatial lag model residuals that look like random noise with no patches." res_ols <- resid(ols) res_slm <- as.numeric((diag(n) - sar_slm_rho * W) %*% y - X %*% fit_lag$par[1:2]) dfr <- data.frame(coords, ols = res_ols, slm = res_slm) lim <- max(abs(c(res_ols, res_slm))) p1 <- tilemap(dfr, dfr$ols, "OLS residuals", paste0("Moran's I = ", round(sar_mi_ols, 3)), diverging = TRUE) + scale_fill_gradient2(low = "#3b4cc0", mid = "white", high = "#b40426", midpoint = 0, limits = c(-lim, lim)) p2 <- tilemap(dfr, dfr$slm, "Lag model residuals", paste0("Moran's I = ", round(sar_mi_slm, 3)), diverging = TRUE) + scale_fill_gradient2(low = "#3b4cc0", mid = "white", high = "#b40426", midpoint = 0, limits = c(-lim, lim)) two_panel(p1, p2) ``` ## How much does it matter? A short simulation One dataset can mislead. We repeat the whole exercise 400 times, comparing the ordinary least squares slope with the lag estimate, and asking how often each ninety-five per cent interval actually covers the true slope. ```{r mc} set.seed(99); B <- 400 ols_b <- slm_b <- ols_cov <- slm_cov <- numeric(B) for (b in 1:B) { yb <- as.numeric(Ainv %*% (X %*% c(b0, b1) + rnorm(n, 0, sig))) o <- lm(yb ~ x); ob <- coef(o)["x"] ose <- summary(o)$coefficients["x", "Std. Error"] ols_b[b] <- ob; ols_cov[b] <- (b1 >= ob - 1.96 * ose & b1 <= ob + 1.96 * ose) fl <- optim(c(coef(o), 0.1, log(sd(resid(o)))), function(par) { be <- par[1:2]; r <- par[3]; s2 <- exp(2 * par[4]) if (r <= rho_lo || r >= rho_hi) return(1e10) Ay <- yb - r * (W %*% yb); e <- Ay - X %*% be -(-n / 2 * log(2 * pi) - n / 2 * log(s2) + sum(log(1 - r * ev)) - sum(e^2) / (2 * s2)) }, method = "BFGS", hessian = TRUE) sb <- fl$par[2]; sse <- sqrt(diag(solve(fl$hessian)))[2] slm_b[b] <- sb; slm_cov[b] <- (b1 >= sb - 1.96 * sse & b1 <= sb + 1.96 * sse) } sar_mc_ols_bias <- 100 * (mean(ols_b) - b1) / b1 sar_mc_ols_cov <- mean(ols_cov) sar_mc_slm_bias <- 100 * (mean(slm_b) - b1) / b1 sar_mc_slm_cov <- mean(slm_cov) round(c(ols_bias_pct = sar_mc_ols_bias, ols_cover = sar_mc_ols_cov, slm_bias_pct = sar_mc_slm_bias, slm_cover = sar_mc_slm_cov), 3) ``` Ordinary least squares overstates the slope by `r round(sar_mc_ols_bias, 0)` per cent on average, and its confidence interval covers the truth in `r round(100 * sar_mc_ols_cov, 0)` per cent of runs: a total failure of inference. The lag model is essentially unbiased, off by `r round(sar_mc_slm_bias, 1)` per cent, with `r round(100 * sar_mc_slm_cov, 0)` per cent coverage, close to the nominal ninety-five. ## What to carry forward For areal data on a lattice, the choice is not whether to model space but how. A spatial lag model treats neighbourhood dependence as a process and rewards you with an interpretable multiplier that separates direct effects from spillover. A spatial error model treats it as a nuisance and simply repairs the standard errors. Fitting the wrong one, or fitting neither, distorts both estimates and inference, as the simulation makes plain. Two honest cautions. Every result here assumes the neighbourhood matrix `W` is correct; a different contiguity rule gives different numbers, and the choice is rarely obvious. And these models assume the spatial dependence is stationary across the map. The [conditional autoregressive models](../conditional-autoregressive-car-models/) tutorial offers a different route to the same goal, built on local conditional distributions rather than a simultaneous system, and [checking a spatial regression](../checking-a-spatial-regression/) shows how to test which specification your data actually support and where even a good spatial model can quietly mislead you. ## Related tutorials - [Spatial autocorrelation and Moran's I](../spatial-autocorrelation-morans-i/) introduces the neighbourhood weights and the statistic reused here. - [Conditional autoregressive (CAR) models](../conditional-autoregressive-car-models/) build areal dependence from local conditional distributions. - [Moran eigenvector spatial filtering](../moran-eigenvector-spatial-filtering/) accounts for autocorrelation inside a generalised linear model. - [Checking a spatial regression](../checking-a-spatial-regression/) tests lag against error and exposes spatial confounding. ## References Anselin, L. (1988). *Spatial Econometrics: Methods and Models*. Kluwer. https://doi.org/10.1007/978-94-015-7799-1 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 Kissling, W. D., & Carl, G. (2008). Spatial autocorrelation and the selection of simultaneous autoregressive models. *Global Ecology and Biogeography*, 17(1), 59-71. https://doi.org/10.1111/j.1466-8238.2007.00334.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