如何为gtsummary多模型回归表添加年份分组的多级表头？

我需要为包含多个回归模型的gtsummary表格生成多级表头，按年份对模型进行分组。现有示例代码如下：

library(dplyr)
library(gtsummary)
library(purrr)

set.seed(92922)
df <- tibble(y_1980 = rbinom(n = 10, size = 1, prob = .4),
             y_1990 = rbinom(n = 10, size = 1, prob = .7),
             x1 = rnorm(10, sd = 1),
             x2 = rnorm(10, sd = 2))

tbls <- c("y_1980 ~ x1", "y_1980 ~ x1 + x2", "y_1990 ~ x1", "y_1990 ~ x1 + x2") %>% 
  map(as.formula) %>% 
  map(glm,
      data = df, 
      family = binomial(link = "logit")) %>% 
  map(tbl_regression, exponentiate = TRUE) %>% 
  map(add_significance_stars, hide_ci = TRUE, hide_p = TRUE, hide_se = FALSE) %>% 
  map(add_glance_table, include = nobs) 

目前我能通过以下代码生成单级表头的表格：

tbls %>% 
  tbl_merge(tab_spanner = c("1980 (1)", "1980 (2)", "1990 (1)", "1990 (2)")) %>% 
  modify_table_body(~.x %>% dplyr::arrange(row_type == "glance_statistic"))

但需要实现两级表头：年份在上，模型编号在下，格式如下：

1980      1990
_______   ______
(1) (2)   (1) (2)

要实现多级表头，可利用tbl_merge()的参数特性，将同一年份的模型绑定到同一主表头下，再设置二级表头显示模型编号，具体代码如下：

tbls %>%
  tbl_merge(
    # 用命名列表指定一级分组表头，键为年份，值为对应模型的列标签
    tab_spanner = list(
      "1980" = c("1980 (1)", "1980 (2)"),
      "1990" = c("1990 (1)", "1990 (2)")
    ),
    # 为每个模型列设置二级表头（仅显示模型编号）
    tab_spanner_label = c("(1)", "(2)", "(1)", "(2)")
  ) %>%
  # 保持原有行排序逻辑，将glance统计行置于底部
  modify_table_body(~.x %>% dplyr::arrange(row_type == "glance_statistic")) %>%
  # 可选：给一级表头添加加粗样式，增强层级区分
  modify_spanner_header(
    all_spanners(),
    ~paste0("**", .x, "**")
  )

代码说明： #

1. tab_spanner命名列表：将同一年份的模型列归到对应年份的一级表头下，实现分组。
2. tab_spanner_label：单独为每个模型列设置二级表头，仅展示模型编号(1)/(2)。
3. modify_spanner_header：可选操作，通过加粗一级表头让多级结构更清晰。

生成的表格会呈现你需要的两级表头结构：年份作为一级标题，模型编号作为二级标题，同一年份的模型列自动归组。

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