R中的频率加权，将结果与Stata进行比较

Question

我正在尝试分析来自明尼苏达大学 IPUMS 数据集的数据，以获取 R 中的 1990 US census。我正在使用survey 包，因为数据是weighted。仅获取家庭数据（并忽略人员变量以保持简单），我正在尝试计算 hhincome (household income) 的平均值。为此，我使用svydesign() 函数和以下代码创建了一个调查设计对象：

> require(foreign)
> ipums.household <- read.dta("/path/to/stata_export.dta")
> ipums.household[ipums.household$hhincome==9999999, "hhincome"] <- NA # Fix missing 
> ipums.hh.design <- svydesign(id=~1, weights=~hhwt, data=ipums.household)
> svymean(ipums.household$hhincome, ipums.hh.design, na.rm=TRUE)
      mean     SE
[1,] 37029 17.365

到目前为止一切顺利。但是，如果我在 Stata（使用 code meant for a different portion of the same dataset）中尝试相同的计算，则会得到不同的标准错误：

use "C:\I\Hate\Backslashes\stata_export.dta"
replace hhincome = . if hhincome == 9999999
(933734 real changes made, 933734 to missing)

mean hhincome [fweight = hhwt] # The code from the link above.

Mean estimation                     Number of obs    = 91746420

--------------------------------------------------------------
             |       Mean   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
    hhincome |   37028.99   3.542749      37022.05    37035.94
--------------------------------------------------------------

而且，看看另一种给这只猫剥皮的方法，survey 的作者，this suggestion 用于频率加权：

expanded.data<-as.data.frame(lapply(compressed.data,
               function(x) rep(x,compressed.data$weights)))

但是，我似乎无法让这段代码工作：

> hh.dataframe <- data.frame(ipums.household$hhincome, ipums.household$hhwt)
> expanded.hh.dataframe <- as.data.frame(lapply(hh.dataframe, function(x) rep(x, hh.dataframe$hhwt)))
Error in rep(x, hh.dataframe$hhwt) : invalid 'times' argument

我似乎无法修复。这可能与this issue有关。

总之：

1. 为什么我在 Stata 和 R 中得不到相同的答案？
2. 哪一个是对的（或者我在这两种情况下都做错了什么）？
3. 假设我的 rep() 解决方案有效，那会复制 Stata 的结果吗？
4. 正确的方法是什么？如果答案允许我使用 plyr 包进行任意计算，而不是仅限于在 survey （svymean()、svyglm() 等）中实现的函数，则表示敬意

更新

因此，在我在这里和通过电子邮件从 IPUMS 收到的出色帮助之后，我正在使用以下代码来正确处理调查权重。我在这里描述以防其他人将来遇到此问题。

初始状态准备

由于 IPUMS 目前没有发布用于将数据导入 R 的脚本，因此您需要从 Stata、@987654351 开始@ 或 SPSS。我现在会坚持使用 Stata。首先从 IPUMS 运行导入脚本。然后在继续添加以下变量之前：

generate strata = statefip*100000 + puma

这将为 240001 形式的每个 PUMA 创建一个唯一整数，前两位数字作为州 fip 代码（在马里兰州为 24），后四位 PUMA id 在每个国家基础。如果您打算使用 R，您可能还会发现运行它也很有帮助

generate statefip_num = statefip * 1

这将创建一个没有标签的附加变量，因为将 .dta 文件导入 R 会应用标签并丢失基础整数。

Stata 和svyset

正如 Keith 所解释的，调查抽样由 Stata 通过调用 svyset 来处理。

对于我现在使用的个人层面分析：

svyset serial [pweight=perwt], strata(strata)

这将权重设置为perwt，对我们在上面创建的变量进行分层，并使用家庭serial 数字来说明聚类。如果我们使用多年，我们可能想尝试

generate double yearserial = year*100000000 + serial

还要考虑纵向聚类。

对于家庭层面的分析（没有年份）：

svyset serial [pweight=hhwt], strata(strata)

应该是不言自明的（尽管我认为在这种情况下连续剧实际上是多余的）。将 serial 替换为 yearserial 将考虑时间序列。

在R

假设您要导入带有上述附加 strata 变量的 .dta 文件并在单个字母处进行分析：

require(foreign)
ipums <- read.dta('/path/to/data.dta')
require(survey)
ipums.design <- svydesign(id=~serial, strata=~strata, data=ipums, weights=perwt)

或者在家庭层面：

ipums.hh.design <- svydesign(id=~serial, strata=~strata, data=ipums, weights=hhwt)

希望有人觉得这很有帮助，非常感谢来自 IPUMS 的 Dwin、Keith 和 Brandon。

Answer 1

1&2) 您从 Lumley 引述的评论写于 2001 年，早于他发表的任何带有调查包的工作，该调查包仅发布了几年。您可能在两种不同的意义上使用“权重”。 （Lumley 在他书中的早期描述了三种可能的感觉。）调查函数 svydesign 使用的是概率权重而不是频率权重。考虑到该数据集的庞大规模，这些似乎不是真正的频率权重，而是概率权重，这意味着调查包结果是正确的，而 Stata 结果是不正确的。如果您不相信，那么调查包提供了 as.svrepdesign() 函数，Lumley 的书中描述了如何从 svydesign-object 创建复制权重向量。

3) 我想是的，但正如 RMN 所说……“那是错误的。”

4) 既然是错误的 (IMO)，那就没有必要了。

Answer 2

您不应该在 Stata 中使用频率权重。这很清楚。如果 IPUMS 没有“复杂”的调查设计，您可以使用：

mean hhincome [pw = hhwt]

或者，为方便起见：

svyset [pw = hhwt]
svy: mean hhincome
svy: regress hhincome `x'

第二个选项的好处是您可以将其用于更复杂的调查设计（通过svyset 上的选项。然后您可以运行大量命令而无需一直键入 [pw...]。

Answer 3

对于无法访问 Stata 或 SAS 的人的轻微添加； （我会把它放在 cmets 但是......） SAScii 库可以使用 SAS 代码文件读入 IPUMS 下载的数据。读取数据的代码来自doc

library(SAScii)
IPUMS.file.location <- "..\\usa_00007dat\\usa_00007.dat"
IPUMS.SAS.read.in.instructions <- "..\\usa_00007dat\\usa_00007.sas"

#store the IPUMS extract as an R data frame!
IPUMS.df <- 
  read.SAScii ( 
    IPUMS.file.location , 
    IPUMS.SAS.read.in.instructions , 
    zipped = F )