问题概述

我正在使用一个数据集，该数据集检查植物生长与各种因素的关系，包括14个不同采样点的海冰范围（SeaIce）。为了分析数据，我使用了带有对数链接函数的Gamma广义线性混合模型（GLMM）。然而，我担心模型残差中的空间自相关，我不确定解决这个问题的最佳方法。
为了解决这个问题，我使用DHARMa包沿着lme4进行模型拟合。

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

您可以在此GitHub repository中找到此模型的数据集。

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

当前模型结构

响应变量

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

预测器

***SeaIce.s：**每个地点每年的海冰范围，主要的预测因素（比例）。
***年龄：**每个测量年份的每株植物的年龄，由于植物年龄和生长之间的生物学关系，作为与海冰的相互作用包括在内。

随机效果

***地点：**代表14个采样地点的系数，这些地点之间的距离不同。我包括一个随机斜率(SeaIce.s | Site)，以说明每个站点的生长-海冰关系中的不同斜率。
***年份：**本研究涵盖1983年至2015年，并作为随机效应(1 | year)纳入。
***个体：**代表测量的每株植物个体的因子。

型号

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

我在this question之后尝试的解决方案之一是将纬度和经度作为固定效果。然而，当我使用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包中的corAR 1（）函数，但它似乎与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

征求建议

我很想知道如何在我目前的建模方法中有效地解决空间自相关问题。