在gtsummary中添加含调查设计效应的ANOVA P值报错求助

我参考gtsummary文档中add_p()和add_difference()的自定义函数章节，尝试传入自定义函数获取包含调查设计效应的ANOVA检验P值，但运行代码时出错。

运行代码 #

rm(list = ls())

library(gtsummary)
library(dplyr)
library(srvyr)
library(survey)


data <- structure(list(id = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25), 
                       strata = c(10, 20, 30, 10, 20, 20, 10, 20, 30, 30, 10, 30, 30, 20, 10, 20, 20, 20, 10, 20, 20, 30, 30, 20, 30), 
                       weight = c(10, 8, 17, 15, 9, 10, 25, 8, 8, 13, 17, 24, 12, 15, 3, 12, 16, 17, 24, 12, 3, 2, 8, 14, 4), 
                       popgroup = c("A", "B", "A", "C", "A", "C", "B", "B", "B", "A", "A", "C", "A", "B", "A", "C", "B", "A", "C", "B", "A", "B", "C", "B", "B"),
                       gender = c("Male", "Female", "Female", "Male", "Female", "Female", "Male", "Male", "Female", "Male", "Female", "Female", "Male", "Female", "Male", "Female", "Female", "Male", "Female", "Female", "Male", "Female", "Male", "Female", "Male"),
                       inc_01 = c(1500, 1200, 130, 500, 750, 2000, 10000, 1500, 1050, 400, 360, 490, 250, 400, 2500, 1300, 800, 540, 690, 520, 600, 700, 700, 600, 400), 
                       inc_02 = c(360, 450, 120, 300, 900, 560, 450, 280, 720, 360, 1000, 900, 530, 820, 640, 520, 130, 140, 150, 650, 240, 130, 200, 300, 500)), 
                  class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA, -25L))


dclus2 <- survey::svydesign(
            id = ~ id,
            strata = ~ strata,
            weights = ~ weight,
            data = data
          )


ttest_common_variance <- function(data, variable = inc_01, by = popgroup, ...) {
                                  model<-svyglm(variable ~ by, design=dclus2)
                                  reg_v <- regTermTest(model,~ by)
                                  reg_v1 <- as.numeric(reg_v$p)
                                  reg_v2 <- tibble( p.value = reg_v1 )
}


tbl_2 <-data |>
  as_survey_design(strata = strata, weights = weight) %>%
  gtsummary::tbl_svysummary(
    by = popgroup,
    type = where(is.numeric) ~ "continuous",
    statistic = list(c(inc_01, inc_02) ~ "{mean} {mean.std.error}"),
    missing = "no",
    digits = list(c(inc_01, inc_02) ~ c(4)),
    include = c(inc_01, inc_02)
  ) |>
  gtsummary::add_p(
    test       = c(inc_01, inc_02) ~ "ttest_common_variance",
    pvalue_fun = function(x) gtsummary::style_pvalue(x, digits = 3)
  ) |>
  as_tibble()

错误信息 #

There was an error in 'add_p()/add_difference()' for variable 'inc_01', p-value omitted:
Error in svyglm.survey.design(variable ~ by, design = dclus2, family = quasibinomial()): all variables must be in design= argument

错误原因是自定义函数直接调用全局环境的dclus2，且未正确处理add_p()传入的符号变量与调查设计对象的关联。修改后的代码如下：

rm(list = ls())

library(gtsummary)
library(dplyr)
library(srvyr)
library(survey)


data <- structure(list(id = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25), 
                       strata = c(10, 20, 30, 10, 20, 20, 10, 20, 30, 30, 10, 30, 30, 20, 10, 20, 20, 20, 10, 20, 20, 30, 30, 20, 30), 
                       weight = c(10, 8, 17, 15, 9, 10, 25, 8, 8, 13, 17, 24, 12, 15, 3, 12, 16, 17, 24, 12, 3, 2, 8, 14, 4), 
                       popgroup = c("A", "B", "A", "C", "A", "C", "B", "B", "B", "A", "A", "C", "A", "B", "A", "C", "B", "A", "C", "B", "A", "B", "C", "B", "B"),
                       gender = c("Male", "Female", "Female", "Male", "Female", "Female", "Male", "Male", "Female", "Male", "Female", "Female", "Male", "Female", "Male", "Female", "Female", "Male", "Female", "Female", "Male", "Female", "Male", "Female", "Male"),
                       inc_01 = c(1500, 1200, 130, 500, 750, 2000, 10000, 1500, 1050, 400, 360, 490, 250, 400, 2500, 1300, 800, 540, 690, 520, 600, 700, 700, 600, 400), 
                       inc_02 = c(360, 450, 120, 300, 900, 560, 450, 280, 720, 360, 1000, 900, 530, 820, 640, 520, 130, 140, 150, 650, 240, 130, 200, 300, 500)), 
                  class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA, -25L))


# 修改后的自定义函数
ttest_common_variance <- function(data, variable, by, ...) {
  # 动态构造模型公式
  formula <- as.formula(paste0(rlang::as_name(variable), " ~ ", rlang::as_name(by)))
  # 使用传入的调查设计对象拟合模型
  model <- svyglm(formula, design = data)
  # 计算检验P值
  reg_v <- regTermTest(model, formula(paste0("~", rlang::as_name(by))))
  tibble(p.value = as.numeric(reg_v$p))
}


tbl_2 <- data |>
  as_survey_design(strata = strata, weights = weight) %>%
  gtsummary::tbl_svysummary(
    by = popgroup,
    type = where(is.numeric) ~ "continuous",
    statistic = list(c(inc_01, inc_02) ~ "{mean} {mean.std.error}"),
    missing = "no",
    digits = list(c(inc_01, inc_02) ~ c(4)),
    include = c(inc_01, inc_02)
  ) |>
  gtsummary::add_p(
    test       = c(inc_01, inc_02) ~ ttest_common_variance,  # 直接传递函数对象
    pvalue_fun = function(x) gtsummary::style_pvalue(x, digits = 3)
  ) |>
  as_tibble()

print(tbl_2)

修改要点 #

1. 自定义函数不再依赖全局环境的dclus2，而是使用add_p()传入的data（即调查设计对象）
2. 用rlang::as_name()将符号变量转为字符串，动态构造模型公式，确保变量能被调查设计对象识别
3. add_p()中直接传递函数对象，而非字符串形式
4. 调整regTermTest的检验公式，确保正确引用分组变量

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