调查包中的分层群集抽样估计

Question

我希望能够从分层抽样设计中估计均值和总数，其中在每个分层中使用了单阶段群集抽样。我相信我已经正确使用“survey”包的“svydesign（）”函数来规定了设计方案。但我不确定如何正确指定分层权重。

下面是示例代码。我使用“weights=”参数提供未经调整的分层权重。我预期“svytotal（）”的估计值和SE将等于分层权重之和（例如，在本示例中为70）乘以“svymean（）”的估计值和SE。但实际上，估计值的差异因子为530（这是计数数据中所有元素的分层权重之和），而SE的差异因子为898（？？？）。我的问题是：（1）我如何向“svydesign（）”提供我的三个分层权重，使其理解？（2）为什么“svytotal（）”和“svymean（）”的估计值和SE没有按同样的因子进行差异？

library(survey)

# example data from a stratified sampling design in which
# single stage cluster sampling is used in each stratum
counts <- data.frame(
  Stratum=rep(c("A", "B", "C"), c(5, 8, 8)), 
  Cluster=rep(1:8, c(3, 2, 3, 2, 3, 2, 3, 3)),
  Element=c(1, 2, 3, 1, 2, 1, 2, 3, 1, 2, 1, 2, 3, 1, 2, 1, 2, 3, 1, 2, 3),
  Count = 1:21
)
# stratum weights
weights <- data.frame(
  Stratum=c("A", "B", "C"),
  W=c(10, 20, 40)
)

# combine counts and weights
both <- merge(counts, weights)

# estimate mean and total count
D <- svydesign(id=~Cluster, strata=~Stratum, weights=~W, data=both)
a <- svymean(~Count, D)
b <- svytotal(~Count, D)

sum(weights$W)  #  70
sum(both$W)     # 530
coef(b)/coef(a) # 530 
SE(b)/SE(a)     # 898.4308

---

第一次更新

我正在添加一个图表来帮助解释我的设计。整个人口是一个已知面积的湖泊（例如，70公顷）。层有已知的面积（10、20和40公顷）。每层分配的簇数不成比例。此外，与可能被抽样的数量相比，这些簇非常小，因此有限人口校正为FPC = 1。

我想计算每单位面积的总体平均值和SE，并且总数等于70倍这个平均值和SE。

---

第二次更新

我编写了代码来从头开始进行计算。 我得到了一个总估计值为920，标准误为61.6。

library(survey)
library(tidyverse)

# example data from a stratified sampling design in which
# single stage cluster sampling is used in each stratum
counts <- data.frame(
  Stratum=rep(c("A", "B", "C"), c(5, 8, 8)),
  Cluster=rep(1:8, c(3, 2, 3, 2, 3, 2, 3, 3)),
  Element=c(1, 2, 3, 1, 2, 1, 2, 3, 1, 2, 1, 2, 3, 1, 2, 1, 2, 3, 1, 2, 3),
  Count = c(5:1, 6:21)
)
# stratum weights
areas <- data.frame(
  Stratum=c("A", "B", "C"),
  A_h=c(10, 20, 40)
)

# calculate cluster means
step1 <- counts %>%
  group_by(Stratum, Cluster) %>%
  summarise(P_hi = sum(Count), m_hi=n())
step2 <- step1 %>%
  group_by(Stratum) %>%
  summarise(
    ybar_h = sum(P_hi) / sum(m_hi),
    n_h = n(),
    sh.numerator = sum((P_hi - ybar_h*m_hi)^2),
    mbar_h = mean(m_hi)
  ) %>%
  mutate(
    S_ybar_h = 1 / mbar_h * sqrt( sh.numerator / (n_h * (n_h-1)) )
  )

# now expand up to strata
step3 <- step2 %>%
  left_join(areas) %>%
  mutate(
    W_h = A_h / sum(A_h)
  ) %>%
  summarise(
    A = sum(A_h),
    ybar_strat = sum(W_h * ybar_h),
    S_ybar_strat = sum(W_h * S_ybar_h / sqrt(n_h))
  ) %>%
  mutate(
    tot = A * ybar_strat,
    S_tot = A * S_ybar_strat
  )

step2
step3

这将产生以下输出：

> step2
# A tibble: 3 x 6
  Stratum ybar_h   n_h sh.numerator   mbar_h S_ybar_h
   <fctr>  <dbl> <int>        <dbl>    <dbl>    <dbl>
1       A    3.0     2         18.0 2.500000 1.200000
2       B    9.5     3        112.5 2.666667 1.623798
3       C   17.5     3         94.5 2.666667 1.488235
> step3
# A tibble: 1 x 5
      A ybar_strat S_ybar_strat   tot   S_tot
  <dbl>      <dbl>        <dbl> <dbl>   <dbl>
1    70   13.14286    0.8800657   920 61.6046

- Jean V. Adams

Answer 1

（对修订后问题的修订答案）

在这种情况下，svytotal 不是您想要的——它是用于被抽样元素的实际人口总数，因此当将人口视为无限大时就没有意义了。整个调查包实际上是为离散、有限的人群设计的，但我们可以绕过它。

我认为您想要获得每个分层的平均值，然后乘以分层权重。要做到这一点，

D <- svydesign(id=~Cluster, strata=~Stratum, data=both)
means<- svyby(~Count, ~Stratum, svymean, design=D)
svycontrast(means, quote(10*A+20*B+40*C))

你会收到一个警告

Warning message:
In vcov.svyby(stat) : Only diagonal elements of vcov() available

这是因为svyby不会返回层次均值之间的协方差。这并不会造成影响，因为各层次实际上是独立样本（这就是分层的含义），所以协方差为零。

- Thomas Lumley

---

这很有帮助。我决定自己计算一遍（请参见我的帖子的第二次更新）。我得到了相同的层平均值和相同的总估计，但标准误差略有不同（你是68.9，我是61.6）。我猜想这可能与你使用而我没有使用的有限人口校正有关。 - Jean V. Adams

Answer 2

svytotal的作用是根据抽样概率来计算权重，因此权重只针对抽样单位进行定义。在svydesign函数中，将这些权重应用于群集和元素（因为是集群抽样），导致总数高达530倍。你需要提供观测权重或足够的信息让svydesign自行计算。如果这是无子抽样的集群抽样，则可以将分层权重分配给各个群集以获得群集权重，然后将其分配到群集内的元素以获得观测权重。或者，如果分层权重是人口中群集数量，则可以使用svydesign函数的fpc参数。

标准误不与点估计按相同比例缩放是因为人口大小未知且必须进行估计。均值是估计的总数除以估计的人口大小，而标准误估计考虑了分母的方差及其与分子的协方差。

- Thomas Lumley

---

我理解你在说权重方面的东西。我需要根据比例计算每个元素的权重（我认为这就是你所说的“观察权重”）。我曾认为 survey 包中可能有内置的方法来实现这一点。关于标准误差的解释既揭示了问题，又令人困惑。在这种情况下，我知道放大因子是什么。它就是一个湖泊的面积。我添加了一些更详细的描述和图表来帮助解释这种情况。 - Jean V. Adams