I'm trying to analyze data from the University of Minnesota IPUMS dataset for the 1990 US census in R. I'm using the survey package because the data is weighted. Just taking the household data (and ignoring the person variables to keep things simple), I am attempting to calculate the mean of hhincome (household income). To do this I created a survey design object using the svydesign() function with the following code:

> 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

So far so good. However, I get a different standard error if I attempt the same calculation in Stata(using 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--------------------------------------------------------------

And, looking at another way to skin this cat, the author of survey, has this suggestion for frequency weighting:

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

However, I can't seem to get this code to work:

> 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

Which I can't seem to fix. This may be related to this issue.

So in sum:

• Why don't I get the same answers in Stata and R?
• Which one is right (or am I doing something wrong in both cases)?
• Assuming I got the rep() solution working, would that replicate Stata's results?
• What's the right way to do it? Kudos if the answer allows me to use the plyr package for doing arbitrary calculations, rather than being limited to the functions implemented in survey(svymean(), svyglm() etc.)