[R] Survey and Stratification
Mark Hempelmann
neo27 at t-online.de
Thu May 26 16:10:41 CEST 2005
Dear WizaRds,
Working through sampling theory, I tried to comprehend the concept of
stratification and apply it with Survey to a small example. My question
is more of theoretic nature, so I apologize if this does not fully fit
this board's intention, but I have come to a complete stop in my efforts
and need an expert to help me along. Please help:
age<-matrix(c(rep(1,5), rep(2,3), 1:8, rep(3,5), rep(4,3), rep(5,5),
rep(3,3), rep(15,5), rep(12,3), 23,25,27,21,22, 33,27,29), ncol=6, byrow=F)
colnames(age)<-c("stratum", "id", "weight", "nh", "Nh", "y")
age<-as.data.frame(age)
## create survey design object
age.des1<-svydesign(ids=~id, strata=~stratum, weight=~Nh, data=age)
svymean(~y, age.des1)
## gives mean 25.568, SE 0.9257
age.des2<-svydesign(ids=~id, strata=~stratum, weight=~I(nh/Nh), data=age)
svymean(~y, age.des2)
## gives mean 25.483, SE 0.9227
age.des3<-svydesign(ids=~id, strata=~stratum, weight=~weight, data=age)
svymean(~y, age.des3)
## gives mean 26.296, SE 0.9862
age.des4<-svydesign(ids=~id, strata=~stratum, data=age)
svymean(~y, age.des4)
## gives mean 25.875, SE 0.9437
age.des3 is the only estimator I am able to compute per hand correctly.
It is stratified random sampling with inverse probablility weighting
with weight= nh/Nh ## sample size/ stratum size.
Basically, I thought the option weight=~Nh as well as weight=~I(nh/Nh)
would result in the same number, but it does not. I am reading
Thompson(02), Cochran(77) and of course Lumley on his Survey package,
but I can't find my mistake.
I thought the Hansen-Hurwitz estimator per stratum offers the right numbers:
p1=5/15, p2=3/12, so y1.total=1/5*(3*118), y2.total=1/3*(4*89) and the
stratified estimator with this design should be:
1/27(y1.total+y2.total), obviously wrong. How on earth do I get the
numbers Survey is calculating?
I am very sorry to bother you with this problem, however, I didn't find
anybody who was willing to help me.
Thank you so much
Mark
More information about the R-help
mailing list