Clase de 2022/01/05

Regresion Logistica en Machine Learning

La regresión logística es otra técnica que el machine learning ha tomado prestada del campo de la estadística. Es una forma poderosa de modelar un resultado binomial con una o más variables explicativas. Mide la relación entre la variable dependiente categórica y una o más variables independientes estimando probabilidades usando una función logística, que es la distribución logística acumulativa.

Regresión logística en R con glm

Hoy estudiaremos un ejemplo de una regresión logística binaria. Trabajaremos con el paquete ISLR, que proporcionará el conjunto de datos, y la función glm(), que generalmente usa el metodo de modelos lineales generalizados para ajustar el modelo de regresión logística.

¿Por qué usamos glm?

Dado a que es una generalización flexible de la regresión lineal ordinaria (que usa MCO), y permite variables de respuesta que tienen distribución de errores distintos de una distribución normal. Usar una regresión lineal en el caso de variables binarias, es decir, utilizando un modelo que pronostique la probabilidad de elegir entre sí/no (una variable de Bernoulli) es incluso menos apropiado, ya que las probabilidades están limitadas en ambos extremos (deben estar entre 0 y 1). Hay que imaginar, por ejemplo, un modelo que pronostique la probabilidad de que una persona determinada vaya a la playa en función de la temperatura. Un modelo razonable podría predecir, por ejemplo, que un cambio en 10 grados hace que una persona tenga el doble de probabilidades de ir o no ir a la playa. Pero, ¿qué significa “el doble” en términos de probabilidad? No puede significar literalmente duplicar el valor de probabilidad (por ejemplo, 50% se convierte en 100%, 75% se convierte en 150%, etc.). Más bien, es la razón de oportunidades la que se duplica: de una razón de oportunidades 2:1, a una razón de oportunidades 4:1, a una razón de oportunidades 8:1, etc. Tal modelo es un log-odds o un modelo logístico, que es lo que veremos hoy.

Librerías

Lo primero es cargar las librerias y usando el paquete ISLR, obtendremos la base de datos que vamos a usar. ISLR provee los datos usadas en el libro “An Introduction to Statistical Learning with Applications in R”.

#install.packages("ISLR", dependencies = TRUE)
#install.packages("Amelia", dependencies = TRUE)
#install.packages("mlbench", dependencies = TRUE)
#install.packages("corrplot", dependencies = TRUE)
require(ISLR) # Collection of data-sets used in the book 'An Introduction to Statistical Learning with Applications in R'.
## Loading required package: ISLR
library(Amelia) # A tool that "multiply imputes" missing data in a single cross-section (such as a survey), from a time series (like variables collected for each year in a country), or from a time-series-cross-sectional data set.
## Loading required package: Rcpp
## ## 
## ## Amelia II: Multiple Imputation
## ## (Version 1.8.0, built: 2021-05-26)
## ## Copyright (C) 2005-2022 James Honaker, Gary King and Matthew Blackwell
## ## Refer to http://gking.harvard.edu/amelia/ for more information
## ##
library(mlbench) # A collection of artificial and real-world machine learning benchmark problems, including, e.g., several data sets from the UCI repository.
library(corrplot) # Provides a visual exploratory tool on correlation matrix that supports automatic variable reordering to help detect hidden patterns among variables.
## corrplot 0.92 loaded
library(lattice) # A powerful and elegant high-level data visualization system inspired by Trellis graphics, with an emphasis on multivariate data.
library(caret) # Misc functions for training and plotting classification and regression models.
## Loading required package: ggplot2
#Si tienen problemas con Rcpp usen el siguiente codigo #Sys.setenv(PATH="%PATH%;C:/Rtools/gcc-4.6.3/bin;c:/Rtools/bin")

Hoy usaremos la base de datos Smarket, la cual contiene infomacion diaria de los rendimientos porcentuales del índice bursátil S&P 500 entre 2001 y 2005.

names(Smarket)
## [1] "Year"      "Lag1"      "Lag2"      "Lag3"      "Lag4"      "Lag5"     
## [7] "Volume"    "Today"     "Direction"

Explorando la data

head(Smarket)
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  
##            
##            
##            
##

Primeras visualizaciones

La visualización de datos es quizás la forma más rápida y útil de resumir y aprender más sobre los datos. Comenzará explorando las variables numéricas individualmente.

Los histogramas proporcionan un gráfico de barras de una variable numérica dividida en contenedores con la altura que muestra el número de instancias que caen en cada contenedor. Son útiles para obtener una indicación de la distribución de un atributo.

par(mfrow=c(1,8))
for(i in 1:8) {
    hist(Smarket[,i], main=names(Smarket)[i])
}

Es extremadamente difícil de ver, pero la mayoría de las variables muestran una distribución normal.

Otra forma de ver la distribución de los datos de una manera diferente es utilizando diagramas de caja y bigotes. El cuadro captura el 50% medio de los datos, la línea muestra la mediana y los bigotes de las gráficas muestran la extensión razonable de los datos. Cualquier punto fuera de los bigotes es un buen candidato para outliers.

par(mfrow=c(1,8))
for(i in 1:8) {
    boxplot(Smarket[,i], main=names(Smarket)[i])
}

Como vemos, los Lags y Today tienen un rango similar. De lo contrario, no hay señales de valores atípicos.

Los datos faltantes tienen un gran impacto en el modelado. Por lo tanto, vamos a usar un gráfico de missings para tener una idea rápida de la cantidad de datos faltantes. El eje x muestra atributos y el eje y muestra instancias. Las líneas horizontales indican datos faltantes para una instancia, los bloques verticales representan datos faltantes para un atributo.

missmap(Smarket, col=c("blue", "red"), legend=FALSE)

¡Tenemos suerte! ¡No faltan datos!

Comencemos a calcular la correlación entre cada par de variables numéricas. Estas correlaciones por pares se pueden representar en una gráfica de matriz de correlación para tener una idea de qué variables se mueven juntas.

correlations <- cor(Smarket[,1:8])
corrplot(correlations, method="circle")

Se utilizó una representación de puntos donde el azul representa la correlación positiva y el rojo negativo. Cuanto mayor sea el punto, mayor será la correlación. Puede ver que la matriz es simétrica y que las diagonales están perfectamente correlacionadas positivamente porque muestra la correlación de cada variable consigo misma. Desafortunadamente, ninguna de las variables está correlacionada entre sí.

Hagamos una gráfica de los datos. Hay una función pairs() que traza las variables de Smarket en una matriz de diagrama de dispersión. En este caso, Direction, su respuesta binaria, es el indicador de color:

pairs(Smarket, col=Smarket$Direction)

Parece que no hay mucha correlación aquí. La variable de clase se deriva de la variable Today, por lo que Up y Down parece hacer una división. Aparte de eso, no hay mucho que hacer.

Echemos un vistazo a la distribución de densidad de cada variable desglosada por valor de Direction. Al igual que la matriz de gráficos de dispersión anterior, el gráfico de densidad por Dirección puede ayudar a ver la separación de Up y Down. También puede ayudar a comprender la superposición en los valores de Direction para una variable.

x <- Smarket[,1:8]
y <- Smarket[,9]
scales <- list(x=list(relation="free"), y=list(relation="free"))
featurePlot(x=x, y=y, plot="density", scales=scales)

Puede ver que los valores de Direction se superponen para todas estas variables, lo que significa que es difícil predecir Up o Down en función de solo una o dos variables.

Construcción del modelo de regresión logística

Ahora llamamos a la función glm.fit(). El primer argumento que se entrega a esta función es una fórmula de R (lo que se encuentra dentro de los paréntesis). En este caso, la fórmula indica que Direction es la respuesta, mientras que las variables Lag y Volume son los predictores.

Sin embargo, en este caso, debemos dejar en claro que deseamos ajustar un modelo de regresión logística. Ésto se resuelve estableciendo el argumento family y luego indicando binomial. De esta manera, le dicimos a glm() que se ajuste a modelo de regresión logística en lugar de uno de los muchos otros modelos en los que podemos aplicar este algoritmo (glm).

# Logistics Regression
glm.fit <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume, data = Smarket, family = binomial)

A continuación, vamos a hacer un summary(), para entender como es el ajuste:

summary(glm.fit)
## 
## 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

Como podemos ver, summary() devuelve la estimación, los errores estándar, la puntuación z y los valores p de cada uno de los coeficientes. Parece que ninguno de los coeficientes es significativo aquí. También le da la desviación nula (la desviación solo para la media) y la desviación residual (la desviación para el modelo con todos los predictores). Hay una diferencia muy pequeña entre los 2, junto con 6 grados de libertad.

Ahora asignamos el resultado predict() de glm.fit() a glm.probs, agregando el argumento type y fijandolo como "response". Esto hará predicciones sobre los datos de entrenamiento que usamos para ajustar el modelo y entregará un vector de probabilidades ajustadas.

Veamos a las primeras 5 probabilidades, están muy cerca del 50%:

glm.probs <- predict(glm.fit,type = "response")
glm.probs[1:5]
##         1         2         3         4         5 
## 0.5070841 0.4814679 0.4811388 0.5152224 0.5107812

Ahora vamos a hacer una predicción de si el mercado subirá o bajará en función de los retrasos y otros predictores. En particular, convertiré las probabilidades en clasificaciones estableciendo un umbral en 0.5. Para hacer esto, utilizaremos un comando ifelse().

glm.pred <- ifelse(glm.probs > 0.5, "Up", "Down")

glm.pred es un vector de verdaderos y falsos (dato lógico o booleano). Si glm.probs es mayor que 0.5, glm.pred lo llamará "Up"; de lo contrario, lo llamará "False".

Aquí, adjuntamos el dataframe Smarket y creamos una tabla de glm.pred, que son los altibajos de Direction. También tomas la media de éstos.

attach(Smarket)
table(glm.pred,Direction)
##         Direction
## glm.pred Down  Up
##     Down  145 141
##     Up    457 507
# Mean
mean(glm.pred == Direction)
## [1] 0.5216

De la tabla, las instancias en las diagonales son donde obtiene la clasificación correcta, y fuera de las diagonales es donde se equivoca. Parece que cometimos muchos errores. La media da una proporción de 0,52 (accuracy).

Creación de muestras de prueba y entrenamiento

¿Cómo podemos mejorar el modelo? Dividiremos los datos en un conjunto de entrenamiento y un conjunto de prueba, esta estrategia es una practica adecuada.

# Make training and test set
train = Year<2005
glm.fit <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume,
               data = Smarket,
               family = binomial,
               subset = train)

glm.probs <- predict(glm.fit,
                    newdata = Smarket[!train,],
                    type = "response")

glm.pred <- ifelse(glm.probs > 0.5, "Up", "Down")

Veamos este fragmento de código en detalle:

  • train es igual al año menor que 2005. Para todo el año menor que 2005, obtendremos True; de lo contrario, obtendremos False. (¿considera usted esto una buena forma de dividir train y test?) Luego, reajustamos el modelo con glm.fit (), excepto que el subconjunto es igual a 'Train', lo que significa que solo se ajusta a los datos del año anterior a 2005. Luego usamos la función predict() nuevamente para glm.probs para predecir los datos restantes en un año mayor o igual a 2005. Para los nuevos datos, usamos Smarket indexado por !Train (!Train es True si el año es mayor o igual a 2005). Establecemos type igual a "response" para predecir las probabilidades. *Finalmente, usamos la función ifelse() nuevamente con glm.pred para predecir las variables como Up y Down (¿de lo contrario, cual sería el resultado esperado?).

Ahora creamos una nueva variable para almacenar un nuevo subconjunto Test y lo llama Direction.2005. La variable de respuesta sigue siendo Direction. Hacemos una tabla y calculamos la media en este nuevo conjunto de prueba:

Direction.2005 = Smarket$Direction[!train]
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

Ja, lo hizo peor que en el caso anterior. ¿Cómo pudo pasar esto?

Resolver el overfitting (sobreajuste)

Bueno, es posible que hayamos sobreajustado los datos. Para solucionar este problema, ajustaremos un modelo más pequeño y utilizaremos Lag1, Lag2, Lag3 como predictores, dejando de lado todas las demás variables. El resto del código es el mismo.

# Fit a smaller model
glm.fit = glm(Direction ~ Lag1 + Lag2 + Lag3, data = Smarket, family = binomial, subset = train)
glm.probs = predict(glm.fit, newdata = Smarket[!train,], type = "response")
glm.pred = ifelse(glm.probs > 0.5, "Up", "Down")
table(glm.pred, Direction.2005)
##         Direction.2005
## glm.pred Down  Up
##     Down   39  31
##     Up     72 110
mean(glm.pred == Direction.2005)
## [1] 0.5912698

Bueno, obtuvimos una tasa de clasificación del 59%, no tan mal. El uso del modelo más pequeño parece funcionar mejor.

Por último, haremos un summary() a glm.fit para ver si hay cambios significativos.

summary(glm.fit)
## 
## Call:
## glm(formula = Direction ~ Lag1 + Lag2 + Lag3, family = binomial, 
##     data = Smarket, subset = train)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.338  -1.189   1.072   1.163   1.335  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.032230   0.063377   0.509    0.611
## Lag1        -0.055523   0.051709  -1.074    0.283
## Lag2        -0.044300   0.051674  -0.857    0.391
## Lag3         0.008815   0.051495   0.171    0.864
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1383.3  on 997  degrees of freedom
## Residual deviance: 1381.4  on 994  degrees of freedom
## AIC: 1389.4
## 
## Number of Fisher Scoring iterations: 3

Nada se volvió significativo, al menos los valores P son mejores, lo que indica una mejora en la predicción.