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.

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

Boosting

Consiste en ajustar secuencialmente múltiples modelos sencillos, llamados weak learners, de forma que cada modelo aprende de los errores del anterior. Como valor final, al igual que en bagging, se toma la media de todas las predicciones (variables continuas) o la clase más frecuente (variables cualitativas). Tres de los métodos de boosting más empleados son AdaBoost, Gradient Boosting y * Stochastic Gradient Boosting*.

library(gbm)
## Loaded gbm 2.1.5
boston.boost<- gbm(medv~., data = boston.train, distribution = "gaussian",
                   n.trees = 10000, shrinkage = 0.01, interaction.depth = 8)
summary(boston.boost)

##             var    rel.inf
## lstat     lstat 42.4684267
## rm           rm 25.0005083
## dis         dis  9.9861621
## nox         nox  4.1284482
## age         age  3.9688560
## crim       crim  3.6989371
## black     black  3.0748729
## tax         tax  2.6329781
## ptratio ptratio  1.6973551
## indus     indus  1.3915215
## chas       chas  1.0904641
## rad         rad  0.7089181
## zn           zn  0.1525519

Notar que debemos especificar el tipo de distribución en distribution = “gaussian” al trabajar en árboles de regresión. Por defecto esto trabaja en distribución Bernoulli para clasificación binaria. n.trees corresponde a l número de árboles a ajustar. A mayor n.trees podría haber mayor sobreajuste. shrinkage = es un argumento para decidir cuanta contribución aporta cada nuevo árbol al modelo. Finalmente interaction.depth = es cuantas separaciones tendrá cada árbol. La mejor manera de ajustar estos parámetros es mediante Validación cruzada y grid search (que se verá en la próxima ayudantía).

Al aplicar en la base test

boston.boost.pred.test<- predict(boston.boost, boston.test, n.trees = 10000)
mean((boston.test$medv-boston.boost.pred.test)^2)
## [1] 9.723984

Notar que el error cambia dependiendo de la cantidad de árboles. Al igual que en los casos anteriores, una sobreproducción de árboles podría llevar a sobreajuste del modelo.

ntree<- seq(100, 10000, 100)
predmat<- predict(boston.boost, newdata = boston.test, n.trees = ntree)
err<- apply((predmat-boston.test$medv)^2, 2, mean)
plot(ntree, err, type = 'l', col=2, lwd=2, xlab = "n.trees", ylab = "Test MSE")
abline(h=min(test.err), lty=2)

Notar que la línea punteada muestra el mejor caso de error de Random Forest que se vio anteriormente.