R: subsetting survey design fails after raking

Question

all. I hope you can help me with a problem I encountered while analyzing data from a survey. I have actually done all the analysis, but after that obtained more reliable data about the population distribution of some variables, so I decided to adjust weights by means of post-stratification to compare results. I followed instructions here and all seemed to go well. I did some calculations with the new survey design object (called r0) that I obtained after calling rake() on the unweighted survey design, and those too went well.

But they I needed to do some calculations using only a subset of data, so I subset the survey design in the following manner:

r1 <- subset(r0, subset = inlf==1)

This seemed to go well too, but then I called the following, trying to obtain means of the variable "employed" (which is a factor with levels 0 and 1) for women and men (female==0 and female==1):

emp_fem <- svyby(~employed, by = ~female, r1, svymean, vartype = c("se", "ci"))

... and I got a bunch of NAs where the means and CI values should be, and NaNs where the SE values should be.

I examined my survey design object r1, and discovered that it has the same number of observations as r0, when it should have fewer (note, a call to identical (r0, r1) does return FALSE, but I can't see the difference between the two objects). I have done the same kind of subsetting with the original survey design (before post-stratification), and obtained the expected smaller number of observations. It seems that the subset function did not work after post-stratification for some reason! So I figured the NAs and NaNs are coming from the fact that NA values that should have been excluded after subsetting are still in the survey design, and of course svymean can't deal with them. Does anyone have any ideas about the cause of this problem or how I could resolve it?

I updated both R and the survey package just now, so they are not out of date.

UPDATE: Here is a minimal reproducible example. This is not my actual data, but it reproduces the error:

# Minimal dataset
emp_df <- data.frame(female = c(1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0),
                     decade = c(2, 1, 3, 3, 2, 4, 4, 4, 3, 4, 3, 1, 2, 4, 3, 2, 4, 2, 3, 3), 
                     inlf =     c(1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1),
                     employed = c(1, 0, 1, NA, 1, 0, 1, 0, 1, 1, 0, 0, NA, NA, 1, 0, 1, NA, 0, 1))

# An unweighted survey design
u0 <- svydesign(id=~1, data = emp_df)

# A calculation with unweighted data
inlf_fem_u <- svyby(~inlf, by = ~female, u0, svymean, vartype = c("se", "ci"))

# An unweighted survey design subset 
u1 <- subset(u0, subset = inlf==1)

# A calculation with unweighted data subset 
emp_fem_u <- svyby(~employed, by = ~female, u1, svymean, vartype = c("se", "ci"))

# Weight adjustment

# A dataframe for population gender distribution 
femalepop <- data.frame(female = c(0, 1),
                     Freq = nrow(emp_df) * c(0.5023, 0.4977))

# A dataframe for population age distribution by decade
decadepop <- data.frame(decade = c("1", "2", "3", "4"),
                         Freq = nrow(emp_df) * c(0.1318, 0.2890, 0.3045, 0.2747))

# Rake to obtain weights
r0 <- rake(design = u0, 
           sample.margins = list(~female, ~decade), 
           population.margins = list(femalepop, decadepop))

# A calculation with weighted data
inlf_fem <- svyby(~inlf, by = ~female, r0, svymean, vartype = c("se", "ci"))

# A weighted data subset 
r1 <- subset(r0, subset = inlf==1)

# A calculation with weighted data subset
emp_fem <- svyby(~employed, by = ~female, r1, svymean, vartype = c("se", "ci"))

Calculations using 'u0', 'u1' and 'r0' use the expected number of observations, but I can't get emp_fem because r1 still uses NA values that it should not include at all - it still has 20 observations whereas it should have 16.

Answer 1

There are two issues here. First, you need na.rm=TRUE to drop the NA values, just as you would with mean() in base R.

Second, when you subset a survey design object there is no particular reason to expect it will get smaller. Mathematically, the correct implementation of subsetting is to set the weights to zero for observations that are not in the subpopulation, which does not change the number of rows in the object.

In some situations, R can reduce the size of the object (to save memory) when subsetting. A replicate-weights design object can just be subsetted in a straightforward way. A svydesign object that has not been post-stratified, raked, or calibrated can be subsetted and have additional information added to allow standard error calculations to reconstruct the missing zero-weight rows. For svydesign objects with raking/calibration and for twophase objects, no such memory-saving is implemented.

It's not obvious whether dim should report the size of the actual object or the size of the subset. The argument for the former is that whether anything is actually removed is just a memory optimisation and so shouldn't be visible to the user. The argument for the latter is that it's confusing if dim lies to you about the size of the object. Currently it reports the size of the object.