Este texto se basa en los siguientes materiales:

> library(tidyverse)
> library(broom)
> library(car)

Introducción

En este capítulo presentamos la regresión logística como herramienta de modelado en aquellas situaciones en las que hay una variable de respuesta categórica con dos niveles, por ejemplo, sí y no. La regresión logística es un tipo de modelo lineal generalizado (GLM) para variables de respuesta donde la regresión múltiple regular no funciona muy bien. Los GLM se pueden considerar como un enfoque de modelado de dos etapas. Primero modelamos la variable de respuesta utilizando una distribución de probabilidad, como la distribución binomial o de Poisson. Segundo, modelamos el parámetro de la distribución usando una colección de predictores y una forma especial de regresión múltiple. En última instancia, la aplicación de un GLM se sentirá muy similar a la regresión múltiple, incluso si algunos de los detalles son diferentes.

El problema: autopercepción de clase

Vamos, como siempre, a trabajar con la ENES que contiene una pregunta sobre autoposicionamiento de clases. Sin embargo, la misma solamente se realizó a los PSH (principales sostenes económicos del hogar) y a les cónyuges de los mismos. Por ello, será necesario trabajar no con toda la población de personas sino solamente con estos últimos.

> psh_cony <- read_rds("./data/ENES_psh_cony.rds")

Pero primero, tenemos que recodificar las variables de percecpión de clase MS_v260 y MS_v261 a una sola variable que vamos a llamar percep_clase y luego a esta, dictomizarla. A su vez, vamos a tratar de replicar el esquema que utilizó Rodolfo Elbert en este artículo. Allí, realizó una desagregación de los trabajadores por cuenta propia en dos categorías:

  • informales: aquellos TCP que trabajan en tareas no calificadas y operativas
  • pequeña burguesía formal: aquellos TCP que trabajan en tareas profesionales y técnicas

Dado que de esta forma, quedan muy pocos casos en la cateogría “Pequeña burguesía”, la vamos a unificar con la categoría de “Empleadores” bajo el rótulo: “Clases propietarias”.

> psh_cony <- psh_cony %>%
+   mutate(
+     v109 = if_else(v109 == 'Varón', 'Masculino', 'No masculino'),
+     percep_clase = case_when(
+            !is.na(MS_v260) ~ MS_v260,
+            is.na(MS_v260) ~ MS_v261
+            )
+     )
> 
> psh_cony <- psh_cony %>%
+   mutate(
+     salud = case_when(
+       v134a %in% c("Obra social", "Prepaga", "Prepaga a través de Obra Social") ~ 'Tiene cobertura',
+       TRUE ~  'No tiene cobertura'
+       ),
+     
+     percep_obrera = as.factor(case_when(
+             percep_clase == 'Clase baja' | percep_clase == 'Clase obrera' ~ 'Obrera',
+             is.na(percep_clase) ~ NA_character_,
+             TRUE ~ 'No obrera')),
+     
+     percep_obrera_orig = as.factor(case_when(
+             MS_v260 == 'Clase baja' | MS_v260 == 'Clase obrera' ~ 'Obrera',
+             is.na(MS_v260) ~ NA_character_,
+             TRUE ~ 'No obrera')),
+     CNO_calif_agg = case_when(
+            CNO_calif == 'Profesional' | CNO_calif == 'Técnica' ~ 'Prof./técn',
+            CNO_calif == 'Operativa' | CNO_calif == 'No calificada' ~ 'Op./no calif.',
+          ),
+     nivel_ed_agg = case_when(
+                   nivel_ed == 'Menores de 5 años' | 
+                   nivel_ed == 'Sin instrucción (incluye nunca asistió o sólo asistió a sala de 5)' |
+                   nivel_ed == 'Primaria/EGB incompleto' | 
+                   nivel_ed == 'Primaria/EGB completo' |
+                   nivel_ed == 'Educación especial' | nivel_ed == 'NS/NR'~ '0_Bajo',
+      
+                 nivel_ed == 'Secundario/Polimodal incompleto' | 
+                 nivel_ed == 'Secundario/Polimodal completo' ~ '1_Medio',
+     
+                 nivel_ed == 'Terciario incompleto' | 
+                 nivel_ed == 'Terciario completo' |
+                 nivel_ed == 'Universitario incompleto' | 
+                 nivel_ed == 'Universitario completo' ~ '2_Alto'),
+          
+          class_eow_agg = case_when(
+          EOW_class == 'Managers' | EOW_class == 'Supervisores' ~ 'Manag./superv.',
+          EOW_class == 'Trabajadores' ~ 'Trabajadores',
+          EOW_class == 'Pequeña burguesía' & CNO_calif_agg == 'Prof./técn' ~ 'Propietarias',
+          EOW_class == 'Pequeña burguesía' & CNO_calif_agg == 'Op./no calif.' ~ 'TCP informales',
+          
+          EOW_class == 'Empleadores' ~ 'Propietarias'),
+    )

Modelando la probabilidad de un evento

La regresión logística es un modelo lineal generalizado donde el resultado es una variable categórica de dos niveles. El resultado, \(Y_{i}\), toma el valor 1 (en nuestra aplicación, esto representa percepción de clase obrera) con probabilidad \(p_{i}\) y valor 0 con probabilidad \(1-p_{i}\). Debido a que cada observación tiene un “contexto” ligeramente diferente, por ejemplo, posición objetiva de clase diferente, niveles educativos diferentes, etc., la probabilidad \(p_{i}\) será diferente para cada observación. En última instancia, es esta probabilidad la que modelamos en relación con las variables predictoras: examinaremos qué características de la persona se asocian a una percepción obrera.

Notación para un modelo de regresión logística

La variable de resultado para un GLM (acrónimo de “generalized linera model”) se denota por \(Y_{i}\), donde el índice \(i\) se utiliza para representar la i-ésima observación. En la ENES, \(Y_{i}\) se utilizará para representar si la persona \(i\) tiene percepción de clase obrera (\(Y_{i} = 1\)) o no (\(Y_{i} = 0\)).

Las variables predictoras son representadas de la siguiente manera: \(X_{1,i}\) es el valor de la variable 1 en el caso \(i\), \(X_{2,i}\) es el valor de la variable 2 en el caso \(i\), y así sucesivamente.

\[transformacion(p_{i}) = \beta_{0} + \beta_{1}X_{1,i} + \beta_{2}X_{2,i} + ... + + \beta_{k}X_{k,i}\]

Queremos elegir una transformación en la ecuación que tenga sentido práctico y matemático. Por ejemplo, queremos una transformación que haga que el rango de posibilidades del lado izquierdo de la ecuación sea igual al rango de posibilidades del lado derecho. Si no hubiera transformación para esta ecuación, el lado izquierdo solo podría tomar valores entre 0 y 1, pero el lado derecho podría tomar valores fuera de este rango.

Por ejemplo, en nuestro caso, sabemos que la variable percep_obrera asume solamente dos valores. ¿Cómo se distribuyen los casos según los ingresos individuales?

> psh_cony %>%
+     drop_na(percep_obrera) %>%
+     mutate(percep_obrera = if_else(percep_obrera == "Obrera", 1, 0)) %>%
+     ggplot(aes(y = percep_obrera, x = ITI)) + geom_smooth(method = "lm", se = FALSE) +
+     geom_point() + ylim(0, 1) + theme_minimal() + labs(x = "Ingreso total individual",
+     y = "Percep. obrera")

En azul vemos la recta de regresión lineal tomando como variable dependiente la autopercepción y como independiente los ingresos. ¿Qué problemas se podrían presentar?

Una transformación común para evitar que los valores estimados caigan fuera del rango 0-1 (rango admimitdo por una probabilidad) es calcular el logit que puede ser escrito de la siguiente manera:

\[logit(p_{i}) = log_e\left(\frac{p_{i}}{1-p_{i}}\right)\] Los valores de la transformación logit podemos verlos en el gráfico a continuación.

> logit <- function(x) {
+     log(x/(1 - x))
+ }
> 
> 
> logit(seq(0.01, 0.99, 0.01)) %>%
+     as_tibble() %>%
+     ggplot(aes(y = seq(0.01, 0.99, 0.01), x = value)) + geom_line() + geom_hline(yintercept = 1,
+     linetype = "dashed") + geom_hline(yintercept = 0, linetype = "dashed") + ylim(-0.01,
+     1.01) + theme_minimal() + labs(x = "Logit", y = "Prob")

En nuestro ejemplo, vamos a usar unas 8 variables predictoras, por lo cual \(k=8\). Si bien la elección precisa de una función logit no es intuitiva, se basa en la teoría que sustenta los modelos lineales generalizados, que está más allá del alcance de este curso. Afortunadamente, una vez que ajustamos un modelo usando software, empezaremos a sentir que estamos de vuelta en el contexto de la regresión múltiple, incluso si la interpretación de los coeficientes es un poco más compleja.

A continuación escribimos nuevamente la ecuación que relaciona a los predictores con \(Y_{i}\) usando la transformación logit sobre \(p_{i}\):

\[log_e\left(\frac{p_{i}}{1-p_{i}}\right) = \beta_{0} + \beta_{1}X_{1,i} + \beta_{2}X_{2,i} + ... + + \beta_{k}X_{k,i}\]

Para convertir valores en la escala de regresión logística a la escala de probabilidad, necesitamos transformar hacia atrás y luego resolver para \(p_{i}\):

\[log_e\left(\frac{p_{i}}{1-p_{i}}\right) = \beta_{0} + \beta_{1}X_{1,i} + \beta_{2}X_{2,i} + ... + + \beta_{k}X_{k,i} \\ \frac{p_{i}}{1-p_{i}} = e^{ \beta_{0} + \beta_{1}X_{1,i} + \beta_{2}X_{2,i} + ... + + \beta_{k}X_{k,i}} \\ p_{i} = (1-p_{i}) e^{ \beta_{0} + \beta_{1}X_{1,i} + \beta_{2}X_{2,i} + ... + + \beta_{k}X_{k,i}} \\ p_{i} = e^{ \beta_{0} + \beta_{1}X_{1,i} + \beta_{2}X_{2,i} + ... + \beta_{k}X_{k,i}} - p_{i} \times e^{ \beta_{0} + \beta_{1}X_{1,i} + \beta_{2}X_{2,i} + ... + + \beta_{k}X_{k,i}} \\ p_{i} + p_{i} \times e^{ \beta_{0} + \beta_{1}X_{1,i} + \beta_{2}X_{2,i} + ... + \beta_{k}X_{k,i}} = e^{ \beta_{0} + \beta_{1}X_{1,i} + \beta_{2}X_{2,i} + ... + \beta_{k}X_{k,i}} \\ p_{i} (1 + e^{ \beta_{0} + \beta_{1}X_{1,i} + \beta_{2}X_{2,i} + ... + \beta_{k}X_{k,i}}) = e^{ \beta_{0} + \beta_{1}X_{1,i} + \beta_{2}X_{2,i} + ... + \beta_{k}X_{k,i}} \\ p_{i} = \frac{e^{ \beta_{0} + \beta_{1}X_{1,i} + \beta_{2}X_{2,i} + ... + \beta_{k}X_{k,i}}}{ 1 + e^{ \beta_{0} + \beta_{1}X_{1,i} + \beta_{2}X_{2,i} + ... + \beta_{k}X_{k,i}}}\]

Al igual que con la mayoría de los problemas de datos aplicados, sustituimos los parámetros (\(\beta_{k}\)) por las estimaciones puntuales (\(\hat{\beta_{k}}\)).

Un ejemplo mínimo para empezar

Empecemos con un ejemplo simple. Queremos ver si entre las diferentes posiciones en el hogar (cónyuge y PSH) existen diferentes percepciones de clase. Así, podemos entrenar nuestra regresión logística. Para eso, usamos la función glm() que, como verán, tiene la misma sintaxis que lm().

> glm_1 <- psh_cony %>%
+     drop_na(percep_obrera) %>%
+     glm(percep_obrera ~ v111, family = "binomial", data = .)
> 
> tidy(glm_1)
## # A tibble: 2 × 5
##   term        estimate std.error statistic   p.value
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)  -0.656     0.0237    -27.7  2.30e-168
## 2 v111Cónyuge  -0.0835    0.0393     -2.13 3.35e-  2

De esta forma, entrenamos una regresión logística para la autopercepción de clase obrera y, como variable regresora, la posicion en el hogar: v111 toma valor 1 cuando la persona es cónyuge y 0 cuando la persona es el PSH. La ecuación del modelo es:

\[log_e\left(\frac{P(obrera)_{i}}{1-P(obrera)_{i}}\right) = -0.656 - 0.084 \times v111_{cony, i}\]

Si se considera una persona (de les cónyuges o PSH) elegida al azar y se trata del PSH, entonces v111 toma valor 0 y el lado derecho de la ecuación del modelo es igual a -0.249. Resolviendo para \(p(obrera)_{i}\):

\[\frac{e^{-0.656}}{1 + e^{-0.656}} = 0.3416\]

Es decir que la probabilidad estimada de autopercepción obrera es de \(\hat{p}(obrera)_{i}=0.3416\). Así como etiquetamos un valor estimado de \(y_{i}\) con un “sombrero” en regresión de una variable y múltiple, hacemos lo mismo para esta probabilidad

En cambio, si se trata de una persona que es la cónyuge, v111 toma valor 1. Así, \[p(obrera)_{i} = \frac{e^{- 0.656 - 0.084 \times 1}}{1 + e^{-0.656 - 0.084 \times 1}} = 0.3230\] Si bien sabemos que existe cierta diferencia (aunque muy pequeña) en la percepción de clase, nos gustaría tener en cuenta muchas variables diferentes a la vez para comprender cómo cada una de las diferentes características de la persona se vinculan a la percepción obrera.

Ahora bien… podemos hacer esta cuenta volviendo a utilizar la función predict. Solamente tenemos que generar un dataframe o un tibble que contenga los predictores y sus valores:

> tibble(v111 = "Cónyuge") %>%
+     predict(glm_1, ., type = "response")
##         1 
## 0.3231066

Modelo logístico con muchas variables

Ajustemos el modelo de regresión logística con 5 predictores. Al igual que la regresión múltiple, el resultado se puede presentar en una tabla de resumen, que se muestra en la tabla siguiente.

> glm_2 <- psh_cony %>%
+     filter(estado == "Ocupado") %>%
+     glm(percep_obrera ~ class_eow_agg + v109 + v111 + v108 + nivel_ed_agg, data = .,
+         family = "binomial")
> 
> tidy(glm_2)
## # A tibble: 9 × 5
##   term                        estimate std.error statistic   p.value
##   <chr>                          <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)                  0.0649    0.126       0.513 6.08e-  1
## 2 class_eow_aggPropietarias   -0.178     0.114      -1.56  1.20e-  1
## 3 class_eow_aggTCP informales  0.564     0.0900      6.26  3.75e- 10
## 4 class_eow_aggTrabajadores    0.419     0.0775      5.41  6.36e-  8
## 5 v109No masculino             0.163     0.0568      2.88  4.03e-  3
## 6 v111Cónyuge                 -0.256     0.0625     -4.09  4.31e-  5
## 7 v108                        -0.00997   0.00199    -5.02  5.16e-  7
## 8 nivel_ed_agg1_Medio         -0.608     0.0568    -10.7   9.95e- 27
## 9 nivel_ed_agg2_Alto          -1.51      0.0708    -21.3   1.12e-100

Al igual que con la regresión múltiple, podríamos recortar algunas variables del modelo. Aquí usaremos una estadística llamada criterio de información de Akaike (AIC), que es análoga al \(R^2\) en regresión múltiple. AIC es un método popular de selección de modelos utilizado en muchas disciplinas y es elogiado por su énfasis en la incertidumbre y parsimonia del modelo. AIC selecciona un “mejor” modelo clasificando los modelos de mejor a peor según sus valores de AIC. En el cálculo del AIC de un modelo, se aplica una penalización por incluir variables adicionales. Esta penalización por la complejidad añadida del modelo intenta lograr un equilibrio entre el ajuste insuficiente (muy pocas variables en el modelo) y el ajuste excesivo (demasiadas variables en el modelo). Cuando se usa AIC para la selección de modelos, los modelos con un valor de AIC más bajo se consideran “mejores”. Recuerde que al usar \(R^2\) ajustado en su lugar, seleccionamos modelos con valores más altos. Es importante tener en cuenta que AIC proporciona información sobre la calidad de un modelo en relación con otros modelos, pero no proporciona información sobre la calidad general de un modelo.

En este caso, posiblemente una primera varible candidata a eliminar podría ser la variable asociada al género. No obstante, veremos cuestiones vinculadas a la selección de modelos dentro de algunas clases.

La variable v109 tiene solamente dos valores (masculino y no masculino)… basándose en los resultados del modelo, ¿qué nos dice esa variable acerca de la autopercepción de clase?

El coeficiente mostrado corresponde al nivel de no masculino es positivo El mismo refleja una ganancia en la proporción de autopercepción obrera para aquellas personas que no son masculinas. El modelo sugiere que las personas no masculinas tienden a tener una percepción de clase obrera mayor que las personas masculinas.

Si se fijan en el coeficiente asociado a la variable v111 ha cambiado (aunque ha mantenido el signo). En la mayoría de los datos provenientes de estudios observacionales -como una encuesta- y no experimentales, es común que las estimaciones puntuales cambien un poco y, a veces, mucho, dependiendo de qué otras variables se incluyan en el modelo.

Actividad 1

Use el modelo anterior para estimar la probabilidad de tener autopercepción de clase obrera que presenta una persona masculina, principal sostén del hogar, de 38 años de edad, de nivel educativo medio y cuya posición de clase objetiva es pequeña burguesía formal (clases propietarias).

Podemos empezar escribiendo la ecuación completa del modelo:

\[log_e\left(\frac{P(obrera)_{i}}{1-P(obrera)_{i}}\right) = \\ 0.064888 - 0.177762 \times class\_eow\_agg_{prop} + \\ 0.563875 \times class\_eow\_agg_{tcp\ informal} + \\ 0.419077 \times class\_eow\_agg_{trabajador} - \\ 0.255518 \times v111_{conyugye} - \\ 0.163277 \times v109{no\ masculino} - \\ 0.009968 \times v108 - \\ 0.608011 \times nivel\_ed\_agg_{medio} - \\ 1.507580 \times nivel\_ed\_agg_{alto} \]

Ahora, podemos reemplazar los valores correspondentes a cada variable

\[log_e\left(\frac{P(obrera)_{i}}{1-P(obrera)_{i}}\right) = \\ 0.064888 - 0.177762 \times 1 + \\ 0.563875 \times 0 + \\ 0.419077 \times 0 - \\ 0.255518 \times 0 - \\ 0.163277 \times 0 - \\ 0.009968 \times 38 - \\ 0.608011 \times 1 - \\ 1.507580 \times 0 = -1.099677 \] Ahora podemos resolver para \(\hat{p}\) para obtener las probabilidades de este perfil de presentar una autopercepción de clase obrera.

\[\frac{e^{-1.099677}}{1 + e^{-1.099677}} = 0.2498\]

Es decir, un 25% de chances de presentar autopercepción obrera.

Hagámoslo con predict(). Es importante pasar el mismo nombre de las variables predictoras y escribir las categorías de la misma forma en que entraron al modelo. Si no va a arrojar un error:

> tibble(v109 = "Masculino", v111 = "PSH", v108 = 38, nivel_ed_agg = "1_Medio", class_eow_agg = "Propietarias") %>%
+     predict(glm_2, ., type = "response")
##         1 
## 0.2498008