Moran eigenvector spatial filtering

spatial

regression

autocorrelation

glm

Turn a spatial weights matrix into orthogonal map patterns, add them to a Poisson GLM to remove residual autocorrelation, and see where filtering falls short.

Author

Tidy Ecology

Published

2026-08-13

The SAR and CAR models are built for Gaussian responses. Ecological areal data are often counts or presences, where you would rather stay inside a Poisson or binomial generalised linear model. Fitting a full autoregressive likelihood for those families is awkward. Moran eigenvector spatial filtering offers a lighter alternative: it distils the neighbourhood structure into a handful of orthogonal map patterns, then drops them into an ordinary GLM as extra covariates. The spatial autocorrelation is soaked up by the patterns, and everything else about the model stays familiar.

This tutorial builds the eigenvector maps by hand, confirms that each one’s Moran’s I equals a scaled eigenvalue, and uses forward selection to clean the residual autocorrelation out of a Poisson regression. It also does something many treatments avoid: it measures, across many simulated datasets, how well filtering actually controls the false positives that spatial autocorrelation causes. The answer is encouraging but not triumphant.

Map patterns from a weights matrix

Take the binary adjacency matrix A, centre it on both sides with C = I - 11'/n, and find the eigenvectors of C A C. Those eigenvectors, the Moran eigenvector maps (MEM), are mutually orthogonal spatial patterns. The one with the largest eigenvalue is the broadest gradient across the map; later ones describe progressively finer structure; and eigenvectors with negative eigenvalues encode negative autocorrelation, the checkerboards. The tidy fact that makes them useful is that the Moran’s I of eigenvector k is exactly (n / sum(A)) times its eigenvalue.

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
}
one <- rep(1, n); Cn <- diag(n) - outer(one, one) / n
eg  <- eigen(Cn %*% A %*% Cn, symmetric = TRUE)
keep <- abs(eg$values) > 1e-8
E <- eg$vectors[, keep]                       # the eigenvector maps
moranE <- (n / sum(A)) * eg$values[keep]      # Moran's I of each map
broad <- which(moranE > 0)[order(-moranE[moranE > 0])]  # broad patterns first
mem_npos <- length(broad)

A quick check confirms the eigenvalue identity: the analytic Moran’s I of a map agrees with a direct computation.

k <- broad[1]
c(from_eigenvalue = moranE[k], from_moran_test = moran_test(E[, k], A)$I)

from_eigenvalue from_moran_test 
        1.02337         1.02337

The two agree to machine precision, and there are 189 maps carrying positive autocorrelation out of 399 in total.

Two grid maps side by side, each a Moran eigenvector map. The left panel is a broad-scale pattern, a smooth gradient across the whole grid, labelled with a high Moran's I. The right panel is a much finer pattern with many small alternating patches, labelled with a lower Moran's I.

A spatial red herring

Here is the trap that filtering is meant to defuse. We simulate a spatially structured covariate that has no real effect on the response, alongside an unmeasured smooth field that does drive the counts. Because both the covariate and the response carry spatial structure, a naive analysis is prone to see a relationship that is not there.

set.seed(4233)
x   <- as.numeric(scale(1.1 * sin(pi * coords$col / 11) -
                        0.9 * cos(pi * coords$row / 12) + rnorm(n, 0, 0.5)))
fsp <- as.numeric(scale(1.2 * sin(pi * coords$col / 8 + 0.6) +
                        1.0 * cos(pi * coords$row / 7 - 0.4)))  # unmeasured field
b0 <- 1.4; b1_true <- 0
set.seed(20)
y <- rpois(n, exp(b0 + b1_true * x + 0.85 * fsp))
c(mean = mean(y), zeros = sum(y == 0), max = max(y))

   mean   zeros     max 
 5.7875 30.0000 27.0000

Fit the count model that ignores space and inspect both the covariate and the residual autocorrelation.

g0 <- glm(y ~ x, family = poisson)
s0 <- summary(g0)$coefficients
mem_base_b1    <- s0["x", "Estimate"]
mem_base_se    <- s0["x", "Std. Error"]
mem_base_z     <- s0["x", "z value"]
mem_base_p     <- s0["x", "Pr(>|z|)"]
mem_base_moran <- moran_test(residuals(g0, type = "pearson"), A)$I
round(c(b1 = mem_base_b1, se = mem_base_se, z = mem_base_z,
        p = mem_base_p, resid_moran = mem_base_moran), 4)

         b1          se           z           p resid_moran 
    -0.0738      0.0209     -3.5326      0.0004      0.7178

The covariate looks convincingly significant, with a z of -3.53 and a p-value of 4.1^{-4}, yet its true effect is zero. The Pearson residuals carry a Moran’s I of 0.72, a blazing signal of leftover spatial structure. This is the red herring that Diniz-Filho et al. (2003) warned about: autocorrelation manufactures false confidence, and can even flip the sign of a relationship.

Filtering it out

Forward selection adds eigenvector maps one at a time, each time choosing the map most correlated with the current residuals, and stops once the residual autocorrelation is no longer significant. We restrict the candidates to the broad-scale maps, where a smooth field lives.

cand <- broad[1:80]; sel <- integer(0); r <- resid(lm(y ~ x))
for (s in 1:25) {
  if (moran_test(r, A)$p > 0.05) break
  cc <- abs(cor(E[, setdiff(cand, sel)], r))
  sel <- c(sel, setdiff(cand, sel)[which.max(cc)])
  r <- resid(lm(y ~ x + E[, sel]))
}
gm <- glm(y ~ x + E[, sel], family = poisson)
sm <- summary(gm)$coefficients
mem_nsel      <- length(sel)
mem_filt_b1   <- sm["x", "Estimate"]
mem_filt_se   <- sm["x", "Std. Error"]
mem_filt_z    <- sm["x", "z value"]
mem_filt_p    <- sm["x", "Pr(>|z|)"]
mem_filt_moran <- moran_test(residuals(gm, type = "pearson"), A)$I
round(c(n_maps = mem_nsel, b1 = mem_filt_b1, se = mem_filt_se,
        z = mem_filt_z, p = mem_filt_p, resid_moran = mem_filt_moran), 4)

     n_maps          b1          se           z           p resid_moran 
    18.0000      0.0047      0.0233      0.2005      0.8411      0.0359

Adding 18 maps pulls the residual Moran’s I down from 0.72 to 0.036, no longer significant, and the spurious effect evaporates: the covariate’s p-value climbs from 4.1^{-4} to 0.84. The GLM machinery never changed. We simply gave it enough spatial patterns to explain the clustering that the covariate was otherwise blamed for.

Two panels tracking the forward selection. The left panel is a line chart showing residual Moran's I falling steadily from about 0.72 towards near zero as more eigenvector maps are added. The right panel shows the covariate p-value rising from far below 0.05 up across the 0.05 line and settling around 0.5 as the spurious signal is absorbed.

Does filtering really fix the inference?

Removing autocorrelation from one dataset is reassuring but not proof. The sterner test is to repeat the whole procedure on many datasets where the covariate genuinely has no effect, and count how often each method wrongly declares it significant. We compare the naive Poisson GLM with adaptive filtering, selecting maps afresh for every dataset.

select_mem <- function(y, x, cand, maxK = 25) {
  sel <- integer(0); r <- resid(lm(y ~ x))
  for (s in 1:maxK) {
    if (moran_test(r, A)$p > 0.05) break
    cc <- abs(cor(E[, setdiff(cand, sel)], r))
    sel <- c(sel, setdiff(cand, sel)[which.max(cc)]); r <- resid(lm(y ~ x + E[, sel]))
  }
  sel
}
mkx <- function(seed) { set.seed(seed)
  as.numeric(scale(1.1 * sin(pi * coords$col / 11) -
                   0.9 * cos(pi * coords$row / 12) + rnorm(n, 0, 0.5))) }
mkf <- function(seed) { set.seed(seed)
  as.numeric(scale(1.2 * sin(pi * coords$col / 8 + 0.6) +
                   1.0 * cos(pi * coords$row / 7 - 0.4))) }
set.seed(505); B <- 300; naive <- adapt <- logical(B); nsel <- integer(B)
for (b in 1:B) {
  xb <- mkx(10000 + b); fb <- mkf(20000 + b); yb <- rpois(n, exp(b0 + 0 * xb + 0.85 * fb))
  naive[b] <- summary(glm(yb ~ xb, family = poisson))$coefficients["xb", "Pr(>|z|)"] < 0.05
  se <- select_mem(yb, xb, broad[1:80]); nsel[b] <- length(se)
  adapt[b] <- summary(glm(yb ~ xb + E[, se], family = poisson))$coefficients["xb", "Pr(>|z|)"] < 0.05
}
mem_mc_naive <- mean(naive); mem_mc_adapt <- mean(adapt); mem_mc_medk <- median(nsel)
round(c(naive = mem_mc_naive, filtered = mem_mc_adapt, median_maps = mem_mc_medk), 3)

      naive    filtered median_maps 
      0.903       0.113      20.000

The naive GLM is a disaster: it flags the null covariate as significant 90 per cent of the time, not the nominal five. Adaptive filtering, using a median of 20 maps, cuts that to 11 per cent. That is a large improvement, roughly a fivefold reduction, but it is honestly still more than triple the nominal rate. Filtering helps a great deal; it does not fully cure the problem.

This matches what the literature found the hard way. Bini et al. (2009) compared spatial and non-spatial regression across dozens of real datasets and concluded that the shifts in coefficients are largely idiosyncratic, with no reliable rule for when a spatial method will behave. Eigenvector filtering reduces the red herring without eliminating it, and the number of maps you keep is a genuine tuning choice that trades false positives against power.

What to carry forward

Moran eigenvector filtering is the pragmatic choice when your response is a count or a presence and a full SAR or CAR likelihood is more than you want to fit. It converts the neighbourhood into orthogonal map patterns, and a modest selection of them, added as covariates, strips most of the spatial autocorrelation out of the residuals. That whitening is real and worth having.

Keep two things in view. The maps are patterns, not mechanisms; they describe where the response clusters, not why, so they support honest inference rather than causal claims. And the control they buy is partial: false-positive rates fall sharply but not to nominal, and the choice of how many maps to keep matters. When the covariate you care about is itself spatially smooth, those same maps can start absorbing its signal, which is the confounding problem that checking a spatial regression confronts directly.

References

Bini, L. M., Diniz-Filho, J. A. F., Rangel, T. F. L. V. B., Akre, T. S. B., et al. (2009). Coefficient shifts in geographical ecology: an empirical evaluation of spatial and non-spatial regression. Ecography, 32(2), 193-204. https://doi.org/10.1111/j.1600-0587.2009.05717.x

Diniz-Filho, J. A. F., Bini, L. M., & Hawkins, B. A. (2003). Spatial autocorrelation and red herrings in geographical ecology. Global Ecology and Biogeography, 12(1), 53-64. https://doi.org/10.1046/j.1466-822X.2003.00322.x

Dray, S., Legendre, P., & Peres-Neto, P. R. (2006). Spatial modelling: a comprehensive framework for principal coordinate analysis of neighbour matrices (PCNM). Ecological Modelling, 196(3-4), 483-493. https://doi.org/10.1016/j.ecolmodel.2006.02.015

Griffith, D. A., & Peres-Neto, P. R. (2006). Spatial modeling in ecology: the flexibility of eigenfunction spatial analyses. Ecology, 87(10), 2603-2613.

--- title: "Moran eigenvector spatial filtering" description: "Turn a spatial weights matrix into orthogonal map patterns, add them to a Poisson GLM to remove residual autocorrelation, and see where filtering falls short." date: 2026-08-13 categories: [spatial, regression, autocorrelation, glm] image: thumbnail.png --- The [SAR](../spatial-autoregressive-models-sar/) and [CAR](../conditional-autoregressive-car-models/) models are built for Gaussian responses. Ecological areal data are often counts or presences, where you would rather stay inside a Poisson or binomial generalised linear model. Fitting a full autoregressive likelihood for those families is awkward. Moran eigenvector spatial filtering offers a lighter alternative: it distils the neighbourhood structure into a handful of orthogonal *map patterns*, then drops them into an ordinary GLM as extra covariates. The spatial autocorrelation is soaked up by the patterns, and everything else about the model stays familiar. This tutorial builds the eigenvector maps by hand, confirms that each one's Moran's I equals a scaled eigenvalue, and uses forward selection to clean the residual autocorrelation out of a Poisson regression. It also does something many treatments avoid: it measures, across many simulated datasets, how well filtering actually controls the false positives that spatial autocorrelation causes. The answer is encouraging but not triumphant. ## Map patterns from a weights matrix Take the binary adjacency matrix `A`, centre it on both sides with `C = I - 11'/n`, and find the eigenvectors of `C A C`. Those eigenvectors, the *Moran eigenvector maps* (MEM), are mutually orthogonal spatial patterns. The one with the largest eigenvalue is the broadest gradient across the map; later ones describe progressively finer structure; and eigenvectors with negative eigenvalues encode negative autocorrelation, the checkerboards. The tidy fact that makes them useful is that the Moran's I of eigenvector `k` is exactly `(n / sum(A))` times its eigenvalue. ```{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")) } memmap <- function(coords, value, title, subtitle = NULL) { df <- data.frame(col = coords$col, row = coords$row, v = value) ggplot(df, aes(col, row, fill = v)) + geom_tile() + coord_equal() + scale_fill_gradient2(low = "#3b4cc0", mid = "white", high = "#b40426", midpoint = 0) + 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)) } 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)) } ``` ```{r mem} 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 } one <- rep(1, n); Cn <- diag(n) - outer(one, one) / n eg <- eigen(Cn %*% A %*% Cn, symmetric = TRUE) keep <- abs(eg$values) > 1e-8 E <- eg$vectors[, keep] # the eigenvector maps moranE <- (n / sum(A)) * eg$values[keep] # Moran's I of each map broad <- which(moranE > 0)[order(-moranE[moranE > 0])] # broad patterns first mem_npos <- length(broad) ``` A quick check confirms the eigenvalue identity: the analytic Moran's I of a map agrees with a direct computation. ```{r mem-check} k <- broad[1] c(from_eigenvalue = moranE[k], from_moran_test = moran_test(E[, k], A)$I) ``` The two agree to machine precision, and there are `r mem_npos` maps carrying positive autocorrelation out of `r n - 1` in total. ```{r fig1-maps} #| echo: false #| fig-alt: "Two grid maps side by side, each a Moran eigenvector map. The left panel is a broad-scale pattern, a smooth gradient across the whole grid, labelled with a high Moran's I. The right panel is a much finer pattern with many small alternating patches, labelled with a lower Moran's I." kb <- broad[1]; kf <- broad[55] two_panel( memmap(coords, E[, kb], "A broad-scale map", paste0("Moran's I = ", round(moranE[kb], 2))), memmap(coords, E[, kf], "A fine-scale map", paste0("Moran's I = ", round(moranE[kf], 2))) ) ``` ## A spatial red herring Here is the trap that filtering is meant to defuse. We simulate a spatially structured covariate that has *no* real effect on the response, alongside an unmeasured smooth field that does drive the counts. Because both the covariate and the response carry spatial structure, a naive analysis is prone to see a relationship that is not there. ```{r simulate} set.seed(4233) x <- as.numeric(scale(1.1 * sin(pi * coords$col / 11) - 0.9 * cos(pi * coords$row / 12) + rnorm(n, 0, 0.5))) fsp <- as.numeric(scale(1.2 * sin(pi * coords$col / 8 + 0.6) + 1.0 * cos(pi * coords$row / 7 - 0.4))) # unmeasured field b0 <- 1.4; b1_true <- 0 set.seed(20) y <- rpois(n, exp(b0 + b1_true * x + 0.85 * fsp)) c(mean = mean(y), zeros = sum(y == 0), max = max(y)) ``` Fit the count model that ignores space and inspect both the covariate and the residual autocorrelation. ```{r base-glm} g0 <- glm(y ~ x, family = poisson) s0 <- summary(g0)$coefficients mem_base_b1 <- s0["x", "Estimate"] mem_base_se <- s0["x", "Std. Error"] mem_base_z <- s0["x", "z value"] mem_base_p <- s0["x", "Pr(>|z|)"] mem_base_moran <- moran_test(residuals(g0, type = "pearson"), A)$I round(c(b1 = mem_base_b1, se = mem_base_se, z = mem_base_z, p = mem_base_p, resid_moran = mem_base_moran), 4) ``` The covariate looks convincingly significant, with a z of `r round(mem_base_z, 2)` and a p-value of `r signif(mem_base_p, 2)`, yet its true effect is zero. The Pearson residuals carry a Moran's I of `r round(mem_base_moran, 2)`, a blazing signal of leftover spatial structure. This is the red herring that Diniz-Filho et al. (2003) warned about: autocorrelation manufactures false confidence, and can even flip the sign of a relationship. ## Filtering it out Forward selection adds eigenvector maps one at a time, each time choosing the map most correlated with the current residuals, and stops once the residual autocorrelation is no longer significant. We restrict the candidates to the broad-scale maps, where a smooth field lives. ```{r filter} cand <- broad[1:80]; sel <- integer(0); r <- resid(lm(y ~ x)) for (s in 1:25) { if (moran_test(r, A)$p > 0.05) break cc <- abs(cor(E[, setdiff(cand, sel)], r)) sel <- c(sel, setdiff(cand, sel)[which.max(cc)]) r <- resid(lm(y ~ x + E[, sel])) } gm <- glm(y ~ x + E[, sel], family = poisson) sm <- summary(gm)$coefficients mem_nsel <- length(sel) mem_filt_b1 <- sm["x", "Estimate"] mem_filt_se <- sm["x", "Std. Error"] mem_filt_z <- sm["x", "z value"] mem_filt_p <- sm["x", "Pr(>|z|)"] mem_filt_moran <- moran_test(residuals(gm, type = "pearson"), A)$I round(c(n_maps = mem_nsel, b1 = mem_filt_b1, se = mem_filt_se, z = mem_filt_z, p = mem_filt_p, resid_moran = mem_filt_moran), 4) ``` Adding `r mem_nsel` maps pulls the residual Moran's I down from `r round(mem_base_moran, 2)` to `r round(mem_filt_moran, 3)`, no longer significant, and the spurious effect evaporates: the covariate's p-value climbs from `r signif(mem_base_p, 2)` to `r round(mem_filt_p, 2)`. The GLM machinery never changed. We simply gave it enough spatial patterns to explain the clustering that the covariate was otherwise blamed for. ```{r fig2-results} #| echo: false #| fig-alt: "Two panels tracking the forward selection. The left panel is a line chart showing residual Moran's I falling steadily from about 0.72 towards near zero as more eigenvector maps are added. The right panel shows the covariate p-value rising from far below 0.05 up across the 0.05 line and settling around 0.5 as the spurious signal is absorbed." sel2 <- integer(0); rr <- resid(lm(y ~ x)) tr_moran <- moran_test(rr, A)$I tr_p <- summary(glm(y ~ x, family = poisson))$coefficients["x", "Pr(>|z|)"] for (s in seq_along(sel)) { sel2 <- c(sel2, sel[s]); rr <- resid(lm(y ~ x + E[, sel2])) tr_moran <- c(tr_moran, moran_test(rr, A)$I) tr_p <- c(tr_p, summary(glm(y ~ x + E[, sel2], family = poisson))$coefficients["x", "Pr(>|z|)"]) } steps <- 0:length(sel) g1 <- ggplot(data.frame(step = steps, moran = tr_moran), aes(step, moran)) + geom_hline(yintercept = 0, colour = "grey70", linetype = 2) + geom_line(colour = "#275139", linewidth = 0.9) + geom_point(colour = "#275139", size = 1.5) + labs(title = "Residuals whiten as maps are added", subtitle = "forward selection path", x = "eigenvector maps added", y = "residual Moran's I") + theme_te() + theme(axis.text = element_text(size = 8), panel.grid.major.y = element_line(colour = "grey92")) g2 <- ggplot(data.frame(step = steps, p = tr_p), aes(step, p)) + geom_hline(yintercept = 0.05, colour = "grey40", linetype = 2) + geom_line(colour = "#c98a1a", linewidth = 0.9) + geom_point(colour = "#c98a1a", size = 1.5) + annotate("text", x = max(steps) * 0.6, y = 0.09, label = "0.05", size = 3, colour = "grey40") + labs(title = "The spurious effect dies", subtitle = "covariate has no true effect", x = "eigenvector maps added", y = "covariate p-value") + theme_te() + theme(axis.text = element_text(size = 8), panel.grid.major.y = element_line(colour = "grey92")) two_panel(g1, g2) ``` ## Does filtering really fix the inference? Removing autocorrelation from one dataset is reassuring but not proof. The sterner test is to repeat the whole procedure on many datasets where the covariate genuinely has no effect, and count how often each method wrongly declares it significant. We compare the naive Poisson GLM with adaptive filtering, selecting maps afresh for every dataset. ```{r mc} select_mem <- function(y, x, cand, maxK = 25) { sel <- integer(0); r <- resid(lm(y ~ x)) for (s in 1:maxK) { if (moran_test(r, A)$p > 0.05) break cc <- abs(cor(E[, setdiff(cand, sel)], r)) sel <- c(sel, setdiff(cand, sel)[which.max(cc)]); r <- resid(lm(y ~ x + E[, sel])) } sel } mkx <- function(seed) { set.seed(seed) as.numeric(scale(1.1 * sin(pi * coords$col / 11) - 0.9 * cos(pi * coords$row / 12) + rnorm(n, 0, 0.5))) } mkf <- function(seed) { set.seed(seed) as.numeric(scale(1.2 * sin(pi * coords$col / 8 + 0.6) + 1.0 * cos(pi * coords$row / 7 - 0.4))) } set.seed(505); B <- 300; naive <- adapt <- logical(B); nsel <- integer(B) for (b in 1:B) { xb <- mkx(10000 + b); fb <- mkf(20000 + b); yb <- rpois(n, exp(b0 + 0 * xb + 0.85 * fb)) naive[b] <- summary(glm(yb ~ xb, family = poisson))$coefficients["xb", "Pr(>|z|)"] < 0.05 se <- select_mem(yb, xb, broad[1:80]); nsel[b] <- length(se) adapt[b] <- summary(glm(yb ~ xb + E[, se], family = poisson))$coefficients["xb", "Pr(>|z|)"] < 0.05 } mem_mc_naive <- mean(naive); mem_mc_adapt <- mean(adapt); mem_mc_medk <- median(nsel) round(c(naive = mem_mc_naive, filtered = mem_mc_adapt, median_maps = mem_mc_medk), 3) ``` The naive GLM is a disaster: it flags the null covariate as significant `r round(100 * mem_mc_naive, 0)` per cent of the time, not the nominal five. Adaptive filtering, using a median of `r mem_mc_medk` maps, cuts that to `r round(100 * mem_mc_adapt, 0)` per cent. That is a large improvement, roughly a fivefold reduction, but it is honestly still more than triple the nominal rate. Filtering helps a great deal; it does not fully cure the problem. This matches what the literature found the hard way. Bini et al. (2009) compared spatial and non-spatial regression across dozens of real datasets and concluded that the shifts in coefficients are largely idiosyncratic, with no reliable rule for when a spatial method will behave. Eigenvector filtering reduces the red herring without eliminating it, and the number of maps you keep is a genuine tuning choice that trades false positives against power. ## What to carry forward Moran eigenvector filtering is the pragmatic choice when your response is a count or a presence and a full SAR or CAR likelihood is more than you want to fit. It converts the neighbourhood into orthogonal map patterns, and a modest selection of them, added as covariates, strips most of the spatial autocorrelation out of the residuals. That whitening is real and worth having. Keep two things in view. The maps are patterns, not mechanisms; they describe where the response clusters, not why, so they support honest inference rather than causal claims. And the control they buy is partial: false-positive rates fall sharply but not to nominal, and the choice of how many maps to keep matters. When the covariate you care about is itself spatially smooth, those same maps can start absorbing its signal, which is the confounding problem that [checking a spatial regression](../checking-a-spatial-regression/) confronts directly. ## Related tutorials - [Spatial lag and spatial error models (SAR)](../spatial-autoregressive-models-sar/) fit the full autoregressive likelihood for Gaussian responses. - [Conditional autoregressive (CAR) models](../conditional-autoregressive-car-models/) build the same dependence from local conditionals. - [Checking a spatial regression](../checking-a-spatial-regression/) tests which spatial model fits and exposes spatial confounding. - [Spatial cross-validation for distribution models](../spatial-cross-validation-sdm/) shows the validation side of the same autocorrelation problem. ## References Bini, L. M., Diniz-Filho, J. A. F., Rangel, T. F. L. V. B., Akre, T. S. B., et al. (2009). Coefficient shifts in geographical ecology: an empirical evaluation of spatial and non-spatial regression. *Ecography*, 32(2), 193-204. https://doi.org/10.1111/j.1600-0587.2009.05717.x Diniz-Filho, J. A. F., Bini, L. M., & Hawkins, B. A. (2003). Spatial autocorrelation and red herrings in geographical ecology. *Global Ecology and Biogeography*, 12(1), 53-64. https://doi.org/10.1046/j.1466-822X.2003.00322.x Dray, S., Legendre, P., & Peres-Neto, P. R. (2006). Spatial modelling: a comprehensive framework for principal coordinate analysis of neighbour matrices (PCNM). *Ecological Modelling*, 196(3-4), 483-493. https://doi.org/10.1016/j.ecolmodel.2006.02.015 Griffith, D. A., & Peres-Neto, P. R. (2006). Spatial modeling in ecology: the flexibility of eigenfunction spatial analyses. *Ecology*, 87(10), 2603-2613.