R调用tbl_uvregression报错C栈使用率接近上限的解决方案

问题描述 #

有一个包含2946条观测、600个变量的数据框，需要针对其中599个变量生成单因素logistic回归结果表，使用gtsummary包的tbl_uvregression()函数代码如下：

RAPOA_labelled[,-1] %>%    #移除ID列
  tbl_uvregression(
    method = glm,
    y = GIR.2cat,          #因变量
    method.args = list(family = binomial),
    exponentiate = TRUE,
    pvalue_fun = ~style_pvalue(.x, digits = 3)
  ) %>%
  add_nevent() %>%         #添加结局事件发生数
  bold_p() %>%             #p<0.05加粗
  bold_labels()

运行时触发报错：

Error: C stack usage 7971168 is too close to the limit.

查询C栈信息结果：

> Cstack_info()
   size    current   direction eval_depth 
7969177      12800           1          2 

报错原因 #

tbl_uvregression()在设计上没有针对数百级别的批量变量做性能优化，处理大量变量时内部嵌套调用会超过R默认的C栈容量上限，触发溢出报错。

解决方案 #

方案1：分批次调用gtsummary函数（保留原包格式化能力） #

将599个变量拆分为多组，每组50-60个变量分别生成回归表，再合并结果：

# 提取自变量列表（排除ID和因变量）
all_pred <- colnames(RAPOA_labelled)[!colnames(RAPOA_labelled) %in% c("ID", "GIR.2cat")]
# 拆分为10组
pred_groups <- split(all_pred, cut(seq_along(all_pred), 10, labels = FALSE))
# 逐组生成单因素回归表
tbl_list <- lapply(pred_groups, function(group_vars) {
  RAPOA_labelled[, c("GIR.2cat", group_vars)] %>%
    tbl_uvregression(
      method = glm,
      y = GIR.2cat,
      method.args = list(family = binomial),
      exponentiate = TRUE,
      pvalue_fun = ~style_pvalue(.x, digits = 3)
    )
})
# 合并所有表格并添加格式
final_uv_table <- tbl_stack(tbl_list) %>%
  add_nevent() %>%
  bold_p() %>%
  bold_labels()

方案2：直接使用循环批量拟合（性能更高） #

如果不需要gtsummary的实时格式化能力，可以使用自定义循环批量拟合回归，你提供的循环代码已经可以实现需求，也可以优化为purrr的map写法提升代码整洁度，生成的结果数据框可以直接导出为csv，或者传入as_flextable()等函数生成发表格式表格：

library(purrr)
library(broom)
all_pred <- colnames(datos_rapoa_gir)[!colnames(datos_rapoa_gir) %in% c("ID", "GIR.2cat")]
univars <- map_dfr(all_pred, function(var) {
  mod <- glm(reformulate(var, response = "GIR.2cat"), 
             family = binomial, data = datos_rapoa_gir, na.action = na.omit)
  tidy_mod <- tidy(mod, conf.int = TRUE, exponentiate = FALSE)
  tidy_mod$OR <- exp(tidy_mod$estimate)
  tidy_mod$LowIC <- exp(tidy_mod$conf.low)
  tidy_mod$HighIC <- exp(tidy_mod$conf.high)
  return(tidy_mod[, c("term", "estimate", "std.error", "p.value", "OR", "LowIC", "HighIC")])
})

方案3：调整R栈上限（不推荐） #

仅临时解决问题，变量进一步增加时仍会报错，Windows可在R快捷方式属性的目标栏末尾添加--max-ppsize=100000，macOS/Linux可在终端执行ulimit -s 16384后启动R。 #

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