Regresión Logística
Pablo Gaitán
2022-03-29
Limpieza de datos
Empezamos por cargar paquetes, datos, y por crear el ID
library(tidyverse)
library(skimr)
library(sjPlot)
library(sjlabelled)
library(sjmisc)
library(broom)
library(effectsize)
library(performance)
library(see)
library(ggmosaic)
library(modelbased)
library(readr)
TSDEM <- read_csv("Datos/TSDEM.csv")
TSDEM1 <- TSDEM %>%
mutate(HOGAR= as.character(HOGAR)) %>%
mutate(id_personal= paste(UPM_DIS, VIV_SEL, HOGAR, N_REN, sep=".")) # ID
TENBIARE <- read_csv("Datos/TENBIARE.csv")
TENBIARE1 <- TENBIARE %>%
mutate(HOGAR= as.character(HOGAR)) %>%
mutate(id_personal= paste(UPM_DIS, VIV_SEL, HOGAR, N_REN, sep=".")) # ID
enbiare <- TENBIARE1 %>% # nuevo objeto bd
left_join(TSDEM1, by= "id_personal")Recuerden que estaré siguiendo estos dos tutoriales. Unos sobre la interpretación de probabilidad, momios y razón de momios y el otro más directamente sobre regresión logística, ambos cortesía del OARC de UCLA: siempre muy recomendable.
Depresión
La variable dependiente es depresión medida con CESD-7.
enbiare <- enbiare %>%
mutate(PD2_6= case_when(
PD2_6== 0 ~ 3,
PD2_6== 1 ~ 2,
PD2_6== 2 ~ 1,
PD2_6== 3 ~ 0))
enbiare <- enbiare %>%
mutate(d_total = PD2_1 + PD2_2 + PD2_3 + PD2_4 + PD2_5 + PD2_6 + PD2_7)
summary(enbiare$d_total)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 1.000 3.000 4.373 7.000 21.000Así se ve su histograma
enbiare %>%
ggplot(aes(x= d_total))+
geom_histogram(binwidth = 1, fill= "darkorange", color="white")+
geom_vline(xintercept = mean(enbiare$d_total), color= "blue")+
geom_vline(xintercept = 9, color= "black")+ # punto de corte
theme_minimal()
Pero a nosotros nos importa convertirla a dummy (dep_level):
Sin síntomas= 0
Con síntomas = 1
enbiare <- enbiare %>%
mutate(dep_level= if_else(d_total >= 9, 1, 0))
enbiare <- enbiare %>%
mutate(dep_level2= ordered(dep_level, levels=c("0", "1"), labels= c("Sin Síntomas", "Con Síntomas")))
table(enbiare$dep_level2, useNA = "ifany")##
## Sin Síntomas Con Síntomas
## 26050 5116prop.table(table(enbiare$dep_level2, useNA = "ifany"))*100 # ordinal##
## Sin Síntomas Con Síntomas
## 83.58468 16.41532prop.table(table(enbiare$dep_level, useNA = "ifany"))*100 # numerica##
## 0 1
## 83.58468 16.41532La prevalencia es de 16.4%. Esto viene de la probabilidad de tener depresión:
5116/ (26050 + 5116) = 0.164
5116/ (26050 + 5116) ## [1] 0.1641532¿Cuál es el odds? fórmula= p/ 1-p
1- 0.164## [1] 0.8360.164/ (1- 0.164)## [1] 0.1961722Y ahora el log de los odds
log(0.1961722)## [1] -1.628762Así se ve el modelo nulo de la regresión logística. Y ahí nuestro log odds!
glm(dep_level ~ 1, data = enbiare, family = "binomial")##
## Call: glm(formula = dep_level ~ 1, family = "binomial", data = enbiare)
##
## Coefficients:
## (Intercept)
## -1.628
##
## Degrees of Freedom: 31165 Total (i.e. Null); 31165 Residual
## Null Deviance: 27830
## Residual Deviance: 27830 AIC: 27830Si este es el primer resultado que tenemos, ¿cómo volvemos?
exp(-1.628)/(1+exp(-1.628))## [1] 0.1641045Y ahí está nuestro 16%!
Sexo
Empecemos por un control dicotómico
enbiare <- enbiare %>%
mutate(mujer= case_when(
SEXO== 1 ~ 0,
SEXO== 2 ~ 1))
enbiare <- enbiare %>%
mutate(mujer_ord= ordered(mujer, levels = c("0", "1" ), labels=c("Hombre", "Mujer")))
table(enbiare$mujer_ord, useNA = "ifany")##
## Hombre Mujer
## 14255 16911prop.table(table(enbiare$mujer_ord, enbiare$dep_level2), margin = 1)*100##
## Sin Síntomas Con Síntomas
## Hombre 88.32690 11.67310
## Mujer 79.58725 20.41275chisq.test(table(enbiare$mujer_ord, enbiare$dep_level2))##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(enbiare$mujer_ord, enbiare$dep_level2)
## X-squared = 429.96, df = 1, p-value < 2.2e-16enbiare %>%
ggplot()+
geom_mosaic(aes(x = product(mujer_ord), fill = dep_level2))+
theme_classic()
Ahora veamos los odds combinados
addmargins(table(enbiare$mujer_ord, enbiare$dep_level2))##
## Sin Síntomas Con Síntomas Sum
## Hombre 12591 1664 14255
## Mujer 13459 3452 16911
## Sum 26050 5116 31166¿Cuáles son los odds de que un hombre tenga depresión y los odds de que una mujer la tenga?
(1664/ 14255) / (12591/ 14255) # odds dep hombres## [1] 0.13215791664/12591 # odds dep hombres## [1] 0.1321579(3452/ 16911) / (13459/ 16911) # odds dep mujeres## [1] 0.25648273452 / 13459 # odds dep mujeres## [1] 0.2564827¿Cuál es el odds ratio (OR) hombres a mujeres ?
(1664 / 12591) / (3452 / 13459) # OR## [1] 0.51527030.1321579 / 0.2564827 # OR## [1] 0.5152702¿Cómo interpretamos un odds ratio (OR) cuando es menor a 1?
1- 0.515## [1] 0.485Los momios de que un hombre tenga depresión son 48% menos a que una mujer tenga depresión.
Y así se ve el log odds
log(0.5152702)## [1] -0.6630639Ahora veamos la regresión logística la cual da el output en log odds
mod1 <- glm(dep_level ~ mujer, data = enbiare, family = "binomial")
summary(mod1)##
## Call:
## glm(formula = dep_level ~ mujer, family = "binomial", data = enbiare)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.6757 -0.6757 -0.4983 -0.4983 2.0726
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.02376 0.02608 -77.59 <2e-16 ***
## mujer 0.66306 0.03232 20.52 <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: 27831 on 31165 degrees of freedom
## Residual deviance: 27390 on 31164 degrees of freedom
## AIC: 27394
##
## Number of Fisher Scoring iterations: 4Eso se ve bien, para ya sabemos que nos es muy difícil interpretar logs. Vamos a pasarlo a OR:
coef(mod1)## (Intercept) mujer
## -2.0237579 0.6630637exp(coef(mod1)) # Estos son los momios## (Intercept) mujer
## 0.1321579 1.9407290Los momios de que una mujer tenga depresión son 94% mayores a que un hombre tenga depresión
Salió de aquí.
0.2564827/ 0.1321579 # OR## [1] 1.9407291.94 - 1## [1] 0.94No pierdan de vista qué valor va arriba, en el numerador, porque ese guía la interpretación. Veamos la interpretación inversa.
0.1321579 / 0.2564827## [1] 0.5152702 1 - 0.515## [1] 0.485En este caso debemos restar el 1.
Eso significa que los momios de que un hombre tenga depresión son 48.5% menores a que una mujer tenga depresión.
Y en probabilidad?
exp(0.6630637)/(1+exp(0.6630637))## [1] 0.6599483Las mujeres tienen un 65.9% mayor probabilidad de tener depresión que los varones (sin ajustar)
Préstamos
enbiare <- enbiare %>%
mutate_at(
vars(starts_with("PF1_")),
funs(case_when(
.== 1 ~ 1,
.== 2 ~ 0)))## Warning: `funs()` was deprecated in dplyr 0.8.0.
## Please use a list of either functions or lambdas:
##
## # Simple named list:
## list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`:
## tibble::lst(mean, median)
##
## # Using lambdas
## list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.enbiare <- enbiare %>%
mutate(prest_cont = PF1_1 + PF1_2 + PF1_3 + PF1_4 + PF1_5 + PF1_6)
table(enbiare$prest_cont, useNA = "ifany")##
## 0 1 2 3 4 5 6
## 19609 3385 2563 2156 1927 1064 462prop.table(table(enbiare$prest_cont, useNA = "ifany"))*100##
## 0 1 2 3 4 5 6
## 62.917923 10.861195 8.223705 6.917795 6.183020 3.413977 1.482385Vamos a ver sus asociaciones: ¿ven un patrón ascendente?
enbiare %>%
group_by(prest_cont) %>%
summarise(medias= mean(dep_level)) %>%
ungroup()## # A tibble: 7 × 2
## prest_cont medias
## <dbl> <dbl>
## 1 0 0.119
## 2 1 0.188
## 3 2 0.208
## 4 3 0.260
## 5 4 0.297
## 6 5 0.303
## 7 6 0.338Vean cómo voy a graficar esta tabla
enbiare %>%
group_by(prest_cont) %>%
summarise(medias= mean(dep_level)) %>%
ggplot(aes(x= prest_cont, y= medias)) +
geom_bar(stat='identity', fill= "skyblue" )+
geom_hline(yintercept = 0.16, color= "red")+ # media de depresion
theme_classic()
Una gráfica pero viendo la media de préstamos
enbiare %>%
ggplot( aes(x= dep_level2, y=prest_cont)) +
geom_boxplot()+
geom_jitter(color="gray", size=0.4, alpha=0.3) +
stat_summary(fun.y="mean", color= "red")+
theme_classic()## Warning: `fun.y` is deprecated. Use `fun` instead.## Warning: Removed 2 rows containing missing values (geom_segment).
Vean cómo, cada préstamo adicional, va subiendo la probabilidad de estar en el cuadro de arriba. Ese aumento en probabilidad lo indica la línea azul.
enbiare %>%
ggplot(aes(x=prest_cont, y=dep_level))+
geom_point() +
geom_jitter(color="gray", size=0.4, alpha=0.3)+
stat_smooth(method="glm", se=TRUE, method.args = list(family=binomial))+
theme_classic()## `geom_smooth()` using formula 'y ~ x'
Sí vemos asociaciones, pero ahora el reto es examinar si esas asociaciones se mantienen una vez que tenemos los controles. Ahora agrupamos por sexo….
enbiare %>%
ggplot(aes(x=prest_cont, y=dep_level))+
geom_point() +
geom_jitter(color="gray", size=0.4, alpha=0.3)+
stat_smooth(method="glm", se=TRUE, method.args = list(family=binomial))+
facet_wrap(~ mujer_ord) +
theme_classic()## `geom_smooth()` using formula 'y ~ x'
Ahora voy a ir más rápido.
Educación
enbiare <- enbiare %>%
mutate(educ= case_when(
NIVEL== "00" ~ 0,
NIVEL== "01" ~ 0,
NIVEL== "02" ~ 0, # Primaria
NIVEL== "03" ~ 1,
NIVEL== "04" ~ 1,
NIVEL== "05" ~ 1,
NIVEL== "06" ~ 1, # Prepa
NIVEL== "07" ~ 2,
NIVEL== "08" ~ 2,
NIVEL== "09" ~ 2,
NIVEL== "10" ~ 2, # Superior
NIVEL== "99" ~ NA_real_))
enbiare <- enbiare %>%
mutate(educ_ord= ordered(educ, levels = c("0", "1", "2"), labels=c("Primaria", "Preparatoria", "Superior")))
table(enbiare$educ_ord, useNA = "ifany")##
## Primaria Preparatoria Superior <NA>
## 7658 15898 7549 61prop.table(table(enbiare$educ_ord, enbiare$dep_level2), margin = 1)*100##
## Sin Síntomas Con Síntomas
## Primaria 77.10891 22.89109
## Preparatoria 84.18040 15.81960
## Superior 88.89919 11.10081chisq.test(table(enbiare$educ_ord, enbiare$dep_level2))##
## Pearson's Chi-squared test
##
## data: table(enbiare$educ_ord, enbiare$dep_level2)
## X-squared = 393.56, df = 2, p-value < 2.2e-16enbiare %>%
ggplot()+
geom_mosaic(aes(x = product(educ_ord), fill = dep_level2))+
theme_classic()
Desempleo
enbiare <- enbiare %>%
mutate(desempleo= case_when(
PE15== 1 ~ 1, # sí
PE15== 2 ~ 0, # Sí, lo recuperó
PE15== 3 ~ 0)) # No
enbiare <- enbiare %>%
mutate(desempleo_ord= ordered(desempleo, levels = c("0", "1"),
labels=c("No", "Sí")))
table(enbiare$desempleo_ord, useNA = "ifany")##
## No Sí
## 26866 4300prop.table(table(enbiare$desempleo_ord, enbiare$dep_level2), margin = 1)*100##
## Sin Síntomas Con Síntomas
## No 84.63113 15.36887
## Sí 77.04651 22.95349chisq.test(table(enbiare$desempleo_ord, enbiare$dep_level2))##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(enbiare$desempleo_ord, enbiare$dep_level2)
## X-squared = 154.86, df = 1, p-value < 2.2e-16enbiare %>%
ggplot()+
geom_mosaic(aes(x = product(desempleo_ord), fill = dep_level2))+
theme_classic()
Regresión logística
Lo primero, como habrán notado ya, es que cambia el comando en dos puntos clave. El primero es que ya no corremos un modelo de regresión lineal (“lm”) sino uno de regresión logística que se inscribe en los generalized linear models: “glm”. Al transformar la var dependiente (dicotómica - > odds -> odds ratio -> log odds ratio), ya no nos basamos en la distribución normal, o gaussiana sino que usamos una distribución binomial - la cual muestra la distribución de probabilidad de una serie de éxitos de eventos. Si sólo hablamos de un evento, se usa la distribución Bernoulli, si son una serie de eventos (como en las regresiones) entonces usamos la distribución binomial. También cambiamos el método de estimación; pasamos de mínimos cuadrados ordinarios (OLS) con el comando lm a uno de máxima verosimilitud.
Veamos ahora un regresión logística múltiple
mod2 <- glm(dep_level ~ prest_cont + mujer + educ + desempleo, data = enbiare, family = "binomial")
summary(mod2)##
## Call:
## glm(formula = dep_level ~ prest_cont + mujer + educ + desempleo,
## family = "binomial", data = enbiare)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3357 -0.6362 -0.5116 -0.4047 2.3812
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.967494 0.035880 -54.835 < 2e-16 ***
## prest_cont 0.228485 0.008976 25.454 < 2e-16 ***
## mujer 0.647771 0.033004 19.627 < 2e-16 ***
## educ -0.403496 0.023255 -17.351 < 2e-16 ***
## desempleo 0.313546 0.042456 7.385 1.52e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 27776 on 31104 degrees of freedom
## Residual deviance: 26178 on 31100 degrees of freedom
## (61 observations deleted due to missingness)
## AIC: 26188
##
## Number of Fisher Scoring iterations: 4Este output está en log odds. Algo que queda claro es que todos los coeficientes son estadísticamente significativos (vean los p-values). Nuestro problema es la interpretación porque es poco informativo para los mortales que no pensamos en logaritmos. Veamos sólo el coef de préstamos:
- Por cada préstamo adicional (i.e., unidad adicional en X1), el logaritmo de los momios de depresión aumenta en 0.228.
A nosotros lo que nos gusta son los odds ratio:
exp(coef(mod2))## (Intercept) prest_cont mujer educ desempleo
## 0.1398067 1.2566941 1.9112767 0.6679810 1.3682687Ahora sí vamos a mirar cada coeficiente:
El intercepto ahora nos nos interesa, pero podríamos sacar la probabilidad condicional de tener depresión.
Prestamos: Por cada préstamo adicional, los momios de tener depresión aumentan en un 25%, controlando por sexo, educación y desempleo.
Sexo: Los momios de las mujeres de tener depresión son 91% mayores que los de los hombres, controlando por préstamos, educación y desempleo.
Eduación. Por cada nivel eduactivo adicional, los momios de tener depresión disminuyen en un 33% (1- 0.667), controlando por préstamos, sexo, y desempleo.
Desempleo: Los momios de las personas desempleadas de tener depresión son 36% mayores que los de los empleados, controlando por préstamos, sexo, y educación.
Y, mejor aún, si nos da los intervalos de confianza con los que construimos los aviones de las gráficas.
exp(cbind(OR =coef(mod2), confint(mod2)))## Waiting for profiling to be done...## OR 2.5 % 97.5 %
## (Intercept) 0.1398067 0.1302684 0.1499440
## prest_cont 1.2566941 1.2347508 1.2789751
## mujer 1.9112767 1.7918378 2.0393563
## educ 0.6679810 0.6381742 0.6990898
## desempleo 1.3682687 1.2585597 1.4864845Graficas del modelo
Empecemos por los coeficientes con sus intervalos de confianza (“los aviones”). Noten que ya están en Odds Ratios porque eso es lo que se suele reportar. De igual manera, recuerden que la línea de referencia es el 1. Si lo quisieran en en log odds, tendrían que poner el argumento “transform = NULL” , pero no creo que sea necesario.
plot_model(mod2)
Y fíjense cómo lo puedo combinar con ggplot para mejorar la presentación y hacerla más profesional.
plot_model(mod2, vline.color = "red", show.values = TRUE, value.offset = .3)+
labs(
y= "Razón de momios",
x= "Variables",
fill= "",
title= "Relación entre depresión y préstamos",
subtitle= "Ejercicio de clase",
caption= "Realizado con el paquete sjplots")+
theme_classic()
Ahora vamos a ver los efectos marginales con nuestra variable de tratamiento.
Vean cómo nos la de en… probabilidades! Y con su elegante sombreado indicando los intervalos de confianza
plot_model(mod2, type = "pred", terms = "prest_cont")+
theme_classic()
Y sí, adivinaron, podemos meter otras variables para contar historias con más textura.
plot_model(mod2, type = "pred", terms = c("prest_cont", "mujer"))+
theme_classic()
plot_model(mod2, type = "pred", terms = c("prest_cont", "educ", "mujer"))+
theme_classic()
Predicciones del modelo
Usemos el poder de la predicción… Ahora metemos una nueva base de datos; una en la cual desconocemos el valor de la var dependiente (i.e. depresión) y justamente, a queremos anticipar o predecir. Noten la composición, la cual ayuda a comparar probabilidades de acuerdo a perfiles de personas específicas. (Ignoren la “L”, al usar DataPasta y pegarla como tibble se añade la L, que significa integer). Estamos tratando de modelar una mujer, de educación prepa y empleada. La única diferencia entre ellas 7 es el número de préstamos.
nueva_bd <- tibble::tribble(
~prest_cont, ~mujer, ~educ, ~desempleo,
0L, 1L, 1L, 0L,
1L, 1L, 1L, 0L,
2L, 1L, 1L, 0L,
3L, 1L, 1L, 0L,
4L, 1L, 1L, 0L,
5L, 1L, 1L, 0L,
6L, 1L, 1L, 0L
)Y vean qué chulada… Estas son las probabilidades que predice el modelo. Probabilidades de las que nos gustan, que van de 0-1. El numero de la primera fila indica la fila de la nueva base que ingresamos. Es la mujer 1, luego la mujer 2, etc.; abajo aparece su probabilidad correspondiente.
predicciones <- predict(mod2, newdata = nueva_bd, type = "response")
predicciones## 1 2 3 4 5 6 7
## 0.1514571 0.1832123 0.2199000 0.2615815 0.3080439 0.3587495 0.4128221Una mujer empleada y con preparatoria, sin préstamos, tiene un 15% de probabilidad de tener depresión. Si una mujer equivalente pide un préstamo, ahora su probabilidad es de 18%. Si pide dos préstamos, es de 23%. Si pide 6 es de 36%!
Esta función nos permite calcular los valores para toda la muestra.
pred1 <- estimate_response(mod2)## `estimate_response()` is deprecated.
## Please use `estimate_expectation()` (for conditional expected values) or
## `estimate_prediction()` (for individual case predictions) instead.head(pred1)## Model-based Expectation
##
## dep_level | prest_cont | mujer | educ | desempleo | Predicted | SE | 95% CI | Residuals
## ---------------------------------------------------------------------------------------------------
## 0.00 | 0.00 | 0.00 | 2.00 | 0.00 | 0.06 | 2.21e-03 | [0.05, 0.06] | -0.06
## 0.00 | 1.00 | 1.00 | 0.00 | 0.00 | 0.25 | 5.42e-03 | [0.24, 0.26] | -0.25
## 0.00 | 0.00 | 1.00 | 2.00 | 0.00 | 0.11 | 3.34e-03 | [0.10, 0.11] | -0.11
## 0.00 | 0.00 | 1.00 | 0.00 | 0.00 | 0.21 | 5.18e-03 | [0.20, 0.22] | -0.21
## 0.00 | 0.00 | 1.00 | 1.00 | 0.00 | 0.15 | 3.05e-03 | [0.15, 0.16] | -0.15
## 0.00 | 0.00 | 1.00 | 0.00 | 0.00 | 0.21 | 5.18e-03 | [0.20, 0.22] | -0.21
##
## Variable predicted: dep_levelPrecisión en la estimación
Voy a añadir un tema que en realidad debe de tomar más tiempo y más cuidado. Pero es final de semestre y prefiero que al menos lo vean ustedes más a profundidad después. La ENBIARE, como la mayoría de las encuestas del INEGI, es una encuesta con un diseño muestral complejo. Esto significa que van a aleatorizando por niveles (ej estados, UPM, etc.). Y por razones que exceden este espacio, ponderan algunas observaciones más que otras. Es decir, algunas personas con perfiles raros, por ser más difíciles de encontrar, valen más y, en cambio, perfiles de personas muy comunes, valen un poco menos porque hay muchos en la muestra. al final lo que logra el inegi con estos ponderadores muestrales es que la composición de la muestra sea similar a la composición de la población mexicana. en términos prácticos, al usar restos datos, si queremos obtener lo mismo que el INEGI, entonces SIEMPRE deberíamos estimar tanto descriptivos como modelos utilizando estos ponderadores.
En esta base la variable con los ponderadores muestrales es la de FAC_ELE, también conocida como factores de expansión. si una persona tiene un factor de 122, su respuesta equivale a la de 122 personas. Si una persona tiene un factor de 2,710 , su respuesta equivale a 2,710.
summary(enbiare$FAC_ELE)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 122 949 1779 2710 3248 45955Si vemos las observaciones en la base, tenemos 5,116 personas con síntomas de depresión.
table(enbiare$dep_level2)##
## Sin Síntomas Con Síntomas
## 26050 5116Pero… si lo estimamos con factores de expansión… tenemos 12, 967, 919. Increíble, yo sé.
enbiare %>%
count(dep_level2, wt = FAC_ELE)## # A tibble: 2 × 2
## dep_level2 n
## <ord> <dbl>
## 1 Sin Síntomas 71482017
## 2 Con Síntomas 12967919Ahora comparemos descriptivos. La muestra sola nos da 16.4%
enbiare %>%
summarise(mean(dep_level)*100)## # A tibble: 1 × 1
## `mean(dep_level) * 100`
## <dbl>
## 1 16.4Pero… con los pesos ya nos da 15.35% Se parece pero no es igual. Esta es la prevalencia que deberíamos reportar. Quien vaya a trabajar con descriptivos de bases del INEGI le recomiendo este tutorial de svyr y la bibliografía asociada, en especial Lumley.
enbiare %>%
summarise(weighted.mean(dep_level, FAC_ELE)*100)## # A tibble: 1 × 1
## `weighted.mean(dep_level, FAC_ELE) * 100`
## <dbl>
## 1 15.4Volviendo a nuestra regresión… lo primero es decirle a R que es un diseño complejo. Recuerden que el sufijo .x es por cómo las pegamos. Idealmente, no debe llevar eso.
library(srvyr)##
## Attaching package: 'srvyr'## The following object is masked from 'package:stats':
##
## filterenbiare_survey <- enbiare %>%
as_survey_design(ids = UPM_DIS.x, # cluster
weights = FAC_ELE, #factor
strata = EST_DIS.x) # estratoAhora ya tenemos un espejo de nuestra base, pero con diseño complejo. Van estadísticos como deberían de calcularse.
enbiare_survey %>%
summarize(survey_mean(dep_level, vartype = "ci"))## # A tibble: 1 × 3
## coef `_low` `_upp`
## <dbl> <dbl> <dbl>
## 1 0.154 0.148 0.160Por sexo
enbiare_survey %>%
group_by(mujer) %>%
summarize(survey_mean(dep_level, vartype = "ci"))## # A tibble: 2 × 4
## mujer coef `_low` `_upp`
## <dbl> <dbl> <dbl> <dbl>
## 1 0 0.107 0.100 0.114
## 2 1 0.195 0.186 0.204Por prestamos
enbiare_survey %>%
group_by(prest_cont) %>%
summarize(survey_mean(dep_level, vartype = "ci"))## # A tibble: 7 × 4
## prest_cont coef `_low` `_upp`
## <dbl> <dbl> <dbl> <dbl>
## 1 0 0.108 0.102 0.115
## 2 1 0.185 0.166 0.205
## 3 2 0.196 0.172 0.220
## 4 3 0.264 0.233 0.294
## 5 4 0.287 0.256 0.318
## 6 5 0.281 0.243 0.319
## 7 6 0.346 0.283 0.410Y, finalmente, nuestra regresión con ponderadores
library(survey)## Loading required package: grid## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## Loading required package: survival##
## Attaching package: 'survey'## The following object is masked from 'package:graphics':
##
## dotchartmod3 <- svyglm(dep_level ~ prest_cont + mujer + educ + desempleo, design = enbiare_survey, family = "binomial")## Warning in eval(family$initialize): non-integer #successes in a binomial glm!summary(mod3)##
## Call:
## svyglm(formula = dep_level ~ prest_cont + mujer + educ + desempleo,
## design = enbiare_survey, family = "binomial")
##
## Survey design:
## Called via srvyr
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.14466 0.05394 -39.761 < 2e-16 ***
## prest_cont 0.24638 0.01291 19.085 < 2e-16 ***
## mujer 0.71247 0.04964 14.352 < 2e-16 ***
## educ -0.34369 0.03485 -9.863 < 2e-16 ***
## desempleo 0.31838 0.06452 4.934 8.35e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 0.9958239)
##
## Number of Fisher Scoring iterations: 4Se parecen, pero no son iguales. ele efecto de educacion…
exp(coef(mod2))## (Intercept) prest_cont mujer educ desempleo
## 0.1398067 1.2566941 1.9112767 0.6679810 1.3682687exp(coef(mod3))## (Intercept) prest_cont mujer educ desempleo
## 0.1171080 1.2793904 2.0390246 0.7091457 1.3748929plot_models(mod2, mod3)
Calidad del modelo
Recuerden que todos los modelos tienen supuestos que deben cumplirse. En el caso de la regresión logística, al tener una variable dependiente dicotómica, los supuestos son un poco diferentes a los de la lineal.
Los supuestos básicos son:
1- Independencia de errores 2- Ausencia de Colinearidad 3- Ausencia de outliers
Les sugiero leer este artículo para profundizar en los supuestos y formas de atenderlos: Logistic regression: a brief primer.
Veamos la función global del paquete performance, que ya detecta que es un glm
check_model(mod2)
No se ve muy bien. Lo primero que nos dice es que no hay homogeneidad de varianza. Las predicciones son mejores para un grupo que para el otro. No es un modelo muy consistente. Más adelante lo exploramos un poco más. Esto nos pega en el supuesto 1. La cuarta gráfica también explora esto de otra manera y también nos dice que los residuales altos no son normales. Muy probablemente nos falta incluir variables que sean buenas para detectar a quienes suelen tener depresión.
La gráfica 2 habla del supuesto 2 y se ve bien.
La gráfica 3 también nos dice que no tenemos valores influyentes altos-problemáticos.
Con esto yo concluiría que el problema es que no estamos controlando por suficientes variables. Es un modelo incompleto.
Veamos una prueba formal de si tenemos homogeneidad de varianzas: nop.
check_homogeneity(mod2)## Warning: Variances differ between groups (Bartlett Test, p = 0.000).Veamos los residuales. No está ideal, pero creo que puede aguantar. Vean la nota
residuales <- binned_residuals(mod2)## Warning: About 92% of the residuals are inside the error bounds (~95% or higher would be good).residuales
Terminemos con una medición global del modelo. La idea es hacer una prueba de hipótesis donde comparamos nuestro modelo teórico (la regresión) con un modelo empírico (los datos). Si el valor p es mayor a 0.05 entonces no hay diferencia significativa, que es lo que queremos.
Valores p bajos nos dicen que nuestros datos cuentan una historia distinta a nuestro modelo.
performance_hosmer(mod2, n_bins = 4)## # Hosmer-Lemeshow Goodness-of-Fit Test
##
## Chi-squared: 2.004
## df: 2
## p-value: 0.367## Summary: model seems to fit well.Esto se ve bien, pero la verdad es que este test es muy sensible a qué tan estrictos queremos ser. Al partirlo en 8 es más estricto porque son 8 partes que deben ser similares.
performance_hosmer(mod2, n_bins = 8)## # Hosmer-Lemeshow Goodness-of-Fit Test
##
## Chi-squared: 14.895
## df: 6
## p-value: 0.021## Summary: model does not fit well.Precisión del modelo
Ahora veamos unos estadísticos más interesantes. ¿Qué tan buena puntería tiene nuestro modelo?
performance_accuracy(mod2)## # Accuracy of Model Predictions
##
## Accuracy: 66.89%
## SE: 0.74%-points
## Method: Area under CurveLe atinamos al 67% de nuestras predicciones. ¿Alto o bajo?
Es mejor que el 1
performance_accuracy(mod1)## # Accuracy of Model Predictions
##
## Accuracy: 57.58%
## SE: 0.63%-points
## Method: Area under CurveNormalmente estas métricas nos ayudan a decidirnos qué modelos son mejores que otros. Y lo calibramos al incluir más o menos variables.
La curva de ROC nos ayuda a decidir entre modelos. Queremos la curva que se acerque más a la zona de arriba a la derecha de la gráfica.
plot(performance_roc(mod2, mod1))
Vamos hacerlo un poco más estricto. Vamos a restringir 20% de nuestros datos y que con eso nos diga la predicción. Esta técnica se llama “cross-validation”. Aquí más infor para quien quiera explorarlo. https://www.r-bloggers.com/2019/10/evaluating-model-performance-by-building-cross-validation-from-scratch/
performance_accuracy(mod2, method = "cv",k = 5, verbose = TRUE)## # Accuracy of Model Predictions
##
## Accuracy: 66.86%
## SE: 0.42%-points
## Method: Area under CurveVamos a probar otro método, el de remuestreo o bootstap que vimos el semestre pasado. Usemos 100 submuestras para ver cómo le va.
performance_accuracy(mod2, method = "boot", n= 100, verbose = TRUE)## # Accuracy of Model Predictions
##
## Accuracy: 66.83%
## SE: 0.41%-points
## Method: Area under CurveSeguimos en 2 de 3. Ahora usemos uno todavía más estricto, que ya no es como un volado, sino que parte de la estructura del modelo pero sin predictores.
El resultado nos dice que sí mejoramos con nuestros controles pero que tampoco es dramático el desempeño. Tiene que ser arriba del 50% para que sea mejor que un volado, pero los predictores no nos están ayudando tanto.
performance_pcp(mod2, ci = 0.95, method = "Herron", verbose = TRUE)## # Percentage of Correct Predictions from Logistic Regression Model
##
## Full model: 74.07% [73.58% - 74.56%]
## Null model: 72.56% [72.06% - 73.05%]
##
## # Likelihood-Ratio-Test
##
## Chi-squared: 1598.357
## df: 4.000
## p-value: 0.000En suma, estas métricas nos ayudan a elegir entre varios modelos, por ejemplo, con 5 o con 6 predictores. Pero también nos recuerdan qué tanto podemos confiar en ellos cuando queramos hacer predicciones. Los modelos ayudan a explicar y a predecir. Y para ambas funciones necesitamos estas medidas de calidad.