如何高效将同一glmer模型拟合到数据多子集并提取系数？

问题背景 #

需要按变量的不同水平拆分数据集，为每个子集拟合结构一致的复杂逻辑斯蒂混合效应模型（glmer），并高效提取各模型的系数。示例模型结构如下：

mod1 <- glmer(outcome ~
                 var1*var0 +
                 var2*var0 +
                 var3+
                 var4+
                 (1+var0|subj),
               family=binomial("logit"),
               control=glmerControl(calc.derivs=FALSE,
                                    optimizer="bobyqa",
                                    optCtrl=list(maxfun=2e5)),
               data=subset(mydata, var5 == T))

以sleepstudy数据集为例，已创建group和sex变量，目标是按sex的两个水平拆分数据，分别拟合模型并提取系数：

data(sleepstudy)
sleepstudy$group <- car::recode(sleepstudy$Subject,
   "c('308','309','310','330','331','332','333','334') = 'group1'; 
     c('335','337','349','350','351','352','369','370','371','372') = 'group2'")
sleepstudy$sex <- car::recode(sleepstudy$Subject,
                   "c('308','310','331','333','335','349','351','369','371') = 'm';
   c('309','330','332','334','337','350','352','370','372') = 'f'")

# 单模型示例
model <- glmer(group~Days+Reaction+(1|Subject), 
                family=binomial("logit"),
                data=sleepstudy)

高效解决方案 #

方法1：tidyverse工具链（dplyr + purrr + broom.mixed） #

代码简洁，输出格式规整，便于后续分析：

library(tidyverse)
library(lme4)
library(broom.mixed)

# 按sex分组、拟合模型、提取固定效应系数
coef_by_sex <- sleepstudy %>%
  group_nest(sex) %>%  # 按sex嵌套数据集
  mutate(
    # 批量拟合模型，可直接添加control参数
    model = map(data, ~glmer(group~Days+Reaction+(1|Subject), 
                            family=binomial("logit"),
                            data=.x,
                            control=glmerControl(calc.derivs=FALSE,
                                                 optimizer="bobyqa",
                                                 optCtrl=list(maxfun=2e5)))),
    # 提取系数及统计量
    coefs = map(model, tidy, effects = "fixed")
  ) %>%
  unnest(coefs) %>%
  select(sex, term, estimate, std.error, statistic, p.value)  # 保留需要的列

print(coef_by_sex)

• 若需提取随机效应，将tidy函数的effects参数改为"all"或"random"即可。

方法2：Base R实现 #

适合习惯base R语法的场景：

library(lme4)

# 按sex拆分数据集
split_data <- split(sleepstudy, sleepstudy$sex)

# 定义模型拟合函数
fit_glmer <- function(df) {
  glmer(group~Days+Reaction+(1|Subject), 
        family=binomial("logit"),
        data=df,
        # 添加控制参数
        control=glmerControl(calc.derivs=FALSE,
                             optimizer="bobyqa",
                             optCtrl=list(maxfun=2e5)))
}

# 批量拟合模型
models <- lapply(split_data, fit_glmer)

# 提取固定效应系数并整理为数据框
coef_df <- do.call(rbind, lapply(models, function(m) {
  tibble::enframe(fixef(m), name = "term", value = "estimate")
}))
# 添加分组标识
coef_df$sex <- rep(names(models), each = length(fixef(models[[1]])))

print(coef_df)

注意事项 #

• 模型收敛问题：若遇到拟合警告，可调整glmerControl中的优化参数（如增大maxfun、更换optimizer）。
• 复杂模型适配：两种方法都支持任意复杂度的glmer模型，只需将模型公式和参数对应替换即可。

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