##############chap 2: simple random sampling#################

#########Example 2.5,2.6 and 2.10####################################
###data agpop.dat and agsrs.dat The U.S government conducts a census of agriculture every five years,
##collecting data on all farms in the 50 states. The census of agriculture provides data on number of farms,
##the total acreage devoted to farms, farm size, yield of different crops, and a wide variety of other agricuture measures for 
#N=3078 counties and county-equivalent in the United States. agpop.dat is the population data and agsrs.dat is an simple random sample of 
#size 300 from agpop.dat

install.packages("survey")
install.packages("SDaA")
library(survey) #survey package
library(SDaA) #data sets package

#############Analysis of a simple random sample agsrs.dat
data(agsrs)
agsrs[1:10,]
agsrs10<-agsrs[1:10,]
fix(agsrs10)
n<-nrow(agsrs)
dsrs <- svydesign(id = ~1, fpc=rep(3078,300), data = agsrs)
smean<-svymean(~acres92,dsrs)
smean
confint(smean, level=.95,df=n-1) #without df, default is using normal distribution

#calculate mean and standard error by direct formula
ybar<-mean(agsrs$acres92)
hatvybar<-(1-300/3078)*var(agsrs$acres92)/300
seybar<-sqrt(hatvybar)
seybar
#calculate confidence interval by direct formula using t distribution
smean2<-as.data.frame(smean)
Mean_CI <- c(smean2[1] - qt(.975, n-1)*smean2[2], smean2[1] + qt(.975, n-1)*smean2[2])
names(Mean_CI) <- c("lower", "upper")
Mean_CI


#calculate total
stotal<-svytotal(~acres92,dsrs)
stotal
confint(stotal, level=.95,df=n-1)
s<-3078*smean2[1]
v2<-3078*smean2[2]


######calculate proportion
agsrsp<-agsrs
for(i in 1:n){if (agsrsp[i,3]<2e+05) agsrsp[i,3]=1
              if (agsrsp[i,3]>=2e+05) agsrsp[i,3]=0} 
dsrsp <- svydesign(id = ~1, fpc=rep(3078,300), data = agsrsp)
smeanp<-svymean(~acres92,dsrsp)
smeanp
confint(smeanp, level=.95,df=n-1)



###########calculate the population parameters using agpop.dat
data(agpop)
agpop[1:10,]
n2<-nrow(agpop)
n2
agpopmean<-mean(agpop$acres92)
agpopmean
agpoptotal<-sum(agpop$acres92)
agpoptotal
stdybar<-sqrt((1-n/3078)*var(agpop$acres92)/n)
stdybar