Capítulo 9 Regresión logística
9.1 Resumen
En este instructivo, estimamos e interpretamos modelos de regresión logística para variables dependientes categóricas.
Vamos a utilizar las siguientes librerías:
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.0.3
## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.3.2 v purrr 0.3.4
## v tibble 3.0.6 v dplyr 1.0.4
## v tidyr 1.1.2 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0
## Warning: package 'tibble' was built under R version 4.0.3
## Warning: package 'dplyr' was built under R version 4.0.3
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(ggeffects) # efectos en modelos de regresion
library(sjPlot) # tablas de regresion
library(haven) # datos en formato .dta
9.2 Repaso
Estimar modelos de regresión lineal:
- Operacionalización de modelos teóricos para evaluar hipótesis sobre relaciones entre variables.
- Interpretar coeficientes: dirección, magnitud, significancia.
- Regresión múltiple para controlar por otros factores.
9.3 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). Una mejor opción es utilizar la regresión logística o logit. En este instructivo, veremos cómo podemos usar ambos.
9.4 Datos
Pero primero, necesitamos datos. Carguemos unos datos directamente de la web. Se trata de información sobre la admisión a programas de posgrado universitario en una univerdiad de Estados Unidos:
<- read_csv("https://stats.idre.ucla.edu/stat/data/binary.csv") dat
## Parsed with column specification:
## cols(
## admit = col_double(),
## gre = col_double(),
## gpa = col_double(),
## rank = col_double()
## )
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 cuatro variables:
admit
: dummy de admisión (0=no admitio; 1=admitido).gre
: nota en el examen GRE. Puede tomar valores de 200 a 800.gpa
: promedio crédito acumulado de pregrado. Puede tomar valores de 0 a 4.rank
: ránquin de la universidad de pregrado. Hay 4 categorías, de mayor (1) a menor (4) calidad.
Estamos interesados en la probabilidad de que un individuo sea admitido a la universidad (admit = 1
). Convertimos dos variables categóricas -amit
y rank
- a tipo factor, porque pese a que conceptualmente son categóricas, están codificadas como numéricas en R. Hacemos esto con la funciones case_when()
(para recodificar números como categorías) y as.factor()
(un cambio de clase simple):
<- dat %>%
dat mutate(
admit_fct = case_when(
== 0 ~ "No", # usar los valores de la variable admit
admit == 1 ~ "Sí"
admit
),admit_fct = as.factor(admit_fct),
rank_fct = as.factor(rank)
)
Ya que 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",
y = "Número de observaciones")
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
Veamos una de estas relaciones como una gráfica de dispersión con una línea de tendencia (¡estimada por ggplot2
con 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. Busquemos una alternativa.
9.5 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
):
<- lm(
mpl ~ gre + gpa + rank_fct,
admit data = dat
)
Usemos summary()
para ver un resumen de los resultados del modelo:
summary(mpl)
##
## Call:
## lm(formula = admit ~ gre + gpa + rank_fct, 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 *
## rank_fct2 -0.1623653 0.0677145 -2.398 0.0170 *
## rank_fct3 -0.2905705 0.0702453 -4.137 4.31e-05 ***
## rank_fct4 -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 es similar a OLS, 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 a la universidad. 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%. ¿¡Cómo!?
La interpretación del intercepto apunta a uno de los problemas principales del MPL: el modelo puede arrojar “probabilidades” por debajo de 0. Veamos otros 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 predict() %>% # extrae las predicciones o "fitted values" del modelo
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. 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.
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
## # x = gre
##
## x | Predicted | SE | 95% CI
## -------------------------------------
## 220 | 0.36 | 0.10 | [0.17, 0.56]
## 295 | 0.40 | 0.09 | [0.23, 0.56]
## 365 | 0.43 | 0.08 | [0.28, 0.57]
## 435 | 0.46 | 0.07 | [0.32, 0.59]
## 510 | 0.49 | 0.06 | [0.37, 0.61]
## 585 | 0.52 | 0.06 | [0.41, 0.63]
## 655 | 0.55 | 0.06 | [0.44, 0.66]
## 800 | 0.61 | 0.07 | [0.47, 0.75]
##
## Adjusted for:
## * gpa = 3.39
## * rank_fct = 1
##
##
## $gpa
##
## # Predicted values of admit
## # x = gpa
##
## x | Predicted | SE | 95% CI
## --------------------------------------
## 2.00 | 0.30 | 0.11 | [0.09, 0.52]
## 2.50 | 0.38 | 0.08 | [0.22, 0.54]
## 3.00 | 0.46 | 0.06 | [0.34, 0.58]
## 3.50 | 0.54 | 0.06 | [0.43, 0.65]
## 4.00 | 0.62 | 0.07 | [0.48, 0.75]
##
## Adjusted for:
## * gre = 587.70
## * rank_fct = 1
##
##
## $rank_fct
##
## # Predicted values of admit
## # x = rank_fct
##
## x | Predicted | SE | 95% CI
## -----------------------------------
## 1 | 0.52 | 0.06 | [0.41, 0.63]
## 2 | 0.36 | 0.04 | [0.29, 0.43]
## 3 | 0.23 | 0.04 | [0.15, 0.31]
## 4 | 0.20 | 0.05 | [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(
terms = c("gpa", "rank_fct[1, 4]"), # mantiene gre constante en la media
%>%
) plot() +
scale_y_continuous(labels = scales::percent) +
labs(title = "Valores esperados: probabilidad de ser admitido",
x = "GPA", y = "Pr(Y=1)", color = "Rank pregrado")
## 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.6 Regresión logística
La alternativa al MPL es la regresión logística o el modelo logit.12. A diferencia de OLS, la regresión logística, está diseñada explícitamente para modelar variables dependientes categóricas (específicamente binarias) y se estima via máxima verosimilitud (o maximum likelihood estimation, MLE), no mínimos cuadrados ordinarios.13 El logit está diseñadao para ver la asociación entre una serie de variables independientes y la probabilidad de que Y=1 dado X (Pr(Y|X)=1).
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:
<- glm( # modelos lineales generalizados estimados por MLE
logit ~ gre + gpa + rank_fct, # formula
admit_fct data = dat, # datos
family = binomial() # tipo de modelo/distribucion
)
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( # modelos lineales generalizados estimados por MLE
~ gre + gpa + rank, # formula
admit data = dat, # datos
family = "gaussian" # tipo de modelo/distribucion
%>%
) summary() # ver coeficientes
##
## Call:
## glm(formula = admit ~ gre + gpa + rank, family = "gaussian",
## data = dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.6617 -0.3417 -0.1947 0.5061 0.9556
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.1824127 0.2169695 -0.841 0.4010
## gre 0.0004424 0.0002101 2.106 0.0358 *
## gpa 0.1510402 0.0633854 2.383 0.0176 *
## rank -0.1095019 0.0237617 -4.608 5.48e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.1978669)
##
## Null deviance: 86.677 on 399 degrees of freedom
## Residual deviance: 78.355 on 396 degrees of freedom
## AIC: 493.07
##
## Number of Fisher Scoring iterations: 2
Volviendo al modelo logit, 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. Los log-odds son el logaritmo de los odds:14 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 14=0.25 (y los log-odds son log(0.25)=−1.39), pero la probabilidad de que llegue tarde es 15=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_fct, 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 *
## rank_fct2 -0.675443 0.316490 -2.134 0.032829 *
## rank_fct3 -1.340204 0.345306 -3.881 0.000104 ***
## rank_fct4 -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 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. Por defecto, tab_model()
nos muestra odds en los coeficientes:
tab_model(
mpl, logit,dv.labels = c("MPL", "Logit")
)
MPL | Logit | |||||
---|---|---|---|---|---|---|
Predictors | Estimates | CI | p | Odds Ratios | CI | p |
(Intercept) | -0.26 | -0.68 – 0.17 | 0.231 | 0.02 | 0.00 – 0.17 | <0.001 |
gre | 0.00 | 0.00 – 0.00 | 0.042 | 1.00 | 1.00 – 1.00 | 0.038 |
gpa | 0.16 | 0.03 – 0.28 | 0.015 | 2.23 | 1.17 – 4.32 | 0.015 |
rank_fct [2] | -0.16 | -0.30 – -0.03 | 0.017 | 0.51 | 0.27 – 0.94 | 0.033 |
rank_fct [3] | -0.29 | -0.43 – -0.15 | <0.001 | 0.26 | 0.13 – 0.51 | <0.001 |
rank_fct [4] | -0.32 | -0.48 – -0.17 | <0.001 | 0.21 | 0.09 – 0.47 | <0.001 |
Observations | 400 | 400 | ||||
R2 / R2 adjusted | 0.100 / 0.089 | 0.102 |
Los resultados son prácticamente idénticos… ¿por qué hacer una regresión logística, entonces? Porque el logit nos permite pasar entre log-odds, odds y probabilidades de forma tal que podemos hacer mejores predicciones.
9.6.1 Coeficientes como probabilidades
La mejor forma de interpretar los resultados de un modelo de regresión logística (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:
logistic(α)=exp(α)1+exp(α)
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 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
9.6.1.1 Visualizar probabilidades
Estos resultados son más interpretables de forma gráfica. Podemos usar ggpredict()
nuevamente, especificando las variables independientes cuyos efectos queremos ver usando terms =
:
%>%
logit ggpredict(terms = c("gpa", "rank_fct")) %>%
plot() +
labs(title = "Valores esperados: probabilidad de ser admitido",
color = "Rank", 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 que la asociación entre X (GPA) y Y (probabilidad de ser admitido) no es constante, sino que depende del valor de X.
- 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!
9.6.2 Evaluar el modelo
¿Cómo evaluar un modelo de regresión logística? Miremos dos opciones.
9.6.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.
9.6.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. Recordemos que para cada observación en nuestros datos, sabemos si fue admitido o no. Además, podemos calcular la probabilidad de ser admitido a partir de los resultados del modelo.
Para facilitar la comparación entre probabilidades (probabilidad de ser admitido o no) y valores de 0 o 1 (admitido o no), asumimos que si el modelo estimó 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”:
<- as_tibble(logit$model)
modelo_aug <- modelo_aug %>%
modelo_aug mutate(fitted = logit$fitted.values,
residuals = logit$residuals,
std_residuals = rstandard(logit),
cooks_d = cooks.distance(logit))
<- modelo_aug %>%
preds mutate(
model_prob = plogis(fitted), # probabilidad de ser admitido
model_pred = if_else(
> 0.5, "Sí", "No" # convertir a binario (clasificar)
model_prob
)
) preds
## # A tibble: 400 x 10
## admit_fct gre gpa rank_fct fitted residuals std_residuals cooks_d
## <fct> <dbl> <dbl> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 No 380 3.61 3 0.173 -1.21 -0.621 0.000579
## 2 Sí 660 3.67 3 0.292 3.42 1.58 0.00450
## 3 Sí 800 4 1 0.738 1.35 0.788 0.00145
## 4 Sí 640 3.19 4 0.178 5.61 1.87 0.0132
## 5 No 520 2.93 4 0.118 -1.13 -0.505 0.000302
## 6 Sí 760 3 2 0.370 2.70 1.43 0.00631
## 7 Sí 560 2.98 1 0.419 2.39 1.33 0.00535
## 8 No 400 3.08 2 0.217 -1.28 -0.704 0.000610
## 9 Sí 540 3.39 3 0.201 4.98 1.80 0.00573
## 10 No 700 3.92 2 0.518 -2.07 -1.22 0.00268
## # ... with 390 more rows, and 2 more variables: model_prob <dbl>,
## # model_pred <chr>
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 (una matriz de confusión):
%>%
preds group_by(admit_fct, model_pred) %>% # agrupar
summarize(casos = n()) %>%
ungroup() %>%
mutate(prop = casos/sum(casos)) %>%
pivot_wider(admit_fct, names_from = model_pred, values_from = prop) # organizar
## `summarise()` has grouped output by 'admit_fct'. You can override using the `.groups` argument.
## # A tibble: 2 x 2
## admit_fct Sí
## <fct> <dbl>
## 1 No 0.682
## 2 Sí 0.318
El modelo predice correctamente (No-No y Sí-Sí) el 71% de las observaciones. Sin embargo, tenmos aproximadamente un 24% de falsos negativos. Si nuestra pretensión es construir un modelo que identifique correctamente los admitidos, claramente podríamos mejorar. Esto lo haríamos a) agregando variables, b) transformando las existentes de forma tal que haya una relación más clara con la variable dependiente o c) cambiando el modelo.
9.7 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.8 Ejercicios
9.8.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
disponibles en https://github.com/josefortou/lab-book/tree/master/data) para intentar responder estas preguntas.
9.8.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):
<- read_dta("data/nafta.dta") %>%
nafta zap_label()
nafta
## # A tibble: 434 x 4
## vote democrat pcthispc cope93
## <dbl> <dbl> <dbl> <dbl>
## 1 0 1 8 92
## 2 0 1 1 92
## 3 0 1 40 100
## 4 1 0 0 42
## 5 0 0 22 8
## 6 1 0 5 0
## 7 1 1 66 83
## 8 1 1 14 83
## 9 0 1 0 100
## 10 0 1 3 92
## # ... with 424 more rows
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.8.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)=λ(β0+β1Dem.+β2COPE+β3Dem.×COPE+β4Porc. latino)
Donde λ es la función de distribución acumulada (CDF) logística.
9.8.1.3 Estimación
Estime el modelo usando la función glm()
.
9.8.1.4 Resultados
Presente los resultados en una tabla de regresión:
Presente los resultados utilizando gráficas:
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.↩︎