MicroDatosEs, EPA y modelos mixtos
José Luis Cañadas Reche
Agosto de 2018
Librerías
list.of.packages <- c("MicroDatosEs", "lme4","RCurl", "dplyr", "sjPlot", "ggplot2", "classInt", "reshape2", "sf")
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)##
## 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## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source'))`## Linking to GEOS 3.6.2, GDAL 2.2.3, proj.4 4.9.3## [[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] TRUEDescargamos la EPA
¡ El INE tiene un ftp accesible!
file_location <- "ftp://www.ine.es/temas/epa/datos_2t18.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
Seleccionamos solo las variables que nos interesan
# 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
}## # A tibble: 6 x 7
## prov edad nforma aoi factorel nforma3 gedad
## <chr> <chr> <chr> <fct> <dbl> <fct> <fct>
## 1 Álava de 45 A 49… Educación s… o 22506 Est. Universi… De 35 a 54
## 2 Álava de 45 A 49… Educación s… o 22506 Est. Universi… De 35 a 54
## 3 Álava de 35 A 39… Educación s… o 23231 Est. Universi… De 35 a 54
## 4 Álava de 35 A 39… Educación p… i 23231 Est primarios… De 35 a 54
## 5 Álava de 5 A 9 a… <NA> <NA> 23414 <NA> "15 años o…
## 6 Álava de 5 A 9 a… <NA> <NA> 23414 <NA> "15 años o…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: 161,506 x 8
## prov edad nforma aoi factorel nforma3 gedad llave
## <chr> <chr> <chr> <fct> <dbl> <fct> <fct> <chr>
## 1 Álava de 45 A … Educación … o 22506 Est. Univer… De 35 a … 1
## 2 Álava de 45 A … Educación … o 22506 Est. Univer… De 35 a … 1
## 3 Álava de 35 A … Educación … o 23231 Est. Univer… De 35 a … 1
## 4 Álava de 35 A … Educación … i 23231 Est primari… De 35 a … 1
## 5 Álava de 5 A 9… <NA> <NA> 23414 <NA> "15 años… 1
## 6 Álava de 5 A 9… <NA> <NA> 23414 <NA> "15 años… 1
## 7 Álava de 45 A … Educación … o 20807 Est. Univer… De 35 a … 1
## 8 Álava de 40 A … Educación … o 20807 Est. Univer… De 35 a … 1
## 9 Álava de 0 A 4… <NA> <NA> 20808 <NA> "15 años… 1
## 10 Álava de 40 A … Educación … o 20146 Est. Univer… De 35 a … 1
## # ... with 161,496 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.
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
## [[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
Para leer el mapa utilizo la función st_read del paquete sf que es un paquete nacido bajo el paraguas del R-Consortium. Este paquete permite leer shapefiles y guardar su geometría en un objeto de tipo data.frame, lo cual facilita mucho su tratamiento.
## Reading layer `spain_provinces_ind_2' from data source `/home/jose/Dropbox/Public/shapefile_spain/spain_provinces_ind_2.shp' using driver `ESRI Shapefile'
## Simple feature collection with 116 features and 10 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -389564.9 ymin: 3826792 xmax: 1127057 ymax: 4859444
## epsg (SRID): NA
## proj4string: NA## [1] "sf" "data.frame"## Simple feature collection with 6 features and 10 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -13952 ymin: 4589271 xmax: 603260 ymax: 4859444
## epsg (SRID): NA
## proj4string: NA
## OBJECTID AREA PERIMETER P20099_ P20099_ID NOMBRE99
## 1 1 7.980748e+09 1032580.0333 2 1 Coruña (A)
## 2 2 1.979346e+09 284275.3736 3 12 Guipúzcoa
## 3 3 1.241556e+03 425.0905 4 0 <NA>
## 4 4 1.963521e+07 21810.3954 5 15 Cantabria
## 5 5 1.559045e+10 821443.1246 6 16 León
## 6 6 1.399732e+10 994635.4371 7 19 Burgos
## Shape_Leng Shape_Area PROV COM geometry
## 1 1032579.9729 7.980748e+09 15 12 POLYGON ((91555 4756380, 90...
## 2 284275.4018 1.979346e+09 20 16 POLYGON ((603213 4794695, 6...
## 3 425.0838 1.241235e+03 <NA> <NA> POLYGON ((568823.4 4793617,...
## 4 21810.4169 1.963520e+07 39 06 POLYGON ((479862 4789888, 4...
## 5 821443.0056 1.559045e+10 24 07 POLYGON ((187928 4685080, 1...
## 6 994635.4222 1.399732e+10 09 07 POLYGON ((497093 4759050, 4...## [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 116## [1] "sf" "data.frame"## [1] "Coruña (A)" "Guipúzcoa"
## [3] "Cantabria" "León"
## [5] "Burgos" "Palencia"
## [7] "Vizcaya" "Pontevedra"
## [9] "Palencia" "Palencia"
## [11] "Palencia" "Burgos"
## [13] "Zaragoza" "Burgos"
## [15] "Pontevedra" "Ourense"
## [17] "Pontevedra" "Rioja (La)"
## [19] "Coruña (A)" "Girona"
## [21] "Burgos" "Burgos"
## [23] "Girona" "Pontevedra"
## [25] "Pontevedra" "Navarra"
## [27] "Navarra" "Pontevedra"
## [29] "Pontevedra" "Valladolid"
## [31] "Zamora" "Valladolid"
## [33] "Girona" "Barcelona"
## [35] "Palencia" "Soria"
## [37] "Valladolid" "Palencia"
## [39] "Lleida" "Segovia"
## [41] "Tarragona" "Salamanca"
## [43] "Ávila" "Madrid"
## [45] "Castellón / Castelló" "Madrid"
## [47] "Guadalajara" "Cáceres"
## [49] "Toledo" "Balears (Illes)"
## [51] "Valencia / València" "Balears (Illes)"
## [53] "Balears (Illes)" "Ciudad Real"
## [55] "Ciudad Real" "Badajoz"
## [57] "Albacete" "Balears (Illes)"
## [59] "Balears (Illes)" "Balears (Illes)"
## [61] "Balears (Illes)" "Balears (Illes)"
## [63] "Balears (Illes)" "Murcia"
## [65] "Córdoba" "Jaén"
## [67] "Huelva" "Sevilla"
## [69] "Granada" "Almería"
## [71] "Murcia" "Murcia"
## [73] "Córdoba" "Málaga"
## [75] "Palmas (Las)" "Palmas (Las)"
## [77] "Palmas (Las)" "Ceuta"
## [79] "Santa Cruz de Tenerife" "Palmas (Las)"
## [81] "Santa Cruz de Tenerife" "Melilla"
## [83] "Santa Cruz de Tenerife" "Palmas (Las)"
## [85] "Santa Cruz de Tenerife" "Lugo"
## [87] "Asturias" "Cantabria"
## [89] "Vizcaya" "Navarra"
## [91] "Álava" "Huesca"
## [93] "Lleida" "Barcelona"
## [95] "Tarragona" "Teruel"
## [97] "Guadalajara" "Cuenca"
## [99] "Valencia / València" "Alicante /Alacant"
## [101] "Cádiz"## [1] "15" "20" "39" "24" "09" "34" "48" "36" "34" "34" "34" "09" "50" "09"
## [15] "36" "32" "36" "26" "15" "17" "09" "09" "17" "36" "36" "31" "31" "36"
## [29] "36" "47" "49" "47" "17" "08" "34" "42" "47" "34" "25" "40" "43" "37"
## [43] "05" "28" "12" "28" "19" "10" "45" "07" "46" "07" "07" "13" "13" "06"
## [57] "02" "07" "07" "07" "07" "07" "07" "30" "14" "23" "21" "41" "18" "04"
## [71] "30" "30" "14" "29" "35" "35" "35" "51" "38" "35" "38" "52" "38" "35"
## [85] "38" "27" "33" "39" "48" "31" "01" "22" "25" "08" "43" "44" "19" "16"
## [99] "46" "03" "11"# pongo en minúsculas las variables en prov.map@data
colnames(prov.map) <- tolower(colnames(prov.map))Hago 6 intervalos “bonitos” de la tasa de paro estimada por el modelo mixto utilizando la función classIntervals y con la función findInterval se lo asigno a cada combinación de provincia, edad y nivel de estudios del data frame
## 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]
## 119 155 103 59 22 10Para poder unir los polígonos con el mapa necesitamos que el código de provincia se llame igual en ambos datasets.
## # A tibble: 6 x 9
## provincia prov gedad nforma3 n n_parado tasa_estim_dire…
## <fct> <fct> <fct> <fct> <int> <int> <dbl>
## 1 Álava 01 De 1… Est pr… 1 0 0
## 2 Álava 01 De 1… Est. S… 68 13 0.180
## 3 Álava 01 De 1… Est. U… 51 5 0.115
## 4 Álava 01 De 3… Est pr… 7 1 0.122
## 5 Álava 01 De 3… Est. S… 167 20 0.125
## 6 Álava 01 De 3… Est. U… 231 14 0.0621
## # ... with 2 more variables: tasa_estim_mixto <dbl>, intervalos <fct>Mapa con ggplot2
## Warning: Column `prov` joining factor and character vector, coercing into
## character vector# Utilizamos geom_sf para la geometría
ggplot(datos_pint) +
geom_sf(aes( fill = intervalos),
colour = "grey30") +
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 2T 2018\npor edad y estudios\n con modelo mixto")
Y voilá, como siempre tenemos que a mayor edad menor tasa de paro, a mayor nivel de estudios menor tasa de paro y qué si vives por el norte mucho mejor para tus perspectivas laborales.