使用nonnest2对lme4非嵌套线性混合模型执行Vuong测试遇内存问题

问题背景 #

用lme4拟合了两个非嵌套混合效应模型，计划通过nonnest2包的Vuong检验（1989）比较模型优劣，但运行时频繁出现**“无法分配76.4 Gb大小的向量”**错误或脚本直接终止。模型保存为RDS仅5.4MB，参考官方文档未找到可用示例，代码如下：

library(lme4)
library(nonnest2)
library(merDeriv)

m_1 <- lmer(LogRT ~ surprisal_S + covariate + (1|SubjectID) + (1 | UniqueWordID), data=df, REML=FALSE)

m_2 <- lmer(LogRT ~ surprisal_T + covariate + (1|SubjectID) + (1 | UniqueWordID), data=df, REML=FALSE)

# follow the docs
vcl <- function(obj) vcov(obj, full=TRUE)
vuongtest(m_1, m_2, vc1=vcl, vc2=vcl, nested=FALSE)

核心原因 #

代码中指定vcov(obj, full=TRUE)会生成包含所有随机效应分组水平的全协方差矩阵，当SubjectID或UniqueWordID的分组水平数量较多（比如上万级）时，矩阵规模会呈平方级增长，直接耗尽内存。

可行解决方案 #

1. 使用固定效应协方差矩阵替代全矩阵 #

Vuong检验的核心逻辑可基于固定效应参数的协方差矩阵实现，避免生成庞大的随机效应协方差矩阵。修改协方差函数：

# 仅提取固定效应的协方差矩阵
vcl_fixed <- function(obj) vcov(obj)
vuongtest(m_1, m_2, vc1=vcl_fixed, vc2=vcl_fixed, nested=FALSE)

2. 验证数据规模 #

先检查分组变量的水平数量：

length(unique(df$SubjectID))
length(unique(df$UniqueWordID))

如果任一变量水平数超过1000，全协方差矩阵的大小会达到百万级以上，必然触发内存溢出，此时优先使用方案1。

3. 小样本测试代码逻辑 #

随机抽取小比例数据（比如10%的受试者）验证代码是否能正常运行，排除语法或其他逻辑错误：

# 随机抽取10%的SubjectID
sample_subj <- sample(unique(df$SubjectID), size=0.1*length(unique(df$SubjectID)))
df_sample <- df[df$SubjectID %in% sample_subj, ]

# 重新拟合模型并运行检验
m_1_sample <- lmer(LogRT ~ surprisal_S + covariate + (1|SubjectID) + (1 | UniqueWordID), data=df_sample, REML=FALSE)
m_2_sample <- lmer(LogRT ~ surprisal_T + covariate + (1|SubjectID) + (1 | UniqueWordID), data=df_sample, REML=FALSE)

vuongtest(m_1_sample, m_2_sample, vc1=vcl_fixed, vc2=vcl_fixed, nested=FALSE)

4. 手动实现Vuong检验（备选） #

如果nonnest2包的内存问题无法解决，可手动计算核心统计量（大样本下结果与官方实现近似）：

library(merDeriv)

# 获取每个观测的对数似然值
ll1 <- logLikVec(m_1)
ll2 <- logLikVec(m_2)
ll_diff <- ll1 - ll2

# 计算Vuong统计量与p值
n <- length(ll_diff)
vuong_stat <- mean(ll_diff) / (sd(ll_diff)/sqrt(n))
p_val <- 2 * pnorm(-abs(vuong_stat))

# 输出结果
cat("Vuong统计量：", round(vuong_stat, 3), "\n")
cat("双侧p值：", round(p_val, 4), "\n")

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