如何在R语言线性混合模型中设置Q1*2001为参考组？

问题说明 #

需要将线性混合模型的参考类别指定为Q1*2001，以下是现有代码及示例数据：

原代码 #

library(lme4)       # linear mixed-effects models
library(lmerTest)   # test for linear mixed-effects models
library(gtsummary)

names(trajectories)

m1 <- lmer(distance ~ (as.factor(year))* Quintile + (1 | id), data = trajectories)

summary(m1)

m2 <- tbl_regression (m1,
                include= c("as.factor(year):Quintile"),
                pvalue_fun = ~ style_pvalue(.x, digits = 3),
                estimate_fun = function(x) gtsummary::style_sigfig(x, digits = 3),
                add_estimate_to_reference_rows = TRUE,
                )%>% 
  bold_p() %>%  
  bold_labels()

print(m2)

尝试的组合变量代码 #

m1 <- lmer(distance ~ yr_qun + (1 | id), data = trajectories)

summary(m1)

contrasts(trajectories$yr_qun) <- contr.treatment(50, base = 1)

m2 <- tbl_regression (m1,
                include= c("yr_qun"),
                pvalue_fun = ~ style_pvalue(.x, digits = 3),
                estimate_fun = function(x) gtsummary::style_sigfig(x, digits = 3),
                add_estimate_to_reference_rows = TRUE,
                )%>% 
  bold_p() %>%  
  bold_labels() %>% 
as_kable()

print(m2)

示例数据 #

trajectories <- structure(list(id = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L,
7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 10L,
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L), year = c(2001L,
2002L, 2003L, 2004L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L,
2001L, 2002L, 2003L, 2004L, 2005L, 2006L, 2007L, 2008L, 2009L,
2010L, 2001L, 2002L, 2003L, 2004L, 2005L, 2006L, 2007L, 2008L,
2009L, 2010L, 2001L, 2002L, 2003L, 2004L, 2005L, 2006L, 2007L,
2008L, 2009L, 2010L, 2001L, 2002L, 2003L, 2004L, 2005L, 2006L,
2007L, 2008L, 2009L, 2010L, 2001L, 2002L, 2003L, 2004L, 2005L,
2006L, 2001L, 2002L, 2003L, 2004L, 2005L, 2006L, 2007L, 2008L,
2009L, 2010L, 2004L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L,
2001L, 2002L, 2003L, 2004L, 2005L, 2006L, 2007L, 2008L, 2009L,
2010L), distance = c(15, 20, 21.5, 23, 21, 21.5, 24, 25.5, 20.5,
24, 21, 20, 21.5, 23, 21, 21.5, 24, 25.5, 20.5, 24, 21, 20, 21.5,
23, 21, 21.5, 24, 25.5, 20.5, 24, 21, 20, 21.5, 23, 21, 21.5,
24, 25.5, 20.5, 24, 21, 20, 21.5, 23, 21, 21.5, 24, 25.5, 20.5,
24, 21, 20, 21.5, 23, 21, 21.5, 21, 20, 21.5, 23, 21, 21.5, 24,
25.5, 20.5, 24, 23, 21, 21.5, 24, 25.5, 20.5, 24, 15, 20, 21.5,
23, 21, 21.5, 24, 25.5, 20.5, 24), age = c(8L, 9L, 10L, 11L,
12L, 13L, 14L, 15L, 16L, 17L, 11L, 12L, 13L, 14L, 15L, 16L, 17L,
18L, 19L, 20L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L,
24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 6L, 7L, 8L,
9L, 10L, 11L, 12L, 13L, 14L, 15L, 9L, 10L, 11L, 12L, 13L, 14L,
6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 18L, 19L, 20L,
21L, 22L, 23L, 24L, 28L, 40L, 41L, 42L, 43L, 44L, 45L, 46L, 47L,
48L), Quintile = structure(c(5L, 2L, 3L, 3L, 2L, 2L, 4L, 2L,
5L, 5L, 1L, 4L, 2L, 5L, 4L, 3L, 3L, 4L, 3L, 3L, 1L, 3L, 1L, 2L,
1L, 5L, 2L, 4L, 1L, 4L, 3L, 2L, 5L, 3L, 4L, 4L, 3L, 1L, 4L, 3L,
4L, 1L, 4L, 4L, 5L, 1L, 5L, 2L, 2L, 2L, 3L, 5L, 3L, 3L, 4L, 1L,
3L, 1L, 1L, 5L, 2L, 4L, 1L, 3L, 2L, 4L, 1L, 3L, 4L, 5L, 3L, 3L,
1L, 2L, 2L, 3L, 1L, 2L, 3L, 5L, 5L, 2L, 5L), levels = c("Q1",
"Q2", "Q3", "Q4", "Q5"), class = "factor"), sex = structure(c(2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L), levels = c("F", "M"), class = "factor")), class = "data.frame", row.names = c(NA,
-83L))

解决方案 #

方法一：调整单个因子的参考水平 #

通过分别设置year和Quintile的参考水平，让它们的交互项默认参考组为Q1*2001：

library(lme4)       
library(lmerTest)   
library(gtsummary)

# 重新设置year的参考水平为2001
trajectories$year <- factor(trajectories$year, levels = c(2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010))
# 重新设置Quintile的参考水平为Q1（原本已是Q1，可省略，若需确认则保留）
trajectories$Quintile <- factor(trajectories$Quintile, levels = c("Q1", "Q2", "Q3", "Q4", "Q5"))

# 拟合模型，交互项参考组自动为Q1*2001
m1 <- lmer(distance ~ year * Quintile + (1 | id), data = trajectories)
summary(m1)

# 生成回归表格
m2 <- tbl_regression(m1,
                     include= c("year:Quintile"),
                     pvalue_fun = ~ style_pvalue(.x, digits = 3),
                     estimate_fun = function(x) gtsummary::style_sigfig(x, digits = 3),
                     add_estimate_to_reference_rows = TRUE) %>% 
  bold_p() %>%  
  bold_labels()

print(m2)

方法二：创建组合因子并指定参考组 #

将year和Quintile合并为单个因子，直接指定2001_Q1作为参考组，更直观可控：

# 创建组合因子，分隔符为下划线
trajectories$yr_qun <- interaction(trajectories$year, trajectories$Quintile, sep = "_")
# 将2001_Q1设为参考组
trajectories$yr_qun <- relevel(factor(trajectories$yr_qun), ref = "2001_Q1")

# 拟合模型
m1 <- lmer(distance ~ yr_qun + (1 | id), data = trajectories)
summary(m1)

# 生成回归表格
m2 <- tbl_regression(m1,
                     include= c("yr_qun"),
                     pvalue_fun = ~ style_pvalue(.x, digits = 3),
                     estimate_fun = function(x) gtsummary::style_sigfig(x, digits = 3),
                     add_estimate_to_reference_rows = TRUE) %>% 
  bold_p() %>%  
  bold_labels() %>% 
  as_kable()

print(m2)

注意事项 #

• 方法一中，确保year的第一个水平是2001，Quintile的第一个水平是Q1，交互项的参考组会自动对应两者的组合。
• 方法二中，interaction()生成的组合格式为年份_Quintile，若修改分隔符，参考组名称需同步调整。
• 两种方法均可实现需求，可根据分析习惯选择。

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