使用Sandwich、merDeriv配合lme4::glmer计算稳健标准误报错求助

需要为lme4包glmer()函数拟合的广义线性混合模型计算稳健标准误，此前使用sandwich+merDeriv组合的代码可正常运行，当前执行时触发报错。
原运行代码：

a = lme4::glmer(gholocaust ~ cg06131859_KYNU + IC + gender + age+ (1|cID), 
     data = data, 
     family = "binomial", nAGQ = 5L)
sandwich::sandwich(a, bread. = bread.glmerMod, meat. = meat(a, level = 2))

报错信息 #

Error in (function (cond) :
error in evaluating the argument 'x' in selecting a method for function 't': subscript out of bounds

分步验证结果 #

• 单独执行glmer()模型拟合代码无报错，模型正常收敛，输出如下：

Generalized linear mixed model fit by maximum likelihood (Adaptive Gauss-Hermite Quadrature, nAGQ = 5) ['glmerMod']
 Family: binomial  ( logit )
Formula: gholocaust ~ cg06131859_KYNU + IC + gender + age + (1 | cID)
   Data: data
     AIC      BIC   logLik deviance df.resid 
136.2090 151.4046 -62.1045 124.2090       87 
Random effects:
 Groups Name        Std.Dev.
 cID    (Intercept) 0.5122  
Number of obs: 93, groups:  cID, 47
Fixed Effects:
    (Intercept)  cg06131859_KYNU               IC        genderאשה              age  
        2.53816         -5.72758         -0.83361          0.61304          0.00775  

• 单独调用merDeriv的bread.glmerMod()计算面包矩阵可正常运行，输出为6×6矩阵（包含5个固定效应参数、1个随机效应方差参数）：

(Intercept) cg06131859_KYNU           IC   genderאשה          age cov_cID.(Intercept)
[1,]  729.518988    -384.5030066 -320.5693995 -17.1856860 -10.10491538         -3.45021545
[2,] -384.503007     530.5340435  181.3728027 -16.2948357   0.57913609        -34.97224468
[3,] -320.569399     181.3728027  284.3496470   4.7488551   0.23081025         17.61985923
[4,]  -17.185686     -16.2948357    4.7488551  11.0148593   0.50889828          5.55635382
[5,]  -10.104915       0.5791361    0.2308103   0.5088983   0.32781282          0.02613785
[6,]   -3.450215     -34.9722447   17.6198592   5.5563538   0.02613785         57.24944251

1. S3方法未注册：代码中没有显式加载merDeriv包，仅直接调用bread.glmerMod，导致merDeriv为glmerMod类编写的meat.glmerMod()方法没有被注册到R的S3方法调度表里。执行meat(a, level=2)时，R会调用sandwich包的默认meat()函数，该函数不识别glmerMod对象的结构，计算出的meat矩阵仅包含固定效应对应的5行5列，和bread返回的6×6矩阵维度不匹配，矩阵运算时触发下标越界。
2. 参数传递方式错误：sandwich()函数的meat.参数要求传入计算meat矩阵的函数对象，原代码直接传入meat(a, level=2)的求值结果，没有给sandwich内部做参数对齐、维度校验的空间，进一步放大了维度不匹配的问题。

1. 显式加载merDeriv包，确保包内置的glmerMod类对应的bread、meat S3方法被正确注册，不要仅用::调用单个函数。
2. 修正sandwich()的参数传递逻辑，将meat对应的方法函数传入，把集群层级参数level=2作为额外参数传给sandwich()，由函数内部调度正确的方法计算、对齐矩阵维度。

修正后的可运行代码：

# 加载所有依赖包，注册S3方法
library(lme4)
library(sandwich)
library(merDeriv)
# 拟合模型
a <- glmer(gholocaust ~ cg06131859_KYNU + IC + gender + age + (1|cID), 
           data = data, 
           family = "binomial", 
           nAGQ = 5L)
# 计算集群稳健方差-协方差矩阵（按cID集群，对应level=2）
robust_vcov <- sandwich(
  a,
  bread. = bread.glmerMod,
  meat. = meat.glmerMod,
  level = 2
)
# 提取稳健标准误直接取对角元开平方即可
robust_se <- sqrt(diag(robust_vcov))

如果执行上述代码仍报错，将merDeriv更新到0.2-4及以上版本、sandwich更新到3.0-2及以上版本即可，旧版本存在两类对象方法适配的已知bug。

内容的提问来源于stack exchange，提问作者Tsachi Ein Dor