Capítulo 9 Logit

En la base de datos Default,²⁵ la variable default pertenece a una de dos categorías, Yes o No. En vez de modelar esta respuesta \(Y\) directamente, la regresión logística modela la probabilidad de que \(Y\) pertenezca a una categoría. Por lo tanto, si se busca estimar la probabilidad de default, se puede transformar la variable original en una variable binaria \(1\) (Yes); \(0\) (No).

En general, la estrategia consiste en:

\(Y\) variable binaria \(0/1\).

\(X\) vector de \(K\) predictores.

Se debe construir un modelo para:

\(p = Pr(Y = 1 \mid X)\)

Como vimos en la subsección 6.6.3.1 el clasificador de Bayes:

\(\hat{Y} = 1\) si \(p > 0,5\)
\(\hat{Y} = 0\) si \(p \le 0,5\)

Recordar que en el caso particular de la probabilidad de default puede resultar relativamente más costoso para un banco clasificar a un deudor malo en la categoría No default que a un deudor bueno como Si default (pérdida asimétrica). Por lo tanto, podría elegirse un umbral más bajo (ej. \(p > 0,3\)) para la categoría Si default.

9.1 Modelo logit

\[\begin{equation} \tag{9.1} p = \frac{e^z}{1+e^z} \end{equation}\]

donde \(z \equiv X \beta\) con \(\beta\) un vector de \(K\) coeficientes.

El panel izquierdo de la Figura 9.1 muestra la probabilidad estimada usando regresión lineal donde algunos valores de probabilidad son negativos (también podría haber valores mayores a \(1\) para valores elevados de la variable balance, ambos extremos inconsistentes con la noción de probabilidad). Los puntos naranja indican los valores \(0/1\) codificados por defecto (No o Si). El panel derecho exhibe la probabilidad de incumplimiento estimada mediante regresión logística donde todos los valores \(\in (0,1)\).

Figura 9.1: Regresión lineal vs. logística

9.1.1 Interpretacion de coeficientes en el modelo logit

Dada la forma funcional del modelo logit, primero se introduce el concepto de odds ratio (\(\frac{p}{1 - p}\)) para luego derivar la interpretación de un coeficiente. El se define odds ratio como la probabilidad de que un evento suceda en relación a que el mismo evento no suceda. Por ejemplo, en promedio \(1\) de cada \(5\) personas con un odds de \(1/4\) no pagará su deuda, ya que la probabilidad de default \(p(X) = 0,2\) implica un odds de \(\frac{0,2}{1 − 0,2} = 0,25 = 1/4\).

En el modelo logit:

\[\begin{equation} \tag{9.2} (1 - p) = \frac{1}{1+e^z} \end{equation}\]

Entonces el odds ratio:

\[\begin{equation} \tag{9.3} \frac{p}{1 - p} = e^z \end{equation}\]

Si se toma el logaritmo del odds ratio queda:

\[\begin{equation} \tag{9.4} ln(\frac{p}{1 - p}) = z = X \beta \end{equation}\]

El lado izquierdo se llama log odds o logit. El modelo de regresión logística (9.1) tiene un logit que es lineal en \(X\). Por lo tanto, el \(k-ésimo\) coeficiente es:

\[\begin{equation} \tag{9.5} \beta_k = \frac{\partial ln(\frac{p}{1 - p})}{\partial X_k} \end{equation}\]

Aumentar \(X_k\) en una unidad cambia el log odds en \(\beta_k\). Dado que la relación entre \(p(X)\) y \(X\) en (9.1) no es lineal, \(\beta_k\) no corresponde al cambio en \(p(X)\) asociado con el aumento de una unidad en \(X\). La cantidad que \(p(X)\) cambia debido a un cambio de una unidad en \(X\) depende del valor actual de \(X\). Pero independientemente del valor de \(X\), si \(\beta_k\) es positivo, entonces el aumento de \(X\) se asocia con el aumento de \(p(X)\), y si \(\beta_k\) es negativo, el aumento de \(X\) se asocia con la disminución de \(p(X)\).²⁶

En este caso, \(\hat{{\beta}}\) se estima por máxima verosimilitud (no es una forma cerrada como \(MCO\)).

\[\begin{equation} \tag{9.6} \hat{\theta} = argmax_{\theta} \prod_{i = 1}^{n} L(\theta; y_i) \end{equation}\]

Luego se reemplazan los coeficientes estimados en (9.7) para obtener las probabilidades estimadas.

\[\begin{equation} \tag{9.7} \hat{p_i} = \frac{e^{X_i \hat{\beta}}}{1+e^{X_i \hat{\beta}}} \end{equation}\]

9.2 Aplicacion practica

Comenzaremos examinando algunas estadísticas y gráficos de los datos Smarket, que forman parte de la biblioteca ISLR2. Esta base de datos consiste en rendimientos porcentuales para el índice bursátil \(S\&P 500\) para más de \(1.250\) días, desde principios de \(2001\) hasta finales de \(2005\). Para cada fecha, registra los rendimientos porcentuales de los cinco días hábiles anteriores, de Lag1 a Lag5. También registra el volume (el número de acciones negociadas el día anterior, en miles de millones), Today (el rendimiento porcentual en la fecha en cuestión) y direction (si el mercado estaba Up o Down en este día). El objetivo es predecir la “dirección” (una respuesta cualitativa) usando las otras características.

library(ISLR2)
names(Smarket)

## [1] "Year"      "Lag1"      "Lag2"      "Lag3"      "Lag4"      "Lag5"     
## [7] "Volume"    "Today"     "Direction"

dim(Smarket)

## [1] 1250    9

summary(Smarket)

##       Year           Lag1                Lag2                Lag3          
##  Min.   :2001   Min.   :-4.922000   Min.   :-4.922000   Min.   :-4.922000  
##  1st Qu.:2002   1st Qu.:-0.639500   1st Qu.:-0.639500   1st Qu.:-0.640000  
##  Median :2003   Median : 0.039000   Median : 0.039000   Median : 0.038500  
##  Mean   :2003   Mean   : 0.003834   Mean   : 0.003919   Mean   : 0.001716  
##  3rd Qu.:2004   3rd Qu.: 0.596750   3rd Qu.: 0.596750   3rd Qu.: 0.596750  
##  Max.   :2005   Max.   : 5.733000   Max.   : 5.733000   Max.   : 5.733000  
##       Lag4                Lag5              Volume           Today          
##  Min.   :-4.922000   Min.   :-4.92200   Min.   :0.3561   Min.   :-4.922000  
##  1st Qu.:-0.640000   1st Qu.:-0.64000   1st Qu.:1.2574   1st Qu.:-0.639500  
##  Median : 0.038500   Median : 0.03850   Median :1.4229   Median : 0.038500  
##  Mean   : 0.001636   Mean   : 0.00561   Mean   :1.4783   Mean   : 0.003138  
##  3rd Qu.: 0.596750   3rd Qu.: 0.59700   3rd Qu.:1.6417   3rd Qu.: 0.596750  
##  Max.   : 5.733000   Max.   : 5.73300   Max.   :3.1525   Max.   : 5.733000  
##  Direction 
##  Down:602  
##  Up  :648  
##            
##            
##            
##

pairs(Smarket)

La función cor() produce una matriz que contiene todas las correlaciones por pares entre los predictores.

cor(Smarket[, -9])

##              Year         Lag1         Lag2         Lag3         Lag4
## Year   1.00000000  0.029699649  0.030596422  0.033194581  0.035688718
## Lag1   0.02969965  1.000000000 -0.026294328 -0.010803402 -0.002985911
## Lag2   0.03059642 -0.026294328  1.000000000 -0.025896670 -0.010853533
## Lag3   0.03319458 -0.010803402 -0.025896670  1.000000000 -0.024051036
## Lag4   0.03568872 -0.002985911 -0.010853533 -0.024051036  1.000000000
## Lag5   0.02978799 -0.005674606 -0.003557949 -0.018808338 -0.027083641
## Volume 0.53900647  0.040909908 -0.043383215 -0.041823686 -0.048414246
## Today  0.03009523 -0.026155045 -0.010250033 -0.002447647 -0.006899527
##                Lag5      Volume        Today
## Year    0.029787995  0.53900647  0.030095229
## Lag1   -0.005674606  0.04090991 -0.026155045
## Lag2   -0.003557949 -0.04338321 -0.010250033
## Lag3   -0.018808338 -0.04182369 -0.002447647
## Lag4   -0.027083641 -0.04841425 -0.006899527
## Lag5    1.000000000 -0.02200231 -0.034860083
## Volume -0.022002315  1.00000000  0.014591823
## Today  -0.034860083  0.01459182  1.000000000

Como era de esperar, las correlaciones entre las variables rezagadas y los rendimientos de hoy son cercanas a cero. En otras palabras, parece haber poca correlación entre los rendimientos de hoy y los rendimientos de días anteriores. La única correlación sustancial es entre Year y volume. Al graficar los datos, que están ordenados cronológicamente, vemos que el volume aumenta con el tiempo. En otras palabras, el número promedio de acciones negociadas diariamente aumentó de 2001 a 2005.

attach(Smarket)
plot(Volume)

A continuación, ajustaremos un modelo de regresión logística para predecir la direction utilizando Lag1 hasta Lag5 y volume. La función glm() se puede utilizar para muchos tipos de modelos lineales generalizados, incluida la regresión logística. La sintaxis de la función glm() es similar a la de lm(), excepto que debemos utilizar el argumento family = binomial para decirle a R que ejecute una regresión logística.

glm.fits <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume,
    data = Smarket, family = binomial)
summary(glm.fits)

## 
## Call:
## glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + 
##     Volume, family = binomial, data = Smarket)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.446  -1.203   1.065   1.145   1.326  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.126000   0.240736  -0.523    0.601
## Lag1        -0.073074   0.050167  -1.457    0.145
## Lag2        -0.042301   0.050086  -0.845    0.398
## Lag3         0.011085   0.049939   0.222    0.824
## Lag4         0.009359   0.049974   0.187    0.851
## Lag5         0.010313   0.049511   0.208    0.835
## Volume       0.135441   0.158360   0.855    0.392
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1731.2  on 1249  degrees of freedom
## Residual deviance: 1727.6  on 1243  degrees of freedom
## AIC: 1741.6
## 
## Number of Fisher Scoring iterations: 3

El valor \(p\) más pequeño está asociado con Lag1. El coeficiente negativo sugiere que si el mercado tuvo un rendimiento positivo ayer, es menos probable que suba hoy. Sin embargo, con un valor de \(0,145\), el valor de \(p\) sigue siendo relativamente grande, por lo que no hay pruebas claras de una asociación real entre Lag1 y direction.

La función predict() se puede utilizar para predecir la probabilidad de que el mercado suba, dados los valores de los predictores. La opción type = "response" le dice a R que genere probabilidades de la forma \(P(Y=1|X)\). Si no se proporciona ninguna base de datos a la función predict(), se calculan las probabilidades para los datos de entrenamiento que se usaron para ajustar el modelo de regresión logística. Se imprimen las primeras diez probabilidades estimadas. Sabemos que estos valores corresponden a la probabilidad de que el mercado suba, en lugar de que baje, porque la función contrasts() indica que R ha creado una variable ficticia con un 1 para Up.

glm.probs <- predict(glm.fits, type = "response")
glm.probs[1:10]

##         1         2         3         4         5         6         7         8 
## 0.5070841 0.4814679 0.4811388 0.5152224 0.5107812 0.5069565 0.4926509 0.5092292 
##         9        10 
## 0.5176135 0.4888378

contrasts(Direction)

##      Up
## Down  0
## Up    1

Para hacer una predicción sobre si el mercado subirá o bajará en un día en particular, debemos convertir estas probabilidades pronosticadas en etiquetas de clase, Up o Down. Los siguientes dos comandos crean un vector de predicciones de clase basado en si la probabilidad prevista de un aumento del mercado es mayor o menor que \(0,5\).

glm.pred <- rep('Down', 1250)
glm.pred[glm.probs > .5] = 'Up'

El primer comando crea un vector de \(1.250\) elementos Down. La segunda línea transforma en Up todos los elementos para los que la probabilidad predicha de un aumento del mercado supera los \(0,5\). Dadas estas predicciones, la función table() se puede usar para producir una matriz de confusión (subsección 6.6.4) para determinar cuántas observaciones se clasificaron correcta o incorrectamente. Al ingresar dos vectores cualitativos, R creará una tabla de \(2 \times 2\) con recuentos del número de veces que ocurrió cada combinación. Por ejemplo, pronosticó Up y el mercado aumentó, predijo Up y el mercado disminuyó, etc.

Los elementos de la diagonal de la matriz de confusión indican predicciones correctas (accuracy), mientras que los fuera de la diagonal representan predicciones incorrectas. Por lo tanto, el modelo predijo correctamente que el mercado subiría en \(507\) días y que bajaría en \(145\) días, para un total de \(507+145 = 652\) predicciones correctas. La función mean() se puede utilizar para calcular la fracción de días en los que la predicción fue correcta. En este caso, la regresión logística predijo correctamente el movimiento del mercado \(52,2\%\) del tiempo.

table(glm.pred, Direction)

##         Direction
## glm.pred Down  Up
##     Down  145 141
##     Up    457 507

(507 + 145) / 1250

## [1] 0.5216

mean(glm.pred == Direction)

## [1] 0.5216

La función confusionMatrix() de la librería caret calcula la matriz de confusión con estadísticas completas.

library('caret')
confusionMatrix(factor(glm.pred), factor(Direction), positive = 'Up')

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Down  Up
##       Down  145 141
##       Up    457 507
##                                           
##                Accuracy : 0.5216          
##                  95% CI : (0.4935, 0.5496)
##     No Information Rate : 0.5184          
##     P-Value [Acc > NIR] : 0.4216          
##                                           
##                   Kappa : 0.0237          
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.7824          
##             Specificity : 0.2409          
##          Pos Pred Value : 0.5259          
##          Neg Pred Value : 0.5070          
##              Prevalence : 0.5184          
##          Detection Rate : 0.4056          
##    Detection Prevalence : 0.7712          
##       Balanced Accuracy : 0.5116          
##                                           
##        'Positive' Class : Up              
##

A primera vista, parece que el modelo de regresión logística funciona un poco mejor que un modelo aleatorio. Sin embargo, este resultado es engañoso porque entrenamos y testeamos el modelo en el mismo conjunto de observaciones de \(1.250\). En otras palabras, \(100\%-52.2\%=47.8\%\), es la tasa de error de entrenamiento. Como hemos visto anteriormente, la tasa de error de entrenamiento suele ser demasiado optimista: tiende a subestimar la tasa de error de test.

Para evaluar mejor la precisión del modelo de regresión logística, podemos ajustar el modelo usando parte de los datos y luego examinar qué tan bien predice los datos retenidos. Esto producirá una tasa de error más realista, en el sentido de que en la práctica estamos interesados en el rendimiento de nuestro modelo, no en los datos que usamos para ajustar el modelo, sino en los días futuros para los que se desconocen los movimientos del mercado.

Para implementar esta estrategia, primero crearemos un vector correspondiente a las observaciones de \(2001\) a \(2004\). Luego usaremos este vector para crear una base de datos retenidos de observaciones de \(2005\).

train <- (Year < 2005)
Smarket.2005 <- Smarket[!train, ]
dim(Smarket.2005)

## [1] 252   9

Direction.2005 <- Direction[!train]

El objeto train es un vector de \(1.250\) elementos, correspondientes a las observaciones de la base de datos. Los elementos del vector que corresponden a observaciones que ocurrieron antes de \(2005\) se establecen en VERDADERO, mientras que los que corresponden a observaciones en \(2005\) se establecen en FALSO. El objeto train es un vector booleano, ya que sus elementos son VERDADERO y FALSO. Los vectores booleanos se pueden utilizar para obtener un subconjunto de filas o columnas de una matriz. Por ejemplo, el comando Smarket[train, ] elegirá una submatriz de la base de datos del mercado de valores, correspondiente solo a las fechas anteriores a \(2005\), ya que son aquellos para los que los elementos de train son VERDADERO. El símbolo ! (negación) se puede utilizar para invertir todos los elementos de un vector booleano. Es decir, !train es un vector similar a train, excepto que los elementos que son VERDADERO en train se cambian a FALSO en !train, y los elementos que son FALSO en train se cambia a VERDADERO en !train. Por lo tanto, Smarket[!train, ] produce una submatriz de los datos del mercado de valores que contiene solo las observaciones para las cuales train es FALSO, es decir, las observaciones con fechas en \(2005\) (\(252\) observaciones).

Ahora ajustamos un modelo de regresión logística usando solo el subconjunto de las observaciones que corresponden a fechas anteriores a \(2005\), usando el argumento subset. Luego obtenemos las probabilidades pronosticadas de que el mercado de valores suba para cada uno de los días de la base de test, es decir, para los días de \(2005\).

glm.fits <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume,
    data = Smarket, family = binomial, subset = train)

glm.probs <- predict(glm.fits, Smarket.2005,  type = "response")

Tener en cuenta que hemos entrenado y testeado el modelo en dos bases de datos completamente separadas: el entrenamiento se realizó solo con las fechas anteriores a \(2005\) y el test se realizó solo con las fechas de \(2005\). Finalmente, calculamos las predicciones para \(2005\) y las comparamos con los movimientos reales del mercado durante ese período de tiempo.

glm.pred <- rep("Down", 252)
glm.pred[glm.probs > .5] <- "Up"
table(glm.pred, Direction.2005)

##         Direction.2005
## glm.pred Down Up
##     Down   77 97
##     Up     34 44

mean(glm.pred == Direction.2005)

## [1] 0.4801587

mean(glm.pred != Direction.2005)

## [1] 0.5198413

La notación != significa que no es igual a (o distinto), por lo que el último comando calcula la tasa de error en los datos de test. Los resultados son bastante decepcionantes: la tasa de error de test es de \(52 \%\). Por supuesto, este resultado no es tan sorprendente, dado que, en general, no se esperaría poder utilizar los rendimientos de días anteriores para predecir el rendimiento futuro del mercado.

Bibliografia

Wooldridge, Jeffrey. 2012. Introductory Econometrics: A Modern Approach. Vol. 5th edition. South-Western College Publishing.

Esta base pertenece a la librería ISLR2.↩︎
Para la interpretación de coeficientes en la práctica usualmente se analizan los efectos marginales. Ver (Wooldridge 2012) capítulo \(17\).↩︎