ggeffects::ggpredict处理nlme::lme模型无法返回总体水平预测区间

分析存在自相关性的重复测量数据时，需使用nlme包的lme()函数拟合带有correlation = corAR1()自相关结构的线性混合效应模型，调用ggeffects包的ggpredict()函数、设置type = "re"参数获取总体水平预测区间（PI）时返回结果不符合预期，但等价设定的lme4::lmer()模型可正常返回正确结果，具体异常如下：

• 对无自相关结构的基础nlme::lme()模型调用ggpredict()时，结果仅返回自变量x的取值，完全缺失预测值与95%预测区间，即便显式在terms参数中指定"Subject [0]"也无改善
• 为lme()模型添加corAR1自相关结构后，ggpredict()可返回预测值与区间，但结果仅针对Subject因子的第一个水平（Subject=3），无法得到总体水平（Subject=0）的预测结果

library(lme4)
library(nlme)
library(ggeffects)

Data <- data.frame(
  Subject = factor(c(3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 
                     4, 4, 4, 4, 4, 4, 4, 4, 4, 
                     5, 5, 5, 5, 5, 5, 5, 5, 5, 
                     6, 6, 6, 6, 6, 6, 6, 6, 
                     7, 7, 7, 7, 7, 7, 7, 7, 
                     8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 
                     9, 9, 9, 9, 9, 9, 9, 9, 
                     13, 13, 13, 13, 13, 13, 13, 13, 
                     14, 14, 14, 14, 14, 14, 14, 14, 14, 
                     19, 19, 19, 19, 19, 19, 19)), 
  x = c(20.0, 28.5, 38.0, 47.5, 57.0, 66.5, 76.0, 85.5, 95.0, 100.0, 
           21.0, 31.5, 42.0, 53.0, 63.0, 73.5, 84.0, 95.0, 100.0, 
           20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 
           22.0, 33.0, 44.0, 56.0, 67.0, 78.0, 89.0, 100.0, 
           21.5, 32.0, 43.0, 54.0, 65.0, 76.0, 86.5, 100.0, 
           20.0, 29.0, 38.5, 48.5, 58.0, 67.5, 77.0, 87.0, 96.5, 100.0, 
           23.0, 33.0, 44.0, 56.0, 67.0, 78.0, 89.0, 100.0, 
           23.5, 34.5, 46.5, 57.5, 69.5, 80.5, 92.5, 100.0, 
           20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 
           25.0, 37.5, 50.0, 62.5, 75.0, 87.5, 100.0)/100,
  y = c(1.10, 1.00, 1.25, 1.60, 1.40, 1.20, 2.50, 4.60, 6.80, 10.40, 
          0.90, 1.00, 0.75, 0.90, 1.10, 1.70, 4.35, 9.95, 11.45, 
          1.20, 0.70, 1.30, 1.40, 0.70, 1.25, 2.30, 4.30, 8.20, 
          1.55, 1.15, 0.95, 1.10, 1.90, 3.25, 7.20, 14.30, 
          1.85, 2.00, 1.70, 2.00, 2.35, 3.30, 7.30, 12.10, 
          2.20, 1.95, 1.15, 1.55, 1.65, 3.00, 4.45, 9.05, 13.75, 15.85, 
          1.55, 1.20, 1.35, 1.60, 1.65, 4.70, 6.45, 10.80, 
          1.00, 0.90, 1.00, 1.10, 1.60, 3.60, 8.05, 12.30, 
          0.85, 1.00, 1.05, 1.00, 1.35, 2.00, 3.65, 6.75, 13.10, 
          2.25, 2.35, 2.40, 2.80, 4.90, 8.15, 13.50)
)

Model.lme4 <- lmer(
  y ~ x + (1 | Subject),
  data = Data
)

# 拟合不带自相关结构的lme模型
Model.nlme <- lme(
  fixed = y ~ x,
  random = ~ 1 | Subject, 
  data = Data,
)

# lmer()模型可正常返回预期结果
ggpredict(
  Model.lme4,
  terms = c("x [all]"),
  type = "re",
)

# 对无自相关结构的lme()模型调用时，预测值和预测区间缺失
ggpredict(
  Model.nlme,
  terms = c("x [all]"),
  type = "re",
)

添加correlation = corAR1()参数后，无预测值的问题仍然存在；显式调用terms = c("x [all]", "Subject [0]")也无法解决。当添加目标自相关结构后，可返回预测值与预测区间，但结果仅对应Subject因子的第一个水平：

Model.nlme <- lme(
  fixed = y ~ x,
  random = ~ 1 | Subject, 
  correlation = corAR1(form = ~ x | Subject),
  data = Data,
)

ggpredict(
  Model.nlme,
  terms = c("x [all]"),
  type = "re",
)

该问题不属于代码编写错误，是ggeffects包对带相关结构的nlme::lme模型的type = "re"参数适配不完善导致：

• 无自相关结构的基础lme模型返回空预测值，是旧版本ggeffects对nlme对象的方差分量提取逻辑存在bug，无法正确识别随机效应方差时就会返回空预测结果
• 添加corAR1结构后默认返回第一个受试者的结果，是因为模型存在自相关结构时，ggpredict()默认将分组变量作为条件变量取参考水平，没有正确识别type = "re"对应的总体平均（随机效应设为0）的计算逻辑

方案1：更新ggeffects并使用predict_response()接口 #

新版本ggeffects已优化对nlme模型的支持，推荐使用官方推荐的predict_response()作为预测接口，显式指定总体水平参数：

# 先将ggeffects更新到最新版本
# install.packages("ggeffects")
library(ggeffects)

predict_response(
  Model.nlme,
  terms = "x [all]",
  margin = "empirical",
  type = "re",
  condition = c(Subject = 0)
)

方案2：手动调用nlme原生预测函数计算 #

如果包函数适配仍有问题，可以直接调用predict.lme()手动计算总体水平预测区间，逻辑完全可控：

# 构建预测数据集
newdat <- data.frame(x = sort(unique(Data$x)))

# level=0对应总体水平预测，不叠加个体随机效应
pred <- predict(
  Model.nlme,
  newdata = newdat,
  level = 0,
  se.fit = TRUE
)

# 计算预测区间标准误：固定效应不确定性+随机效应方差+残差方差
se_pred <- sqrt(
  pred$se.fit^2 + 
    as.numeric(VarCorr(Model.nlme)[1, "Variance"]) + 
    Model.nlme$sigma^2
)

# 组装结果
result <- data.frame(
  x = newdat$x,
  Predicted = pred$fit,
  CI_low = pred$fit - 1.96 * se_pred,
  CI_high = pred$fit + 1.96 * se_pred
)

方案3：改用glmmTMB包拟合带AR1结构的模型 #

glmmTMB包支持AR1自相关结构的线性混合模型拟合，且ggeffects对glmmTMB对象的适配成熟度远高于带相关结构的nlme对象：

library(glmmTMB)
# 拟合带AR1结构的混合模型
Model.glmm <- glmmTMB(
  y ~ x + ar1(x + 0 | Subject),
  data = Data
)
# 直接调用ggpredict即可正常返回总体水平预测区间
ggpredict(Model.glmm, terms = "x [all]", type = "re")

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