MicroDatosEs, EPA y modelos mixtos
José Luis Cañadas Reche
Junio de 2018
Librerías
list.of.packages <- c("MicroDatosEs","rgdal", "lme4","RCurl", "dplyr", "sjPlot", "ggplot2", "broom", "classInt", "reshape2")
new.packages <-
list.of.packages[!(list.of.packages %in% installed.packages()[, "Package"])]
if (length(new.packages))
install.packages(new.packages)
# Importar libreria
lapply(list.of.packages,
require,
character.only = TRUE,
quietly = TRUE)## rgdal: version: 1.3-1, (SVN revision 747)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: /usr/share/gdal/2.2
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: (autodetected)
## Linking to sp version: 1.2-7##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source'))`## [[1]]
## [1] TRUE
##
## [[2]]
## [1] TRUE
##
## [[3]]
## [1] TRUE
##
## [[4]]
## [1] TRUE
##
## [[5]]
## [1] TRUE
##
## [[6]]
## [1] TRUE
##
## [[7]]
## [1] TRUE
##
## [[8]]
## [1] TRUE
##
## [[9]]
## [1] TRUE
##
## [[10]]
## [1] TRUEDescargamos la EPA
¡ El INE tiene un ftp accesible!
file_location <- "ftp://www.ine.es/temas/epa/datos_1t18.zip"
td <- tempdir()
tf <- tempfile(tmpdir = td, fileext = ".zip")
# descargar el fichero
download.file(file_location, tf)
# Nombre del fichero descomprimido
fname <- unzip(tf, list = TRUE)$Name[1]
# Descomprimir el fichero en directorio temporal
unzip(tf,
files = fname,
exdir = td,
overwrite = TRUE)
# fpath is the full path to the extracted file
fpath <- file.path(td, fname)Leemos la EPA usando la función epa2005del paquete MicroDatosEs
epa <- epa2005(fpath)
names(epa) <- tolower(names(epa))Seleccionamos solo las variables que nos interesan
epa <- subset(epa, select = c(prov, edad, nforma, aoi, factorel))
DT::datatable(epa[1:100,])# Recode data
recodificacion <- function (dat) {
dat$aoi[grepl("^Inactivos", dat$aoi)] <- "i"
dat$aoi[grepl("[O-o]cupados", dat$aoi)] <- "o"
dat$aoi[grepl("^Parados", dat$aoi)] <- "p"
dat$aoi <- factor(dat$aoi)
dat$nforma3 <- dat$nforma
dat$nforma3[dat$nforma == "Analfabetos" |
dat$nforma == "Educación primaria" |
dat$nforma == "Educación primaria incompleta"] <-
"Est primarios o menos"
dat$nforma3[dat$nforma == "Educación superior"] <-
"Est. Universitarios"
dat$nforma3[dat$nforma == "Primera etapa de educación secundaria" |
dat$nforma == "Segunda etapa de educación secundaria, orientación general" |
dat$nforma == "Segunda etapa de educación secundaria, orientación profesional"] <-
"Est. Secundarios"
dat$nforma3 <- factor(dat$nforma3)
dat$gedad <- dat$edad
dat$gedad[dat$edad == "de 0 A 4 años" |
dat$edad == "de 5 A 9 años" |
dat$edad == "de 10 A 15 años"] <- "15 años o menos "
dat$gedad[dat$edad == "de 16 A 19 años" |
dat$edad == "de 20 A 24 años" |
dat$edad == "de 25 A 29 años" |
dat$edad == "de 30 A 34 años"] <- "De 16 a 34"
dat$gedad[dat$edad == "de 35 A 39 años" |
dat$edad == "de 40 A 44 años" |
dat$edad == "de 45 A 49 años" |
dat$edad == "de 50 A 54 años"] <- "De 35 a 54"
dat$gedad[dat$edad == "de 55 A 59 años" |
dat$edad == "de 60 A 64 años" |
dat$edad == "65 o más años"] <- "55 o más"
dat$gedad <-
factor(dat$gedad,
levels = c("15 años o menos ", "De 16 a 34", "De 35 a 54", "55 o más"))
dat
}epa <- recodificacion(epa)
head(epa)## # A tibble: 6 x 7
## prov edad nforma aoi factorel nforma3 gedad
## <chr> <chr> <chr> <fct> <dbl> <fct> <fct>
## 1 Álava de 40 A 44 años Educación superior o 22223 Est. Uni… De 35…
## 2 Álava de 40 A 44 años Educación superior o 22223 Est. Uni… De 35…
## 3 Álava de 35 A 39 años Educación superior o 23566 Est. Uni… De 35…
## 4 Álava de 35 A 39 años Educación primaria i 23566 Est prim… De 35…
## 5 Álava de 5 A 9 años <NA> <NA> 22461 <NA> "15 a…
## 6 Álava de 5 A 9 años <NA> <NA> 22461 <NA> "15 a…Añadimos el código de provincia para poder unir posteriormente con un shapefile.
En los ficheros internos del paquete MicroDatosEs vienen metadatos útiles que nos sirven para este propósito.
metadatos_epa <- system.file("metadata", "epa_mdat2.txt",
package = "MicroDatosEs")
epa_mdat2 <- read.table(metadatos_epa, header = T, sep = "\t",
fileEncoding = "UTF-8", stringsAsFactors = FALSE)
head(epa_mdat2)## var tipo nulo llave valor
## 1 CICLO N NA <NA>
## 2 CCAA D NA 1 Andalucía
## 3 CCAA D NA 2 Aragón
## 4 CCAA D NA 3 Asturias (Principado de)
## 5 CCAA D NA 4 Baleares (Islas)
## 6 CCAA D NA 5 Canariasepa_mdat2 <- epa_mdat2[epa_mdat2$var=="PROV",]
epa <- epa %>%
left_join(epa_mdat2[, c("llave","valor")], by = c("prov" = "valor"))
epa## # A tibble: 160,875 x 8
## prov edad nforma aoi factorel nforma3 gedad llave
## <chr> <chr> <chr> <fct> <dbl> <fct> <fct> <chr>
## 1 Álava de 40 A 44 años Educació… o 22223 Est. Univ… De 35 … 1
## 2 Álava de 40 A 44 años Educació… o 22223 Est. Univ… De 35 … 1
## 3 Álava de 35 A 39 años Educació… o 23566 Est. Univ… De 35 … 1
## 4 Álava de 35 A 39 años Educació… i 23566 Est prima… De 35 … 1
## 5 Álava de 5 A 9 años <NA> <NA> 22461 <NA> "15 añ… 1
## 6 Álava de 5 A 9 años <NA> <NA> 22461 <NA> "15 añ… 1
## 7 Álava de 40 A 44 años Educació… o 22413 Est. Univ… De 35 … 1
## 8 Álava de 40 A 44 años Educació… o 22413 Est. Univ… De 35 … 1
## 9 Álava de 0 A 4 años <NA> <NA> 21084 <NA> "15 añ… 1
## 10 Álava de 40 A 44 años Educació… o 26193 Est. Univ… De 35 … 1
## # ... with 160,865 more rowsCálculo de la tasa de paro
Eliminamos los menores de 16 años y los inactivos.
# creo variable numérica para las provincias
epa$prov <- factor(epa$prov)
epa$prov.n <- epa$llave
epa$prov.n <- factor(ifelse(nchar(as.character(epa$prov.n))==1,
paste0("0",epa$prov.n),
as.character(epa$prov.n)))
epa$llave <- NULL
# eliminar menores de 16 años
epa <- epa[epa$gedad!= "15 años o menos ", ]
# eliminar inactivos
epa <- epa[ epa$aoi != "i", ]
# eliminar niveles que no se usan
epa$gedad <- droplevels( epa$gedad)
levels(epa$gedad)## [1] "De 16 a 34" "De 35 a 54" "55 o más"Modelo de efectos aleatorios
La forma usual de calcular la tasa de paro es dividiendo el número de personas en paro entre el total de personas activas. Si hacemos esto teniendo en cuenta la provincia, la formación y el grupo de edad se obtienen resultados no creíbles. Esto es debido a la falta de personas (en la EPA) en esas categorías.
epa %>%
group_by(prov, prov.n, gedad, nforma3) %>%
summarise(
n = n(),
n_parado = sum(aoi=="p"),
tasa_estim_directa = weighted.mean(aoi == "p", factorel)
) %>% arrange(n, n_parado) %>%
DT::datatable()¿Nos creemos que la tasa de paro en Zamora en el grupo de edad de 16 a 34 años de edad son estudios primarios es del 16% y en Alava en ese mismo grupo de edad es del 100%, cuándo en la encuesta sólo hay 1 parado sobre 3 encuestados en Zamora y 2 parados sobre 2 encuestados en Álava?
Los modelos mixtos permiten hacer “partial pooling” compartiendo información entre los grupos, y obtener estimaciones utilizando el estimador BLUP). Un excelente libro para aprender más sobre esto y otras cosas es uno de Gelman enlace
Ajustamos un modelo mixto muy sencillo, dónde consideramos como efectos aleatorios la provincia, el grupo de edad y la educación.
fit_mixto <- glmer(aoi == "p" ~ (1 | prov.n) + (1 | gedad) + (1 | nforma3), family = binomial, data = epa)epa$tasa_estim_mixto <- fitted(fit_mixto)tasa.paro <- epa %>% group_by(prov, prov.n, gedad, nforma3) %>%
summarise(
n = n(),
n_parado = sum(aoi=="p"),
tasa_estim_directa = weighted.mean(aoi == "p", factorel),
tasa_estim_mixto = mean(tasa_estim_mixto)) %>% ungroup()
names(tasa.paro)[1:2] <- c("provincia","prov")
tasa.paro %>%
arrange(n, n_parado) %>%
DT::datatable()Gráfico de los efectos aleatorios del modelo
plot_model(fit_mixto, type = "re", grid = FALSE, sort = "sort.all")## [[1]]
##
## [[2]]
##
## [[3]]
p <- ggplot(tasa.paro, aes( x = tasa_estim_mixto, fill = nforma3))
p + geom_density(alpha = 0.4) +facet_grid(~gedad) 
Extra. Mapa
Lectura shapefile. Se puede descargar el shapefile del ine, del paquete SIANE o de mi dropbox público. En este último caso, me he permitido la licencia de acercar un poco las Islas Canarias. enlace
# usar ruta absoluta en readOGR
prov.map <- readOGR(dsn= "/home/jose/Dropbox/Public/shapefile_spain/spain_provinces_ind_2.shp",layer="spain_provinces_ind_2")## OGR data source with driver: ESRI Shapefile
## Source: "/home/jose/Dropbox/Public/shapefile_spain/spain_provinces_ind_2.shp", layer: "spain_provinces_ind_2"
## with 116 features
## It has 10 fieldsprov.map$OBJECTID## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
## [18] 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
## [35] 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
## [52] 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
## [69] 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85
## [86] 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102
## [103] 103 104 105 106 107 108 109 110 111 112 113 114 115 116prov.map <- prov.map[!is.na(prov.map$NOMBRE99),]
head(prov.map@data, n = 2)## OBJECTID AREA PERIMETER P20099_ P20099_ID NOMBRE99 Shape_Leng
## 0 1 7980747647 1032580.0 2 1 Coruña (A) 1032580.0
## 1 2 1979346035 284275.4 3 12 Guipúzcoa 284275.4
## Shape_Area PROV COM
## 0 7980747663 15 12
## 1 1979346073 20 16class(prov.map)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"# nombre
prov.map$NOMBRE99## [1] Coruña (A) Guipúzcoa Cantabria
## [4] León Burgos Palencia
## [7] Vizcaya Pontevedra Palencia
## [10] Palencia Palencia Burgos
## [13] Zaragoza Burgos Pontevedra
## [16] Ourense Pontevedra Rioja (La)
## [19] Coruña (A) Girona Burgos
## [22] Burgos Girona Pontevedra
## [25] Pontevedra Navarra Navarra
## [28] Pontevedra Pontevedra Valladolid
## [31] Zamora Valladolid Girona
## [34] Barcelona Palencia Soria
## [37] Valladolid Palencia Lleida
## [40] Segovia Tarragona Salamanca
## [43] Ávila Madrid Castellón / Castelló
## [46] Madrid Guadalajara Cáceres
## [49] Toledo Balears (Illes) Valencia / València
## [52] Balears (Illes) Balears (Illes) Ciudad Real
## [55] Ciudad Real Badajoz Albacete
## [58] Balears (Illes) Balears (Illes) Balears (Illes)
## [61] Balears (Illes) Balears (Illes) Balears (Illes)
## [64] Murcia Córdoba Jaén
## [67] Huelva Sevilla Granada
## [70] Almería Murcia Murcia
## [73] Córdoba Málaga Palmas (Las)
## [76] Palmas (Las) Palmas (Las) Ceuta
## [79] Santa Cruz de Tenerife Palmas (Las) Santa Cruz de Tenerife
## [82] Melilla Santa Cruz de Tenerife Palmas (Las)
## [85] Santa Cruz de Tenerife Lugo Asturias
## [88] Cantabria Vizcaya Navarra
## [91] Álava Huesca Lleida
## [94] Barcelona Tarragona Teruel
## [97] Guadalajara Cuenca Valencia / València
## [100] Alicante /Alacant Cádiz
## 52 Levels: Álava Albacete Alicante /Alacant Almería Asturias ... Zaragoza# pero es mejor quedarse con variable del código
prov.map$PROV## [1] 15 20 39 24 09 34 48 36 34 34 34 09 50 09 36 32 36 26 15 17 09 09 17
## [24] 36 36 31 31 36 36 47 49 47 17 08 34 42 47 34 25 40 43 37 05 28 12 28
## [47] 19 10 45 07 46 07 07 13 13 06 02 07 07 07 07 07 07 30 14 23 21 41 18
## [70] 04 30 30 14 29 35 35 35 51 38 35 38 52 38 35 38 27 33 39 48 31 01 22
## [93] 25 08 43 44 19 16 46 03 11
## 52 Levels: 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 ... 52# pongo en minúsculas las variables en prov.map@data
colnames(prov.map@data) <- tolower(colnames(prov.map@data))Para pintar con ggplot es necesario convertir prov.map a un tipo de objeto adecuado.
prov.broom <- tidy(prov.map, region="prov")
head(prov.broom)## long lat order hole piece group id
## 1 537005 4770915 1 FALSE 1 01.1 01
## 2 537134 4770789 2 FALSE 1 01.1 01
## 3 536957 4770526 3 FALSE 1 01.1 01
## 4 537472 4769642 4 FALSE 1 01.1 01
## 5 537033 4769624 5 FALSE 1 01.1 01
## 6 536893 4769460 6 FALSE 1 01.1 01intervalos <- classIntervals(tasa.paro$tasa_estim_mixto, 6, style = "pretty")
intervalos## style: pretty
## [0,0.1) [0.1,0.2) [0.2,0.3) [0.3,0.4) [0.4,0.5) [0.5,0.6) [0.6,0.7]
## 94 152 115 64 26 15 2tasa.paro$intervalos <- as.factor(findInterval(tasa.paro$tasa_estim_mixto, intervalos$brks, all.inside=TRUE))Para poder unir los polígonos con el mapa necesitamos que el código de provincia se llame igual en ambos datasets.
head(tasa.paro)## # A tibble: 6 x 9
## provincia prov gedad nforma3 n n_parado tasa_estim_direc…
## <fct> <fct> <fct> <fct> <int> <int> <dbl>
## 1 Álava 01 De 16 a… Est primarios… 2 2 1
## 2 Álava 01 De 16 a… Est. Secundar… 65 12 0.186
## 3 Álava 01 De 16 a… Est. Universi… 50 4 0.121
## 4 Álava 01 De 35 a… Est primarios… 6 1 0.149
## 5 Álava 01 De 35 a… Est. Secundar… 175 27 0.156
## 6 Álava 01 De 35 a… Est. Universi… 232 11 0.0450
## # ... with 2 more variables: tasa_estim_mixto <dbl>, intervalos <fct>ggplot(tasa.paro) +
geom_map(aes(map_id = prov, fill = intervalos),
map = prov.broom,
colour = "grey30") +
expand_limits(x = prov.broom$long, y = prov.broom$lat) +
facet_grid(gedad ~ nforma3) +
scale_fill_brewer(
palette = "Reds",
labels = c(
"Menos del 10%",
"[10% - 20%)",
"[20% - 30%)",
"[30% - 40%)",
"[40% - 50%)",
"[50% - 60%)",
"[60% - 70%)",
"[70% - 80%]")) +
scale_x_continuous(breaks = NULL) +
scale_y_continuous(breaks = NULL) +
theme(
axis.text.y = element_blank(),
axis.text.x = element_blank(),
plot.title = element_text(face = "bold", size = rel(1.2)),
legend.text = element_text(
face = "bold",
color = "black",
size = rel(1.2)),
panel.background = element_rect(fill = "#EEEDF6"),
strip.text = element_text(face = "bold", size = rel(1))) +
labs(list(x = "", y = "", fill = "")) +
ggtitle("Tasa de paro 1T 2018\npor edad y estudios\n con modelo mixto")