生成负二项模型三向交互表：R语言工具调用遇阻求助

模型背景 #

拟合的负二项最终模型包含2个显著三向交互项，模型代码如下：

Final_model_acetotal_demo_17_20 <-  glm.nb(acetotal ~ (group + fed_pov_level + gender + race_eth + metro)^2 
            + group:race_eth:metro
            + group:race_eth:fed_pov_level, 
            data=NSCH2_2017_2020)

遇到的问题 #

尝试多种工具生成美观表格均未达到预期：

• 使用tbl_regression(Final_model_acetotal_demo_17_20, exponentiate = T, intercept = T)：程序长时间运行无输出
• 使用tab_model(Final_model_acetotal_demo_17_20, ci_method="wald")：提示模型矩阵秩亏（参数groupDHH:race_ethOther:metroYes无法估计），随后程序无响应
• 使用texreg::knitreg(...)：仅能生成优于控制台的简易表格，美观度不足

可行解决方案 #

1. 先解决秩亏问题（根源性处理） #

秩亏是导致多数工具卡顿/报错的核心原因，先排查并修复：

• 检查交互项对应分组是否存在零频数单元格：

# 生成三变量交叉频数表
table(NSCH2_2017_2020$group, NSCH2_2017_2020$race_eth, NSCH2_2017_2020$metro)

若存在零频数组，可合并小众类别（比如将race_eth的低频次类别合并为"Other"），或删除对应观测后重新拟合模型，再尝试表格生成。

2. 用gtsummary进阶参数优化运行效率 #

tbl_regression卡顿多因默认处理全量变量，可通过参数限制输出范围：

library(gtsummary)
Final_model_acetotal_demo_17_20 %>%
  tbl_regression(
    exponentiate = TRUE,
    intercept = TRUE,
    # 仅保留核心变量与目标交互项，减少计算量
    include = c(group:race_eth:metro, group:race_eth:fed_pov_level, group, race_eth, metro, fed_pov_level),
    # 自定义变量标签提升可读性
    label = list(
      group ~ "Group",
      race_eth ~ "Race/Ethnicity",
      metro ~ "Metro Status",
      fed_pov_level ~ "Poverty Level"
    )
  ) %>%
  add_significance_stars() %>%
  as_gt() %>%
  gt::tab_options(table.width = gt::pct(100)) # 适配页面宽度

3. 使用modelsummary生成高质量表格 #

modelsummary对复杂交互项支持更友好，运行效率更高：

library(modelsummary)
modelsummary(Final_model_acetotal_demo_17_20,
             exponentiate = TRUE,
             conf_level = 0.95,
             stars = TRUE,
             # 自定义交互项名称，避免默认的冗长标签
             coef_rename = c(
               "groupDHH" = "Group: DHH",
               "race_ethOther" = "Race/Ethnicity: Other",
               "metroYes" = "Metro: Yes",
               "groupDHH:race_ethOther:metroYes" = "Group: DHH × Race: Other × Metro: Yes",
               "groupDHH:race_ethOther:fed_pov_levelHigh" = "Group: DHH × Race: Other × Poverty: High"
             ),
             output = "gt") # 输出为gt格式，可直接导出为HTML/Word

若仍有秩亏提示，可使用稳健标准误并跳过无法估计的参数：

modelsummary(Final_model_acetotal_demo_17_20,
             exponentiate = TRUE,
             conf_level = 0.95,
             stars = TRUE,
             vcov = "HC3",
             coef_omit = "groupDHH:race_ethOther:metroYes", # 跳过无法估计的参数
             output = "interaction_table.html") # 直接导出HTML文件

4. 手动提取结果定制表格 #

若上述工具仍有问题，可手动提取交互项结果，用gt包定制：

library(broom)
library(gt)

# 提取模型的OR值、置信区间与P值
model_results <- tidy(Final_model_acetotal_demo_17_20, exponentiate = TRUE, conf.int = TRUE)

# 筛选出三向交互项
interaction_results <- model_results[grepl(":", model_results$term) & stringr::str_count(model_results$term, ":") == 2, ]

# 生成定制化表格
gt(interaction_results) %>%
  cols_label(
    term = "Interaction Term",
    estimate = "OR",
    conf.low = "95% CI Lower",
    conf.high = "95% CI Upper",
    p.value = "P-value"
  ) %>%
  fmt_number(columns = c(estimate, conf.low, conf.high), decimals = 2) %>%
  fmt_pvalue(columns = p.value, decimals = 3) %>%
  tab_header(title = "Three-Way Interaction Results of Negative Binomial Model") %>%
  opt_table_lines() %>%
  gtsave("custom_interaction_table.html") # 导出为HTML

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