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

Question

我试图从明尼苏达大学IPUMS数据集的数据分析1990年美国人口普查在R.我正在使用该survey包,因为数据是加权的.只是拿家庭数据(并忽略人员变量以保持简单),我试图计算hhincome(家庭收入)的平均值.为此,我使用具有以下代码的函数创建了一个调查设计对象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(使用代码表示同一数据集的不同部分),我会得到不同的标准错误:

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,有这个频率加权的建议:

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

我似乎无法修复.这可能与此问题有关.

总而言之:

1. 为什么我没有拿到在相同的答案Stata和R？
2. 哪一个是正确的(或者我在两种情况下做错了什么)？
3. 假设我得到了rep()解决方案,那会复制Stata结果吗？
4. 什么是正确的方法呢？荣誉如果答案让我用plyr包做任意的计算,而不是仅限于实现的功能survey(svymean(),svyglm()等等)

更新

所以在我通过电子邮件收到我和IPUMS提供的出色帮助后,我正在使用以下代码来正确处理调查权重.我在这里描述以防其他人将来遇到这个问题.

最初的Stata准备

由于IPUMS目前没有公布脚本导入其数据导入R,您需要从启动Stata,SAS或SPSS.我Stata现在坚持.首先从IPUMS运行导入脚本.然后在继续添加以下变量之前:

generate strata = statefip*100000 + puma

这为每个PUMA表单240001 创建一个唯一的整数,前两位数字作为状态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

假设您使用上面解释.dta的附加strata变量导入文件并在单个字母处进行分析:

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 ...].