###################hw 2######################################
##prob 1
library(survey)
data(api)
mean(apipop$enroll,na.rm=TRUE)
sum(apipop$enroll,na.rm=TRUE)

##stratified sample, stratified by school type (elementary, middle and high school)
dstrat<-svydesign(id=~1,strat=~stype,weights=~pw,data=apistrat,fpc=~fpc)
summary(dstrat)
svymean(~enroll,dstrat)
confint(svymean(~enroll,dstrat))
svytotal(~enroll,dstrat)
confint(svytotal(~enroll,dstrat))

##one-stage cluster sample, clusted by district
dclus1<-svydesign(id=~dnum,weights=~pw,data=apiclus1,fpc=~fpc)
summary(dclus1)
svymean(~enroll,dclus1)
confint(svymean(~enroll,dclus1))
svytotal(~enroll,dclus1)
confint(svytotal(~enroll,dclus1))

##two-stage cluster sample, schools within district clusters
dclus2<-svydesign(id=~dnum+snum,data=apiclus2,fpc=~fpc1+fpc2)
summary(dclus2)
svymean(~enroll,dclus2,na.rm=TRUE)
confint(svymean(~enroll,dclus2,na.rm=TRUE))
svytotal(~enroll,dclus2,na.rm=TRUE)
confint(svytotal(~enroll,dclus2,na.rm=TRUE))

##SRS
dsrs<-svydesign(id=~1,data=apisrs,fpc=~fpc)
summary(dsrs)
svymean(~enroll,dsrs)
confint(svymean(~enroll,dsrs))
svytotal(~enroll,dsrs)
confint(svytotal(~enroll,dsrs))

##prob 2
##hw##
library(survey)
library(SDaA)
data(coots)
csize<-coots$csize
volume<-coots$volume
oweight<-csize/2
ybar<-sum(oweight*volume)/sum(oweight)

##treat the cluster sample as stratified random sample
dstr <- svydesign(id = ~1, strata = ~clutch, weight=oweight, data = coots)
smean<-svymean(~volume, dstr)
smean

dcount<-function(){
     
     data<-coots
     volume<-coots$volume
     ln<-length(data$clutch)
     a<-table(data$clutch)
     rn<-row.names(a)
     ll<-length(rn)
     thetahat<-rep(0,ll)
     

    for(i in 1:ll){
        data<-coots
        csize<-data$csize
        oweight<-csize/2
       for(j in 1:ln){
            if(rn[i]==data$clutch[j]){
                        oweight[j]<-0
                      
                         }
                     }

       for(k in 1:ln){
            
           oweight[k]<-oweight[k]*ll/(ll-1)
         
           }

thetahat[i]<-sum(oweight*volume)/sum(oweight)
vjk<-(ll-1)/ll*sum((thetahat-ybar)^2)

}
list(thetahat,vjk)
}