Nest data en R

Pese al título, esta entrada solo pretende intentar explicar las funcionalidades de una data tipo nest. Las próximas líneas presentan un marco de datos, que es una forma tradicional de una base de datos (no requiere explicación) y un nest data, en esta última agrupamos por variable/clase haciendo un collapse o summarise de la data donde variables/clases pasan a ser el id de la data. Posteriormente, las próximas celdas funcionan como lista, por ende, en ella podemos guardar diversos tipos de objetos, sin importar su dimensión y recuperar estos objetos adelante, sin necesidad de crear nuevos objetos y manteniendo una mejor organización.  

# Marco de datos

+ head(iris)

  Sepal.Length Sepal.Width Petal.Length Petal.Width Species

1          5.1         3.5          1.4         0.2  setosa

2          4.9         3.0          1.4         0.2  setosa

3          4.7         3.2          1.3         0.2  setosa

4          4.6         3.1          1.5         0.2  setosa

5          5.0         3.6          1.4         0.2  setosa

6          5.4         3.9          1.7         0.4  setosa

# nest data

> iris %>%

+   group_by(Species) %>%

+   nest()

# A tibble: 3 x 2

  Species    data            

  <fct>      <list>          

1 setosa     <tibble [50 x 4]>

2 versicolor <tibble [50 x 4]>

3 virginica  <tibble [50 x 4]>

Ahora se utiliza la función vectorizada (lapply) para estimar un modelo de regresión que se aplicara para cada una de los grupos de especies disponibles den la data:

nestdata <- iris %>%

  group_by(Species) %>%

  nest() %>%

  mutate(suma = lapply(data, function(x) sum(x[[1]])),

 MCO_model = lapply(data, function(x) lm(Sepal.Length ~ Petal.Length, data = x)))

nestdata

nestdata

# A tibble: 3 x 4

  Species    data              suma      MCO_model

  <fct>      <list>            <list>    <list>  

1 setosa     <tibble [50 x 4]> <dbl [1]> <lm>    

2 versicolor <tibble [50 x 4]> <dbl [1]> <lm>    

3 virginica  <tibble [50 x 4]> <dbl [1]> <lm>

# equivalente

lm(Sepal.Length ~ Petal.Length, data = iris[Specie== “setosa”    ,])

lm(Sepal.Length ~ Petal.Length, data = iris[Specie== “versicolor”,])

lm(Sepal.Length ~ Petal.Length, data = iris[Specie== “virginica” ,])

Al imprimir la estructura de la data nest se puede notar que en la primera columna se guarda en cada celda una lista con la información de la data correspondiente a cada variable y en la columna MCO_model se guardan los resultados de los modelos marcados como equivalentes:

nestdata$data[[1]]

# A tibble: 50 x 4

   Sepal.Length Sepal.Width Petal.Length Petal.Width

          <dbl>       <dbl>        <dbl>       <dbl>

 1          5.1         3.5          1.4         0.2

 2          4.9         3            1.4         0.2

 3          4.7         3.2          1.3         0.2

 4          4.6         3.1          1.5         0.2

 5          5           3.6          1.4         0.2

 6          5.4         3.9          1.7         0.4

 7          4.6         3.4          1.4         0.3

 8          5           3.4          1.5         0.2

 9          4.4         2.9          1.4         0.2

10          4.9         3.1          1.5         0.1

# ... with 40 more rows

nestdata$MCO_model[[1]]

Call:

lm(formula = Sepal.Length ~ Petal.Length, data = x)

Coefficients:

 (Intercept)  Petal.Length 

      4.2132        0.5423

También las funciones pueden modificarse para agregar otros elementos:

dataness <- iris %>%

  group_by(Species) %>%

  nest() %>%

  mutate(MCO_model = lapply(data, function(x) lm(Sepal.Length ~ Petal.Length, data = x)),

         rediduo =  lapply(data, function(x) lm(Sepal.Length ~ Petal.Length, data = x)$residuals))

dataness

# A tibble: 3 x 4

  Species    data              MCO_model rediduo  

  <fct>      <list>            <list>    <list>   

1 setosa     <tibble [50 x 4]> <lm>      <dbl [50]>

2 versicolor <tibble [50 x 4]> <lm>      <dbl [50]>

3 virginica  <tibble [50 x 4]> <lm>      <dbl [50]>

Note que al hacer un unnest de la data original, podemos recuperar la data original, pero con el residuo correspondiente a cada regresión anexado:

resids <- unnest(dataness, rediduo, data)

resids

# A tibble: 150 x 6

   Species rediduo Sepal.Length Sepal.Width Petal.Length Petal.Width

   <fct>     <dbl>        <dbl>       <dbl>        <dbl>       <dbl>

 1 setosa   0.128           5.1         3.5          1.4         0.2

 2 setosa  -0.0724          4.9         3            1.4         0.2

 3 setosa  -0.218           4.7         3.2          1.3         0.2

 4 setosa  -0.427           4.6         3.1          1.5         0.2

 5 setosa   0.0276          5           3.6          1.4         0.2

 6 setosa   0.265           5.4         3.9          1.7         0.4

 7 setosa  -0.372           4.6         3.4          1.4         0.3

 8 setosa  -0.0266          5           3.4          1.5         0.2

 9 setosa  -0.572           4.4         2.9          1.4         0.2

10 setosa  -0.127           4.9         3.1          1.5         0.1

# ... with 140 more rows

resids %>%

  ggplot(aes(x=Petal.Length, y=rediduo, group = Species)) +

  geom_point() +

  geom_hline(yintercept=0,col="red") +

  facet_wrap(~Species, nrow = 1, scales = "free_x")

Referencias

-          An, S. (2017). Linear regression and residuals. In: http://rstudio-pubs-static.s3.amazonaws.com/310380_07eebbdd4fac48aba8399c738ae0c31d.html.

-          Barrett, Tyson. List-columns in data.table: Nesting and unnesting data 2 tables and vectors. Utah State University.

-          Dancho, Matt (2018). Forecasting Time Series Groups in the tidyverse. In: https://cran.rstudio.com/web/packages/sweep/vignettes/SW01_Forecasting_Time_Series_Groups.html.

Generando matrices markovianas de transición en el mercado laboral dominicano con la Encuesta Continua de Fuerza de Trabajo (ENCFT)

En el próximo código se crean las probabilidades de transición de la PET entre tres estados: ocupados, desocupados e inactivos. Note que cargamos todas (18 bases) las ENCFT entre 2014-3 y 2018-4, lo que permite considerar todas las bases dentro de un bucle sin necesidad de repetir el procedimiento base por base. Todas las bases se guardan en un objeto único llamado miembros.

rm(list = ls(all.names = TRUE))

library(dplyr)

# Dinámina del desempleo

setwd("C:/Users/Nery Ramirez/Dropbox/Dinamica ML/Data")

# Carga todas las bases

miembros <- lapply(Sys.glob("miembros*.txt"), read.delim)

bases<-sum(lengths(lengths(miembros)))

[1] 18

Podemos verificar la lista miembro, donde hemos cargado todas las bases de datos. Ahora procedemos a crear dos matrices: una para guardar los nombres de las bases con la que trabajamos, lo que nos permitirá leer mejor los resultados (matrixNombres [18x1] 18 es el número de bases cargadas) y otra donde se guardan los coeficientes (resultFinalDesempleo [36x1]).

La idea general, como estamos trabajando con 3 estados, es tener 4 matrices 3x3, como explicamos adelante, convertirlas en vector y terminar con una matriz 36x18. Ahora creamos un for anidado, donde la general es para cada par de base levantadas en periodos continuos, se crean 4 matrices 3x3 para capturar las probabilidades de transición de un estado a otro.

La primera parte del bucle anidado [for (jbase in 1:(bases-1)){] permite cargar las bases de forma continuas (ejemplo: 2014-1 y 2014-2), note que va desde la primera base hasta la penúltima, porque dentro se obtiene la segunda base como el índice del bucle más 1 [[[jbase+1]]], a estas bases la llamamos base 1 y base 2.

Posteriormente, combinando los nombres de casa encuesta (data1 y data2) se crea un vector [1x1] conteniendo la combinación de nombre de las encuestas correspondientes ("20143-20144").

(Perdonando por no usar mutate de tidyverse, dicho procedimiento se hace de la manera tradicional utilizando el operador $) ahora se crea una variable cualitativa con la información de los tres estados. Note que solo debemos indicarle a la variable dicotómica de ocupados (1=ocupados) cuales son los inactivos del mercado laboral:

0 = empleados

1=desempleado

2=inactivos

table(data1$stado1)

    0     1     2

 8963   509 10755

Aunque solo se utiliza una vez, posteriormente de deja en comentario una parte donde se pretende validar visualmente si no se había cometido algún error con los códigos, al programar tareas este paso es sumamente importante, porque evita cometer errores graves. Aunque lo mejor siempre es tratar de establecer algún procedimiento de validación.

Luego creamos la base3 combinando base1 y base2 por el [by="ID_PERSONA"], posteriormente nos quedamos solo con las personas con edad de trabajar en la primera encuesta [dplyr::filter(PET.x==1)].

Luego se crean 4 matrices [3x3] que se resetean con cada vuelta del bucle:

MatrixMarkPro: guarda probabilidades de transición sin considerar el factor de expansión.

MatrixMarkPer: guarda el número de personas según la condición de transferencia entre un estado y otro, sin considerar el factor de expansión.  

MatrixMarkProW: guarda probabilidades de transición considerando el factor de expansión. 

MatrixMarkPerW: guarda el número de personas según la condición de transferencia entre un estado y otro, considerando el factor de expansión.  

Finalmente usamos un for doble o anidado, dentro del for anterior, para completar estas matrices, note que la idea es muy sencilla, recorriendo el número de estados en cada bucle [0:sk] para crear probabilidades para construir la primera matriz, en la posición [1,1] estamos verificando la probabilidad de que un individuo este en el estado cero y continúe en ese estado, porque estos son los dos primeros valores que asumen los índices del bucle.

for (s1 in 0:sk) {

  for (s2 in 0:sk) {

    # sin factor de expansión

    MatrixMarkPro[s1+1,s2+1]  <- sum(stado1==s1 & stado2==s2)/sum(stado1==s1)

    MatrixMarkPer[s1+1,s2+1]  <- sum(stado1==s1 & stado2==s2)

   

    # Usando factor de expansión

    MatrixMarkProW[s1+1,s2+1] <- sum(FACTOR_EXPANSION.x[stado1==s1 & stado2==s2])/sum(FACTOR_EXPANSION.x[stado1==s1])

    MatrixMarkPerW[s1+1,s2+1] <- sum(FACTOR_EXPANSION.x[stado1==s1 & stado2==s2])

  }

}

Asignémosle valores a los índices del bucle en la primera y segunda corrida del bucle, para plantearlo más claro:

Corrida 1 del blucle S1 = 0 S2 = 0 Por tanto, estamos verificando la probabilidad [en la matriz MatrixMarkPro] de permanecer en el estado cero (de permanecer empleados). Para ello sumamos la cantidad de individuos que tienen este estado en ambas bases [sum(stado1==0 & stado2==0)]  y se divide entre el total de PET con el estado 0 en la primera encuesta [sum(stado1==0)]. En el caso de contar personas, solo eliminamos la división entre sum(stado1==0) de la cantidad de personas que pasa de un estado a otro. Por último, para considerar el factor de expansión repetimos el ejercicio anterior, pero sumando en vez de estados (un 1 por cada observación) se suma el factor de expansión correspondiente.

Corrida 2 del blucle S1 = 0 S2 = 1 Por tanto, estamos verificando la probabilidad [en la matriz MatrixMarkPro] de pasar del estado 0 al 1, empleado a desempleado. Para ello sumamos la cantidad de individuos que tienen el estado cero en la base 0 y 1 en la base 1 [sum(stado1==0 & stado2==1)]  y se divide entre el total de PET con el estado 0 en la primera encuesta [sum(stado1==0)]. En el caso de contar personas, solo eliminamos la división entre sum(stado1==0) de la cantidad de personas que pasa de un estado a otro. Por último, para considerar el factor de expansión repetimos el ejercicio anterior, pero sumando en vez de estados (un 1 por cada observación) se suma el factor de expansión correspondiente.

Las matrices creadas en el paso anterior se transforman en vectores de forma que podemos colocar una matriz debajo de otra, correspondiendo cada observación a los coeficientes estimados para base de datos, con el fin de poder organizar todas las matrices en una tabla plana:

> MatrixMarkPro

           [,1]       [,2]       [,3]

[1,] 0.95254403 0.01190476 0.03555121

[2,] 0.30352304 0.42005420 0.27642276

[3,] 0.07665222 0.01814735 0.90520043

> as.vector(MatrixMarkPro)

[1] 0.95254403 0.30352304 0.07665222 0.01190476 ...

L hubiese sido mucho más eficiente guardar estas matrices en una lista y luego convertirlas en una tabla plana, es bien fácil usando la función apply, hace el código más corto y entendible y mantiene disponible toda la información generada. Esto lo vi luego de crear el programa y me dio flojera corregirlo.

# ::::::::::::::::::::::::::::::::::::::::::

# NOTA: correr el bucle anidado desde aquí

# ::::::::::::::::::::::::::::::::::::::::::

matrixNombres <-c()

skk <- 3 #considerados estados

resultFinalDesempleo <- matrix(NA, nrow=1, ncol=skk*skk*4) 

for (jbase in 1:(bases-1)){

data1 <- as.data.frame(miembros[[jbase]])

data2 <- as.data.frame(miembros[[jbase+1]])

namess<-paste("d",data1$TRIMESTRE[1],data2$TRIMESTRE[2],sep = "_")

# crea tres estados 0=empleado; 1=empleado 3=inactivo

data1$stado1 <- data1$DESOCUPADO

data1$stado1[data1$PEA==0]<-2

data2$stado2 <- data2$DESOCUPADO

data2$stado2[data2$PEA==0]<-2

# Validar

# View(cbind(esta=data1$stado1, desocupa=data1$DESOCUPADO, pea=data1$PEA, data2$stado2, data2$DESOCUPADO, data2$PEA))

# https://rpubs.com/williamsurles/293454

data3 <- inner_join(data1,data2,by="ID_PERSONA") %>%

         dplyr::filter(PET.x==1)

attach(data3)

# :::::::::::::: Automatico

sk<-skk -1   #length(table(t1))  #Estados: Emmpelo, desempleo

MatrixMarkPro <- matrix(NA, sk+1, sk+1) # Probabilidades

MatrixMarkPer <- matrix(NA, sk+1, sk+1) # Cuenta personas

MatrixMarkProW <- matrix(NA, sk+1, sk+1) # Probabilidades

MatrixMarkPerW <- matrix(NA, sk+1, sk+1) # Cuenta personas

for (s1 in 0:sk) {

  for (s2 in 0:sk) {

    MatrixMarkPro[s1+1,s2+1]  <- sum(stado1==s1 & stado2==s2)/sum(stado1==s1)

    MatrixMarkPer[s1+1,s2+1]  <- sum(stado1==s1 & stado2==s2)

    

    # Usando factor de expansión

    MatrixMarkProW[s1+1,s2+1] <- sum(FACTOR_EXPANSION.x[stado1==s1 & stado2==s2])/sum(FACTOR_EXPANSION.x[stado1==s1])

    MatrixMarkPerW[s1+1,s2+1] <- sum(FACTOR_EXPANSION.x[stado1==s1 & stado2==s2])

  }

} #Cierra creación de matrices para el periodo

dinamicaMatrix <- c(as.vector(MatrixMarkPro),

                              as.vector(MatrixMarkPer),

                              as.vector(MatrixMarkProW),

                              as.vector(MatrixMarkPerW))

resultFinalDesempleo <- rbind(resultFinalDesempleo,dinamicaMatrix)

matrixNombres<-c(matrixNombres,paste(data1$TRIMESTRE[1],data2$TRIMESTRE[1],sep = "-"))

print(paste(data1$TRIMESTRE[1],data2$TRIMESTRE[1],sep = "-"))

 }

#Cierra el bucle de cargar base de datos

Resultados

En el caso de las probabilidades de transición obtenidas para el periodo [1] "20183-20184" y colocadas debajo, se verifica que en el periodo un 95% de las personas que estaban ocupadas permaneció ocupada pera el periodo siguiente, el 1.15% paso a ser desempleado y 3.3% paso a estar inactivo. El 45% de los desempleados continúo desempleado, el 29% paso a estar inactiva y el 26% encontró empleo. El 90% de los inactivos continuaron en la misma posición y 7.4% encontró trabajo, 2% comenzó a buscar empleo.

> MatrixMarkPro            [,1]       [,2]       [,3] [1,] 0.95254403 0.01190476 0.03555121 [2,] 0.30352304 0.42005420 0.27642276 [3,] 0.07665222 0.01814735 0.90520043

> MatrixMarkProW            [,1]       [,2]       [,3] [1,] 0.95516926 0.01152740 0.03330334 [2,] 0.26040401 0.45056803 0.28902796 [3,] 0.07445017 0.01930278 0.90624705

Graficando las columnas de esta matriz de resultados fines, podemos estudiar la evolución histórica de las matrices de transición (L lo siento por no terminar el proceso en R). Se observa como durante el perioro la probabilidad de permanecer empleado se ha incrementado, mientras se reduce la probabiildad de perder el empleo o quedar desempleado.