在R调查包中使用“adjust”选项来处理孤立的电源单元：为什么当其他分层数据发生变化时，单个电源单元分层方差不会改变？

Question

我有一些调查数据，采用分层简单随机抽样设计，其中一些层包含仅一个抽样单位（即使层人口数量可能>1）。在R调查包中，这些被称为“孤独的PSU” (http://r-survey.r-forge.r-project.org/survey/exmample-lonely.html)。有几种处理这种情况的选项，我感兴趣的是“调整”选项。 options(survey.lonely.psu="adjust") 的文档说明：

“单个PSU层的数据以样本总均值为中心，而不是层均值。”

在尝试此选项时，我预计如果我更改另一个层中的数据，则单个PSU层的方差将发生变化，但事实并非如此。以下是一个小例子：

library(survey)
options(survey.lonely.psu="adjust")

# sample 1
dat1 <- data.frame(N = c(3, 3, 2), h = c(1, 1, 2), y = c(2, 6, 15))
survey1 <- svydesign(~1, fpc = ~N, strata = ~h, data = dat1)
svyby(~y, by = ~h, design = survey1, FUN = svytotal)

请注意结果中分层2的标准误差，这是单一PSU层：

  h  y        se
1 1 12  3.464102
2 2 30 21.213203

现在，如果我像这样更改分层1中的数据

# sample 2
(dat2 <- data.frame(N = c(3, 3, 2), h = c(1, 1, 2), y = c(200, 600, 15)))
(survey2 <- svydesign(~1, fpc = ~N, strata = ~h, data = dat2))
svyby(~y, by = ~h, design = survey2, FUN = svytotal)

对于层级1，结果按预期发生变化，但是层级2的标准误差仍然相同。

  h    y       se
1 1 1200 346.4102
2 2   30  21.2132

我是否误解了文档的意思，还是这可能是一个错误？

顺便说一下，这是我的sessionInfo:

R version 3.1.3 (2015-03-09)
Platform: i386-w64-mingw32/i386 (32-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1
other attached packages:
[1] survey_3.30-3

编辑

我对收到的初始答案的理解是，当使用svyby函数对数据进行子集处理时，方差调整不会生效。然而，当我比较分层总体的总数方差与总体方差时，似乎就像没有孤立PSU时一样，总体方差只是独立抽样分层方差的方差：

> vcov(svyby(~y, by = ~h, design = survey1, FUN = svytotal))
   1   2
1 12   0
2  0 450

> vcov(svytotal(~y, survey1))
    y
y 462

似乎如果在合并所有数据时进行某种中心化到总平均值的处理，后面的方差应该是不同的。

作为相关问题，这里比较了使用svyby估算均值和总数时的情况：

> svyby(~y, by = ~h, design = survey1, FUN = svymean)
  h  y    se
1 1  4 1.155
2 2 15 0.000

> svyby(~y, by = ~h, design = survey1, FUN = svytotal)
  h  y     se
1 1 12  3.464
2 2 30 21.213

我对在估计总体时为什么会对层2（其中包含一个孤立的PSU）进行方差估计感到困惑，但在估计平均值时却不这样做感到困惑。

- Bryan

Answer 1

对复杂设计的一个子集进行分析基本上相当于将采样权重设为零，以便于不在子集之外的观察结果不会被计算入“总均值”中。

- Thomas Lumley

---

谢谢Thomas。也许使用svyby不是理解adjust选项的最佳方法。尽管如此，我编辑了我的问题，以显示分层方差之和（其中第2层有一个孤立的PSU）等于总体方差。如果进行了相对于总平均值的居中处理，它们似乎不应该是等价的。此外，当使用svyby时，我很困惑为什么第2层均值没有方差，但总体有方差。 - Bryan

我认为svytotal()函数中的"adjust"选项存在一个错误，但是svymean()函数没有。我在这里尝试解释一下：https://dev59.com/5o_ea4cB1Zd3GeqPIwjA#73588156 - bschneidr

Answer 2

使用TSL方差公式撤销您的结果后，我们得到层2方差= 21.213203 ^ 2 = 450。除以抽样分数，您将获得s ^ 2 = 900。 900是30 ^ 2。除非我漏掉了什么，否则该软件包似乎假定子组人口平均值为0。这对于线性化平均值是正确的，但不适用于总数。

- Jules

---

这是正确的。调查包已经有一个关于“adjust”选项的错误已经有一段时间了。你是对的，它只是将孤立的PSU居中到0，然后将平方加到估计的协方差矩阵中。这是一个错误。它应该做的是将孤立的PSU居中到所有分层平均PSU的平均值。这里有一个可重现的例子，演示了这个错误： https://github.com/bschneidr/r-forge-survey-mirror/issues/5 - bschneidr

最终写了一篇关于这个 bug 的博客文章。调查包的行为对于 svymean() 是有意义的，但对于 svytotal() 来说肯定是一个 bug。影响函数的总体均值和加权样本均值对于均值、比率等都是零，但总数是一个例外。 https://www.practicalsignificance.com/posts/bugs-with-singleton-strata/ - bschneidr

Answer 3

我认为这些设置只有在单例行与其他记录某种方式组合后才真正开始生效？ 因此，隔离单例记录的 svyby() 不会捕获其他记录的方差，但包括单例记录和其他记录的 svymean() 将根据设置而表现不同。 请注意每个 svymean 调用的标准误差不同：

options(digits=22)
library(survey)
options(survey.lonely.psu="adjust")


# sample 1
dat1 <- data.frame(h = c(1, 1, 2 ), w=1:3 , y = c(2, 6, 15))
survey1 <- svydesign(~1,strata=~h, w=~w,data = dat1)
svyby(~y, ~h,survey1,svymean)
svymean(~y, survey1)

# sample 2
dat1 <- data.frame(h = c(1, 1, 2 ), w=1:3 , y = c(200, 600, 15))
survey1 <- svydesign(~1,strata=~h, w=~w,data = dat1)
svyby(~y, ~h,survey1,svymean)
svymean(~y, survey1)


options(survey.lonely.psu="average")

# sample 1
dat1 <- data.frame(h = c(1, 1, 2 ), w=1:3 , y = c(2, 6, 15))
survey1 <- svydesign(~1,strata=~h, w=~w,data = dat1)
svyby(~y, ~h,survey1,svymean)
svymean(~y, survey1)

# sample 2
dat1 <- data.frame(h = c(1, 1, 2 ), w=1:3 , y = c(200, 600, 15))
survey1 <- svydesign(~1,strata=~h, w=~w,data = dat1)
svyby(~y, ~h,survey1,svymean)
svymean(~y, survey1)


options(survey.lonely.psu="remove")

# sample 1
dat1 <- data.frame(h = c(1, 1, 2 ), w=1:3 , y = c(2, 6, 15))
survey1 <- svydesign(~1,strata=~h, w=~w,data = dat1)
svyby(~y, ~h,survey1,svymean)
svymean(~y, survey1)

# sample 2
dat1 <- data.frame(h = c(1, 1, 2 ), w=1:3 , y = c(200, 600, 15))
survey1 <- svydesign(~1,strata=~h, w=~w,data = dat1)
svyby(~y, ~h,survey1,svymean)
svymean(~y, survey1)

- Anthony Damico

---

谢谢Anthony。你的例子和我的有点不同，因为你的权重被错误地指定了（应该是w=c(3/2, 3/2, 2/1)），你的调查设计是带替换的，因为没有包括fpc，而且你看的是平均数而不是总数。尽管如此，你提出了一个很好的观点，这让我比较了svymeans和svytotals，有和没有svyby。但我仍然有点困惑，所以我会尝试编辑我的原始问题以进一步说明。 - Bryan

@Bryan 去掉权重，不同的 lonely.psu 选项结果仍然不同。Lumley博士的答案是正确的。 - Anthony Damico

Answer 4

简短回答

TL;DR: 当使用svytotal()时，调整选项在调用'adjust'时会出现错误，但对于svymean()则没有问题。本博客文章详细介绍了这个问题： https://www.practicalsignificance.com/posts/bugs-with-singleton-strata/

详细回答

目前，使用'adjust'选项时，survey包只是将孤立的PSU居中到0，然后将平方加到估计的协方差矩阵中。这对于svymean()、svyratio()和许多其他函数来说是完全合理的，但对于svytotal()来说是一个错误。

为什么呢？

关于如何使用影响函数估计方差的背景知识

R包使用基于影响函数的线性化方法来估计方差。本博客文章解释了这意味着什么，并提供了一些参考文献：

https://www.practicalsignificance.com/posts/survey-covariances-using-influence-functions/#how-does-this-work

本质上，对于感兴趣的统计量，我们通过将其表示为影响函数的加权总和并估计该加权总和的方差来计算其方差。对于诸如均值或比率之类的统计量，影响函数的加权总和始终为0。总数是一种特殊情况，其中影响函数等于原始变量，因此它们的总和通常不为零。

survey.lonely.psu = 'adjust' 的作用

“adjust”选项的文档说，它通过“重新定位”孤立PSU来工作。如果PSU不孤立，我们可以通过计算PSU总数与其层中平均PSU总数之间的平方差来得到其方差贡献的合理估计。但是，由于它是其层中唯一的PSU，因此该差异为零。理论上，使用survey.lonely.psu ='adjust'，我们改为通过计算PSU总数与所有层中平均PSU总数之间的平方差来估计其方差贡献。

调查包实际执行的操作以及这对svytotal()的问题所在

因为调查包使用基于影响函数的线性化方法，所以在所有层中PSU总平均数对于除总数以外的每个统计量都是零。我认为由于这个事实，调查包被编写成当用户指定survey.lonely.psu = 'adjust'时，孤立PSU的方差贡献仅通过平方PSU总数（即减去零）来考虑。这很聪明，适用于均值、比率等等。但它在总数方面失败了，因为其影响函数的总和具有唯一的非零属性。

https://www.practicalsignificance.com/posts/bugs-with-singleton-strata/

- bschneidr