如何将MuMIn的model.avg对象转为含相同系数的lme4模型

核心思路 #

由于averaging类对象无法被部分lme4兼容的函数识别，我们需要手动构造一个lme4模型对象，替换其中的固定效应系数、随机效应参数（可选）和残差方差，使其匹配模型平均的结果。注意：MuMIn的dredge默认固定随机效应结构不变，仅筛选固定效应，因此模型平均通常仅针对固定效应，随机效应可选择沿用原始结构或对候选模型的随机参数做平均。

具体步骤与代码 #

1. 提取模型平均的关键参数 #

从model.avg结果中提取平均后的固定效应系数，若需要平均随机效应，还需从候选模型中提取相关参数：

library(lme4)
library(MuMIn)
options(na.action = na.fail)

# 原始模型与模型平均流程
m1 <- lmer(distance ~ age * Sex + (1|Subject), data = Orthodont)
d1 <- dredge(m1)
a1 <- model.avg(d1)

# 提取包含所有固定项的平均系数（full=TRUE确保保留所有候选固定效应）
fixed_coef <- coef(a1, full = TRUE)

# 若需要平均随机效应，先提取所有候选模型
candidate_models <- get.models(d1, subset = TRUE)

2. 创建模板lme4模型 #

构建一个包含所有固定效应项的模板模型，确保其结构与模型平均的full模型一致：

# 模板模型需包含所有固定效应项（与fixed_coef的变量对应）
m2 <- lmer(distance ~ age + Sex + age:Sex + (1|Subject), data = Orthodont)

3. 替换固定效应系数 #

直接修改lme4模型对象的@beta属性，替换为平均后的固定效应：

# 替换固定效应系数
m2@beta <- fixed_coef

4. 处理随机效应（可选） #

若需要对随机效应的条件模式（即ranef结果）或方差参数做平均，可按以下方式替换：

# 平均随机效应条件模式（每个Subject的截距）
ranef_list <- lapply(candidate_models, function(x) ranef(x)$Subject)
avg_ranef <- colMeans(do.call(rbind, ranef_list))
m2@u <- as.matrix(avg_ranef)  # 替换模型的随机效应条件模式

# 平均随机效应方差参数（theta）
avg_theta <- mean(sapply(candidate_models, function(x) x@theta))
m2@theta <- avg_theta

# 平均残差方差（sigma）
avg_sigma <- mean(sapply(candidate_models, function(x) x@sigma))
m2@sigma <- avg_sigma

5. 更新方差协方差矩阵 #

修改参数后，需重新计算模型的方差协方差矩阵，避免后续函数调用出错：

m2@vcov <- vcov(m2)

注意事项 #

• 手动修改lme4的S4对象属于非官方操作，部分依赖模型内部结构的函数可能出现异常，建议测试关键功能（如predict、summary）。
• 若仅需要固定效应的平均结果，可跳过随机效应的平均步骤，直接保留模板模型的随机效应参数。
• 务必确保模板模型的固定效应项与fixed_coef的变量完全匹配，否则系数赋值会因维度不匹配报错。

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