如何在R中动态子集化调查设计时迭代svytable对象

我有一份调查数据，需要针对ethnicity、sex、age三个变量，同时实现变量遍历和按因子水平子集化调查设计对象，最终获取各分组下的比例列联表。

我找到了两个接近目标的实现方案，但无法将它们结合以满足需求，尝试过程中多次出现无法解决的错误。另外，srvyr包的相关方案运行速度极慢，因此仅考虑使用survey包解决问题。

可复现代码 #

library(dplyr)
library(survey)

df <- structure(list(USOC_Wave = c(2, 6, 9, 5, 11, 11, 2, 2, 4, 3, 
10, 8, 11, 11, 11, 5, 6, 10, 6, 1, 9, 4, 9, 4, 11, 12, 5, 2, 
10, 11, 7, 5, 4, 11, 6, 10, 9, 13, 3, 7, 5, 10, 8, 7, 6, 12, 
12, 1, 12, 5), ethnicity = c("White", "White", "White", "Asian", 
"Asian", "White", "White", "White", "White", "White", "White", 
"Asian", "White", "White", "White", "White", "White", "Asian", 
"White", "White", NA, "Asian", "Asian", "White", "White", "White", 
"White", "White", "White", "White", "White", "White", "Asian", 
"White", "White", "Asian", "White", "White", "White", "White", 
"White", "White", "White", "White", "White", "White", "Asian", 
"White", "White", "Mixed"), sex = c("Men", "Men", "Men", "Men", 
"Men", "Women", "Men", "Men", "Women", "Men", "Women", "Women", 
"Women", "Women", "Women", "Women", "Women", "Men", "Women", 
"Women", "Men", "Men", "Men", "Men", "Women", "Women", "Women", 
"Men", "Women", "Women", "Women", "Men", "Women", "Women", "Women", 
"Men", "Women", "Women", "Women", "Men", "Women", "Women", "Women", 
"Women", "Women", "Men", "Men", "Men", "Women", "Women"), age = c("16-29", 
"30-64", "30-64", "30-64", "30-64", "30-64", "16-29", "65+", 
"30-64", "30-64", "30-64", "30-64", "30-64", "65+", "30-64", 
"30-64", "30-64", "30-64", "65+", "30-64", "16-29", "65+", "16-29", 
"16-29", "16-29", "30-64", "30-64", "30-64", "65+", "16-29", 
"30-64", "30-64", "65+", "30-64", "30-64", "30-64", "16-29", 
"16-29", "30-64", "30-64", "30-64", "30-64", "30-64", "30-64", 
"16-29", "30-64", "65+", "65+", "30-64", "30-64"), strata = c(2902, 
3165, 3069, 2108, 3943, 2683, 2521, 3175, 3232, 3256, 42, 3401, 
2326, 2108, 701, 2074, 1, 5122, 12, 2721, 5122, 3991, 3717, 3157, 
2311, 101, 2717, 118, 2425, 2584, 2523, 2222, 2400, 2729, 2199, 
3361, 10, 2427, 2151, 2584, 2327, 2, 2750, 3297, 2363, 114, 2750, 
2574, 2843, 4121), psu = c(3804, 4330, 4138, 2215, 38089, 3365, 
3041, 4350, 4464, 4512, 156, 11187, 2651, 2216, 1672, 2147, 3, 
52063, 47, 3441, 52086, 40537, 26666, 4314, 2621, 403, 3433, 
458, 2849, 3168, 3045, 2443, 2800, 3458, 2397, 9013, 31, 2854, 
2302, 3168, 2653, 6, 3500, 4593, 2725, 447, 3499, 3148, 3686, 
46785), weight_cs = c(2.80231904983521, 0, 0.950280964374542, 
0.28423735499382, 0.300251632928848, 0.766829490661621, 2.18452429771423, 
0.680638015270233, 0.224062830209732, 2.74595475196838, 0.718028843402863, 
0.340109616518021, 2.88688373565674, 1.17885708808899, 0.620745718479156, 
1.20946884155273, 0.57785838842392, 0.305908054113388, 0.727640688419342, 
1.17930126190186, 0, 0.623862087726593, 0.372526079416275, 0, 
1.3677384853363, 2.87374138832092, 1.31425619125366, 0.462548196315765, 
1.18157768249512, 0.814507722854614, 1.21053576469421, 2.14700984954834, 
0.449016481637955, 1.151535987854, 0.790829658508301, 0.359708696603775, 
3.43058443069458, 0.309507787227631, 1.17791354656219, 1.70297181606293, 
0.741691768169403, 1.52170836925507, 0, 0.989463746547699, 1.34024882316589, 
0.842447340488434, 0.869455099105835, 0.846965670585632, 3.40495872497559, 
0.730816066265106), unpaid = c(2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 1, 
1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 1, 2, 1)), row.names = c(NA, 
-50L), class = c("tbl_df", "tbl", "data.frame"))

uos_design <- survey::svydesign(id= ~psu, strata= ~strata, survey.lonely.psu="adjust",
                              weights= ~weight_cs, data=df)

现有方案 #

方案1：单个变量的动态子集化 #

该方案能针对单个变量，按其因子水平子集化调查对象并生成比例列联表，但需要手动为每个变量重复编写代码，无法自动化遍历多个变量：

groups <- unique(df$ethnicity) # 获取唯一水平

tablefun <- function(i){svytable(~USOC_Wave+unpaid+ethnicity,
                                 design = subset(uos_design, ethnicity == i)) %>%
    prop.table(margin = 1)  
}

results <- do.call(rbind, lapply(groups, tablefun)) # 转换为数据框

方案2：遍历多个变量生成列联表 #

该方案能遍历指定的多个变量生成列联表，但没有按变量的因子水平进行子集化：

lapply(names(df[c("ethnicity","sex","age")]), function(x){
  svytable(bquote(~.(as.name(x)) +USOC_Wave + unpaid), 
           design = uos_design) %>%
    prop.table(margin = 1)
} )

我尝试过用mapply替代lapply、在外层添加循环等方法，但均未成功实现两者的结合。

解决方案 #

通过嵌套lapply结合bquote动态构建表达式，同时实现变量遍历和因子水平子集化：

# 定义需要处理的目标变量
target_vars <- c("ethnicity", "sex", "age")

# 嵌套遍历：先遍历变量，再遍历每个变量的非NA因子水平
all_results <- lapply(target_vars, function(var) {
  # 获取当前变量的非NA唯一水平
  levels <- unique(df[[var]])
  levels <- levels[!is.na(levels)]
  
  # 遍历每个水平，生成子集化的比例表
  var_results <- lapply(levels, function(level) {
    # 动态构建子集化条件
    subset_expr <- bquote(.(as.name(var)) == .(level))
    # 子集化调查设计对象
    sub_design <- subset(uos_design, eval(subset_expr))
    # 动态构建列联表公式并计算比例，最后转为数据框
    svytable(bquote(~USOC_Wave + unpaid + .(as.name(var))), design = sub_design) %>%
      prop.table(margin = 1) %>%
      as.data.frame()
  })
  
  # 合并当前变量的所有分组结果
  do.call(rbind, var_results)
})

# 给每个变量的结果命名，方便后续调用
names(all_results) <- target_vars

# 如果需要将所有变量的结果合并为一个大的数据框
combined_results <- do.call(rbind, all_results)

代码说明 #

• 外层lapply遍历三个目标变量，内层lapply遍历每个变量的非NA因子水平，实现双重遍历
• 使用bquote动态构建子集条件和列联表公式，解决变量名的动态引用问题
• 最终结果可以按变量单独存储，也可以合并成一个完整的数据框，满足不同的后续分析需求

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