Regresión logística múltiple

A partir de los datos de la European Social Survey (ESS) tratmos de responder qué variables demográficas (edad, género, haber nacido en el país o no y años de educación completados) explican el sentirse cercano a un partido político. Trataremos de establecer un orden de importancia de los predictores.

Concepto

Nuestro modelo de regresión logística, también llamado modelo logit, es un modelo de respuesta o elección binaria que permite estimar la probabilidad de la variable cualitativa binaria clsprty en función de las covariables (las demás variables independientes). Las observaciones se clasifican en el grupo según se sientan cercanas o no a un partido político dependiendo del valor que tomen las variables empleadas como predictoras.

Carga datos y adecuación variables

Dado que el objetivo de la actividad es la implementación del modelo en sí, no se realiza el debido análisis exploratorio de datos (Junto al procesamiento de datos, son claves para cualquier proyecto de ciencia de datos)

Descripción de variables

clsparty variable binaria sobre el sentimiento de cercanía (1) o no (0) de un partido político

agea variable edad

gndr variable factor género femenino (1) o masculino (0)

brncntr variable factor si ha nacido en el país (1) o no (0)

eduyrs años de educación completado

# Importamos librerias
library(dplyr)
library(foreign)
library(MASS)
library(rcompanion)

# Importamos el dataset
ESS8ES <- read.spss("ESS8ES.sav", reencode='utf-8', to.data.frame = TRUE)


# Seleccion de las variables
data <- ESS8ES %>% 
  dplyr::select(clsprty, agea, gndr, brncntr, eduyrs)


# Recodificacion de las variables
data <- data %>% mutate(clsprty = recode(clsprty,
                                          No = 0,
                                          Yes = 1))

data <- data %>% mutate(gndr = recode(gndr,
                             Male = 0,
                             Female = 1))

data <- data %>% mutate(brncntr = recode(brncntr,
                            Yes = 1,
                            No = 0))


# Asignamos a cada variable las etiqueta correspondiente
data$clsprty <- factor(data$clsprty)
data$gndr <- factor(data$gndr)
data$brncntr <- factor(data$brncntr)
data$eduyrs <- as.numeric(data$eduyrs)
data$agea <- as.numeric(data$agea)

Modelo general

Comenzamos generando un modelo completo introduciendo todas las variables como predictores.

La expresión de la ecuación de regresión logística es la siguiente:

\[\begin{align*} y = Pr(y = 1 | x) = \frac{1}{1 + e^{-(a+b_{1}x_{1}+b_{2}x_{2}+...+b_{j}x_{j})}} = \frac{1}{1 + e^{a + \sum_{j=1}^{p}b_{j}x_{j}}} \end{align*}\]

log_model <- glm(clsprty ~ agea + gndr + 
                        brncntr + eduyrs,
                      data = data,
                      family = 'binomial')
summary(log_model)

## 
## Call:
## glm(formula = clsprty ~ agea + gndr + brncntr + eduyrs, family = "binomial", 
##     data = data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9287  -1.2594   0.8399   1.0342   1.6795  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.580917   0.257116  -6.149 7.81e-10 ***
## agea         0.016885   0.002992   5.644 1.66e-08 ***
## gndr1       -0.123025   0.095898  -1.283 0.199537    
## brncntr1     0.575401   0.154285   3.729 0.000192 ***
## eduyrs       0.067357   0.010399   6.477 9.36e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2531.7  on 1869  degrees of freedom
## Residual deviance: 2460.2  on 1865  degrees of freedom
##   (88 observations deleted due to missingness)
## AIC: 2470.2
## 
## Number of Fisher Scoring iterations: 4

En la salida del modelo, observamos los coeficientes del modelo, la desviación estándar de las observaciones, el estadístico z y los p-valores asociados.

Acorde al modelo, la variable gndr tiene un p-value=0.19. Se mejora el siguiente modelo excluyendo este predictor pues no es estadísticamente significativo y se situaría en la región de rechazo ya que el valor se encuentra por debajo del nivel de significación del 99%.Las demás variables están positivimante relacionadas con la variable binaria.

Descarte variables

Generamos un segundo modelo, descartando la variable gndr

# Modelo descartando variable gndr
log_model_desc <- glm(clsprty ~ agea + 
                        brncntr + eduyrs,
                      data = data,
                      family = 'binomial')
summary(log_model_desc)

## 
## Call:
## glm(formula = clsprty ~ agea + brncntr + eduyrs, family = "binomial", 
##     data = data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9002  -1.2627   0.8425   1.0331   1.6524  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.641555   0.252646  -6.497 8.17e-11 ***
## agea         0.016800   0.002988   5.623 1.87e-08 ***
## brncntr1     0.579942   0.154185   3.761 0.000169 ***
## eduyrs       0.067181   0.010387   6.468 9.93e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2531.7  on 1869  degrees of freedom
## Residual deviance: 2461.9  on 1866  degrees of freedom
##   (88 observations deleted due to missingness)
## AIC: 2469.9
## 
## Number of Fisher Scoring iterations: 4

Una vez generado el segundo modelo, donde las 3 variables son estadísticamente significativas realizamos la evaluación de ambos

Evaluación de modelos

AIC

En la salida de ambos modelos se incluye el criterio de información de Akaike (AIC), cuya regla es escoger los modelos con los valores más bajos.

\[\begin{align*} AIC = -2(ln \space verosimilitud \space - \space número \space de \space parámetros) \end{align*}\]

donde ln verosimilitud es el máximo valor de la función de verosimilitud para el modelo estimado.

El segundo modelo obtiene valores más bajos.

Likelihood

El cociente de verosimilitud (likelihood) usa la diferencia entre la probabilidad de obtener los valores observados con el modelo logístico creado y las probabilidades de hacerlo con un modelo sin relación entre las variables (modelo nulo, sin predictores). Para ello calcula la significancia de la diferencia de residuos entre ambos modelos.

# Ratio de verosimilitud primer modelo
nagelkerke(log_model)$Likelihood.ratio.tes

##  Df.diff LogLik.diff  Chisq    p.value
##       -4     -81.275 162.55 4.1487e-34

# Ratio de verosimilitud segundo modelo
nagelkerke(log_model_desc)$Likelihood.ratio.tes

##  Df.diff LogLik.diff Chisq    p.value
##       -3     -80.452 160.9 1.1698e-34

El segundo modelo mejora significativamente (el contraste p-valor es menor en comparación al primer modelo) la verosimilitud del modelo.

La diferencia (Log-Likelihood), aunque relativamente alta a la del primer modelo, indica que el modelo ha mejorado.

Transformación logarítmica

Para evitar los problemas que supondría modelar nuestras variables mediante un modelo de regresión lineal por mínimos cuadrados \(\beta_{0} + \beta_{1}x\) al obtener valores de Y menores que 0 o mayores que 1, lo que entra en contradicción con el hecho de que las probabilidades logit siempre están dentro del rango [0,1], se usa la función sigmoide. Para voles de x muy grandes positivos, el valor de \(e^{-x}\) es aproximadamente 0 por lo que el valor de la función sigmoide es 1. Para valores de x muy grandes negativos, el valor \(e^{-x}\) tiende a infinito pot lo que valor de la función sigmoide es 0.

\[\begin{align*} Función \space sigmoide = \sigma(x) = \frac{1}{1 + e^{-x}} \end{align*}\]

Se sustituye la x de la ecuación por la función lineal \((\beta_{0} + \beta_{1}x)\), donde Pr(y=1|x) se puede interpretar como: la probabilidad de que la variable cualitativa Y adquiera el valor 1 (se siente cercano en el caso de nuestro modelo) dado los predictores. Por lo que si la probabilidad de que suceda (se siente cercano) es P, la probabilidad de que suceda (no se siente cercano) es igual a 1 menos la probabilidad P.

La función sigmoide puede ajustarse de forma sencilla con métodos de regresión lineal si se emplea su versión logarítimica, obteniendo lo que se conoce como LOG of ODDs. Se define ODDs o la razón de probabilidad de verdadero, como la ratio de probabilidad entre la probabilidad de evento verdadero y la probabilidad de evento falso.

\[\begin{align*} odds = \frac{P}{1 - P} = \frac{Probabilidad \space de \space que \space ocurra \space el \space suceso }{Probabilidad \space de \space que \space no \space ocurra \space el \space suceso} \end{align*}\]

La transformación de probabilidades a ODDs es monotónica, si la probabilidad aumenta también lo hacen los ODDs, y viceversa. El rango de valores que pueden tomar los ODDs es de [0,\(\infty\)]; Dado que el valor de una probabilidad está acotada entre [0,1] se recurre a una transformación logit que consiste en el logaritmo natural de los ODDs. Esto permite convertir el rango de probabilidad previamente limitado [0,1] a \([-\infty, +\infty\)].

Aplicando la transformación logarítmica nos queda la transformación logit que permite identificar el modelo en forma lineal y aditiva:

\[\begin{align*} log(\frac{P}{1 - P} = a + bx) \end{align*}\]

El resultado de los coeficientes de regresión logística expresados mediante el logaritmo de odds de nuestro siguiente modelo:

# Logaritmo de odds del segundo modelo
exp(coef(log_model_desc))

## (Intercept)        agea    brncntr1      eduyrs 
##   0.1936786   1.0169424   1.7859341   1.0694894

Las variables que mayor incremento producen en la ratio de odds es brnctr. Para establecer un orden de importancia de las variables respecto a las demás usamos el z-valor, de esta manera establecemos el orden en el que producen mayor cambio para atribuir la cercanía o no a un partido.

La variable gndr es la que menor incidencia tiene en la cercanía política.

# Ratio de odds y sus intervalos de 
# confianza
exp(cbind(OR = coef(log_model_desc),
          confint.default(log_model_desc)))

##                    OR    2.5 %    97.5 %
## (Intercept) 0.1936786 0.118040 0.3177854
## agea        1.0169424 1.011005 1.0229146
## brncntr1    1.7859341 1.320148 2.4160621
## eduyrs      1.0694894 1.047937 1.0914847

Comparación de las predicciones con las observaciones

Se emplea un threshold de 0.7. Si la probabilidad de sentimiento de cercanía a un partido político es superior a 0.7 se asigna al nivel 1 (sí se siente), si es menor se asigna al nivel 0 (no se siente)

predicciones <- ifelse(test = log_model_desc$fitted.values > 0.7, yes = 1, no = 0)


matriz_confusion <- table(log_model_desc$model$clsprty, predicciones, 
                          dnn=c("observaciones", "predicciones"))

matriz_confusion

##              predicciones
## observaciones   0   1
##             0 708  59
##             1 933 170

fourfoldplot(matriz_confusion, color = c("#CC6666", "#99CC99"),
             conf.level = 0, margin = 1, main = "Matriz de confusión")

Conclusión

El modelo es capaz de clasificar correctamente \(\frac{708+170}{708+170+933+59}=0.469\)(46.9%) de las observaciones de entrenamiento. Si se realiza en detalle cómo se distribuye el error, se aprecia que el modelo solo ha sido capaz de identificar 170 observaciones de 933 de los ciudadanos que se sienten cercanos a un partido político.