library(tree)
datos = read.table('coches.txt',header=T)
datos = na.omit(datos)
datos$origen = factor(datos$origen, labels = c("USA","Europa","Japon"))
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(datos)
train = sample(1:ndat,ndat/2) # la mitad de los datos para entranamiento
datos_train = datos[train,] # trainning set
datos_test = datos[-train,] # test set
\[ \hat f(x) = \sum _{b=1}^{B} \lambda \hat f_b(x), \quad r_i \rightarrow 0 \] - El parámetro \(\lambda\) controla la velociad del proceso.
library(gbm)
## Loaded gbm 2.1.4
set.seed(321)
boost1 = gbm(consumo ~ ., data = datos_train, distribution = "gaussian", n.trees = 5000, interaction.depth = 4)
summary(boost1)
## var rel.inf
## peso peso 39.3648351
## cc cc 22.7062692
## cv cv 18.0349291
## ano ano 11.9272781
## acel acel 5.8688626
## cilindros cilindros 1.3570329
## origen origen 0.7407932
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.01896061
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.461189
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.004823553
yp_test_boost2 = predict(boost2, newdata = datos_test, n.trees = 5000)
( MSE_test_boost2 = mean((y_test - yp_test_boost2)^2) )
## [1] 2.610173