Installing Packages
#install.packages("readxl")
#install.packages("survey")
#install.packages("sampling")
#install.package("srvyr")
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
library(srvyr)
## Warning: package 'srvyr' was built under R version 4.5.3
##
## Attaching package: 'srvyr'
## The following object is masked from 'package:stats':
##
## filter
Reading Data
base<-readRDS(file = "Marco.rds")
Creating the Sampling Frame of Schools
colegios=base%>%group_by(COLE_AREA_UBICACION,COLE_COD_MCPIO_UBICACION,COLE_COD_DANE_ESTABLECIMIENTO)%>%
summarise(x=n(),y=mean(PUNT_GLOBAL))
## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by COLE_AREA_UBICACION,
## COLE_COD_MCPIO_UBICACION, and COLE_COD_DANE_ESTABLECIMIENTO.
## ℹ Output is grouped by COLE_AREA_UBICACION and COLE_COD_MCPIO_UBICACION.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(COLE_AREA_UBICACION, COLE_COD_MCPIO_UBICACION,
## COLE_COD_DANE_ESTABLECIMIENTO))` for per-operation grouping
## (`?dplyr::dplyr_by`) instead.
Calculating the Absolute Frequencies of the School’s Location Area
Variable
table(colegios$COLE_AREA_UBICACION)
##
## RURAL URBANO
## 3375 6935
Ordering by the Stratification Variable: School Location Area
colegios=colegios[order(colegios$COLE_AREA_UBICACION),]
Neyman Allocation
n=500
(sy=colegios%>%group_by(COLE_AREA_UBICACION)%>%summarise(sy=sd(x)))
## # A tibble: 2 × 2
## COLE_AREA_UBICACION sy
## <chr> <dbl>
## 1 RURAL 23.2
## 2 URBANO 63.2
(ne1=round(n*nrow(colegios)*sy$sy/sum(nrow(colegios)*sy$sy),0))
## [1] 134 366
sy$neyman=ne1
Proportional Allocation
(Ne=colegios%>%group_by(COLE_AREA_UBICACION)%>%summarise(Ne=n()))
## # A tibble: 2 × 2
## COLE_AREA_UBICACION Ne
## <chr> <int>
## 1 RURAL 3375
## 2 URBANO 6935
(ne2=round(n*Ne$Ne/nrow(colegios),0))
## [1] 164 336
sy$Ne=Ne$Ne
sy$prop=ne2
Allocation Proportional to the x-Total
(te=colegios%>%group_by(COLE_AREA_UBICACION)%>%summarise(te=sum(x)))
## # A tibble: 2 × 2
## COLE_AREA_UBICACION te
## <chr> <int>
## 1 RURAL 85745
## 2 URBANO 460467
(ne3=round(n*te$te/sum(te$te),0))
## [1] 78 422
sy$te=te$te
sy$stotal=ne3
Creating the “Ne” Column
colegios$Ne=ifelse(colegios$COLE_AREA_UBICACION=="RURAL",3375,6935)
Selecting the Stratified Sample
set.seed(123)
s1=strata(colegios, c("COLE_AREA_UBICACION"), ne1, method="srswor")
m1=colegios[s1$ID_unit,]
m1$pik1=s1$Prob
Creating the “dsgn” Objec
dsgn_str=svydesign(id=~1,data=m1, strata=~COLE_AREA_UBICACION ,fpc=~Ne, weights=~1/pik1)
Calculating the Mean Estimator under STSI (Simple Random Sampling
Design in all Strata)
(est_str=svymean(~y,dsgn_str))
## mean SE
## y 246.9 1.5905
salida=function(est,alpha){
est=as.data.frame(est)
names(est)[2]="se"
est$cv=100*(est$se/est[,1])
est$ic_low=est[,1]-qnorm(1-alpha/2)*est$se
est$ic_upp=est[,1]+qnorm(1-alpha/2)*est$se
return(round(est,2))
}
alpha=0.05
(tabla_str=salida(est_str,alpha))
## mean se cv ic_low ic_upp
## y 246.9 1.59 0.64 243.78 250.01
Calculating the Mean Estimator under STSI by Stratum
(est_str1=svyby(~y,~COLE_AREA_UBICACION,dsgn_str,svymean))
## COLE_AREA_UBICACION y se
## RURAL RURAL 220.7838 2.533839
## URBANO URBANO 254.0899 1.904719
est_str1=as.data.frame(est_str1)
est_str1=est_str1[,-1]
alpha=0.05
(tabla_str1=salida(est_str1,alpha))
## y se cv ic_low ic_upp
## RURAL 220.78 2.53 1.15 215.82 225.75
## URBANO 254.09 1.90 0.75 250.36 257.82
Creating a New Data Frame with a New Stratification Variable
colegios1=base%>%group_by(COLE_CARACTER,COLE_COD_MCPIO_UBICACION,COLE_COD_DANE_ESTABLECIMIENTO)%>%
summarise(x=n(),y=mean(PUNT_GLOBAL))
## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by COLE_CARACTER, COLE_COD_MCPIO_UBICACION,
## and COLE_COD_DANE_ESTABLECIMIENTO.
## ℹ Output is grouped by COLE_CARACTER and COLE_COD_MCPIO_UBICACION.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(COLE_CARACTER, COLE_COD_MCPIO_UBICACION,
## COLE_COD_DANE_ESTABLECIMIENTO))` for per-operation grouping
## (`?dplyr::dplyr_by`) instead.
Calculating the Absolute Frequencies of the “school character”
Variable
(t=as.data.frame(table(colegios1$COLE_CARACTER)))
## Var1 Freq
## 1 409
## 2 - 27
## 3 ACADÉMICO 5838
## 4 NO APLICA 353
## 5 TÉCNICO 1286
## 6 TÉCNICO/ACADÉMICO 2276
names(t)=c("COLE_CARACTER","Ne")
colegios1=merge(colegios1,t)
Selecting some Categories of the “school character” Variable
colegios1=subset(colegios1,colegios1$COLE_CARACTER=="ACADÉMICO"|
colegios1$COLE_CARACTER=="TÉCNICO"|
colegios1$COLE_CARACTER=="TÉCNICO/ACADÉMICO")
Ordering by the Stratification Variable: “school character”
colegios1=colegios1[order(colegios1$COLE_CARACTER),]
colegios1=colegios1[!is.na(colegios1$COLE_CARACTER),]
Neyman Allocation
n=400
(sy1=colegios1%>%group_by(COLE_CARACTER)%>%summarise(sy1=sd(x)))
## # A tibble: 3 × 2
## COLE_CARACTER sy1
## <chr> <dbl>
## 1 ACADÉMICO 50.7
## 2 TÉCNICO 45.6
## 3 TÉCNICO/ACADÉMICO 71.9
(ne11=round(n*nrow(colegios1)*sy1$sy1/sum(nrow(colegios1)*sy1$sy1),0))
## [1] 121 109 171
Selecting the Sample
s11=strata(colegios1, c("COLE_CARACTER"), ne11, method="srswor")
m11=colegios1[s11$ID_unit,]
m11$pik1=s11$Prob
Creating the “dsgn” Objec
dsgn_str1=svydesign(id=~1,data=m11, strata=~COLE_CARACTER ,fpc=~m11$Ne, weights=~1/pik1)
Calculating the Mean Estimator under STSI (Simple Random Sampling
Design in all Strata)
(est_str2=svymean(~y,dsgn_str1))
## mean SE
## y 238.98 1.6314
alpha=0.05
(tabla_str2=salida(est_str2,alpha))
## mean se cv ic_low ic_upp
## y 238.98 1.63 0.68 235.78 242.18
Calculating the Mean estimator under STSI by Stratum
(est_str3=svyby(~y,~COLE_CARACTER,dsgn_str1,svymean))
## COLE_CARACTER y se
## ACADÉMICO ACADÉMICO 246.9753 3.355436
## TÉCNICO TÉCNICO 236.6103 2.928007
## TÉCNICO/ACADÉMICO TÉCNICO/ACADÉMICO 239.1219 2.211636
est_str3=as.data.frame(est_str3)
est_str3=est_str3[,-1]
alpha=0.05
(tabla_str3=salida(est_str3,alpha))
## y se cv ic_low ic_upp
## ACADÉMICO 246.98 3.36 1.36 240.40 253.55
## TÉCNICO 236.61 2.93 1.24 230.87 242.35
## TÉCNICO/ACADÉMICO 239.12 2.21 0.92 234.79 243.46
Creating a New Data Frame with a New Stratification Variable:
“school schedule”
colegios2=base%>%group_by(COLE_JORNADA,COLE_COD_MCPIO_UBICACION,COLE_COD_DANE_ESTABLECIMIENTO)%>%
summarise(x=n(),y=mean(PUNT_GLOBAL))
## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by COLE_JORNADA, COLE_COD_MCPIO_UBICACION,
## and COLE_COD_DANE_ESTABLECIMIENTO.
## ℹ Output is grouped by COLE_JORNADA and COLE_COD_MCPIO_UBICACION.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(COLE_JORNADA, COLE_COD_MCPIO_UBICACION,
## COLE_COD_DANE_ESTABLECIMIENTO))` for per-operation grouping
## (`?dplyr::dplyr_by`) instead.
Calculating the Absolute Frequencies of the “school schedule”
Variable
(t1=as.data.frame(table(colegios2$COLE_JORNADA)))
## Var1 Freq
## 1 COMPLETA 2584
## 2 MAÑANA 5570
## 3 NOCHE 1081
## 4 SABATINA 1185
## 5 TARDE 1181
## 6 UNICA 1107
names(t1)=c("COLE_JORNADA","Ne")
colegios2=merge(colegios2,t1)
Ordering by the Stratification Variable: “school schedule”
colegios2=colegios2[order(colegios2$COLE_JORNADA),]
Neyman Allocation
n=400
(sy2=colegios2%>%group_by(COLE_JORNADA)%>%summarise(sy2=sd(x)))
## # A tibble: 6 × 2
## COLE_JORNADA sy2
## <chr> <dbl>
## 1 COMPLETA 31.3
## 2 MAÑANA 43.1
## 3 NOCHE 32.7
## 4 SABATINA 43.1
## 5 TARDE 39.0
## 6 UNICA 56.8
(ne12=round(n*nrow(colegios2)*sy2$sy2/sum(nrow(colegios2)*sy2$sy2),0))
## [1] 51 70 53 70 63 92
Selecting the Stratified Sample
s2=strata(colegios2, c("COLE_JORNADA"), ne12, method="srswor")
m2=colegios2[s2$ID_unit,]
m2$pik1=s2$Prob
Creating the “dsgn” Objec
dsgn_str2=svydesign(id=~1,data=m2, strata=~COLE_JORNADA ,fpc=~m2$Ne, weights=~1/pik1)
Calculating the Mean Estimator under STSI (Simple Random Sampling
Design in all Strata)
(est_str4=svymean(~y,dsgn_str2))
## mean SE
## y 233.84 1.3548
alpha=0.05
(tabla_str4=salida(est_str4,alpha))
## mean se cv ic_low ic_upp
## y 233.84 1.35 0.58 231.18 236.49
Calculating the Mean Estimator under STSI by Stratum
(est_str5=svyby(~y,~COLE_JORNADA,dsgn_str2,svymean))
## COLE_JORNADA y se
## COMPLETA COMPLETA 269.7504 5.315654
## MAÑANA MAÑANA 243.8771 3.808833
## NOCHE NOCHE 210.5204 2.891820
## SABATINA SABATINA 203.1628 2.169995
## TARDE TARDE 238.3909 3.339034
## UNICA UNICA 250.4649 2.629792
est_str5=as.data.frame(est_str5)
est_str5=est_str5[,-1]
alpha=0.05
(tabla_str5=salida(est_str5,alpha))
## y se cv ic_low ic_upp
## COMPLETA 269.75 5.32 1.97 259.33 280.17
## MAÑANA 243.88 3.81 1.56 236.41 251.34
## NOCHE 210.52 2.89 1.37 204.85 216.19
## SABATINA 203.16 2.17 1.07 198.91 207.42
## TARDE 238.39 3.34 1.40 231.85 244.94
## UNICA 250.46 2.63 1.05 245.31 255.62
Creating a New data Frame with a Two Stratification variables:
“school schedule” and “school’s location area”
colegios3=base%>%group_by(COLE_JORNADA,COLE_AREA_UBICACION,
COLE_COD_MCPIO_UBICACION,COLE_COD_DANE_ESTABLECIMIENTO)%>%
summarise(x=n(),y=mean(PUNT_GLOBAL))
## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by COLE_JORNADA, COLE_AREA_UBICACION,
## COLE_COD_MCPIO_UBICACION, and COLE_COD_DANE_ESTABLECIMIENTO.
## ℹ Output is grouped by COLE_JORNADA, COLE_AREA_UBICACION, and
## COLE_COD_MCPIO_UBICACION.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(COLE_JORNADA, COLE_AREA_UBICACION,
## COLE_COD_MCPIO_UBICACION, COLE_COD_DANE_ESTABLECIMIENTO))` for per-operation
## grouping (`?dplyr::dplyr_by`) instead.
Creating the “estrato” variable
colegios3$estrato=paste(colegios3$COLE_JORNADA,colegios3$COLE_AREA_UBICACION,sep="-")
(t2=as.data.frame(table(colegios3$estrato)))
## Var1 Freq
## 1 COMPLETA-RURAL 864
## 2 COMPLETA-URBANO 1759
## 3 MAÑANA-RURAL 2087
## 4 MAÑANA-URBANO 3525
## 5 NOCHE-RURAL 129
## 6 NOCHE-URBANO 953
## 7 SABATINA-RURAL 234
## 8 SABATINA-URBANO 953
## 9 TARDE-RURAL 145
## 10 TARDE-URBANO 1038
## 11 UNICA-RURAL 337
## 12 UNICA-URBANO 776
names(t2)=c("estrato","Ne")
colegios3=merge(colegios3,t2)
Selecting the Stratified Sample
n=c(rep(40,12))
s3=strata(colegios3, c("COLE_JORNADA","COLE_AREA_UBICACION"), n, method="srswor")
m3=colegios3[s3$ID_unit,]
m3$pik1=s3$Prob
Creating the “dsgn” Objec
dsgn_str3=svydesign(id=~1,data=m3, strata=~estrato ,fpc=~m3$Ne, weights=~1/pik1)
Calculating the Mean Estimator under STSI (Simple Random Sampling
Design in all strata)
(est_str6=svymean(~y,dsgn_str3))
## mean SE
## y 214.58 1.227
alpha=0.05
(tabla_str6=salida(est_str6,alpha))
## mean se cv ic_low ic_upp
## y 214.58 1.23 0.57 212.18 216.99
Calculating the Mean Estimator under STSI by Stratum
(est_str7=svyby(~y,~estrato,dsgn_str3,svymean))
## estrato y se
## COMPLETA-RURAL COMPLETA-RURAL 239.4032 6.970808
## COMPLETA-URBANO COMPLETA-URBANO 282.1714 4.089806
## MAÑANA-RURAL MAÑANA-RURAL 222.6982 4.749680
## MAÑANA-URBANO MAÑANA-URBANO 244.7605 4.687515
## NOCHE-RURAL NOCHE-RURAL 198.2983 2.219854
## NOCHE-URBANO NOCHE-URBANO 203.4842 2.967289
## SABATINA-RURAL SABATINA-RURAL 195.1094 2.086668
## SABATINA-URBANO SABATINA-URBANO 210.9208 2.458919
## TARDE-RURAL TARDE-RURAL 215.5231 3.365092
## TARDE-URBANO TARDE-URBANO 240.6460 3.321648
## UNICA-RURAL UNICA-RURAL 235.3736 4.827345
## UNICA-URBANO UNICA-URBANO 254.9312 3.394311
est_str7=as.data.frame(est_str7)
est_str7=est_str7[,-1]
alpha=0.05
(tabla_str7=salida(est_str7,alpha))
## y se cv ic_low ic_upp
## COMPLETA-RURAL 239.40 6.97 2.91 225.74 253.07
## COMPLETA-URBANO 282.17 4.09 1.45 274.16 290.19
## MAÑANA-RURAL 222.70 4.75 2.13 213.39 232.01
## MAÑANA-URBANO 244.76 4.69 1.92 235.57 253.95
## NOCHE-RURAL 198.30 2.22 1.12 193.95 202.65
## NOCHE-URBANO 203.48 2.97 1.46 197.67 209.30
## SABATINA-RURAL 195.11 2.09 1.07 191.02 199.20
## SABATINA-URBANO 210.92 2.46 1.17 206.10 215.74
## TARDE-RURAL 215.52 3.37 1.56 208.93 222.12
## TARDE-URBANO 240.65 3.32 1.38 234.14 247.16
## UNICA-RURAL 235.37 4.83 2.05 225.91 244.83
## UNICA-URBANO 254.93 3.39 1.33 248.28 261.58