如何在R中提取并报告glmer模型的随机斜率及置信区间？

核心逻辑 #

个体的随机斜率是固定效应值 + 对应个体的随机效应偏差，但不能直接将固定效应的置信限与随机效应的condsd简单组合——因为固定效应置信区间针对的是总体平均水平的变异，而随机效应是个体相对于总体的偏差，两者的变异来源需要结合模型参数的分布来综合计算。最可靠的方式是通过模拟模型参数的后验分布，生成每个个体斜率的分布后提取置信区间。

具体实现步骤 #

1. 拟合模型（沿用你的代码） #

library(lme4)
library(dplyr)
library(tidyr)
library(arm) # 用于模拟参数后验分布

# 拟合二项式混合效应模型
mod <- glmer(dependent~pred1+pred2+(pred1+pred2|ID), data = df, family = "binomial")

2. 整理基础效应值 #

先把固定效应和个体随机效应整理成结构化数据：

# 提取固定效应点估计
fixed_df <- tibble(
  term = names(fixef(mod)),
  fixef = fixef(mod)
)

# 提取每个ID的随机效应偏差
random_df <- ranef(mod)$ID %>%
  rownames_to_column(var = "ID") %>%
  rename_with(~gsub("\\(Intercept\\)", "intercept", .x)) %>% # 重命名截距项方便后续处理
  pivot_longer(cols = -ID, names_to = "term", values_to = "condval")

3. 模拟参数分布计算置信区间 #

通过模拟模型参数的后验样本，生成每个个体斜率的分布，再取分位数得到置信区间：

# 模拟1000组模型参数样本（样本量可根据稳定性需求调整）
sim_params <- sim(mod, n.sim = 1000)

# 整理模拟的固定效应样本
sim_fixed <- t(sim_params@fixef) %>%
  as.data.frame() %>%
  rownames_to_column(var = "sim_id") %>%
  pivot_longer(cols = -sim_id, names_to = "term", values_to = "sim_fixef")

# 整理模拟的随机效应样本（每个ID的随机偏差）
sim_random <- sim_params@ranef$ID %>%
  lapply(function(x) as.data.frame(t(x)) %>% rownames_to_column(var = "ID")) %>%
  bind_rows(.id = "sim_id") %>%
  rename_with(~gsub("\\(Intercept\\)", "intercept", .x)) %>%
  pivot_longer(cols = -c(sim_id, ID), names_to = "term", values_to = "sim_condval")

# 合并模拟数据，计算每个样本下的个体斜率
sim_slopes <- sim_fixed %>%
  inner_join(sim_random, by = c("sim_id", "term")) %>%
  mutate(sim_slope = sim_fixef + sim_condval)

# 计算每个ID-变量组合的斜率点估计（中位数）和95%置信区间
slope_ci <- sim_slopes %>%
  group_by(ID, term) %>%
  summarise(
    rnd_slope = median(sim_slope),
    lower_ci = quantile(sim_slope, 0.025),
    upper_ci = quantile(sim_slope, 0.975),
    .groups = "drop"
  )

# 合并手动计算的斜率与模拟得到的置信区间
final_results <- random_df %>%
  inner_join(fixed_df, by = "term") %>%
  mutate(manual_slope = fixef + condval) %>%
  inner_join(slope_ci, by = c("ID", "term"))

4. 额外处理（可选） #

如果需要将对数优势比尺度的斜率转换为概率尺度，可应用plogis()函数：

final_results <- final_results %>%
  mutate(
    rnd_slope_prob = plogis(rnd_slope),
    lower_ci_prob = plogis(lower_ci),
    upper_ci_prob = plogis(upper_ci)
  )

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