Regresión Logística

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.

Ejemplo y Base de Datos

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"

Ejemplo: Uso de anticonceptivos modernos

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

Probabilidades

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\]

Odds o posibilidades

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\).

Odds Ratio o Ventaja Relativa

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.

OR como medida de asociación

Como medida de asociación:

Odds ratio y probabilidades

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%

Logaritmo natural o Neperiano

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...\)

Transformaciones logarítmicas y regresión

Las transformaciones logarítmicas con útiles si tomamos en cuenta que:

Modelo de Regresión Logística Binomial

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.

Cálculo de los coeficientes del modelo logístico binomial

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

Interpretación de los resultados

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.

Interpretación de coeficientes

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.

Modelo de Regresión Logística

##         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)\]

Desvianza y ajuste del modelo (1)

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.

Cálculo de los coeficientes del modelo logístico binomial

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

Desvianza y ajuste del modelo (2)

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

Pseudo \(R^2\)

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.

Regresión logística Múltiple

## 
## 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

Cálculo de probabilidades

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

Prueba de significancia del modelo

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

Prueba de Significancia del modelo - hacia atrás

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

Pseudo \(R^2\) para el modelo 2

PseudoR2(modelo2, c("CoxSnell", "Nagel"))
##   CoxSnell Nagelkerke 
## 0.01451020 0.01978328

Sintaxis en R para esta lección: Función de tabulación cruzada

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)}
}

Carga de datos y tablas de contingencia

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

Sintaxis modelos logísticos

# 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"))