Cap # 8, FIABILIDAD - Oscar Padilla

8.3.1 Fitting Classification Trees

La libreria Tree se utiliza para construir arboles de clasificacion y regresion.

library(tree)

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() 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 ")

Usamos la funcion data.frame() para unir High con el resto de los datos de Carseats.

Carseats =data.frame(Carseats ,High)

Usamos la funcion tree() para ajustar un arbol de clasificacion con el fin de predecir High usando todas las variables excepto Sale.

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.

summary (tree.carseats )

Classification tree:
tree(formula = High ~ . - Sales, data = Carseats)
Variables actually used in tree construction:
[1] "ShelveLoc"   "Price"       "Income"      "CompPrice"   "Population"  "Advertising"
[7] "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 propiedades mas atractivas de los arboles es que se pueden mostrar graficamente. Usamos la funcion plot() para mostrar la estructura del arbol, y la funcion text() para mostrar las etiquetas de nodo. La tasa de error de entrenamiento es del 9%. 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)

Se indican mediante asteriscos las ramas que conducen a los nodos terminales.

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 ) *

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

(86+57) /200

[1] 0.715

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.

set.seed (3)
cv.carseats =cv.tree(tree.carseats ,FUN=prune.misclass )
names(cv.carseats )

[1] "size"   "dev"    "k"      "method"

cv.carseats

$size
[1] 21 19 14  9  8  5  3  2  1

$dev
[1] 74 76 81 81 75 77 78 85 81

$k
[1] -Inf  0.0  1.0  1.4  2.0  3.0  4.0  9.0 18.0

$method
[1] "misclass"

attr(,"class")
[1] "prune"         "tree.sequence"

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)

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    97    25
     Yes   20    58

(94+60) /200

[1] 0.77

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

(86+62) /200

[1] 0.74

8.3.2 Fitting Regression Trees

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

En el contexto de un arbol de regresion, la desviacion es simplemente la suma de los errores cuadrados para el arbol.

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')

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)

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

La 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 Bagging and Random Forests

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.

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 se deben considerar los 13 predictores para cada division del arbol; en otras palabras, se debe realizar el empaquetamiento.

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

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

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, con 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

Usando la funcion importance(), podemos ver la importancia de cada variable.

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

Se pueden producir graficos de estas medidas de importancia utilizando la funcion varImpPlot().

varImpPlot (rf.boston )

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.5

set.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 )

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")

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

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, con λ = 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] 18.33455