【问题标题】：Using the `survey` package in R to estimate covariances between overlapping domains使用 R 中的 `survey` 包来估计重叠域之间的协方差
【发布时间】：2021-09-28 17:26:00
【问题描述】：

问题：

假设我们有测量两个重叠领域（30-50 岁和 40-60 岁）的疫苗接种状态的调查数据，我们想要估计这两个群体的疫苗接种率。显然，由于重叠（40-50 岁的人），域估计值将是相关的。

我们如何使用survey 包来估计来自诸如这些重叠域的估计的协方差？

理想情况下，是否可以使用svyby() 或svybys() 来执行此操作？

示例数据：

以下是基于survey 包中内置的api 数据集的示例。这是加利福尼亚州学校的多阶段样本。变量stype 指示学校是小学、初中还是高中。变量api00 是一个数值变量，它总结了每所学校在 2000 年的标准化考试成绩。

对于这个例子，假设我们要比较两个重叠域之间api00 的平均值：（1）小学和初中； (2) 初中和高中。

# Create a survey design object ----
  library(survey)
  data(api)
  
  dclus2 <- svydesign(id=~dnum+snum,
                      fpc=~fpc1+fpc2,
                      data=apiclus2)
  
# Add indicator variables for the domains ----

  dclus2 <- transform(dclus2,
                      E_or_M = stype %in% c("E", "M"),
                      M_or_H = stype %in% c("M", "H"))
  
# Add variables for outcome by domain indicators ----
  
  dclus2 <- transform(dclus2,
                      api00_E_or_M = api00 * E_or_M,
                      api00_M_or_H = api00 * M_or_H)

# Estimate means by domain ----
  
  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

【问题讨论】：

标签： r survey

【解决方案1】：

您可以使用带有选项covmat = TRUE 的svyratio() 函数来获取这些估计值及其方差-协方差矩阵。

  estimates <- svyratio(numerator = ~ api00_E_or_M + api00_M_or_H,
                        denominator = ~ E_or_M + M_or_H,
                        design = dclus2, covmat = TRUE)

唯一的挑战是，除了您不关心的估算值之外，输出还包含您关心的估算值。在下面的输出中，我们只关心对角线元素。

> estimates$ratio
#>               E_or_M    M_or_H
#> api00_E_or_M 682.0563 1851.9722
#> api00_M_or_H 229.7408  623.8102

estimates_of_interest <- diag(estimates$ratio)
print(estimates_of_interest)
#> [1] 682.0563 623.8102

同样，您通过调用 vcov() 获得的方差-协方差矩阵包含我们不关心的比率的估计方差和协方差，如下所示。

  full_vcov <- vcov(estimates)
  
  print(full_vcov)
#>                     api00_E_or_M/E_or_M api00_M_or_H/E_or_M api00_E_or_M/M_or_H
#> api00_E_or_M/E_or_M          1047.29150            34.07083            5189.788
#> api00_M_or_H/E_or_M            34.07083          2743.16185          -22346.457
#> api00_E_or_M/M_or_H          5189.78822        -22346.45657          211988.476
#> api00_M_or_H/M_or_H           882.76208          -109.79708            5981.838
#>                     api00_M_or_H/M_or_H
#> api00_E_or_M/E_or_M            882.7621
#> api00_M_or_H/E_or_M           -109.7971
#> api00_E_or_M/M_or_H           5981.8380
#> api00_M_or_H/M_or_H            909.3903

为了得到我们真正关心的方差-协方差矩阵，我们必须像这样对其进行子集化：


  vcov_of_estimates <- full_vcov[c("api00_E_or_M/E_or_M", "api00_M_or_H/M_or_H"),
                                 c("api00_E_or_M/E_or_M", "api00_M_or_H/M_or_H")]
  
  print(vcov_of_estimates)
#>                     api00_E_or_M/E_or_M api00_M_or_H/M_or_H
#> api00_E_or_M/E_or_M           1047.2915            882.7621
#> api00_M_or_H/M_or_H            882.7621            909.3903

这篇博文解释了统计细节： https://www.practicalsignificance.com/posts/how-correlated-are-survey-estimates-from-overlapping-groups/

【讨论】：