线性混合模型回归参考组定义求助：逐年对比距离时间影响

核心思路 #

无需手动构建yr_qun交互项，直接让R自动处理年份（year）与暴露分组（qun）的交互关系，通过重新指定因子水平的顺序，让每个年份对应的low组成为该年份内的参考基准。

步骤与代码示例 #

1. 准备模拟数据 #

先生成符合研究场景的模拟数据（可替换为你的真实数据）：

# 模拟数据
set.seed(123)
n <- 1000
year <- rep(2001:2005, each = n/5)
qun <- sample(c("low", "medium", "high"), n, replace = TRUE)
outcome <- rnorm(n, mean = ifelse(qun == "low", 0, ifelse(qun == "medium", 0.5, 1)) + (year-2001)*0.1, sd = 1)
id <- rep(1:(n/10), each = 10) # 个体ID用于混合效应
df <- data.frame(year = factor(year), qun = factor(qun), outcome = outcome, id = factor(id))

2. 重新编码交互因子，指定参考组 #

将year和qun合并为新因子，调整水平顺序让每个年份的low组排在最前面（成为参考组）：

# 创建年份-分组合并因子
df$yr_qun <- interaction(df$year, df$qun, sep = "_")

# 自定义因子水平顺序：每个年份的low组优先
levels_order <- unlist(lapply(unique(df$year), function(y) {
  c(paste(y, "low", sep = "_"), paste(y, c("medium", "high"), sep = "_"))
}))
df$yr_qun <- factor(df$yr_qun, levels = levels_order)

3. 拟合线性混合模型 #

使用lme4包拟合模型，此时每个年份的low组会自动作为该年份的参考组：

library(lme4)
model <- lmer(outcome ~ yr_qun + (1|id), data = df)
summary(model)

查看输出结果时，2001_medium、2001_high是相对于2001_low的估计值，2002_medium、2002_high是相对于2002_low的估计值，以此类推。

可选：用emmeans验证组间对比 #

如果需要更直观的组间差异结果，可借助emmeans包输出每个年份内的分组对比：

library(emmeans)
emm <- emmeans(model, ~ qun | year)
pairs(emm) # 输出每个年份中medium/high组与low组的对比结果

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