library(MASS)
library(randomForest)
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
library(ggplot2)
## 
## Attaching package: 'ggplot2'
## The following object is masked from 'package:randomForest':
## 
##     margin
library(tibble)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:randomForest':
## 
##     combine
## The following object is masked from 'package:MASS':
## 
##     select
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
data("Boston")
index <- sample(nrow(Boston),nrow(Boston)*0.70) # separar train/test
boston.train <- Boston[index,]
boston.test <- Boston[-index,]
rm(Boston)

El número de columnas es ncol(Boston)

El data set tiene 14 columnas: 1 variables respuesta y 13 predictores.

Random Forest

Random Forest es un derivado del Bagging, y mejora significativamente la predicción. La idea es que, aleatoriamente, se selecciona un subconjunto de predictores como variables candidatas en cada separación del árbol. La razón de esto es adaptarse a las correlaciónes naturales entre las variables y los árboles, así se reduce la varianza al agregar árboles.

Uso de Random Forest

library(randomForest)
boston.rf<- randomForest(medv~., data = boston.train, importance=TRUE)
boston.rf
## 
## Call:
##  randomForest(formula = medv ~ ., data = boston.train, importance = TRUE) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 4
## 
##           Mean of squared residuals: 12.59328
##                     % Var explained: 85.11

Por defecto, el modelo de árbol de regresión utiliza un tercio de las variables como predictoras en cada iteración. Mientras que utiliza la raíz de las variables para las el problema de clasificación poor defecto. En ambos casos se puede cambiar modificando la opción mtry=. También se puede especificar el número de árboles con ntree=. por defecto se utilizan 500. El argumento importance=TRUE nos permite ver la importancia de cada variable en el modelo.

tabla<- rownames_to_column(data.frame(boston.rf$importance),"Variable")
tabla
##    Variable   X.IncMSE IncNodePurity
## 1      crim  9.9613927     1935.2968
## 2        zn  0.9728210      239.4470
## 3     indus  5.9674937     1671.9717
## 4      chas  0.2550211      118.8992
## 5       nox  9.5223864     1946.2755
## 6        rm 32.7600070     8143.5607
## 7       age  3.8001590      795.7722
## 8       dis  6.4876023     1756.4042
## 9       rad  1.8822130      325.2100
## 10      tax  4.7348743     1133.1837
## 11  ptratio  8.3986984     2444.9074
## 12    black  1.7853023      630.8243
## 13    lstat 60.0653950     8206.1475
tabla %>% ggplot(aes(y=X.IncMSE,x=Variable))+geom_bar(stat="identity")+
  coord_flip()+theme_classic()+ 
    geom_text(aes(x=Variable,y=X.IncMSE,label=X.IncMSE),hjust=0)

El Random Forest guarda todos los errores del OOB para cada ntree desde 1 a 500. podemos graficarlo para ver cómo el error OOB cambia según el número de árboles ntree.

plot(boston.rf$mse, type='l', col=2, lwd=2, xlab = "ntree", ylab = "OOB Error")

Podemos ver el error al aplicar el modelo en la base test

boston.rf.pred<- predict(boston.rf, boston.test)

mean((boston.test$medv-boston.rf.pred)^2)
## [1] 8.454668
rm(boston.rf)

Como se dijo previamente, el número de variables candidatas en cada iteración o separación es la raíz de las variables totales: Si son 13 variables, m~4. También podemos especificar este número en el argumento mtry=. Ahora podemos graficar las diferencias del error en OOB y el testing y como éstos cambian según el argumento mtry=

oob.err<- rep(0, 13)
test.err<- rep(0, 13)
for(i in 1:13){
  fit<- randomForest(medv~., data = boston.train, mtry=i)
  oob.err[i]<- fit$mse[500]
  test.err[i]<- mean((boston.test$medv-predict(fit, boston.test))^2)
  cat(i, " ")
}
## 1  2  3  4  5  6  7  8  9  10  11  12  13

En rojo: test error En azul: OOB error

matplot(cbind(test.err, oob.err), pch=15, col = c("red", "blue"), type = "b", ylab = "MSE", xlab = "mtry")