如何从嵌套数据的多组线性模型中提取Wald p值

问题背景 #

我有16个二元候选变量，想要观察每个变量对15个结果变量的“贡献”，且每个候选变量对应的协变量组合不同。相关变量定义如下：

Candidates = c("C1","C2",...)
covariates = c("age","sex"...)
outcomes = c("Cancer","Heart disease",...)

现有操作流程 #

1. 按结果变量对数据进行嵌套处理：

ndat <- data |> group_by(outcome) |> nest()

2. 为每个候选变量定义对应模型（协变量随候选变量变化）：

model1 = function(df){
lm(value ~ C1 + age + sex, data = df)}

model2 = function(df){
lm(value ~ C2 + sex + ethnicity, data = df)}

3. 将所有模型应用到嵌套数据中：

ndat <- ndat |>
mutate(m1 = map(data, model1),
m2 = map(data, model2),
...
)

此时ndat是一个分组tibble：第一列为结果变量，其余列对应不同模型，列内元素是各结果对应的线性模型列表。

核心需求 #

需要批量提取每个模型中候选变量的Wald p值，最终实现：每列对应一个候选变量，每行对应一个结果变量，单元格为对应候选变量在该结果模型中的p值。

尝试过的方法及问题 #

• 使用broom::glance仅能返回整体模型指标，无法提取单个变量的p值；
• 单个提取可通过summary(ndat$model1[[1]])$coefficients[,4][[1]]实现，但无法规模化处理16个候选变量+15个结果变量的场景。

复现示例（mtcars数据集） #

注：示例代码未加载tidyverse包，因此出现函数未找到的错误

data <- mtcars |>
  mutate(c1  = rbinom(nrow(mtcars),prob=0.05, size = 1),
         c2 = rbinom(nrow(mtcars), prob = 0.1, size =1),
         c3 = rbinom(nrow(mtcars), prob = 0.5, size = 1))
#> Error in mutate(mtcars, c1 = rbinom(nrow(mtcars), prob = 0.05, size = 1), : could not find function "mutate"


candidates <- c("c1","c2","c3")
covars <- c("disp","hp","drat","wt")
outcomes <- c("mpg","qsec")

outcome_cols <- names(data)[names(data) %in% outcomes]

dat_long <- data |>
  pivot_longer(cols=all_of(outcome_cols), names_to = "outcome", values_to = "value")
#> Error in pivot_longer(data, cols = all_of(outcome_cols), names_to = "outcome", : could not find function "pivot_longer"

dat_n <- dat_long |>
  group_by(cyl) |>
  nest()
#> Error in nest(group_by(dat_long, cyl)): could not find function "nest"

c_models <- c("c1_mod","c2_mod","c3_mod")

c1_mod <- function(df){
  lm(value ~ c1 + disp + hp, data = df)
}

c2_mod <- function(df){
  lm(value ~ c2 + disp + drat, data = df)
}

c3_mod <- function(df){
  lm(value ~ c3 + drat + wt, data = df)
}

dat_n <- dat_n |>
  mutate(c1 = map(data, c1_mod),
         c2 = map(data, c2_mod),
         c3 = map(data, c3_mod))
#> Error in mutate(dat_n, c1 = map(data, c1_mod), c2 = map(data, c2_mod), : could not find function "mutate"

# 单个模型提取整体指标可行（需加载tidyverse）
glancec1 <- dat_n |>
  mutate(glance = map(c1, broom::glance)) |>
  unnest(glance) |>
  select(cyl, BIC, adj.r.squared) |>
  mutate(model = "c1")
#> Error in mutate(select(unnest(mutate(dat_n, glance = map(c1, broom::glance)), : could not find function "mutate"

# 尝试批量处理所有模型失败
glances <- data.frame()

for (i in c_models) {
  print(i)
  glance <- dat_n |>
    mutate(glance = map(i, broom::glance)) |>
    unnest(glance) |>
    select(cyl, BIC, adj.r.squared) |>
    mutate(model = as.character(i))
  
  glances <- bind_rows(glances, glance)
}
#> [1] "c1_mod"
#> Error in mutate(select(unnest(mutate(dat_n, glance = map(i, broom::glance)), : could not find function "mutate"

^{Created on 2023-06-23 with reprex v2.0.2}

解决方案 #

第一步：修复示例代码基础错误 #

首先必须加载tidyverse和broom包，否则会出现函数未找到的报错：

library(tidyverse)
library(broom)

第二步：批量提取候选变量的Wald p值 #

方法1：针对少量候选变量的直接提取

定义一个提取指定变量p值的函数，然后手动为每个模型列生成对应p值列：

# 定义提取指定变量p值的函数
extract_p <- function(model, var_name) {
  tidy(model) |>
    filter(term == var_name) |>
    pull(p.value)
}

# 批量处理每个模型列
ndat <- ndat |>
  mutate(
    p_C1 = map_dbl(m1, extract_p, var_name = "C1"),
    p_C2 = map_dbl(m2, extract_p, var_name = "C2"),
    # 按此格式添加所有16个候选变量的提取语句
  )

# 保留结果变量和p值列
result <- ndat |>
  select(outcome, starts_with("p_"))

方法2：针对大量候选变量的自动化处理

如果候选变量和模型列较多，用循环实现自动化批量处理，避免手动重复代码：

# 定义模型列名和对应候选变量名
model_cols <- paste0("m", 1:16)
candidate_vars <- paste0("C", 1:16)

# 批量生成p值列
for (i in seq_along(model_cols)) {
  ndat <- ndat |>
    mutate(!!paste0("p_", candidate_vars[i]) := map_dbl(!!sym(model_cols[i]), extract_p, var_name = candidate_vars[i]))
}

# 提取最终结果
result <- ndat |>
  select(outcome, starts_with("p_"))

针对mtcars示例的完整可运行代码

library(tidyverse)
library(broom)

# 生成数据
data <- mtcars |>
  mutate(c1  = rbinom(nrow(mtcars),prob=0.05, size = 1),
         c2 = rbinom(nrow(mtcars), prob = 0.1, size =1),
         c3 = rbinom(nrow(mtcars), prob = 0.5, size = 1))

candidates <- c("c1","c2","c3")
outcomes <- c("mpg","qsec")

# 转换为长格式并嵌套
dat_long <- data |>
  pivot_longer(cols=all_of(outcomes), names_to = "outcome", values_to = "value")

dat_n <- dat_long |>
  group_by(outcome) |>
  nest()

# 批量存储模型函数
model_funs <- list(
  c1 = function(df) lm(value ~ c1 + disp + hp, data = df),
  c2 = function(df) lm(value ~ c2 + disp + drat, data = df),
  c3 = function(df) lm(value ~ c3 + drat + wt, data = df)
)

# 批量应用模型到嵌套数据
dat_n <- dat_n |>
  mutate(
    across(all_of(names(model_funs)), ~map(data, model_funs[[cur_column()]]))
  )

# 定义提取p值的函数
extract_p <- function(model, var_name) {
  tidy(model) |>
    filter(term == var_name) |>
    pull(p.value)
}

# 批量提取每个候选变量的p值
dat_n <- dat_n |>
  mutate(
    across(all_of(names(model_funs)), ~map_dbl(.x, extract_p, var_name = cur_column()), .names = "p_{.col}")
  )

# 查看最终结果
dat_n |>
  select(outcome, starts_with("p_"))

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