Logistic Regression

En R la regresión logística está incluida dentro de la función glm() que engloba distintos tipos de generalized linear models. Para especificar que se quiere una regresión logística se tiene que indicar en el argumento family=binomial.

En primer lugar se prueba a incluir todas las variables disponibles como predictores excepto año.

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

## 
## 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

Para conocer como se ha codificado la variable respuesta (Direction) se recurre a la función contrast().

contrasts(Smarket$Direction)

##      Up
## Down  0
## Up    1

El modelo ha considerado como 1 el nivel Up y cómo 0 el nivel Down.

De entre los predictores introducidos, el que menor p-value tiene es Lag1 (0.145). Dado que su coeficiente estimado es negativo, si el mercado subió el día anterior, es menos probable que suba hoy. Sin embargo, el p-value es considerablemente alto por lo que no hay una clara evidencia de que exista relación entre Lag1 y Direction.

Una vez generado el modelo se tiene que evaluar cómo de bueno es prediciendo los resultados. Para ello se tienen en cuenta dos tipos de error:

• Trainig error: % de errores que comete el modelo al intentar predecir las observaciones que se han utilizado para generarlo.

• Test error: % de errores que comete el modelo al intentar predecir nuevas observaciones que no se han empleado en la creación del modelo. Este es el % de error que más interesa y que por lo generar siempre es mayor que el Training Error.

La función predict() de R permite predecir la probabilidad de que la variable respuesta pertenezca al nivel de referencia (en este caso 1 = Up), dados unos valores del predictor o predictores. Por defecto devuelve el valor Log of ODDs ya que es esta la variable que se modela en la regresión logística. Si se emplea el argumento type=response se devuelve directamente el valor de la probabilidad. Si a la función predict() no se le pasa ningún dataframe como argumento, devuelve directamente las predicciones para los datos empleados en su creación (training data).

predicciones <- predict(object = modelo_logistico, type = "response")
head(predicciones)

##         1         2         3         4         5         6 
## 0.5070841 0.4814679 0.4811388 0.5152224 0.5107812 0.5069565

El siguiente paso una vez predichas las probabilidades es utilizarlas para determinar a qué clase de la variable cualitativa pertenece cada observación. En este caso se considera como límite de decisión p=0.5. Si p>0.5 la observación se clasifica en el nivel 1(Up) de lo contrario se clasifica en el nivel 0 (Down).

prediccion <- data.frame(probabilidad = predicciones,
                         clase = rep(NA,length(predicciones)))
prediccion[prediccion$probabilidad < 0.5,"clase"] <- "Down"
prediccion[prediccion$probabilidad > 0.5,"clase"] <- "Up"
head(prediccion)

Mediante la función table() se puede obtener una matriz que muestre los aciertos, falsos positivos y falsos negativos:

table(clase_predicha = prediccion$clase, clase_real = Smarket$Direction)

##               clase_real
## clase_predicha Down  Up
##           Down  145 141
##           Up    457 507

paste( "% de acierto:", mean(prediccion$clase == Smarket$Direction))

## [1] "% de acierto: 0.5216"

paste( "% de error:",mean(prediccion$clase != Smarket$Direction))

## [1] "% de error: 0.4784"

El modelo logístico ha sido capaz de predecir correctamente el 52.2% de las observaciones, solo un poco mejor de lo que cabría esperar por azar (50%). El training error es de 47.8%. A efectos prácticos, es mucho más interesante conocer el test error. Un modo de hacerlo consiste en emplear solo un subconjunto de los datos disponibles para crear el modelo y con el resto de datos estudiar la precisión de las predicciones. En este caso se emplean las observaciones pertenecientes al intervalo 2001-2004 a modo de traning data y se evaluará el error de predicción con las observaciones del 2005.

train_data <- Smarket[Smarket$Year < 2005,]
test_data <- Smarket[!(Smarket$Year < 2005),]

# Se crea el modelo de regresión logística
modelo <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume ,
              data = train_data, family = "binomial")

# Se realizan las predicciones para el set de datos no empleado en la creación
# del modelo
predicciones <- predict(object = modelo,newdata = test_data, type = "response")
head(predicciones)

##       999      1000      1001      1002      1003      1004 
## 0.5282195 0.5156688 0.5226521 0.5138543 0.4983345 0.5010912

# Se considera como threshold de clasificación el 0.5
predicciones[predicciones > 0.5] <- "Up"
predicciones[predicciones != "Up"] <- "Down"

Finalmente se genera una tabla cruzada para identificar las predicciones correctas y erróneas.

table(clase_predicha = predicciones, clase_real = test_data$Direction)

##               clase_real
## clase_predicha Down Up
##           Down   77 97
##           Up     34 44

# % de aciertos
paste("% de acierto:", mean(predicciones == test_data$Direction))

## [1] "% de acierto: 0.48015873015873"

# % de errores
paste("% de error:", mean(predicciones != test_data$Direction))

## [1] "% de error: 0.51984126984127"

El % de aciertos empleando un nuevo conjunto de datos es del 47% y el test error del 52%. El modelo no se puede considerar útil ya que el error es superior al esperado por azar (50%). Esto no es de extrañar ya que el p-value de los predictores introducidos en el modelo es muy alto, indicando que no hay apenas relación entre ellos y la variable dependiente.

Una posible forma de mejorar el modelo es incluir únicamente los predictores que tienen menor p-value: Lag1, Lag2

# Se crea el modelo de regresión logística
modelo <- glm(Direction ~ Lag1 + Lag2, data = train_data, family = "binomial")

# Se realizan las predicciones para el set de datos no empleado en la creación
# del modelo
predicciones <- predict(object = modelo, newdata = test_data, type = "response")
head(predicciones)

##       999      1000      1001      1002      1003      1004 
## 0.5098275 0.5208237 0.5332635 0.5260574 0.5072103 0.5061388

# Se considera como threshold de clasificación el 0.5
predicciones[predicciones > 0.5]   <- "Up"
predicciones[predicciones != "Up"] <- "Down"

table(clase_predicha = predicciones, clase_real = test_data$Direction)

##               clase_real
## clase_predicha Down  Up
##           Down   35  35
##           Up     76 106

# % de aciertos
paste("% de acierto:", mean(predicciones == test_data$Direction))

## [1] "% de acierto: 0.55952380952381"

# % de errores
paste("% de error:", mean(predicciones != test_data$Direction))

## [1] "% de error: 0.44047619047619"

Se ha conseguido una pequeña mejora, esta vez el modelo es capaz de predecir correctamente la dirección del mercado en un 56% de las ocasiones, pasando a ser ligeramente mejor que lo esperado por azar. Aun así lo recomendable sería probarlo con nuevos datos para confirmarlo.

Además de estudiar el test error global, es interesante estudiar cómo se reparte este error entre falsos negativos y falsos positivos. Puede ocurrir que un modelo sea mucho mejor prediciendo en una dirección que en otra. Para el modelo aquí generado, la predicción de subida (Up) tiene un porcentaje de acierto de 106/(106+76)=58 y la de bajada de 35/(35+35)=50. Por lo tanto, es más seguro apostar por los días en los que el modelo indica subida que por los que indica bajada.

test_data$prediccion <- (predicciones)
test_data$acierto <- ifelse(test = test_data$Direction == test_data$prediccion,
                            yes = "Si", no = "No")
ggplot(data = test_data, aes(x = Lag1, y = Lag2, color = Direction,
                             shape = acierto)) +
  geom_point() +
  scale_shape_manual(values = c("No" = 4, "Si" = 19 )) +
  theme_bw()

Linear Discriminant Analysis (LDA)

El método LDA se puede emplear para variables cualitativas con dos o más niveles. En R se crean los modelos LDA con la función lda() del paquete MASS.

library(MASS)
modelo_lda <- lda(Direction ~ Lag1 + Lag2, data = train_data)
modelo_lda

## Call:
## lda(Direction ~ Lag1 + Lag2, data = train_data)
## 
## Prior probabilities of groups:
##     Down       Up 
## 0.491984 0.508016 
## 
## Group means:
##             Lag1        Lag2
## Down  0.04279022  0.03389409
## Up   -0.03954635 -0.03132544
## 
## Coefficients of linear discriminants:
##             LD1
## Lag1 -0.6420190
## Lag2 -0.5135293

Las prior probabilities of groups se estiman a partir de la proporción de observaciones de cada clase presentes en el set de datos empleado para crear el modelo. La media de los grupos se corresponden con la media de cada predictor dentro de cada una de las clases (en este caso Up y Down) y se emplean como estimación de μk del modelo. Se observa una ligera tendencia a que en promedio ambos predictores sean negativos cuando el mercado sube (Up) y positivo cuando baja (Down).

Se calculan las predicciones para los datos del año 2005:

predicciones_lda <- predict(object = modelo_lda, test_data)
#La función predict() aplicada a un modelo LDA devuelve una lista con 3 objetos.
# class: contiene la clasificación predicha por el modelo.
# posterior: las probabilidades posteriores de que la observación pertenezca a
# cada clase. Se asigna la observación a la clase con mayor probabilidad 
# posterior.
head(predicciones_lda$class, n = 3)

## [1] Up Up Up
## Levels: Down Up

head(predicciones_lda$posterior, n = 3)

##           Down        Up
## 999  0.4901792 0.5098208
## 1000 0.4792185 0.5207815
## 1001 0.4668185 0.5331815

Evaluación del modelo:

table(clase_predicha = predicciones_lda$class,
      clase_real =  test_data$Direction)

##               clase_real
## clase_predicha Down  Up
##           Down   35  35
##           Up     76 106

# % de aciertos
paste("% de acierto:", mean(predicciones_lda$class == test_data$Direction))

## [1] "% de acierto: 0.55952380952381"

# % de errores
paste("% de error:", mean(predicciones_lda$class != test_data$Direction))

## [1] "% de error: 0.44047619047619"

El modelo tiene un test error del 44%, es decir, que acierta en un 56% de las ocasiones. Está ligeramente por encima de los aciertos esperados por azar. Los resultados son muy similares a los obtenidos mediante regresión logística.

test_data$prediccion <- predicciones_lda$class
test_data$acierto <- ifelse(test = test_data$Direction == test_data$prediccion,
                            yes = "Si", no = "No")
ggplot(data = test_data, aes(x = Lag1, y = Lag2, color = Direction,
                             shape = acierto)) +
  geom_point() +
  scale_shape_manual(values = c("No" = 4, "Si" = 19 )) + 
  theme_bw()

Quadratic Discriminant Analysis (QDA)

El método QDA se puede emplear para variables cualitativas con dos o más niveles. En R se crean los modelos QDA con la función qda() del paquete MASS. La sintaxis y el tipo de output es igual que el de la función lda() descrita anteriormente.

library(MASS)
modelo_qda <- qda(Direction ~ Lag1 + Lag2, data = train_data)
modelo_qda

## Call:
## qda(Direction ~ Lag1 + Lag2, data = train_data)
## 
## Prior probabilities of groups:
##     Down       Up 
## 0.491984 0.508016 
## 
## Group means:
##             Lag1        Lag2
## Down  0.04279022  0.03389409
## Up   -0.03954635 -0.03132544

Predicciones con nuevos datos (2005):

predicciones_qda <- predict(object = modelo_qda, test_data)
head(predicciones_qda$class)

## [1] Up Up Up Up Up Up
## Levels: Down Up

head(predicciones_qda$posterior)

##           Down        Up
## 999  0.4873243 0.5126757
## 1000 0.4759011 0.5240989
## 1001 0.4636911 0.5363089
## 1002 0.4739253 0.5260747
## 1003 0.4903426 0.5096574
## 1004 0.4913561 0.5086439

Evaluación del modelo:

table(clase_predicha = predicciones_qda$class,
      clase_real =  test_data$Direction)

##               clase_real
## clase_predicha Down  Up
##           Down   30  20
##           Up     81 121

# % de aciertos
paste("% de acierto:", mean(predicciones_qda$class == test_data$Direction))

## [1] "% de acierto: 0.599206349206349"

# % de errores
paste("% de error:", mean(predicciones_qda$class != test_data$Direction))

## [1] "% de error: 0.400793650793651"

El modelo obtenido mediante quadratic discriminant analysis tiene un test error del 40%, es capaz de predecir correctamente casi un 60% de las veces. Este modelo es superior a los obtenidos por LDA y regresión logística, indicando que el modelo cuadrático captura mejor la verdadera relación entre predictores y variable dependiente que los modelos lineales.

test_data$prediccion <- predicciones_qda$class
test_data$acierto <- ifelse(test = test_data$Direction == test_data$prediccion,
                            yes = "Si", no = "No")
ggplot(data = test_data, aes(x = Lag1, y = Lag2, color = Direction,
                             shape = acierto)) +
  geom_point() +
  scale_shape_manual(values = c("No" = 4, "Si" = 19 )) + 
  theme_bw()

K-Nearest-Neighbors

En R, el método de predicción no paramétrico K-NN se aplica mediante la función knn() del paquete class. Esta función tiene un funcionamiento distinto a las dos anteriores. En lugar de generar un modelo que después se interroga con nuevos valores de los predictores, en un único paso se evalúa los datos de training, se pasan los valores de test y se devuelven las predicciones. Requiere 4 argumentos obligatorios:

• train: Una matriz (que no dataframe) que contenga el valor de los predictores para las observaciones que se van a emplear como trainig data.

• test: Una matriz (que no dataframe) que contenga el valor de los predictores para las nuevas observaciones cuya variable respuesta se quiere predecir, test data.

• cl: un vector que contenga la verdadera clasificación de las observaciones empleadas en el training data.

• K: el número de observaciones vecinas empleadas por el modelo. Cuantas más observaciones vecinas participen en una predicción, más se atenúa la influencia que tiene cada una de ellas por separado. Es por lo tanto el parámetro que determina la flexibilidad del modelo. Cuantos menor sea k, mayor la flexibilidad pero, a su vez, mayor el riesgo de overfitting.

La función knn() requiere fraccionar el data set en 3 partes. Es importante que se mantenga el orden para que los valores del argumento train se correspondan con los del argumento cl.

Dado que la escala de ambos predictores empleados en este ejercicio es la misma, no es necesario estandarizarlos.

library(class)
# En caso de empate entre los vecinos, se elige uno aleatoriamente.
# Esto influye en la reproducibilidad.
set.seed(1)
prediccion_knn <- knn(train = train_data[, c("Lag1", "Lag2")],
                      test = test_data[, c("Lag1", "Lag2")] ,
                      cl = train_data[,"Direction"], k = 3 )
head(prediccion_knn)

## [1] Down Down Down Up   Up   Down
## Levels: Down Up

Evaluación del modelo:

table(clase_predicha = prediccion_knn, clase_real = test_data$Direction)

##               clase_real
## clase_predicha Down Up
##           Down   48 55
##           Up     63 86

# % de aciertos
paste("% de acierto:", mean(prediccion_knn == test_data$Direction))

## [1] "% de acierto: 0.531746031746032"

# % de errores
paste("% de error:", mean(prediccion_knn != test_data$Direction))

## [1] "% de error: 0.468253968253968"

El modelo por K-NN empleando k=3 es capaz de predecir correctamente el 53% de las veces. Test error de 47%.

test_data$prediccion <- prediccion_knn
test_data$acierto <- ifelse(test = test_data$Direction == test_data$prediccion,
                            yes = "Si", no = "No")
ggplot(data = test_data, aes(x = Lag1, y = Lag2, color = Direction,
                             shape = acierto)) +
  geom_point() +
  scale_shape_manual(values = c("No" = 4, "Si" = 19 )) + 
  theme_bw()

Conclusión

De entre todos los métodos empleados con este set de datos, el QDA es el que menor test error ha mostrado y por lo tanto el que mejor recoge la relación entre los predictores Lag1 y Lag2, y la variable respuesta Direction.