Installing Packages

#install.packages("survey")
#install.packages("sampling")
#install.packages("writexl")

Loading Packages

library(readxl)
## Warning: package 'readxl' was built under R version 4.5.3
library(survey)
## Warning: package 'survey' was built under R version 4.5.3
## Loading required package: grid
## Loading required package: Matrix
## Loading required package: survival
## 
## Attaching package: 'survey'
## The following object is masked from 'package:graphics':
## 
##     dotchart
library(sampling)
## Warning: package 'sampling' was built under R version 4.5.3
## 
## Attaching package: 'sampling'
## The following objects are masked from 'package:survival':
## 
##     cluster, strata

Reading Data

base<-readRDS(file = "Marco.rds")

Creating the Bogotá Dataset

base_bogota<-subset(base,base$COLE_COD_MCPIO_UBICACION==11001)

Creating the Academic School Variable

base_bogota$academico=ifelse(base_bogota$COLE_CARACTER=="ACADÉMICO",1,0)

Sample Selection under Bernoulli Sampling

s.ber=function(base,pi,seed){
  set.seed(seed)
  base$u=runif(N)
  base$m=ifelse(base$u<pi,1,0)
  base=subset(base,base$m==1)
  base$pik=rep(pi,nrow(base))
  return(base)
}

e.ber.prop=function(variable,muestra,N){
  for(i in 1:length(table(muestra[,c(variable)]))){
    muestra$ind<-ifelse(muestra[,c(variable)]==unique(muestra[,c(variable)])[i],1,0)
    muestra.i=subset(muestra,muestra[,c(variable)]==unique(muestra[,c(variable)])[i])
    p.i=sum(muestra$ind*muestra$Fexp)/N
    se.i=1/N*sqrt((1/pi)*((1/pi)-1)*sum(muestra$ind^2))
    tab.i=data.frame(Categoria=paste0(variable,unique(muestra[,c(variable)])[i]),p_ber=round(p.i,4),se_ber=round(se.i,4))
    ifelse(i==1,(tabla=tab.i),(tabla=rbind(tabla,tab.i)))
  }
  return(tabla)
}

n=500
N=nrow(base_bogota)
seed=123
pi=n/N
muestra_be=s.ber(base_bogota,pi,seed)

Creating the “dsgn” Objec

muestra_be$Fexp=1/muestra_be$pik
dsgn_Be=svydesign(id=~1,data=muestra_be,weights=~Fexp,fpc=~rep(N,nrow(muestra_be)))

Total Estimation Under the Bernoulli Design

(est_Be=svytotal(~academico,dsgn_Be,deff=T))
##             total      SE   DEff
## academico 68419.2  1367.6 1.0001

Direct Computation of the Total Estimator

(total_est=sum(muestra_be$academico*muestra_be$Fexp))
## [1] 68419.23

Direct Computation of the Standard Error of the Total Estimator

(se=sqrt((1/pi)*((1/pi)-1)*sum(muestra_be$academico^2)))
## [1] 3356.514

The Alternative Estimator of the Total

(t_a=(N*sum(muestra_be$academico))/nrow(muestra_be))
## [1] 69110.34
(se_a=sqrt(N^2*((1-pi)/nrow(muestra_be))*(1-(1/nrow(muestra_be)))*var(muestra_be$academico)))
## [1] 1379.932
est_Be=as.data.frame(est_Be)


(Tabla_BE=data.frame(metodo=c("Survey","t_pi","t.alternativo"),t_est=c(est_Be[1,1],total_est,t_a),se=c(est_Be[1,2],se,se_a)))
##          metodo    t_est       se
## 1        Survey 68419.23 1367.556
## 2          t_pi 68419.23 3356.514
## 3 t.alternativo 69110.34 1379.932
(CVm1<-se/est_Be[1,1])
## [1] 0.04905804
(CVm2<-se_a/t_a)
## [1] 0.01996708

Estimation of the Global Mean Score

(est1_Be=svymean(~PUNT_GLOBAL,dsgn_Be,deff=T))
##                 mean       SE   DEff
## PUNT_GLOBAL 265.5616   2.0894 1.0001

Direct Computation of the Mean Estimator

(media_Be=sum(muestra_be$PUNT_GLOBAL*muestra_be$Fexp)/N)
## [1] 262.906

Direct Computation of the Standard Error of the Mean Estimator

(se_1=1/N*sqrt((1/pi)*((1/pi)-1)*sum(muestra_be$PUNT_GLOBAL^2)))
## [1] 11.96087

Poisson Sampling

Sample Selection under Poisson Sampling

s.pois=function(base,pi_k,seed){
  set.seed(seed)
  base$u=runif(N)
  base$pik=pi_k
  base$m=ifelse(base$u<base$pik,1,0)
  base=subset(base,base$m==1)
  return(base)
}

Sample Selection Under the Poisson Design

n=500
N=nrow(base_bogota)
seed=123
x=rnorm(nrow(base_bogota),10,1)
pi_k=n*x/sum(x)
muestra_po=s.pois(base_bogota,pi_k,seed)

Creating the “dsgn” Objec for the Poisson Design

muestra_po$Fexp=1/muestra_po$pik
dsgn_po=svydesign(id=~1,data=muestra_po,weights=~Fexp,fpc=~rep(N,nrow(muestra_po)))

Estimator of the Total Under the Poisson Design

(est_po=svytotal(~academico,dsgn_po,deff=T))
##             total      SE   DEff
## academico 70080.6  1390.7 1.0726

##Direct Computation of the Total Estimator Under the Poisson Design

total_est=sum(muestra_po$academico*muestra_po$Fexp)

Direct Computation of the Standard Error of the Total Estimator Under the Poisson Design

(se_po=sqrt(sum(1/muestra_po$pik*(1/muestra_po$pik-1)*muestra_po$academico^2)))
## [1] 3416.836

The Alternative Estimator of the Total Under the Poisson Design

(t_a_po=N*sum(muestra_po$academico/muestra_po$pik)/sum(muestra_po$Fexp))
## [1] 69918.83