R Code for Systematic Sampling

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")

Reading the Bogotá Dataset

base_bogota=readRDS(file = "Marco Bogota.RDS")

Creating the Academic School Variable

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

Sample Selection under Systematic Sampling

s.sis=function(base,n,seed){
  N=nrow(base)
  a=round(N/n,0)
  set.seed(seed)
  u=round(runif(1,1,a),0)
  s=seq(u,u+a*(n-1),a) 
  base=base[s,]
  base$pik=1/a
  return(base)
}

n=700
N=nrow(base_bogota)
seed=123
muestra.sis=s.sis(base_bogota,n,seed)

Creating the “dsgn” Objec

muestra.sis$Fexp=1/muestra.sis$pik
dsgn.sis=svydesign(id=~1,data=muestra.sis,weights=~Fexp,fpc=~rep(N,n))

Total Estimation Under Systematic Sampling

(est.sis=svytotal(~academico,dsgn.sis,deff=T))

##             total      SE DEff
## academico 69502.0  1136.4    1

Direct Computation of the Total Estimator

f=n/N
total=sum(muestra.sis$academico*muestra.sis$Fexp)
se=sqrt(N^2*(1-f)/n*var(muestra.sis$academico))
est.sis=as.data.frame(est.sis)
data.frame(metodo=c("Survey","Formula"),t=c(est.sis[1,1],total),se=c(est.sis[1,2],se))

##    metodo     t       se
## 1  Survey 69502 1136.370
## 2 Formula 69502 1139.562

Estimation of the Mean

(est1.sis=svymean(~PUNT_GLOBAL,dsgn.sis,deff=T))

##                 mean       SE DEff
## PUNT_GLOBAL 266.0071   1.8255    1

Direct Computation of the Mean Estimator

media=sum(muestra.sis$PUNT_GLOBAL*muestra.sis$Fexp)/N
se=1/N*sqrt(N^2*(1-f)/n*var(muestra.sis$PUNT_GLOBAL)*n/(n-1))
est1.sis=as.data.frame(est1.sis)
data.frame(metodo=c("Survey","Formula"),media=c(est1.sis[1,1],media),se=c(est1.sis[1,2],se))

##    metodo    media       se
## 1  Survey 266.0071 1.825476
## 2 Formula 265.2621 1.826781

Estimating the Proportion

(est2.sis=svymean(~DESEMP_INGLES,dsgn.sis))

##                     mean     SE
## DESEMP_INGLESA- 0.307143 0.0174
## DESEMP_INGLESA1 0.304286 0.0173
## DESEMP_INGLESA2 0.234286 0.0160
## DESEMP_INGLESB+ 0.035714 0.0070
## DESEMP_INGLESB1 0.118571 0.0122

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
salida(est2.sis,alpha)

##                 mean   se    cv ic_low ic_upp
## DESEMP_INGLESA- 0.31 0.02  5.66   0.27   0.34
## DESEMP_INGLESA1 0.30 0.02  5.69   0.27   0.34
## DESEMP_INGLESA2 0.23 0.02  6.81   0.20   0.27
## DESEMP_INGLESB+ 0.04 0.01 19.57   0.02   0.05
## DESEMP_INGLESB1 0.12 0.01 10.27   0.09   0.14