Fiabilidad, Laboratorio 1
Laboratorio 1: Arboles de Decision
Ajustando Los Arboles de Clasificacion
Para ello empesaremos utilizando la libreria tree de R
Para este ejercicio usaremos el data set Carseats, en donde la propiedad Sales es continua y la transformaremos en binaria, creando la variable High utilizando la funcion ifelse de R, donde esta sera Yes si las ventas son mayores a 8 y No para lo contrario:
library(ISLR)
car_seats <- Carseats
car_seats$High <- ifelse(Carseats$Sales <= 8, 0,1)
car_seats$ShelveLoc <- as.factor(car_seats$ShelveLoc)
car_seats$Urban <- as.factor(car_seats$Urban)
car_seats$US <- as.factor(car_seats$US)
car_seatsAhora usando la funcion tree para poder ajustar un arbol para clasificar la variable High utilizando todas las variables del dataset excluyendo Sales:
tree.carseats =tree (High~.-Sales, car_seats )
summary(tree.carseats)
Regression tree:
tree(formula = High ~ . - Sales, data = car_seats)
Variables actually used in tree construction:
[1] "ShelveLoc" "Price" "Income" "Advertising" "CompPrice" "Age"
Number of terminal nodes: 19
Residual mean deviance: 0.08725 = 33.24 / 381
Distribution of residuals:
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.93750 -0.06250 -0.02000 0.00000 0.06667 0.98000 Para los arboles de clasificacion las deviacion reportada esta dada por la siguiente formula:
\[\begin{aligned} -2\Sigma_m\Sigma_k n_{mk}log\hat{p}_{mk} \end{aligned}\]Donde \(n_{mk}\) es el numero de observaciones en el m avo nodo terminal que terneces a la k clase. Una pequena desviacion indica al arbol que provee un buen ajuste para la data. La desviacion residual media reportada es simplemente la desviacion dividido \(n-|T_0|\), el cual es el caso para \(400 - 19 = 28\). Una de las ventajas de la funcion tree es que nos permite graficarla utilizando plot:
plot(tree.carseats)
text(tree.carseats, pretty = 0)
Esto es un indicador de que la clasificacion mas importante es Shelveloc, ahora veamos si simplemente imprimimos la variable tree.carseats nos da las clasificaciones de los arboles:
tree.carseats
node), split, n, deviance, yval
* denotes terminal node
1) root 400 96.7600 0.41000
2) ShelveLoc: Bad,Medium 315 67.5100 0.31110
4) Price < 92.5 46 9.7390 0.69570
8) Income < 57 10 2.1000 0.30000 *
9) Income > 57 36 5.6390 0.80560 *
5) Price > 92.5 269 49.8100 0.24540
10) Advertising < 13.5 224 33.5000 0.18300
20) CompPrice < 124.5 96 5.6250 0.06250 *
21) CompPrice > 124.5 128 25.4300 0.27340
42) Price < 109.5 21 4.2860 0.71430
84) ShelveLoc: Bad 5 0.0000 0.00000 *
85) ShelveLoc: Medium 16 0.9375 0.93750 *
43) Price > 109.5 107 16.2600 0.18690
86) Price < 126.5 42 9.3330 0.33330
172) Age < 49.5 20 4.8000 0.60000
344) CompPrice < 137 14 3.4290 0.42860 *
345) CompPrice > 137 6 0.0000 1.00000 *
173) Age > 49.5 22 1.8180 0.09091 *
87) Price > 126.5 65 5.4460 0.09231
174) CompPrice < 147.5 50 0.9800 0.02000 *
175) CompPrice > 147.5 15 3.3330 0.33330
350) CompPrice < 151.5 8 1.8750 0.62500 *
351) CompPrice > 151.5 7 0.0000 0.00000 *
11) Advertising > 13.5 45 11.1100 0.55560
22) Age < 54.5 25 4.0000 0.80000 *
23) Age > 54.5 20 3.7500 0.25000
46) CompPrice < 122.5 10 0.0000 0.00000 *
47) CompPrice > 122.5 10 2.5000 0.50000
94) Price < 125 5 0.0000 1.00000 *
95) Price > 125 5 0.0000 0.00000 *
3) ShelveLoc: Good 85 14.7500 0.77650
6) Price < 142.5 73 8.6300 0.86300
12) Income < 34.5 13 3.2310 0.53850
24) Price < 109.5 6 0.0000 1.00000 *
25) Price > 109.5 7 0.8571 0.14290 *
13) Income > 34.5 60 3.7330 0.93330 *
7) Price > 142.5 12 2.2500 0.25000 *Para poder evaluar que el arbol es bueno para predecir medirlo con trainging data y test data:
set.seed(2)
train=sample(1:nrow(car_seats), 200)
car_seats.test <- car_seats[-train,]
high.test <- car_seats$High[-train]
tree.carseats <- tree(as.factor(High)~.-Sales, car_seats, subset=train)
tree.pred <- predict(tree.carseats, car_seats.test, type="class")
table(tree.pred, high.test) high.test
tree.pred 0 1
0 86 27
1 30 57Ahora veremos si podando el arbol podemos encontrar mejores resultados, para ello ustaremo las fucion cv.tree, la cual hace una validacion para verificar el nivel optimo de complejidad, el costo de esta operacion va de la mano de la seleccion de sequencias de arboles que usemos. Para ello especificamos el parametro FUN=prune.misclass de manera que le indicamos que utilize la tasa de errores en la clasificacion como criterio para podar el arbol:
set.seed(3)
cv.carseats <- cv.tree(tree.carseats, FUN=prune.misclass)
names(cv.carseats)[1] "size" "dev" "k" "method"cv.carseats
$size
[1] 19 17 14 13 9 7 3 2 1
$dev
[1] 55 55 53 52 50 56 69 65 80
$k
[1] -Inf 0.0000000 0.6666667 1.0000000 1.7500000 2.0000000 4.2500000 5.0000000 23.0000000
$method
[1] "misclass"
attr(,"class")
[1] "prune" "tree.sequence"Con esto vemos que el arbol con menos errores de clasificacion es el numero 9 con unicamente 50 errores de clasificacion. Ahora grafiquemos esto:
par(mfrow=c(1,2))
plot(cv.carseats$size, cv.carseats$dev, type="b")
plot(cv.carseats$k, cv.carseats$dev, type="b")
Ahora podemos aplicar prune.misclass() para podar el arbol hasta dejarlo unicamente con 9 nodos.
prune.carseats <- prune.misclass(tree.carseats, best=9)
plot(prune.carseats)
text(prune.carseats, pretty=0)
Veamos ahora que tambien predice:
tree.pred <- predict(prune.carseats, car_seats.test, type="class")
tbl <- table(tree.pred, high.test)
tbl high.test
tree.pred 0 1
0 94 24
1 22 60(tbl[1,1] + tbl[2,2])/sum(tbl)[1] 0.77Con esto vemos que el 77% de las observaciones del test fueron acertadas, por lo que no solo simplificamos el arbol si no que mejoramos la presion del arbol. Finalmete hagamos la prueba utilizando un arbol con mayor complejidad:
prune.carseats <- prune.misclass(tree.carseats, best=15)
plot(prune.carseats)
text(prune.carseats, pretty=0)
tree.pred <- predict(prune.carseats, car_seats.test, type="class")
tbl <- table(tree.pred, high.test)
tbl high.test
tree.pred 0 1
0 86 22
1 30 62(tbl[1,1] + tbl[2,2])/sum(tbl)[1] 0.74Ajustando Arboles de Regresion
Ahora ajustaremos un arbol de regresion al dataset Boston:
library(MASS)
set.seed(1)
train <- sample(1:nrow(Boston), nrow(Boston)/2)
tree.boston <- tree(medv~., Boston, subset=train)
summary(tree.boston)
Regression tree:
tree(formula = medv ~ ., data = Boston, subset = train)
Variables actually used in tree construction:
[1] "lstat" "rm" "dis"
Number of terminal nodes: 8
Residual mean deviance: 12.65 = 3099 / 245
Distribution of residuals:
Min. 1st Qu. Median Mean 3rd Qu. Max.
-14.10000 -2.04200 -0.05357 0.00000 1.96000 12.60000 Ahora grafiquemos nuestro arbol
plot(tree.boston)
text(tree.boston, pretty=0)
Ahora veamos si podemos podar un poco este arbol:
cv.boston <- cv.tree(tree.boston)
plot(cv.boston$size, cv.boston$dev, type="b")
Como vemos de la grafica aca el arbol mas complejo es el que mejor prediccion tiene, aunque pareciera no haber mucha diferencia con el de 7 nodos. De cualquier manera haremos la prueba de podarlo:
prune.boston <- prune.tree(tree.boston, best=5)
plot(prune.boston)
text(prune.boston, pretty=0)
Finalmente utilizando el arbol sin podar haremos pruebas de prediccion para mediar la presicion del arbol:
yhat <- predict(tree.boston, newdata = Boston[-train,])
boston.test <- Boston[-train, "medv"]
plot(yhat, boston.test)
abline(0,1)
mean((yhat-boston.test)^2)[1] 25.04559Esto nos dice que el MSE es de 25, su raiz cuadrada es de 5.005, por lo que podemos determinar que el precio promedio de las casas en Boston es $5005.
Bagging y Random Forest
Ahora utilizando el paquete randomForest de R generaremos un arbol para el dataset de Boston, donde bagging es un caso especial de random forest donde \(m=p\) y utilizamos el mismo paquete para ambos casos:
library(randomForest)
set.seed(1)
bag.boston <- randomForest(medv~., data=Boston, subset = train, mtry=13, importance=TRUE)
bag.boston
Call:
randomForest(formula = medv ~ ., data = Boston, mtry = 13, importance = TRUE, subset = train)
Type of random forest: regression
Number of trees: 500
No. of variables tried at each split: 13
Mean of squared residuals: 11.15723
% Var explained: 86.49Al utilizar el parametro mtry=13 le estamos indicando que todos los parametros deben ser utilizados por lo que el proceso de bagging se asegura, ahora veamos como este arbol predice:
yhat.bag <- predict(bag.boston, newdata = Boston[-train,])
plot(yhat.bag, boston.test)
abline(0,1)
mean((yhat.bag-boston.test)^2)[1] 13.50808Como vemos el MSE nos quedo a la mitad de lo que obtuvimos con un arbol podado, ahora utilizando el paremtro ntree de randomForest podemos utilizar varios arboles:
bag.boston <- randomForest(medv~., data=Boston, subset=train, mtry=13, ntree=25)
yhat.bag <- predict(bag.boston, newdata = Boston[-train,])
mean((yhat.bag-boston.test)^2)[1] 13.94835Ahora veamos lo que sucede cuando reducimos el valor de mtry
set.seed(1)
rf.boston <- randomForest(medv~., data=Boston, subset=train, mtry=6, importance = TRUE)
yhat.rf <- predict(rf.boston, newdata = Boston[-train,])
mean((yhat.rf-boston.test)^2)[1] 11.66454Como vemos aca random Forest nos devolvio un MSE aun menor, por lo que esta tecnica funciono mejor que el Bagging, ahora veamos la importancia del arbol:
importance(rf.boston) %IncMSE IncNodePurity
crim 12.132320 986.50338
zn 1.955579 57.96945
indus 9.069302 882.78261
chas 2.210835 45.22941
nox 11.104823 1044.33776
rm 31.784033 6359.31971
age 10.962684 516.82969
dis 15.015236 1224.11605
rad 4.118011 95.94586
tax 8.587932 502.96719
ptratio 12.503896 830.77523
black 6.702609 341.30361
lstat 30.695224 7505.73936Este output nos da dos varaibles, primero el porcentaje de incremento del MSE. El segundo es una medida de impuresa de la variable, la cual aumenta las particiones de los valores de la variable. Veamos ahora la importancia de las variables:
varImpPlot(rf.boston)
Esta grafica nos termina de confirmar que las variables rm y lstat son las de mayor importancia en el modelo.
Bosting
Aqui usaremos el paquete gbm de R, lo cual nos permite ajustar los arboles de regresion para el dataset de Boston. Al tratarse de un problema de regresion utilizaremos el parametro distribution="gaussian", si se tratara de un problema de clasificacion utilizariamos distribution="bernoulli", el parametro n.trees=5000 indica que queremos utilizar 5000 arboles y interaction.depth=4 limitamos a solo 4 nodos de profundidad en cada uno:
library(gbm)
set.seed(1)
boost.boston <- gbm(medv~., data=Boston[-train,], distribution="gaussian", n.trees = 5000, interaction.depth = 4)
summary(boost.boston)
Nuevamente confirmamos por mucho querm y lstat son las variables con mayor importancia en el modelo, ahora intentemos graficar por partes:
par(mfrow=c(1,2))
plot(boost.boston, i="rm")
plot(boost.boston, i="lstat")
Ahora utilizaremos el modelo boosteado para predecir:
yhat.boston <- predict(boost.boston, newdata = Boston[-train,], n.trees = 5000)
mean((yhat.boston-boston.test)^2)[1] 9.03777e-07Veemos que en este caso el MSE es extremadamente bajo, ahora veamos si ajustamos los parametros \(\lambda\) del modelo
boost.boston <- gbm(medv~., data=Boston[-train,], distribution="gaussian", n.trees = 5000, interaction.depth = 4, shrinkage = 0.2, verbose = F)
yhat.boston <- predict(boost.boston, newdata = Boston[-train,], n.trees = 5000)
mean((yhat.boston-boston.test)^2)[1] 5.186409e-12Veemos esto aun mejora el modelo grandemente.