Lab 8. Árboles de decisión
1. Ajuste de árboles de clasificación
La libreria tree se utiliza para construir árboles de clasificación y regresión.
library(tree)Primero usamos árboles de clasificación para analizar el conjunto de datos Carseats. En estos datos, Sales es una variable continua, por lo que comenzamos por recodificarla como variable binaria. Usamos la función ifelse() para crear una variable, llamada High, que adquiere un valor de Yes si la variable Sales excede 8, y toma un valor de No de lo contrario.
library(ISLR)
attach(Carseats)
View(Carseats)
summary(Carseats)## Sales CompPrice Income Advertising
## Min. : 0.000 Min. : 77 Min. : 21.00 Min. : 0.000
## 1st Qu.: 5.390 1st Qu.:115 1st Qu.: 42.75 1st Qu.: 0.000
## Median : 7.490 Median :125 Median : 69.00 Median : 5.000
## Mean : 7.496 Mean :125 Mean : 68.66 Mean : 6.635
## 3rd Qu.: 9.320 3rd Qu.:135 3rd Qu.: 91.00 3rd Qu.:12.000
## Max. :16.270 Max. :175 Max. :120.00 Max. :29.000
## Population Price ShelveLoc Age Education
## Min. : 10.0 Min. : 24.0 Bad : 96 Min. :25.00 Min. :10.0
## 1st Qu.:139.0 1st Qu.:100.0 Good : 85 1st Qu.:39.75 1st Qu.:12.0
## Median :272.0 Median :117.0 Medium:219 Median :54.50 Median :14.0
## Mean :264.8 Mean :115.8 Mean :53.32 Mean :13.9
## 3rd Qu.:398.5 3rd Qu.:131.0 3rd Qu.:66.00 3rd Qu.:16.0
## Max. :509.0 Max. :191.0 Max. :80.00 Max. :18.0
## Urban US
## No :118 No :142
## Yes:282 Yes:258
##
##
##
## High<-ifelse(Sales <=8,"No","Yes")Finalmente, usamos la función data.frame() para fusionar High con el resto datos Carseats y utilizamos la función as.factor para cambiar la variable de tipo character a factor. Un factor es un tipo de datos específico a R. Puede ser descrito como un dato numérico representado por una etiqueta, en este caso Yes y No.
Carseats<-data.frame(Carseats,High)
class(Carseats$High)## [1] "character"Carseats$High <- as.factor(Carseats$High)
class(Carseats$High)## [1] "factor"Ahora usamos la función tree() para ajustar un árbol de clasificación para predecir High usando todas las variables menos Sales. La sintaxis de la función tree() es bastante similar a la de la función lineal lm().
tree.carseats<-tree(High~.-Sales, data=Carseats)
#se quitanlas ventas porque ya estan incluidas en highLa función summary () enumera las variables que se usan como nodos internos en el árbol, el número 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"
## [6] "Advertising" "Age" "US"
## Number of terminal nodes: 27
## Residual mean deviance: 0.4575 = 170.7 / 373
## Misclassification error rate: 0.09 = 36 / 400Vemos que la tasa de error de entrenamiento es del 9%. Una de las propiedades más atractivas de los árboles es que se pueden mostrar gráficamente. Usamos la función plot () para mostrar la estructura del árbol y la función text () para mostrar las etiquetas de los nodos. El argumento pretty = 0 le indica a R que incluye los nombres de categoría para cualquier predictador cualitativo, en lugar de simplemente mostrar una letra para cada categoría.
plot(tree.carseats)
text(tree.carseats, pretty=0, cex=0.5)
El indicador más importante de Sales parece ser la ubicación de ShelveLoc, ya que la primera rama diferencia Good ubicaciones de las Ubicaciones de Bad y Medium. Si simplemente escribimos el nombre del objeto del árbol, R imprime la salida correspondiente a cada rama del árbol. R Muestra el criterio de división (por ejemplo, Price <92.5), luego el número de observaciones en esa rama, la desviación, la predicción general de la rama (Sí o No) y la fracción de observaciones en esa rama que toman valores de Sí y No. Las ramas que conducen a nodos terminales se indican mediante asteriscos.
tree.carseats## node), 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 ) *Para evaluar adecuadamente el rendimiento de un árbol de clasificación en estos datos, debemos estimar el error de prueba en lugar de simplemente calcular el error de entrenamiento.
Lo que tenemos que hacer es divir las observaciones en un conjunto de entrenamiento y un conjunto de prueba, como ya lo hemos hecho en otros métodos, construimos el árbol usando el conjunto de entrenamiento y evaluamos su desempeño en los datos de la prueba. La función predict() se puede usar para este propósito. En el caso de un árbol de clasificación, el argumento type=“class” instruye a R para que devuelva la predicción de clase real.
set.seed(123)
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 87 13
## Yes 35 65(87+65)/200## [1] 0.76El 76% de las predicciones se realizaron correctamente. A continuación, consideramos si podar el árbol podría conducir a mejores resultados.
Para este fin la función cv.tree() realiza una validación cruzada para que determine el nivel óptimo de complejidad del árbol; Se utiliza poda de complejidad de costos para seleccionar una secuencia de árboles para su consideración. Usamos el argumento FUN=prune.misclass para indicar que queremos que la tasa de error de clasificación guíe el proceso de validación cruzada y poda. La función cv.tree() informa el número de nodos terminales de cada árbol considerado (tamaño), así como la tasa de error correspondiente y el valor del parámetro de costo-complejidad utilizado.
set.seed(123)
cv.carseats<-cv.tree(tree.carseats,FUN=prune.misclass)
#usamos misclas para...
names(cv.carseats)## [1] "size" "dev" "k" "method"cv.carseats## $size
## [1] 20 18 16 13 8 3 2 1
##
## $dev
## [1] 74 75 75 74 75 75 91 93
##
## $k
## [1] -Inf 0.000000 0.500000 1.333333 2.000000 3.000000 13.000000
## [8] 20.000000
##
## $method
## [1] "misclass"
##
## attr(,"class")
## [1] "prune" "tree.sequence"Noten que dev corresponde a la tasa de error de validación cruzada. ¿Cuál es el nodo del árbol con el valor mínimo de tasa de error? Trazamos la tasa de error en función del tamaño.
plot(cv.carseats$size ,cv.carseats$dev ,type="b")
Aplicamos la función prune.misclass() para podar el árbol.
prune.carseats<-prune.misclass(tree.carseats,best=5)
plot(prune.carseats)
text(prune.carseats,pretty=0)
#si se cumple es para la izquierda
#si no se cumple es para la derecha¿Qué tan bien funciona este árbol podado en el conjunto de datos de prueba? Una vez más, aplicamos la función predict ()
tree.pred<-predict(prune.carseats,Carseats.test,type="class")
table(tree.pred,High.test)## High.test
## tree.pred No Yes
## No 88 12
## Yes 34 66#3
(85+58)/200## [1] 0.715#8
(88+66)/200 #77% entonces es el que tiene mas precision## [1] 0.77#13
(86+65)/200## [1] 0.755Ahora el 77% las observaciones de prueba están clasificadas correctamente, por lo el proceso de poda produjo un árbol más interpretable y mejoró un poco la precisión de la clasificación.
Si aumentamos el valor best, obtenemos un árbol podado más grande con menor precisión de clasificación:
prune.carseats<-prune.misclass(tree.carseats, best=13)
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 13
## Yes 36 65(86+65)/200## [1] 0.7552. Ajuste de árboles de regresión
Aquí ajustamos un árbol de regresión al conjunto de datos de Boston. Primero, creamos un conjunto de entrenamiento y ajustamos el árbol a los datos de entrenamiento
library(MASS)
library(tree)
set.seed(123)
train <- sample(1:nrow(Boston), nrow(Boston)/2)
Boston.data.test<-Boston[-train, ]
tree.boston<-tree(medv~., data=Boston, subset=train) #el valor medio de las casas: medv
summary(tree.boston)##
## Regression tree:
## tree(formula = medv ~ ., data = Boston, subset = train)
## Variables actually used in tree construction:
## [1] "rm" "lstat" "crim" "tax"
## Number of terminal nodes: 8
## Residual mean deviance: 12.98 = 3179 / 245
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -11.86000 -1.99500 0.03714 0.00000 1.86800 16.45000Observe que la salida de summary () indica que solo cuatro de las variables se han utilizado para construir el árbol. En el contexto de un árbol de regresión, la desviación es simplemente la suma de los errores al cuadrado del árbol. Ahora trazamos el árbol.
plot(tree.boston)
text(tree.boston, pretty=0)
La variable rm representa el número medio de habitaciones por vivienda, lstat mide el porcentaje de individuos con un estado socioeconómico más bajo. El árbol indica que entre mayor sea el número de habitaciones las casas son más caras, así como los valores más bajos de lstat corresponden a las casas más caras. Ahora utilizamos la función cv.tree() para ver si podar el árbol mejorará el rendimiento
cv.boston<-cv.tree(tree.boston)
plot(cv.boston$size, cv.boston$dev, type="b")
En este caso, el árbol más complejo se selecciona mediante validación cruzada. Sin embargo, si deseamos podar el árbol, podríamos hacerlo de la siguiente manera, utilizando la función prune.tree ():
prune.boston<-prune.tree(tree.boston, best=4)
plot(prune.boston)
text(prune.boston, pretty=0)
De acuerdo con los resultados de la validación cruzada, utilizamos las predicciones de árbol sin podar para hacer predicciones en el conjunto de prueba
tree.boston## node), split, n, deviance, yval
## * denotes terminal node
##
## 1) root 253 20510.0 22.55
## 2) rm < 6.941 219 8862.0 20.13
## 4) lstat < 14.805 141 3385.0 23.37
## 8) lstat < 9.95 82 2169.0 25.30
## 16) crim < 1.38393 77 908.5 24.66
## 32) rm < 6.594 55 341.3 23.23 *
## 33) rm > 6.594 22 173.6 28.23 *
## 17) crim > 1.38393 5 743.0 35.16 *
## 9) lstat > 9.95 59 488.6 20.69 *
## 5) lstat > 14.805 78 1313.0 14.27
## 10) tax < 551.5 35 273.4 17.46 *
## 11) tax > 551.5 43 390.3 11.66 *
## 3) rm > 6.941 34 2085.0 38.15
## 6) rm < 7.437 22 607.1 33.55 *
## 7) rm > 7.437 12 161.7 46.58 *boston.pred<-predict(tree.boston, newdata=Boston.data.test)
boston.test<-Boston[-train,"medv"]
sqrt(mean((boston.pred-boston.test)^2))## [1] 4.848972La raíz cuadrada del MSE, indica que este modelo conduce a predicciones de prueba que están dentro del valor medio real de la vivienda para el suburbio.
3. Bagging and Random Forests
Aquí aplicamos Bagging y Random Forests a los datos Boston, usando la librería randomForest en R. Los resultados exactos obtenidos en esta sección pueden depender de la versión de R y de la versión del paquete randomForest() instalado en su computadora. Recuerde que el Bagging es simplemente un caso especial de un Random Forests con \(m = p\). Por lo tanto, la función randomForest () se puede utilizar para realizar Random Forests y Bagging. Realizamos Bagging como sigue
library(randomForest)## randomForest 4.7-1.1## Type rfNews() to see new features/changes/bug fixes.set.seed(123)
bag.boston<-randomForest(medv~.,data=Boston,subset=train,mtry=13,importance=TRUE)
#bagging vs random forest
#investigar alae
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: 13.59049
## % Var explained: 83.23El argumento mtry= 13 indica que los 13 predictores deben ser considerados para cada división del árbol. ¿Qué tan bien se desempeña este modelo sobre el set de prueba?
boston.pred.bag <- predict(bag.boston, newdata=Boston.data.test)
plot(boston.pred.bag, boston.test)+
abline(0,1)
## integer(0)mean((boston.pred.bag-boston.test)^2)## [1] 14.82658¿Cual es el MSE asociado con el árbol de regresión en bagging? Podríamos cambiar el número de árboles cultivados por RandomForest() usando el argumento siguiente ntree:
set.seed(123)
bag.boston<-randomForest(medv~.,data=Boston,subset=train,mtry=13,ntree=25)
boston.pred.bag <- predict(bag.boston, newdata=Boston.data.test)
mean((boston.pred.bag-boston.test)^2)## [1] 14.55574El crecimiento de un Random Forests procede exactamente de la misma manera, excepto que usamos un valor menor de su argumento. Por defecto, randomForest() usa variables \(p / 3\) cuando construye un bosque aleatorio de árboles de regresión y \(\sqrt{p}\) variables cuando construye un bosque aleatorio de árboles de clasificación. Aquí usamos mtry = 6.
set.seed(123)
rf.boston<-randomForest(medv~.,data=Boston,subset=train,mtry=6,importance =TRUE)
boston.pred.rf <- predict(rf.boston ,newdata=Boston.data.test)
mean((boston.pred.rf-boston.test)^2)## [1] 12.76183¿Cual es el valor del conjunto de prueba MSE? ¿se mejoró?; Usando la función importance (), podemos ver la importancia de cada variable y se pueden generar gráficos de estas medidas importantes utilizando la función varImpPlot().
importance(rf.boston)## %IncMSE IncNodePurity
## crim 12.697705 1439.85189
## zn 1.948040 54.09570
## indus 6.850607 596.16940
## chas 2.404661 75.58075
## nox 10.815893 703.32507
## rm 32.755672 7305.23476
## age 10.222574 515.34594
## dis 13.608222 1112.04788
## rad 4.600246 131.35233
## tax 9.239586 508.45391
## ptratio 9.690679 717.36922
## black 5.722377 301.14534
## lstat 29.112077 6743.11775varImpPlot(rf.boston)
Se muestran dos medidas de variables de importancia. El primero se basa en la disminución media de la precisión en las predicciones de las muestras fuera de bagging cuando una variable dada se excluye del modelo. Esta última es una medida de la disminución total en la impureza de los nodos que resulta de las divisiones sobre esa variable, promediada sobre todos los árboles. En el caso de los árboles de regresión, la impureza del nodo se mide por el RSS de entrenamiento, y para los árboles de clasificación por la desviación.
Los resultados indican que en todos los árboles considerados en el bosque aleatorio, el nivel de riqueza de la comunidad (lstat) y el tamaño de la casa (rm) son por mucho las dos variables más importantes para el valor medio de la casa.
4. Boosting
Aquí usamos la libreria gbm, y dentro de la función gbm(), para ajustar los árboles de regresión al conjunto de datos Boston. Ejecutamos gbm() con la opción de distribution=“gaussiano” ya que este es un problema de regresión; si se tratara de un problema de clasificación binaria, utilizaríamos distribution=“bernoulli”. El argumento n.trees=5000 indica que queremos 5000 árboles, y la opción interaction.depth=4 limita la profundidad de cada árbol. Tome en cuenta que para la función gbm no es posible utilizar el argumento subset=train, así que desde el argumento del conjunto de datos hay que seleccionar sólo las observaciones de entrenamiento.
library(gbm)## Loaded gbm 2.1.9## This version of gbm is no longer under development. Consider transitioning to gbm3, https://github.com/gbm-developers/gbm3set.seed(123)
boost.boston<-gbm(medv~.,data=Boston[train,],
distribution="gaussian",n.trees=5000,interaction.depth=4)La función summary () produce un gráfico y las estadísticas de influencia relativa.
summary(boost.boston)
## var rel.inf
## lstat lstat 36.9047193
## rm rm 31.6492543
## dis dis 9.9236656
## crim crim 4.7636681
## age age 3.1727355
## ptratio ptratio 2.8884337
## black black 2.8779828
## tax tax 2.8055114
## nox nox 2.2500602
## indus indus 1.0427958
## chas chas 0.8189095
## rad rad 0.7834895
## zn zn 0.1187742Vemos que lstat y rm son por mucho las variables más importantes. También podemos producir gráficos de dependencia parcial para estas dos variables. Estas gráficas de dependencia parcial ilustran el efecto marginal de las variables seleccionadas sobre la respuesta después de integrar las otras variables. En este caso, como podríamos esperar, los precios medios de las casas están aumentando con rm y disminuyendo con lstat.
par(mfrow=c(1,2))
plot(boost.boston, i="rm")
plot(boost.boston, i="lstat")
5. Resuelve
Este problema tiene que ver con el conjunto de datos OJ que forma parte del paquete ISLR.
Crea un conjunto de entrenamiento que contenga una muestra aleatoria de 800 observaciones y un conjunto de pruebas que contenga las observaciones restantes. Utilice una semilla de 555.
Ajusta un árbol a los datos de entrenamiento, con Purchase como la respuesta y las otras variables como predictores. Utilice la función summary() para producir estadísticas resumidas sobre el árbol, y describa los resultados obtenidos. ¿Cuál es la tasa de error de entrenamiento? ¿Cuántos nodos terminales tiene el árbol?
Crea una gráfica del árbol, e interpreta los resultados.
Predice la respuesta en los datos de la prueba, y producir una matriz de confusión comparando las etiquetas de la prueba con las etiquetas de la prueba pronosticada. ¿Cuál es la tasa de error de la prueba?
Aplica la función cv.tree() al conjunto de entrenamiento para determinar el tamaño óptimo del árbol.
Produce un gráfico con el tamaño del árbol en el eje x y la tasa de error de clasificación cruzada en el eje y.
¿Qué tamaño de árbol corresponde a la menor tasa de error de clasificación validada cruzada?
Produce un árbol podado que corresponda al tamaño óptimo de árbol obtenido mediante validación cruzada. Si la validación cruzada no conduce a la selección de un árbol podado, crea un árbol podado con cinco nodos terminales.