Spatial Regression: Lag and Error Models

Big Idea

The GLS chapter gave us one tool for handling spatially autocorrelated residuals: specify the correlation structure with a variogram and let GLS down-weight correlated observations. That works well. But there is a whole family of spatial regression models that take a different approach, and understanding them will round out your toolkit and sharpen your thinking about what spatial structure in a regression actually means.

GLS models the error correlation structure using a variogram – a continuous function of distance. The SAR (Simultaneous AutoRegressive) family instead uses an explicit spatial weights matrix: the same neighbor structure you built for Moran’s I and LISA. Two SAR models matter most for applied environmental work. The spatial error model says that the errors – the stuff your predictors don’t explain – are spatially correlated via a neighbor process. The spatial lag model says something more provocative: that the response variable at each location is directly influenced by the response at neighboring locations. These two models answer different questions, and it matters which one you fit.

Note

This module follows directly from the Regression: GLS chapter and assumes you are comfortable with the OLS-to-GLS workflow there. The spatial weights machinery is the same as in the LISA chapter.

Packages

Code 
library(spatialreg)
library(spdep)
library(sf)
library(tidyverse)

spatialreg (Bivand et al. 2021) carries the main functions for fitting SAR models. It was spun off from spdep (Bivand 2026) a few years ago, which means older code you find online might load these functions from spdep directly. If you see spdep::lagsarlm() in old tutorials, that is the same function. Load spatialreg explicitly. We also use sf (Pebesma 2026) for spatial data structures and tidyverse (Wickham 2023) for data wrangling.

Two Models, Two Questions

Let’s be concrete about what these models actually say.

A note on notation before we dive in. In what follows, \(y\) is a vector of \(n\) response values (one per location), \(X\) is the design matrix, \(\beta\) is the vector of regression coefficients, and \(W\) is the spatial weights matrix you already know from Moran’s I and LISA. The operation \(Wy\) computes a weighted average of \(y\) at neighboring locations – that is the spatial lag.

You will also notice that the same Greek letters appear here with different meanings than they had in earlier chapters. In the kriging chapter, \(\lambda\) was the vector of interpolation weights. Here it is a scalar measuring how strongly errors are correlated across neighbors. In the GLS chapter, \(\rho\) appeared in the code as the parameter we used to generate autocorrelated errors. Here \(\rho\) is the spatial autoregression coefficient on the response itself. This is not a conspiracy. Notation is there to serve the math, not to rule it, and the same letters get recycled constantly across fields, papers, and packages. Treat each symbol as local to what you are reading.

Spatial error model (SEM):

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

The residuals \(u\) follow a spatial autoregressive process. A residual at location \(i\) is a function of the residuals at its neighbors (via \(\lambda W\)) plus independent noise \(\varepsilon\). The coefficient \(\lambda\) (lambda) captures how strongly residuals are autocorrelated across neighbors. The coefficients \(\beta\) have the same interpretation as in OLS: a one-unit increase in \(x\) is associated with a \(\beta_1\) change in \(y\), on average, holding everything else constant.

Ecologically, the SEM is the right choice when space is a nuisance. Your predictors don’t capture all the spatial structure in \(y\), so the leftovers are autocorrelated. Classic example: you’re modeling tree growth as a function of soil nitrogen, but soil moisture, microclimate, and a dozen other unmeasured things also drive growth and are themselves spatially structured. The SEM absorbs that shared unmeasured spatial structure into the error term and gives you clean, valid inference on the nitrogen-growth relationship.

Spatial lag model (SLM):

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

Now the response variable \(y\) at each location appears on both sides of the equation. It is a function of \(Wy\) – the weighted average of \(y\) at neighboring locations – plus your predictors \(X\) and independent noise \(\varepsilon\). The coefficient \(\rho\) (rho) is a scalar that captures how strongly \(y\) at one location is pulled toward \(y\) at its neighbors.

This is a fundamentally different claim from the SEM. The SLM says there is a spatial spillover process: \(y\) at one location genuinely influences \(y\) at neighboring locations. Classic examples are disease spread (infection rate in one area influences neighboring areas), invasive species dispersal (abundance at one site drives colonization of nearby sites), or any system where diffusion, movement, or social contagion is part of the process. If the biology or ecology predicts that kind of mechanism, the lag model is the one to reach for.

One important caution: the \(\beta\) coefficients from the SLM are not total marginal effects. More on this below.

Choosing Between Them: LM Tests

There is a diagnostic test for choosing the appropriate model. The Rao score tests in spdep::lm.RStests() take an OLS fit and a weights matrix and test for evidence of spatial error vs spatial lag structure in the residuals. You get five test statistics from test = "all":

  • RSerr: tests for spatial error autocorrelation. Significant means the error model looks useful.
  • RSlag: tests for spatial lag structure. Significant means the lag model looks useful.
  • adjRSerr: robust test for error, controlling for the possibility of a lag. Use this when both RSerr and RSlag are significant.
  • adjRSlag: robust test for lag, controlling for error structure. Same situation.
  • SARMA: joint test for both.

The decision tree is simple. If only one of RSerr or RSlag is significant, fit that model. If both are significant, look at the robust tests and fit whichever of adjRSerr or adjRSlag is more significant. If neither is significant, OLS is fine.

Toy Example

We will build simulated data following the same spirit as the GLS chapter: a response \(y\) that is a function of \(x\) with spatially autocorrelated errors. We know going in that the spatial error model is the right answer. We will run the LM tests, confirm they point to the error model, and then fit both models to see what each gives us.

Data

Code 
n        <- 75
easting  <- runif(n, 0, 100)
northing <- runif(n, 0, 100)
points   <- cbind(easting, northing)

# Build spatially autocorrelated errors via a SAR process
dnb   <- dnearneigh(x = points, d1 = 0, d2 = 150)
dsts  <- nbdists(nb = dnb, coords = points)
p     <- 2.25
idw   <- lapply(dsts, function(x) x^-p)
wList <- nb2listw(neighbours = dnb, glist = idw, style = "W")

inv     <- spatialreg::invIrW(x = wList, rho = 0.75)
epsilon <- inv %*% rnorm(n)
epsilon <- scale(epsilon[,1])[,1]

# Response variable: y is a function of x plus autocorrelated noise
x  <- rnorm(n)
B0 <- 3
B1 <- 0.5
y  <- B0 + B1*x + epsilon

dat <- data.frame(easting, northing, y, x)

We also need a k-nearest-neighbor weights matrix for the SAR models. Same approach as the LISA chapter.

Code 
nb_k8 <- knn2nb(knearneigh(points, k = 8))
W <- nb2listw(nb_k8, style = "W")

OLS and Residual Check

Start with OLS.

Code 
ols <- lm(y ~ x, data = dat)
dat$olsResids <- residuals(ols)
moran.test(dat$olsResids, W)

    Moran I test under randomisation

data:  dat$olsResids  
weights: W    

Moran I statistic standard deviate = 7.0468, p-value = 9.152e-13
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.343949179      -0.013513514       0.002573207

Moran’s I is significantly positive on the OLS residuals. We know we have a problem. Run the LM tests.

LM Tests

Code 
lm.RStests(ols, listw = W, test = "all")

    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = y ~ x, data = dat)
test weights: W

RSerr = 39.288, df = 1, p-value = 3.657e-10


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = y ~ x, data = dat)
test weights: W

RSlag = 31.251, df = 1, p-value = 2.267e-08


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = y ~ x, data = dat)
test weights: W

adjRSerr = 8.5111, df = 1, p-value = 0.00353


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = y ~ x, data = dat)
test weights: W

adjRSlag = 0.47438, df = 1, p-value = 0.491


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = y ~ x, data = dat)
test weights: W

SARMA = 39.763, df = 2, p-value = 2.321e-09

Both RSerr and RSlag will be significant – a spatial error process can produce apparent lag structure because the autocorrelated errors bleed into the response. Look at the robust tests. adjRSerr is significant; adjRSlag is not. That points to the error model, which is the right answer since that is how we built the data.

Spatial Error Model

Code 
semFit <- errorsarlm(y ~ x, data = dat, listw = W)
summary(semFit)

Call:errorsarlm(formula = y ~ x, data = dat, listw = W)

Residuals:
      Min        1Q    Median        3Q       Max 
-1.832670 -0.516352 -0.011576  0.503105  2.692372 

Type: error 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value  Pr(>|z|)
(Intercept) 3.096759   0.336247  9.2098 < 2.2e-16
x           0.482327   0.086319  5.5877 2.301e-08

Lambda: 0.72049, LR test value: 23.41, p-value: 1.3092e-06
Asymptotic standard error: 0.10165
    z-value: 7.0877, p-value: 1.3636e-12
Wald statistic: 50.235, p-value: 1.3636e-12

Log likelihood: -94.04644 for error model
ML residual variance (sigma squared): 0.66157, (sigma: 0.81337)
Number of observations: 75 
Number of parameters estimated: 4 
AIC: 196.09, (AIC for lm: 217.5)

The coefficient table gives \(\hat\beta_0\) and \(\hat\beta_1\), interpreted exactly like OLS coefficients. The additional parameter is \(\lambda\), the spatial autocorrelation in the errors. A value near 1 means strong spatial error dependence; near 0 means the error model adds little over OLS. Given how we generated the data (\(\rho = 0.75\)), lambda should be substantial, and the coefficient estimates should land near our known values of \(\beta_0 = 3\) and \(\beta_1 = 0.5\).

Check the residuals.

Code 
dat$semResids <- residuals(semFit)
moran.test(dat$semResids, W)

    Moran I test under randomisation

data:  dat$semResids  
weights: W    

Moran I statistic standard deviate = 1.0261, p-value = 0.1524
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.038829368      -0.013513514       0.002601947

Moran’s I on the SEM residuals should be near zero. The spatial structure has been absorbed into \(\lambda\) and what is left looks like independent noise.

Spatial Lag Model

For comparison, fit the lag model on the same data.

Code 
slmFit <- lagsarlm(y ~ x, data = dat, listw = W)
summary(slmFit)

Call:lagsarlm(formula = y ~ x, data = dat, listw = W)

Residuals:
      Min        1Q    Median        3Q       Max 
-2.068155 -0.557181  0.051375  0.542870  2.708741 

Type: lag 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value  Pr(>|z|)
(Intercept)  1.14051    0.38178  2.9874  0.002814
x            0.43158    0.09069  4.7589 1.947e-06

Rho: 0.62184, LR test value: 18.187, p-value: 2.0028e-05
Asymptotic standard error: 0.11873
    z-value: 5.2375, p-value: 1.6276e-07
Wald statistic: 27.431, p-value: 1.6276e-07

Log likelihood: -96.65795 for lag model
ML residual variance (sigma squared): 0.72882, (sigma: 0.85371)
Number of observations: 75 
Number of parameters estimated: 4 
AIC: 201.32, (AIC for lm: 217.5)
LM test for residual autocorrelation
test value: 10.136, p-value: 0.0014539

The lag model estimates \(\rho\), the spatial autoregression of the response. Because the data were generated under an error process, this is the wrong model. It fits reasonably (the models are related) but the LM tests already pointed us elsewhere, and the AIC confirms it.

Code 
AIC(semFit, slmFit)
       df      AIC
semFit  4 196.0929
slmFit  4 201.3159

The SEM should win. This is a good workflow to internalize: LM tests first, then model fitting, then AIC as a sanity check.

Now let’s run this on data we didn’t build.

Bird Species Richness in Mexico

The birdRichnessMexico dataset has 200 point locations across Mexico with bird species richness (nSpecies) and six WorldClim bioclimatic variables extracted at each location. Maps use tmap (Tennekes 2018). Two predictors have strong ecological justification: mean annual precipitation (map) and temperature annual range (tempRange). Wetter areas support more species; thermally variable areas support fewer. That’s the energy-water hypothesis and Rapoport’s rule operating across the Mexican landscape.

Code 
library(tmap)
birds_sf <- readRDS("data/birdRichnessMexico.rds")
coords_birds <- st_coordinates(birds_sf)

nb_birds <- knn2nb(knearneigh(coords_birds, k = 8))
W_birds  <- nb2listw(nb_birds, style = "W")

Map the richness surface first.

Code 
tmap_mode("view")
tm_shape(birds_sf) +
  tm_symbols(col = "nSpecies", palette = "viridis",
             title.col = "Species", size = 0.4)

OLS and Residual Check

Code 
ols_birds <- lm(nSpecies ~ map + tempRange, data = birds_sf)
summary(ols_birds)

Call:
lm(formula = nSpecies ~ map + tempRange, data = birds_sf)

Residuals:
     Min       1Q   Median       3Q      Max 
-136.424  -13.629   -2.138   15.146   92.269 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 213.383856  15.079562  14.151  < 2e-16 ***
map           0.078438   0.004528  17.324  < 2e-16 ***
tempRange    -1.758331   0.473349  -3.715 0.000265 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 27.74 on 197 degrees of freedom
Multiple R-squared:  0.8009,    Adjusted R-squared:  0.7989 
F-statistic: 396.3 on 2 and 197 DF,  p-value: < 2.2e-16

R² of 0.80. Both predictors are significant with the expected signs. Now check the residuals.

Code 
birds_sf$ols_resids <- residuals(ols_birds)
moran.test(birds_sf$ols_resids, W_birds)

    Moran I test under randomisation

data:  birds_sf$ols_resids  
weights: W_birds    

Moran I statistic standard deviate = 12.42, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.397871539      -0.005025126       0.001052372

Moran’s I is 0.40, strongly significant. The model explains 80% of the variance but the residuals carry substantial spatial signal. Map them to see where the model is over- and under-predicting.

Code 
tm_shape(birds_sf) +
  tm_symbols(col = "ols_resids",
             palette = "-RdBu",
             midpoint = 0,
             title.col = "Residual",
             size = 0.4)

The pattern is not random. There are regions where richness is consistently higher than the climate model predicts and others where it’s consistently lower. That’s what Moran’s I is flagging, and the map shows you where.

LM Tests

Code 
lm.RStests(ols_birds, listw = W_birds, test = "all")

    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = nSpecies ~ map + tempRange, data = birds_sf)
test weights: W_birds

RSerr = 140.81, df = 1, p-value < 2.2e-16


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = nSpecies ~ map + tempRange, data = birds_sf)
test weights: W_birds

RSlag = 33.448, df = 1, p-value = 7.321e-09


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = nSpecies ~ map + tempRange, data = birds_sf)
test weights: W_birds

adjRSerr = 107.73, df = 1, p-value < 2.2e-16


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = nSpecies ~ map + tempRange, data = birds_sf)
test weights: W_birds

adjRSlag = 0.37184, df = 1, p-value = 0.542


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = nSpecies ~ map + tempRange, data = birds_sf)
test weights: W_birds

SARMA = 141.18, df = 2, p-value < 2.2e-16

Both RSerr and RSlag are significant. Go to the robust tests. adjRSerr is large and stays significant; adjRSlag collapses to near zero. The error model is the right fit.

Spatial Error Model

Code 
sem_birds <- errorsarlm(nSpecies ~ map + tempRange,
                        data = birds_sf, listw = W_birds)
summary(sem_birds)

Call:errorsarlm(formula = nSpecies ~ map + tempRange, data = birds_sf, 
    listw = W_birds)

Residuals:
      Min        1Q    Median        3Q       Max 
-133.6497  -12.8050    0.7134   11.8756   60.3406 

Type: error 
Coefficients: (asymptotic standard errors) 
               Estimate  Std. Error z value  Pr(>|z|)
(Intercept) 245.0812346  20.1854247 12.1415 < 2.2e-16
map           0.0727603   0.0052084 13.9699 < 2.2e-16
tempRange    -2.7528487   0.6849927 -4.0188  5.85e-05

Lambda: 0.72462, LR test value: 74.723, p-value: < 2.22e-16
Asymptotic standard error: 0.059692
    z-value: 12.139, p-value: < 2.22e-16
Wald statistic: 147.36, p-value: < 2.22e-16

Log likelihood: -909.5156 for error model
ML residual variance (sigma squared): 478.86, (sigma: 21.883)
Number of observations: 200 
Number of parameters estimated: 5 
AIC: 1829, (AIC for lm: 1901.8)

Lambda is 0.72, highly significant. There is strong spatial autocorrelation in the errors above and beyond what precipitation and temperature range explain. Check the residuals.

Code 
birds_sf$sem_resids <- residuals(sem_birds)
moran.test(birds_sf$sem_resids, W_birds)

    Moran I test under randomisation

data:  birds_sf$sem_resids  
weights: W_birds    

Moran I statistic standard deviate = 0.43109, p-value = 0.3332
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.008848150      -0.005025126       0.001035655

Moran’s I drops to near zero. The spatial structure has been absorbed into lambda. Compare AIC.

Code 
AIC(ols_birds, sem_birds)
          df      AIC
ols_birds  4 1901.754
sem_birds  5 1829.031

A drop of about 73 units. The error model is substantially better.

What Changed

The coefficients shift after accounting for spatial error dependence. The effect of tempRange gets stronger: OLS estimated -1.76 species per degree of range; the SEM estimates -2.75. OLS was underestimating the temperature effect because spatially autocorrelated errors were confounding the estimate. The precipitation coefficient shifts less, from 0.078 to 0.073, which makes sense given how dominant that predictor is.

What does lambda = 0.72 mean ecologically? Precipitation and temperature range explain most of the richness gradient, but there’s substantial residual spatial structure they don’t capture. Some of it is almost certainly habitat heterogeneity: the Sierra Madre Occidental and Oriental are centers of endemism and provide topographic complexity that flat climate variables miss. Some may reflect historical biogeography – proximity to the Neotropical source pool in southern Mexico leaves a signature that current climate alone doesn’t explain. The SEM doesn’t tell you what the missing variable is. It tells you that something spatially structured is out there and that ignoring it was biasing your inference.

A Note on Interpreting Lag Model Coefficients

In an OLS or SEM, \(\hat\beta_1 = 0.5\) means: a one-unit increase in \(x\) is associated with a 0.5-unit increase in \(y\). That is a total effect.

In a SLM, a change in \(x_i\) affects \(y_i\) directly, but that change in \(y_i\) then propagates to neighboring locations because their \(y\) values depend on \(y_i\) via the lag term. Those ripples feed back. The \(\hat\beta_1\) from lagsarlm is only the direct local effect, before the spatial propagation plays out. The full marginal effect, including spillovers, requires computing impacts from the fitted model.

Code 
impacts(slmFit, listw = W, R = 500)
Impact measures (lag, exact):
           Direct  Indirect    Total
x dy/dx 0.4626068 0.6786596 1.141266

The output decomposes effects into direct (same-location), indirect (spillover to neighbors), and total. For most environmental science applications, the total effects are what you want to report if you are using a lag model. If the indirect effects are small relative to the direct effects, the distinction matters less in practice, but check.

The lag model is making a more specific and more mechanistic claim about the data-generating process than the error model. The error model is the cleaner choice when you just want valid inference on your predictors and space is nuisance variation.

GLS vs. SAR: When to Use Which

They are often interchangeable in practice for the spatial error case.

GLS uses a variogram, a continuous function of distance, to model the error correlation structure. It captures the range, sill, and nugget of the spatial dependence and is very flexible. SAR models use a discrete neighbor structure: you define who counts as a neighbor and the autocorrelation operates only among neighbors. GLS tends to work better when you have a clear distance-decay structure and a reasonable sample size to fit the variogram. SAR works better when your data are naturally defined by adjacency or neighborhood (watershed catchments, habitat patches, gridded climate data) and when you want to lean on the weights matrix you have already built for exploratory work like LISA.

The spatial lag model has no real GLS equivalent. If the biology or ecology predicts a spillover or diffusion process, fit the lag model. Just be careful with the coefficients.

Wrapping Up

GLS and SAR are both corrections to the same problem: OLS residuals that carry spatial signal the model hasn’t explained. Which one you reach for depends on your data and your question. For most field-based environmental science, GLS is the workhorse. It’s flexible, connects directly to the variogram, and the output is easy to interpret. SAR is worth knowing when the spatial structure is better described by a neighbor graph than a continuous distance decay, or when there’s a credible spillover process in the response itself. That’s a real mechanism sometimes – think about how disease spreads, or how animal behavior in one patch influences neighboring patches. But for lead in floodplain soil, space is probably just a proxy for flood exposure. Use GLS and sleep well.

Exercises

The Meuse River data contains several candidate predictors for soil lead contamination. river_dist_m is the distance from the river channel in meters and ffreq is a flooding frequency class (a factor with three levels: 1 = most frequent, 3 = least frequent). Both are plausible drivers of lead deposition: proximity to the channel means exposure to recent flood sediment, and flooding frequency determines how much sediment accumulates.

Code 
meuse2      <- readRDS("data/meuse2.Rds")
meuse.grid2 <- readRDS("data/meuse.grid2.Rds")
meuse_sf <- st_as_sf(meuse2, coords = c("x", "y"), crs = 28992)
meuse_sf$log_lead <- log(meuse_sf$lead)
# get river distance and flooding frequency from the grid
covars_sf <- st_as_sf(meuse.grid2, coords = c("x","y"), crs = 28992) %>%
  select(ffreq, river_dist_m)
meuse_sf <- st_join(meuse_sf, covars_sf, join = st_nearest_feature) %>%
  mutate(ffreq = factor(ffreq))

  1. Build a k-nearest-neighbor weights matrix with \(k = 8\).

  2. Fit OLS of log_lead ~ river_dist_m + ffreq. Interpret the coefficients – do the signs make ecological sense? Run Moran’s I on the residuals.

  3. Run the LM tests. Both RSerr and RSlag will be significant. Look at the robust tests to break the tie. Which model do they recommend?

  4. Fit both the SEM and the SLM. Check the residuals of each using Moran’s I. One model cleans up the spatial structure; the other doesn’t. Which is which, and what does that tell you?

  5. Compare AIC across OLS, SEM, and SLM. Does AIC agree with the LM tests?

  6. Compare coefficients on river_dist_m and ffreq between OLS and the SEM. The point estimates shift and the standard errors change. What does the shift in standard errors tell you about what OLS was doing wrong?

  7. Is the spatial error model or the spatial lag model more ecologically defensible for soil lead? What process would the lag model imply, and is that a plausible mechanism for how lead ends up in floodplain soil?

  8. In the bird richness example, lambda = 0.72 after accounting for precipitation and temperature range. What do you think is driving the residual spatial structure? Name at least two candidate variables not in the model and explain why they might be spatially structured across Mexico.

Further Reading

Dormann et al. (2007) sections 3.3 and 3.4 cover spatial error models and related approaches. Bivand (2008) chapter 9 goes deep on the SAR model family and the spatialreg functions; it’s a technical chapter, but worth returning to once you’ve worked through the examples here.