Para construir árbles de regresión y clasificación utilizaremos la libreria tree() acompñado del dataset Carseats para predecir la variable sales, que es continua, pero la convertiremos a forma binaria en la variable High indicado Yes si sales es mayor a 8 y luego agregamos esta nueva varible al dataset original.

library(tree)
library(ISLR)
attach(Carseats)
High=ifelse (Sales <=8,"No","Yes")

Carseats =data.frame(Carseats ,High)
tree.carseats =tree(High~.-Sales ,Carseats )

luego con la función summary() veremos las variables utilizadas para calcular los nodos, la cantidad de nodos terminales y el ratio de error del árbol.

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 / 400

Una de las características más atractivas de los árboles es que podes mostrarlos graficamente, podemos graficar la estructura con plot() y mostrar las etiquetas con text().

plot(tree.carseats )
text(tree.carseats ,pretty =0)

Parece que la variable mas importante para Sales es Shelving location, ya que desde la primer rama distingue las ubicaciones.

Si solo llamamos al arbol por su nombre nos mostrará el criterio que se utiliza para sus cortes, la etiqueta que asignará cada rama y la cantidad de registros que encajan.

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 poder evaluar el árbolo no podemos hacerlo solo con un set de datos, entonces lo dividiremos en prueba y entrenamiento como lo hemos hecho anteriormente y comprabaremos la exactitud del modelo, para esto usaremos la funcion predict() junto con el parámetro type=“class” que le indica que nos devuelva las etiquetas de clsificación. Esto nos devuelve un 71.5% de predicciones correctas.

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  104  33
##       Yes  13  50

ahora probaremos podar el árbol porque consideramos que eso puede arrojar mejores resultados. La función cv.tree() realiza una validación cruzada para buscar un equilibrio de entre complejidad y costos. El parámetro FUN=prune.misclass indica que queremos usar el ratio de error de clasificación para guiar la validación cruzada y la poda en lugar del parámetro estándar de la función cv.tree() que es la desviación. La función nos devuelve la cantidad de nodos terminales que es size, el error correspondiente, y el valor de costo-complejidad.

set.seed (3)
cv.carseats = cv.tree(tree.carseats , FUN = prune.misclass)
names(cv.carseats)
## [1] "size"   "dev"    "k"      "method"

tomemeos en cuenta que a pesar que el parámetro devuelto se llama dev, en realidad contiene el ratio de error de la validación cruzada. El árbol con 9 nodos terminales es el de menor error.

par(mfrow = c(1, 2))
plot(cv.carseats$size , cv.carseats$dev , type ="b")
plot(cv.carseats$k , cv.carseats$dev , type ="b")

ahora usamos la función prune.misclass() para podar el árbol

prune.carseats =prune.misclass (tree.carseats ,best =9)
plot(prune.carseats )
text(prune.carseats ,pretty =0)

ahora utilizemos el set de prueba

tree.pred = predict (prune.carseats , Carseats.test , type ="class")
table(tree.pred , High.test)
##          High.test
## tree.pred No Yes
##       No  97  25
##       Yes 20  58

conseguimos el 77% de exactitud, quiere decir que la poda no solo permite comprender mejor el modelo sino que también mejoró la exactitud.

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  102  30
##       Yes  15  53

8.3.2 Fitting Regression Trees

Acá usaremos el dataset Boston. Primero divimos los datos y ajustamos el modelo.

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] "rm"    "lstat" "crim"  "age"  
## Number of terminal nodes:  7 
## Residual mean deviance:  10.38 = 2555 / 246 
## Distribution of residuals:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -10.1800  -1.7770  -0.1775   0.0000   1.9230  16.5800

notemos que solo 3 variables han sido utilizadas para construir el árbol. Veamoslo de forma gráfica.

plot(tree.boston )
text(tree.boston ,pretty =0)

la variable lstat mide el porcenaje de individuos con bajo nivel socioeconómico. El modelo dice que a menor nivel socioeconomico poseen casas mas cara. El árbol arroja una media de $46,400 para casas grandes en los suburbios donde los residentes tienen un nivel socioeconómico entre rm>=7.437 y lstat<9.715.

Probemos podar y utilizar validación cruzada para elegir el mejor resultado.

cv.boston =cv.tree(tree.boston )
plot(cv.boston$size ,cv.boston$dev ,type='b')

prune.boston =prune.tree(tree.boston ,best =5)
plot(prune.boston )
text(prune.boston ,pretty =0)

el árbol mas complejo fue seleccionado y la poda no fue tan eficiente por lo que nos quedamos con el resultado de la validación cruzada

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] 35.28688

El error mínimo cuadrado del árbol es de 25.05. La raíz cuadrada de este está alrededor de 5.005, indicando que el modelo arroja resultados de $5,005 alrededor de la media verdadera de hogares.

8.3.3 Bagging and Random Forests

Acá usaremos el paquete randomForest. Los resultados acá pueden variar según la versión de R y del paquete que estemos utilizando. Bagging es un caso especial de árbol aleatorio con m = p y usaremos la función randomForest() para realizar ambos.

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.39601
##                     % Var explained: 85.17

el argumento mtry=13 indica que todos los predictores deben ser utilizados, es decir que se utilice bagging

yhat.bag = predict (bag.boston ,newdata =Boston [-train ,])
plot(yhat.bag , boston.test)
abline (0,1)

mean(( yhat.bag -boston.test)^2)
## [1] 23.59273

el error mínimo cuadrado es de 13.16, que es casi la mitad del anterior. Podemos cambiar la cantidad de árboles a usar con el parámetro 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] 23.66716

hacer crecer un árbol aleatorio funciona en la misma forma pero con un mtry más pequeño, por defecto la función randomForest() utiliza \(p/3\) para los de regresion y \(\sqrt{p}\) cuando es de clasificacion. Nosotros usaremos 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] 19.62021

el error mínimo cuadrado bajó a 11.31, es decir mejoró. Con importance() podemos ver la importancia de cada predictor

importance (rf.boston )
##           %IncMSE IncNodePurity
## crim    16.697017    1076.08786
## zn       3.625784      88.35342
## indus    4.968621     609.53356
## chas     1.061432      52.21793
## nox     13.518179     709.87339
## rm      32.343305    7857.65451
## age     13.272498     612.21424
## dis      9.032477     714.94674
## rad      2.878434      95.80598
## tax      9.118801     364.92479
## ptratio  8.467062     823.93341
## black    7.579482     275.62272
## lstat   27.129817    6027.63740
varImpPlot (rf.boston )

Boosting

Acá usaremos el paquete gbm

library (gbm)
## Loaded gbm 2.1.5
set.seed (1)
boost.boston = gbm(
  medv~.,
  data = Boston [train , ],
  distribution =
   "gaussian",
  n.trees = 5000 ,
  interaction.depth = 4
)

summary() nos devuelve la influencia de las variables en el árbol

summary (boost.boston )

vemos que lstat y rm son las variables mas importantes.

par(mfrow =c(1,2))
plot(boost.boston ,i="rm")

plot(boost.boston ,i="lstat")

ahora usamos boosting para predecir en el set de pruebas

yhat.boost = predict (boost.boston , newdata = Boston [-train ,], n.trees =
                        5000)
mean((yhat.boost - boston.test) ^ 2)
## [1] 18.84709

El error mínimo cuadrado es 11.8 que es similar al de árbol aleatorio y mejor que con bagging. Podemos intentar con un gamma = 0.2 para buscar mejores resultados.

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] 18.33455

el resultado si fue mejor.

libera memoria

rm()
gc()
##           used (Mb) gc trigger  (Mb) max used  (Mb)
## Ncells 1727021 92.3    2588025 138.3  2588025 138.3
## Vcells 4001755 30.6    8388608  64.0  6102030  46.6