[R] Post Stratification

[Thomas Lumley]

On Sun, 18 Jun 2006, Mark Hempelmann wrote:

> Dear WizaRds,
>
> 	having met some of you in person in Vienna,  I think even more fondly
> of this community and hope to continue on this route. It was great
> talking with you and learning from you. Thank you. I am trying to work
> through an artificial example in post stratification. This is my dataset:
>
> library(survey)
> age <- data.frame(id=1:8, stratum=rep( c("S1","S2"),c(5,3)),
> weight=rep(c(3,4),c(5,3)), nh=rep(c(5,3),c(5,3)),
> Nh=rep(c(15,12),c(5,3)), y=c(23,25,27,21,22, 77,72,74) )
>
> pop.types <- table(stratum=age$stratum)

No, you need population totals here, that's the whole point of 
post-stratification.

> age.post <- svydesign(ids=~1, strata=NULL, data=age, fpc=~Nh)

Since the sampling isn't stratified you can't give stratum-specific fpc. 
This should give an error message. I will add one.

> post <- postStratify(design=age.post, strata=~stratum, population=pop.types)

Yes, but pop.types is defined wrongly

>
> pop.types <- data.frame(stratum = c("S1","S2"), Freq = c(15, 12))

This is the correct definition



> However, I compute:
> Nh=c(15,12); nh=c(5,3); sh=by(age$y, age$stratum, var); N=sum(Nh)
> # Mean estimator
> y.bar=by(age$y, age$stratum, mean) ## 23.6; 74.33
> estimator=1/N*sum(Nh*y.bar) ## 46.14815
> # Variance estimator
> vari=1/N^2*sum(Nh*(Nh-nh)*sh/nh)
> sqrt(vari) ##	.7425903
>
> and with Taylor expansion .7750118
>

This is partly the problem with giving stratum-specific fpc and partly 
because "survey" is using a less elegant but more general estimator (as is 
often the case).  The estimator is the one due to Valliant, described in 
the reference on the help page.

The primary difference from the estimator you given is that R does not try 
to compute finite population corrections for each post-stratum.  Another 
more subtle difference is that R's formula is estimating the conditional 
variance, conditional on the estimated weights, rather than the 
unconditional variance.


So, to do the computations correctly

> library(survey)
> age <- data.frame(id=1:8, stratum=rep( c("S1","S2"),c(5,3)),
+ weight=rep(c(3,4),c(5,3)), nh=rep(c(5,3),c(5,3)),
+ Nh=rep(c(15,12),c(5,3)), y=c(23,25,27,21,22, 77,72,74) )
>
> age$N<-rep(27,8)
>
> pop.types <- data.frame(stratum = c("S1","S2"), Freq = c(15, 12))
>

The sampling is unstratified from a population of size 27=15+12

> age.design <- svydesign(ids=~1, data=age, fpc=~N)
> age.post<-postStratify(age.design, strata=~stratum, population=pop.types)
>
> svymean(~y, age.design)
     mean     SE
y 42.625 7.8163
> svymean(~y, age.post)
     mean     SE
y 46.148 0.6737
>

Now, post-stratification adjusts the weights to sum to 15 and 12 for the 
two post-strata, so the mean estimate is a weighted average of the two 
post-stratum-specific means, weighted by 15/27 and 12/27. This is what you 
found.

The variance is the variance of the mean of the residuals after 
subtracting the post-stratum-specific mean.  That is, we can get the 
post-stratified estimates by

Adjust the weights
> age.design2 <- svydesign(ids=~1, data=age, fpc=~N,weight=~weight)
Create a 'residuals' variable
> age.design2<-update(age.design2, yres=y-ifelse(stratum=="S1", 
23.6,74.33333333))

Post-stratified mean is 46.148
> svymean(~y, age.design2)
     mean     SE
y 46.148 8.2317

Standard error of post-stratified mean is 0.6737
> svymean(~yres, age.design2)
            mean     SE
yres 1.4815e-09 0.6737


The formula for the variance of the estimated total is

> (n/(n-1))*((N-sum(nh))/N)*sum(yr2*(Nh/nh)^2)
[1] 330.915

which is formula (11) of the Valliant reference. The variance of the
mean is the same thing divided by N^2:

> (n/(n-1))*((N-sum(nh))/N)*sum(yr2*(Nh/nh)^2)/(N*N)
[1] 0.45393

where yr2<-by(yres^2,age$stratum,sum).  These agree exactly with R's 
estimates
> vcov(svytotal(~y, design=age.post))
         y
y 330.915
> vcov(svymean(~y, design=age.post))
         y
y 0.45393


 	-thomas



Thomas Lumley			Assoc. Professor, Biostatistics
tlumley at u.washington.edu	University of Washington, Seattle