R中的频率加权，与STATA比较结果 [英] Frequency weighting in R, comparing results with Stata
问题描述
我正在尝试分析来自明尼苏达大学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:IHateBackslashesstata_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使用了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有关。
总而言之：
- 为什么我在
Stata和R中没有得到相同的答案？ - 哪一个是正确的(或者我在两种情况下都做错了什么)？
- 假设我使
rep()解决方案起作用，该解决方案是否会复制Stata的结果？ - 正确的方法是什么？如果答案允许我使用
plyr包进行任意计算，而不是仅限于survey(svymean()，svyglm()等)中实现的函数，则值得称赞。)
更新
因此，在我从这里和通过电子邮件从IPUMS获得了极好的帮助之后，我使用以下代码来正确处理调查权重。我在此进行说明，以备将来其他人遇到此问题。初始状态准备
由于IPUM当前不发布将其数据导入R的脚本，因此您需要从Stata、SAS或SPSS开始。目前，我将继续使用Stata。首先从IPUMS运行导入脚本。然后在继续之前添加以下变量：
generate strata = statefip*100000 + puma
PUMA创建一个唯一的整数，前两位是州FIP码(在马里兰州为24)，后四位是在每个州唯一的id。如果您要使用R，您可能也会发现运行此命令也很有帮助
generate statefip_num = statefip * 1
这将创建一个不带标签的附加变量，因为将.dta文件导入R会应用标签并丢失基础整数。
状态和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。
推荐答案
1&；2)您从Lumley引用的评论写于2001年，早于他发表的任何关于调查包的工作，而调查包只发布了几年。你可能在两种不同的意义上使用"权重"。(Lumley在书的开头描述了三种可能的感觉。)调查函数svyDesign使用概率权重而不是频率权重。考虑到数据集的巨大规模，这些似乎不是真正的频率权重，而是概率权重，这将意味着调查包结果是正确的，而STATA结果是错误的。如果您不信服，那么调查包提供了函数.svepDesign()，Lumley的书使用该函数描述了如何从svyDesign-对象创建复制权重向量。
3)我想是的，但正如RMN所说……"这是错误的。"
4)因为这是错误的(国际海事组织)，所以没有必要。
这篇关于R中的频率加权，与STATA比较结果的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
