Muestreo_Lumey.R

Muestreo_Lumey.R

#-----------------TUTORIAL DE MUESTREOS EN R--------------
R version 3.2.1 (2015-06-18) -- "World-Famous Astronaut"
Copyright (C) 2015 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)
#---------------BASADO EN EL MANUAL DE LUMEY--------------
rm(list=ls())
#Introducción
#El presente tutorial tiene como objetivo principal el introducir al
#usuario en el muestreo simple y complejo en el ambiente de programación
#analítica R.Se cubren los siguientes temas:
#    La selección de tamaños de muestra simples(base y dplyr),
#    Selección sistemática(SamplinUtil),
#    Cálculo de tamaños de muestra(SamplinUtil),
#    Cálculo de intervalos de confianza(DescTools),
#    Selección de muestras con PPT (pps),
#    Muestreo complejo con el paquete survey.
#    Comparación de los resultados arrojados en R con
#información oficial de la Encuesta de Hogares del INEC.


#Ejercicios
#Selección de muestras simples
#Muestreo irrestricto Aleatorio(MIA):
#El R Base contiene la útil función sample en la que se
#pueden obtener muestras aleatorias con o sin reemplazo
#de manera sencilla.

#Vamos a cargar un dataset para ejemplificar y a cuantificar sus dimensiones:

crime<-data.frame(crimtab)
dim(crime)

## [1] 924   3

#Como podemos ver, el dataset crime contiene 924 filas y 3 columnas.
#Seguidamente seleccionamos una muestra de 30 casos sin reemplazo:

#Selección de la muestra

#Tamaño de la muestra
n<-30

muestramia<- sample(1:nrow(crime),size=n,replace=FALSE)
muestramia

##  [1]  64 658 551 887 546 509 639 654 328 659 453 258 279 550 876 676 121
## [18] 170 695 610  55 824 138  96   3 266 256 776 445 638

#Asignar los elementos de la muestra al data frame de datos
crimemuestramia<- crime[muestramia, ]

head(crimemuestramia)

##     Var1   Var2 Freq
## 64  11.5 144.78    0
## 658 12.1 180.34    4
## 551  9.8 175.26    0
## 887  9.8 195.58    0
## 546 13.5 172.72    0
## 509  9.8 172.72    0

#Selección de muestras simples con dplyr
#Una forma más sencilla de obtener una muestra es con el paquete dplyr
#Este paquete es sumamente útil para el tratamiento de datos, adicionalmente
#contiene una función para obtener muestras simples de un data frame:

library(dplyr)

#Muestra sin reemplazo
crimemuestramia2<- crime %>%
  sample_n(size=n,replace=FALSE)

head(crimemuestramia2)

##     Var1   Var2 Freq
## 252 13.5 154.94    0
## 537 12.6 172.72    5
## 203 12.8  152.4    0
## 119 12.8 147.32    0
## 702 12.3 182.88    0
## 649 11.2 180.34    0

#Muestra con pesos
crimemuestramia3<- crime %>%
  sample_n(size=n,weight=Freq)

head(crimemuestramia3)

##     Var1   Var2 Freq
## 314 11.3 160.02   24
## 402 11.7  165.1   37
## 177 10.2  152.4    2
## 616 12.1  177.8   10
## 486 11.7 170.18   45
## 625   13  177.8    1

#Muestra con una proporción de casos
crimemuestramia4<- crime %>%
  sample_frac(0.05)

head(crimemuestramia4);dim(crimemuestramia4)

##     Var1   Var2 Freq
## 260 10.1 157.48    0
## 409 12.4  165.1    2
## 137 10.4 149.86    1
## 473 10.4 170.18    0
## 259   10 157.48    2
## 279   12 157.48    1

## [1] 46  3

#Selección de muestras sistemáticas(Instalación del paquete SamplingUtil)
#Para el ejemplo del muestreo sistemático utilizaremos la función sys.sample
#del paquete SamplingUtil.Para instalar este paquete, se debe inicialmente
#instalar el paquete devtools ya que no se encuentra en CRAN.Las instrucciones
#son las siguientes:
#Instalar devtools:install.packages("devtools")
library(devtools)
#Instalar SamplingUtils:install_github("DFJL/SamplingUtil")
#De esta manera podemos proceder a generar la muestra sistemática:

#Cargar libreria:

library(SamplingUtil)

muestrasis<- sys.sample(N=nrow(crime),n=30)
muestrasis

##  [1]  12  42  72 102 132 162 192 222 252 282 312 342 372 402 432 462 492
## [18] 522 552 582 612 642 672 702 732 762 792 822 852 882

#Asignar los elementos de la muestra al data frame de datos
crimemuestrasis<- crime[muestrasis, ]

head(crimemuestrasis)

##     Var1   Var2 Freq
## 12  10.5 142.24    0
## 42  13.5 142.24    0
## 72  12.3 144.78    0
## 102 11.1 147.32    0
## 132  9.9 149.86    0
## 162 12.9 149.86    0


#Tamaños de muestra simples

#En este tutorial veremos tres variantes para cálculo de tamaño muestral
#asumiendo MIA.Para esto, utilizaremos la data de la ENAHO 2013 y la función
#nsize() incluída en el paquete SamplingUtil.

#Como variable continua se utilizará la Cantidad de personas del hogar (TamHog),
#para el ejemplo de proporciones, se creará la variable Proporción de Hogares
#Unipersonales.

#Error relativo para variables continuas.
#se desea conocer el tamaño de muestra para que el error relativo no supere el 5%,
#con un alpha de 5%.

#Generando df de ENAHO2013 a nivel de hogar

df<- ENAHO2013 %>%  #Define nuevo data frame
  mutate(TamHog=as.numeric(TamHog), #Convierte TamHog a numérica
         phu=ifelse(TamHog>1,0,1)) %>% #Crea variable de prop Hog. Unipersonales.
  select(FACTOR:ZONA,TamHog,phu,-REGION) %>% #Selecciona variables de interés.
  group_by(SEGMENTO,CUESTIONARIO,HOGAR,ZONA) %>% #Genera esquema de agrupación.
  summarise(TamHog=mean(TamHog), #Cálculo de variables sumarizadas.
            phu=mean(phu)) %>%
  mutate(id=paste0(SEGMENTO,CUESTIONARIO,HOGAR,ZONA)) #Genera id único para uso posterior.


#Cálculo del tamaño de muestra con error relativo

#Error relativo
r<-0.05

nsizeR<- nsize(x=df$TamHog,r=r,alpha=0.05)

nsizeR

## $n0
## [1] 782
## 
## $n_adjust
## [1] 731

#------Error absoluto para variables continuas.
#Cálculo del tamaño de muestra con error absoluto
#Igualar el error absoluto con el error relativo
abs<-mean(df$TamHog)*r
nsizeA<- nsize(x=df$TamHog,abs=abs,alpha=0.05)
nsizeA

## $n0
## [1] 782
## 
## $n_adjust
## [1] 731

#------Error relativo para variables dicotómicas.
#Cálculo de tamaño de muestra para proporción
nsizeP<- nsize(x=df$phu,abs=0.02,alpha=0.05)

nsizeP

## $n0
## [1] 924
## 
## $n_adjust
## [1] 854

#---------Tamaños de muestra estratificados
#-----Asignación proporcional.

#Cálculo de las proporciones por estrato

Estratos<- df %>%
  select(ZONA,TamHog) %>%
  group_by(ZONA) %>%
  summarise(n=n(),
            s=sd(TamHog)) %>%
  mutate(p=n/sum(n))

#-----------------Estratos

## Source: local data frame [2 x 4]
## 
##     ZONA    n        s         p
## 1 Urbana 4790 3.149710 0.4269543
## 2  Rural 6429 1.945832 0.5730457

#Asignación de la muestra proporcional a los estratos
nsizeProp<-nstrata(n=nsizeR[[2]],wh=Estratos[,4],method="proportional")

nsizeProp

##     p
## 1 313
## 2 419

#---------------Asignación óptima.

#Asignación de la muestra óptima a los estratos

#Costo de entrevista por estrato
ch<-c(5,10)

nsizeOpt<-nstrata(n=nsizeR[[2]],wh=Estratos[,4],sh=Estratos[,3],ch,method="optimal")
nsizeOpt

##     p
## 1 461
## 2 271

#-------------Asignación de Neyman.

#Asignación de la muestra óptima a los estratos(asume costos iguales)

nsizeNeyman<-nstrata(n=nsizeR[[2]],wh=Estratos[,4],sh=Estratos[,3],method="neyman")
nsizeNeyman

##     p
## 1 400
## 2 332

#Cálculo de intervalos de confianza asumiendo MIA

#En R existen múltiples posibilidades de calcular los intervalos de
#confianza.Se utilizará el paquete DescTools para este propósito.Este 
#paquete contiene muchas funciones misceláneas de las tareas estadísticas típicas.

#Calcularemos el intervalo de confianza para la media del Tamaño del hogar,
#con el tamaño de muestra obtenida en el ejercicio anterior y posteriormente
#se comparará el error obtenido(suponiendo que el valor verdadero es el de total
#de la encuesta)

#Muestra sin reemplazo
muestra<- data.frame(df) %>%
  sample_n(size=nsizeR[[2]],replace=FALSE)

crimemuestramia2<- crime %>%
  sample_n(size=n,replace=FALSE)

#Carga de paquete DescTools

library(DescTools)

#Intervalos de confianza 95% clásicos

ICTamHog<- MeanCI(x=muestra$TamHog, trim = 0,conf.level = 0.95, na.rm = FALSE)

ICTamHog

##     mean   lwr.ci   upr.ci 
## 3.393981 3.277623 3.510339

#Media del total de la encuesta
mean(df$TamHog)

## [1] 3.547134

#Diferencia relativa

difR<-paste0(abs(round((ICTamHog[1]-mean(df$TamHog))/mean(df$TamHog),3))*100,"%")


#De esta manera, con este ejercicio comprobamos que el valor estimado cumple con
#los requisitos del error relativo deseado en el cálculo del tamaño muestral,
#ya que la diferencia fue de 4.3% inferior al 5 % deseado.
#Selección de muestra con ppt(Paquete pps)

#Hemos abarcado aspectos básicos del muestreo de elementos.A partir de ahora
#veremos algunos ejemplos de muestreo complejo.

#Existen diversos paquetes para selecciones de muestras complejas en R, como el
#sampling o el samplingbook.En este caso se utilizará el pps que se concentra en
#muestreo mediante ppt y permite la selección sistemática estratificada.

#En el siguiente ejercicio se seleccionará una muestra de UPM utilizando como
#tamaño la cantidad de hogares, estratificada por Zona,utilizando como base
#el tamaño de muestra ya calculado para la variable Tamaño del hogar.

#Se debe obtener una muestra proporcional para cada estrato y seleccionar
#en cada conglomerado una submuestra de 5 hogares.
#Generar el dataframe de UPMS y su respectivo tamaño

UPMs<- df %>%
  mutate(TamHog=as.numeric(TamHog)) %>%
  select(SEGMENTO,ZONA,TamHog) %>%
  group_by(ZONA,SEGMENTO) %>%
  summarise(Mta=n()) %>%
  arrange(desc(ZONA),SEGMENTO)

head(UPMs)

## Source: local data frame [6 x 3]
## Groups: ZONA
## 
##     ZONA SEGMENTO Mta
## 1 Urbana        2   7
## 2 Urbana        4   9
## 3 Urbana        5   9
## 4 Urbana        8  13
## 5 Urbana        9  10
## 6 Urbana       10   9

#Se asigna proporcionalmente la muestra a cada estrato

#Recordemos la distribución proporcional de los estratos:
nsizeProp

##     p
## 1 313
## 2 419

#Cálculo del número de UPMs por estrato
b<-5
aurbano<- ceiling(nsizeProp[1,1]/b)
arural<- ceiling(nsizeProp[2,1]/b)

#Unir en un solo objeto para facilitar input de función
asizeppt<-rbind(aurbano,arural)
asizeppt

##         [,1]
## aurbano   63
## arural    84

#Selección de la muestra con PPT(sistemática)
library(pps)
muestrappt<-ppssstrat(sizes=UPMs$Mta,stratum=UPMs$ZONA,n=asizeppt[,1])
muestraUPMs<-UPMs[muestrappt,]

#Muestra de la selección
head(muestraUPMs)

## Source: local data frame [6 x 3]
## Groups: ZONA
## 
##     ZONA SEGMENTO Mta
## 1 Urbana       10   9
## 2 Urbana       25   8
## 3 Urbana       44   9
## 4 Urbana       58  13
## 5 Urbana       82  11
## 6 Urbana       97   9

tail(muestraUPMs)

## Source: local data frame [6 x 3]
## Groups: ZONA
## 
##    ZONA SEGMENTO Mta
## 1 Rural     1083  10
## 2 Rural     1101  10
## 3 Rural     1114  14
## 4 Rural     1125  10
## 5 Rural     1135  11
## 6 Rural     1146  11

#Verificar selección por estrato
Freq(muestraUPMs$ZONA)

##    level freq  perc cumfreq cumperc
## 1 Urbana   63 0.429      63   0.429
## 2  Rural   84 0.571     147   1.000

#Procedimiento para seleccionar submuestras en cada cluster de tamaño fijo

#Listado de segmentos seleccionados
segs<- unique(muestraUPMs$SEGMENTO)
sample<- list()

#Identifica el número de columna del id(se requiere para el ciclo)

idcolnum<-which( colnames(df)=="id" )

#Ciclo para seleccionar la muestra en cada segmento
for(i in 1:length(segs)){
  subsample<-data.frame(df[which(df$SEGMENTO==segs[i]),])
  sample[[i]]<-sys.sample(nrow(subsample),b)
 #Previene de los números que se pasan del total del segmento
  #if(sample[[i]][b]>nrow(subsample)){
  #  sample[[i]][b]<-1
  #  }
 #Asigna en cada elemento de las submuestras los elem. seleccionados
 sample[[i]]<-subsample[unlist(sample[[i]]),idcolnum]
 }

#Genera el data frame de ids seleccionados(ya que estaban en una lista)
sampledf<-data.frame(id=unlist(sample))

#Uniendo el data frame de datos con la muestra seleccionada mendiante la llave creada
muestrappt<-inner_join(unique(sampledf),unique(df),by="id")

#Verificar el procedimiento
head(muestrappt)

##           id SEGMENTO CUESTIONARIO HOGAR   ZONA TamHog phu
## 1  3111Rural        3           11     1  Rural      6   0
## 2 1011Urbana       10            1     1 Urbana      2   0
## 3 1021Urbana       10            2     1 Urbana      6   0
## 4 1031Urbana       10            3     1 Urbana      2   0
## 5 1041Urbana       10            4     1 Urbana      5   0
## 6 1051Urbana       10            5     1 Urbana      3   0

Freq(muestrappt$ZONA)

##    level freq  perc cumfreq cumperc
## 1 Urbana  317 0.428     317   0.428
## 2  Rural  424 0.572     741   1.000

#------------------------------------Diseños de muestreo Complejos
#---------Especificación diseños muestrales complejos(paquete survey)

#Un paquete robusto para el cálculo correcto de estadísticos en diseños complejos
#es el survey.A continuación se ejemplifica la función svydesign, la cuál establece
#los diferentes parámetros para poder calcular los errores muestrales adecuadamente.
#Inicialmente utilizaremos la muestra con ppt seleccionada en el ejercicio anterior.

#Generar probabilidades

#Recordemos las medidas de tamaño por UPM:
UPMs

## Source: local data frame [1,119 x 3]
## Groups: ZONA
## 
##      ZONA SEGMENTO Mta
## 1  Urbana        2   7
## 2  Urbana        4   9
## 3  Urbana        5   9
## 4  Urbana        8  13
## 5  Urbana        9  10
## 6  Urbana       10   9
## 7  Urbana       11  12
## 8  Urbana       12   9
## 9  Urbana       16   6
## 10 Urbana       18   8
## ..    ...      ... ...

#Cálculo de prob en Primera etapa

#Estrato Urbano

UPMU<- UPMs %>%
  filter(ZONA=="Urbana") %>%
  group_by(ZONA) %>%
  mutate(fa=Mta/(sum(Mta)/aurbano),
         fb=b/Mta)

#Estrato Rural

UPMR<- UPMs %>%
  filter(ZONA=="Rural") %>%
  group_by(ZONA) %>%
  mutate(fa=Mta/(sum(Mta)/arural),
         fb=b/Mta)

#Unir ambos estratos

UPMs <- rbind(UPMU,UPMR)

#Incluir probabilidades en data frame de muestra
muestrappt2<- left_join(muestrappt,UPMs,by=c("SEGMENTO","ZONA"))


#Cargar paquete de muestreo complejo
library(survey)

#Establecer el diseño muestral
design1<-svydesign(data=muestrappt2,id=~SEGMENTO+id,
                   strata=~ZONA,pps="brewer",prob=~fa+fb)

#Resumen del diseño
summary(design1)

## Stratified 2 - level Cluster Sampling design (with replacement)
## With (151, 741) clusters.
## svydesign(data = muestrappt2, id = ~SEGMENTO + id, strata = ~ZONA, 
##     pps = "brewer", prob = ~fa + fb)
## Probabilities:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.06533 0.06533 0.06533 0.06551 0.06576 0.06576 
## First-level Stratum Sizes: 
##            Rural Urbana
## obs          424    317
## design.PSU    88     63
## actual.PSU    88     63
## Data variables:
##  [1] "SEGMENTO"     "ZONA"         "id"           "CUESTIONARIO"
##  [5] "HOGAR"        "TamHog"       "phu"          "Mta"         
##  [9] "fa"           "fb"

#---------------------Estadísticos en diseños muestrales complejos

#Calculando estadísticos de TamHog

TamHogmc<- svymean(~TamHog,design=design1,deff=TRUE)

#Media del Tamaño del hogar
TamHogmc

##           mean      SE   DEff
## TamHog 3.45115 0.10447 1.1018

#Efecto del diseño
deff(TamHogmc)

##   TamHog 
## 1.101787

#Intervalo de confianza
confint(TamHogmc)

##           2.5 %   97.5 %
## TamHog 3.246387 3.655921

#Coeficiente de variación
cv(TamHogmc)

##            TamHog
## TamHog 0.03027248

#Diferencia relativa

difRc<-paste0(abs(round((TamHogmc[1]-mean(df$TamHog))/mean(df$TamHog),3))*100,"%")
difRc
#De esta manera, con este ejercicio comprobamos que el valor estimado cumple
#con los requisitos del error relativo deseado en el cálculo del tamaño muestral,
#ya que la diferencia fue de 2.7% inferior al 5 % deseado.
#Ejemplo con diseño complejo real

#Finalmente, se desarrollará un ejemplo con un diseño complejo real,utilizando la ENAHO 2012.

#El diseño muestral se detalla a continuación:
#Diseño muestral: diseño probabilístico, estratificado y bietápico.
#En la primera etapa se seleccionaron segmentos censales o unidades primarias de muestreo (UPM)
#con probabilidad proporcional al tamaño (PPT), y en la segunda etapa se seleccionaron las
#viviendas o unidades secundarias de muestreo (USM) con probabilidades iguales de selección
#dentro de cada segmento, mediante muestreo sistemático con arranque aleatorio.

#Dominios de estudio: las seis regiones de planificación y las zonas urbana y rural.

#Objetos para generar variables(debido a que el Dataframe de ENAHO2012 tiene las labels del SPSS y
# no los valores)

profesionales<- c("Profesionales científicos e intelectuales")
tecnicos<- c("Técnicos y profesionales de nivel medio")

#Generando df de ENAHO2012 a nivel de hogar

df12<- ENAHO2012 %>%  #Define nuevo data frame
  mutate(pobre=ifelse(np=="No pobre",0,1), #Crea variable de pobreza.
         Profesional=ifelse(OcupFuerzaTrab==profesionales,1,0), #Crea flag de Profesional.
         Tecnico=ifelse(OcupFuerzaTrab==tecnicos,1,0)) %>%  #Crea flag de tecnico.
  select(FACTOR:ZONA,pobre,Profesional,Tecnico,ipcn) %>% #Selecciona variables de interés.
  group_by(SEGMENTO,CUESTIONARIO,HOGAR,ZONA,REGION) %>% #Genera esquema de agrupación.
  summarise(Factor=mean(FACTOR),
            pobre=mean(pobre),
            Profesionales=sum(Profesional),
            Tecnicos=sum(Tecnico),
            Ingreso=mean(ipcn)) %>%
  mutate(id=paste0(SEGMENTO,CUESTIONARIO,HOGAR,ZONA,REGION)) #Genera id único para uso posterior.

#Establecer el diseño muestral

design2<-svydesign(data=df12,id=~SEGMENTO+id,
                   strata=~REGION+ZONA,pps="brewer",weights=~Factor)

summary(design2)

## Stratified 2 - level Cluster Sampling design (with replacement)
## With (1120, 11374) clusters.
## svydesign(data = df12, id = ~SEGMENTO + id, strata = ~REGION + 
##     ZONA, pps = "brewer", weights = ~Factor)
## Probabilities:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.001403 0.008065 0.015380 0.016030 0.021280 0.062500 
## First-level Stratum Sizes: 
##            Central Chorotega Pacífico Central Brunca Huetar Atlántica
## obs           4940      1220             1144   1598             1218
## design.PSU     468       116              124    160              128
## actual.PSU     468       116              124    160              128
##            Huetar Norte
## obs                1254
## design.PSU          124
## actual.PSU          124
## Data variables:
##  [1] "SEGMENTO"      "CUESTIONARIO"  "HOGAR"         "ZONA"         
##  [5] "REGION"        "Factor"        "pobre"         "Profesionales"
##  [9] "Tecnicos"      "Ingreso"       "id"

#Calculando estadísticos de Pobreza por hogar nacionales

#Porcentaje de pobreza
pobremc<- svymean(~pobre,design=design2,deff=TRUE,na.rm=TRUE)

#Media de pobreza
pobremc

##            mean        SE  DEff
## pobre 0.2067336 0.0063786 2.836

#Efecto del diseño
deff(pobremc)

##    pobre 
## 2.835952

#Intervalo de confianza
confint(pobremc)

##           2.5 %    97.5 %
## pobre 0.1942318 0.2192355

#Coeficiente de variación
cv(pobremc)

##           pobre
## pobre 0.0308542

#Total de pobres
pobretmc<-svytotal(~pobre,design=design2, deff=TRUE,na.rm=TRUE)

pobretmc

##        total     SE   DEff
## pobre 279514  10334 4.0715

#Efecto del diseño
deff(pobretmc)

##    pobre 
## 4.071516

#Intervalo de confianza
confint(pobretmc)

##          2.5 %   97.5 %
## pobre 259260.8 299767.2

#Coeficiente de variación
cv(pobretmc)

##            pobre
## pobre 0.03696945

#Estimación de quantiles(Ingreso per cápita)
ingresoq<- svyquantile(~Ingreso,quantile= seq(0.2,1,0.2),design=design2, deff=TRUE,na.rm=TRUE)

ingresoq

##             0.2      0.4      0.6    0.8        1
## Ingreso 85889.4 149602.9 245042.2 433385 13059042

#Cálculo de razón(Profesionales/Técnicos)
ptratio<- svyratio(~Profesionales, ~Tecnicos, design=design2, deff=TRUE,na.rm=TRUE)

ptratio

## Ratio estimator: svyratio.survey.design2(~Profesionales, ~Tecnicos, design = design2, 
##     deff = TRUE, na.rm = TRUE)
## Ratios=
##               Tecnicos
## Profesionales 1.257501
## SEs=
##                Tecnicos
## Profesionales 0.1431679

#Efecto del diseño
deff(ptratio)

## [1] 2.278801

#Intervalo de confianza
confint(ptratio)

##                            2.5 %   97.5 %
## Profesionales/Tecnicos 0.9768973 1.538105

#Coeficiente de variación
cv(ptratio)

##                Tecnicos
## Profesionales 0.1138511

#Estimaciones por estrato(Zona)

pobremZ<-svyby(~pobre, ~ZONA, design2, svymean,deff=TRUE,na.rm=TRUE)

#Media de pobreza por zona
pobremZ

##          ZONA     pobre          se DEff.pobre
## Urbana Urbana 0.1834762 0.008193071   2.200074
## Rural   Rural 0.2714831 0.007898585   2.071340

#Efecto del diseño
deff(pobremZ)

## [1] 2.200074 2.071340

#Intervalo de confianza
confint(pobremZ)

##            2.5 %    97.5 %
## Urbana 0.1674180 0.1995343
## Rural  0.2560022 0.2869640

#Coeficiente de variación
cv(pobremZ)

##     Urbana      Rural 
## 0.04465469 0.02909420

#Estimaciones por estrato(Región)

pobremR<-svyby(~pobre,~REGION, design2, svymean,deff=TRUE,na.rm=TRUE)

#Media de pobreza por zona
pobremR

##                            REGION     pobre          se DEff.pobre
## Central                   Central 0.1613506 0.008443248   2.602328
## Chorotega               Chorotega 0.3356425 0.018022067   1.795762
## Pacífico Central Pacífico Central 0.2688588 0.022885659   3.090166
## Brunca                     Brunca 0.3332134 0.021407389   3.339586
## Huetar Atlántica Huetar Atlántica 0.2480790 0.020845431   2.861739
## Huetar Norte         Huetar Norte 0.2320398 0.017439694   2.156538

#Efecto del diseño
deff(pobremR)

## [1] 2.602328 1.795762 3.090166 3.339586 2.861739 2.156538

#Intervalo de confianza
confint(pobremR)

##                      2.5 %    97.5 %
## Central          0.1448021 0.1778990
## Chorotega        0.3003199 0.3709651
## Pacífico Central 0.2240038 0.3137139
## Brunca           0.2912557 0.3751712
## Huetar Atlántica 0.2072227 0.2889353
## Huetar Norte     0.1978586 0.2662210

#Coeficiente de variación
cv(pobremR)

##          Central        Chorotega Pacífico Central           Brunca 
##       0.05232859       0.05369424       0.08512146       0.06424527 
## Huetar Atlántica     Huetar Norte 
##       0.08402738       0.07515820

#Estimaciones combinadas

pobremZR<-svyby(~pobre,~REGION+ZONA, design2, svymean,deff=TRUE,na.rm=TRUE)
pobremZR