Generalised least squares for spatial data in R

GLS

spatial

regression

ecology tutorial

When regression residuals are spatially autocorrelated, OLS standard errors mislead. Fit a GLS model in nlme, pick a correlation structure by AIC, and compare.

Author

Tidy Ecology

Published

2026-05-06

Modified

2026-11-05

Corrected 5 November 2026: the post said the spread of OLS slopes across replicate landscapes is close to three times the standard error a single fit reports, and that OLS scatters more than GLS, without computing either. A new chunk now repeats the simulation two hundred times and reports the ratio, the coverage of both sets of intervals and the GLS spread.

Corrected 14 October 2026: the text said the GLS intervals cover the true values because their standard errors are wider; it now says that GLS also makes the estimate more precise and that here the intervals cover the truth mainly because the estimates moved towards it. The variogram caption no longer blames the noisy far points on a shortage of plot pairs, since every point averages about the same number of pairs; the farthest pairs draw on fewer distinct plots.

Corrected 27 September 2026: the variogram text said the semivariance levels off at the sill with the fitted curve through the points; the empirical points in fact overshoot the sill at middle distances and fall away at the largest ones.

Ordinary least squares assumes the residuals are independent. Field data rarely cooperates: plots close together tend to share unmeasured conditions (soil, drainage, microclimate), so their residuals look alike. The post on Moran’s I showed how to detect that dependence. This post is the remedy. Generalised least squares lets the residuals carry a spatial correlation structure, which fixes the problem ordinary regression hides: standard errors far smaller than the real spread of the estimate.

We will build a synthetic landscape so the truth is known, fit it both ways, and compare what each fit claims to know.

Setup: packages, colours and plot theme

library(nlme)
library(ggplot2)

paper <- "#f5f4ee"; ink <- "#16241d"; body <- "#2c3a31"
forest <- "#275139"; label <- "#46604a"; faint <- "#5d6b61"; line <- "#dad9ca"
base_theme <- theme_minimal(base_size = 13) +
  theme(plot.background = element_rect(fill = paper, colour = NA),
        panel.background = element_rect(fill = paper, colour = NA),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_line(colour = line, linewidth = 0.3),
        axis.title = element_text(colour = body), axis.text = element_text(colour = faint),
        legend.title = element_text(colour = body, size = 10),
        legend.text = element_text(colour = faint, size = 9))

A landscape with correlated residuals

The helper draws a Gaussian random field: correlated noise over space, built from an exponential covariance and a Cholesky factor. The covariate (canopy openness) is itself spatially structured, as real environmental gradients are. The response (understory cover) has a genuine canopy effect with a slope of 0.5 and an intercept of 10, plus spatially autocorrelated residual noise.

gp_field <- function(coords, range, nugget = 0) {
  d <- as.matrix(dist(coords))
  S <- exp(-d / range)
  diag(S) <- 1 + nugget
  as.vector(t(chol(S)) %*% rnorm(nrow(coords)))
}

set.seed(204)
n  <- 150
co <- data.frame(xc = runif(n, 0, 10), yc = runif(n, 0, 10))
canopy <- 50 + 15 * scale(gp_field(co, range = 3))[, 1]
cover  <- 10 + 0.5 * canopy + 8 * gp_field(co, range = 2.5, nugget = 0.15)
df <- data.frame(co, canopy = canopy, cover = cover)

The symptom

Fit ordinary least squares (gls with independent errors is identical to lm, and fitting it this way keeps every model on the same footing for later comparison) and map its residuals across the landscape.

m_ols <- gls(cover ~ canopy, data = df, method = "REML")
df$ols_resid <- residuals(m_ols, type = "response")
summary(m_ols)$tTable

                 Value  Std.Error   t-value      p-value
(Intercept) 24.5295648 1.97538230 12.417629 1.040608e-24
canopy       0.2733377 0.03785188  7.221244 2.511352e-11

ggplot(df, aes(xc, yc, colour = ols_resid)) +
  geom_point(size = 3.4, alpha = 0.95) +
  scale_colour_gradient2(low = "#b5534e", mid = "#e7e4d6", high = "#1d5b4e",
                         midpoint = 0, name = "OLS\nresidual") +
  coord_equal() +
  labs(x = "Easting", y = "Northing") +
  base_theme + theme(legend.position = "right")

Scatter of plot locations coloured by OLS residual. Positive residuals form a green cluster through the centre, negative residuals form red clusters at the lower left and right, showing spatial structure. — OLS residuals by plot location. Nearby plots share similar residuals: green patches sit above the fit, red patches below.

If the residuals were independent, this map would be a random wash of colours. Instead the greens and reds clump: positive residuals run through the centre, negative ones gather at the edges. A formal Moran’s I test (covered in the earlier post) would confirm it, but the eye is enough here. The independence assumption behind the OLS standard errors does not hold, so the p-values cannot be trusted.

The fix: a correlation structure on the residuals

gls from nlme fits the same linear model but lets you specify how residuals covary. corExp says the correlation between two plots decays exponentially with the distance between them; nugget = TRUE adds a share of variance that is not spatial, soaking up measurement error and fine-scale noise.

m_exp <- gls(cover ~ canopy, data = df, method = "REML",
             correlation = corExp(form = ~ xc + yc, nugget = TRUE))
summary(m_exp)$tTable

                 Value  Std.Error   t-value      p-value
(Intercept) 13.1905550 3.77229633  3.496691 6.221073e-04
canopy       0.4627138 0.04444089 10.411894 2.175866e-19

Compare the two fits against the values that generated the data (intercept 10, slope 0.5). OLS put the slope at 0.27 and the intercept at 24.5; GLS puts the slope at 0.46 and the intercept at 13.2. Resist the obvious reading. One landscape cannot tell a biased method from a noisy one, so repeat the simulation: two hundred fresh landscapes from the same recipe, each fitted both ways.

set.seed(205)
reps <- t(replicate(200, {
  co_r  <- data.frame(xc = runif(n, 0, 10), yc = runif(n, 0, 10))
  can_r <- 50 + 15 * scale(gp_field(co_r, range = 3))[, 1]
  cov_r <- 10 + 0.5 * can_r + 8 * gp_field(co_r, range = 2.5, nugget = 0.15)
  d_r   <- data.frame(co_r, canopy = can_r, cover = cov_r)
  o <- summary(lm(cover ~ canopy, data = d_r))$coefficients["canopy", 1:2]
  g <- tryCatch(   # an occasional fit fails to converge; it is dropped
    summary(gls(cover ~ canopy, data = d_r,
                correlation = corExp(form = ~ xc + yc, nugget = TRUE)))$tTable["canopy", 1:2],
    error = function(e) c(NA_real_, NA_real_))
  c(ols = o[[1]], ols_se = o[[2]], gls = g[[1]], gls_se = g[[2]])
}))
g_ok   <- !is.na(reps[, "gls"])
covers <- function(est, se) mean(abs(est - 0.5) < qt(0.975, n - 2) * se)
rep_sum <- c(ols_mean   = mean(reps[, "ols"]),
             ols_sd     = sd(reps[, "ols"]),
             ols_se     = mean(reps[, "ols_se"]),
             ols_cover  = covers(reps[, "ols"], reps[, "ols_se"]),
             gls_sd     = sd(reps[g_ok, "gls"]),
             gls_cover  = covers(reps[g_ok, "gls"], reps[g_ok, "gls_se"]),
             gls_failed = sum(!g_ok))
round(rep_sum, 3)

  ols_mean     ols_sd     ols_se  ols_cover     gls_sd  gls_cover gls_failed 
     0.496      0.117      0.041      0.505      0.056      0.934      2.000

OLS is not biased here: across the replicates its slope averages 0.50, the true 0.5, because the correlated noise is independent of canopy. What the correlation destroys is precision. The OLS slopes scatter with a standard deviation of 0.117, about 2.8 times the standard error a single OLS fit reports (0.041 on average), so the nominal 95 per cent interval covers the true slope in only 50 per cent of the landscapes. In the landscape above, the estimate landed at 0.27 while the fit insists it knows the slope to within 0.04. The reason is that OLS treated all 150 plots as independent evidence, when clusters of nearby plots are partly carrying the same information. Telling the model that, through the correlation structure, shrinks the effective sample size and widens the standard errors (the intercept’s nearly doubles, from 2.0 to 3.8). It also weights the plots by the fitted correlation instead of treating them as independent, which makes the GLS estimate more precise than the OLS one: across the replicates the GLS slopes scatter with a standard deviation of 0.056, and its intervals cover the true slope in 93 per cent of the landscapes. In the single landscape above, the GLS intervals cover the truth mainly because the estimates moved towards it, not because the intervals widened: the slope’s standard error grows only from 0.038 to 0.044, while the estimate moves from 0.27 to 0.46. Under spatial autocorrelation the OLS point estimate is fine on average, though it scatters more than the GLS one; the confidence attached to it is not.

The fitted structure has interpretable parts:

coef(m_exp$modelStruct$corStruct, unconstrained = FALSE)

    range    nugget 
2.3940610 0.1041068

The range, 2.4, is the distance over which correlation falls to about a third: plots closer than that are meaningfully linked, plots much farther apart are effectively independent. The nugget, 0.10, says about a tenth of the residual variance is non-spatial noise. A residual variogram shows the same thing: semivariance climbs steeply over the first few distance units. It does not settle neatly on the sill, though: the empirical points overshoot the fitted sill (the dotted line) at middle distances and fall away at the largest ones, so the fitted exponential is a smooth summary of the initial rise rather than a close fit to every point.

vg  <- Variogram(m_exp, resType = "response")
cs  <- coef(m_exp$modelStruct$corStruct, unconstrained = FALSE)
sig2 <- summary(m_exp)$sigma^2
curve_df <- data.frame(dist = seq(0.01, max(vg$dist), length.out = 200))
curve_df$gamma <- sig2 * (1 - (1 - cs[["nugget"]]) * exp(-curve_df$dist / cs[["range"]]))

ggplot() +
  geom_point(data = vg, aes(dist, variog), colour = forest, size = 2.4, alpha = 0.8) +
  geom_line(data = curve_df, aes(dist, gamma), colour = "#c9a13f", linewidth = 1) +
  geom_hline(yintercept = sig2, linetype = "dotted", colour = faint) +
  labs(x = "Distance between plots", y = "Semivariance of residuals") +
  base_theme

Semivariance of residuals against distance. Points rise from about 20 at the shortest distance to a peak near 75 to 78 at distances of about 4 to 5, well above the dotted sill near 60, then fall to about 30 by distance 10. The fitted gold exponential curve passes close to the first few points, sits below the points at middle distances and approaches the sill only near distance 10. — Empirical semivariogram of the residuals (points) with the fitted exponential model (line). Structure is clear out to the range; at large distances the empirical estimate is noisy: every point averages about the same number of plot pairs, since the distances are binned by quantiles, but the farthest pairs can only join plots towards opposite edges of the landscape, so they draw on fewer distinct plots.

Choosing the structure

corExp is one option among several. corGaus decays as a Gaussian (smoother fields), corSpher reaches independence at a finite distance. Fit the candidates and compare them, along with the independent-errors model, by AIC.

m_gau <- update(m_exp, correlation = corGaus(form = ~ xc + yc, nugget = TRUE))
m_sph <- update(m_exp, correlation = corSpher(form = ~ xc + yc, nugget = TRUE))
AIC(m_ols, m_exp, m_gau, m_sph)

      df       AIC
m_ols  3 1014.4773
m_exp  5  911.6535
m_gau  5  920.2874
m_sph  5  910.9282

Two things stand out. Any spatial structure beats the independent-errors model by about 100 AIC units: an overwhelming preference for modelling the autocorrelation at all. Among the spatial structures the gap is small. Exponential and spherical are within one AIC unit of each other (indistinguishable on this evidence), and Gaussian trails by about nine. The practical lesson is that the decision that matters is whether to model spatial correlation, not which decay shape you choose; pick the best by AIC, but do not agonise between exponential and spherical when they tie. Comparing correlation structures this way uses REML, which is the default; switch to method = "ML" only when you are comparing models with different fixed effects.

Why it is worth the trouble

The cost of ignoring autocorrelation is not abstract. Simulate the null case: a predictor with no real effect on the response, but both the predictor and the residuals spatially structured, fitted both ways one hundred and fifty times, counting how often each method declares the slope significant at 0.05. With a true slope of zero, every rejection is a false positive.

set.seed(771)
nsim <- 150; nn <- 100
cs2 <- data.frame(xc = runif(nn, 0, 10), yc = runif(nn, 0, 10))
dd  <- as.matrix(dist(cs2))
Lx <- t(chol({ S <- exp(-dd / 3);   diag(S) <- 1;        S }))
Le <- t(chol({ S <- exp(-dd / 2.5); diag(S) <- 1 + 0.15; S }))

res <- t(replicate(nsim, {
  xf <- 50 + 15 * scale(Lx %*% rnorm(nn))[, 1]
  yf <- 10 + 0 * xf + 8 * (Le %*% rnorm(nn))[, 1]   # true slope is zero
  p_ols <- summary(lm(yf ~ xf))$coefficients["xf", "Pr(>|t|)"]
  p_gls <- tryCatch(
    summary(gls(yf ~ xf, data = data.frame(cs2, xf, yf),
                correlation = corExp(form = ~ xc + yc, nugget = TRUE)))$tTable["xf", "p-value"],
    error = function(e) NA_real_)
  c(ols = p_ols, gls = p_gls)
}))
ok <- !is.na(res[, "gls"])
c(OLS = mean(res[ok, "ols"] < 0.05), GLS = mean(res[ok, "gls"] < 0.05))

       OLS        GLS 
0.33557047 0.08724832

td <- data.frame(method = factor(c("OLS", "GLS"), levels = c("OLS", "GLS")),
                 rate = c(mean(res[ok, "ols"] < 0.05), mean(res[ok, "gls"] < 0.05)))
ggplot(td, aes(method, rate, fill = method)) +
  geom_col(width = 0.55) +
  geom_hline(yintercept = 0.05, linetype = "dashed", colour = "#b5534e", linewidth = 0.6) +
  geom_text(aes(label = sprintf("%.1f%%", 100 * rate)), vjust = -0.5, colour = body, size = 4.4) +
  annotate("text", x = 1.5, y = 0.058, label = "nominal 5%",
           colour = "#b5534e", size = 3.3, hjust = 0.5) +
  scale_fill_manual(values = c(OLS = "#b5534e", GLS = "#275139")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.12)),
                     labels = function(x) paste0(x * 100, "%")) +
  labs(x = NULL, y = "False-positive rate (true slope = 0)") +
  base_theme + theme(legend.position = "none")

Two bars: OLS false-positive rate about 34 per cent, well above a dashed line at 5 per cent; GLS about 9 per cent, near the line. — False-positive rate when the slope is truly zero. OLS rejects far more than the nominal 5 per cent; GLS stays close to it.

Ordinary least squares calls a non-existent effect significant about a third of the time, nearly seven times the rate it advertises. GLS pulls that back close to 5 per cent. The handful of points by which it overshoots is the usual price of estimating the correlation parameters from the same data; it is a rounding error next to the OLS inflation.

A few practical notes. gls is for continuous, roughly Gaussian responses. For counts or presence-absence with spatial dependence, the analogue is a spatial term inside a mixed model or a spatial smooth inside a GAM. Start with corExp and a nugget, read the range to judge the scale of dependence, and always compare against the independent-errors model so you can see what the correction buys you.

References

Pinheiro and Bates 2000 Mixed-Effects Models in S and S-PLUS, Springer (10.1007/b98882).

Zuur, Ieno, Walker, Saveliev and Smith 2009 Mixed Effects Models and Extensions in Ecology with R, Springer (10.1007/978-0-387-87458-6).

Dormann et al. 2007 Ecography 30(5):609-628 (10.1111/j.2007.0906-7590.05171.x).

Legendre 1993 Ecology 74(6):1659-1673 (10.2307/1939924).

--- title: "Generalised least squares for spatial data in R" description: "When regression residuals are spatially autocorrelated, OLS standard errors mislead. Fit a GLS model in nlme, pick a correlation structure by AIC, and compare." date: "2026-05-06 14:00" date-modified: "2026-11-05" categories: [R, GLS, spatial, regression, ecology tutorial] image: thumbnail.png image-alt: "A landscape map of regression residuals where nearby plots share similar values, forming green and red clusters that reveal spatial autocorrelation." --- *Corrected 5 November 2026: the post said the spread of OLS slopes across replicate landscapes is close to three times the standard error a single fit reports, and that OLS scatters more than GLS, without computing either. A new chunk now repeats the simulation two hundred times and reports the ratio, the coverage of both sets of intervals and the GLS spread.* *Corrected 14 October 2026: the text said the GLS intervals cover the true values because their standard errors are wider; it now says that GLS also makes the estimate more precise and that here the intervals cover the truth mainly because the estimates moved towards it. The variogram caption no longer blames the noisy far points on a shortage of plot pairs, since every point averages about the same number of pairs; the farthest pairs draw on fewer distinct plots.* *Corrected 27 September 2026: the variogram text said the semivariance levels off at the sill with the fitted curve through the points; the empirical points in fact overshoot the sill at middle distances and fall away at the largest ones.* Ordinary least squares assumes the residuals are independent. Field data rarely cooperates: plots close together tend to share unmeasured conditions (soil, drainage, microclimate), so their residuals look alike. The [post on Moran's I](../spatial-autocorrelation-morans-i/) showed how to detect that dependence. This post is the remedy. Generalised least squares lets the residuals carry a spatial correlation structure, which fixes the problem ordinary regression hides: standard errors far smaller than the real spread of the estimate. We will build a synthetic landscape so the truth is known, fit it both ways, and compare what each fit claims to know. ```{r} #| code-fold: true #| code-summary: "Setup: packages, colours and plot theme" #| results: hide #| message: false #| warning: false library(nlme) library(ggplot2) paper <- "#f5f4ee"; ink <- "#16241d"; body <- "#2c3a31" forest <- "#275139"; label <- "#46604a"; faint <- "#5d6b61"; line <- "#dad9ca" base_theme <- theme_minimal(base_size = 13) + theme(plot.background = element_rect(fill = paper, colour = NA), panel.background = element_rect(fill = paper, colour = NA), panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = line, linewidth = 0.3), axis.title = element_text(colour = body), axis.text = element_text(colour = faint), legend.title = element_text(colour = body, size = 10), legend.text = element_text(colour = faint, size = 9)) ``` ## A landscape with correlated residuals The helper draws a Gaussian random field: correlated noise over space, built from an exponential covariance and a Cholesky factor. The covariate (canopy openness) is itself spatially structured, as real environmental gradients are. The response (understory cover) has a genuine canopy effect with a slope of 0.5 and an intercept of 10, plus spatially autocorrelated residual noise. ```{r} gp_field <- function(coords, range, nugget = 0) { d <- as.matrix(dist(coords)) S <- exp(-d / range) diag(S) <- 1 + nugget as.vector(t(chol(S)) %*% rnorm(nrow(coords))) } set.seed(204) n <- 150 co <- data.frame(xc = runif(n, 0, 10), yc = runif(n, 0, 10)) canopy <- 50 + 15 * scale(gp_field(co, range = 3))[, 1] cover <- 10 + 0.5 * canopy + 8 * gp_field(co, range = 2.5, nugget = 0.15) df <- data.frame(co, canopy = canopy, cover = cover) ``` ## The symptom Fit ordinary least squares (gls with independent errors is identical to `lm`, and fitting it this way keeps every model on the same footing for later comparison) and map its residuals across the landscape. ```{r} m_ols <- gls(cover ~ canopy, data = df, method = "REML") df$ols_resid <- residuals(m_ols, type = "response") summary(m_ols)$tTable ``` ```{r} #| fig-cap: "OLS residuals by plot location. Nearby plots share similar residuals: green patches sit above the fit, red patches below." #| fig-alt: "Scatter of plot locations coloured by OLS residual. Positive residuals form a green cluster through the centre, negative residuals form red clusters at the lower left and right, showing spatial structure." ggplot(df, aes(xc, yc, colour = ols_resid)) + geom_point(size = 3.4, alpha = 0.95) + scale_colour_gradient2(low = "#b5534e", mid = "#e7e4d6", high = "#1d5b4e", midpoint = 0, name = "OLS\nresidual") + coord_equal() + labs(x = "Easting", y = "Northing") + base_theme + theme(legend.position = "right") ``` If the residuals were independent, this map would be a random wash of colours. Instead the greens and reds clump: positive residuals run through the centre, negative ones gather at the edges. A formal Moran's I test (covered in the [earlier post](../spatial-autocorrelation-morans-i/)) would confirm it, but the eye is enough here. The independence assumption behind the OLS standard errors does not hold, so the p-values cannot be trusted. ## The fix: a correlation structure on the residuals `gls` from nlme fits the same linear model but lets you specify how residuals covary. `corExp` says the correlation between two plots decays exponentially with the distance between them; `nugget = TRUE` adds a share of variance that is not spatial, soaking up measurement error and fine-scale noise. ```{r} m_exp <- gls(cover ~ canopy, data = df, method = "REML", correlation = corExp(form = ~ xc + yc, nugget = TRUE)) summary(m_exp)$tTable ``` Compare the two fits against the values that generated the data (intercept 10, slope 0.5). OLS put the slope at `r sprintf("%.2f", summary(m_ols)$tTable["canopy", "Value"])` and the intercept at `r sprintf("%.1f", summary(m_ols)$tTable["(Intercept)", "Value"])`; GLS puts the slope at `r sprintf("%.2f", summary(m_exp)$tTable["canopy", "Value"])` and the intercept at `r sprintf("%.1f", summary(m_exp)$tTable["(Intercept)", "Value"])`. Resist the obvious reading. One landscape cannot tell a biased method from a noisy one, so repeat the simulation: two hundred fresh landscapes from the same recipe, each fitted both ways. ```{r} set.seed(205) reps <- t(replicate(200, { co_r <- data.frame(xc = runif(n, 0, 10), yc = runif(n, 0, 10)) can_r <- 50 + 15 * scale(gp_field(co_r, range = 3))[, 1] cov_r <- 10 + 0.5 * can_r + 8 * gp_field(co_r, range = 2.5, nugget = 0.15) d_r <- data.frame(co_r, canopy = can_r, cover = cov_r) o <- summary(lm(cover ~ canopy, data = d_r))$coefficients["canopy", 1:2] g <- tryCatch( # an occasional fit fails to converge; it is dropped summary(gls(cover ~ canopy, data = d_r, correlation = corExp(form = ~ xc + yc, nugget = TRUE)))$tTable["canopy", 1:2], error = function(e) c(NA_real_, NA_real_)) c(ols = o[[1]], ols_se = o[[2]], gls = g[[1]], gls_se = g[[2]]) })) g_ok <- !is.na(reps[, "gls"]) covers <- function(est, se) mean(abs(est - 0.5) < qt(0.975, n - 2) * se) rep_sum <- c(ols_mean = mean(reps[, "ols"]), ols_sd = sd(reps[, "ols"]), ols_se = mean(reps[, "ols_se"]), ols_cover = covers(reps[, "ols"], reps[, "ols_se"]), gls_sd = sd(reps[g_ok, "gls"]), gls_cover = covers(reps[g_ok, "gls"], reps[g_ok, "gls_se"]), gls_failed = sum(!g_ok)) round(rep_sum, 3) ``` OLS is not biased here: across the replicates its slope averages `r sprintf("%.2f", rep_sum[["ols_mean"]])`, the true 0.5, because the correlated noise is independent of canopy. What the correlation destroys is precision. The OLS slopes scatter with a standard deviation of `r sprintf("%.3f", rep_sum[["ols_sd"]])`, about `r sprintf("%.1f", rep_sum[["ols_sd"]] / rep_sum[["ols_se"]])` times the standard error a single OLS fit reports (`r sprintf("%.3f", rep_sum[["ols_se"]])` on average), so the nominal 95 per cent interval covers the true slope in only `r sprintf("%.0f", 100 * rep_sum[["ols_cover"]])` per cent of the landscapes. In the landscape above, the estimate landed at `r sprintf("%.2f", summary(m_ols)$tTable["canopy", "Value"])` while the fit insists it knows the slope to within `r sprintf("%.2f", summary(m_ols)$tTable["canopy", "Std.Error"])`. The reason is that OLS treated all 150 plots as independent evidence, when clusters of nearby plots are partly carrying the same information. Telling the model that, through the correlation structure, shrinks the effective sample size and widens the standard errors (the intercept's nearly doubles, from `r sprintf("%.1f", summary(m_ols)$tTable["(Intercept)", "Std.Error"])` to `r sprintf("%.1f", summary(m_exp)$tTable["(Intercept)", "Std.Error"])`). It also weights the plots by the fitted correlation instead of treating them as independent, which makes the GLS estimate more precise than the OLS one: across the replicates the GLS slopes scatter with a standard deviation of `r sprintf("%.3f", rep_sum[["gls_sd"]])`, and its intervals cover the true slope in `r sprintf("%.0f", 100 * rep_sum[["gls_cover"]])` per cent of the landscapes. In the single landscape above, the GLS intervals cover the truth mainly because the estimates moved towards it, not because the intervals widened: the slope's standard error grows only from `r sprintf("%.3f", summary(m_ols)$tTable["canopy", "Std.Error"])` to `r sprintf("%.3f", summary(m_exp)$tTable["canopy", "Std.Error"])`, while the estimate moves from `r sprintf("%.2f", summary(m_ols)$tTable["canopy", "Value"])` to `r sprintf("%.2f", summary(m_exp)$tTable["canopy", "Value"])`. Under spatial autocorrelation the OLS point estimate is fine on average, though it scatters more than the GLS one; the confidence attached to it is not. The fitted structure has interpretable parts: ```{r} coef(m_exp$modelStruct$corStruct, unconstrained = FALSE) ``` The range, `r sprintf("%.1f", coef(m_exp$modelStruct$corStruct, unconstrained = FALSE)[["range"]])`, is the distance over which correlation falls to about a third: plots closer than that are meaningfully linked, plots much farther apart are effectively independent. The nugget, `r sprintf("%.2f", coef(m_exp$modelStruct$corStruct, unconstrained = FALSE)[["nugget"]])`, says about a tenth of the residual variance is non-spatial noise. A residual variogram shows the same thing: semivariance climbs steeply over the first few distance units. It does not settle neatly on the sill, though: the empirical points overshoot the fitted sill (the dotted line) at middle distances and fall away at the largest ones, so the fitted exponential is a smooth summary of the initial rise rather than a close fit to every point. ```{r} #| fig-cap: "Empirical semivariogram of the residuals (points) with the fitted exponential model (line). Structure is clear out to the range; at large distances the empirical estimate is noisy: every point averages about the same number of plot pairs, since the distances are binned by quantiles, but the farthest pairs can only join plots towards opposite edges of the landscape, so they draw on fewer distinct plots." #| fig-alt: "Semivariance of residuals against distance. Points rise from about 20 at the shortest distance to a peak near 75 to 78 at distances of about 4 to 5, well above the dotted sill near 60, then fall to about 30 by distance 10. The fitted gold exponential curve passes close to the first few points, sits below the points at middle distances and approaches the sill only near distance 10." vg <- Variogram(m_exp, resType = "response") cs <- coef(m_exp$modelStruct$corStruct, unconstrained = FALSE) sig2 <- summary(m_exp)$sigma^2 curve_df <- data.frame(dist = seq(0.01, max(vg$dist), length.out = 200)) curve_df$gamma <- sig2 * (1 - (1 - cs[["nugget"]]) * exp(-curve_df$dist / cs[["range"]])) ggplot() + geom_point(data = vg, aes(dist, variog), colour = forest, size = 2.4, alpha = 0.8) + geom_line(data = curve_df, aes(dist, gamma), colour = "#c9a13f", linewidth = 1) + geom_hline(yintercept = sig2, linetype = "dotted", colour = faint) + labs(x = "Distance between plots", y = "Semivariance of residuals") + base_theme ``` ## Choosing the structure `corExp` is one option among several. `corGaus` decays as a Gaussian (smoother fields), `corSpher` reaches independence at a finite distance. Fit the candidates and compare them, along with the independent-errors model, by AIC. ```{r} m_gau <- update(m_exp, correlation = corGaus(form = ~ xc + yc, nugget = TRUE)) m_sph <- update(m_exp, correlation = corSpher(form = ~ xc + yc, nugget = TRUE)) AIC(m_ols, m_exp, m_gau, m_sph) ``` Two things stand out. Any spatial structure beats the independent-errors model by about 100 AIC units: an overwhelming preference for modelling the autocorrelation at all. Among the spatial structures the gap is small. Exponential and spherical are within one AIC unit of each other (indistinguishable on this evidence), and Gaussian trails by about nine. The practical lesson is that the decision that matters is whether to model spatial correlation, not which decay shape you choose; pick the best by AIC, but do not agonise between exponential and spherical when they tie. Comparing correlation structures this way uses REML, which is the default; switch to `method = "ML"` only when you are comparing models with different fixed effects. ## Why it is worth the trouble The cost of ignoring autocorrelation is not abstract. Simulate the null case: a predictor with no real effect on the response, but both the predictor and the residuals spatially structured, fitted both ways one hundred and fifty times, counting how often each method declares the slope significant at 0.05. With a true slope of zero, every rejection is a false positive. ```{r} set.seed(771) nsim <- 150; nn <- 100 cs2 <- data.frame(xc = runif(nn, 0, 10), yc = runif(nn, 0, 10)) dd <- as.matrix(dist(cs2)) Lx <- t(chol({ S <- exp(-dd / 3); diag(S) <- 1; S })) Le <- t(chol({ S <- exp(-dd / 2.5); diag(S) <- 1 + 0.15; S })) res <- t(replicate(nsim, { xf <- 50 + 15 * scale(Lx %*% rnorm(nn))[, 1] yf <- 10 + 0 * xf + 8 * (Le %*% rnorm(nn))[, 1] # true slope is zero p_ols <- summary(lm(yf ~ xf))$coefficients["xf", "Pr(>|t|)"] p_gls <- tryCatch( summary(gls(yf ~ xf, data = data.frame(cs2, xf, yf), correlation = corExp(form = ~ xc + yc, nugget = TRUE)))$tTable["xf", "p-value"], error = function(e) NA_real_) c(ols = p_ols, gls = p_gls) })) ok <- !is.na(res[, "gls"]) c(OLS = mean(res[ok, "ols"] < 0.05), GLS = mean(res[ok, "gls"] < 0.05)) ``` ```{r} #| fig-cap: "False-positive rate when the slope is truly zero. OLS rejects far more than the nominal 5 per cent; GLS stays close to it." #| fig-alt: "Two bars: OLS false-positive rate about 34 per cent, well above a dashed line at 5 per cent; GLS about 9 per cent, near the line." td <- data.frame(method = factor(c("OLS", "GLS"), levels = c("OLS", "GLS")), rate = c(mean(res[ok, "ols"] < 0.05), mean(res[ok, "gls"] < 0.05))) ggplot(td, aes(method, rate, fill = method)) + geom_col(width = 0.55) + geom_hline(yintercept = 0.05, linetype = "dashed", colour = "#b5534e", linewidth = 0.6) + geom_text(aes(label = sprintf("%.1f%%", 100 * rate)), vjust = -0.5, colour = body, size = 4.4) + annotate("text", x = 1.5, y = 0.058, label = "nominal 5%", colour = "#b5534e", size = 3.3, hjust = 0.5) + scale_fill_manual(values = c(OLS = "#b5534e", GLS = "#275139")) + scale_y_continuous(expand = expansion(mult = c(0, 0.12)), labels = function(x) paste0(x * 100, "%")) + labs(x = NULL, y = "False-positive rate (true slope = 0)") + base_theme + theme(legend.position = "none") ``` Ordinary least squares calls a non-existent effect significant about a third of the time, nearly seven times the rate it advertises. GLS pulls that back close to 5 per cent. The handful of points by which it overshoots is the usual price of estimating the correlation parameters from the same data; it is a rounding error next to the OLS inflation. A few practical notes. `gls` is for continuous, roughly Gaussian responses. For counts or presence-absence with spatial dependence, the analogue is a spatial term inside a [mixed model](../glmm-nested-counts-pseudoreplication/) or a spatial smooth inside a [GAM](../gam-species-response-curves/). Start with `corExp` and a nugget, read the range to judge the scale of dependence, and always compare against the independent-errors model so you can see what the correction buys you. ## References Pinheiro and Bates 2000 Mixed-Effects Models in S and S-PLUS, Springer (10.1007/b98882). Zuur, Ieno, Walker, Saveliev and Smith 2009 Mixed Effects Models and Extensions in Ecology with R, Springer (10.1007/978-0-387-87458-6). Dormann et al. 2007 Ecography 30(5):609-628 (10.1111/j.2007.0906-7590.05171.x). Legendre 1993 Ecology 74(6):1659-1673 (10.2307/1939924). ## Related tutorials - [Separable space-time covariance, tested](../separable-space-time-covariance/) - [Stream networks and tail-up covariance](../stream-networks-and-tail-up-covariance/) - [Conditional autoregressive (CAR) models](../conditional-autoregressive-car-models/) - [The variogram and the GP kernel](../the-variogram-and-the-gp-kernel/)