Práctica 8. Modelos generalizados: La regresión logística
AE
20/01/2020
Previo
¡Recuerda poner el directorio!
setwd("/Users/anaescoto/Dropbox/2020/R_invierno")Volvemos a importar la base completa del cuestionario sociodemográfico de la Encuesta Nacional de Ocupación y Empleo, trimestre III de 2019. O si guardaste tu ambiente, vuévelo a cargar.
library(haven)
SDEMT319 <- read_dta("./datos/SDEMT319.dta")Nuevos paquetes
#install.packages("ResourceSelection", dependencies = TRUE)
#install.packages("lmtest", dependencies = TRUE)Librerías
library(tidyverse)── Attaching packages ───────────────────────────────────────────── tidyverse 1.2.1 ──✔ ggplot2 3.2.1 ✔ purrr 0.3.3
✔ tibble 2.1.3 ✔ dplyr 0.8.3
✔ tidyr 1.0.0 ✔ stringr 1.4.0
✔ readr 1.3.1 ✔ forcats 0.4.0── Conflicts ──────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()library(sjlabelled)
Attaching package: 'sjlabelled'The following object is masked from 'package:forcats':
as_factorThe following object is masked from 'package:dplyr':
as_labelThe following objects are masked from 'package:haven':
as_factor, read_sas, read_spss, read_stata, write_sas,
zap_labelsIntro
En esta práctica vamos a revisar los elementos básicos para la regresión logística. El proceso en R para todos los modelos generalizados se parece mucho. Por tanto, no será difícil que luego puedas utilizar otras funciones de enlace. Vamos a hacer una sub-base de nuestras posibles variables explicativas. Esto es importante porque sólo podemos comparar modelos con la misma cantidad de observaciones. Intentaremos predecir el acceso a seguridad social
mydata<- SDEMT319 %>%
filter(clase2==1 & anios_esc<99 & eda<98) %>%
select("ing_x_hrs","t_loc","eda", "sex", "fac", "anios_esc", "imssissste")¿Cómo está codificada la variable “imssissste”
table(mydata$imssissste)
1 2 3 4 5
62321 11726 1171 102347 1309 table(as_label(mydata$imssissste))
No aplica Imss
0 62321
Issste Otras instituciones
11726 1171
No recibe atención médica No especificado
102347 1309 Vamos a volver dicotómica (0,1) nuestra variable [y de paso vemos cómo se recodifica en R]
mydata$y_binomial<-NA
mydata$y_binomial[mydata$imssissste<4]<-1
mydata$y_binomial[mydata$imssissste>3]<-0
newdata <- na.omit(mydata)Regresión Logística
Un solo predictor
modelo0<-glm(y_binomial ~ eda, family = binomial("logit"), data=newdata, na.action=na.exclude)
summary(modelo0)
Call:
glm(formula = y_binomial ~ eda, family = binomial("logit"), data = newdata,
na.action = na.exclude)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.193 -1.056 -0.943 1.281 1.626
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.1928322 0.0140398 13.73 <2e-16 ***
eda -0.0130982 0.0003385 -38.69 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 243432 on 178873 degrees of freedom
Residual deviance: 241912 on 178872 degrees of freedom
AIC: 241916
Number of Fisher Scoring iterations: 4confint(modelo0)Waiting for profiling to be done... 2.5 % 97.5 %
(Intercept) 0.16531856 0.22035380
eda -0.01376199 -0.01243502Predicción de probabilidades
Para predecir la probabilidad, primero chequemos el rango de nuestra variabe explicativa
range(newdata$eda)[1] 12 97Hacemos un vector con los valores que queremos predecir
xeda <- 15:97Vamos a utilizar el comando “predict” para predecir los valores. Podemos el argumento “response” para que nos dé el logito
y_logito <- predict(modelo0, list(eda = xeda))
y_prob<- predict(modelo0, list(eda = xeda), type= "response")
results_m0<-cbind(y_logito, y_prob, xeda)
results_m0<-as.data.frame(results_m0)Hoy podemos graficar
library(ggplot2)
ggplot(data=results_m0, aes(x=xeda, y=y_prob)) +
geom_point()
##Coeficientes exponenciados
Para interpretar mejor los coeficientes suelen exponenciarse y hablar de las veces que aumentan o disminuyen los momios con respecto a la unidad como base
exp(coef(modelo0))(Intercept) eda
1.2126793 0.9869872 Agregando una variable
modelo1<-glm(y_binomial ~ eda + as_label(t_loc), family = binomial("logit"), data=newdata, na.action=na.exclude)
summary(modelo1)
Call:
glm(formula = y_binomial ~ eda + as_label(t_loc), family = binomial("logit"),
data = newdata, na.action = na.exclude)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.3638 -1.0924 -0.6755 1.1602 2.1566
Coefficients:
Estimate
(Intercept) 0.6002679
eda -0.0143507
as_label(t_loc)Localidades de 15 000 a 99 999 habitantes -0.5198754
as_label(t_loc)Localidades de 2 500 a 14 999 habitantes -0.9557794
as_label(t_loc)Localidades menores de 2 500 habitantes -1.5457069
Std. Error
(Intercept) 0.0152367
eda 0.0003533
as_label(t_loc)Localidades de 15 000 a 99 999 habitantes 0.0147286
as_label(t_loc)Localidades de 2 500 a 14 999 habitantes 0.0162260
as_label(t_loc)Localidades menores de 2 500 habitantes 0.0177512
z value Pr(>|z|)
(Intercept) 39.40 <2e-16
eda -40.62 <2e-16
as_label(t_loc)Localidades de 15 000 a 99 999 habitantes -35.30 <2e-16
as_label(t_loc)Localidades de 2 500 a 14 999 habitantes -58.90 <2e-16
as_label(t_loc)Localidades menores de 2 500 habitantes -87.08 <2e-16
(Intercept) ***
eda ***
as_label(t_loc)Localidades de 15 000 a 99 999 habitantes ***
as_label(t_loc)Localidades de 2 500 a 14 999 habitantes ***
as_label(t_loc)Localidades menores de 2 500 habitantes ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 243432 on 178873 degrees of freedom
Residual deviance: 230323 on 178869 degrees of freedom
AIC: 230333
Number of Fisher Scoring iterations: 4confint(modelo1)Waiting for profiling to be done... 2.5 %
(Intercept) 0.57041722
eda -0.01504345
as_label(t_loc)Localidades de 15 000 a 99 999 habitantes -0.54876554
as_label(t_loc)Localidades de 2 500 a 14 999 habitantes -0.98763481
as_label(t_loc)Localidades menores de 2 500 habitantes -1.58059781
97.5 %
(Intercept) 0.63014436
eda -0.01365858
as_label(t_loc)Localidades de 15 000 a 99 999 habitantes -0.49102975
as_label(t_loc)Localidades de 2 500 a 14 999 habitantes -0.92402889
as_label(t_loc)Localidades menores de 2 500 habitantes -1.51101285Este modelo tiene coeficientes que deben leerse “condicionados”. Es decir, en este caso tenemos que el coeficiente asociado a la eda, mantiene constante el valor del “t_loc”" y viceversa.
Bondad de Ajuste
Devianza
La devianza es una medida de la bondad de ajuste de los modelos lineales generalizados. O más bien, es una medida de la no-bondad del ajust, puesto que valores más altos indican un peor ajuste.
R nos da medidas de devianza: la devianza nula y la desviación residual. La devianza nula muestra qué tan bien la variable de respuesta se predice mediante un modelo que incluye solo la intersección (gran media).
Prueba de Verosimilitud
image
¿Cómo saber si ha mejorado nuestro modelo? Podemos hacer un test que compare las devianzas(tendría la misma lógica que nuestra prueba F del modelo lineal). Para esto tenemos que instalar un paquete “lmtest”
library(lmtest)Loading required package: zoo
Attaching package: 'zoo'The following objects are masked from 'package:base':
as.Date, as.Date.numericlrtest0<-lrtest(modelo0, modelo1)
lrtest0Likelihood ratio test
Model 1: y_binomial ~ eda
Model 2: y_binomial ~ eda + as_label(t_loc)
#Df LogLik Df Chisq Pr(>Chisq)
1 2 -120956
2 5 -115162 3 11589 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Como puedes ver, el resultado muestra un valor p muy pequeño (<.001). Esto significa que agregar el t_loc al modelo lleva a un ajuste significativamente mejor sobre el modelo original.
Podemos seguir añadiendo variables sólo “sumando” en la función
modelo2<-glm(y_binomial ~ eda + as_label(t_loc) + as_label(sex), family = binomial("logit"), data=newdata, na.action=na.exclude)
summary(modelo2)
Call:
glm(formula = y_binomial ~ eda + as_label(t_loc) + as_label(sex),
family = binomial("logit"), data = newdata, na.action = na.exclude)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.3772 -1.0930 -0.6736 1.1596 2.1473
Coefficients:
Estimate
(Intercept) 0.6305440
eda -0.0143513
as_label(t_loc)Localidades de 15 000 a 99 999 habitantes -0.5207803
as_label(t_loc)Localidades de 2 500 a 14 999 habitantes -0.9585720
as_label(t_loc)Localidades menores de 2 500 habitantes -1.5535996
as_label(sex)Mujer -0.0711363
Std. Error
(Intercept) 0.0158406
eda 0.0003532
as_label(t_loc)Localidades de 15 000 a 99 999 habitantes 0.0147316
as_label(t_loc)Localidades de 2 500 a 14 999 habitantes 0.0162341
as_label(t_loc)Localidades menores de 2 500 habitantes 0.0177917
as_label(sex)Mujer 0.0101256
z value Pr(>|z|)
(Intercept) 39.806 < 2e-16
eda -40.627 < 2e-16
as_label(t_loc)Localidades de 15 000 a 99 999 habitantes -35.351 < 2e-16
as_label(t_loc)Localidades de 2 500 a 14 999 habitantes -59.047 < 2e-16
as_label(t_loc)Localidades menores de 2 500 habitantes -87.322 < 2e-16
as_label(sex)Mujer -7.025 2.14e-12
(Intercept) ***
eda ***
as_label(t_loc)Localidades de 15 000 a 99 999 habitantes ***
as_label(t_loc)Localidades de 2 500 a 14 999 habitantes ***
as_label(t_loc)Localidades menores de 2 500 habitantes ***
as_label(sex)Mujer ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 243432 on 178873 degrees of freedom
Residual deviance: 230274 on 178868 degrees of freedom
AIC: 230286
Number of Fisher Scoring iterations: 4confint(modelo2)Waiting for profiling to be done... 2.5 %
(Intercept) 0.59951015
eda -0.01504404
as_label(t_loc)Localidades de 15 000 a 99 999 habitantes -0.54967626
as_label(t_loc)Localidades de 2 500 a 14 999 habitantes -0.99044335
as_label(t_loc)Localidades menores de 2 500 habitantes -1.58856965
as_label(sex)Mujer -0.09098502
97.5 %
(Intercept) 0.66160438
eda -0.01365932
as_label(t_loc)Localidades de 15 000 a 99 999 habitantes -0.49192881
as_label(t_loc)Localidades de 2 500 a 14 999 habitantes -0.92680569
as_label(t_loc)Localidades menores de 2 500 habitantes -1.51882617
as_label(sex)Mujer -0.05129317Y podemos ver si introducir esta variable afectó al ajuste global del modelo
lrtest1<-lrtest(modelo1, modelo2)
lrtest1Likelihood ratio test
Model 1: y_binomial ~ eda + as_label(t_loc)
Model 2: y_binomial ~ eda + as_label(t_loc) + as_label(sex)
#Df LogLik Df Chisq Pr(>Chisq)
1 5 -115162
2 6 -115137 1 49.403 2.085e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Test Hosmer-Lemeshow Goodness of Fit “GOF”
El teste Homer-Lemeshow se calcula sobre los datos una vez que las observaciones se han segmentado en grupos basados en probabilidades predichas similares. Este teste examina si las proporciones observadas de eventos son similares a las probabilidades predichas de ocurrencia en subgrupos del conjunto de datos, y lo hace con una prueba de chi cuadrado de Pearson.
¡Ojo! No queremos rechazar la hipótesis nula. LoLa hipótesis nula sostiene que el modelo se ajusta a los datos y en el siguiente ejemplo rechazamos H0.
library(ResourceSelection)ResourceSelection 0.3-5 2019-07-22hoslem.test(newdata$imssissste, fitted(modelo2))
Hosmer and Lemeshow goodness of fit (GOF) test
data: newdata$imssissste, fitted(modelo2)
X-squared = 5330330, df = 8, p-value < 2.2e-16No obstante, esta prueba ha sido criticada. Checa la postura de Paul Allison https://statisticalhorizons.com/hosmer-lemesho Es un problema que tenemos en muestras grandes. Casi siempre preferimos el enfoque de la devianza.
Tabla de modelos estimados
library(stargazer)#stargazer(modelo0, modelo1,modelo2, type = 'latex', header=FALSE)stargazer(modelo0, modelo1,modelo2,
type = 'text', header=FALSE)
===============================================================================================
Dependent variable:
--------------------------------------
y_binomial
(1) (2) (3)
-----------------------------------------------------------------------------------------------
eda -0.013*** -0.014*** -0.014***
(0.0003) (0.0004) (0.0004)
as_label(t_loc)Localidades de 15 000 a 99 999 habitantes -0.520*** -0.521***
(0.015) (0.015)
as_label(t_loc)Localidades de 2 500 a 14 999 habitantes -0.956*** -0.959***
(0.016) (0.016)
as_label(t_loc)Localidades menores de 2 500 habitantes -1.546*** -1.554***
(0.018) (0.018)
as_label(sex)Mujer -0.071***
(0.010)
Constant 0.193*** 0.600*** 0.631***
(0.014) (0.015) (0.016)
-----------------------------------------------------------------------------------------------
Observations 178,874 178,874 178,874
Log Likelihood -120,956.200 -115,161.700 -115,137.000
Akaike Inf. Crit. 241,916.400 230,333.500 230,286.100
===============================================================================================
Note: *p<0.1; **p<0.05; ***p<0.01Para sacar los coeficientes exponenciados
stargazer(modelo0, modelo1,modelo2,
type = 'text', header=FALSE,
apply.coef = exp,
apply.se = exp)
===============================================================================================
Dependent variable:
--------------------------------------
y_binomial
(1) (2) (3)
-----------------------------------------------------------------------------------------------
eda 0.987 0.986 0.986
(1.000) (1.000) (1.000)
as_label(t_loc)Localidades de 15 000 a 99 999 habitantes 0.595 0.594
(1.015) (1.015)
as_label(t_loc)Localidades de 2 500 a 14 999 habitantes 0.385 0.383
(1.016) (1.016)
as_label(t_loc)Localidades menores de 2 500 habitantes 0.213 0.211
(1.018) (1.018)
as_label(sex)Mujer 0.931
(1.010)
Constant 1.213 1.823* 1.879*
(1.014) (1.015) (1.016)
-----------------------------------------------------------------------------------------------
Observations 178,874 178,874 178,874
Log Likelihood -120,956.200 -115,161.700 -115,137.000
Akaike Inf. Crit. 241,916.400 230,333.500 230,286.100
===============================================================================================
Note: *p<0.1; **p<0.05; ***p<0.01También la librería “sjPlot” tiene el comando “plot_model()” (instala el comando si no lo tienes)
library(sjPlot)Registered S3 methods overwritten by 'survival':
method from
nobs.coxph insight
nobs.survreg insightplot_model(modelo1)Waiting for profiling to be done...
Por default nos da los coeficientes exponenciados.
plot_model(modelo1, type="eff")Can't compute marginal effects, 'effects::Effect()' returned an error.
Reason: contrasts can be applied only to factors with 2 or more levels
You may try 'ggpredict()' or 'ggemmeans()'.
Can't compute marginal effects, 'effects::Effect()' returned an error.
Reason: non-conformable arguments
You may try 'ggpredict()' or 'ggemmeans()'.NULL¿Cómo saber lo que tiene esos gráficos? Es bueno guardar siempre estos resultados en un objeto. Este objeto es una lista de dos listas
get<-plot_model(modelo1, type="eff")Can't compute marginal effects, 'effects::Effect()' returned an error.
Reason: contrasts can be applied only to factors with 2 or more levels
You may try 'ggpredict()' or 'ggemmeans()'.
Can't compute marginal effects, 'effects::Effect()' returned an error.
Reason: non-conformable arguments
You may try 'ggpredict()' or 'ggemmeans()'.get$edaNULLget$eda$dataNULLget$t_locNULLget$t_loc$dataNULL