Capítulo 9 Regresión logística

9.1 Resumen

En este capítulo:

Estimamos e interpretamos modelos de regresión logística para variables dependientes categóricas binarias.

Principales conceptos: regresión logística, modelo de probabilidad lineal. Funciones clave: glm().

9.1.1 Librerías

Vamos a utilizar las siguientes librerías:

library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --

## v ggplot2 3.3.3     v purrr   0.3.4
## v tibble  3.1.2     v dplyr   1.0.6
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   1.4.0     v forcats 0.5.1

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(ggeffects) # efectos en modelos de regresion
library(modelsummary) # tablas de regresion
library(broom) # organizar resultados de modelos
library(haven) # leer archivos .dta

theme_set(theme_light())

9.1.2 Datos

Como ejemplo, vamos a cargar directamente de la web una base de datos con información sobre la admisión a programas de posgrado universitario en una universidad de Estados Unidos:

dat <- read_csv("https://stats.idre.ucla.edu/stat/data/binary.csv")

## 
## -- Column specification --------------------------------------------------------
## cols(
##   admit = col_double(),
##   gre = col_double(),
##   gpa = col_double(),
##   rank = col_double()
## )

Al final, para un taller, debemos descargar el siguiente archivos de datos y guardarlo en la carpeta /data de nuestro proyecto de RStudio:

Datos de votos por el NAFTA: link. Para descargar, hacer click en “Download” o “View raw”.

9.2 Variables dependientes categóricas

No todos los fenómenos que estudiamos se pueden expresar como variables numéricas, ordenadas de menor a mayor y en las que cada intervalo (cada aumento de una unidad) es constante. En ocasiones, nos interesa estudiar cualidades o categorías, como el tipo de régimen, la adopción de una política pública o la presencia de guerra civil en un país.

El modelo de regresión lineal por mínimos cuadrados ordinarios (OLS) no es ideal para modelar variables dependientes categóricas. Sin embargo, sí lo podemos utilizar, una aplicación que llamamos el modelo de probabilidad lineal (MPL). Típicamente, una mejor opción es utilizar otro tipo de regresión conocido como regresión logística o el modelo logit. En este instructivo, veremos cómo podemos usar ambos.

Empecemos por mirar los datos sobre admisiones a posgrados:

dat

## # A tibble: 400 x 4
##    admit   gre   gpa  rank
##    <dbl> <dbl> <dbl> <dbl>
##  1     0   380  3.61     3
##  2     1   660  3.67     3
##  3     1   800  4        1
##  4     1   640  3.19     4
##  5     0   520  2.93     4
##  6     1   760  3        2
##  7     1   560  2.98     1
##  8     0   400  3.08     2
##  9     1   540  3.39     3
## 10     0   700  3.92     2
## # ... with 390 more rows

Tenemos 400 observaciones, que corresponden a 400 candidatos a posgrado, y cuatro variables:

admit: dummy de admisión (0=no admitido; 1=admitido).
gre: nota en el examen GRE del candidato. Puede tomar valores de 200 a 800.
gpa: promedio crédito acumulado de pregrado del candidato. Puede tomar valores de 0 a 4.
rank: ránquin de la universidad de pregrado del candidato. Hay 4 categorías, de mayor (1) a menor (4) calidad.

Estamos interesados en entender los efectos de la historia académica de un candidato sobre la probabilidad de que un individuo sea admitido a la universidad (admit = 1). Convertimos dos variables categóricas -admit y rank- a tipo factor, porque pese a que conceptualmente son categóricas, están codificadas como numéricas en R:

dat <- dat %>%
  mutate(
    admit_fct = factor(admit, levels = c(0, 1), labels = c("No", "Sí")),
    rank = factor(rank)
  )

Es importante primero hacer un análisis exploratorio de nuestros datos. Como ya recodificamos y limpiamos un poco los datos, veamos la distribución de la variable dependiente con una gráfica de frecuencias:

dat %>%
  ggplot(aes(x = admit_fct)) +
  geom_bar() +
  labs(x = "Admitido a posgrado",
       y = "Número de candidatos")

O quizás nos interese ver si hay diferencias entre candidatos que vienen de distintos tipos de universidades:

dat %>%
  ggplot(aes(x = admit_fct, fill = rank)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent) +
  scale_fill_brewer() +
  labs(x = "Admitido a posgrado",
       y = "Porcentaje de candidatos",
       fill = "Ránquin universidad \nde pregrado")

Y veamos la relaciones entre esta variable y las dos variables numéricas en la base de datos:

dat %>%
  group_by(admit_fct) %>%
  summarize(
    obs = n(),
    media_gre = mean(gre),
    media_gpa = mean(gpa)
  )

## # A tibble: 2 x 4
##   admit_fct   obs media_gre media_gpa
##   <fct>     <int>     <dbl>     <dbl>
## 1 No          273      573.      3.34
## 2 Sí          127      619.      3.49

Hasta ahora, todo bien… Pero veamos una de estas relaciones como una gráfica de dispersión con una línea de tendencia (¡estimada automáticamente por ggplot2 via OLS!):

dat %>%
  ggplot(aes(x = gre, y = admit)) +
  geom_point(alpha = 0.25, size = 2) +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  scale_y_continuous(breaks = c(0, 1), labels = c("No", "Sí")) +
  labs(x = "Puntaje GRE", y = "Admisión")

## `geom_smooth()` using formula 'y ~ x'

Algo no está bien… La línea de tendencia de geom_smooth() es estimada usando OLS - el modelo hace lo posible por ajustarse a los datos, pero tratar de conectar las observaciones no admitidas y sus puntajes con las de estudiantes admitidos es difícil por la naturaleza discreta de la variable dependiente. Exploremos las alternativas.

9.3 El model de probabilidad lineal (MPL)

Recordemos: nos interesa ver la probabilidad de que un estudiante sea admitido, dadas sus notas en el examen GRE, su promedio crédito acumulado (GPA) y el ránquin de su universidad de pregrado (1-4).

Podemos hacer una primera aproximación estimando un modelo de probabilidad lineal (MPL). Usamos lm() como si estuviéramos estimando cualquier otro modelo de regresión lineal - ¡porque eso es lo que estamos haciendo! Noten que la variable dependiente en la fórmula es la versión numérica de la misma (admit) y no el factor que creamos arriba (admit_fct):

mpl <- lm(admit ~ gre + gpa + rank,
          data = dat)
summary(mpl)

## 
## Call:
## lm(formula = admit ~ gre + gpa + rank, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7022 -0.3288 -0.1922  0.4952  0.9093 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.2589102  0.2159904  -1.199   0.2314    
## gre          0.0004296  0.0002107   2.038   0.0422 *  
## gpa          0.1555350  0.0639618   2.432   0.0155 *  
## rank2       -0.1623653  0.0677145  -2.398   0.0170 *  
## rank3       -0.2905705  0.0702453  -4.137 4.31e-05 ***
## rank4       -0.3230264  0.0793164  -4.073 5.62e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4449 on 394 degrees of freedom
## Multiple R-squared:  0.1004, Adjusted R-squared:  0.08898 
## F-statistic: 8.795 on 5 and 394 DF,  p-value: 6.333e-08

La interpretación de los resultados es similar a una regresión por OLS con una variable numérica, excepto que conceptualizamos los coeficientes como cambios en la probabilidad de que la variable dependiente sea igual a 1 (\(Pr(Y|X) = 1\)) asociados a cambios en las variables independientes. En este caso, se trata de la probabilidad de que admit sea igual 1, o sea, de ser admitido al posgrado. Así, los resultados indican que:

Un incremento de 1 unidad en la nota del examen GRE, se asocia con un incremento de 0.04% en la probabilidad de ser amitido.
Un incremento de 1 unidad en el promedio crédito acumulado, se asocia con un incremento de 16% en la probabilidad de ser amitido.
Comparado con estudiantes de universidades en la categoría 1 del ránquin, estudiantes de la categoría 2 tienen 16% menos probabilidad de ser admitidos.
Comparado con estudiantes de universidades en la categoría 1 del ránquin, estudiantes de la categoría 3 tienen 29% menos probabilidad de ser admitidos.
Comparado con estudiantes de universidades en la categoría 1 del ránquin, estudiantes de la categoría 4 (la más baja) tienen 32% menos probabilidad de ser admitidos.
Finalmente, el intercepto nos dice que la probabilidad de ser admitido para un estudiante de una universidad en la categoría 1 del ránquin, con notas de 0 en el GRE y promedio crédito acumulado (GPA) de 0 es de -25%. Wait…what?

La interpretación del intercepto apunta a uno de los problemas principales del MPL: el modelo puede arrojar “probabilidades” por debajo de 0 o por encima de 1, lo cual no tiene sentido, ya que estas deben estar entre 0 y 1. Veamos esto con otro ejemplo. Calculemos el valor esperado de \(Y\) (la probabilidad esperada de ser admitido) para un estudiante de una universidad en la categoría 4 y los valores más bajos de GRE y GPA observados en los datos. Reemplazamos en la ecuación del modelo y extraemos los coeficientes usando la función coef():

coef(mpl)[[1]] + 
  coef(mpl)[[2]]*min(dat$gre) + 
  coef(mpl)[[3]]*min(dat$gpa) + 
  coef(mpl)[[4]]*0 +
  coef(mpl)[[5]]*0 +
  coef(mpl)[[6]]*1

## [1] -0.1359216

¡Este pobre estudiante tiene una probabilidad de -13% de ser admitido! Y esto no es un caso aislado, sino que pasa frecuentemente con el MPL. Para ver qué tan frecuentemente ocurre con nuestro modelo, miremos la distribución de las predicciones que estima. Para varias observaciones, el modelo predice probabilidades negativas (por debajo de 0):

mpl %>%
  # extraemos las predicciones o "fitted values" del modelo
  predict() %>% 
  as_tibble() %>% 
  ggplot(aes(x = value)) +
  geom_histogram(binwidth = 0.01) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
  scale_x_continuous(labels = scales::percent) +
  labs(x = expression(hat(Y)), y = "Número de observaciones")

Esto no significa que el MPL sea inútil. Estas probabilidades negativas (y en otras ocasiones, por encima de 1.0 o 100%) solo suceden para valores extremos de las variables independientes. Además, las podemos interpretar como que la probabilidad de que el evento ocurra es extremadamente baja (o extremadamente alta, en el caso de que estén por encima de 1). Finalmente, el MPL es utilizado frecuentemente en la literatura académica de experimentos o evaluación técnica de políticas públicas con experimentos pues da resultados similares a la regresión logística, es más fácil de interpretar y se ajusta mejor a algunas nociones de causalidad y efectos.

Pese a sus limitaciones, las predicciones de un modelo como este tienen valor. Veamos los efectos de las distintas variables independientes usando ggpredict() de ggeffects:

mpl %>% 
  ggpredict()

## $gre
## # Predicted values of admit
## 
## gre | Predicted |       95% CI
## ------------------------------
## 220 |      0.36 | [0.17, 0.56]
## 295 |      0.40 | [0.23, 0.56]
## 365 |      0.43 | [0.28, 0.57]
## 435 |      0.46 | [0.32, 0.59]
## 510 |      0.49 | [0.37, 0.61]
## 585 |      0.52 | [0.41, 0.63]
## 655 |      0.55 | [0.44, 0.66]
## 800 |      0.61 | [0.47, 0.75]
## 
## Adjusted for:
## *  gpa = 3.39
## * rank =    1
## 
## $gpa
## # Predicted values of admit
## 
##  gpa | Predicted |       95% CI
## -------------------------------
## 2.00 |      0.30 | [0.09, 0.52]
## 2.50 |      0.38 | [0.22, 0.54]
## 3.00 |      0.46 | [0.34, 0.58]
## 3.50 |      0.54 | [0.43, 0.65]
## 4.00 |      0.62 | [0.48, 0.75]
## 
## Adjusted for:
## *  gre = 587.70
## * rank =      1
## 
## $rank
## # Predicted values of admit
## 
## rank | Predicted |       95% CI
## -------------------------------
## 1    |      0.52 | [0.41, 0.63]
## 2    |      0.36 | [0.29, 0.43]
## 3    |      0.23 | [0.15, 0.31]
## 4    |      0.20 | [0.09, 0.30]
## 
## Adjusted for:
## * gre = 587.70
## * gpa =   3.39
## 
## attr(,"class")
## [1] "ggalleffects" "list"        
## attr(,"model.name")
## [1] "."

O para un “perfil” específico:

mpl %>% 
  ggpredict(
    # mantenemos gre constante en la media, variamos las otras dos
    terms = c("gpa", "rank[1, 4]"), 
  ) %>%
  plot() +
  scale_y_continuous(labels = scales::percent) +
  labs(title = "Valores esperados: probabilidad de ser admitido",
       x = "GPA", y = "Pr(Y=1)", color = "Ránquin")

## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

Vemos que la probabilidad de ser admitido aumenta a medida que sube el promedio crédito acumulado o GPA, pero que los estudiantes que provenienen de universidades con un mejor ránquin tienen una ventaja comparados con los de universidades menos prestigiosas.

9.4 Regresión logística

La alternativa al MPL es la regresión logística o el modelo logit.¹¹. A diferencia de OLS, la regresión logística, está diseñada explícitamente para modelar variables dependientes categóricas binarias y se estima via máxima verosimilitud (o maximum likelihood estimation, MLE), en vez de mínimos cuadrados ordinarios/OLS.¹² El logit está diseñado para ver modelar la asociación entre una serie de variables independientes y la probabilidad de que \(Y = 1\) dado unas \(X\).

Estimemos el mismo modelo que trabajamos arriba, pero esta vez usando regresión logística. En R, usamos la funciónm glm() y especificamos que el modelo usa una distribución binomial, lo cual le indica a R que es una regresión logística:

logit <- glm(admit_fct ~ gre + gpa + rank, 
             data = dat, 
             # tipo de modelo/distribucion
             family = "binomial")

Aparte: podemos estimar un modelo de regresión lineal con glm() y obtener los mismos resultados que con lm() - comparen con los resultados del MPL:

glm(admit ~ gre + gpa + rank,
    data = dat, 
    family = "gaussian") %>%
  summary() # ver coeficientes

## 
## Call:
## glm(formula = admit ~ gre + gpa + rank, family = "gaussian", 
##     data = dat)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.7022  -0.3288  -0.1922   0.4952   0.9093  
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.2589102  0.2159904  -1.199   0.2314    
## gre          0.0004296  0.0002107   2.038   0.0422 *  
## gpa          0.1555350  0.0639618   2.432   0.0155 *  
## rank2       -0.1623653  0.0677145  -2.398   0.0170 *  
## rank3       -0.2905705  0.0702453  -4.137 4.31e-05 ***
## rank4       -0.3230264  0.0793164  -4.073 5.62e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.1979062)
## 
##     Null deviance: 86.677  on 399  degrees of freedom
## Residual deviance: 77.975  on 394  degrees of freedom
## AIC: 495.12
## 
## Number of Fisher Scoring iterations: 2

Volviendo al modelo logit, los coeficientes que vemos en los resultados del modelo son “log-odds” y, por el momento, solo interpretamos la dirección de la relación (el signo del coeficiente, positivo o negativo) y la significancia estadística de los mismos, no su magnitud. Los log-odds son el logaritmo de los odds:¹³ la tasa o razón entre el número veces que algo sucede y el número de veces que no sucede. Mientras, probabilidad es la tasa de casos positivos a casos posibles. Si yo llego tarde a clase 1 de cada 5 veces, los odds de yo llegar tarde son \(\frac{1}{4} = 0.25\) (y los log-odds son \(\log(0.25) = -1.39\)), pero la probabilidad de que llegue tarde es \(\frac{1}{5} = 0.20\). Son formas distintas de expresar ideas muy similares.

En ese orden de ideas, los resultados del modelo logístico que estimamos corroboran lo encontrado con el MPL - aumentos en el GRE y GPA se asocian con aumentos la probabilidad de admisión, mientras que los estudiantes de universidades de universidades mejor ranqueadas tienen una mayor probabilidad de admisión que los demás:

summary(logit)

## 
## Call:
## glm(formula = admit_fct ~ gre + gpa + rank, family = "binomial", 
##     data = dat)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6268  -0.8662  -0.6388   1.1490   2.0790  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.989979   1.139951  -3.500 0.000465 ***
## gre          0.002264   0.001094   2.070 0.038465 *  
## gpa          0.804038   0.331819   2.423 0.015388 *  
## rank2       -0.675443   0.316490  -2.134 0.032829 *  
## rank3       -1.340204   0.345306  -3.881 0.000104 ***
## rank4       -1.551464   0.417832  -3.713 0.000205 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 499.98  on 399  degrees of freedom
## Residual deviance: 458.52  on 394  degrees of freedom
## AIC: 470.52
## 
## Number of Fisher Scoring iterations: 4

Podemos convertir los log-odds a odds (también llamados odds ratios) usando la función exponencial exp(). En este caso, valores por encima de 1 indican que hay una asociación positiva (aumentos en la variable independiente se asocian con aumentos en \(Pr(Y) = 1\)) y valores por debajo de 1 indican una asociación negativa (reducciones en \(Pr(Y) = 1\)):

coef(logit) %>%
  as_tibble() %>%
  mutate(odds = exp(value))

## # A tibble: 6 x 2
##      value   odds
##      <dbl>  <dbl>
## 1 -3.99    0.0185
## 2  0.00226 1.00  
## 3  0.804   2.23  
## 4 -0.675   0.509 
## 5 -1.34    0.262 
## 6 -1.55    0.212

Veamos una tabla resumen comparando ambos modelos -el MPL y el logit:

modelsummary(
  list(
    "MPL" = mpl,
    "Logit" = logit
  ),
  stars = TRUE
)

## Warning: In version 0.8.0 of the `modelsummary` package, the default significance markers produced by the `stars=TRUE` argument were changed to be consistent with R's defaults.
## This warning is displayed once per session.

	MPL	Logit
(Intercept)	-0.259	-3.990***
	(0.216)	(1.140)
gre	0.000*	0.002*
	(0.000)	(0.001)
gpa	0.156*	0.804*
	(0.064)	(0.332)
rank2	-0.162*	-0.675*
	(0.068)	(0.316)
rank3	-0.291***	-1.340***
	(0.070)	(0.345)
rank4	-0.323***	-1.551***
	(0.079)	(0.418)
Num.Obs.	400	400
R2	0.100
R2 Adj.	0.089
AIC	495.1	470.5
BIC	523.1	494.5
Log.Lik.	-240.560	-229.259
F	8.795
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

Por defecto, modelsummary nos muestra log-odds en los coeficientes, pero podemos exponenciarlos y presentar odds-ratios en cambio:

modelsummary(
  list(
    "MPL" = mpl,
    "Logit" = logit
  ),
  stars = TRUE, 
  exponentiate = TRUE
)

	MPL	Logit
(Intercept)	-0.259	0.019***
	(0.216)	(1.140)
gre	0.000*	1.002*
	(0.000)	(0.001)
gpa	0.156*	2.235*
	(0.064)	(0.332)
rank2	-0.162*	0.509*
	(0.068)	(0.316)
rank3	-0.291***	0.262***
	(0.070)	(0.345)
rank4	-0.323***	0.212***
	(0.079)	(0.418)
Num.Obs.	400	400
R2	0.100
R2 Adj.	0.089
AIC	495.1	470.5
BIC	523.1	494.5
Log.Lik.	-240.560	-229.259
F	8.795
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

En últimas, los resultados parecen ser prácticamente idénticos en términos de dirección y significancia estadística. ¿Por qué hacer una regresión logística, entonces? Porque el modelo logit nos permite pasar entre log-odds, odds y probabilidades de forma tal que podemos hacer mejores predicciones y mostrar “efectos” más efectivamente.

9.4.1 Coeficientes como probabilidades

La mejor forma de interpretar los resultados de un modelo de regresión logística (la forma más fácil y más informativa) pasa por convertir los log-odds a probabilidades. Esto lo hacemos pasando los log-odds por la función logística, de ahí el nombre:

\[ \text{logistic}(\alpha) = \frac{\exp(\alpha)}{1 + \exp(\alpha)} \]

En R, podríamos calcular probabilidades usando la función logística así:

exp(coeficiente) / (1 + exp(coeficiente))

Sin embargo, ya existe la función plogis, que es el “link” o vínculo entre probabilidades y log-odds. Gráficamente, la función logística se ve así (y describe mejor la relación entre una variable numérica y una categórica que el OLS que hicimos arriba):

tibble(x = -15:15, y = plogis(x)) %>%
  ggplot(aes(x, y)) +
  geom_line(size = 1, color = "red")

Volviendo a nuestro modelo de regresión, la probabilidad que resulta de pasar el coeficiente (en log-odds) por la función logística es la probabilidad de que \(Y\) sea igual a 1 \(Pr(Y) = 1\) en la mitad de la curva de la función logística:

coef(logit) %>%
  as_tibble() %>%
  mutate(
    odds = exp(value),
    prob = plogis(value)
  )

## # A tibble: 6 x 3
##      value   odds   prob
##      <dbl>  <dbl>  <dbl>
## 1 -3.99    0.0185 0.0182
## 2  0.00226 1.00   0.501 
## 3  0.804   2.23   0.691 
## 4 -0.675   0.509  0.337 
## 5 -1.34    0.262  0.207 
## 6 -1.55    0.212  0.175

Otra forma de interpretar los coeficientes de una regresión logística en términos de probabilidades es usando la regla de “dividir por 4”. Tomamos el coeficiente de la regresión, lo dividimos por 4 y el resultado es una aproximación al “efecto” máximo de \(X\) sobre la probabilidad de que \(Y\) ocurra:

coef(logit)/4

##   (Intercept)           gre           gpa         rank2         rank3 
## -0.9974947683  0.0005661064  0.2010093873 -0.1688607320 -0.3350509791 
##         rank4 
## -0.3878659192

Así, por ejemplo, un incremento de 1 unidad en la variable gpa se asocia con un incremento aproximado del 20% en la probabilidad de ser admitido al posgrado. En este caso, vemos una diferencia entre este valor y el que obtuvimos con el MPL, que era de aproximadamente 15.5%.

9.4.1.1 Visualizar probabilidades

Pese a que podemos usar log-odds, odds-ratios y probabilidades, los resultados de una regresión logística son más interpretables de forma gráfica. Vamos a usar ggpredict() nuevamente, especificando las variables independientes cuyos efectos queremos ver usando terms =:

logit %>%
  ggpredict(terms = c("gpa", "rank")) %>%
  plot() +
  labs(title = "Valores esperados: probabilidad de ser admitido",
  color = "Ránquin", x = "GPA", y = "Pr(Y=1)")

## Data were 'prettified'. Consider using `terms="gpa [all]"` to get smooth plots.

Una representación gráfica del modelo de regresión logística nos permite evidenciar dos elementos muy importantes en su interpretación:

La pendiente de la línea cambia levemente a lo largo de su trayecto, indicando que se trata de una curva -y no una recta como en OLS- y que la asociación entre \(X\) (GPA) y \(Y\) (probabilidad de ser admitido) no es constante, sino que depende del valor de \(X\) misma.
La pendiente de la curva para cada grupo definido por la variable categórica \(Z\) (ránquin) también es distinta, lo cual indica que la asociación entre \(X\) y \(Y\) depende de las demás variables incluidas en el modelo. En otras palabras, ¡hay una interacción inherente en una regresión logística, aunque no la hayamos especificado!

9.4.2 Evaluar el modelo

Hemos estimado e interpretado los resultados de un modelo de regresión logística. Ahora, ¿cómo lo evaluamos? Miremos dos opciones.

9.4.2.1 Variables omitidas

Esta es una evaluación tanto teórica, como metodológica. Al igual que en la regresión lineal por OLS, para que uno modelo de regresión logística tenga coeficientes insesgados, es necesario que hayamos incluido todas las variables relevantes (causales). Así mismo, es importante que no hayamos incluido otras variables irrelevantes. En últimas, debemos evitar el sesgo por variables omitidas. Si queremos hacer inferencia causal, este punto es crítico.

9.4.2.2 Predicciones y falsos positivos

Otra forma común de evaluar un modelo de regresión logística es comparando las predicciones del modelo con los valores reales. Vale aclarar que en este caso estamos minimizando la importancia de la estimación de un efecto causal y, más bien, enfocándonos en la capacidad predictiva del modelo. Estos son dos objetivos muy distintos.

Recordemos que para cada candidato en nuestros datos, sabemos si fue admitido o no a un programa posgrado (el valor de \(Y\), admit). Además, una vez estimamos un modelo de regresión logística, estimamos la probabilidad esperada de ser admitido para cada candidato.

Para facilitar la comparación entre probabilidades (probabilidad de ser admitido o no) y valores de 0 o 1 (admitido o no), debemos imponer un umbral a partir del cual clasificamos nuestras predicciones de la admisión de cada candidato como “admitido” o “no admitido”. Por ejemplo, vamos a asumir que si el modelo estima una probabilidad de admisión por encima del 50% (0.5) para una observación, la predicción del modelo es que ese individuo será “admitido” al programa de posgrado. Usando la función augment de broom podemos obtener los valores esperados (predicciones) de \(Y\) para cada observación que entró al modelo:

preds <- logit %>% 
  # agregamos las predicciones a los datos
  # especificamos que queremos probabilidades 
  augment(type.predict = "response") %>% 
  # seleccionamos unas columnas para simplificar
  select(admit_fct, gre, gpa, rank, prob_admit = .fitted) %>% 
  # aplicamos nuestra regla para clasificar predicciones como admitido
  mutate(pred_admit = if_else(prob_admit > 0.5, "Sí", "No"))
preds

## # A tibble: 400 x 6
##    admit_fct   gre   gpa rank  prob_admit pred_admit
##    <fct>     <dbl> <dbl> <fct>      <dbl> <chr>     
##  1 No          380  3.61 3          0.173 No        
##  2 Sí          660  3.67 3          0.292 No        
##  3 Sí          800  4    1          0.738 Sí        
##  4 Sí          640  3.19 4          0.178 No        
##  5 No          520  2.93 4          0.118 No        
##  6 Sí          760  3    2          0.370 No        
##  7 Sí          560  2.98 1          0.419 No        
##  8 No          400  3.08 2          0.217 No        
##  9 Sí          540  3.39 3          0.201 No        
## 10 No          700  3.92 2          0.518 Sí        
## # ... with 390 more rows

Para ver el desempeño del modelo, podemos ver:

Positivos verdaderos: datos dicen sí, modelo predice sí.
Falsos positivos: datos dicen no, modelo predice sí.
Negativos verdaderos: datos dicen no, modelo predice no.
Falsos negativos: datos dicen sí, modelo predice no.

Organicemos una tabla cruzada con esa información (también conocida como una matriz de confusión):

preds %>%
  group_by(admit_fct, pred_admit) %>%
  summarize(casos = n()) %>%
  ungroup() %>%
  mutate(prop = casos/sum(casos)) %>%
  pivot_wider(admit_fct, names_from = pred_admit, values_from = prop)

## `summarise()` has grouped output by 'admit_fct'. You can override using the `.groups` argument.

## # A tibble: 2 x 3
##   admit_fct    No     Sí
##   <fct>     <dbl>  <dbl>
## 1 No        0.635 0.0475
## 2 Sí        0.242 0.075

El que estimamos modelo predice correctamente (No-No y Sí-Sí) aproximadamente el 71% de las observaciones. Sin embargo, tenemos aproximadamente un 24% de falsos negativos y 5% de falsos positivos. Si nuestra pretensión es construir un modelo que identifique correctamente los admitidos, claramente podríamos mejorar. Esto lo haríamos a) agregando variables relevantes, b) transformando las existentes de forma tal que haya una relación más clara con la variable dependiente o c) cambiando el modelo completamente.

9.5 Conclusión

Con frecuencia queremos explicar, entender, describir o predecir variable dependientes categóricas binarias. Para esto, podemos estimar un modelo de probabilidad lineal, pero es mejor utilizar una regresión logística. Un modelo de regresión logística nos permite estimar la probabilidad de que un evento ocurra dadas unas variables independientes. Así mismo, nos permite cómo cambios en las variables independientes se asocian con cambios en la probabilidad de que ocurra un evento. Finalmente, podemos evaluar estos modelos según su capacidad de predecir correctamente.

9.6 Ejercicios

9.6.1 Explicando el voto nominal

¿Por qué algunos legisladores votan a favor de un proyecto de ley? ¿Cómo se relaciona la identificación partidista con la probabilidad de que un legislador vote “sí”? Utilicemos datos (el archivo nafta.dta disponible en https://github.com/josefortou/lab-book/tree/master/data/nafta.dta) para intentar responder estas preguntas usando modelos de regresión logística.

9.6.1.1 Los datos

Miremos cómo votaron los representates a la cámara de Estados Unidos frente al Tratado de Libre Comercio de América del Norte (TLCAN o NAFTA):

nafta <- read_dta("data/nafta.dta") %>%
  zap_label()
nafta %>% 
  count(vote)

## # A tibble: 2 x 2
##    vote     n
##   <dbl> <int>
## 1     0   200
## 2     1   234

234 votaron a favor y 200 en contra. Además de la variable de voto a favor/en contra (vote que toma los valores “Yes” y “No”), tenemos tres variables más:

democrat: dummy de membrecía partidista (Republicano o Demócrata).
pcthispc: porcentaje de la población que es de origen latino en el distrito de donde viene el representante.
cope93: indicador COPE de posiciones pro-sindicato de cada legislador (0=anti-sindicatos, 100=pro-sindicatos).

Recodifiquemos las variables categóricas:

nafta <- nafta %>%
  mutate(
    vote = factor(vote, levels = c(0, 1), labels = c("No", "Yes")),
    democrat = factor(democrat, levels = c(0, 1), labels = c("Rep./other", "Dem."))
  )

Ahora, veamos cuántos representantes votaron a favor o en contra en cada partido:

nafta %>% 
  count(vote, democrat)

## # A tibble: 4 x 3
##   vote  democrat       n
##   <fct> <fct>      <int>
## 1 No    Rep./other    47
## 2 No    Dem.         153
## 3 Yes   Rep./other   133
## 4 Yes   Dem.         101

9.6.1.2 Hipótesis

Tenemos las siguientes hipótesis:

Mayores valores COPE (más pro-sindicalismo) se asocian con una menor probabilidad de votar por NAFTA.
El efecto de COPE sobre la probabilidad de votar por NAFTA depende de la afiliación partidista: el efecto será mayor para Demócratas que para Republicanos.

Además, controlamos por una explicación alternativa: es posible que una mayor población latina en el distrito de origen se asocie con una mayor probabilidad de votar por NAFTA por razones de lobby, presión política y deseos de reelección.

Así que estimamos el siguiente modelo:

\[ Pr(Y = 1|X) = \lambda(\beta_0 + \beta_1\text{Dem.} + \beta_2\text{COPE} + \beta_3\text{Dem.}\times\text{COPE} + \beta_4\text{Porc. latino}) \]

Donde \(\lambda\) es la función de distribución acumulada (CDF) logística.

9.6.1.3 Estimación

Estime el modelo usando la función glm().

9.6.1.4 Resultados

Presente los resultados en una tabla de regresión:

Presente resultados que evalúen las hipótesis utilizando gráficas de probabilidades esperadas/predicciones:

Otros modelos como el probit son similares, pero menos populares y con frecuencia equivalentes en su interpretación o resultados.↩︎
Mientras que OLS busca los parámetros que minimizan una distancia -los cuadrados de los residuos-, MLE busca los parámetros que maximizan la probabilidad o likelihood de observar los datos que tenemos.↩︎
En español, a veces se usan los términos “oportunidades”, “probabilidades” o “momios”. Porque no hay una terminología equivalente en español, usamos odds.↩︎