La regresión logística es una técnica de estimación estadística que se emplea cuando la variable dependiente es una variable categórica o cualitativa.
Existen dos tipos de modelos de regresión logística:
En esta lección nos concentraremos en las características de la regresión logística binomial.
Como ejemplo para esta lección, utilizaremos un subconjunto de variables de la Encuesta Nacional de Demografía y Salud 2017, realizada por el Instituto Nacional de Estadística e Informática. La base de datos y el libro de códigos pueden descargarse del siguiente enlace: endes17s
load("endes17s.rda")
names(endes17s)
## [1] "CASEID" "edad" "urbano" "tamano" "lengua_indig"
## [6] "year_educ" "nivedu" "relac_jefeh" "jefe_hogar" "exp_medios"
## [11] "nse" "nse_ptje" "seguro" "FACTOREXP" "auto_etnic"
## [16] "auto_etnic2" "region" "hijos_tot" "hijos_vivos" "past_anticon"
## [21] "hoy_anticon" "hoy_anticon2" "condon" "ecivil" "activ_sex"
## [26] "problem_emb"
A partir de los datos correspondientes a las mujeres que han sido activas sexualmente en el mes anterior de la encuesta, se ha calculado la siguiente tabla cruzada:
Tabla 1: Uso de métodos anticonceptivos modernos según área urbana o rural (No. de casos)
## urbano
## hoy_anticon2 Urbana Rural TOTALf
## No 5012 2675 7687
## Sí 9499 3361 12860
## TOTALc 14511 6036 20547
Una probabilidad es una medida que nos indica qué tan posible es que ocurra un evento y se calcula mediante dividiendo las veces de que en que ocurre un evento específico sobre el total de eventos que se observan:
\[\pi = \frac{N_{(x=1)}}{N}\] Las probabilades oscilan entre 0 y 1. De acuerdo con los datos de la tabla 1, la probabilidad de encontrar a una mujer sexualmente activa que use métodos anticonceptivos modernos es de:
\[\pi = \frac{12860}{20547} = 0.623\]
Un odds es otra manera de representar la misma idea que una probabilidad, mide la posibilidad de que ocurra un evento en relación a la posibilidad de que ese evento no ocurra:
\[Odds = \frac{\pi}{1-\pi}\]
En nuestro ejemplo, los Odds para el mismo evento serían:
\[Odds = \frac{12860/20547}{7687/20547} = 1.67\] Ello nos indica que por cada mujer que no usa métodos modernos hay 1.67 mujeres que sí lo hacen. Otra forma de decirlo: es 1.67 veces más probable encontrar a una mujer que use métodos modernos que a una que no lo haga. A diferencia de las probabilidades, los Odds varían entre 0 y \(+\infty\).
Los Odds Ratio (OR) pueden interpretarse como una medida de asociación entre dos variables. Mide la posibilidad de que un evento ocurra dependiendo de la presencia o ausencia de otra condición, por ejemplo: la posibilidad de que se usen métodos anticonceptivos modernos en función del área de residencia urbana:
\[OR = \frac{9499/5012}{3361/2675} = 1.51\] En este ejemplo la posibilidad de usar métodos anticonceptivos modernos respecto de no usarlos son 1.51 veces mayores en zonas urbanas que en zonas rurales.
Como medida de asociación:
Es posible convertir un Odd Ratio en una probabilidad de la siguiente manera:
\[\pi = \frac{OR}{OR + 1}\] Por ejemplo:
\[\pi = \frac{1.51}{1.51 + 1}=0.602\] Un OR = 1.51 equivale a una probabilidad de 0.601 o de 60.2%
Un logaritmo es el exponente o potencia a la cual debe elevarse una base para obtener un número determinado:
\[b^x = n \Leftrightarrow log_b(n) = x\]
El logaritmo Neperiano o natural tiene como base el número e o constante de Neper, cuyo valor aproximado es \(e=2.7182818284590452354...\)
Las transformaciones logarítmicas con útiles si tomamos en cuenta que:
El modelo de regresión logística binomial se representa mediante la siguiente ecuación:
\[log_e(\frac{\pi}{1-{\pi}}) = b_0 + b_1X_1\] En este caso, el coeficiente \(b_1\) mide el efecto de \(X_1\) en el logaritmo natural del Odds de la variable dependiente del evento que estamos analizando.
En nuestro ejemplo, sería el efecto del área de residencia, urbana o rural, en los odds de usar métodos anticonceptivos modernos.
modelo1 <- glm(hoy_anticon2~urbano, binomial,
data = endes17s[endes17s$activ_sex=="Ultimo mes", ])
summary(modelo1)
##
## Call:
## glm(formula = hoy_anticon2 ~ urbano, family = binomial, data = endes17s[endes17s$activ_sex ==
## "Ultimo mes", ])
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4581 -1.2758 0.9206 0.9206 1.0821
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.63935 0.01746 36.62 <2e-16 ***
## urbanoRural -0.41106 0.03124 -13.16 <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: 27168 on 20546 degrees of freedom
## Residual deviance: 26996 on 20545 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 27000
##
## Number of Fisher Scoring iterations: 4
El modelo resultante para el uso de métodos anticonceptivos modernos entre mujeres en edad fértil que son activas sexualmente se puede representar de la siguiente manera:
\[log_e(\frac{\pi}{1-\pi}) = 0.64 - (0.41*Urbana)\] Como se aprecia en los resultados, en el caso de la variable “Urbana”, la categoría de control es “Zona urbana”, por lo que el coficiente \(b_1\) mide el efecto de vivir en una zona rural en el uso de métodos anticonceptivos modernos. En este caso, el coeficiente \(b_0\) sería el logaritmo natural de las posibilidades de uso de anticonceptivos modernos en zonas urbanas.
Dado que: \(\pi = \frac{OR}{OR + 1}\)
Podemos decir que: \(\pi = \frac{e^{b_0+b_1X_1}}{{e^{b_0+b_1X_1}}+1}\)
Por lo tanto en zonas urbanas: \(\pi = \frac{e^{0.64}}{{e^{0.64}}+1}= 0.655\)
Y en zonas rurales: \(\pi = \frac{e^{0.64-0.41}}{{e^{0.64-0.41}}+1}= 0.557\)
como se aprecia, la ecuación de regresión logística en nuestro ejemplo modela los resultados de la tabla de contingencia que vimos al inicio de la presentación.
## Urbana Rural TOTAL
## No 34.5 44.3 37.4
## Sí 65.5 55.7 62.6
## TOTAL 100.0 100.0 100.0
## Nvalid 14511.0 6036.0 20547.0
\[log_e(\frac{\pi}{1-\pi}) = 0.64 - (0.41*Urbana)\]
La desvianza puede interpretarse de manera análoga a la Suma Total de Cuadrados de una regresión lineal. En el caso de la regresión logística, se reportan dos desvianzas, una para un modelo “nulo”, que incluye únicamente la intercepción y ninguna variable independiente; otra que es la desvianza residual, es decir la que queda luego de haber incluido las variables independientes del modelo. Cuanto mayor sea la diferencia entre ambas desvianzas, mejor será el ajuste del modelo.
modelo1 <- glm(hoy_anticon2~urbano, binomial,
data = endes17s[endes17s$activ_sex=="Ultimo mes", ])
summary(modelo1)
##
## Call:
## glm(formula = hoy_anticon2 ~ urbano, family = binomial, data = endes17s[endes17s$activ_sex ==
## "Ultimo mes", ])
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4581 -1.2758 0.9206 0.9206 1.0821
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.63935 0.01746 36.62 <2e-16 ***
## urbanoRural -0.41106 0.03124 -13.16 <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: 27168 on 20546 degrees of freedom
## Residual deviance: 26996 on 20545 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 27000
##
## Number of Fisher Scoring iterations: 4
Es posible usar una prueba de Chi Cuadrado para evaluar si la diferencia entre las desvianzas es estadísticamente significativa:
anova(modelo1, test = "Chisq")
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: hoy_anticon2
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 20546 27168
## urbano 1 172.08 20545 26996 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Otra medida que indica el “ajuste” de un modelo de regresión logística binomial es el “Pseudo \(R^2\)”. Existen diversos tipos de Pseudo \(R^2\), entre los más conocidos está el Cox-Snell \(R^2\) y el Nagelkerke \(R^2\). El R tiene diversas funciones que permiten calcularlos:
library(DescTools)
PseudoR2(modelo1, c("CoxSnell", "Nagel"))
## CoxSnell Nagelkerke
## 0.008339999 0.011370799
En principio, cuanto mayor es el Pseudo \(R^2\) mejor explica el modelo a la variable dependiente.
##
## Call:
## glm(formula = hoy_anticon2 ~ urbano + year_educ + nse, family = binomial,
## data = endes17s[endes17s$activ_sex == "Ultimo mes", ])
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5859 -1.3448 0.8911 0.9420 1.2154
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.050587 0.056060 0.902 0.36686
## urbanoRural -0.139618 0.044692 -3.124 0.00178 **
## year_educ 0.034039 0.004445 7.658 1.89e-14 ***
## nseMedio bajo 0.252854 0.047961 5.272 1.35e-07 ***
## nseMedio 0.293689 0.057343 5.122 3.03e-07 ***
## nseMedio alto 0.160673 0.062549 2.569 0.01021 *
## nseAlto 0.162417 0.069825 2.326 0.02002 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 27168 on 20546 degrees of freedom
## Residual deviance: 26867 on 20540 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 26881
##
## Number of Fisher Scoring iterations: 4
Con la siguiente sintaxis podemos calcular las probabilidades estimadas por nuestro modelo de que mujeres con diferentes características usen métodos anticonceptivos modernos:
new.data <- data.frame(urbano = c("Urbana", "Rural"), year_educ = 10, nse = "Bajo")
predict(modelo2, new.data, type = "response")
## 1 2
## 0.5965173 0.5625104
new.data <- data.frame(urbano = "Urbana", year_educ = c(10, 17), nse = "Bajo")
predict(modelo2, new.data, type = "response")
## 1 2
## 0.5965173 0.6523185
new.data <- data.frame(urbano = "Urbana", year_educ = 10, nse = c("Bajo", "Alto"))
predict(modelo2, new.data, type = "response")
## 1 2
## 0.5965173 0.6349221
La prueba de Chi Cuadrado evalúa si hay una reducción significativa de la desvianza por el hecho de incluir una variable adicional en el modelo de regresión.
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: hoy_anticon2
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 20546 27168
## urbano 1 172.081 20545 26996 < 2.2e-16 ***
## year_educ 1 88.698 20544 26907 < 2.2e-16 ***
## nse 4 39.547 20540 26867 5.371e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Se puede también hacer una prueba de Chi Cuadrado de la desvianza para evaluar el efecto de quitar una variable
drop1(modelo2, test = "Chisq")
## Single term deletions
##
## Model:
## hoy_anticon2 ~ urbano + year_educ + nse
## Df Deviance AIC LRT Pr(>Chi)
## <none> 26867 26881
## urbano 1 26877 26889 9.732 0.00181 **
## year_educ 1 26926 26938 58.626 1.906e-14 ***
## nse 4 26907 26913 39.547 5.371e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
PseudoR2(modelo2, c("CoxSnell", "Nagel"))
## CoxSnell Nagelkerke
## 0.01451020 0.01978328
tab_cruzada <- function(tabla, pc="s"){
tab.0 <- tabla
if(pc=="s"){
tab.1 <- addmargins(tab.0, 1, FUN = list(list(TOTALc = sum)))
tab.1 <- addmargins(tab.1, 2, FUN = list(list(TOTALf = sum)))
return(tab.1)}
if (pc=="col"){
tab.2 <- addmargins(tab.0, 2, FUN = list(list(TOTAL = sum)))
tab.2 <- round(prop.table(tab.2, 2)*100,2)
tab.2 <- round(addmargins(tab.2, 1, FUN = list(list(TOTAL = sum))),1)
tab.0 <- tab.0
tab.0 <- addmargins(tab.0, 2, FUN = list(list(TOTAL = sum)))
tab.0 <- addmargins(tab.0, 1, FUN = list(list(Nvalid = sum)))
tab.2 <- rbind(tab.2, tail(tab.0, 1))
return(tab.2)}
if (pc=="fila"){
tab.2 <- addmargins(tab.0, 1, FUN = list(list(TOTAL = sum)))
tab.2 <- round(prop.table(tab.2, 1)*100,2)
tab.2 <- round(addmargins(tab.2, 2, FUN = list(list(TOTAL = sum))),1)
tab.0 <- tab.0
tab.0.1 <- addmargins(tab.0, 1, FUN = list(list(TOTAL = sum)))
tab.0.1 <- addmargins(tab.0.1, 2, FUN = list(list(Nvalid = sum)))
tab.2 <- cbind(tab.2, Nvalid=tab.0.1[, "Nvalid"])
return(tab.2)}
}
load("endes17s.rda")
names(endes17s)
# Función para tabulación cruzada
tabla1 <- xtabs(~hoy_anticon2+urbano,
data = endes17s[endes17s$activ_sex=="Ultimo mes", ])
tab_cruzada(tabla1, pc="s") # Tabla valores absolutos
tab_cruzada(tabla1, pc="col") # Tabla en % verticales
# Modelo1
modelo1 <- glm(hoy_anticon2~urbano, binomial,
data = endes17s[endes17s$activ_sex=="Ultimo mes", ])
summary(modelo1)
anova(modelo1, test = "Chisq")
library(DescTools)
PseudoR2(modelo1, c("CoxSnell", "Nagel"))
#Modelo 2
modelo2 <-glm(hoy_anticon2~urbano + year_educ + nse, binomial,
data = endes17s[endes17s$activ_sex=="Ultimo mes", ])
summary(modelo2)
anova(modelo2, test = "Chisq")
drop1(modelo2, test = "Chisq")
PseudoR2(modelo2, c("CoxSnell", "Nagel"))