Diapositivas del día 2

Se pueden descargar de este link (si no funciona, pégalo en el browser) https://www.dropbox.com/s/qpvh7hsf4ja8ssn/D%C3%ADa2.pptx?dl=0

Previo: preparación de nuestra base de datos

Cargamos la base de nuevo

*No hay necesidad de instalar los comandos de nuevo, pero es necesario para escribir esta libreta de instrucciones, per si checar nuestro directorio

install.packages("haven", repos = "http://cran.us.r-project.org", dependencies = TRUE)
## 
## The downloaded binary packages are in
##  /var/folders/fr/mw1x21js54367mjdhqsjfwqm0000gn/T//RtmpD3hUat/downloaded_packages
setwd("~/Dropbox/TeTra2018/CursoR")

Con la librería haven, abrimos la base

library(haven)
sdemt218 <- read_dta("SDEMT218.dta")

Volvemos a hacer nuestro filtro, hoy incluimos quitar los missings de la variable de ocupación, que será nuestra variable multinivel

filtro<-sdemt218$clase2==1 & sdemt218$anios_esc!=99 & sdemt218$ing_x_hrs>0
base<-sdemt218[filtro,]

Volvemos a hacer nuestra variable logaritmica

base$log_ing_x_hrs<-log(base$ing_x_hrs)

Modelos Multinivel

Repaso

Hacemos un modelo lineal, para repasar.

lineal0 <- lm(log_ing_x_hrs ~ anios_esc, data = base)
summary( lineal0 )
## 
## Call:
## lm(formula = log_ing_x_hrs ~ anios_esc, data = base)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.8067 -0.3737 -0.0191  0.3711  4.3131 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.7490158  0.0050084   548.9   <2e-16 ***
## anios_esc   0.0649206  0.0004609   140.8   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.663 on 114964 degrees of freedom
## Multiple R-squared:  0.1472, Adjusted R-squared:  0.1471 
## F-statistic: 1.984e+04 on 1 and 114964 DF,  p-value: < 2.2e-16

Vamos a guardar los coeficientes, porque los queremos graficar

install.packages("tidyverse", repos = "http://cran.us.r-project.org", dependencies = TRUE)
## 
## The downloaded binary packages are in
##  /var/folders/fr/mw1x21js54367mjdhqsjfwqm0000gn/T//RtmpD3hUat/downloaded_packages
library(broom)

lineal0
## 
## Call:
## lm(formula = log_ing_x_hrs ~ anios_esc, data = base)
## 
## Coefficients:
## (Intercept)    anios_esc  
##     2.74902      0.06492
coefplot <- tidy(lineal0)

Usando ggplot los vamos a graficar.

library(ggplot2)
q <-ggplot(base, aes(x =anios_esc, y = log_ing_x_hrs) ) + geom_point() + theme_minimal() +geom_abline( intercept = coefplot$estimate[1], slope = coefplot$estimate[2])
q

Modelo sin intercepto

En el cuestionario sociodemográfico hay una clasificación de ocupaciones (versión resumida del SINCO)

Categorías
  • 0 No aplica
  • 1 Profesionales, técnicos y trabajadores del arte
  • 2 Trabajadores de la educación
  • 3 Funcionarios y directivos
  • 4 Oficinistas
  • 5 Trabajadores industriales artesanos y ayudantes
  • 6 Comerciantes
  • 7 Operadores de transporte
  • 8 Trabajadores en servicios personales
  • 9 Trabajadores en protección y vigilancia
  • 10 Trabajadores agropecuarios
  • 11 No especificado

    • Vamos a ver cómo se comporta el logaritmo de los ingreso. Podemos forzar a que no haya intercepto, y tenemos básicamente una estimación de medias. (R2 es alto porque hay colinealidad)

    base$ocup<-as.factor(base$c_ocu11c)
lineal1 <- lm(log_ing_x_hrs ~ ocup - 1, data = base)
summary(lineal1)
    ## 
## Call:
## lm(formula = log_ing_x_hrs ~ ocup - 1, data = base)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.3540 -0.3595 -0.0193  0.3429  4.8084 
## 
## Coefficients:
##        Estimate Std. Error t value Pr(>|t|)    
## ocup1  3.952762   0.006378  619.80   <2e-16 ***
## ocup2  4.234475   0.010511  402.88   <2e-16 ***
## ocup3  4.167961   0.015953  261.27   <2e-16 ***
## ocup4  3.623236   0.006402  565.95   <2e-16 ***
## ocup5  3.363021   0.003514  956.92   <2e-16 ***
## ocup6  3.197970   0.004692  681.55   <2e-16 ***
## ocup7  3.371700   0.008282  407.11   <2e-16 ***
## ocup8  3.253216   0.004590  708.74   <2e-16 ***
## ocup9  3.486425   0.020988  166.12   <2e-16 ***
## ocup10 2.880596   0.006971  413.25   <2e-16 ***
## ocup11 3.440881   0.148798   23.12   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6486 on 114955 degrees of freedom
## Multiple R-squared:  0.9651, Adjusted R-squared:  0.9651 
## F-statistic: 2.893e+05 on 11 and 114955 DF,  p-value: < 2.2e-16

    Modelo multinivel

    1. Versión más simple, movemos el intercepto entre los grupos

    1 | group Intercepto aleatorio con media fija entre grupos

    install.packages("lme4", repos = "http://cran.us.r-project.org", dependencies = TRUE)
    ## 
## The downloaded binary packages are in
##  /var/folders/fr/mw1x21js54367mjdhqsjfwqm0000gn/T//RtmpD3hUat/downloaded_packages
    library(lme4)
    ## Loading required package: Matrix
    mlm <-lmer( log_ing_x_hrs  ~ 1 +  (1|ocup), data = base)

summary(mlm)
    ## Linear mixed model fit by REML ['lmerMod']
## Formula: log_ing_x_hrs ~ 1 + (1 | ocup)
##    Data: base
## 
## REML criterion at convergence: 226795.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -9.7969 -0.5544 -0.0298  0.5288  7.4135 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ocup     (Intercept) 0.1775   0.4213  
##  Residual             0.4207   0.6486  
## Number of obs: 114966, groups:  ocup, 11
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   3.5439     0.1277   27.75
    coef(mlm)
    ## $ocup
##    (Intercept)
## 1     3.952669
## 2     4.234046
## 3     4.167068
## 4     3.623218
## 5     3.363033
## 6     3.198013
## 7     3.371767
## 8     3.253250
## 9     3.486567
## 10    2.880777
## 11    3.452305
## 
## attr(,"class")
## [1] "coef.mer"

    No nos da los intervalos de confianza pero se los podemos pedir

    confint(mlm)
    ## Computing profile confidence intervals ...
    ##                 2.5 %    97.5 %
## .sig01      0.2781389 0.6530917
## .sigma      0.6459536 0.6512564
## (Intercept) 3.2827353 3.8046427

    También les podemos pedir que nos dé los estimados de los efectos aleatorios de esos intercpetos por grupo

    ranef(mlm)
    ## $ocup
##    (Intercept)
## 1   0.40878557
## 2   0.69016285
## 3   0.62318451
## 4   0.07933500
## 5  -0.18084970
## 6  -0.34586992
## 7  -0.17211618
## 8  -0.29063285
## 9  -0.05731558
## 10 -0.66310594
## 11 -0.09157782

    Y también podemos pedir los efectos fijos de la ecuación, que no cambian a través de los grupos

    fixef(mlm)
    ## (Intercept) 
##    3.543883
    1. Incluimos el predictor de los años de escolaridad

    1 | group Intercepto aleatorio con media fija entre grupos, pero con un predictor a nivel individual

    mlm1 <-lmer( log_ing_x_hrs  ~  anios_esc +  (1|ocup), data = base)
summary(mlm1)
    ## Linear mixed model fit by REML ['lmerMod']
## Formula: log_ing_x_hrs ~ anios_esc + (1 | ocup)
##    Data: base
## 
## REML criterion at convergence: 221740.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.1638  -0.5582  -0.0274   0.5424   7.2589 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ocup     (Intercept) 0.09722  0.3118  
##  Residual             0.40256  0.6345  
## Number of obs: 114966, groups:  ocup, 11
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) 3.1053285  0.0950266   32.68
## anios_esc   0.0387577  0.0005383   71.99
## 
## Correlation of Fixed Effects:
##           (Intr)
## anios_esc -0.064

    Intervalos y resultados por grupos

    confint(mlm1)
    ## Computing profile confidence intervals ...
    ##                  2.5 %     97.5 %
## .sig01      0.20544828 0.48392293
## .sigma      0.63188705 0.63707440
## (Intercept) 2.91091789 3.29905056
## anios_esc   0.03770488 0.03981566
    ranef(mlm1)
    ## $ocup
##    (Intercept)
## 1   0.28091803
## 2   0.50686391
## 3   0.49234288
## 4   0.01955873
## 5  -0.08390156
## 6  -0.28588355
## 7  -0.09302706
## 8  -0.18165758
## 9  -0.06322703
## 10 -0.47892704
## 11 -0.11305975
    1. Hoy también dejaremos que las pendientes de los años de escolaridad cambien entre las ocupaciones

    Para efectos ilustrativos correremos este código para ver cómo se estiman las pendientes en 11 regresiones separadas en cada grupo

    beta.hat <- list()
alpha.hat <- list()
for(i in 1:11){
  unit.lm <- lm(log_ing_x_hrs  ~  anios_esc , data = base[base$c_ocu11c == i,] )
  beta.hat[i] <- coef(unit.lm)[2]
  alpha.hat[i] <- coef(unit.lm)[1]
}
beta.hat <- as.numeric(beta.hat)
alpha.hat <- as.numeric(alpha.hat)

    Este código nos guardó las 11 pendientes. Hoy las vamos a graficar, los coeficientes de los años de escolaridad

    hist(beta.hat)

    También podemos graficar los 11 interceptos de los modelos individuales

    hist(alpha.hat)

    ¿Cómo introducir esto a nuestro modelo?

    mlm2 <-lmer( log_ing_x_hrs  ~  anios_esc +  (1 + anios_esc|ocup), data = base)
summary(mlm2)
    ## Linear mixed model fit by REML ['lmerMod']
## Formula: log_ing_x_hrs ~ anios_esc + (1 + anios_esc | ocup)
##    Data: base
## 
## REML criterion at convergence: 220986.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.1895  -0.5485  -0.0260   0.5428   7.1911 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  ocup     (Intercept) 0.0458785 0.21419       
##           anios_esc   0.0003964 0.01991  -0.22
##  Residual             0.3998152 0.63231       
## Number of obs: 114966, groups:  ocup, 11
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 2.981411   0.068159  43.742
## anios_esc   0.046123   0.006253   7.376
## 
## Correlation of Fixed Effects:
##           (Intr)
## anios_esc -0.279

    Intervalos y resultados por grupos

    #confint(mlm2)
ranef(mlm2)
    ## $ocup
##    (Intercept)    anios_esc
## 1   0.15643142  0.009628791
## 2   0.19193500  0.019989495
## 3   0.12763093  0.025926304
## 4  -0.18104010  0.017883219
## 5   0.09677424 -0.013808539
## 6  -0.27105338  0.003808557
## 7   0.28281720 -0.034567840
## 8   0.09612104 -0.025465585
## 9  -0.06819787  0.003895420
## 10 -0.33898153 -0.009783750
## 11 -0.09243697  0.002493929
    fixef(mlm2)
    ## (Intercept)   anios_esc 
##  2.98141112  0.04612347
    1. Variables de nivel superior

    Supongamos que el nivel de participación de las mujeres en las ocupaciones hace una diferencia en sus ingresos por hora.

    Chequemos si esto varía con un tabulado de sexo, guardándolo en un objeto

    freq.sexo<-table(base$c_ocu11c,base$sex)
freq.sexo
    ##     
##          1     2
##   1   5940  4403
##   2   1452  2356
##   3   1057   596
##   4   4423  5841
##   5  25269  8791
##   6   8897 10210
##   7   6065    68
##   8   7616 12350
##   9    852   103
##   10  7905   753
##   11    11     8

    Podemos usar la función prop.table para ver los resultados, usando la opción 1, porque queremos la proporción por filas

    prop.sexo<-prop.table(freq.sexo,1)
prop.sexo
    ##     
##               1          2
##   1  0.57430146 0.42569854
##   2  0.38130252 0.61869748
##   3  0.63944344 0.36055656
##   4  0.43092362 0.56907638
##   5  0.74189665 0.25810335
##   6  0.46564086 0.53435914
##   7  0.98891244 0.01108756
##   8  0.38144846 0.61855154
##   9  0.89214660 0.10785340
##   10 0.91302841 0.08697159
##   11 0.57894737 0.42105263

    Como es un objeto podemos convertirlo en un data frame y hacerle un merge

    prop.sexo<-as.data.frame(prop.sexo)
#Me quedo con los datos de las mujeres
prop.sexo<-prop.sexo[prop.sexo$Var2==2,]

    Hoy creamos la variable c_ocu11c, y ponemos nombre de las variables de nuestra mini-dataframe

    names(prop.sexo)<-c( "c_ocu11c", "sex_ind", "part_mujer")
    base2<-merge(base, prop.sexo, by="c_ocu11c")

    Hoy sí nuestro modelo

    mlm3 <-lmer( log_ing_x_hrs  ~  anios_esc + part_mujer + (1 + anios_esc|ocup), data = base2)
summary(mlm3)
    ## Linear mixed model fit by REML ['lmerMod']
## Formula: log_ing_x_hrs ~ anios_esc + part_mujer + (1 + anios_esc | ocup)
##    Data: base2
## 
## REML criterion at convergence: 220986.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.1895  -0.5485  -0.0260   0.5428   7.1911 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  ocup     (Intercept) 0.0519675 0.22796       
##           anios_esc   0.0003972 0.01993  -0.26
##  Residual             0.3998152 0.63231       
## Number of obs: 114966, groups:  ocup, 11
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  2.94111    0.13641  21.561
## anios_esc    0.04616    0.00627   7.362
## part_mujer   0.10947    0.31981   0.342
## 
## Correlation of Fixed Effects:
##            (Intr) ans_sc
## anios_esc  -0.161       
## part_mujer -0.848 -0.006

    Intervalos y resultados por grupos

    #confint(mlm3)
ranef(mlm3)
    ## $ocup
##    (Intercept)    anios_esc
## 1   0.15052437  0.009564155
## 2   0.16893956  0.019679706
## 3   0.12951665  0.025817267
## 4  -0.20316515  0.017854262
## 5   0.10883134 -0.013848593
## 6  -0.28932548  0.003777061
## 7   0.32210137 -0.034627727
## 8   0.06880723 -0.025513599
## 9  -0.04361194  0.004171383
## 10 -0.30846593 -0.009792004
## 11 -0.10415201  0.002918088
    fixef(mlm3)
    ## (Intercept)   anios_esc  part_mujer 
##  2.94111148  0.04616239  0.10946496

    Modelos con tres niveles

    Supongamos que los individuos están anidados en ocupaciones y que esas ocupaciones están anidadas en entidades.

    mlm4 <-lmer( log_ing_x_hrs  ~  anios_esc  +  (1|ent) + (1|ocup), data = base2)
summary(mlm4)
    ## Linear mixed model fit by REML ['lmerMod']
## Formula: log_ing_x_hrs ~ anios_esc + (1 | ent) + (1 | ocup)
##    Data: base2
## 
## REML criterion at convergence: 215749.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.2814  -0.5606  -0.0354   0.5371   7.6091 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ent      (Intercept) 0.02174  0.1474  
##  ocup     (Intercept) 0.09726  0.3119  
##  Residual             0.38157  0.6177  
## Number of obs: 114966, groups:  ent, 32; ocup, 11
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) 3.1216078  0.0985136   31.69
## anios_esc   0.0371671  0.0005268   70.56
## 
## Correlation of Fixed Effects:
##           (Intr)
## anios_esc -0.060

    Intervalos y resultados por grupos

    #confint(mlm4)
ranef(mlm4)
    ## $ent
##    (Intercept)
## 1  -0.05039768
## 2   0.19524380
## 3   0.28600172
## 4  -0.14241680
## 5   0.10920082
## 6   0.08754683
## 7  -0.29803975
## 8   0.13511651
## 9   0.08265613
## 10 -0.05894533
## 11  0.02440909
## 12 -0.22627447
## 13 -0.05361581
## 14  0.14481986
## 15 -0.08696546
## 16  0.07891327
## 17 -0.13480445
## 18  0.07047623
## 19  0.22404620
## 20 -0.23145790
## 21 -0.14317160
## 22  0.10680685
## 23  0.06103538
## 24  0.02179712
## 25  0.16175025
## 26  0.16464432
## 27 -0.03722892
## 28  0.02900936
## 29 -0.22042268
## 30 -0.12264191
## 31 -0.13067663
## 32 -0.04641436
## 
## $ocup
##    (Intercept)
## 1   0.27783311
## 2   0.53040771
## 3   0.49018701
## 4   0.01212813
## 5  -0.09236378
## 6  -0.27852936
## 7  -0.09745646
## 8  -0.18794053
## 9  -0.05589967
## 10 -0.44795059
## 11 -0.15041561
    fixef(mlm4)
    ## (Intercept)   anios_esc 
##  3.12160775  0.03716709