Tarea
Luis Saucedo
3 de mayo de 2019
Ejercicio 8
Buscar predecir la variable “Sales” usando árboles de regresión y aproximaciones relacionadas, tratando la variable respuesta como variable numérica.
- Dividir en conjunto de prueba y de entrenamiento.
set.seed(101010)
train=sample(1:nrow(Carseats), 300) #El libro usa la mitad para train, yo usaré el 75 % que es 300
Carseats.test<-Carseats[-train, "Sales"]
Carseats.train<-Carseats[train, "Sales"]El conjunto de prueba tiene el 75 % de las observaciones, que serían 300. El conjunto test tiene 100 observaciones.
- Ajustar un árbol de regresión al conjunto de entrenamiento. Graficar el árbol e interpretar los resultados. Reportar el MSE.
Ajusto el árbol con el conjunto de entrenamiento y obtengo un ábol de tamaño 15, que utiliza las variables: “price” indica el precio que cobra la compañía por los asientos de carro en cada lugar. “Age” es la edad promedio en la zona. “income” es el nivel de ingresos en la comunidad (en miles de dólares). “CompPrice” es el precio que le cobran al competidor por cada zona. Por último, la variable respuesta “Sales” es el número de ventas (en miles) por cada zona, por lo tanto las ramas de los árboles van a ser el promedio de unidades vendidas (en miles) por cada zona y por cada rama. El MSE obtenido es: 5.01
Carseats.train<-Carseats[train, "Sales"]
simple_tree<-tree(Sales~ ., data=Carseats, subset=train)
#plot(simple_tree); text(simple_tree, pretty = 0) #graficando árbol
###MSE Test
yhat<-predict(simple_tree, newdata=Carseats[-train,])#árbol, datos a predecir
mean((yhat-Carseats.test)*(yhat-Carseats.test)) #MSE 5.01 Test error
errors[1,2]<-mean((yhat-Carseats.test)^2)
#MSE Train
yhat_t<-predict(simple_tree, newdata=Carseats[train,]) #árbol, datos a predecir
mean((yhat_t-Carseats.train)^2) #MSE 2.39 Test error
errors[1,1]<-round(mean((yhat_t-Carseats.train)^2), 2)
Fig 1. Árbol de regresión de tamaño 15, usando la base de datos Carseats y tomando como variable respuesta la variable “Sales”.
Al ser un árbol de regresión es fácil interpretar los resultados. La primera división del árbol es sobre la variable “ShelveLoc”, que indica la calidad de ubicación de estanterías para los asientos de automóvil en cada sitio. La rama de la izquiera son las observaciones con calidad mala y/o media y la derecha son las de calidad buena.
Primero notar que donde hay más ventas es en donde el promedio es igual a 8, 9 ó 10. Donde hay menos ventas el promedio es: 3 y 4 y ventas regulares: 6 y 7 (en miles).
El lugar donde en promedio hay 3000 unidades de ventas, es decir con menos ventas, es donde la calidad es mala y el precio es mayor o igual que 143.5. Además la media de la variable “price” es: 116, entonces podría considerar que 143.5 es un precio alto. Así, tiene sentido que haya pocas ventas, ya que el precio es alto y la calidad para guardar los asientos es mala. Por otro lado, mala calidad, precio arriba de la media, y zona donde habitan personas grandes tiene, en promedio, 4000 ventas, también tiene sentido ya que la gente tiene una calidad media para guardar los asientos, pero es gente de mayor edad, por lo que casi no tienen hijos.
Veamos dónde hay más ventas. Precio menor que 106.5, edad menor que 60.5 y calidad media es donde hay 10000 en promedio. Son asientos baratos, gente joven que tiene calidad media para guardarlos. Por otro lado, si la edad es mayor o igual que 60.5, el precio menor que 80.5 y calidad media y/o baja hay muchas ventas, esto dice que la gente mayor compra asientos baratos aunque tengan un lugar de calidad media o baja para guardarlos.
Viendo ahora la primer rama derecha del árbol, es la calidad buena para guardarlos; donde hay más ventas es cuando el precio es menor que 109.5, o sea, muchas ventas cuando la gente tiene un buen lugar para guardar los asientos y el precio del asiento está por debajo de la media. Además si el precio está entre 109.5 y 135 tenemos muchas ventas.
Las variables “shelveLoc” y “Price” arrojan bloques heterogéneos, yo diría que tienen importancia mayor que las otras variables, más adelante lo veremos.
- Usar Cross-Validation para determinar el nivel óptimo de la complejidad del árbol. Indicar si el árbol podado mejora el MSE.
El método Cross-Validation da un gráfico en el cual se puede ver el error por tamaño de árbol. Con este gráfico me apoyo para seleccionar el tamaño de árbol que minimice el MSE y ver si mejora en comparación al árbol sin podar.
set.seed(101010)
cv_tree<-cv.tree(simple_tree)
#plot(cv_tree$size, cv_tree$dev, type="b", xlab="Tamaño del árbol", ylab="MSE", col="blue") #El error sigue disminuyendo, pero voy a podar a 7 y con 11
prune_simple_tree<-prune.tree(simple_tree, best= 7)
#plot(prune_simple_tree); text(prune_simple_tree, pretty=0) #Graficando árbol podado
######## MSE Test
yhat1<-predict(prune_simple_tree, newdata=Carseats[-train,])#árbol, datos a predecir
mean((yhat1-Carseats.test)^2) #MSE 5.36-11 5.32-7 No mejora el árbol podado
errors[3,2]<-mean((yhat1-Carseats.test)^2)
######## MSE Train
yhat1_t<-predict(prune_simple_tree, newdata=Carseats[train,])#árbol, datos a predecir
mean((yhat1_t-Carseats.train)^2) #MSE 2.88-11 3.57-7
errors[3,1]<-mean((yhat1_t-Carseats.train)^2)
Fig 2. MSE calculado mediante Cross-Validation en función del tamaño del árbol.
El error disminuye conforme el tamaño del árbol crece, voy a podar con tamaño 11 porque el MSE es menor ahí y con 7 para disminuir el tamaño del árbol, una vez teniendo los árboles podados voy a usar mi conjunto de prueba para ver el rendimiento y compararlos.
| Tamaño | 7 | 11 | 15 |
|---|---|---|---|
| MSE | 5.32 | 5.36 | 5.01 |
El árbol sin podar tiene un MSE de 5.01 el cual es menor que 5.32 y que 5.36. Por lo tanto no mejoraron los árboles podados.
Fig 3. Árbol podado de tamaño 7.
Fig 4. Árbol podado de tamaño 8.
- Usar método “bagging” para analizar los datos. Reportar MSE obtenido y usar la función importance() para determinar las variables más importantes.
En este método el “out-of-bag error” estima al error de prueba, voy a calcular los dos y comparar. De igual manera voy a ver qué variables tienen más importancia.
set.seed(101010)
bagging_tree<-randomForest(Sales~ ., data=Carseats, subset=train, ntree=1000, mtry=10, importance=TRUE)
errors[4,3]<-bagging_tree$mse[1000] ##obb error 2.42
########MSE Test
yhatbg<-predict(bagging_tree, newdata=Carseats[-train,])#árbol, datos a predecir
mean((yhatbg-Carseats.test)^2) #MSE 2.68
errors[4, 2]<- mean((yhatbg-Carseats.test)^2)
########MSE Train
yhatbg_t<-predict(bagging_tree, newdata=Carseats[train,])#árbol, datos a predecir
mean((yhatbg_t-Carseats.train)^2) #MSE 0.42
errors[4, 1]<- mean((yhatbg_t-Carseats.train)^2)
########## importance
#importance(bagging_tree)
#varImpPlot(bagging_tree)| Test | 2.66 |
|---|---|
| OOB | 2.42 |
En este modelo el error de prueba es casi la mitad que los errores de prueba obtenidos con un solo árbol, es de esperarse pues estamos haciendo remuestreo sobre el conjunto de entrenamiento. Además el OOB error es menor que el error de prueba.
Ahora veré la importancia de las variables. Por el primer árbol de regresión ajustado yo diría que las variables más importantes son “ShelveLoc” y “Price”, ya que los valores en los nodos son heterogéneos.
Fig 5. Izquiera: incremento del MSE al excluir a la variable de la izquierda. Derecha: incremento de pureza al excluir a la variable de la izquierda.
Y en efecto, al excluir las variables “ShelveLoc” ó “Price” tanto el MSE y la pureza incrementan, lo que indica que estas son variables importantes.
- Usar “random forests” para analizar los datos. Reportar MSE. Usar función de importancia y describir el efecto del valor de variables utilizadas en cada división, en la tasa de error obtenida.
Primero voy a hacer “random forest”. Aquí lo interesante es ver cómo cambia el MSE al variar el valor \(m<p\), donde p es el número total de predictores; el valor default que se usa es: el número de predictores entre 3, en este caso es 3. En el cual obtuve:
set.seed(101010)
random_tree<-randomForest(Sales~ ., data=Carseats, subset=train, importance=TRUE, ntree=1000, mtry=3) #m=p/3=3
errors[5, 3]<-random_tree$mse[1000] ##obb error 2.7
########MSE Test
yhatrf<-predict(random_tree, newdata=Carseats[-train,])#árbol, datos a predecir
mean((yhatrf-Carseats.test)^2) #MSE 2.91
errors[5,2]<-mean((yhatrf-Carseats.test)^2)
########MSE Train
yhatrf_t<-predict(random_tree, newdata=Carseats[train,])#árbol, datos a predecir
mean((yhatrf_t-Carseats.train)^2) #MSE .52
errors[5,1]<-mean((yhatrf_t-Carseats.train)^2) | Test | 2.92 |
|---|---|
| OOB | 2.74 |
Luego voy a determinar el efecto de m, número de predictores, sobre la tasa de error para escoger un mejor valor.
#####Excoger el mejor valor para m<p, tal que MSE es el menor
mse<-c()
for(i in 1:10){
mse[i]<-best_tree(i)
}
x<-seq(1:10)
data<-data.frame(x, mse)
#ggplot(data, aes(x=x, y=mse))+
# geom_line(linetype="longdash", size=0, color="blue")+
# geom_point(shape=15, color="purple", size=5)+
# theme_bw()+
# theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
# xlab("Valor de m")+
# ylab("MSE")
Fig 6. Gráfica del MSE según el número de predictores usados para hacer el random forest.
El valor propuesto es 7. Voy a ajustar otro árbol con este número de predictores y calcular su tasa de error y su OBB.
################ si m=7, el MSE es el menor
set.seed(101010)
random_tree_best<-randomForest(Sales~ ., data=Carseats, subset=train, importance=TRUE, ntree=1000, mtry=7) #m=7
errors[6, 3]<-random_tree_best$mse[1000] ##obb error 2.39
########MSE Test
yhatrfb<-predict(random_tree_best, newdata=Carseats[-train,])#árbol, datos a predecir
mean((yhatrfb-Carseats.test)^2) #MSE 2.64
errors[6,2]<-mean((yhatrfb-Carseats.test)^2)
########MSE Train
yhatrfb_t<-predict(random_tree_best, newdata=Carseats[train,])#árbol, datos a predecir
mean((yhatrfb_t-Carseats.train)^2) #MSE .42
errors[6,1]<-mean((yhatrfb_t-Carseats.train)^2)| Tamaño | m=7 | m=3 |
|---|---|---|
| OOB | 2.39 | 2.74 |
| Test | 2.64 | 2.92 |
En efecto, usando m=7 hay una tasa de error y OBB menor que usando m=3. Ahora la importancia, como el modelo con m=7 tiene tasas de error menor voy a ver la importancia sobre este.
Fig 7. Gráfica del MSE según el número de predictores usados para hacer el random forest.
De nuevo, las variables con más importancia son: “ShelveLoc” y “Price”.
Por último, una tabla con las tasas de errores calculadas. En la cual, random forest con m=7 es la que tiene menor tasa de error de prueba y menor OOB.
| Train Errors | Test Errors | Out-of-Bag Errors | |
|---|---|---|---|
| Single Tree | 2.39 | 5.01 | - |
| Pruned Tree 11 | 2.88 | 5.36 | - |
| Pruned Tree 7 | 3.57 | 5.32 | - |
| Bagged Tree | 0.42 | 2.68 | 2.42 |
| RandomForestDef | 0.56 | 2.92 | 2.74 |
| RandomForestBes | 0.42 | 2.64 | 2.39 |
Ejercicio 7.
- Crear plot mostrando las tasas de error de prueba resultantes de hacer randomforest para valores de “mtry” y “ntree”.
Voy a usar los valores 4, 5, 6, 7, 8 y 13. Primero divido en train y test, luego voy a ajustar los randomforests a cada uno de los valores de m, haciendo variar el número de árboles.
#Test y train
set.seed(010101)
train<-sample(1:506, 506/2) #El lirbo entrena con una mitad y prueba con la otra.
Boston.test<-Boston[-train, "medv"]
##############funcion
MSErf<-function(m, n){
set.seed(1000)
rf<-randomForest(medv~ ., data=Boston, subset=train, importance=TRUE, ntree=n, mtry=m)
yhat<-predict(rf, newdata=Boston[-train,]) #árbol, datos a predecir
mse<-mean((yhat-Boston.test)^2)
print(mse)
}
######Graficando
mse13<-c()
for(i in 1:500){
mse13[i]<-MSErf(13, i)
}
mse7<-c()
for(i in 1:500){
mse7[i]<-MSErf(7, i)
}
mse4<-c()
for(i in 1:500){
mse4[i]<-MSErf(4, i)
}
mse10<-c()
for(i in 1:500){
mse10[i]<-MSErf(10, i)
}
mse6<-c()
for(i in 1:500){
mse6[i]<-MSErf(6, i)
}
#plot(x, mse7, type="l", col="red", xlab="ntrees", ylab="MSE")
#lines(x, mse13, type="l", col="purple")
#lines(x, mse4, type="l", col="yellow")
#lines(x, mse10, type="l", col="black")
#lines(x, mse6, type="l", col="green")
#legend("topright", legend = c("mse13-mor", "mse10-negr", "mse4-amar", "mse7-rojo", "mse6-verd"))
Fig 8. Gráfica del MSE según el número de predictores elegidos en función del numero de árboles.
Al hacer bagging obtenemos el MSE más alto (m=13). Cuando m toma los valores 5, 6 y 7 el MSE se parece mucho y es menor que en los demás casos, pues cuando m=4 ó m=8 el MSE es mayor que el de m=5, 6 y 7.