含时间序列数据的混合效应模型空间自相关处理及lme报错解决

Hey there! Let's tackle this problem step by step—you're already ahead of the game by identifying spatial autocorrelation with variograms and Moran's I, so great job on that front. Below's a breakdown of how to fix your lme error and effectively address spatial autocorrelation in your spatiotemporal dataset:

First: Fixing the lme Model Error #

The cannot have zero distances in "corSpatial" error stems directly from your data structure: multiple observations (different years) share identical (x,y) coordinates for the same region. This creates zero-distance pairs in the spatial distance matrix, which corSpatial functions can't handle (it leads to singular covariance matrix issues). Here are 3 actionable fixes:

1. Use a Combined Spatiotemporal Correlation Structure #

Since you have both spatial and temporal data, incorporate time into the correlation structure to avoid zero-distance pairs. This lets the model capture both spatial and temporal autocorrelation:

mod.cor <- lme(FLkm ~ Monsoon.Precip + Monsoon.Temp,
               correlation = corGaus(form = ~ x + y + Year | State, nugget = TRUE),
               data = NE1,
               random = ~1|State)

The ~x + y + Year | State formula means the model calculates correlation based on both spatial coordinates and year, within each state. This ensures no two observations have identical "distance" values (even for the same region, years differ).

2. Add Tiny Noise to Spatial Coordinates #

If you only want to model spatial autocorrelation, add a minuscule amount of random noise to your (x,y) coordinates to eliminate exact overlaps. Make sure the noise scale is far smaller than your actual spatial resolution so it doesn't bias results:

# Add tiny jitter to coordinates (1e-6 is negligible for most spatial datasets)
NE1$x_jitter <- NE1$x + rnorm(nrow(NE1), mean = 0, sd = 1e-6)
NE1$y_jitter <- NE1$y + rnorm(nrow(NE1), mean = 0, sd = 1e-6)

# Refit the model with jittered coordinates
mod.cor <- lme(FLkm ~ Monsoon.Precip + Monsoon.Temp,
               correlation = corGaus(form = ~ x_jitter + y_jitter, nugget = TRUE),
               data = NE1,
               random = ~1|State)

3. Switch to gls with Combined Spatial-Temporal Correlation #

If the state-level random effect isn't critical, use gls (Generalized Least Squares) which offers more flexible correlation structure combinations. For example, pair spatial correlation with temporal autocorrelation:

library(nlme)
mod.gls <- gls(FLkm ~ Monsoon.Precip + Monsoon.Temp,
               correlation = corAR1(form = ~1|RegionID) * corGaus(form = ~x+y, nugget = TRUE),
               data = NE1)

Here corAR1 models first-order temporal autocorrelation within each region, while corGaus handles spatial autocorrelation.

Other Methods to Account for Spatial Autocorrelation #

Beyond fixing the lme model, the approaches you mentioned can be adapted to your spatiotemporal data with these details:

1. OLS + Spatial Lag/Error Models (SAR/SEM) #

OLS ignores spatial autocorrelation, but you can extend it to spatial regression models while controlling for time effects with year fixed variables:

library(spdep)
# Build a spatial weight matrix (k-nearest neighbors example)
coords <- NE1[,c("x","y")]
knn <- knearneigh(coords, k = 5)
w_nb <- knn2nb(knn)
w_mat <- nb2listw(w_nb, style = "W")

# Add year fixed effects and fit a Spatial Lag Model (SAR)
NE1$Year <- as.factor(NE1$Year)
mod.sar <- lagsarlm(FLkm ~ Monsoon.Precip + Monsoon.Temp + Year,
                    data = NE1,
                    listw = w_mat)

# Or fit a Spatial Error Model (SEM) if autocorrelation is in residuals
mod.sem <- errorsarlm(FLkm ~ Monsoon.Precip + Monsoon.Temp + Year,
                      data = NE1,
                      listw = w_mat)

2. Spatiotemporal Geographically Weighted Regression (GTWR) #

Standard GWR works for cross-sectional data, but GTWR extends it to account for both spatial and temporal heterogeneity:

library(GTWR)
# Convert year to a numeric time variable (e.g., 2000 = 0, 2018 = 18)
NE1$Time <- NE1$Year - 2000
# Create spatiotemporal coordinates
coords_t <- cbind(NE1$x, NE1$y, NE1$Time)

# Fit GTWR (use cross-validation to select optimal bandwidth)
mod.gtwr <- GTWR(formula = FLkm ~ Monsoon.Precip + Monsoon.Temp,
                 data = NE1,
                 coords = coords_t,
                 bw = 80, # Adjust this with cv.gtwr()
                 adaptive = TRUE)

GTWR lets regression coefficients vary across both space and time, making it ideal for capturing local spatiotemporal patterns.

3. Panel Spatial Autoregression (PSAR) #

Since you have a panel dataset (590 regions × 19 years), PSAR combines panel model strengths (controlling for individual/time effects) with spatial autocorrelation handling:

library(plm)
library(spdep)
# Convert data to panel format (replace "RegionID" with your actual region identifier)
pdata <- pdata.frame(NE1, index = c("RegionID", "Year"))

# Fit a panel spatial lag model with random effects
mod.psar <- spml(FLkm ~ Monsoon.Precip + Monsoon.Temp,
                 data = pdata,
                 listw = w_mat,
                 model = "random",
                 effect = "individual")

Final Notes #

Choose the method based on your research goals:

• Stick with lme if you need to retain state-level random effects.
• Use GTWR if you want to explore local spatiotemporal variation in coefficients.
• Opt for SAR/SEM or PSAR if your priority is correcting for spatial autocorrelation in model estimates.

内容的提问来源于stack exchange，提问作者Alice