library(survey)
source("c:/courses/st446/rcode/confintt.r")

# Cluster sample with unequal-size clusters (Figure 10)

M0 = 400
N = 49
n =  8
y <- c(401,337,279,321,280,171,187,216)
clusterid <- c(1,2,3,4,5,6,7,8)

# Method 1: SRS of clusters using cluster ratios

fpc1 <- c(rep(N,n))
Mvec <- c(16,16,8,8,8,4,4,4)

ratio_ttl  <- M0*y      
ratio_mn   <- y

Fig10 <- data.frame(cbind(fpc1,ratio_ttl,ratio_mn,Mvec))
Fig10

# Create the sampling design
dsgn10 <- svydesign(data=Fig10,id=~1,fpc=~fpc1)

# Method 1: Estimation of the population mean

estmean1 <- svyratio(~ratio_mn,~Mvec,design=dsgn10)
confint.t(estmean1,tdf=n-1,level=.95)

# Method 1: Estimation of the population total

esttotal1 <- svyratio(~ratio_ttl,~Mvec,design=dsgn10)
confint.t(esttotal1,tdf=n-1,level=.95)


# Method 2: SRS of clusters using cluster totals

fpc2 <- c(rep(N,8))
wgt2a <- c(rep(N/n,n))
wgt2b <- wgt2a/M0
#wgt2b <- c(rep(N/(n*M0),n))

Fig10a <- data.frame(cbind(clusterid,y,wgt2a,wgt2b,fpc2))
Fig10a
dsgn10a <-svydesign(ids=~clusterid,weights=~wgt2a,fpc=~fpc2,data=Fig10a)
dsgn10b <-svydesign(ids=~clusterid,weights=~wgt2b,fpc=~fpc2,data=Fig10a)

# Method 2: Estimation of population total

esttotal2 <- svytotal(~y,design=dsgn10a)
print(esttotal2,digits=15)
confint.t(esttotal2,level=.95,tdf=n-1)

# Method 2: Estimation of population mean

estmean2 <- svytotal(~y,design=dsgn10b)
print(estmean2,digits=15)
confint.t(estmean2,level=.95,tdf=n-1)


# Method 3: Sampling proportional to size with replacement

p <- Mvec/M0
t  <- y/p
t2 <- t/M0

Fig10c <- data.frame(t,t2)
Fig10c
dsgn10c <-svydesign(ids=~1,data=Fig10c)

# Method 3: Estimation of population total

esttotal3 <- svymean(~t,design=dsgn10c)
print(esttotal3,digits=15)
confint.t(esttotal3,level=.95,tdf=n-1)

# Method 3: Estimation of population mean

estmean3 <- svymean(~t2,design=dsgn10c)
print(estmean3,digits=15)
confint.t(estmean3,level=.95,tdf=n-1)