LAB_SEC_4.6
Milton Godinez
3/6/2019
ISLR laboratorio sección 4.6
library (ISLR)## Warning: package 'ISLR' was built under R version 3.4.4library (MASS)
library(corrplot)## Warning: package 'corrplot' was built under R version 3.4.4## corrplot 0.84 loadedlibrary(RColorBrewer)## Warning: package 'RColorBrewer' was built under R version 3.4.4Conociendo el set de datos
names(Smarket )## [1] "Year" "Lag1" "Lag2" "Lag3" "Lag4" "Lag5"
## [7] "Volume" "Today" "Direction"summary (Smarket )## Year Lag1 Lag2
## Min. :2001 Min. :-4.922000 Min. :-4.922000
## 1st Qu.:2002 1st Qu.:-0.639500 1st Qu.:-0.639500
## Median :2003 Median : 0.039000 Median : 0.039000
## Mean :2003 Mean : 0.003834 Mean : 0.003919
## 3rd Qu.:2004 3rd Qu.: 0.596750 3rd Qu.: 0.596750
## Max. :2005 Max. : 5.733000 Max. : 5.733000
## Lag3 Lag4 Lag5
## Min. :-4.922000 Min. :-4.922000 Min. :-4.92200
## 1st Qu.:-0.640000 1st Qu.:-0.640000 1st Qu.:-0.64000
## Median : 0.038500 Median : 0.038500 Median : 0.03850
## Mean : 0.001716 Mean : 0.001636 Mean : 0.00561
## 3rd Qu.: 0.596750 3rd Qu.: 0.596750 3rd Qu.: 0.59700
## Max. : 5.733000 Max. : 5.733000 Max. : 5.73300
## Volume Today Direction
## Min. :0.3561 Min. :-4.922000 Down:602
## 1st Qu.:1.2574 1st Qu.:-0.639500 Up :648
## Median :1.4229 Median : 0.038500
## Mean :1.4783 Mean : 0.003138
## 3rd Qu.:1.6417 3rd Qu.: 0.596750
## Max. :3.1525 Max. : 5.733000uso de la función cor() para obtener la matriz de coeficientes de correlacion entre todas las variables del dataset, tomar en cuenta que todas las variables deben ser numéricas para el cálculo
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.000000000forma grafica de matriz de correlacion
M <-cor(Smarket [,-9])
corrplot(M, type="upper", order="hclust",
col=brewer.pal(n=8, name="RdYlBu"))
se ve que la relación más fuerte es entre las variables Year y Volume.
Al graficar los datos se ve un crecimiento sobre volumen en los últimos años.
attach (Smarket )
plot(Volume )
Regresión logística
Se realizará un modelo de regresión logística para predecir Direction utilizando las variables de Lag1 a Lag5 y Volume. La sintaxis para la función glm() es similar a lm() con la diferencia que se debe pasar el parámetro family=binomial para indicar a R que es un modelo de regresión y no cualquier modelo lineal.
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: 3la función coef() nos permite acceder solo a los coeficientes de este modelo.
coef(glm.fits)## (Intercept) Lag1 Lag2 Lag3 Lag4
## -0.126000257 -0.073073746 -0.042301344 0.011085108 0.009358938
## Lag5 Volume
## 0.010313068 0.135440659También podemos utilizar la función summary() para acceder a cierta información relevante como los valors p.
summary (glm.fits)$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.3924004La funcion predict() se puede utilizar para predicir la probabilidad de que el mercado suba según los valores de los predictores. El parámetro type=“response” hace que devuelva las probabilidades en la forma P(Y = 1|X).
glm.probs =predict(glm.fits, type ="response")
glm.probs [1:10]## 1 2 3 4 5 6 7
## 0.5070841 0.4814679 0.4811388 0.5152224 0.5107812 0.5069565 0.4926509
## 8 9 10
## 0.5092292 0.5176135 0.4888378Los valores que devuelve la función corresponden a la probabildad de subir del mercado, porque la función contrasts() nos indica que R ha creado una variable genérica con el valor 1 para Up.
contrasts(Direction)## Up
## Down 0
## Up 1para poder predecir necesitamos convertir esas probabilidades en etiquetas que nos indiquen si sube o baja y eso se logra con las siguientes instrucciones.
La primer linea crea un vector de 1,250 valores idnicando “Down”, la segunda linea convierte a “Up” todos los registros cuya probabilidad predicha es mayor a 0.5.
glm.pred=rep ("Down" ,1250)
glm.pred[glm.probs > .5]="Up"Luego de eso podemos utilizar la función table() para generar nuestra matriz de confusión y así saber los valores que acertamos y los fallidos.
table(glm.pred, Direction )## Direction
## glm.pred Down Up
## Down 145 141
## Up 457 507la diagonal son las predicciones correctas y lo que se encuentra fuera de la diagonal son errores en la prediccion. La función mean() nos ayudar a calcular el porcentaje de aciertos que tuvimos.
mean(glm.pred== Direction )## [1] 0.5216en este caso usamos todo el dataset para entrenar el modelo, ahora usaremos los datos de 2001 a 2004 para entrenar el modelo y luego usaremos el 2005 para comprobar la efectividad.
train =(Year <2005)
Smarket.2005 = Smarket [!train ,]
dim(Smarket.2005)## [1] 252 9Direction.2005= Direction [!train]ahora ajustamos un modelo con los datos anteriores a 2005 utilizando el argumento “subset”
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")por último comparamos los resultados obtenidos con los datos reales
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 44mean(glm.pred== Direction.2005)## [1] 0.4801587mean(glm.pred!= Direction.2005)## [1] 0.5198413con esto obtenes resultados peores a una asignación aleatoria.
Si recordamos la variable “Lag1” tenía el valor p mas bajo por lo que podemos probar a retirar las variables que no aportan al modelo y volver a ajustarlo.
Probaremos de nuevo utilizando solo las variables “Lag1” y “Lag2” que son las que parecen tener mas potencial para aportar al modelo.
glm.fits = glm(
Direction∼Lag1 + Lag2 ,
data = Smarket ,
family = binomial ,
subset = train
)
glm.probs = predict (glm.fits, Smarket.2005 , type = "response")
glm.pred = rep ("Down" , 252)
glm.pred[glm.probs > .5] = "Up"
table(glm.pred, Direction.2005)## Direction.2005
## glm.pred Down Up
## Down 35 35
## Up 76 106mean(glm.pred== Direction.2005)## [1] 0.5595238106/(106+76)## [1] 0.5824176los datos han mejorado muy poca dando una efetividad general de 56% pero si nos fijamos solo en la predicción de “Up” la efectividad es de 58% por lo que podriamos adoptar la estrategia de solo utilizar el modelo cuando predice un aumento e ignorarlo cuando predice la baja.
A continuación veremos como utilizar el modelo si quisieramos predecir solo para ciertos días con los valores conocidos.
predict (glm.fits,
newdata = data.frame(Lag1 = c(1.2 , 1.5) , Lag2 = c(1.1 ,-0.8)),
type = "response")## 1 2
## 0.4791462 0.4960939Linear Discriminant Analysis (LDA)
Para ajustar este modelo utilizamos la función lda() que tiene un estructura muy similar a lm() y glm() paro sin el parámetro “family”.
Ajustaremos con los datos anteriores a 2005.
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.5135293plot(lda.fit)
los resultados indican ˆπ1 = 0.492 y ˆπ2 = 0.508 es decir que el 49.2% de de las observaciones de entrenamiento corresponden a días en los que el mercado bajó. En los resultados también se nos muestra las medias de las variables usadas y vemos que los 2 días previos tienen una media negativa cuando el mercado sube y son positivos cuando el mercado tiende a la baja.
La función predict() devuelve los valores: “class” nos indica la predicción del modelo sobre el movimiento del mercado. El segundo, “posterior”, es una matriz cuya n-ésima columna indica la probabilidad de que x observación pertenezca a n-ésima clase.
lda.pred=predict (lda.fit , Smarket.2005)
names(lda.pred)## [1] "class" "posterior" "x"podemos ver que la predicción entre la regresión lineal y LDA son casi idénticas
lda.class = lda.pred$class
table(lda.class , Direction.2005)## Direction.2005
## lda.class Down Up
## Down 35 35
## Up 76 106si aplicamos un corte a los probabilidades en 0.5 y asignamos las etiquetas nosotros mismos podremos replicar los resultados obtenidos anteriormente.
sum(lda.pred$posterior[ ,1] >=.5)## [1] 70sum(lda.pred$posterior[,1]<.5)## [1] 182lda.pred$posterior[1:20 ,1]## 999 1000 1001 1002 1003 1004 1005
## 0.4901792 0.4792185 0.4668185 0.4740011 0.4927877 0.4938562 0.4951016
## 1006 1007 1008 1009 1010 1011 1012
## 0.4872861 0.4907013 0.4844026 0.4906963 0.5119988 0.4895152 0.4706761
## 1013 1014 1015 1016 1017 1018
## 0.4744593 0.4799583 0.4935775 0.5030894 0.4978806 0.4886331lda.class [1:20]## [1] Up Up Up Up Up Up Up Up Up Up Up Down Up Up
## [15] Up Up Up Down Up Up
## Levels: Down Upsi queremos utilizar un corte diferente a 0.5, por ejemplo 0.9
sum(lda.pred$posterior[,1]>.9)## [1] 0y vemos que ningún día de 2005 alcanza ese corte
Quadratic Discriminant Analysis (QDA)
Para ajustar un modelo de este tipo utilizamos la función qda() que tiene un uso similar a las anteriores.
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.03132544El resultado nos muestra una media de grupos pero no contiene los coeficiones de los descriminantes lineares porque QDA utiliza un cuadratico en lugar de un linear. La función predict() funciona de la misma forma.
qda.class = predict (qda.fit , Smarket.2005)$class
table(qda.class , Direction.2005)## Direction.2005
## qda.class Down Up
## Down 30 20
## Up 81 121mean(qda.class == Direction.2005)## [1] 0.5992063Este modelo consigue un 60% de efectividad superando a los anteriores, se aconseja probarlo en en un dataset mas grande para estar más seguros.
K-Nearest Neighbors (KNN)
Este modelo se configura de forma diferente pero mantiene la estructura de primero ajustar un modelo y luego utilizar para buscar una predicción. La función que utilizaremos es knn() y requiere los siguientes parametros.
- Una matriz con los predictores de la matriz para entrenamiento etiquetada “train.X”
- Una matriz con los predictores de la matriz que deseamos utilizar para predecir, etiquetada como “test.X”
- Un vector que contenga las etiquetas de clases por las observaciones de entrenamiento, llamada train.Direction.
- Un valor “k” que es el numero de vecinos que se utilizaran para la clasificación.
Formaremos las matrices mencionadas con la función cbind()
library (class)
train.X=cbind(Lag1 ,Lag2)[train ,]
test.X=cbind (Lag1 ,Lag2)[!train ,]
train.Direction =Direction [train]Con esto podemos realizar el ajuste del modelo, para este ejemplo estableceremos la semilla que utilizará R para que los resultados puedan ser replicables.
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(83+43) /252## [1] 0.5los resultados utilizando k=1 no son efectivos, ya que solo el 50% de las observaciones se predicen correctamente. Utilizar un K mayor podría ser de mejora para el modelo, probaremos con 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 87mean(knn.pred == Direction.2005)## [1] 0.5357143Los resultados tuvieron muy poca mejora, por lo que parece que QDA tiene el mejor desempeño entre los métodos examinados.