在 R 调查包中使用"adjust"孤独 PSU 选项:当其他地层中的数据发生变化时,为什么单个 PSU 地层方差不会改变?



我有来自分层简单随机抽样设计的调查数据,其中一些地层只包含单个抽样单位(即使地层人口规模可能大于1)。在R调查包(http://r-survey.r-forge.r-project.org/survey/exmample-lonely.html)中,这些被称为"孤独的psu"。有几个选项可以处理这种情况,我感兴趣的是"调整"选项。

options(survey.lonely.psu="adjust")的文档说明

"单psu地层的数据集中在样本grand处而不是阶层的平均数。"

在实验这个选项时,我曾期望如果我改变另一个层的数据,我的单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

现在如果我像这样改变strata 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

我是否误解了文档的含义,或者这可能是一个bug?

仅供参考,这是我的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层(其中包含一个lonley PSU)的方差,而在估计平均值时却没有。

分析一个复杂设计的子集基本上等同于分析将子集外的观测值的抽样权值设置为零。这意味着这些观测值不会进入"大均值"。

使用TSL方差公式撤回您的结果,我们有地层2方差= 21.213203^2 = 450。除以抽样分数,得到s²= 900。900是30^2。除非我遗漏了什么,否则这个包似乎假设子群总体均值是0。这对线性化的平均值是正确的,但对总数不是。

我相信这些设置只有在单例行以某种方式与其他记录组合时才真正开始生效?因此,隔离单例记录的svyby()不会捕获其他记录的方差,但是包括单例记录和其他记录的svymean()将根据设置表现不同。注意每个svymean调用的SE不同:

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)

简短回答

TL;DR:当使用svytotal()时,调查包目前有一个"调整"选项的错误,但svymean()工作正常。这篇博文将详细介绍:https://www.practicalsignificance.com/posts/bugs-with-singleton-strata/

长回答

目前,使用"调整"选项,调查包只是将孤独的PSU集中在0周围,然后将平方添加到估计的协方差矩阵中。对于svymean()svyratio()和许多其他函数来说,这是完全合理的,但对于svytotal()来说,这是一个错误。

为什么?

调查包如何使用影响函数估计方差的背景信息

R包使用基于影响函数的线性化方法估计方差。这篇博文解释了这是什么意思,并提供了一些参考:

https://www.practicalsignificance.com/posts/survey-covariances-using-influence-functions/这是如何工作的

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

survey.lonely.psu = 'adjust'应该做什么

"调整"的文档Option说它通过"recentering"孤独的电源模块。如果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/

最新更新