r - 来自调查包的分层整群抽样估计

Question

我想从分层抽样设计中估计均值和总数，其中每个层都使用单阶段整群抽样。我相信我已经使用 svydesign() 正确指定了设计的功能调查包裹。但我不确定如何正确指定层权重。

示例代码如下所示。我使用 weights= 提供未调整的层权重争论。我预计估计和 SE 来自 svytotal()将等于层权重的总和(在示例中为 70)乘以估计值和来自 svymean() 的 SE .相反，估计值相差 530 倍(这是计数数据中所有元素的层权重之和)，而 SE 相差 898 倍(???)。我的问题是 (1) 如何将我的 3 个层权重提供给 svydesign()以它理解的方式，以及 (2) 为什么不是来自 svytotal() 的估计值和 SE和 svymean()因相同的因素而不同？

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，以及等于该平均值和 SE 的 70 倍的总数。



第二次更新

我编写了代码来从头开始计算。我用 se 61.6 得到了 920 的总估计值。

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

Answer 1

svytotal正在做我认为应该在这里做的事情:权重基于采样概率，因此它们仅针对采样单位定义。 svydesign call 将这些权重应用于集群和(因为集群采样)元素，使总数高出 530 倍。您需要为 svydesign 提供观察权重或足够的信息自己计算它们。如果这是没有子抽样的聚类抽样，您可以将层权重划分到聚类上以获得聚类权重，并将其划分到聚类内的元素上以获得观察权重。或者，如果层权重是总体中的聚类数，则可以使用 fpc论据 svydesign
SE 的缩放方式与点估计不同，这是因为人口规模未知，必须进行估计。均值是估计总数除以估计人口规模，SE 估计考虑了分母的方差及其与分子的协方差。