1 Estimación de los árboles de regresión

library(tree)
datos = read.table('coches.txt',header=T)
t1 = tree(consumo ~ .,data=datos)
plot(t1)
text(t1, cex = 0.75)

summary(t1)
## 
## Regression tree:
## tree(formula = consumo ~ ., data = datos)
## Variables actually used in tree construction:
## [1] "cc"   "cv"   "ano"  "peso"
## Number of terminal nodes:  10 
## Residual mean deviance:  2.142 = 816 / 381 
## Distribution of residuals:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -6.50000 -0.95830 -0.04938  0.00000  0.93480  5.50000

Los regresores que son influyentes en el consumo son los que se han empleado en las particiones.

print(t1)
## node), split, n, deviance, yval
##       * denotes terminal node
## 
##  1) root 391 5911.00 11.230  
##    2) cc < 3482 234  856.20  8.709  
##      4) cv < 84.5 129  217.00  7.659  
##        8) ano < 76.5 48   54.31  8.688 *
##        9) ano > 76.5 81   81.80  7.049 *
##      5) cv > 84.5 105  322.00 10.000  
##       10) ano < 78.5 72  169.10 10.610  
##         20) peso < 921 48   69.92  9.958 *
##         21) peso > 921 24   37.83 11.920 *
##       11) ano > 78.5 33   67.33  8.667 *
##    3) cc > 3482 157 1359.00 14.980  
##      6) cv < 143.5 82  358.00 13.110  
##       12) peso < 1242 60  185.00 12.500 *
##       13) peso > 1242 22   89.86 14.770 *
##      7) cv > 143.5 75  399.90 17.030  
##       14) peso < 1453.5 46  132.80 16.070 *
##       15) peso > 1453.5 29  157.20 18.550  
##         30) ano < 73.5 22   85.09 19.360 *
##         31) ano > 73.5 7   12.00 16.000 *

Lo que devuelve la tabla es:

  • Numero del nodo
  • split: criterio para hacer la partición del nodo
  • n: numero de datos que componen el nodo.
  • deviance: = RSS = \(\sum(y_i - \hat{y}_i)^2\).
n = nrow(datos)
( deviance_root = sum( (datos$consumo - mean(datos$consumo))^2 ) )
## [1] 5910.742
  • yval: predicción del nodo = \(\bar{y}\)
mean(datos$consumo)
## [1] 11.22762
  • un asterisco * para indicar un nodo terminal u hoja.

2 Parámetros del árbol

t2 = tree(consumo ~ .,data=datos, control = tree.control(nobs=nrow(datos),mincut = 20,minsize = 40,mindev = 0.05))
plot(t2)
text(t2, cex=.75)

control:

  • minsize: tamaño mínimo del nodo para que se divida en dos (por defecto, minsize = 10).
  • mincut: si al dividir un nodo, uno de los nodos hijo tiene un tamaño inferior a éste, no se divide el nodo (por defecto, mincut = 5). Ojo, mincut ha de ser menor que minsize/2.
  • mindev: para que un nodo se divida, la deviance del nodo tiene que ser mayor que mindev por la deviance del nodo raiz.
print(t2)
## node), split, n, deviance, yval
##       * denotes terminal node
## 
## 1) root 391 5911.0 11.230  
##   2) cc < 3482 234  856.2  8.709  
##     4) cv < 84.5 129  217.0  7.659 *
##     5) cv > 84.5 105  322.0 10.000 *
##   3) cc > 3482 157 1359.0 14.980  
##     6) cv < 143.5 82  358.0 13.110 *
##     7) cv > 143.5 75  399.9 17.030 *
  • si queremos obtener un árbol que se ajusta perféctamente a los datos hay que utilizar mindev = 0 y minsize = 2. Estos árboles predicen muy mal porque son muy sensibles a ligeros cambios en los datos.
t3 = tree(consumo ~ ., data=datos, control = tree.control(nobs=nrow(datos),minsize = 2,mindev = 0.0))
plot(t3)

#text(t3, cex=.75)

3 Podado

t3_prune = prune.tree(t3, best = 8)
plot(t3_prune)
text(t3_prune, cex=.75)

4 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(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

Entrenamiento del modelo

t3 <- tree(consumo ~ ., data = datos_train)
plot(t3)
text(t3)

4.1 Error del modelo

y_train_p <- predict(t3, newdata = datos_train)
y_train = datos_train$consumo
plot(y_train,y_train_p)
abline(0,1) # como pintamos y vs yp, la relacion perfecta deberia ser una recta a 45º (m=1)

( MSE_train = mean((y_train-y_train_p)^2) ) # error cuadratico medio en los datos de training
## [1] 1.703067

4.2 Error de predicción

# prediccion del consumo con los datos de test
y_test_p = predict(t3, newdata = datos_test)

# error del test set
y_test = datos_test$consumo
(MSE_test = mean((y_test-y_test_p)^2))
## [1] 3.114328

El error es mucho mayor con los datos de validación!