Laboratorio 8.3 - Arboles de decision
Cesar Tinoco - 13003387
8.3.1 Clasificando Arboles de Clasificacion
La libreria Tree se utiliza para construir arboles de clasificacion y regresion. Ahora cargaremos la libreria.
library(tree)Primero usamos classification tree para analizar el conjunto de datos de Carseats. Sale es una variable continua, por lo que comenzamos por recodificarla como una variable binaria. Usamos la función ifelse() para crear una variable, llamada Alto, que toma un valor de Yes si la variable Sale supera 8, y toma un valor de No en caso contrario.
library (ISLR)
attach (Carseats )
High=ifelse (Sales <=8," No"," Yes ")Finalmente, usamos la funcion data.frame() para fusionar High con el resto de los datos de Carseats.
Carseats =data.frame(Carseats ,High)Ahora usamos la funcion tree() para ajustar un arbol de clasificacion con el fin de predecir High usando todas las variables excepto Sale. La sintaxis de la función tree() es bastante similar a la de la funcion lm().
tree.carseats =tree(High∼.-Sales ,Carseats )La funcion de summary() enumera las variables que se utilizan como nodos internos en el arbol, el numero de nodos terminales y la tasa de error (entrenamiento).
summary (tree.carseats )
Classification tree:
tree(formula = High ~ . - Sales, data = Carseats)
Variables actually used in tree construction:
[1] "ShelveLoc" "Price" "Income" "CompPrice" "Population" "Advertising" "Age"
[8] "US"
Number of terminal nodes: 27
Residual mean deviance: 0.4575 = 170.7 / 373
Misclassification error rate: 0.09 = 36 / 400 Vemos que la tasa de error de entrenamiento es del 9%. Una de las propiedades mas atractivas de los arboles es que se ppueden mostrar graficamente. Usamos la funcion plot() para mostrar la estructura del arbol, y la funcion text() para mostrar las etiquetas de nodo. El argumento pretty = 0 indica a R que ha incluido los nombres de categoria para cualquier predictor cualitativo, en lugar de simplemente mostrar una letra para cada categoria.
plot(tree.carseats )
text(tree.carseats ,pretty =0)
Las ramas que conducen a los nodos terminales se indican mediante asteriscos.
tree.carseatsnode), split, n, deviance, yval, (yprob)
* denotes terminal node
1) root 400 541.500 No ( 0.59000 0.41000 )
2) ShelveLoc: Bad,Medium 315 390.600 No ( 0.68889 0.31111 )
4) Price < 92.5 46 56.530 Yes ( 0.30435 0.69565 )
8) Income < 57 10 12.220 No ( 0.70000 0.30000 )
16) CompPrice < 110.5 5 0.000 No ( 1.00000 0.00000 ) *
17) CompPrice > 110.5 5 6.730 Yes ( 0.40000 0.60000 ) *
9) Income > 57 36 35.470 Yes ( 0.19444 0.80556 )
18) Population < 207.5 16 21.170 Yes ( 0.37500 0.62500 ) *
19) Population > 207.5 20 7.941 Yes ( 0.05000 0.95000 ) *
5) Price > 92.5 269 299.800 No ( 0.75465 0.24535 )
10) Advertising < 13.5 224 213.200 No ( 0.81696 0.18304 )
20) CompPrice < 124.5 96 44.890 No ( 0.93750 0.06250 )
40) Price < 106.5 38 33.150 No ( 0.84211 0.15789 )
80) Population < 177 12 16.300 No ( 0.58333 0.41667 )
160) Income < 60.5 6 0.000 No ( 1.00000 0.00000 ) *
161) Income > 60.5 6 5.407 Yes ( 0.16667 0.83333 ) *
81) Population > 177 26 8.477 No ( 0.96154 0.03846 ) *
41) Price > 106.5 58 0.000 No ( 1.00000 0.00000 ) *
21) CompPrice > 124.5 128 150.200 No ( 0.72656 0.27344 )
42) Price < 122.5 51 70.680 Yes ( 0.49020 0.50980 )
84) ShelveLoc: Bad 11 6.702 No ( 0.90909 0.09091 ) *
85) ShelveLoc: Medium 40 52.930 Yes ( 0.37500 0.62500 )
170) Price < 109.5 16 7.481 Yes ( 0.06250 0.93750 ) *
171) Price > 109.5 24 32.600 No ( 0.58333 0.41667 )
342) Age < 49.5 13 16.050 Yes ( 0.30769 0.69231 ) *
343) Age > 49.5 11 6.702 No ( 0.90909 0.09091 ) *
43) Price > 122.5 77 55.540 No ( 0.88312 0.11688 )
86) CompPrice < 147.5 58 17.400 No ( 0.96552 0.03448 ) *
87) CompPrice > 147.5 19 25.010 No ( 0.63158 0.36842 )
174) Price < 147 12 16.300 Yes ( 0.41667 0.58333 )
348) CompPrice < 152.5 7 5.742 Yes ( 0.14286 0.85714 ) *
349) CompPrice > 152.5 5 5.004 No ( 0.80000 0.20000 ) *
175) Price > 147 7 0.000 No ( 1.00000 0.00000 ) *
11) Advertising > 13.5 45 61.830 Yes ( 0.44444 0.55556 )
22) Age < 54.5 25 25.020 Yes ( 0.20000 0.80000 )
44) CompPrice < 130.5 14 18.250 Yes ( 0.35714 0.64286 )
88) Income < 100 9 12.370 No ( 0.55556 0.44444 ) *
89) Income > 100 5 0.000 Yes ( 0.00000 1.00000 ) *
45) CompPrice > 130.5 11 0.000 Yes ( 0.00000 1.00000 ) *
23) Age > 54.5 20 22.490 No ( 0.75000 0.25000 )
46) CompPrice < 122.5 10 0.000 No ( 1.00000 0.00000 ) *
47) CompPrice > 122.5 10 13.860 No ( 0.50000 0.50000 )
94) Price < 125 5 0.000 Yes ( 0.00000 1.00000 ) *
95) Price > 125 5 0.000 No ( 1.00000 0.00000 ) *
3) ShelveLoc: Good 85 90.330 Yes ( 0.22353 0.77647 )
6) Price < 135 68 49.260 Yes ( 0.11765 0.88235 )
12) US: No 17 22.070 Yes ( 0.35294 0.64706 )
24) Price < 109 8 0.000 Yes ( 0.00000 1.00000 ) *
25) Price > 109 9 11.460 No ( 0.66667 0.33333 ) *
13) US: Yes 51 16.880 Yes ( 0.03922 0.96078 ) *
7) Price > 135 17 22.070 No ( 0.64706 0.35294 )
14) Income < 46 6 0.000 No ( 1.00000 0.00000 ) *
15) Income > 46 11 15.160 Yes ( 0.45455 0.54545 ) *Este enfoque conduce a predicciones correctas para alrededor del 71.5% de las ubicaciones en el conjunto de datos de prueba.
set.seed (2)
train=sample (1: nrow(Carseats ), 200)
Carseats.test=Carseats [-train ,]
High.test=High[-train ]
tree.carseats =tree(High∼.-Sales ,Carseats ,subset =train )
tree.pred=predict (tree.carseats ,Carseats.test ,type ="class")
table(tree.pred ,High.test) High.test
tree.pred No Yes
No 86 27
Yes 30 57(86+57) /200[1] 0.715A continuacion, consideramos si la poda del arbol puede llevar a mejores resultados. La función cv.tree() realiza una validacion cruzada para determinar el nivel optimo de complejidad del arbol; la poda de complejidad de costos se utiliza para seleccionar una secuencia de arboles para su consideracion. Utilizamos el argumento FUN=prune.misclass para indicar que queremos que la tasa de error de clasificacion guie el proceso de validacion cruzada y poda, en lugar del valor predeterminado para la función cv.tree(), que es la desviacion. La funcion cv.tree() informa el numero de nodos terminales de cada arbol considerado(tamaño), asi como la tasa de error correspondiente y el valor del parametro de costo-complejidad utilizado (k, que corresponde a α en (8.4)).
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"Tenga en cuenta que, a pesar del nombre, dev corresponde a la tasa de error de validacion cruzada en esta instancia. El arbol con 9 nodos terminales produce la tasa de error de validacion cruzada mas baja, con 50 errores de validacion cruzada. Trazamos la tasa de error como una funcion de tamaño y k.
par(mfrow =c(1,2))
plot(cv.carseats$size ,cv.carseats$dev ,type="b")
plot(cv.carseats$k ,cv.carseats$dev ,type="b")
Ahora aplicamos la funcion prune.misclass() para podar el arbol y obtener el arbol de nueve nodos.
prune.carseats =prune.misclass (tree.carseats ,best =9)
plot(prune.carseats )
text(prune.carseats ,pretty =0)
¿Que tan bien se desempeña este arbol podado en el conjunto de datos de prueba? Una vez mas, aplicamos la funcion predict().
tree.pred=predict (prune.carseats , Carseats.test ,type="class")
table(tree.pred ,High.test) High.test
tree.pred No Yes
No 94 24
Yes 22 60(94+60) /200[1] 0.77Ahora, el 77% de las observaciones de prueba se clasifican correctamente, por lo que el proceso de poda no solo produjo un arbol mas interpretable, sino que tambien mejoro la precision de la clasificacion. Si aumentamos el valor de best, obtenemos un arbol podado mas grande con una menor precision de clasificacion.
prune.carseats =prune.misclass (tree.carseats ,best =15)
plot(prune.carseats )
text(prune.carseats ,pretty =0)
tree.pred=predict (prune.carseats , Carseats.test ,type="class")
table(tree.pred ,High.test) High.test
tree.pred No Yes
No 86 22
Yes 30 62(86+62) /200[1] 0.748.3.2 Ajuste de Arboles de Regresion
Aqui ajustamos un arbol de regresion al conjunto de datos de Boston. Primero, creamos un conjunto de entrenamiento, y ajustamos el arbol a los datos de entrenamiento.
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 Observe que la salida de summary() indica que solo tres de las variables se han utilizado para construir el arbol. En el contexto de un arbol de regresion, la desviacion es simplemente la suma de los errores cuadrados para el arbol. Ahora trazamos el árbol.
plot(tree.boston )
text(tree.boston ,pretty =0)
Ahora usamos la funcion cv.tree() para ver si la poda del arbol mejorar el rendimiento.
cv.boston =cv.tree(tree.boston )
plot(cv.boston$size ,cv.boston$dev ,type='b')
En este caso, el arbol mas complejo se selecciona mediante validacion cruzada. Sin embargo, si deseamos podar el arbol, podriamos hacerlo de la siguiente manera, utilizando la funcion prune.tree()
prune.boston =prune.tree(tree.boston ,best =5)
plot(prune.boston )
text(prune.boston ,pretty =0)
De acuerdo con los resultados de la validacion cruzada, usamos el arbol no afinado para hacer predicciones en el conjunto de pruebas
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.04559La raiz cuadrada del MSE es por lo tanto alrededor de 5.005, indicando que este modelo conduce a predicciones de prueba que estan dentro de aproximadamente $ 5, 005 del valor medio real de la vivienda para el suburbio.
8.3.3 Sacos y Bosques Aleatorios
Aqui aplicamos empaquetamiento y bosques aleatorios a los datos de Boston, usando el paquete randomForest en R. Los resultados exactos obtenidos en esta seccion pueden depender de la version de R y la version del paquete randomForest instalada en la computadora, Recordemos que el embolsado es simplemente un caso especial de un bosque al azar con m = p. Por lo tanto, la funcion randomForest() se puede utilizar para realizar tanto bosques aleatorios como empaquetamiento. Realizamos embolsado como sigue.
library (randomForest)randomForest 4.6-14
Type rfNews() to see new features/changes/bug fixes.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.49El argumento mtry = 13 indica que se deben considerar los 13 predictores para cada division del arbol; en otras palabras, se debe realizar el empaquetamiento. ¿Que tan bien se desempeña este modelo embolsado en el conjunto de prueba?
yhat.bag = predict (bag.boston ,newdata =Boston [-train ,])
plot(yhat.bag , boston.test)
abline (0,1)
mean(( yhat.bag -boston.test)^2)[1] 13.50808El MSE del conjunto de pruebas asociado con el arbol de regresion empaquetado es 13.16, casi la mitad que se obtiene utilizando un arbol unico optimizado. Podriamos cambiar el numero de arboles cultivados por randomForest() usando el argumento ntree.
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.94835El crecimiento de un bosque aleatorio procede exactamente de la misma manera, excepto que usamos un valor mas pequeño del argumento mtry. De forma predeterminada, randomForest() u√s p / 3 variables cuando se crea un bosque aleatorio de arboles de regresion, y p variables cuando se crea un bosque aleatorio de arboles de clasificacion. Aqui usamos mtry = 6.
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.66454El conjunto de prueba MSE es 11.31; esto indica que los bosques aleatorios produjeron una mejora sobre el empaquetamiento en este caso. Usando la funcion importance(), podemos ver la importancia de cada variable.
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.73936Se pueden producir graficos de estas medidas de importancia utilizando la funcion varImpPlot().
varImpPlot (rf.boston )
Los resultados indican que en todos los arboles considerados en el bosque aleatorio, el nivel de riqueza de la comunidad (lstat) y el tamaño de la casa (rm) son, con mucho, las dos variables mas importantes.
8.3.4 Boosting
Aquí usamos el paquete gbm, y dentro de el la funcion gbm(), para ajustar los arboles de regresion aumentados al conjunto de datos de Boston. Ejecutamos gbm() con la opcion distribution = “gaussian” ya que este es un problema de regresion; Si se tratara de un problema de clasificacion binaria, utilizaríamos distribution = “bernoulli”.
library (gbm)Loaded gbm 2.1.5set.seed (1)
boost.boston =gbm(medv∼.,data=Boston [train ,], distribution= "gaussian",n.trees =5000 , interaction.depth =4)La funcion summary() produce un grafico de influencia relativa y tambien genera las estadisticas de influencia relativa.
summary (boost.boston )
Vemos que lstat y rm son las variables más importantes. Podemos Tambien se producir parcelas de dependencia parcial para estas dos variables. En este caso, como cabria esperar, la mediana. Los precios de las casas aumentan con rm y disminuyen con lstat.
par(mfrow =c(1,2))
plot(boost.boston ,i="rm")
plot(boost.boston ,i="lstat")
Ahora usamos el modelo reforzado para predecir medv en el conjunto de prueba:
yhat.boost=predict (boost.boston ,newdata =Boston [-train ,],
n.trees =5000)
mean(( yhat.boost -boston.test)^2)[1] 10.81479La prueba MSE obtenida es 11.8; Similar al MSE de prueba para bosques aleatorios y superior al de ensacado. Si queremos, podemos realizar un aumento con un valor diferente del parametro de contraccion λ en (8.10). El valor predeterminado es 0.001, pero esto se modifica facilmente. Aquí tomamos λ = 0.2.
boost.boston =gbm(medv∼.,data=Boston [train ,], distribution="gaussian",n.trees =5000 , interaction.depth =4, shrinkage =0.2,
verbose =F)
yhat.boost=predict (boost.boston ,newdata =Boston [-train ,],
n.trees =5000)
mean(( yhat.boost -boston.test)^2)[1] 11.51109