使用reformulate/as.formula循环构建Cox模型时tidy_add_n()失效求助

问题背景 #

循环运行调整不同协变量子集的Cox模型，期望输出包含事件数、人年数、风险比（HR）及95%置信区间。直接调用coxph()时tidy_add_n()功能正常；但在循环中使用reformulate或as.formula构建公式时，该函数无法生成n_event和exposure列，反而生成全为缺失值的n列。

失效原因 #

tidy_add_n()（来自broom.helpers包）依赖模型对象中保留的原生公式结构：

• 直接调用coxph()时，公式是原生的Surv(time, status) ~ ...形式，模型对象能正确识别生存分析的响应变量类型，因此tidy_add_n()可正常提取事件数和人年数。
• 用字符串"Surv(time, status)"作为响应传入reformulate时，生成的公式是字符串解析后的形式，模型对象无法正确识别响应变量的生存分析结构，导致tidy_add_n()无法提取所需统计量。

解决方法 #

方案1：构建原生Surv对象的公式（最优解） #

不要用字符串形式的Surv表达式，直接构建包含原生Surv()的公式，让模型对象保留正确的响应变量结构：

test1 <- list(time = c(4,3,1,1,2,2,3), 
              status = c(1,1,1,0,1,1,0), 
              x1 = as.factor(c(0,2,1,1,1,0,0)), 
              x2 = as.factor(c(0,0,0,0,1,1,1)))

# 定义协变量子集
covariates <- c("x1","x2")
# 构建原生Surv公式
formula <- as.formula(paste("Surv(time, status) ~", paste(covariates, collapse = " + ")))

coxph(formula, data = test1) %>%
  tidy_and_attach() %>%
  tidy_add_reference_rows() %>%
  tidy_add_estimate_to_reference_rows() %>%
  tidy_add_header_rows() %>%
  tidy_add_n()

方案2：手动计算事件数和人年数 #

若特殊场景下无法使用原生公式，可从原始数据中手动计算所需统计量，再合并到结果中：

test1 <- list(time = c(4,3,1,1,2,2,3), 
              status = c(1,1,1,0,1,1,0), 
              x1 = as.factor(c(0,2,1,1,1,0,0)), 
              x2 = as.factor(c(0,0,0,0,1,1,1)))

model_data <- test1
# 计算总事件数
total_n_event <- sum(model_data$status)
# 计算总人年数（若需精确统计，可使用survival包的pyears函数）
total_exposure <- sum(model_data$time)

covariates <- c("x1","x2")
reformulate(covariates, response = "Surv(time, status)") %>%
  coxph(data = model_data) %>%
  tidy_and_attach() %>%
  tidy_add_reference_rows() %>%
  tidy_add_estimate_to_reference_rows() %>%
  tidy_add_header_rows() %>%
  # 添加手动计算的统计列
  mutate(n_event = total_n_event, exposure = total_exposure)

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