4.6 Lab: Logistic Regression, LDA, QDA, and KNN

Paula Cazali

4.6.1 The Stock Market Data

library(ISLR)

names(Smarket)

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

summary(Smarket)

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

pairs(Smarket)

Correlaciones entre las variables.

cor(Smarket)

Error in cor(Smarket) : 'x' must be numeric

Da error porque la variable Direction es cualitativa

cor(Smarket[,-9])

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

Como se esperaba las correlaciones entre las variables Lag y Today son muy cercanas a cero. Parece que hay poca correlacion entre el retorno de ese dia y el returno del dia anterior. La unica correlacion sustancial es entre las variables Year y Volume. Con la grafica se puede ver que el Volumen esta incrementando con el tiempo.

NOTA: attach() carga a memoria el dataset, asi se podra referenciar a las columnas de este sin necesidad de poner dataset$columna. Usar solo en este caso.

attach(Smarket)

The following objects are masked from Weekly (pos = 14):

    Direction, Lag1, Lag2, Lag3, Lag4, Lag5, Today, Volume, Year

The following objects are masked from Weekly (pos = 15):

    Direction, Lag1, Lag2, Lag3, Lag4, Lag5, Today, Volume, Year

The following objects are masked from Smarket (pos = 16):

    Direction, Lag1, Lag2, Lag3, Lag4, Lag5, Today, Volume, Year

The following objects are masked from Smarket (pos = 17):

    Direction, Lag1, Lag2, Lag3, Lag4, Lag5, Today, Volume, Year

The following objects are masked from Smarket (pos = 18):

    Direction, Lag1, Lag2, Lag3, Lag4, Lag5, Today, Volume, Year

plot(Volume)

4.6.2 Logistic Regression

Predecir Direction usando todas las variables Lag y Volume. Usar la funcion glm(). La sintaxis es similar a la de la funcion lm(), con la diferencia que necesita el argumento family=binomial para que R corra una regresion logistica en vez de otro tipo de modelo lineal generalizado.

glm.fit <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume,
               data = Smarket,
               family = binomial)
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

Coeficientes del modelo:

coef(glm.fit)

 (Intercept)         Lag1         Lag2         Lag3         Lag4         Lag5       Volume 
-0.126000257 -0.073073746 -0.042301344  0.011085108  0.009358938  0.010313068  0.135440659 

Se puede accesar tambien a informacion especifica del modelo por medio de:

summary(glm.fit)$coef

                Estimate Std. Error    z value  Pr(>|z|)
(Intercept) -0.126000257 0.24073574 -0.5233966 0.6006983
Lag1        -0.073073746 0.05016739 -1.4565986 0.1452272
Lag2        -0.042301344 0.05008605 -0.8445733 0.3983491
Lag3         0.011085108 0.04993854  0.2219750 0.8243333
Lag4         0.009358938 0.04997413  0.1872757 0.8514445
Lag5         0.010313068 0.04951146  0.2082966 0.8349974
Volume       0.135440659 0.15835970  0.8552723 0.3924004

p-values (prueba para que el coeficiente sea significativo) para los coeficientes:

summary(glm.fit)$coef[,4]

(Intercept)        Lag1        Lag2        Lag3        Lag4        Lag5      Volume 
  0.6006983   0.1452272   0.3983491   0.8243333   0.8514445   0.8349974   0.3924004 

La funcion predict() puede ser usada para predecir que probabilidad tienen el mercado de subir, dados los valores de los predictores. Desplegando las primeras 10 probabilidades:

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

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

Se sabe que las probabilidades indican si sube o no ya que la funcion contrasts() indica que una probabilidad de 1 es para que el mercado suba.

contrasts(Direction)

     Up
Down  0
Up    1

Para hacer predicciones en base a si el mercado va a subir o a bajar en un dia en particular se necesita convertir estas probabilidades en labels de clase: Up y Down. Se crea un vector con “Down” y al tener una probabilidad mayor a 0.5 se asigna “Up”.

glm.pred <- rep("Down", 1250)
glm.pred[glm.probs > 0.5] <- "Up"

Ahora se usa la funcion table() la cual produce una matriz de confucion para determinar cuantas observaciones fueron correctamente clasificadas y cuales no.

table(glm.pred, Direction)

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

Los elementos en la diagonal de la matriz de confusion indica las predicciones correctas, mientras los otros valores indican las predicciones incorrectas.

Media de las predicciones acertadas:

(507+145)/1250

[1] 0.5216

mean(glm.pred == Direction)

[1] 0.5216

Por lo que la tasa de error del training es de:

100-52.2

[1] 47.8

Para mejorar la presicion del modelo de regresion logistica se puede usar solo una parte de los datos para el entrenamiento y luego la otra parte para hacer las pruebas.

Se va a crear un training set con las observaciones desde el anio 2001 hasta el 2004. El objeto train es un vector de 1250 elementos que corresponden a las observaciones del dataset. Loe elementos de ese vector que corresponden a las observaciones que ocurrieron antes del 2005 son TRUE y las que no son FALSE. Los vectores booleanos pueden ser usados para hacer subsets de las filas y las columnas de una matriz. Smartek[train,] sera una submatriz del data set de stock market.

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

[1] 252   9

Direction.2005 <- Direction[!train]

Ahora se creara un modelo de regresion logistica ussando solo el subset de las observaciones que corresponden a fechas antes del 2005, usando el argumento subset. Luego se predicen las probabilidades de que el mercado de acciones suba para cada uno de eso dias en el test set.

glm.fit <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume,
               data=Smarket,
               family=binomial,
               subset=train)
glm.probs <- predict(glm.fit, Smarket.2005, type = "response")

Ahora se hacen las predicciones y se despliega la matriz de confusion:

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

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

El promedio de predicciones y la tasa de error del modelo:

mean(glm.pred == Direction.2005)

[1] 0.4801587

mean(glm.pred!=Direction.2005)

[1] 0.5198413

El uso de predictores que no tiene relacion con la respuesta tiende a causar que la tasa de error empeore, ya que los predictores causan un aumento en la varianza sin una disminucion en el bias.

Ahora solo se usaran las variables Lag1 y Lag2 los cuales parecen ser los estadisticamente mas significativos.

glm.fit <- glm(Direction ~ Lag1 + Lag2,
               data = Smarket,
               family = binomial,
               subset = train)
glm.probs <- predict(glm.fit, Smarket.2005, type = "response")
glm.pred <- rep("Down", 252)
glm.pred[glm.probs > 0.5] <- "Up"
table(glm.pred,Direction.2005)

        Direction.2005
glm.pred Down  Up
    Down   35  35
    Up     76 106

mean(glm.pred == Direction.2005)

[1] 0.5595238

106/(106+76)

[1] 0.5824176

Ahora las predicciones tendran un porcentaje mayor: 56% Suponiendo que queremos predecir la direccion del mercado en un dia determinado cuando Lag1 y Lag2 son iguales a 1.2 y 1.1 respectivamente y en otro dia son iguales a 1.5 y -0.8.

predict(glm.fit, newdata = data.frame(Lag1 = c(1.2,1.5), 
                                      Lag2 = c(1.1, -0.8)),
        type = "response")

        1         2 
0.4791462 0.4960939 

4.6.3 Linear Discriminant Analysis

Ahora se usará LDA con el dataset de Smarket. Usaremos lda() de la libreria MASS. Esta funcion se parece a la funcion lm() y glm() excepto por el argumento de family. Usaremos las observaciones de antes del año 2005.

library(MASS)

lda.fit <- lda(Direction ~ Lag1 + Lag2,
               data = Smarket,
               subset = train)
lda.fit

Call:
lda(Direction ~ Lag1 + Lag2, data = Smarket, subset = train)

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

plot(lda.fit)

Los coeficientes del discriminante lineal son una combinacion lineal de Lag1 y Lag2. Si \(???0.642Lag_1???0.514Lag_2\) es mayor el clasificador LDA va a predecir una elevacion en el mercado y si es pequeño va a predecir una declinacion para el mercado. La funcion predict() devuelve una lista de tres elementos. El primer elemento class contiene los predictores de LDA acerca del movimiento del mercado. El segundo elemento posterior es una matriz donde la k-esima columna contiene la probabilidad posterior que corresponde a una observacion que pertenece a la clase k-esima. Finalmente x contiene los discriminates lineales.

lda.pred <- predict(lda.fit, Smarket.2005)
names(lda.pred)

[1] "class"     "posterior" "x"        

lda.class <- lda.pred$class 
table(lda.class,Direction.2005)

         Direction.2005
lda.class Down  Up
     Down   35  35
     Up     76 106

mean(lda.class == Direction.2005)

[1] 0.5595238

Aplicando el threshold de 0.5 para tener las predicciones:

sum(lda.pred$posterior[,1] >= 0.5)

[1] 70

sum(lda.pred$posterior[,1] < 0.5)

[1] 182

lda.pred$posterior[1:20,1]

      999      1000      1001      1002      1003      1004      1005      1006      1007      1008      1009      1010      1011 
0.4901792 0.4792185 0.4668185 0.4740011 0.4927877 0.4938562 0.4951016 0.4872861 0.4907013 0.4844026 0.4906963 0.5119988 0.4895152 
     1012      1013      1014      1015      1016      1017      1018 
0.4706761 0.4744593 0.4799583 0.4935775 0.5030894 0.4978806 0.4886331 

lda.class[1:20]

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

Suponer que queremos predecir que el mercado va a decrementar solamente si estamos certeron que el mercado va a decrecer en ese dia, si la probabilidad posterior es al menos del 90%

sum(lda.pred$posterior[,1] > 0.9)

[1] 0

Por lo que ningun dia del 2005 coincide con ese threshold.

4.6.4 Quadratic Discriminant Analysis

Usaremos la funcion qda() la cual tambien esta en la libreria MASS. La sintaxis de esta funcion es identica a la lda().

qda.fit <- qda(Direction ~ Lag1 + Lag2,
               data = Smarket,
               subset = train)
qda.fit

Call:
qda(Direction ~ Lag1 + Lag2, data = Smarket, subset = train)

Prior probabilities of groups:
    Down       Up 
0.491984 0.508016 

Group means:
            Lag1        Lag2
Down  0.04279022  0.03389409
Up   -0.03954635 -0.03132544

La funcion qda contiene las medias de los grupos. Pero no contiene los coeficientes del discriminante lineal, ya que el clasificador QDA involucra una funcion cuadratica de predictores. La funcion predict() funciona exactamente igual que LDA.

qda.class <- predict(qda.fit, Smarket.2005)$class
table(qda.class, Direction.2005)

         Direction.2005
qda.class Down  Up
     Down   30  20
     Up     81 121

mean(qda.class == Direction.2005)

[1] 0.5992063

Por lo que se puede ver que el modelo QDA es preciso el 60% de las veces incluso sin los datos del 2005. Por lo que se puede decir que el modelo QDA captura de forma mas precisa la relacion de las variables, de mejor manera que LDA y la regresion logistica.

4.6.5 K-Nearest Neighbors

Ahora se usara KNN usando la funcion knn() que es parte de la libreria class. Esta funcion es distinta a las otras ya que no necesita de dos pasos para hacer las predicciones. Necesita cuatro inputs.

1. Una matriz que contenga los predictores asociadss con los datos con el train data. En este ejercicio se llamara train.X.

2. Una matriz que contenga los predictores asociados con los datos con los cuales se quieren hacer predicciones. Esta matriz se llamara test.X.

3. Un vector que contiene las etiquetas de clase para las observaciones del training, train.Direction.

4. Un valor de k, el numero de vecindades cercanas para usar en el clasificador.

Se usara la funcion cbind() para juntar las variables Lag1 y Lag2 en dos matrices:

library(class)
train.X <- cbind(Lag1, Lag2)[train,]
test.X <- cbind(Lag1,Lag2)[!train,]
train.Direction <- Direction[train]

Ahora usamos la funcion knn() para predecir los movimientos del mercado en 2005. Hacemos un set de un random antes de aplicar la funcion, esto asegura que se pueda reproducir siempre el mismo resultado.

set.seed(1)
knn.pred <- knn(train.X, test.X, train.Direction, k = 1)
table(knn.pred, Direction.2005)

        Direction.2005
knn.pred Down Up
    Down   43 58
    Up     68 83

(43+83)/252

[1] 0.5

Las predicciones hechas con \(k=1\) no son muy buenas ya que solo el \(50\%\) de las predicciones son correctas. Ahora se reproduce el modelo usando \(k=3\)

knn.pred <- knn(train.X, test.X, train.Direction, k = 3)
table(knn.pred, Direction.2005)

        Direction.2005
knn.pred Down Up
    Down   48 54
    Up     63 87

mean(knn.pred == Direction.2005)

[1] 0.5357143

Se puede ver que hubo una pequeña mejora. Por lo que se puede concluir que para este analisis el mejor modelo es el QDA.

4.6.6 An Application to Caravan Insurance Data

Finalmente se va a aplicar el metodo KNN al dataset Caravan, que es parte de la libreria ISLR. Este set de datos incluye \(85\) predictores que miden las caracteristicas demograficas de \(5822\) individuos. La variable de respuesta sera Purchase, que indica si un individuo compra o no una poliza de seguro de caravana. Solo el \(6\%\) de las personas compraron un seguro de caravana.

dim(Caravan)

[1] 5822   86

attach(Caravan)

summary(Purchase)

  No  Yes 
5474  348 

348/5822

[1] 0.05977327

Para este problema se va a estandarizar, para que todas las variables tengan una media de \(0\) y una desviacion estandar de \(1\). Esto lo hace la funcion scale().

standardized.X <- scale(Caravan[,-86])
var(Caravan[,1])

[1] 165.0378

var(Caravan[,2])

[1] 0.1647078

var(standardized.X[,1])

[1] 1

var(standardized.X[,2])

[1] 1

Ahora cada columna del dataset standardized.X tiene una media de \(0\) y una desviacion estandar de \(1\). Ahora se va a crear un test data que contenga las primeras \(1000\) observaciones y un training set que contenga el resto.

test <- 1:1000
train.X <- standardized.X[-test,]
test.X <- standardized.X[test,]
train.Y <- Purchase[-test]
test.Y <- Purchase[test]
set.seed(1)
knn.pred <- knn(train.X, test.X, train.Y, k=1)
mean(test.Y != knn.pred)

[1] 0.118

mean(test.Y != "No")

[1] 0.059

Obtenemos la matriz de confusion:

table(knn.pred, test.Y)

        test.Y
knn.pred  No Yes
     No  873  50
     Yes  68   9

9/(68+9)

[1] 0.1168831

Usando ahora \(k=3\)

knn.pred <- knn(train.X, test.X, train.Y, k=3)
table(knn.pred, test.Y)

        test.Y
knn.pred  No Yes
     No  920  54
     Yes  21   5

5/26

[1] 0.1923077

knn.pred <- knn (train.X, test.X, train.Y, k = 5)
table(knn.pred, test.Y)

        test.Y
knn.pred  No Yes
     No  930  55
     Yes  11   4

4/15

[1] 0.2666667

Para comparar tambien se usara una regresion logistica. Usamos el valor de \(0.5\) para la probabilidad del clasificador

glm.fit <- glm(Purchase ~ ., 
               data = Caravan,
               family =binomial,
               subset = -test)

glm.fit: fitted probabilities numerically 0 or 1 occurred

glm.probs <- predict (glm.fit, Caravan[test,], type = "response")
glm.pred <- rep ("No ", 1000)
glm.pred[glm.probs > 0.5] <- "Yes"
table(glm.pred, test.Y)

        test.Y
glm.pred  No Yes
     No  934  59
     Yes   7   0

glm.pred <- rep("No", 1000)
glm.pred[glm.probs > 0.25] <- "Yes"
table(glm.pred, test.Y)

        test.Y
glm.pred  No Yes
     No  919  48
     Yes  22  11

11/(22+11)

[1] 0.3333333

Por lo que podemos predecir de manera correcta en \(33\%\) de estas personas. Lo cual es mucho mejor que adivinar random.