Muestreo estratificado mas y muestreo pipt conglomerados
Para empezar cargemos las librerías requeiridas
library(RCurl) #Para leer archivos de la web colgados en github
library(sampling) #Para selecciones de muestra
library(TeachingSampling) #Para selecciones de muestra y estimación
library(survey) #Para estimaciónd e muestras complejas
library(sqldf) #Realizar consultas en sql
library(ggplot2)
Para empezar leemos la base de datos, es importante indicar cuáles variables son caracter, cuáles son numéricas, la variable evaluados tiene como separador decimal una coma, con gsub se cambia a punto.
options(scipen = 11111)
url <- getURL("https://raw.githubusercontent.com/finiterank/saber_notebooks/master/r-scripts/data/saber2013.csv")
saber2013 <- read.csv(text = url)
clases = rep(NA, 26)
clases[1:5] = rep("character", 5)
clases[6:8] = "factor"
clases[26] = "factor"
saber2013 <- read.csv(text = url, colClasses = clases)
saber2013$evaluados <- as.numeric(gsub(",", "", as.character(saber2013$evaluados)))
En el identificador del colegio se presentan 525 duplicados, por el momento no tengo explicación satisfactoria de por que estos duplicados, para tratamiento de este ejercicio se eliminarán:
table(duplicated(saber2013$codinst))
##
## FALSE TRUE
## 12186 525
saber2013 = saber2013[!duplicated(saber2013$codinst), ]
table(duplicated(saber2013$codinst))
##
## FALSE
## 12186
Se obtiene finalmente una base de datos de 12186 registros. Consideremos un diseño muestral estratificado MAS, se considerará como estrato el cruce de las variables naturaleza y jornada
saber2013$estrato = paste0("Naturaleza: ", saber2013$naturaleza, " y Jornada:",
saber2013$jornada)
niveles = c("Naturaleza: NO OFICIAL y Jornada:COMPLETA U ORDINARIA", "Naturaleza: NO OFICIAL y Jornada:MAÑANA",
"Naturaleza: NO OFICIAL y Jornada:TARDE", "Naturaleza: NO OFICIAL y Jornada:NOCHE",
"Naturaleza: NO OFICIAL y Jornada:SABATINA - DOMINICAL", "Naturaleza: OFICIAL y Jornada:COMPLETA U ORDINARIA",
"Naturaleza: OFICIAL y Jornada:MAÑANA", "Naturaleza: OFICIAL y Jornada:TARDE",
"Naturaleza: OFICIAL y Jornada:NOCHE", "Naturaleza: OFICIAL y Jornada:SABATINA - DOMINICAL")
saber2013$estrato = factor(saber2013$estrato, levels = niveles)
table(saber2013$estrato)
##
## Naturaleza: NO OFICIAL y Jornada:COMPLETA U ORDINARIA
## 1735
## Naturaleza: NO OFICIAL y Jornada:MAÑANA
## 1243
## Naturaleza: NO OFICIAL y Jornada:TARDE
## 272
## Naturaleza: NO OFICIAL y Jornada:NOCHE
## 324
## Naturaleza: NO OFICIAL y Jornada:SABATINA - DOMINICAL
## 362
## Naturaleza: OFICIAL y Jornada:COMPLETA U ORDINARIA
## 2105
## Naturaleza: OFICIAL y Jornada:MAÑANA
## 3536
## Naturaleza: OFICIAL y Jornada:TARDE
## 1294
## Naturaleza: OFICIAL y Jornada:NOCHE
## 873
## Naturaleza: OFICIAL y Jornada:SABATINA - DOMINICAL
## 442
Se desea estimar el promedio del puntaje de matemáticas y el total de evaluados, utilizando un diseño estratificado, el estrato construido anteriormente se utilizará. En la práctica las variables de interés no se disponen, sin emb aro en este ejercicio se ilustra como se selecciona la muestra y como se lleva a cabo el proceso de estimación. En primer lugar procederemos a seleccionar la muestra
cons1 = sqldf("select estrato, count(*) as 'Nh' from saber2013 group by estrato")
cons1$nh = ceiling(0.2 * cons1$Nh)
set.seed(31052014)
indica_estmas = strata(saber2013, stratanames = "estrato", size = cons1$nh,
method = "srswor", description = T)
## Stratum 1
##
## Population total and number of selected units: 3536 347
## Stratum 2
##
## Population total and number of selected units: 1243 249
## Stratum 3
##
## Population total and number of selected units: 1735 65
## Stratum 4
##
## Population total and number of selected units: 272 73
## Stratum 5
##
## Population total and number of selected units: 324 55
## Stratum 6
##
## Population total and number of selected units: 873 421
## Stratum 7
##
## Population total and number of selected units: 1294 708
## Stratum 8
##
## Population total and number of selected units: 2105 175
## Stratum 9
##
## Population total and number of selected units: 442 89
## Stratum 10
##
## Population total and number of selected units: 362 259
## Number of strata 10
## Total number of selected units 2441
mue_estmas = saber2013[indica_estmas$ID_unit, ]
mue_estmas = merge(mue_estmas, cons1)
A continuación se realiza la estimación del promedio de matemáticas, se estima primero el total, y finalmente se divida por N = 12186 (el número de colegios).
cons2 = sqldf("select estrato, max(Nh) as 'Nh', count(*) as 'nh', sum(matematica) as 'sum_mat', variance(matematica) as 'S2_mat',\nsum(evaluados) as 'sum_eval', variance(evaluados) as 'S2_eval'\nfrom mue_estmas group by estrato")
fexp = cons2$Nh/cons2$nh
facvar = ((cons2$Nh^2)/cons2$nh) * (1 - (cons2$nh/cons2$Nh))
cons2$th_mat = cons2$sum_mat * fexp
cons2$Vh_mat = cons2$S2_mat * facvar
cons2$th_eval = cons2$sum_eval * fexp
cons2$Vh_eval = cons2$S2_eval * facvar
cons3 = sqldf("select sum(th_mat) as 't_mat',sum(th_eval) as 't_eval',\nsum(Vh_mat ) as 'Vest_mat ', sum(Vh_eval) as 'Vest_eval' from cons2")
cons3$Prom_mat = cons3$t_mat/nrow(saber2013)
cons3$Vest_prommat = cons3$Vest_mat/(nrow(saber2013)^2)
cons3$cveProm_mat = 100 * sqrt(cons3$Vest_prommat)/cons3$Prom_mat
cons3$cvetot_eval = 100 * sqrt(cons3$Vest_eval)/cons3$t_eval
cons3 = round(cons3, 2)
salida1a = cons3[c(5, 6, 7)]
salida1a$LI95 = salida1a$Prom_mat - 1.96 * sqrt(salida1a$Vest_prommat)
salida1a$LS95 = salida1a$Prom_mat + 1.96 * sqrt(salida1a$Vest_prommat)
salida2a = cons3[c(2, 4, 8)]
salida2a$LI95 = salida2a$t_eval - 1.96 * sqrt(salida2a$Vest_eval)
salida2a$LS95 = salida2a$t_eval + 1.96 * sqrt(salida2a$Vest_eval)
Se obtien las siguientes estimaciones:
salida1a
## Prom_mat Vest_prommat cveProm_mat LI95 LS95
## 1 44.04 0.03 0.4 43.7 44.38
salida2a
## t_eval Vest_eval cvetot_eval LI95 LS95
## 1 584018 222821439 2.56 554761 613276
Con el paquete survey se pueden obtener resultados más inmediatos:
disenoestmas = svydesign(ids = ~1, strata = ~estrato, fpc = ~Nh, data = mue_estmas)
salida1A = svytotal(~matematica, disenoestmas)
salida1A = as.data.frame(salida1A)
salida1A$Prom = (salida1A[, 1])/nrow(saber2013)
salida1A$Var = (salida1A[, 2]^2)/(nrow(saber2013)^2)
salida1A[, 1] = NULL
salida1A[, 1] = NULL
salida1A$cve = round(100 * sqrt(salida1A$Var)/salida1A$Prom, 2)
salida2A = svytotal(~evaluados, disenoestmas)
salida2A = as.data.frame(salida2A)
salida2A$Var = salida2A[, 2]^2
salida2A[, 2] = NULL
salida2A$cve = round(100 * sqrt(salida2A$Var)/salida2A$total, 2)
salida1A
## Prom Var cve
## matematica 44.04 0.03027 0.4
salida2A
## total Var cve
## evaluados 584018 222821439 2.56
Con el paquete survey la estimación del promedio es realizada como una razón, que para un diseño ESTMAS coincide que si asumieran los N conocidos.
salida1C = svymean(~matematica, disenoestmas)
salida1C = as.data.frame(salida1C)
salida1C$cve = round(100 * salida1C[, 2]/salida1C[, 1], 2)
En el segundo punto se realizó un diseño piPT conglomerados, se consideró como conglomerados a los municpios, de los cuales se seleccionarán a 45 de estos, al interior de los congrlomerados seleccionados se realiza un censo de todos los colegios. Se quiere estimar el total de evaluados. Se considera como vaariable auxiliar tx el total de colegio por municipios.
En primer lugar se realiza la selección de la muestra de municipios
CONS1 = sqldf("select cod_municipio, municipio, count(*) as 'tx' from saber2013 group by cod_municipio") #Colegios X municipio
set.seed(31052014)
indica_congpipt = S.piPS(45, CONS1$tx)
muecong = CONS1[indica_congpipt[, 1], ]
muecong = cbind(muecong, indica_congpipt[, 2])
names(muecong)[4] = "piIi"
muecong$municipio = NULL
muestra = merge(saber2013, muecong)
A continuación se realiza el proceso de estimación:
CONS2 = sqldf("select cod_municipio, municipio, max(piIi) as 'piIi',\n sum(evaluados) as 'tyi' from muestra group by cod_municipio")
E.piPS(CONS2$tyi, CONS2$piIi)
## N y
## Estimation 965.65 531273.04014321
## Standard Error 147.93 16033.69134979
## CVE 15.32 3.01797572
## DEFF Inf 0.00005938
Más fácil con la librería Survey, aunque utiliza otro algoritmo (“brewer”) para estimar la varianza, por lo tanto va a diferir ligeramente con respectos a las estimaciones anteriormente encontradas:
muestra$UNOS = 1
disenopiptcong = svydesign(ids = ~cod_municipio + codinst, fpc = ~piIi + UNOS,
data = muestra, pps = "brewer")
salida = svytotal(~evaluados, disenopiptcong)
salida
## total SE
## evaluados 531273 24975
Veamos que tan bueno es el diseño muestral, utilizaremos una función basada en TeachingSampling, para calcular que tan es eficiente es el diseño por conglomerados con respecto al MAS.
Varpipt = function(y, Pik, n) {
y <- cbind(1, y)
y <- as.data.frame(y)
ty = sum(y[, 2])
N = length(y[, 2])
names(y)[1] <- "N"
Total <- matrix(NA, nrow = 4, ncol = dim(y)[2])
rownames(Total) = c("Estimation", "Standard Error", "CVE", "DEFF")
colnames(Total) <- names(y)
for (k in 1:dim(y)[2]) {
bk <- (N * Pik * (1 - Pik))/(N - 1)
P1 <- sum(bk * y[, k]/Pik)
P2 <- sum(bk)
ystar <- Pik * P1/P2
P3 <- bk/(Pik^2)
Vty <- sum(P3 * ((y[, k] - ystar)^2))
CVe <- 100 * sqrt(Vty)/ty
N <- sum(1/Pik)
VMAS <- (N^2) * (1 - (n/N)) * var(y[, k])/(n)
DEFF <- Vty/VMAS
Total[, k] <- c(ty, sqrt(Vty), CVe, DEFF)
}
return(Total)
}
Se obtien los siguientes que se muestran a continuación:
CONS3 = sqldf("select cod_municipio, municipio, count(*) as 'txi' ,sum(evaluados) as 'tyi' from saber2013 group by cod_municipio")
pIi = PikPPS(45, CONS3$txi)
Varpipt(CONS3$tyi, pIi, 45)
## N y
## Estimation 553981.0000 553981.0000000000
## Standard Error 220.4603 18158.2719895611
## CVE 0.0398 3.2777788389
## DEFF Inf 0.0000002101
Veamos la relación entre la variable auxiliar total de colegios por municipios y total de evaluados por municpio:
qplot(txi, tyi, data = CONS3)
Existe una relación alta entre estas dos variables y el intercepto no parece ser significativo con respecto a la magnitud de los datos.