library(survey)
library(sae)
library(dplyr)
library(stratification)
library(sampling)
library(TeachingSampling)
library(readr)
library(readxl)
library(pander)
library(ggplot2)
library(knitr)
library(ggthemes)
library(ltm) #correlaciones
library(corrplot)
library(sp)#visualización y gestión archivos vectoriales
library(rgdal) #leer archivos de geoinformación
library(grDevices)
library(openxlsx) #Exportar archivos excel

La razón del uso de SAE:

“Una vez efectuada una encuesta, con tamaños muestrales establecidos para producir estimaciones fiables a un nivel de agregación determinado, a posteriori se produce una demanda de datos a un nivel más desagregado. Para ello, se desea poder utilizar los datos de dicha encuesta sin incurrir en costes adicionales debidos a un incremento de la muestra (. . . ) No obstante, si la encuesta no se había planificado para estimar a un nivel tan desagregado, el tamaño de muestra en algunas de las áreas puede ser demasiado pequeño, lo cual se traduce en errores de muestreo excesivamente grandes para los estimadores directos de las variables de interés en esas áreas” (Molina,2019).

Debido a la necesidad de obtener estimaciones con poco error muestral en áreas pequeñas donde las principales características son observaciones insuficientes o incluso nula presencia de observaciones, se ha planteado algunos métodos que mejoren estas estimaciones, como es el caso de SAE.

En el actual trabajo, se llevará a cabo la estimación del promedio de los puntajes en matemáticas de la prueba Saber 11 en el año 2018 para cada uno de los municipios del departamento de Cundinamarca. Los métodos utilizados para tal fin son: la técnica de estimación directa, más precisamente el estimador Horvitz Thompson y la técnica de SAE, más precisamente el estimador Battise-Harter y Fuller.

Importamos la base poblacional que incluye los resultados en las distintas materias evaluadas en el ICFES para los departamentos del país con sus respectivos municipios, además de algunas variables socio-económicas. De la base compacta se selecciona nuestro departamento de interés, correspondiente a Cundinamarca.

marco_pob<-base[base$DEPARTAMENTO=="CUNDINAMARCA",]
kable(head(marco_pob,6))

	ID_estud	CODIGO_ICFES	AGSB_NOMBREINSTITUCION	CODIGOMUNICIPIO	NOMBREMUNICIPIO	DEPARTAMENTO	CALENDARIO	NATURALEZA	JORNADA	PERS_GENERO	FINS_ESTRATOVIVIENDAENERGIA	FINS_PERSONASHOGARACTUAL	FINS_CUARTOSHOGARACTUAL	FINS_PISOSHOGAR	FINS_TIENEINTERNET	FINS_TIENECOMPUTADOR	FINS_TIENEAUTOMOVILPARTICULAR	LECTURA_CRITICA_PUNT	MATEMATICAS_PUNT	SOCIALES_CIUDADANAS_PUNT	CIENCIAS_NATURALES_PUNT	INGLES_PUNT	INGLES_DESEM	EVALUADOS
44775	ID44775	8201	INST EDUCATIVA DEPTAL GENERAL CARLOS ALBAN	463	ALBAN	CUNDINAMARCA	Calendario A	Oficial	Mañana	F	3	8	4	3	N	S	N	45	48	49	45	36	A-	48
44776	ID44776	8201	INST EDUCATIVA DEPTAL GENERAL CARLOS ALBAN	463	ALBAN	CUNDINAMARCA	Calendario A	Oficial	Mañana	M	2	3	2	2	N	S	N	44	38	42	49	52	A1	48
44777	ID44777	8201	INST EDUCATIVA DEPTAL GENERAL CARLOS ALBAN	463	ALBAN	CUNDINAMARCA	Calendario A	Oficial	Mañana	F	1	3	2	2	N	N	N	33	47	37	40	52	A1	48
44778	ID44778	8201	INST EDUCATIVA DEPTAL GENERAL CARLOS ALBAN	463	ALBAN	CUNDINAMARCA	Calendario A	Oficial	Mañana	M	3	3	2	4	N	S	N	41	43	55	50	47	A-	48
44779	ID44779	8201	INST EDUCATIVA DEPTAL GENERAL CARLOS ALBAN	463	ALBAN	CUNDINAMARCA	Calendario A	Oficial	Mañana	F	1	4	3	4	N	N	S	41	50	44	49	55	A1	48
44780	ID44780	8201	INST EDUCATIVA DEPTAL GENERAL CARLOS ALBAN	463	ALBAN	CUNDINAMARCA	Calendario A	Oficial	Mañana	F	2	5	2	2	S	S	N	56	51	53	45	48	A-	48

Adicionalmente se importa una base con la muestra previamente seleccionada por un diseño de muestreo tri-etapico: Estratificado-Estratificado-MAS

kable(head(muestra,6))

CODIGOMUNICIPIO	CODIGO_ICFES	estrato_mpio	NOMBREMUNICIPIO	EstratoColegio	ID_estud	AGSB_NOMBREINSTITUCION	DEPARTAMENTO	CALENDARIO	NATURALEZA	JORNADA	PERS_GENERO	FINS_ESTRATOVIVIENDAENERGIA	FINS_PERSONASHOGARACTUAL	FINS_CUARTOSHOGARACTUAL	FINS_PISOSHOGAR	FINS_TIENEINTERNET	FINS_TIENECOMPUTADOR	FINS_TIENEAUTOMOVILPARTICULAR	LECTURA_CRITICA_PUNT	MATEMATICAS_PUNT	SOCIALES_CIUDADANAS_PUNT	CIENCIAS_NATURALES_PUNT	INGLES_PUNT	INGLES_DESEM	EVALUADOS	NI	nI	NII	nII	N_i	n_i
00462	038877	Estrato3	AGUA DE DIOS	GRUPO2	ID173720	I.E.D. TECNICO AGROPECUARIO CALANDAIMA - SEDE PRINCIPAL	CUNDINAMARCA	Calendario A	Oficial	Mañana	F	1	6	3	1	N	N	N	35	42	48	44	56	A1	10	32	32	6	3	10	2
00462	038877	Estrato3	AGUA DE DIOS	GRUPO2	ID173721	I.E.D. TECNICO AGROPECUARIO CALANDAIMA - SEDE PRINCIPAL	CUNDINAMARCA	Calendario A	Oficial	Mañana	F	2	4	3	4	N	S	N	42	53	55	52	53	A1	10	32	32	6	3	10	2
00462	150904	Estrato3	AGUA DE DIOS	GRUPO2	ID486361	I.E.D. ENRIQUE SANTOS MONTEJO - SEDE PRINCIPAL	CUNDINAMARCA	Calendario A	Oficial	Sabatina - Dominical	M	2	3	2	3	N	S	S	58	47	60	57	54	A1	51	32	32	6	3	51	11
00462	150904	Estrato3	AGUA DE DIOS	GRUPO2	ID486336	I.E.D. ENRIQUE SANTOS MONTEJO - SEDE PRINCIPAL	CUNDINAMARCA	Calendario A	Oficial	Sabatina - Dominical	F	2	6	3	2	N	N	N	41	49	42	46	47	A-	51	32	32	6	3	51	11
00462	150904	Estrato3	AGUA DE DIOS	GRUPO2	ID486357	I.E.D. ENRIQUE SANTOS MONTEJO - SEDE PRINCIPAL	CUNDINAMARCA	Calendario A	Oficial	Sabatina - Dominical	F	2	11	3	2	S	S	N	40	39	31	32	51	A1	51	32	32	6	3	51	11
00462	150904	Estrato3	AGUA DE DIOS	GRUPO2	ID486344	I.E.D. ENRIQUE SANTOS MONTEJO - SEDE PRINCIPAL	CUNDINAMARCA	Calendario A	Oficial	Sabatina - Dominical	F	2	4	2	3	S	S	N	48	52	57	51	52	A1	51	32	32	6	3	51	11

dim(muestra)

## [1] 2432   32

length(table(muestra$NOMBREMUNICIPIO))

## [1] 42

length(table(marco_pob$NOMBREMUNICIPIO)[table(marco_pob$NOMBREMUNICIPIO)!=0])

## [1] 116

La muestra seleccionada cuenta con 42 municipios de los 116 que se encuentran en el marco poblacional,lo que significa que en la muestra no se observaron 74 municipios de Cundinamarca.

Estimador Directo

A continuación se mostrará las estimaciones directas del promedio del puntaje de matemáticas con su respectivo error estándar y coeficiente de variación para los dominios observados, que en este estudio corresponden a los distintos municipios de Cundinamarca.

options(survey.lonely.psu="adjust")
diseno<- svydesign(ids = ~CODIGOMUNICIPIO + CODIGO_ICFES + ID_estud,
                   strata = ~estrato_mpio + EstratoColegio,
                   fpc = ~ NI + NII + N_i, data = muestra,
                   nest = T)

Estima_Dir <- svyby(~MATEMATICAS_PUNT, by = ~NOMBREMUNICIPIO, 
                    design = diseno, FUN = svymean)

Estima<-data.frame(round(Estima_Dir[,-1],2),"CV"=round(cv(Estima_Dir),4)*100)
kable(head(Estima,8))

	MATEMATICAS_PUNT	se	CV
AGUA DE DIOS	44.63	1.98	4.43
ALBAN	47.13	0.49	1.04
BOJACÁ	51.82	0.49	0.94
CACHIPAY	47.66	0.42	0.89
CAJICA	51.15	2.00	3.92
CAQUEZA	51.65	1.58	3.07
CHIA	56.87	2.40	4.23
CHOACHI	48.33	0.68	1.40

No se tendrá en cuenta los datos faltantes de las variables categoricas porque luego se necesitará cada modalidad de algunas de estas variables para crear el modelo, o bien se hubieran podido imputar con el uso de alguna función.

marco_pob<-marco_pob[!is.na(marco_pob$PERS_GENERO),]
marco_pob<-marco_pob[!is.na(marco_pob$FINS_TIENEINTERNET),]
marco_pob<-marco_pob[!is.na(marco_pob$FINS_TIENECOMPUTADOR),]
muestra<-muestra[!is.na(muestra$PERS_GENERO),]
muestra<-muestra[!is.na(muestra$FINS_TIENECOMPUTADOR),]
muestra<-muestra[!is.na(muestra$FINS_TIENEINTERNET),]

Antes de llevar a cabo las estimaciones para un modelo de Battese-Harter-Fuller (BHF) se determinará primero qué varibles auxiliares usar en el estudio y para ello se mirarán correlaciones con nuestra variable de interés.

corrplot.mixed(cor(marco_pob[,18:22]),tl.col="black")

Se tomarán inicialmente como variables auxiliares candidatas las correspondientes a: - Naturaleza - Genero - Tiene Internet - Tiene computador - Puntaje Lectura Critica - Puntaje Ciencias Naturales

Para medir el nivel de asociación entre el puntaje de matemáticas (Variable Cuantitativa) y las variables auxiliares categoricas se tendrá en cuenta un anális descriptivo donde se mirará si el promedio de puntaje de matematicas difiere ante las diferentes modalidades de las variables categoricas. Adicionalmente, para una mayor seguridad de la asociación entre estas variables, se hará uso de la función “biserial.cor” de la libreria “ltm” que ayuda a calcular la correlación biserial puntual entre una variable dicotómica y una continua.

Variable NATURALEZA

pander(data.frame(group_by(marco_pob,NATURALEZA)%>%summarise(Promedio=mean(MATEMATICAS_PUNT))))

summarise() ungrouping output (override with .groups argument)

NATURALEZA	Promedio
No oficial	53.15
Oficial	49.95

boxplot(marco_pob$MATEMATICAS_PUNT~marco_pob$NATURALEZA,col="turquoise",ylab = "Puntaje Matemáticas",xlab = "NATURALEZA")

Según resultados gráficos, se puede decir que la naturaleza de la institución sí afecta el puntaje de matemáticas, siendo este mayor en los colegios no oficiales. Así, se puede concluir que existe asociación entre estas dos variables.

Variable GÉNERO

pander(data.frame(group_by(marco_pob,PERS_GENERO)%>%summarise(Promedio=mean(MATEMATICAS_PUNT))))

summarise() ungrouping output (override with .groups argument)

PERS_GENERO	Promedio
F	49.69
M	52.46

boxplot(marco_pob$MATEMATICAS_PUNT~marco_pob$PERS_GENERO ,col="deepskyblue",ylab = "Puntaje Matemáticas",xlab = "GÉNERO")

El género Másculino presenta un mejor promedio en los resultados de matemáticas, por lo que se puede decir que el puntaje de esta área difiere de acuerdo al género del evaluado.

Variable Tiene Internet

pander(data.frame(group_by(marco_pob,FINS_TIENEINTERNET)%>%summarise(Promedio=mean(MATEMATICAS_PUNT))))

summarise() ungrouping output (override with .groups argument)

FINS_TIENEINTERNET	Promedio
N	48.72
S	52.94

boxplot(marco_pob$MATEMATICAS_PUN~marco_pob$FINS_TIENEINTERNET,col="azure2",ylab = "Puntaje Matemáticas",xlab = "INTERNET")

Los evaluados que tienen internet en su casa les fue mejor en la prueba de matemáticas que aquellos que no tienen internet y la diferencia de los dos promedios se puede apreciar considerable.

Variable Tiene Computadora

pander(data.frame(group_by(marco_pob, FINS_TIENECOMPUTADOR)%>%summarise(Promedio=mean(MATEMATICAS_PUNT))))

summarise() ungrouping output (override with .groups argument)

FINS_TIENECOMPUTADOR	Promedio
N	47.83
S	52.29

boxplot(marco_pob$MATEMATICAS_PUNT~marco_pob$FINS_TIENECOMPUTADOR ,col="azure1",ylab = "Puntaje Matemáticas",xlab = "COMPUTADOR")

Los evaluados que poseen computador presentan un promedio en la evaluación de matemáticas mayor que aquellos que no tienen, por lo tanto, podría decirse que hay asociación entre estas dos variables.

Ahora, observamos las correlaciones entre nuestra variable de interés y las variables categoricas candidatas:

Paracor=marco_pob %>% mutate(NATURALEZA=recode(NATURALEZA,"No oficial"=0,"Oficial"=1),
                       PERS_GENERO=recode(PERS_GENERO,"F"=0,"M"=1),
                       FINS_TIENEINTERNET=recode(FINS_TIENEINTERNET,"N"=0,"S"=1),
                       FINS_TIENECOMPUTADOR=recode(FINS_TIENECOMPUTADOR,"N"=0,"S"=1))

a=ltm::biserial.cor(x = Paracor$MATEMATICAS_PUNT,y=Paracor$NATURALEZA)
b=ltm::biserial.cor(x = Paracor$MATEMATICAS_PUNT,y=Paracor$PERS_GENERO)
c=ltm::biserial.cor(x = Paracor$MATEMATICAS_PUNT,y=Paracor$FINS_TIENEINTERNET)
d=ltm::biserial.cor(x = Paracor$MATEMATICAS_PUNT,y=Paracor$FINS_TIENECOMPUTADOR)

pander(data.frame(Variable=c("Naturaleza","Genero","T_Internet","T_compu"),
           Correlacion=round(c(a,b,c,d),4)))

Variable	Correlacion
Naturaleza	0.1576
Genero	-0.1471
T_Internet	-0.2243
T_compu	-0.2184

Se observa que existe correlación entre el Puntaje de matemáticas y cada una de las variables categoricas seleccionadas para el estudio.

Ahora se observará el grado de asociación entre el puntaje de matemáticas y los puntajes de ciencias naturales y Lectura crítica respectivamente.

round(cor(marco_pob$MATEMATICAS_PUNT,marco_pob$CIENCIAS_NATURALES_PUNT),2)

## [1] 0.69

round(cor(marco_pob$MATEMATICAS_PUNT,marco_pob$LECTURA_CRITICA_PUNT),2)

## [1] 0.64

Con correlaciones mayores al 0.64 se puede concluir que existe asociaciones significativas entre estas variables

Modelo Battese-Harter-Fuller

Ya con las variables auxiliares elegidas, se procederá a llevar a cabo el modelo.

Información Auxiliar

Observemos el promedio de cada una de nuestras variables auxiliares por municipio.

# Creación de Dummies de nuestras variables auxiliares de naturaleza categorica
Dummies_Natu <- Domains(marco_pob$NATURALEZA)
Dummies_Gen<-Domains(marco_pob$PERS_GENERO)
Dummies_Int<-Domains(marco_pob$FINS_TIENEINTERNET)
Dummies_Compu<-Domains(marco_pob$FINS_TIENECOMPUTADOR)

# Unión a la base y cambio de nombres
marco_pob<-cbind(marco_pob,Dummies_Natu,Dummies_Gen,Dummies_Int,Dummies_Compu)
names(marco_pob)[27:32]<-c("Femenino","Masculino","No_Int","Si_Int","No_Compu","Si_Compu")

Info_aux<-marco_pob %>% group_by(NOMBREMUNICIPIO) %>%
  summarise(Prom_Cien=mean(CIENCIAS_NATURALES_PUNT),
            Prom_Lec=mean(LECTURA_CRITICA_PUNT),
            Prom_mujer=mean(marco_pob$Femenino),prom_hombre=mean(marco_pob$Masculino),
            Prom_Of=mean(Oficial),Prom_Noof=mean(`No oficial`),
            Prom_Noint=mean(No_Int),Prom_Siint=mean(Si_Int),
            Prom_Nocompu=mean(No_Compu),Prom_Sicompu=mean(Si_Compu),
            N_d=n())

kable(as.data.frame(head(Info_aux)))

NOMBREMUNICIPIO	Prom_Cien	Prom_Lec	Prom_mujer	prom_hombre	Prom_Of	Prom_Noof	Prom_Noint	Prom_Siint	Prom_Nocompu	Prom_Sicompu	N_d
AGUA DE DIOS	49.00000	47.45480	0.5480424	0.4519576	1.0000000	0.0000000	0.6158192	0.3841808	0.3728814	0.6271186	354
ALBAN	48.04054	45.74324	0.5480424	0.4519576	1.0000000	0.0000000	0.8513514	0.1486486	0.5945946	0.4054054	74
ANAPOIMA	49.06215	48.29379	0.5480424	0.4519576	0.8757062	0.1242938	0.6101695	0.3898305	0.4350282	0.5649718	177
ANOLAIMA	49.61081	48.82703	0.5480424	0.4519576	1.0000000	0.0000000	0.7567568	0.2432432	0.5081081	0.4918919	185
ARBELAEZ	52.66071	49.75000	0.5480424	0.4519576	1.0000000	0.0000000	0.8750000	0.1250000	0.0892857	0.9107143	112
BELTRAN	43.40000	42.37778	0.5480424	0.4519576	1.0000000	0.0000000	0.8666667	0.1333333	0.5333333	0.4666667	45

Estimación vía modelo BHF

# Convertimos nuestros dominios a caracteres
marco_pob$NOMBREMUNICIPIO<-as.character(marco_pob$NOMBREMUNICIPIO)
muestra$NOMBREMUNICIPIO<-as.character(muestra$NOMBREMUNICIPIO)
Info_aux$NOMBREMUNICIPIO<-as.character(Info_aux$NOMBREMUNICIPIO)

Tamaños<-Info_aux[,c("NOMBREMUNICIPIO","N_d")] #Intercepto del modelo.
Medias<-Info_aux[,c("NOMBREMUNICIPIO","Prom_Cien","Prom_Lec","Prom_mujer","Prom_Noof","Prom_Siint","Prom_Nocompu")]

set.seed(291119)
BHF<-pbmseBHF(MATEMATICAS_PUNT~CIENCIAS_NATURALES_PUNT+LECTURA_CRITICA_PUNT+NATURALEZA+
                 PERS_GENERO+FINS_TIENECOMPUTADOR+FINS_TIENEINTERNET,dom = NOMBREMUNICIPIO,
              meanxpop = Medias,#medias de las variables auxiliares por dominio
              popnsize = Tamaños,#tamaños de los dominios (de las municipios)
              B=100,#100 replicas bootstrap
              data=muestra)

## 
## Bootstrap procedure with B = 100 iterations starts.
## b = 1 
## b = 2 
## b = 3 
## b = 4 
## b = 5 
## b = 6 
## b = 7 
## b = 8 
## b = 9 
## b = 10 
## b = 11 
## b = 12 
## b = 13 
## b = 14 
## b = 15 
## b = 16 
## b = 17 
## b = 18 
## b = 19 
## b = 20 
## b = 21 
## b = 22 
## b = 23 
## b = 24 
## b = 25 
## b = 26 
## b = 27 
## b = 28 
## b = 29 
## b = 30 
## b = 31 
## b = 32 
## b = 33 
## b = 34 
## b = 35 
## b = 36 
## b = 37 
## b = 38 
## b = 39 
## b = 40 
## b = 41 
## b = 42 
## b = 43 
## b = 44 
## b = 45 
## b = 46 
## b = 47 
## b = 48 
## b = 49 
## b = 50 
## b = 51 
## b = 52 
## b = 53 
## b = 54 
## b = 55 
## b = 56 
## b = 57 
## b = 58 
## b = 59 
## b = 60 
## b = 61 
## b = 62 
## b = 63 
## b = 64 
## b = 65 
## b = 66 
## b = 67 
## b = 68 
## b = 69 
## b = 70 
## b = 71 
## b = 72 
## b = 73 
## b = 74 
## b = 75 
## b = 76 
## b = 77 
## b = 78 
## b = 79 
## b = 80 
## b = 81 
## b = 82 
## b = 83 
## b = 84 
## b = 85 
## b = 86 
## b = 87 
## b = 88 
## b = 89 
## b = 90 
## b = 91 
## b = 92 
## b = 93 
## b = 94 
## b = 95 
## b = 96 
## b = 97 
## b = 98 
## b = 99 
## b = 100

Resultados de estimación:

Estim_mate=BHF$est$eblup
names(Estim_mate)[1]="Mpio"
pander(head(Estim_mate))

Mpio	eblup	sampsize
AGUA DE DIOS	47.63	18
ALBAN	47.31	16
BOJACÃ	49.96	16
CACHIPAY	49.04	18
CAJICA	50.89	37
CAQUEZA	50.81	45

Estimación del error cuadrático medio:

ecm<-BHF$mse # varianza estimada de la estimación
names(ecm)[1]<-"Mpio"
pander(head(ecm))

Mpio	mse
AGUA DE DIOS	0.4845
ALBAN	0.605
BOJACÃ	0.8379
CACHIPAY	0.9026
CAJICA	0.3671
CAQUEZA	0.4381

Comparación de técnicas

CV_BHF<-sqrt(BHF$mse$mse) / BHF$est$eblup$eblup * 100 
pander(head(data.frame(Est_Directo=Estima[,1],CV=Estima[,3],Est_BFH=Estim_mate$eblup,CV_BHF=CV_BHF)))

Est_Directo	CV	Est_BFH	CV_BHF
44.63	4.43	47.63	1.462
47.13	1.04	47.31	1.644
51.82	0.94	49.96	1.832
47.66	0.89	49.04	1.937
51.15	3.92	50.89	1.191
51.65	3.07	50.81	1.303

La tabla nos muestra la comparación entre las estimaciones por método directo y vía modelo BFH junto con sus coeficientes de variación.

plot(Estima_Dir$MATEMATICAS_PUNT, type = "n", ylab = "Estimador Directo", xlab = "Municipio", main="Estimador Directo Vs \n Estimador EBLUP" )
points(Estima_Dir$MATEMATICAS_PUNT, type = "l", col = "brown1", lwd = 2, pch = 1, lty = 1)
points(Estim_mate$eblup, type = "l", col = "coral3", lwd = 2, pch = 4, lty = 2)
legend("top", legend = c("Directo", "EBLUP BHF"), ncol = 2,col = c("brown1", "coral3"), lwd = 1.5,pch = c(1, 4), lty = c(2, 2), cex = 0.8)

plot(Estima$CV, type = "n", ylab = "Coeficiente de Variación",xlab = "Municipio", main="CV de estimador directo Vs \n CV de EBLUP BHF")
points(Estima$CV, type = "l", col = "deepskyblue3", lwd = 2, pch = 1, lty = 1)
points(CV_BHF, type = "l", col = "deeppink", lwd = 2, pch = 4, lty = 2)
legend("top", legend = c("Direct", "EBLUP BHF"), ncol= 2, col = c("deepskyblue3", "deeppink"), lwd = 2,pch = c(1, 4), lty = c(1, 2), cex = 0.8)

Se observa que los coeficientes de variación provenientes del estimador BHF son menores que los coeficientes dados por el estimador directo, lo que significa que el modelo BHF muestra buen desempeño para estimar el promedio en áreas pequeñas, en este caso el puntaje de matemáticas en los municipios de Cundimarca.

Algo muy importante y de mucha utilidad de este método es la capacidad que tiende de estimar dominios que no fueron observados en la muestra.

Estimaciones para dominios externos (74 que no salieron en la muestra)

# 1. Extraer los coeficientes (Betas)
Beta_est <- BHF$est$fit$fixed 
names(Beta_est)[1] <-" Intercepto"
Beta_est <- as.matrix(Beta_est)

# 2.matriz de diseño
# Totales por dominio de las zonas que no estan en la muestra
#matriz de diseño sin intercepto:
Xbar_d <- Info_aux[c("Prom_Cien","Prom_Lec","Prom_mujer","Prom_Noof","Prom_Siint","Prom_Nocompu")]
# creo un vector de unos
Unos <- as.data.frame(as.matrix(rep(1, nrow(Info_aux))))
#matriz de diseño con intercepto
Xbar_d <- cbind(Unos, Xbar_d)
Xbar_d <- as.matrix(Xbar_d) # para poder operar x con beta debe ser matriz
rownames(Xbar_d) <- Info_aux$NOMBREMUNICIPIO # coloco el nombre de las zonas a la matriz diseño
#estimación del puntaje promedio de matemáticas en todas las zonas en la población
Prom_dominios <- Xbar_d %*% Beta_est # estima promedio del puntaje de matemáticas por todos los dominios XB
rownames(Prom_dominios) <- Info_aux$NOMBREMUNICIPIO
Prom_dominios <- as.data.frame(Prom_dominios)
Prom_dominios$domain <- row.names(Prom_dominios) # agrega el nombre las zonas como variable
# Prom_dominios= todas las estimaciones para la población
colnames(Prom_dominios)[1] <- "Ybar_efectosfijos" # el mismo beta genera todas las estimaciones (efecto fijo)

# Conservar los dominios no observados
Prom_dominios_observados <- BHF$est$eblup # estimaciones de la muestra
Prom_dominios <- merge(Prom_dominios, Prom_dominios_observados, by = "domain", all.x = T)
names(Prom_dominios)[1] <- "Municipio"

# Se va a crear un modelo mixto, porque cada municipio tiene su propia ecuación. Se van a estimar betas para cada municipio, para ellos se usará un modelo global o jerarquico (mixto o multinivel)


############ Estimación MSE para dominios no observados ###########

#MODELOS MIXTOS, JERARQUICOS O MULTINIVEL
library(nlme)
#EFECTO ALEATORIO, ECUACIÓN EN CADA ZONA, MODELO DADA CADA ZONA (4 BETAS PARA CADA ZONA)
modelo_mixto <- lme(MATEMATICAS_PUNT~CIENCIAS_NATURALES_PUNT+LECTURA_CRITICA_PUNT+NATURALEZA+
                      PERS_GENERO+FINS_TIENECOMPUTADOR+FINS_TIENEINTERNET, random = ~1 | as.factor(NOMBREMUNICIPIO), data = muestra)

Varest_betaest <- as.matrix(vcov(modelo_mixto)) #tiene en cuenta la matriz de varianzas y covarianzas
#entre las variables, esto induce a un beta diferente en cada zona
sigma2est_u <- BHF[[1]]$fit$refvar #EN modelo_mixto observese que es la misma estimacion

# Identificar los dominios no observados
dominios_noobservados<-unique(marco_pob$NOMBREMUNICIPIO)[!unique(marco_pob$NOMBREMUNICIPIO) %in% unique(muestra$NOMBREMUNICIPIO)]

# Matriz diseño de los no observados
Xbar_d_noobs <- Xbar_d[row.names(Xbar_d) %in% dominios_noobservados,]

# SACAMOS LA ESTIMACIÓN DE LA VARIANZA diag(X*sigma*t(x))+sigma_2
MSE_DominiosNoobservados <- diag((Xbar_d_noobs %*% Varest_betaest %*% t(Xbar_d_noobs)) + sigma2est_u)
MSE_DominiosNoobservados <- as.table(MSE_DominiosNoobservados)
df_MSE_DominiosNoobservados <- as.data.frame(MSE_DominiosNoobservados)
names(df_MSE_DominiosNoobservados) <- c("Municipio", "MSE")# zona y varianza estimada
df_MSE_DominiosNoobservados$ClaseDominio <- "No observado" #no observados en la muestra
pander(head(df_MSE_DominiosNoobservados))

Municipio	MSE	ClaseDominio
ANAPOIMA	0.5974	No observado
ANOLAIMA	0.6206	No observado
ARBELAEZ	0.6264	No observado
BELTRAN	0.6511	No observado
BITUIMA	0.6054	No observado
CABRERA	0.6307	No observado

df_MSE_Dominiosobservados <- BHF$mse #varianza de las estimaciones de la muestra
names(df_MSE_Dominiosobservados) <- c("Municipio", "MSE")
df_MSE_Dominiosobservados$ClaseDominio <- "Observado"

#MSE de los no observados se calcula con modelos mixtos y el de los observados
#con el estimador BLUP
df_MSE_Dominios<- bind_rows(df_MSE_DominiosNoobservados, df_MSE_Dominiosobservados)
df_MSE_Dominios <- df_MSE_Dominios[order(df_MSE_Dominios$Municipio),] # ordena por zona

RESULTADOS FINALES

Resultados <- merge(Prom_dominios, df_MSE_Dominios, by = "Municipio")
Resultados$Yhat_BHF <- ifelse(Resultados$ClaseDominio == "No observado", Resultados$Ybar_efectosfijos,
                              Resultados$eblup)
Resultados$cv <- 100 * sqrt(Resultados$MSE) / Resultados$Yhat_BHF
kable(head(Resultados,16))

Municipio	Ybar_efectosfijos	eblup	sampsize	MSE	ClaseDominio	Yhat_BHF	cv
AGUA DE DIOS	48.10548	47.62561	18	0.4845130	Observado	47.62561	1.4615450
ALBAN	47.17531	47.30562	16	0.6049872	Observado	47.30562	1.6442215
ANAPOIMA	48.63469	NA	NA	0.5974466	No observado	48.63469	1.5892909
ANOLAIMA	48.85597	NA	NA	0.6206198	No observado	48.85597	1.6124830
ARBELAEZ	50.30974	NA	NA	0.6263643	No observado	50.30974	1.5731182
BELTRAN	43.87339	NA	NA	0.6511312	No observado	43.87339	1.8392172
BITUIMA	50.24795	NA	NA	0.6054127	No observado	50.24795	1.5484863
BOJACÁ	50.01691	49.96190	16	0.8378934	Observado	49.96190	1.8321265
CABRERA	46.83873	NA	NA	0.6306814	No observado	46.83873	1.6955083
CACHIPAY	50.03724	49.03642	18	0.9025676	Observado	49.03642	1.9374080
CAJICA	51.52515	50.89118	37	0.3671362	Observado	50.89118	1.1906143
CAPARRAPI	48.83053	NA	NA	0.6331502	No observado	48.83053	1.6295284
CAQUEZA	50.21048	50.81384	45	0.4381172	Observado	50.81384	1.3026062
CARMEN DE CARUPA	48.12381	NA	NA	0.6551235	No observado	48.12381	1.6819055
CHAGUANI	47.65721	NA	NA	0.6360147	No observado	47.65721	1.6734201
CHIA	54.83209	55.51768	150	0.2049386	Observado	55.51768	0.8154184

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\Personal\Desktop\Universidad y Estudios\UNIVERSIDAD\Semestre 2019-II\Muestreo II\3ER CORTE\SAE\DANNIS\Municipios_de_Cundinamarca.shp", layer: "Municipios_de_Cundinamarca"
## with 117 features
## It has 7 fields
## Integer64 fields read as strings:  OBJECTID

Res$CODIGO=as.factor(Res$CODIGO) #codigo que permite realizar el mapa
capa@data=left_join(capa@data,Res,by="CODIGO")
names(capa@data)

##  [1] "OBJECTID"   "CODIGO"     "NOMBRE"     "CODIGO_PRO" "RULEID"    
##  [6] "SHAPE_Leng" "SHAPE_Area" "Municipio"  "Yhat_BHF"   "cve"

spplot(capa,"Yhat_BHF",col.regions=heat.colors(n=100,alpha=0.5,rev = T),
       main="Estimación Promedio de Puntaje Matemáticas por BHF")

spplot(capa,"cve",col.regions=terrain.colors(n=100,alpha=0.5,rev = T),
       main="Coeficiente de Variación de la estimación \ndel Promedio de Puntaje Matemáticas por BHF")

Se obtuvieron buenos resultados debido a que tanto el Error Cuadratico Medio como los coeficientes de variación obtuvieron cantidades considerablemente pequeñas.

Comparación con Promedio Real

prom_real<-as.data.frame(group_by(marco_pob,NOMBREMUNICIPIO) %>% summarise(p=mean(MATEMATICAS_PUNT)))
plot(prom_real$p, type = "n",ylab = "Promedio Real", xlab = "Municipio",main="Promedio Real Vs \n Estimado BHF")
points(prom_real$p, type = "l",col = "cyan3", lwd = 2, pch = 1, lty = 1)
points(Resultados$Yhat_BHF, type = "l", col = "darkorchid1", lwd = 2, pch = 4, lty = 2)
legend("top", legend = c("Real", "Estimado_BHF"), ncol = 2, col = c("cyan3","darkorchid1" ), lwd = 2,pch = c(1, 1), lty = c(2, 8),cex=0.8)

Base<-data.frame(BHF=Resultados$Yhat_BHF,Real=prom_real$p,Mpio=Resultados$Municipio)
ggplot(Base,aes(x=BHF,y=Real))+geom_point(color="deepskyblue")+geom_smooth(color="pink")+theme_light()+ggtitle("Real \nVs BHF") +
  theme(plot.title =element_text(hjust=0.5,face="bold",size=15) ) ->g1
g1 + theme_solarized()+ theme(plot.title = element_text(hjust=0.5,face="bold"))->g2
g2

Se observa que las estimaciones vía modelo Battese del promedio en los puntajes de matemáticas, de tanto los municipios observados en la muestra como los no observados, se ajustan bien a los promedios reales.

Estimación en Areas Pequeñas. Modelo Battese-Harter-Fuller (BHF)

Karen Rojas Ramirez

30/09/2020

Estimador Directo

Modelo Battese-Harter-Fuller

Estimación vía modelo BHF

Resultados de estimación:

Estimación del error cuadrático medio:

Comparación de técnicas

Estimaciones para dominios externos (74 que no salieron en la muestra)

RESULTADOS FINALES

Comparación con Promedio Real