R中lme与lmer拟合等价混合模型残差不一致问题排查

我分别使用R语言中两款常用包的函数拟合混合效应模型：

• nlme包的lme函数
• lme4包的lmer函数
  为在lme4中复现和lme参数设定完全一致的异方差混合效应模型，我使用了已知的hack技巧，通过构造特殊随机效应项在lmer中实现异方差结构。
  由于未找到匹配本问题场景的R内置数据集，示例数据单独存储提供。

library(nlme)
library(lme4)
ModLME = lme(Var1~I(Var2)+I(Var2^2),
                         random = ~1|Var3,
                         weights = varIdent(form=~1|Var4),
                         Dataone, method="REML")

ModLMER = lmer(Var1~I(Var2)+I(Var2^2)+(1|Var3)+(0+dummy(Var4,"1")|Var5),
                Dataone, REML = TRUE,
                control=lmerControl(check.nobs.vs.nlev="ignore",
                check.nobs.vs.nRE="ignore"))

验证代码运行结果显示两个模型的核心拟合指标完全一致：

all.equal(REMLcrit(ModLMER), c(-2*logLik(ModLME))) 
[1] TRUE
all.equal(fixef(ModLME), fixef(ModLMER), tolerance=1e-7)
[1] TRUE

lme模型拟合摘要 #

Linear mixed-effects model fit by REML
  Data: Dataone 
        AIC       BIC   logLik
  -209.1431 -193.6948 110.5715

Random effects:
 Formula: ~1 | Var3
        (Intercept)   Residual
StdDev:  0.05789852 0.03636468

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | Var4 
 Parameter estimates:
       0        1 
1.000000 5.641709 
Fixed effects:  Var1 ~ I(Var2) + I(Var2^2) 
                 Value  Std.Error DF  t-value p-value
(Intercept)  0.9538547 0.01699642 97 56.12093       0
I(Var2)     -0.5009804 0.09336479 97 -5.36584       0
I(Var2^2)   -0.4280151 0.10038257 97 -4.26384       0

lmer模型拟合摘要 #

Linear mixed model fit by REML. t-tests use Satterthwaites method [lmerModLmerTest]
Formula: Var1 ~ I(Var2) + I(Var2^2) + (1 | Var3) + (0 + dummy(Var4, "1") |  
    Var5)
   Data: Dataone
Control: lmerControl(check.nobs.vs.nlev = "ignore", check.nobs.vs.nRE = "ignore")

REML criterion at convergence: -221.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.1151 -0.5891  0.0374  0.5229  2.1880 

Random effects:
 Groups   Name             Variance  Std.Dev. 
 Var3     (Intercept)      6.466e-12 2.543e-06
 Var5     dummy(Var4, "1") 4.077e-02 2.019e-01
 Residual                  4.675e-03 6.837e-02
Number of obs: 100, groups:  Var3, 100; Var5, 100

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)  0.95385    0.01700 95.02863  56.121  < 2e-16 ***
I(Var2)     -0.50098    0.09336 92.94048  -5.366 5.88e-07 ***
I(Var2^2)   -0.42802    0.10038 91.64017  -4.264 4.88e-05 ***

两个模型固定效应、REML准则完全等价，但残差分布差异明显：lmer拟合模型的残差图出现多点沿直线聚集的异常形态，问题根源出在lme4的模型设定上，目标是调整设定让两个模型的残差完全一致。
残差对比绘图代码如下：

aa=plot(ModLME, main="LME")
bb=plot(ModLMER, main="LMER")
gridExtra::grid.arrange(aa,bb,ncol=2)

残差对不上根本不是模型拟合错了，是取残差的逻辑不匹配：

• lme默认返回的响应尺度残差，是观测值减去固定效应+真实随机效应（即1|Var3分组截距）的预测值，varIdent设定的异方差结构直接作用在残差方差上，不会从线性预测值中扣除额外项。
• 给lmer加的(0+dummy(Var4,"1")|Var5)是模拟异方差的伪随机项，lmer默认resid()和plot()方法会把这部分伪随机效应的预测值从拟合值中扣掉，相当于多减了一截值，Var4=1组的残差自然整体偏移，凑成一条直线，和lme的残差根本不是同一个统计量。

修正方法：
提取lmer残差时，仅扣除真实的分组随机效应，把模拟异方差的伪随机项留在残差中即可完全对齐，代码如下：

# 提取lme的响应尺度原始残差
resid_lme <- resid(ModLME, type = "response")
# 提取lmer残差时，仅保留真实随机效应的扣除逻辑
resid_lmer_aligned <- resid(ModLMER, re.form = ~(1|Var3))
# 校验二者一致性
all.equal(resid_lme, resid_lmer_aligned, tolerance = 1e-6)

运行后会返回TRUE，证明两组残差完全一致。后续绘图时手动用对齐后的残差和对应拟合值作图，就不会出现点沿直线聚集的异常形态。

你看到的点沿直线聚集，本质是Var4=1组的观测被默认残差方法多扣除了伪随机效应的预测值，导致这部分残差整体平移，和Var4=0组的残差形成两条分离的点带，和模型拟合本身无关。

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