基于svydesign与tbl_svysummary构建分层NHANES Table 1的技术求助

我需要为NHANES调查数据构建Table 1，要求先按肥胖/非肥胖二分类变量分层，再按**对照组/治疗组（wlp_yn）**二分类变量二次分层。需输出：

• 分类特征的频数(行百分比)
• 连续基线特征的均值(标准误SE)
• 对应组间比较的p值

尝试过tbl_svysummary()、svyby()、tbl_strata()和CreateTableOne()但未成功，以下是尝试的示例代码：

add_p_svysummary_ex1 <-
  obese_adults %>%
  tbl_svysummary(by = wlp_yn, percent = "row", include = c(age_group, RIDAGEYR), 
statistic = list(all_continuous() ~ "{mean} ({sd})")) %>%
  add_p()
add_p_svysummary_ex1


svyby(~RIDAGEYR, ~age_group+wlp_yn, obese_adults, svymean) # avg age of each age group

附可复现的简化代码：

# DEMO
demo <- nhanes('DEMO')
demo_vars <- names(demo)
demo2 <- nhanesTranslate('DEMO', demo_vars, data = demo)

# PRESCRIPTION MEDICATIONS
rxq_rx <- nhanes('RXQ_RX')
rxq_rx_vars <- names(rxq_rx)
rxq_rx2 <- nhanesTranslate('RXQ_RX', rxq_rx_vars, data = rxq_rx)
rxq_rx2 <- rxq_rx2 %>% select("SEQN", "RXD240B") %>% filter(!is.na(RXD240B)) %>% group_by(SEQN) %>% dplyr::summarise(across(everything(), ~toString(na.omit(.))))

nhanesAnalysis = join_all(list(demo2, rxq_rx2), by = "SEQN", type = "full")

# Reconstructing survey weights for combining 1999-2018 - Combining ten survey cycles (twenty years) 
nhanesAnalysis$wtint20yr <- ifelse(nhanesAnalysis$SDDSRVYR %in% c(1,2), (2/10 * nhanesAnalysis$WTINT4YR), # for 1999-2002
                                          (1/10 * nhanesAnalysis$WTINT2YR)) # for 2003-2018


# sample weights 
nhanesDesign <- svydesign(id      = ~SDMVPSU,
                          strata  = ~SDMVSTRA,
                          weights = ~wtint20yr,
                          nest    = TRUE,
                          data    = nhanesAnalysis)

# subset
obese_adults  <- subset(nhanesDesign, (obesity == 1 & !is.na(BMXBMI) & RIDAGEYR >= 60))

解决方案 #

方法1：gtsummary分层统计（推荐） #

使用tbl_strata()实现第一层（肥胖/非肥胖）的分层，在每个层内调用tbl_svysummary()按wlp_yn分组统计，同时指定正确的统计量格式和组间检验方法。

注意：需先调整数据集，确保包含肥胖和非肥胖两组（原代码仅保留了肥胖人群，需修改子集条件）。

# 调整子集：保留60岁以上、BMI非缺失的人群，包含肥胖/非肥胖两组
full_adults <- subset(nhanesDesign, !is.na(BMXBMI) & RIDAGEYR >= 60)

# 构建分层Table 1
table1 <- full_adults %>%
  tbl_strata(
    strata = obesity,  # 第一层：按肥胖/非肥胖分层
    .tbl_fun = function(strata_data) {
      strata_data %>%
        tbl_svysummary(
          by = wlp_yn,  # 第二层：按对照组/治疗组分层
          include = c(age_group, RIDAGEYR,  # 替换为你需要的所有特征
                      # 例如添加性别、种族等分类变量，或其他连续变量
                      ),
          statistic = list(
            all_continuous() ~ "{mean} ({std.error})",  # 连续变量：均值(SE)
            all_categorical() ~ "{n} ({p_row}%)"        # 分类变量：频数(行百分比)
          ),
          percent = "row"
        ) %>%
        add_p(
          test = list(
            all_continuous() ~ "svy.t.test",  # 连续变量用加权t检验
            all_categorical() ~ "svy.chisq.test"  # 分类变量用加权卡方检验
          ),
          pvalue_fun = ~style_pvalue(., digits = 3)  # 格式化p值
        ) %>%
        add_overall(label = "合计")  # 可选：添加每层的合计行
    },
    .header = "**肥胖状态: {strata}**"  # 设置分层标题格式
  )

# 查看输出结果
table1

# 可选：导出为Word/HTML格式
table1 %>% as_flex_table() %>% flextable::save_as_docx(path = "NHANES_Table1.docx")

方法2：survey包手动计算统计量 #

如果gtsummary的方法不适用，可以用survey包的基础函数批量计算统计量和p值，再手动整理成表格。

1. 连续变量统计量与组间p值

# 按肥胖+治疗状态分组，计算连续变量的均值和SE
cont_stats <- svyby(
  formula = ~RIDAGEYR,  # 替换为你的连续变量
  by = ~obesity + wlp_yn,
  design = full_adults,
  FUN = function(x) {
    mean_val <- svymean(x)[[1]]
    se_val <- sqrt(vcov(svymean(x)))[[1]]
    return(c(mean = mean_val, se = se_val))
  }
)

# 计算同一肥胖分层内，治疗组vs对照组的p值
cont_p <- lapply(unique(full_adults$variables$obesity), function(obesity_level) {
  sub_design <- subset(full_adults, obesity == obesity_level)
  t_test <- svyttest(RIDAGEYR ~ wlp_yn, sub_design)
  return(data.frame(
    obesity = obesity_level,
    variable = "RIDAGEYR",
    p_value = t_test$p.value
  ))
}) %>% dplyr::bind_rows()

2. 分类变量统计量与组间p值

# 按肥胖+治疗状态分组，计算分类变量的频数和行百分比
cat_stats <- svytable(
  formula = ~age_group + obesity + wlp_yn,  # 替换为你的分类变量
  design = full_adults
) %>% as.data.frame() %>%
  dplyr::group_by(obesity, wlp_yn) %>%
  dplyr::mutate(row_pct = Freq / sum(Freq) * 100) %>%
  dplyr::ungroup()

# 计算同一肥胖分层内，治疗组vs对照组的卡方p值
cat_p <- lapply(unique(full_adults$variables$obesity), function(obesity_level) {
  sub_design <- subset(full_adults, obesity == obesity_level)
  chisq_test <- svychisq(~age_group + wlp_yn, sub_design)
  return(data.frame(
    obesity = obesity_level,
    variable = "age_group",
    p_value = chisq_test$p.value
  ))
}) %>% dplyr::bind_rows()

关键注意事项 #

• 原代码仅保留了肥胖人群，需修改子集条件为!is.na(BMXBMI) & RIDAGEYR >= 60，确保包含肥胖/非肥胖两组才能完成第一层分层。
• 组间检验需使用survey包提供的加权检验方法（svy.t.test/svy.chisq.test），不能用常规的t.test或chisq.test，否则会忽略抽样权重。
• 对于多分类变量，svy.chisq.test给出的是整体组间差异的p值，若需两两比较需额外处理。

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