Математическое моделирование. Практика 8
Трапезникова Анастасия
28 апреля, 2021
Задание:
Построить две модели для прогноза на основе дерева решений:
1. для непрерывной зависимой переменной;
2. для категориальной зависимой переменной.
Данные и переменные указаны в таблице с вариантами.
Ядро генератора случайных чисел – номер варианта.
Для каждой модели:
1. Указать настроечные параметры метода из своего варианта (например: количество узлов, количество предикторов, скорость обучения).
2. Подогнать модель на обучающей выборке (50% наблюдений). Рассчитать MSE на тестовой выборке.
3. Перестроить модель с помощью метода, указанного в варианте.
4. Сделать прогноз по модели с подобранными в п.3 параметрами на тестовой выборке, оценить его точность и построить график «прогноз-реализация».
Вариант 1
| Номер варианта | Данные | Непрерывный Y | Категориальный Y | Объясняющие переменные | Метод подгонки моделей |
|---|---|---|---|---|---|
| 1 | Boston {MASS} |
\(medv\) | \(high.medv = \begin{cases} \begin{array}{lcl} 1, & если & medv >= 25 \\ 0, & если & medv < 25 \end{array} \end{cases}\) | все остальные | дерево с обрезкой ветвей |
Как сдавать: прислать на почту преподавателя ссылки:
* на html-отчёт с видимыми блоками кода (блоки кода с параметром echo = T), размещённый на rpubs.com.
* на код, генерирующий отчёт, в репозитории на github.com. В текст отчёта включить постановку задачи и ответы на вопросы задания.
Решение:
Деревья решений
Подключаем набор данных Boston
library('GGally') # матричный график разброса ggpairs()
library('tree') # деревья tree()
library(MASS) # набор данных Boston
data(Boston)
head(Boston)## crim zn indus chas nox rm age dis rad tax ptratio black lstat
## 1 0.00632 18 2.31 0 0.538 6.575 65.2 4.0900 1 296 15.3 396.90 4.98
## 2 0.02731 0 7.07 0 0.469 6.421 78.9 4.9671 2 242 17.8 396.90 9.14
## 3 0.02729 0 7.07 0 0.469 7.185 61.1 4.9671 2 242 17.8 392.83 4.03
## 4 0.03237 0 2.18 0 0.458 6.998 45.8 6.0622 3 222 18.7 394.63 2.94
## 5 0.06905 0 2.18 0 0.458 7.147 54.2 6.0622 3 222 18.7 396.90 5.33
## 6 0.02985 0 2.18 0 0.458 6.430 58.7 6.0622 3 222 18.7 394.12 5.21
## medv
## 1 24.0
## 2 21.6
## 3 34.7
## 4 33.4
## 5 36.2
## 6 28.7str(Boston)## 'data.frame': 506 obs. of 14 variables:
## $ crim : num 0.00632 0.02731 0.02729 0.03237 0.06905 ...
## $ zn : num 18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
## $ indus : num 2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
## $ chas : int 0 0 0 0 0 0 0 0 0 0 ...
## $ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
## $ rm : num 6.58 6.42 7.18 7 7.15 ...
## $ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
## $ dis : num 4.09 4.97 4.97 6.06 6.06 ...
## $ rad : int 1 2 2 3 3 3 5 5 5 5 ...
## $ tax : num 296 242 242 222 222 222 311 311 311 311 ...
## $ ptratio: num 15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
## $ black : num 397 397 393 395 397 ...
## $ lstat : num 4.98 9.14 4.03 2.94 5.33 ...
## $ medv : num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...# ?BostonЗагрузим таблицу с данными по стоимости жилья в пригороде Бостона и добавим к ней переменную high.medv – “высокая медианная стоимость” со значениями:
1, если medv >= 25
0, если medv < 25
# новая переменная
high.medv <- ifelse(Boston$medv >= 25, 1, 0)
high.medv <- factor(high.medv, labels = c('yes', 'no'))
Boston$high.medv <- high.medv
# матричные графики разброса переменных
p <- ggpairs(Boston[, c(15, 1:4)], aes(color = high.medv))
suppressMessages(print(p))
p <- ggpairs(Boston[, c(15, 5:8)], aes(color = high.medv))
suppressMessages(print(p))
p <- ggpairs(Boston[, c(15, 9:14)], aes(color = high.medv))
suppressMessages(print(p))
# модель бинарного дерева без переменной medv
tree.boston <- tree(high.medv ~ (.-medv), Boston)
summary(tree.boston)##
## Classification tree:
## tree(formula = high.medv ~ (. - medv), data = Boston)
## Variables actually used in tree construction:
## [1] "rm" "lstat" "crim" "tax" "age" "indus" "nox"
## Number of terminal nodes: 16
## Residual mean deviance: 0.1975 = 96.77 / 490
## Misclassification error rate: 0.03557 = 18 / 506# график результата:
# ветви
plot(tree.boston)
# добавим подписи
text(tree.boston, pretty = 0)
# посмотреть всё дерево в консоли
tree.boston ## node), split, n, deviance, yval, (yprob)
## * denotes terminal node
##
## 1) root 506 580.900 yes ( 0.739130 0.260870 )
## 2) rm < 6.5455 362 165.900 yes ( 0.939227 0.060773 )
## 4) lstat < 9.845 93 91.390 yes ( 0.806452 0.193548 )
## 8) crim < 2.17451 87 69.810 yes ( 0.862069 0.137931 )
## 16) rm < 6.1355 33 0.000 yes ( 1.000000 0.000000 ) *
## 17) rm > 6.1355 54 57.210 yes ( 0.777778 0.222222 )
## 34) tax < 279 18 24.950 no ( 0.500000 0.500000 )
## 68) age < 35.2 7 5.742 no ( 0.142857 0.857143 ) *
## 69) age > 35.2 11 12.890 yes ( 0.727273 0.272727 ) *
## 35) tax > 279 36 20.650 yes ( 0.916667 0.083333 )
## 70) indus < 5.975 23 0.000 yes ( 1.000000 0.000000 ) *
## 71) indus > 5.975 13 14.050 yes ( 0.769231 0.230769 ) *
## 9) crim > 2.17451 6 0.000 no ( 0.000000 1.000000 ) *
## 5) lstat > 9.845 269 41.610 yes ( 0.985130 0.014870 )
## 10) indus < 3.305 5 6.730 yes ( 0.600000 0.400000 ) *
## 11) indus > 3.305 264 23.520 yes ( 0.992424 0.007576 ) *
## 3) rm > 6.5455 144 157.400 no ( 0.236111 0.763889 )
## 6) nox < 0.659 131 115.300 no ( 0.160305 0.839695 )
## 12) rm < 7.0105 71 86.230 no ( 0.295775 0.704225 )
## 24) tax < 267.5 20 0.000 no ( 0.000000 1.000000 ) *
## 25) tax > 267.5 51 69.100 no ( 0.411765 0.588235 )
## 50) age < 35.2 15 11.780 no ( 0.133333 0.866667 ) *
## 51) age > 35.2 36 49.800 yes ( 0.527778 0.472222 )
## 102) nox < 0.435 6 0.000 yes ( 1.000000 0.000000 ) *
## 103) nox > 0.435 30 41.050 no ( 0.433333 0.566667 )
## 206) nox < 0.526 14 11.480 no ( 0.142857 0.857143 ) *
## 207) nox > 0.526 16 19.870 yes ( 0.687500 0.312500 )
## 414) indus < 12.91 8 0.000 yes ( 1.000000 0.000000 ) *
## 415) indus > 12.91 8 10.590 no ( 0.375000 0.625000 ) *
## 13) rm > 7.0105 60 0.000 no ( 0.000000 1.000000 ) *
## 7) nox > 0.659 13 0.000 yes ( 1.000000 0.000000 ) *Теперь построим дерево на обучающей выборке и оценим ошибку на тестовой.
# ядро генератора случайных чисел по номеру варианта
my.seed <- 1
set.seed(my.seed)
# обучающая выборка 50%
train <- sample(1:nrow(Boston), nrow(Boston)*0.5)
# тестовая выборка
boston.test <- Boston[-train,]
high.medv.test <- high.medv[-train]
# строим дерево на обучающей выборке
tree.boston <- tree(high.medv ~ (.-medv), Boston, subset = train)
summary(tree.boston)##
## Classification tree:
## tree(formula = high.medv ~ (. - medv), data = Boston, subset = train)
## Variables actually used in tree construction:
## [1] "rm" "lstat" "nox" "age" "crim" "tax"
## Number of terminal nodes: 12
## Residual mean deviance: 0.1379 = 33.23 / 241
## Misclassification error rate: 0.02767 = 7 / 253# делаем прогноз на тестовой
tree.pred <- predict(tree.boston, boston.test, type = "class")
# матрица неточностей
tbl <- table(tree.pred, high.medv.test)
tbl## high.medv.test
## tree.pred yes no
## yes 159 14
## no 19 61# ACC на тестовой
acc.test <- sum(diag(tbl))/sum(tbl)
names(acc.test)[length(acc.test)] <- 'Boston.class.tree.all'
acc.test## Boston.class.tree.all
## 0.8695652Обобщённая характеристика точности: доля верных прогнозов: 0,869
Теперь обрежем дерево, используя в качестве критерия частоту ошибок классификации. Функция cv.tree() проводит кросс-валидацию для выбора лучшего дерева, аргумент prune.misclass означает, что мы минимизируем ошибку классификации.
set.seed(my.seed)
cv.boston <- cv.tree(tree.boston, FUN = prune.misclass)
# имена элементов полученного объекта
names(cv.boston)## [1] "size" "dev" "k" "method"# сам объект
cv.boston## $size
## [1] 12 9 5 3 2 1
##
## $dev
## [1] 19 15 13 17 23 58
##
## $k
## [1] -Inf 0.00 0.75 1.50 7.00 37.00
##
## $method
## [1] "misclass"
##
## attr(,"class")
## [1] "prune" "tree.sequence"Графики изменения параметров метода по ходу обрезки дерева
# 1. ошибка с кросс-валидацией в зависимости от числа узлов
par(mfrow = c(1, 2))
plot(cv.boston$size, cv.boston$dev, type = "b",
ylab = 'Частота ошибок с кросс-вал. (dev)',
xlab = 'Число узлов (size)')
# размер дерева с минимальной ошибкой
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. ошибка с кросс-валидацией в зависимости от штрафа на сложность
plot(cv.boston$k, cv.boston$dev, type = "b",
ylab = 'Частота ошибок с кросс-вал. (dev)',
xlab = 'Штраф за сложность (k)')
На графике слева видно, что минимум частоты ошибок достигается при числе узлов 5.
Оценим точность дерева с 5 узлами.
# дерево с 5 узлами
prune.boston <- prune.misclass(tree.boston, best = 5)
# визуализация
plot(prune.boston)
text(prune.boston, pretty = 0)
# прогноз на тестовую выборку
tree.pred <- predict(prune.boston, boston.test, type = "class")
# матрица неточностей
tbl <- table(tree.pred, high.medv.test)
tbl## high.medv.test
## tree.pred yes no
## yes 155 13
## no 23 62# ACC на тестовой
acc.test <- c(acc.test, sum(diag(tbl))/sum(tbl))
names(acc.test)[length(acc.test)] <- 'Boston.class.tree.5'
acc.test## Boston.class.tree.all Boston.class.tree.5
## 0.8695652 0.8577075Точность этой модели незначительно снизилась и составляет 0,857.
Увеличив количество узлов, получим более глубокое дерево и немого более точное:
# дерево с 12 узлами
prune.boston <- prune.misclass(tree.boston, best = 12)
# визуализация
plot(prune.boston)
text(prune.boston, pretty = 0)
# прогноз на тестовую выборку
tree.pred <- predict(prune.boston, boston.test, type = "class")
# матрица неточностей
tbl <- table(tree.pred, high.medv.test)
tbl## high.medv.test
## tree.pred yes no
## yes 159 14
## no 19 61# ACC на тестовой
acc.test <- c(acc.test, sum(diag(tbl))/sum(tbl))
names(acc.test)[length(acc.test)] <- 'Carseats.class.tree.12'
acc.test## Boston.class.tree.all Boston.class.tree.5 Carseats.class.tree.12
## 0.8695652 0.8577075 0.8695652# сбрасываем графические параметры
par(mfrow = c(1, 1))Регрессионные деревья
Продолжим использовать набор данных Boston, загрузив его заново
my.seed <- 1
data("Boston")
# ?Boston
str(Boston)## 'data.frame': 506 obs. of 14 variables:
## $ crim : num 0.00632 0.02731 0.02729 0.03237 0.06905 ...
## $ zn : num 18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
## $ indus : num 2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
## $ chas : int 0 0 0 0 0 0 0 0 0 0 ...
## $ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
## $ rm : num 6.58 6.42 7.18 7 7.15 ...
## $ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
## $ dis : num 4.09 4.97 4.97 6.06 6.06 ...
## $ rad : int 1 2 2 3 3 3 5 5 5 5 ...
## $ tax : num 296 242 242 222 222 222 311 311 311 311 ...
## $ ptratio: num 15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
## $ black : num 397 397 393 395 397 ...
## $ lstat : num 4.98 9.14 4.03 2.94 5.33 ...
## $ medv : num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...head(Boston)## crim zn indus chas nox rm age dis rad tax ptratio black lstat
## 1 0.00632 18 2.31 0 0.538 6.575 65.2 4.0900 1 296 15.3 396.90 4.98
## 2 0.02731 0 7.07 0 0.469 6.421 78.9 4.9671 2 242 17.8 396.90 9.14
## 3 0.02729 0 7.07 0 0.469 7.185 61.1 4.9671 2 242 17.8 392.83 4.03
## 4 0.03237 0 2.18 0 0.458 6.998 45.8 6.0622 3 222 18.7 394.63 2.94
## 5 0.06905 0 2.18 0 0.458 7.147 54.2 6.0622 3 222 18.7 396.90 5.33
## 6 0.02985 0 2.18 0 0.458 6.430 58.7 6.0622 3 222 18.7 394.12 5.21
## medv
## 1 24.0
## 2 21.6
## 3 34.7
## 4 33.4
## 5 36.2
## 6 28.7# матричные графики разброса переменных
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))
# обучающая выборка
set.seed(my.seed)
train <- sample(1:nrow(Boston), nrow(Boston)/2) # обучающая выборка -- 50%Построим дерево регрессии для зависимой переменной medv: медианная стоимости домов, в которых живут собственники (тыс. долл.).
# обучаем модель
tree.boston <- tree(medv ~ ., Boston, subset = train)
summary(tree.boston)##
## Regression tree:
## tree(formula = medv ~ ., data = Boston, subset = train)
## Variables actually used in tree construction:
## [1] "rm" "lstat" "crim" "age"
## Number of terminal nodes: 7
## Residual mean deviance: 10.38 = 2555 / 246
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -10.1800 -1.7770 -0.1775 0.0000 1.9230 16.5800# визуализация
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)
В данном случаем минимум ошибки соответствует самому сложному дереву, с 7 узлами. Покажем, как при желании можно обрезать дерево до 6 узлов (ошибка ненамного выше, чем минимальная).
# дерево с 6 узлами
prune.boston = prune.tree(tree.boston, best = 6)
# визуализация
plot(prune.boston)
text(prune.boston, pretty = 0)
Прогноз сначала сделаем по необрезанному дереву, т.к. там ошибка, оцененная по методу перекрёстной проверки, минимальна.
# прогноз по лучшей модели (7 узлов)
yhat <- predict(tree.boston, newdata = Boston[-train, ])
boston.test <- Boston[-train, "medv"]
# график "прогноз -- реализация"
plot(yhat, boston.test)
# линия идеального прогноза
abline(0, 1)
# MSE на тестовой выборке
mse.test <- mean((yhat - boston.test)^2)
names(mse.test)[length(mse.test)] <- 'Boston.regr.tree.7'
mse.test## Boston.regr.tree.7
## 35.28688MSE на тестовой выборке равна 35.28.
Теперь сделаем прогноз по обрезанному дереву.
# прогноз по лучшей модели (6 узлов)
yhat.prune <- predict(prune.boston, newdata = Boston[-train, ])
boston.test <- Boston[-train, "medv"]
# график "прогноз -- реализация"
plot(yhat.prune, boston.test)
# линия идеального прогноза
abline(0, 1)
# MSE обрезанного дерева на тестовой выборке
mse.prune.test <- mean((yhat.prune - boston.test)^2)
names(mse.prune.test)[length(mse.prune.test)] <- 'Boston.regr.tree.6'
mse.prune.test## Boston.regr.tree.6
## 35.16439MSE обрезанного дерева на тестовой выборке равна 35.16
Источники
- Джеймс Г., Уиттон Д., Хасти Т., Тибширани Р. Введение в статистическое обучение с примерами на языке R / пер. с англ. С.Э. Мастицкого. – М.: ДМК Пресс, 2016 – 450 с. Репозиторий с примерами к книге на русском языке: https://github.com/ranalytics/islr-ru