1.-Introducción

La Regresión Logística Simple, desarrollada por David Cox en 1958, es un método de regresión que permite estimar la probabilidad de una variable cualitativa binaria en función de una variable cuantitativa. Una de las principales aplicaciones de la regresión logística es la de clasificación binaria, en el que las observaciones se clasifican en un grupo u otro dependiendo del valor que tome la variable empleada como predictor. Por ejemplo, clasificar a un individuo desconocido como hombre o mujer en función del tamaño de la mandíbula.

Es importante tener en cuenta que, aunque la regresión logística permite clasificar, se trata de un modelo de regresión que modela el logaritmo de la probabilidad de pertenecer a cada grupo. La asignación final se hace en función de las probabilidades predichas.

La existencia de una relación significativa entre una variable cualitativa con dos niveles y una variable continua se puede estudiar mediante otros test estadísticos tales como t-test o ANOVA (un ANOVA de dos grupos es equivalente al t-test). Sin embargo, la regresión logística permite además calcular la probabilidad de que la variable dependiente pertenezca a cada una de las dos categorías en función del valor que adquiera la variable independiente. Supóngase que se quiere estudiar la relación entre los niveles de colesterol y los ataques de corazón. Para ello, se mide el colesterol de un grupo de personas y durante los siguientes 20 años se monitoriza que individuos han sufrido un ataque.

Un t-test entre los niveles de colesterol de las personas que han sufrido ataque vs las que no lo han sufrido permitiría contrastar la hipótesis de que el colesterol y los ataques al corazón están asociados. Si además se desea conocer la probabilidad de que una persona con un determinado nivel de colesterol sufra un infarto en los próximos 20 años, o poder conocer cuánto tiene que reducir el colesterol un paciente para no superar un 50% de probabilidad de padecer un infarto en los próximos 20 años, se tiene que recurrir a la regresión logística.

.- DATOS

# Librerías
#----------------------
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ISLR)
datos <- Default
head(datos)
##   default student   balance    income
## 1      No      No  729.5265 44361.625
## 2      No     Yes  817.1804 12106.135
## 3      No      No 1073.5492 31767.139
## 4      No      No  529.2506 35704.494
## 5      No      No  785.6559 38463.496
## 6      No     Yes  919.5885  7491.559
# Recodificación y selección de variables
# -----------------------------------------
datos <- datos %>% 
         select(default, balance) %>%
         mutate(default = recode(default, 
                                 "No" = 0, 
                                 "Yes" = 1))

head(datos)
##   default   balance
## 1       0  729.5265
## 2       0  817.1804
## 3       0 1073.5492
## 4       0  529.2506
## 5       0  785.6559
## 6       0  919.5885
str(datos)
## 'data.frame':    10000 obs. of  2 variables:
##  $ default: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ balance: num  730 817 1074 529 786 ...
summary(datos)
##     default          balance      
##  Min.   :0.0000   Min.   :   0.0  
##  1st Qu.:0.0000   1st Qu.: 481.7  
##  Median :0.0000   Median : 823.6  
##  Mean   :0.0333   Mean   : 835.4  
##  3rd Qu.:0.0000   3rd Qu.:1166.3  
##  Max.   :1.0000   Max.   :2654.3
# Regresión MCO lineal de default ~ Balance
# -----------------------------------------------
modelo_lineal = lm(default ~balance, data = datos)
summary(modelo_lineal)
## 
## Call:
## lm(formula = default ~ balance, data = datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.23533 -0.06939 -0.02628  0.02004  0.99046 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -7.519e-02  3.354e-03  -22.42   <2e-16 ***
## balance      1.299e-04  3.475e-06   37.37   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1681 on 9998 degrees of freedom
## Multiple R-squared:  0.1226, Adjusted R-squared:  0.1225 
## F-statistic:  1397 on 1 and 9998 DF,  p-value: < 2.2e-16
# Representación de los datos y el modelo
# -----------------------------------------
library(ggplot2)
ggplot(data = datos, aes(x = balance, y = default)) + 
  geom_point(aes(color = as.factor(default)), shape = 1) + 
  geom_smooth(method = "lm", color = "grey20", se = TRUE)+
  theme_bw() + 
  labs(title = "Regresión lineal por MCO", y = "Probabilidad de Default") +
  theme(legend.position = "bottom")
## `geom_smooth()` using formula = 'y ~ x'

En este tipo de modelos con variable dependiente binaria, una regresión lineal puede estimar probabilidades incorrectas, fuera del intervalo [0,1]. Por eso, se deben emplear modelos de probabilidad no lineal como el LOGIT o PROBIT.

# Predicción balance = 100
#-------------------------------
p1 <- predict(object = modelo_lineal, newdata = data.frame(balance = 200))
print(p1)
##           1 
## -0.04921752
# Predicción balance = 10000
#-------------------------------
p2 <- predict(object = modelo_lineal, newdata = data.frame(balance = 10000))
print(p2)
##       1 
## 1.22353

Al emplear funciones no lineales como una función sigmeoide, se estaría estimando una curva que nunca proporcionaría probabilidades fuera del rango [0,1].

# Regresión MCO lineal de default ~ Balance
# -----------------------------------------------
modelo_logistico = glm(default ~balance, data = datos, family = "binomial")
summary(modelo_logistico)
## 
## Call:
## glm(formula = default ~ balance, family = "binomial", data = datos)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.065e+01  3.612e-01  -29.49   <2e-16 ***
## balance      5.499e-03  2.204e-04   24.95   <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: 2920.6  on 9999  degrees of freedom
## Residual deviance: 1596.5  on 9998  degrees of freedom
## AIC: 1600.5
## 
## Number of Fisher Scoring iterations: 8
# Representación de los datos y el modelo
# -----------------------------------------
library(ggplot2)
ggplot(data = datos, aes(x = balance, y = default)) + 
  geom_point(aes(color = as.factor(default)), shape = 1) +
  stat_function(fun = function(x){predict(modelo_logistico,
                                          newdata = data.frame(balance=x),
                                          type = "response")})+
  
  theme_bw()+
  labs(title = "Regresión logística",
       y = "Probabilidad default") +
  theme(legend.position = "bottom")

# Con geom_smooth se puede obtener el gráfico directamente.
#--------------------------------------------------------
ggplot(data = datos, aes(x=balance, y=default)) +
  geom_point(aes(color=as.factor(default)), shape = 1) +
  geom_smooth(method = "glm",
              method.args = list(family = "binomial"),
              color = "grey20",
              se = TRUE) +
  theme_bw() +
  theme(legend.position = "bottom")
## `geom_smooth()` using formula = 'y ~ x'

EJERCICIO

Un estudio quiere establecer un modelo que permita calcular la probabilidad de obtener una matrícula de honor al final del bachillerato en función de la nota que se ha obtenido en matemáticas. La variable matrícula está codificada como 0 si no se tiene matrícula y 1 si se tiene.

matricula <- as.factor(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1,
                         0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1,
                         0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0,
                         0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
                         1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0,
                         1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1,
                         1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1,
                         0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
                         0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0,
                         0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 0,
                         0, 0, 0, 0, 1, 0, 0, 0, 1, 1))
matematicas <- c(41, 53, 54, 47, 57, 51, 42, 45, 54, 52, 51, 51, 71, 57, 50, 43,
                 51, 60, 62, 57, 35, 75, 45, 57, 45, 46, 66, 57, 49, 49, 57, 64,
                 63, 57, 50, 58, 75, 68, 44, 40, 41, 62, 57, 43, 48, 63, 39, 70,
                 63, 59, 61, 38, 61, 49, 73, 44, 42, 39, 55, 52, 45, 61, 39, 41,
                 50, 40, 60, 47, 59, 49, 46, 58, 71, 58, 46, 43, 54, 56, 46, 54,
                 57, 54, 71, 48, 40, 64, 51, 39, 40, 61, 66, 49, 65, 52, 46, 61,
                 72, 71, 40, 69, 64, 56, 49, 54, 53, 66, 67, 40, 46, 69, 40, 41,
                 57, 58, 57, 37, 55, 62, 64, 40, 50, 46, 53, 52, 45, 56, 45, 54,
                 56, 41, 54, 72, 56, 47, 49, 60, 54, 55, 33, 49, 43, 50, 52, 48,
                 58, 43, 41, 43, 46, 44, 43, 61, 40, 49, 56, 61, 50, 51, 42, 67,
                 53, 50, 51, 72, 48, 40, 53, 39, 63, 51, 45, 39, 42, 62, 44, 65,
                 63, 54, 45, 60, 49, 48, 57, 55, 66, 64, 55, 42, 56, 53, 41, 42,
                 53, 42, 60, 52, 38, 57, 58, 65)
datos <- data.frame(matricula, matematicas)
head(datos, 4)
##   matricula matematicas
## 1         0          41
## 2         0          53
## 3         0          54
## 4         0          47
str(datos)
## 'data.frame':    200 obs. of  2 variables:
##  $ matricula  : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ matematicas: num  41 53 54 47 57 51 42 45 54 52 ...
summary(datos)
##  matricula  matematicas   
##  0:151     Min.   :33.00  
##  1: 49     1st Qu.:45.00  
##            Median :52.00  
##            Mean   :52.65  
##            3rd Qu.:59.00  
##            Max.   :75.00
# Tabla de respuesta
# ---------------------------
library(ggplot2)
table(datos$matricula)
## 
##   0   1 
## 151  49
ggplot(data = datos, aes(x = matricula, y = matematicas, color = matricula)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.1)+
  theme_bw()+
  theme(legend.position = "bottom")

Parece existir una diferencia entre la nota de las personas con matrícula y sin matrícula.

Regresión Logística

# MODELO LOGIT
#-----------------------------
modelo <- glm(matricula ~matematicas, data = datos, family = "binomial")
summary(modelo)
## 
## Call:
## glm(formula = matricula ~ matematicas, family = "binomial", data = datos)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -9.79394    1.48174  -6.610 3.85e-11 ***
## matematicas  0.15634    0.02561   6.105 1.03e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 222.71  on 199  degrees of freedom
## Residual deviance: 167.07  on 198  degrees of freedom
## AIC: 171.07
## 
## Number of Fisher Scoring iterations: 5
Odd_beta0 <- exp(-9.79394)
Odd_beta0
## [1] 5.578866e-05
Odd_beta1 <- exp(0.15634)
Odd_beta1
## [1] 1.169224

Odd beta 0: La probabilidad de obtener matrícula sin tener en consideración la nota de matemáticas es:

 P(matricula)/P(no matricula) = 0.00005
 
 P(matricula) = 0.00005*P(no matricula)
 
 P(matricula) < P(no matricula)

Odd beta 1: La probablidad de obtener matrícula teniendo en cuenta o en consideración la nota de matematicas es:

 P(matricula)/P(no matricula) = 1.1692
 
 P(matricula)/P(no matricula) = 1.1692
 
 P(matricula) = 1.1692*P(no matricula)

Se obtiene que, por cada unidad que se incrementa la variable matemáticas, los odds de obtener matrícula se incrementan en promedio 1.169 unidades.

No hay que confundir que si se aumenta la nota de matemáticas en una unidad, lo demás constante, la probabilidad de obtener matrícula aumente en un 1.1692%

La nota media de matemáticas era: Mean: 52.65

z1 = -9.79394 + 0.15634*52.65
prob1 = 1/(1+exp(-z1))
prob1
## [1] 0.1732683
z2 = -9.79394 + 0.15634*(52.65+1)
prob2 = 1/(1+exp(-z2))
prob2
## [1] 0.1968185
diferencialP <- prob2-prob1
diferencialP
## [1] 0.02355016

Se estima que aumentar la nota de matemáticas en una unidad, partiendo del valor medio, en ceteris paribus, provocaría que la probabilidad de obtener una matrícula aumente en un 2.35%.

#intervalo de confianza
#---------------------------
confint(object= modelo, level = 0.95)
## Waiting for profiling to be done...
##                   2.5 %     97.5 %
## (Intercept) -12.9375208 -7.0938806
## matematicas   0.1093783  0.2103937
# Representación sin IC
#-------------------------------

# Codificación 0,1 de la variable respuesta
datos$matricula <- as.character(datos$matricula)
datos$matricula <- as.numeric(datos$matricula)


plot(matricula ~ matematicas, datos, col = "darkblue",
     main = "modelo regresión logística",
     ylab = "P(matrícula=1| matematicas)",
     xlab = "matemáticas", pch = "I")

# type = "response" devuelve las predicciones en forma de probabilidad en lugar de en log_ODDs
curve(predict(modelo, data.frame(matematicas = x), type = "response"),
      col = "firebrick", lwd = 2.5, add = TRUE)

# Gráfico de probabilidad con intervalo de confianza
# ---------------------------------------------------

datos$matricula <- as.character(datos$matricula)
datos$matricula <- as.numeric(datos$matricula)

# Se crea un vector con nuevos valores interpolados en el rango de observaciones.
nuevos_puntos <- seq(from = min(datos$matematicas), to = max(datos$matematicas),
                     by = 0.5)

# Predicciones de los nuevos puntos según el modelo.
predicciones <- predict(modelo, data.frame(matematicas = nuevos_puntos),
                        se.fit = TRUE)

# Mediante la función logit se transforman los log_ODDs a probabilidades.
predicciones_logit <- exp(predicciones$fit) / (1+exp(predicciones$fit))

# logODDs del intervalo se transforman en probabilidades con la función logit.
limite_inferior       <- predicciones$fit - 1.96 * predicciones$se.fit
limite_inferior_logit <- exp(limite_inferior) / (1 + exp(limite_inferior))
limite_superior       <- predicciones$fit + 1.96 * predicciones$se.fit
limite_superior_logit <- exp(limite_superior) / (1 + exp(limite_superior))

# Se crea un dataframe con los nuevos puntos y sus predicciones
datos_curva <- data.frame(matematicas = nuevos_puntos,
                          probabilidad_matricula = predicciones_logit,
                          limite_inferior_logit = limite_inferior_logit, 
                          limite_superior_logit = limite_superior_logit)

ggplot(datos, aes(x = matematicas, y = matricula)) +
      geom_point(aes(color = as.factor(matricula)), shape = "I", size = 3) + 
      geom_line(data = datos_curva, aes(y = probabilidad_matricula),
                color = "firebrick") + 
      geom_line(data = datos_curva, aes(y = limite_inferior_logit),
                linetype = "dashed") + 
      geom_line(data = datos_curva, aes(y = limite_superior_logit),
                linetype = "dashed") + 
      theme_bw() +
      labs(title = "Modelo regresión logística matrícula ~ nota matemáticas",
           y = "P(matrícula = 1 | matemáticas)", y = "matemáticas") + 
      theme(legend.position = "bottom") +
      theme(plot.title = element_text(size = 10))

Evaluación del modelo

A la hora de evaluar la validez y calidad de un modelo de regresión logística, se analiza tanto el modelo en su conjunto como los predictores que lo forman.

Se considera que el modelo es útil si es capaz de mostrar una mejora explicando las observaciones respecto al modelo nulo (sin predictores). El test Likelihood ratio calcula la significancia de la diferencia de residuos entre el modelo de interés y el modelo nulo. El estadístico sigue una distribución chi-cuadrado con grados de libertad equivalentes a la diferencia de grados de libertad de los dos modelos.

# Diferencia de residuos
# En R, un objeto glm almacena la "deviance" del modelo, así como la "deviance" del modelo nulo.
dif_residuos <- modelo$null.deviance - modelo$deviance

# Grados de libertad
df <- modelo$df.null - modelo$df.residual

# p-value
p_value <- pchisq(q = dif_residuos, df = df, lower.tail = FALSE)
paste("Diferencia de residuos:", round(dif_residuos, 4))
## [1] "Diferencia de residuos: 55.6368"
paste("Grados de libertad:", df)
## [1] "Grados de libertad: 1"
paste("p-value:", p_value)
## [1] "p-value: 8.71759108087094e-14"
# El mismo cálculo se puede obtener directamente con:
anova(modelo, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: matricula
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                          199     222.71              
## matematicas  1   55.637       198     167.07 8.718e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En este caso, al ser p-valor prácticamente cero, el modelo obtenido sí es significativo.

library(vcd)
## Loading required package: grid
## 
## Attaching package: 'vcd'
## The following object is masked from 'package:ISLR':
## 
##     Hitters
predicciones <- ifelse(test = modelo$fitted.values > 0.5, yes = 1, no = 0)
matriz_confusion <- table(modelo$model$matricula, predicciones,
                          dnn = c("observaciones", "predicciones"))
matriz_confusion
##              predicciones
## observaciones   0   1
##             0 140  11
##             1  27  22
mosaic(matriz_confusion, shade = T, colorize = T,
       gp = gpar(fill = matrix(c("green", "red", "red", "green"), 2, 2)))

# Fiabilidad
#-------------------------------------------
fiabilidad <- (140+22)/(140+22+11+27)
fiabilidad*100
## [1] 81

El modelo es capaz de predecir con la muestra de entrenamiento correctamente el 81% de los datos o de las observaciones.