从混合模型提取分类系数与p值，生成指定结构结果数据表

需求说明 #

生成包含3列的最终数据表：

1. 暴露分位数（带自定义标签）
2. OR/RR（注：因本例结果为连续变量，实际提取的是回归系数；若为二分类结果，可通过exp()转换为OR）
3. P值（PV）

完整可复现代码 #

set.seed(42) 
n <- 100
dat = data.frame(ID = rep(c(1:25),times=4 ) ,
                 Score = rnorm(n, mean=0.3, sd=0.8))

dat = dat %>%
  group_by(ID)%>%
  dplyr::mutate(exposure1 = rep(c(rnorm(1, mean=6, sd=1.8))),
                exposure2 = rep(c(rnorm(1, mean=3, sd=0.6))),
                age = rep(c(rnorm(1, mean=40, sd=15))))%>%
  ungroup()%>%
  dplyr::mutate(exposure1_quantile = cut(exposure1, breaks = 4, labels = c("Q1","Q2","Q3","Q4")),
                exposure2_quantile = cut(exposure2, breaks = 4, labels = c("Q1","Q2","Q3","Q4")))

# 修正标签定义：改为向量赋值
exposures_var = c("exposure1_quantile","exposure2_quantile")
exposure_var_labels <- c("exposure1 Q1","exposure1 Q2", "exposure1 Q3", 
                         "exposure2 Q1","exposure2 Q2", "exposure2 Q3")
age="age"
outcome = "Score"
exposure_data_table = data.table()

for(i in 1:length(exposures_var)){
  exp = exposures_var[i]
  fixed_effects_formula = as.formula(paste0(outcome, "~",exp,"+",age))
  mixedmodel = lme(fixed = fixed_effects_formula, random = ~1|ID, data=dat, method = "ML")
  
  # 从模型summary提取固定效应的系数与P值
  model_coef <- summary(mixedmodel)$tTable
  
  # 提取当前暴露的哑变量行（排除截距和age）
  exp_rows <- grep(exp, rownames(model_coef))
  
  # 匹配对应的暴露标签
  label_start <- (i-1)*3 +1
  label_end <- label_start +2
  current_labels <- exposure_var_labels[label_start:label_end]
  
  # 组装当前暴露的结果行
  current_result <- data.table(
    `暴露分位数` = current_labels,
    `OR/RR` = round(model_coef[exp_rows, "Value"], 3),
    `PV` = round(model_coef[exp_rows, "p-value"], 4)
  )
  
  exposure_data_table <- rbind(exposure_data_table, current_result)
}

# 查看结果
print(exposure_data_table)

关键修正说明 #

1. 标签定义修正：原代码中exposure_var_labels未正确赋值，改为向量形式存储自定义标签，对应两个暴露的3个哑变量（Q2/Q3/Q4对比Q1）。
2. P值提取：通过summary(mixedmodel)$tTable直接获取固定效应每个系数的P值，anova通常用于检验暴露的整体效应，而非单个哑变量的P值。
3. 标签匹配：按暴露变量的索引，从标签向量中提取对应位置的标签，确保分位数标签与模型系数一一对应。
4. 结果表整理：用data.table直接构建结果行，避免多次rbind的性能问题，同时规范列名与数值格式。

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