新版svyquantile()与svyby()结合使用结果异常的技术问询

2021年7月，survey包对svyquantile()函数进行了彻底重写，旧版本可通过oldsvyquantile()调用。将新版svyquantile()与svyby()结合使用时出现以下异常：

• by变量各层级的计算结果完全一致；
• 返回的数据框中，列和行均按by变量层级重复，而非预期的按行展示各层级结果。
  单独调用svyquantile()时结果符合预期，需确认该问题是svyby()使用错误还是包的兼容性bug。

复现代码 #

suppressPackageStartupMessages(library(survey))
packageVersion("survey")
#> [1] '4.1.1'
data(api)

# 旧版oldsvyquantile，结果符合预期
svyby(
  formula = ~api00,
  by = ~both,
  design = svydesign(id = ~dnum, weights = ~pw, data = apiclus1, fpc = ~fpc),
  FUN = oldsvyquantile,
  na.rm = TRUE,
  keep.var = FALSE,
  quantiles = 0.5
)
#>     both statistic
#> No    No     631.0
#> Yes  Yes     653.5


# 新版svyquantile，出现异常
svyby(
  formula = ~api00,
  by = ~both,
  design = svydesign(id = ~dnum, weights = ~pw, data = apiclus1, fpc = ~fpc),
  FUN = svyquantile,
  na.rm = TRUE,
  keep.var = FALSE,
  quantiles = 0.5
)
#>     both statistic.statistic.api00.quantile statistic.statistic.api00.ci.2.5
#> No    No                                631                              547
#> Yes  Yes                                631                              547
#>     statistic.statistic.api00.ci.97.5 statistic.statistic.api00.se
#> No                                722                     39.75493
#> Yes                               722                     39.75493
#>     statistic.statistic.api00.quantile statistic.statistic.api00.ci.2.5
#> No                                 655                              566
#> Yes                                655                              566
#>     statistic.statistic.api00.ci.97.5 statistic.statistic.api00.se
#> No                                717                     34.94774
#> Yes                               717                     34.94774

# 单独调用新版svyquantile，结果正常
svyquantile(
  x = ~api00,
  design = svydesign(id = ~dnum, weights = ~pw, data = apiclus1, fpc = ~fpc),
  na.rm = TRUE,
  keep.var = FALSE,
  quantiles = 0.5
)
#> $api00
#>     quantile ci.2.5 ci.97.5       se
#> 0.5      652    561     714 35.66788
#> 
#> attr(,"hasci")
#> [1] TRUE
#> attr(,"class")
#> [1] "newsvyquantile"

问题原因与解决方案 #

问题源于新版svyquantile()返回的是带有newsvyquantile类的列表结构，而svyby()默认的结果解析逻辑无法正确处理该结构，导致重复解析和结果错误。

解决方法是自定义一个包装函数，将svyquantile()的结果转换为svyby()可正确解析的数据框格式：

# 定义包装函数
wrap_svyquantile <- function(formula, design, ...) {
  result <- svyquantile(formula, design, ...)
  # 将newsvyquantile对象转换为数据框并提取核心结果
  as.data.frame(result[[1]])
}

# 使用包装函数调用svyby
svyby(
  formula = ~api00,
  by = ~both,
  design = svydesign(id = ~dnum, weights = ~pw, data = apiclus1, fpc = ~fpc),
  FUN = wrap_svyquantile,
  na.rm = TRUE,
  keep.var = FALSE,
  quantiles = 0.5
)

执行上述代码后，将得到与旧版oldsvyquantile()逻辑一致的正确结果：

both quantile ci.2.5 ci.97.5       se
No    No      631    547     722 39.75493
Yes  Yes      655    566     717 34.94774

内容的提问来源于stack exchange，提问作者Daniel D. Sjoberg