Frequency weighting in R, comparing results with Stata

Question

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:

1. Why don't I get the same answers in Stata and R?
2. Which one is right (or am I doing something wrong in both cases)?
3. Assuming I got the rep() solution working, would that replicate Stata's results?
4. 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.)

Update

So after the excellent help I've received here and from IPUMS via email, I'm using the following code to properly handle survey weighting. I describe here in case someone else has this problem in future.

Initial Stata Preparation

Since IPUMS don't currently publish scripts for importing their data into R, you'll need to start from Stata, SAS, or SPSS. I'll stick with Stata for now. Begin by running the import script from IPUMS. Then before continuing add the following variable:

generate strata = statefip*100000 + puma

This creates a unique integer for each PUMA of the form 240001, with first two digits as the state fip code (24 in the case of Maryland) and the last four a PUMA id which is unique on a per state basis. If you're going to use R you might also find it helpful to run this as well

generate statefip_num = statefip * 1

This will create an additional variable without labels, since importing .dta files into R apply the labels and lose the underlying integers.

Stata and svyset

As Keith explained, survey sampling is handled by Stata by invoking svyset.

For an individual level analysis I now use:

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

This sets the weighting to perwt, the stratification to the variable we created above, and uses the household serial number to account for clustering. If we were using multiple years, we might want to try

generate double yearserial = year*100000000 + serial

to account for longitudinal clustering as well.

For household level analysis (without years):

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

Should be self-explanatory (though I think in this case serial is actually superfluous). Replacing serial with yearserial will take into account a time series.

Doing it in R

Assuming you're importing a .dta file with the additional strata variable explained above and analysing at the individual letter:

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

Or at the household level:

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

Hope someone finds this helpful, and thanks so much to Dwin, Keith and Brandon from IPUMS.

Stripling

in case it's of any help, here's the complete R code to analyze the 1990 & 2000 decennial censii – Fantastic

Cool thanks, that's potentially quite useful. I found the 1990 census weirdly less daunting when I first attempted to tackle it. Annoying they seem to have redone all the census tracts so I have no obvious means of creating a time series :( – Stripling

Answer 1

1&2) The comment you cited from Lumley was written in 2001 and predates any of his published work with the survey package which has only been out a few years. You are probably using "weights" in two different senses. (Lumley describes three possible senses early in his book.) The survey function svydesign is using probability weights rather than frequency weights. Seems likely that these are not really frequency weights but rather probability weights, given the massive size of that dataset, and that would mean that the survey package result is correct and the Stata result incorrect. If you are not convinced, then the survey package offers the function as.svrepdesign() with which Lumley's book describes how to create a replicate weight vector from a svydesign-object.

3) I think so, but as RMN said ..."It would be wrong."

4) Since it's wrong (IMO) it's not necessary.

Footstall

Hi and thanks very much for the reply. Interestingly, I get this from STATA – Stripling

Right ok. So Stanford is wrong :). Running mean hhincome [pweight=hhwt] in STATA gives me the same standard error (though again rounded differently). So I need to make sure what type of weights to use. I will double check with the IPUMs. Thanks very much! – Stripling

Oh and I was aware of svyrepdesign but IPUMS explicitly only support replicate sampling in their 2005 and later datasets, and I'm looking at 1990. This is mostly a case of me not knowing how surveys work. Will report back what they tell me, and thanks so much :). – Stripling

Oh and regarding 4: the other problem is how I can use other R commands that aren't part of the survey package (like my beloved ggplot2). The nice thing about STATA in this respect is that the pweight command is essentially a filter I can attach to any other command. – Stripling

If you have a data.frame and want to "replicate" observations, you can use rep and [. df[rep(1:nrow(df), each=df$replic) , ] – Footstall

Hey thanks for the alternative to Lumley's code. Is there a way to do that generically with probability weighting as well? – Stripling

I didn't intend that to be an alternative to Lumley's code, especially when you are working with datasets the size of IPUMS. I just meant that you could do small scale "replications". If you did that with the IPNUMS data, you might get 200 million records. – Footstall

Sorry, I meant an alternative to his 2001 code, not his whole package. – Stripling

Answer 2

You shouldn't be using frequency weights in Stata. That is pretty clear. If IPUMS doesn't have a "complex" survey design, you can just use:

mean hhincome [pw = hhwt]

Or, for convenience:

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

What's nice about the second option is that you can use it for more complex survey designs (via options on svyset. Then you can run lots of commands without having to typ [pw...] all the time.

Encratis

Hi Keith, thanks very much for the STATA approach. It does look considerably more generic than R in the sense that it appears to be a generic conversion factor (a bit like a filter I suppose) that can be combined with just about any other function in STATA. Sadly, it appears I'm stuck using the functions provided by survey within R, and cannot use many of the other nice packages R provides. – Stripling

Hmmm... I've now discovered that Stata won't let me use svy: with any graphing commands, which pushes me back to R. Funny how this all works. Stata does seem to be consistently faster though. – Stripling

Sometimes I just use a svy: .... command, then graph the results as a second step. If you want something really quick, try my program: code.google.com/p/kk-adofiles/source/browse/g/graphbetas.ado (browse for the help file) – Encratis

Answer 3

Slight addition for people who don't have access to Stata or SAS; (I would put this in comments but...) The library SAScii can use the SAS code file to read in the IPUMS downloaded data. The code to read in the data is from the 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 )   

Envoi

IPUMS actually allows you to download extracts as CSVs now - SAScii is unnecessary. but to get the survey object correct, use the dataset-specific examples at asdfree.com :) – Fantastic