On Thu, 26 May 2005, Mark Hempelmann wrote:
> 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)
Ok. Assuming that Nh are the population sizes in each stratum, you have
5/15 sampled in stratum 1 and 3/12 in stratum 2.
This can be specified in a number of ways
You can use
sampling weights of 15/5 and 12/3
sampling probabilities of 5/15 and 3/12
without or without specifiying the finite population correction. The
finite population correction can be specified as 15 and 12 or 5/15 and
3/12, and if the finite population correction is specified the weights are
then optional.
So
d1<-svydesign(ids=~id, strata=~stratum, weight=~I(Nh/nh), data=age)
d2<-svydesign(ids=~id, strata=~stratum, prob=~I(nh/Nh), data=age)
give the with-replacement design (agreeing with your age.des3) and
d3<-svydesign(ids=~id, strata=~stratum, weight=~I(Nh/nh),
fpc=~Nh,data=age)
d4<-svydesign(ids=~id, strata=~stratum, prob=~I(nh/Nh), fpc=~Nh,data=age)
d5<-svydesign(ids=~id, strata=~stratum, weight=~I(Nh/nh),
fpc=~I(nh/Nh),data=age)
d6<-svydesign(ids=~id, strata=~stratum, prob=~I(nh/Nh),
fpc=~I(nh/Nh),data=age)
d7<-svydesign(ids=~id, strata=~stratum, fpc=~Nh,data=age)
d8<-svydesign(ids=~id, strata=~stratum, fpc=~I(nh/Nh),data=age)
all give the without-replacement design. We get> svymean(~y,d1)
mean SE
y 26.296 0.9862> svymean(~y,d2)
mean SE
y 26.296 0.9862> svymean(~y,d3)
mean SE
y 26.296 0.8364> svymean(~y,d4)
mean SE
y 26.296 0.8364> svymean(~y,d5)
mean SE
y 26.296 0.8364> svymean(~y,d6)
mean SE
y 26.296 0.8364> svymean(~y,d7)
mean SE
y 26.296 0.8364> svymean(~y,d8)
mean SE
y 26.296 0.8364
Now, looking at your examples> ## 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
This is wrong: the sampling weight is Nh/nh, not Nh
> age.des2<-svydesign(ids=~id, strata=~stratum, weight=~I(nh/Nh),
data=age)
> svymean(~y, age.des2)
> ## gives mean 25.483, SE 0.9227
This is wrong: the sampling weight is Nh/nh. You need prob=~I(nh/Nh) to
specify sampling fractions.
> age.des3<-svydesign(ids=~id, strata=~stratum, weight=~weight, data=age)
> svymean(~y, age.des3)
> ## gives mean 26.296, SE 0.9862
This is correct and agrees with d1 and d2
> age.des4<-svydesign(ids=~id, strata=~stratum, data=age)
> svymean(~y, age.des4)
> ## gives mean 25.875, SE 0.9437
This is a stratified, unweighted mean, ie mean(age$y).
> 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.
No, it does not. A weight of 3 is not the same as a weight of 1/3. With
the finite population correction it is safe to assume that numbers less
than 1 are sampling fractions and numbers greater than 1 are population
sizes, but this isn't safe when it comes to weights. It is possible that
someone could want to use sampling weights less than 1.
>
> 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.
Since this gives a mean of 7.01 for numbers around 25 it can't be right.
You have divided by sample size twice. You should have
y1.total<-3*118
y2.total<-4*89
You then will get (y1.total+y2.total)/27 to be 26.29630, in agreement
with svymean().
-thomas