如何为tbl_svysummary添加Cohen's D和Cramer's V效应量估计？

问题根源 #

你之前在tbl_summary中使用的效应量函数是针对原始非加权数据编写的，但tbl_svysummary的add_stat传入的是survey::svydesign调查设计对象，而非普通data.frame，直接调用处理原始数据的函数自然会失效。

解决方案：编写加权效应量函数 #

需要基于调查设计对象，使用survey包工具计算加权后的统计量，进而推导效应量：

1. 加权Cohen's D（连续变量） #

手动实现加权版本的Cohen's D，利用svymean和svyvar提取分组加权均值、方差，再结合加权样本量计算合并标准差：

my_weighted_cohens_d <- function(data, variable, by, ...) {
  # 计算分组加权均值与方差
  mean_by <- survey::svymean(as.formula(glue::glue("{variable} ~ {by}")), data = data)
  var_by <- survey::svyvar(as.formula(glue::glue("{variable} ~ {by}")), data = data)
  
  # 提取分组统计量
  m1 <- coef(mean_by)[[1]]
  m2 <- coef(mean_by)[[2]]
  sd1 <- sqrt(coef(var_by)[[1]])
  sd2 <- sqrt(coef(var_by)[[2]])
  
  # 计算分组加权样本量
  n1 <- survey::svytotal(~1, subset(data, data$variables[[by]] == levels(data$variables[[by]])[1]))
  n2 <- survey::svytotal(~1, subset(data, data$variables[[by]] == levels(data$variables[[by]])[2]))
  
  # 计算合并标准差与Cohen's D
  pooled_sd <- sqrt(((n1 - 1)*sd1^2 + (n2 - 1)*sd2^2)/(n1 + n2 - 2))
  d <- abs(m1 - m2)/pooled_sd
  
  round(d, 3)
}

2. 加权Cramer's V（分类变量） #

基于加权列联表计算卡方值，再推导Cramer's V：

my_weighted_cramer_v <- function(data, variable, by, ...) {
  # 生成加权列联表
  tbl <- survey::svytable(as.formula(glue::glue("{variable} ~ {by}")), data = data)
  # 计算加权卡方检验
  chisq_test <- survey::svy_chisq(tbl)
  chi2 <- chisq_test$statistic
  
  # 计算Cramer's V
  n <- sum(tbl)
  min_dim <- min(dim(tbl)) - 1
  v <- sqrt(chi2 / (n * min_dim))
  
  round(v, 3)
}

完整可运行代码 #

library(tidyverse)
library(gtsummary)
library(survey)
library(glue)

# 定义加权效应量函数
my_weighted_cohens_d <- function(data, variable, by, ...) {
  mean_by <- survey::svymean(as.formula(glue::glue("{variable} ~ {by}")), data = data)
  var_by <- survey::svyvar(as.formula(glue::glue("{variable} ~ {by}")), data = data)
  
  m1 <- coef(mean_by)[[1]]
  m2 <- coef(mean_by)[[2]]
  sd1 <- sqrt(coef(var_by)[[1]])
  sd2 <- sqrt(coef(var_by)[[2]])
  
  n1 <- survey::svytotal(~1, subset(data, data$variables[[by]] == levels(data$variables[[by]])[1]))
  n2 <- survey::svytotal(~1, subset(data, data$variables[[by]] == levels(data$variables[[by]])[2]))
  
  pooled_sd <- sqrt(((n1 - 1)*sd1^2 + (n2 - 1)*sd2^2)/(n1 + n2 - 2))
  d <- abs(m1 - m2)/pooled_sd
  
  round(d, 3)
}

my_weighted_cramer_v <- function(data, variable, by, ...) {
  tbl <- survey::svytable(as.formula(glue::glue("{variable} ~ {by}")), data = data)
  chisq_test <- survey::svy_chisq(tbl)
  chi2 <- chisq_test$statistic
  
  n <- sum(tbl)
  min_dim <- min(dim(tbl)) - 1
  v <- sqrt(chi2 / (n * min_dim))
  
  round(v, 3)
}

# 生成带效应量的汇总表
tbl_svysummary_ex1 <- 
  survey::svydesign(~1, data = as.data.frame(Titanic), weights = ~Freq) %>% 
  tbl_svysummary(by = Survived, percent = "row", include = c(Class, Age)) %>% 
  add_p(test = list(all_categorical() ~ "svy.chisq.test")) %>% 
  add_stat(
    fns = list(all_continuous() ~ my_weighted_cohens_d,
               all_categorical() ~ my_weighted_cramer_v)) %>% 
  modify_header(add_stat_1 ~ "**Effect size**")

# 查看结果
tbl_svysummary_ex1

注意事项 #

• 上述Cohen's D函数默认假设分组变量为二分类，若需支持多分类，需调整均值、方差的提取逻辑。
• 所有统计量计算均基于调查加权逻辑，符合复杂抽样数据的分析要求。

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