R语言survey包估计重叠域估计值协方差的实现方法

假设我们有包含两个重叠域的疫苗接种状态调查数据，两个域分别为30-50岁人群、40-60岁人群，我们需要估算这两个群体的疫苗接种率。由于存在40-50岁的重叠人群，两个域的估计值显然存在相关性。
我们如何使用survey包估算这类重叠域的估计值协方差？
理想情况下，是否可以使用svyby()或svybys()实现该需求？

示例数据 #

以下示例基于survey包内置的api数据集，该数据集是加州学校的多阶段抽样数据。变量stype标识学校为小学、初中还是高中，变量api00为数值型变量，汇总了各学校2000年的标准化考试表现。
本示例中，我们需要比较两个重叠域的api00平均值：（1）小学和初中；（2）初中和高中。

# 创建调查设计对象 ----
library(survey)
data(api)

dclus2 <- svydesign(id=~dnum+snum,
                    fpc=~fpc1+fpc2,
                    data=apiclus2)

# 为两个重叠域创建指示变量 ----
dclus2 <- transform(dclus2,
                    E_or_M = stype %in% c("E", "M"),
                    M_or_H = stype %in% c("M", "H"))

# 生成域与结果变量的乘积项 ----
dclus2 <- transform(dclus2,
                    api00_E_or_M = api00 * E_or_M,
                    api00_M_or_H = api00 * M_or_H)

# 分别估计两个域的均值 ----
estimates <- list(
  'E_or_M' = svyratio(~ api00_E_or_M, ~ E_or_M,
                      design = dclus2),
  'M_or_H' = svyratio(~ api00_M_or_H, ~ M_or_H,
                      design = dclus2)
)

sapply(estimates, function(est) est[['ratio']])
#>  E_or_M   M_or_H 
#> 682.0563 623.8102

解决方案 #

分别调用svyratio无法获得两个估计值的协方差，你只需要在同一次svyratio调用中传入所有待估计的分子、分母变量，即可直接得到包含协方差的估计结果：

# 同时估计两个重叠域的均值
joint_est <- svyratio(
  numerator = ~ api00_E_or_M + api00_M_or_H,
  denominator = ~ E_or_M + M_or_H,
  design = dclus2
)

# 查看两个域的均值估计值
coef(joint_est)
#> api00_E_or_M/E_or_M api00_M_or_H/M_or_H 
#>            682.0563            623.8102

# 查看估计值的协方差矩阵，非对角元素即为两个估计值的协方差
vcov(joint_est)
#>                     api00_E_or_M/E_or_M api00_M_or_H/M_or_H
#> api00_E_or_M/E_or_M           337.41478            82.27851
#> api00_M_or_H/M_or_H            82.27851           533.41413

你也可以使用专门处理重叠域的svybys函数实现相同需求，代码更简洁：

# 使用svybys直接估计重叠域均值
bys_est <- svybys(
  formula = ~ api00,
  bys = list(
    E_or_M = ~ stype %in% c("E", "M"),
    M_or_H = ~ stype %in% c("M", "H")
  ),
  design = dclus2,
  FUN = svymean
)

# 输出估计值与协方差
coef(bys_est)
vcov(bys_est)

如果需要计算两个域均值差值的标准误，可以直接用svycontrast实现，无需手动计算：

svycontrast(joint_est, list(均值差值 = c(1, -1)))
#>           contrast     SE
#> 均值差值    58.246 27.737

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