Boosting

1 Lectura de datos

library(tree)
d = read.table('datos/coches.txt',header=T)
d = na.omit(d)
d$origen = factor(d$origen, labels = c("USA","Europa","Japon"))

2 Training set vs Test set

Dividimos los datos en dos partes, una para entrenar el modelo y otra para calcular el error de predicción con datos diferentes de los utilizados para entrenar:

set.seed(321)
ndat = nrow(d)
pos_train = sample(1:ndat,ndat/2) # la mitad de los datos para entranamiento
datos_train = d[pos_train,] # trainning set
datos_test = d[-pos_train,] # test set

3 Boosting

• Se puede utilizar para otros modelos además de los árboles.
• El modelo verifica que: \(y = \hat f(x) + r\).
• Empezar el algoritmo tomando \(\hat f_0(x)=0\) y \(r_0 = y\).
• For b = 1, …, B
  • Estimar el modelo correspondiente \(r_{b-1} = \hat f_b(x) + r_b\) (arbol de regresión en este caso).
  • Actualizar el modelo: \(\hat f_{b}(x) = \hat f_{b-1}(x) + \lambda \hat f_b(x)\)
  • Actualizar los residuos: \(r_b = r_{b-1} - \lambda r_b(x)\)
• El resultado es el modelo boosted será, por tanto:

\[ \hat f(x) = \sum _{b=1}^{B} \lambda \hat f_b(x), \quad r_b \rightarrow 0 \] - El parámetro \(\lambda\) controla la velociad del proceso.

library(gbm)

## Loaded gbm 2.1.5

set.seed(321)
boost1 = gbm(consumo ~ ., data = datos_train, distribution = "gaussian", 
             n.trees = 5000, interaction.depth = 4)

• interaction.depth = 4: Como mucho, cada árbol tendrá cuadro niveles (5 nodos terminales).

4 Importancia de variables

summary(boost1)

##                 var    rel.inf
## peso           peso 31.3580874
## cc               cc 30.3814004
## cv               cv 19.3792893
## ano             ano 11.4855761
## acel           acel  5.9314782
## origen       origen  0.8644751
## cilindros cilindros  0.5996935

5 Predicciones

y_train = datos_train$consumo
yp_train_boost1 = predict(boost1, newdata = datos_train, n.trees = 5000)
( MSE_train_boost1 = mean((y_train - yp_train_boost1)^2) )

## [1] 0.01935361

Es lógico, ya que los modelos con boosting verifican que \(r_i \rightarrow 0\).

y_test = datos_test$consumo
yp_test_boost1 = predict(boost1, newdata = datos_test, n.trees = 5000)
( MSE_test_boost1 = mean((y_test - yp_test_boost1)^2) )

## [1] 2.078316

Que es superior al obtenido con bagging y random forest.

Se puede modificar el parámetro \(\lambda\), que por defecto es 0.001

boost2 = gbm(consumo ~ ., data = datos_train, distribution = "gaussian", 
             n.trees = 5000, interaction.depth = 4, shrinkage = 0.2)

yp_train_boost2 = predict(boost2, newdata = datos_train, n.trees = 5000)
( MSE_train_boost2 = mean((y_train - yp_train_boost2)^2) )

## [1] 0.0102184

yp_test_boost2 = predict(boost2, newdata = datos_test, n.trees = 5000)
( MSE_test_boost2 = mean((y_test - yp_test_boost2)^2) )

## [1] 2.554347