非nlme包广义混合效应模型的半变异函数绘制方法咨询

报错原因 #

nlme包的Variogram()为S3方法，仅适配nlme包原生的lme、gls类模型，未兼容MASS包的glmmPQL和lme4包的glmer模型对象，因此直接调用会出现坐标找不到、距离未指定等不兼容报错。

替代实现方案 #

可以手动提取模型残差，结合gstat包完成半变异函数的计算和绘制，全流程可运行示例代码如下：

步骤1：数据预处理与模型拟合 #

# 加载基础包
library(lme4)
library(MASS)
library(gstat)
library(sp)
data(quakes)

# 修正二分类变量生成逻辑，数值比较不要加引号
quakes$thresh <- ifelse(quakes$mag > 5, 0, 1)

# 拟合glmer二分类模型
model3 <- glmer(as.factor(thresh) ~ depth + (1|stations), data = quakes, family = binomial)

# 拟合glmmPQL二分类模型（可选）
model2 <- glmmPQL(as.factor(thresh) ~ depth, random = ~1|stations, family = binomial, data = quakes)

步骤2：提取残差并计算半变异函数 #

# --- 以glmer模型为例 ---
# 提取标准化皮尔逊残差，匹配原需求normalized类残差
res_df <- data.frame(
  res = residuals(model3, type = "pearson"),
  long = quakes$long,
  lat = quakes$lat
)
# 转换为空间点数据格式
coordinates(res_df) <- ~ long + lat
# 计算半变异函数
semivario <- variogram(res ~ 1, data = res_df)

# --- 如果是glmmPQL模型，仅替换残差提取行即可 ---
# res_df <- data.frame(
#   res = residuals(model2, type = "pearson"),
#   long = quakes$long,
#   lat = quakes$lat
# )

步骤3：绘制半变异图 #

# 绘制半变异函数，支持加平滑线
plot(semivario, pch = 16, col = "steelblue", smooth = TRUE, main = "GLMM残差半变异函数")

如果需要拟合理论半变异模型，可调用gstat包的fit.variogram()函数，输出结构和nlme的Variogram高度一致。

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