Gamma广义线性混合模型(GLMM)中空间自相关的处理问询

我正在使用一个研究植物生长与多种因素（包括14个采样点的海冰范围SeaIce）关系的数据集，采用带log连接函数的Gamma广义线性混合模型(GLMM)进行分析，但担忧模型残差存在空间自相关，不确定最优解决方法。

我使用DHARMa包和lme4包进行模型拟合：

library(DHARMa)
library(lme4)
library(MASS)
library(gstat)
library(dplyr)
library(sf)
library(sp)

数据集可在GitHub仓库获取：

FinalDataset <- read.csv("https://raw.githubusercontent.com/derek-corcoran-barrios/SeaIceQuestion/master/FinalDataset.csv")

当前模型结构 #

响应变量 #

• Growth： 代表每株植物的年生长增量。

预测变量 #

• SeaIce.s： 每个采样点每年的海冰范围（已标准化），为核心关注预测变量。
• age： 每株植物各测量年份的年龄，因植物年龄与生长的生物学关系，与SeaIce.s纳入交互项。

随机效应 #

• Site： 代表14个间距不同的采样点的因子，纳入随机斜率(SeaIce.s | Site)以解释不同采样点生长-海冰关系的斜率差异。
• year： 研究覆盖1983-2015年，作为随机效应(1 | year)纳入。
• Individual： 代表每株被测植物的因子。

模型代码 #

fullmod <- glmer(Growth ~ SeaIce.s * age + (SeaIce.s | Site) + (1 | year) + (1 | Individual), data = FinalDataset, family = Gamma(link = "log"), control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5)))

我的核心问题是如何在模型中合理处理空间自相关。通过以下DHARMa代码检测到显著问题：

res2_null <- simulateResiduals(fullmod)
res3_null <- recalculateResiduals(res2_null, group = FinalDataset$Site)
locs_null <- FinalDataset %>% group_by(Site) %>% summarise(across(c(Latitude, Longitude), mean))
testSpatialAutocorrelation(res3_null, x = locs_null$Longitude, y = locs_null$Latitude)

## 
##  DHARMa Moran's I test for distance-based autocorrelation
## 
## data:  res3_null
## observed = 0.343448, expected = -0.076923, sd = 0.150983, p-value =
## 0.005365
## alternative hypothesis: Distance-based autocorrelation

结果显示存在显著的空间自相关。

已尝试的解决方案 #

尝试方案1 #

参考相关问题，尝试将经纬度作为固定效应纳入模型，但使用DHARMa的testSpatialAutocorrelation()检测仍得到显著的Moran's I值：

LonLatmod <- glmer(Growth ~ SeaIce.s * age + Latitude + Longitude + (SeaIce.s | Site) + (1 | year) + (1 | Individual), data = FinalDataset, family = Gamma(link = "log"), control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5)))

res <-  simulateResiduals(LonLatmod)
res2 <- recalculateResiduals(res, group = FinalDataset$Site)
locs_latlon <- FinalDataset %>% group_by(Site) %>% summarise(across(c(Latitude, Longitude), mean))
testSpatialAutocorrelation(res2, x =  locs_latlon$Longitude, y = locs_latlon$Latitude)

## 
##  DHARMa Moran's I test for distance-based autocorrelation
## 
## data:  res2
## observed = 0.382948, expected = -0.076923, sd = 0.150649, p-value =
## 0.002269
## alternative hypothesis: Distance-based autocorrelation

尝试方案2 #

尝试使用MASS包的glmmPQL()，但即使未添加相关结构也报错：

formula_glmmPQL <- as.formula("Growth ~ SeaIce.s * age")
model_glmmPQL <- glmmPQL(formula_glmmPQL,
                         random = list(~ SeaIce.s|Site, 
                                      ~ 1|year, 
                                      ~1| Individual),
                         data = FinalDataset,
                         na.action=na.omit,
                         family = Gamma(link = "log"))

#Error in pdFactor.pdLogChol(X[[i]], ...) :
#NA/NaN/Inf in foreign function call (arg 3)

尝试方案3 #

尝试使用gstat包的variogram()，但不确定如何解释变异函数形状及如何调整模型：

FinalDatasetSF <- st_as_sf(FinalDataset, coords = c("Longitude", "Latitude"), crs = st_crs(4326))

null_mod <- variogram(log(Growth) ~ 1, FinalDatasetSF)

Abn_fit_null <- fit.variogram(null_mod, model = vgm(1, "Sph", 
                                              700, 1))
plot(null_mod, model=Abn_fit_null)

我还尝试了nlme包的corAR1()函数，但发现它与glmer()模型不兼容。

我尝试了@SarahS提出的解决方案，分别添加和不添加经纬度作为固定效应：

require(glmmTMB)
# Set up the necessary variables
FinalDataset$pos <- numFactor(FinalDataset$Latitude, FinalDataset$Longitude)
FinalDataset$group <- factor(rep(1, nrow(FinalDataset)))
# Fit model
TestA <- glmmTMB(Growth ~ SeaIce.s * age + Latitude + Longitude + (SeaIce.s | Site) + (1 | year) + (1 | Individual) + 1 + exp(pos + 0 | group), data = FinalDataset, family = Gamma(link = "log"))


TestB <- glmmTMB(Growth ~ SeaIce.s * age + Latitude + Longitude + (SeaIce.s | Site) + (1 | year) + (1 | Individual) + 1 + exp(pos + 0 | group), data = FinalDataset, family = Gamma(link = "log"))

但两种模型的空间自相关问题仍未解决：

res <-  simulateResiduals(TestA)
res2 <- recalculateResiduals(res, group = FinalDataset$Site)

locs_latlon <- FinalDataset %>% group_by(Site) %>% summarise(across(c(Latitude, Longitude), mean))

testSpatialAutocorrelation(res2, x =  locs_latlon$Longitude, y = locs_latlon$Latitude)

## 
##  DHARMa Moran's I test for distance-based autocorrelation
## 
## data:  res2
## observed = 0.409445, expected = -0.076923, sd = 0.151346, p-value =
## 0.001311
## alternative hypothesis: Distance-based autocorrelation

以及：

res <-  simulateResiduals(TestB)
res2 <- recalculateResiduals(res, group = FinalDataset$Site)

testSpatialAutocorrelation(res2, x =  locs_latlon$Longitude, y = locs_latlon$Latitude)

## 
##  DHARMa Moran's I test for distance-based autocorrelation
## 
## data:  res2
## observed = 0.409445, expected = -0.076923, sd = 0.151346, p-value =
## 0.001311
## alternative hypothesis: Distance-based autocorrelation

请求建议 #

恳请针对当前建模方法，提供有效处理空间自相关的建议与见解。

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