我正在尝试分析明尼苏达大学IPUMS数据集中1990 US census中1990 US census的数据 . 我正在使用survey包，因为数据是weighted . 只需获取家庭数据（并忽略人员变量以保持简单），我试图计算 hhincome （household income）的平均值 . 为此，我使用svydesign()函数创建了 survey design object ，其代码如下：
> 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有关 .
总而言之：
-
为什么我在 Stata 和 R 得到相同的答案？
-
哪一个是正确的（或者我在两种情况下都做错了什么）？
-
假设我的
rep()解决方案有效，那会复制 Stata 的结果吗？ -
right 这样做的方法是什么？如果答案允许我使用plyr包进行任意计算，而不是仅限于
survey（svymean()，svyglm()等）中实现的功能，那就值得称赞
更新
所以在我通过电子邮件收到我和IPUMS提供的出色帮助后，我正在使用以下代码来正确处理调查权重 . 我在这里描述以防其他人将来遇到这个问题 .
初始Stata准备
由于IPUMS目前不发布用于将其数据导入 R 的脚本，因此您需要从 Stata ， SAS 或 SPSS 开始 . 我现在坚持使用 Stata . 首先从IPUMS运行导入脚本 . 然后在继续添加以下变量之前：
generate strata = statefip*100000 + puma
这为240001形式的每个 PUMA 创建一个唯一的整数，前两个数字作为状态fip代码（在Maryland的情况下为24），后四个为 PUMA id，在每个状态的基础上是唯一的 . 如果你打算使用 R ，你也可能会发现运行它也很有帮助
generate statefip_num = statefip * 1
这将创建一个没有标签的附加变量，因为将 .dta 文件导入 R 会应用标签并丢失基础整数 .
Stata和svyset
正如基思解释的那样，调查抽样由 Stata 通过调用 svyset 来处理 .
对于个人级别分析，我现在使用：
svyset serial [pweight=perwt], strata(strata)
这会将权重设置为 perwt ，即我们在上面创建的变量的分层，并使用家庭 serial 数来计算群集 . 如果我们使用多年，我们可能想尝试
generate double yearserial = year*100000000 + serial
也考虑纵向聚类 .
家庭层面分析（无年）：
svyset serial [pweight=hhwt], strata(strata)
应该是不言自明的（虽然我认为在这种情况下，串行实际上是多余的） . 用 yearserial 替换 serial 将考虑时间序列 .
在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 .
3 回答
1＆2）你引用Lumley的评论是在2001年写的，早于他发表的任何一篇研究报告都是在几年之后才发表的 . 你可能在两种不同的意义上使用“权重” . （Lumley在他的书的早期描述了三种可能的感觉 . ）调查函数svydesign使用概率权重而不是频率权重 . 考虑到该数据集的大小，这似乎不是真正的频率权重而是概率权重，这意味着调查包结果是正确的并且Stata结果不正确 . 如果你不相信，那么调查包提供了函数as.svrepdesign（），Lumley的书描述了如何从svydesign-object创建一个复制权重向量 .
3）我是这么认为的，但正如RMN所说的那样......“这是错误的 . ”
4）因为它是错的（IMO）所以没有必要 .
您不应该在Stata中使用频率权重 . 这很清楚 . 如果IPUMS没有“复杂”的调查设计，您可以使用：
或者，为方便起见：
什么是关于第二个选项的好处是你可以将它用于更复杂的测量设计（通过svyset上的选项 . 然后你可以运行大量命令而无需一直打字[pw ...] .
对无法访问Stata或SAS的人员进行轻微添加; （我会把它放在评论中但是......）库SAScii可以使用SAS代码文件读取IPUMS下载的数据 . 读入数据的代码来自doc