Лабораторная_8
Математическое моделирование
Практика 8
Модели на основе деревьев
В практических примерах ниже показано как:
- строить регрессионные деревья;
- строить деревья классификации;
- делать обрезку дерева;
- использовать бэггинг, бустинг, случайный лес для улучшения качества прогнозирования.
Модели: деревья решений.
Данные: Sales {ISLR}, Boston {MASS}
# Загрузка пакетов
library('tree') # деревья tree()## Warning: package 'tree' was built under R version 3.5.3library('ISLR') # набор данных Carseats## Warning: package 'ISLR' was built under R version 3.5.3library('GGally') # матричный график разброса ggpairs()## Warning: package 'GGally' was built under R version 3.5.3## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 3.5.3library('MASS') # набор данных Boston## Warning: package 'MASS' was built under R version 3.5.3library('randomForest') # случайный лес randomForest()## Warning: package 'randomForest' was built under R version 3.5.3## randomForest 4.6-14## Type rfNews() to see new features/changes/bug fixes.##
## Attaching package: 'randomForest'## The following object is masked from 'package:ggplot2':
##
## marginlibrary('gbm') # бустинг gbm()## Warning: package 'gbm' was built under R version 3.5.3## Loaded gbm 2.1.5data(Boston) # открываем данные
high.medv <- ifelse(Boston$medv <= 25, "0", "1")
# присоединяем к таблице данных
Boston <- cbind(Boston, high.medv)
# ядро генератора случайных чисел
my.seed <- 3
set.seed(my.seed)
train <- sample(1:nrow(Boston), nrow(Boston)/2) # обучающая выборка -- 50%p <- ggpairs(Boston[, c(14, 1:4)])
suppressMessages(print(p))
p <- ggpairs(Boston[, c(14, 5:8)])
suppressMessages(print(p))
p <- ggpairs(Boston[, c(14, 9:13)])
suppressMessages(print(p))
tree.boston <- tree(high.medv ~ ., Boston, subset = train)
summary(tree.boston)##
## Classification tree:
## tree(formula = high.medv ~ ., data = Boston, subset = train)
## Variables actually used in tree construction:
## [1] "medv"
## Number of terminal nodes: 2
## Residual mean deviance: 0 = 0 / 251
## Misclassification error rate: 0 = 0 / 253# визуализация
plot(tree.boston)
text(tree.boston, pretty = 0)
Снова сделаем обрезку дерева в целях улучшения качества прогноза.
# обрезка дерева
cv.boston <- cv.tree(tree.boston)# размер дерева с минимальной ошибкой
plot(cv.boston$size, cv.boston$dev, type = 'b')
opt.size <- cv.boston$size[cv.boston$dev == min(cv.boston$dev)]
abline(v = opt.size, col = 'red', 'lwd' = 2) # соотв. вертикальная прямая
mtext(opt.size, at = opt.size, side = 1, col = 'red', line = 1)
В данном случаем минимум ошибки соответствует дереву с 2 узлами. Сделаем прогноз и рассчитаем MSE.
# обучающая выборка
yhat <- predict(tree.boston, newdata=Boston[-train, ])
# прогноз по лучшей модели (2 узлами)
boston.test <- Boston[-train, "high.medv"]
mse.test <- mean((yhat - train)^2)
mse.test## [1] 84937.05Метод случайного леса
Рассмотрим более сложные методы улучшения качества дерева. Бэггинг – частный случай случайного леса с \(m = p\), поэтому и то, и другое можно построить функцией randomForest().
Для начала используем бэггинг, причём возьмём все 14 предикторов на каждом шаге (аргумент mtry).
# бэггинг с 14 предикторами
boston.test <- Boston[-train,]
set.seed(my.seed)
bag.Boston <- randomForest(high.medv ~ ., data = Boston, subset = train,
mtry = 14, importance = TRUE)
bag.Boston##
## Call:
## randomForest(formula = high.medv ~ ., data = Boston, mtry = 14, importance = TRUE, subset = train)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 14
##
## OOB estimate of error rate: 0.4%
## Confusion matrix:
## 0 1 class.error
## 0 185 0 0.00000000
## 1 1 67 0.01470588# прогноз
yhat.bag <- predict(bag.Boston, newdata = Boston[-train, ])
boston.test <- Boston[-train, "high.medv"]
tbl <- table(yhat.bag, boston.test)
tbl## boston.test
## yhat.bag 0 1
## 0 197 0
## 1 0 56acc.test <- sum(diag(tbl))/sum(tbl)
names(acc.test)[length(acc.test)] <- 'Boston.test'
acc.test## Boston.test
## 1Можно изменить число деревьев с помощью аргумента ntree.
# бэггинг с 13 предикторами и 25 деревьями
bag.Boston <- randomForest(high.medv ~ ., data = Boston, subset = train,
mtry = 13, ntree = 25)
# прогноз
yhat.bag <- predict(bag.Boston, newdata = Boston[-train, ])
boston.test <- Boston[-train, "high.medv"]
tbl <- table(yhat.bag, boston.test)
tbl## boston.test
## yhat.bag 0 1
## 0 197 0
## 1 0 56acc.test <- sum(diag(tbl))/sum(tbl)
names(acc.test)[length(acc.test)] <- 'Boston.test'
acc.test## Boston.test
## 1Критерий в Aсс в обоих случаях равен 1, что говорит о 100%-й точности.
Теперь попробуем вырастить случайный лес. Берём 6 предикторов на каждом шаге.
# обучаем модель
set.seed(my.seed)
rf.boston <- randomForest(high.medv ~ ., data = Boston, subset = train,
mtry = 6, importance = TRUE)
# важность предикторов
importance(rf.boston) # оценки ## 0 1 MeanDecreaseAccuracy MeanDecreaseGini
## crim 3.472723 2.31465973 3.736177 0.48989622
## zn 2.153950 -0.08939633 2.167918 0.18126663
## indus 5.436363 5.10376874 7.331979 3.66296647
## chas 1.297196 0.02667855 1.075199 0.02696158
## nox 1.917966 3.88638408 3.669611 0.69257247
## rm 10.394940 7.51125150 11.929978 21.18831485
## age 3.147285 -0.61109669 2.461373 0.36541658
## dis 3.547443 -1.27326121 2.748210 0.78529972
## rad 2.847373 -0.62961661 2.353177 0.13760783
## tax 6.067422 3.29743435 6.294220 1.80899831
## ptratio 4.153972 3.77337088 5.030666 1.58924664
## black 1.567298 1.95279783 2.462298 0.24947582
## lstat 5.015059 3.05782008 5.868095 5.97438276
## medv 68.683324 64.87955586 78.745276 62.37088266varImpPlot(rf.boston) # графики
По полученным данным можно сделать вывод о том, что наибольшее влияние в модели оказывают такие показатели, как medv (средняя стоимость домов, занимаемых владельцами, в 1000 долларов) и rm (среднее количество комнат в доме).