Математическое моделирование. Практика 8
Трапезникова Анастасия
28 апреля, 2021
Задание:
Построить две модели для прогноза на основе дерева решений:
1. для непрерывной зависимой переменной;
2. для категориальной зависимой переменной.
Данные и переменные указаны в таблице с вариантами.
Ядро генератора случайных чисел – номер варианта.
Для каждой модели:
1. Указать настроечные параметры метода из своего варианта (например: количество узлов, количество предикторов, скорость обучения).
2. Подогнать модель на обучающей выборке (50% наблюдений). Рассчитать MSE на тестовой выборке.
3. Перестроить модель с помощью метода, указанного в варианте.
4. Сделать прогноз по модели с подобранными в п.3 параметрами на тестовой выборке, оценить его точность и построить график «прогноз-реализация».
Вариант 9
| Номер варианта | Данные | Непрерывный Y | Категориальный Y | Объясняющие переменные | Метод подгонки моделей |
|---|---|---|---|---|---|
| 9 | Auto {ISLR} |
mpg |
\(high.mpg = \begin{cases} \begin{array}{lcl} 1, & если & mpg >= 29 \\ 0, & если & mpg < 29 \end{array} \end{cases}\) | остальные, кроме name |
дерево с обрезкой ветвей |
Как сдавать: прислать на почту преподавателя ссылки:
* на html-отчёт с видимыми блоками кода (блоки кода с параметром echo = T), размещённый на rpubs.com.
* на код, генерирующий отчёт, в репозитории на github.com. В текст отчёта включить постановку задачи и ответы на вопросы задания.
Решение:
Деревья решений
Подключаем набор данных Auto
library('ISLR') # набор данных Auto
library('GGally') # матричный график разброса ggpairs()
library('tree') # деревья tree()
data(Auto)Загрузим таблицу с данными по параметрам автомобилей и добавим к ней переменную high.mpg – “высокий расход топлива” со значениями:
1, если mpg >= 29
0, если mpg < 29
str(Auto)## 'data.frame': 392 obs. of 9 variables:
## $ mpg : num 18 15 18 16 17 15 14 14 14 15 ...
## $ cylinders : num 8 8 8 8 8 8 8 8 8 8 ...
## $ displacement: num 307 350 318 304 302 429 454 440 455 390 ...
## $ horsepower : num 130 165 150 150 140 198 220 215 225 190 ...
## $ weight : num 3504 3693 3436 3433 3449 ...
## $ acceleration: num 12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
## $ year : num 70 70 70 70 70 70 70 70 70 70 ...
## $ origin : num 1 1 1 1 1 1 1 1 1 1 ...
## $ name : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...head(Auto)## mpg cylinders displacement horsepower weight acceleration year origin
## 1 18 8 307 130 3504 12.0 70 1
## 2 15 8 350 165 3693 11.5 70 1
## 3 18 8 318 150 3436 11.0 70 1
## 4 16 8 304 150 3433 12.0 70 1
## 5 17 8 302 140 3449 10.5 70 1
## 6 15 8 429 198 4341 10.0 70 1
## name
## 1 chevrolet chevelle malibu
## 2 buick skylark 320
## 3 plymouth satellite
## 4 amc rebel sst
## 5 ford torino
## 6 ford galaxie 500# новая переменная
high.mpg <- ifelse(Auto$mpg >= 29, 1, 0)
high.mpg <- factor(high.mpg, labels = c('yes', 'no'))
Auto$high.mpg <- high.mpg
# матричные графики разброса переменных
p <- ggpairs(Auto[, c(10, 1:5)], aes(color = high.mpg))
suppressMessages(print(p))
p <- ggpairs(Auto[, c(10, 6:8)], aes(color = high.mpg))
suppressMessages(print(p))
# модель бинарного дерева без переменных mpg и name
tree.auto <- tree(high.mpg ~ (.-name-mpg), Auto)
summary(tree.auto)##
## Classification tree:
## tree(formula = high.mpg ~ (. - name - mpg), data = Auto)
## Variables actually used in tree construction:
## [1] "horsepower" "year" "weight" "displacement"
## Number of terminal nodes: 10
## Residual mean deviance: 0.2449 = 93.56 / 382
## Misclassification error rate: 0.05102 = 20 / 392# график результата:
# ветви
plot(tree.auto)
# добавим подписи
text(tree.auto, pretty = 0)
# посмотреть всё дерево в консоли
tree.auto ## node), split, n, deviance, yval, (yprob)
## * denotes terminal node
##
## 1) root 392 451.500 yes ( 0.737245 0.262755 )
## 2) horsepower < 84.5 132 161.900 no ( 0.303030 0.696970 )
## 4) year < 76.5 51 67.350 yes ( 0.627451 0.372549 )
## 8) weight < 2091 22 23.580 no ( 0.227273 0.772727 )
## 16) displacement < 90.5 15 7.348 no ( 0.066667 0.933333 ) *
## 17) displacement > 90.5 7 9.561 yes ( 0.571429 0.428571 ) *
## 9) weight > 2091 29 14.560 yes ( 0.931034 0.068966 ) *
## 5) year > 76.5 81 52.220 no ( 0.098765 0.901235 )
## 10) weight < 2255 58 0.000 no ( 0.000000 1.000000 ) *
## 11) weight > 2255 23 29.720 no ( 0.347826 0.652174 ) *
## 3) horsepower > 84.5 260 91.110 yes ( 0.957692 0.042308 )
## 6) year < 79.5 227 12.850 yes ( 0.995595 0.004405 )
## 12) year < 78.5 208 0.000 yes ( 1.000000 0.000000 ) *
## 13) year > 78.5 19 7.835 yes ( 0.947368 0.052632 ) *
## 7) year > 79.5 33 40.490 yes ( 0.696970 0.303030 )
## 14) displacement < 137 10 12.220 no ( 0.300000 0.700000 )
## 28) displacement < 115.5 5 6.730 yes ( 0.600000 0.400000 ) *
## 29) displacement > 115.5 5 0.000 no ( 0.000000 1.000000 ) *
## 15) displacement > 137 23 17.810 yes ( 0.869565 0.130435 ) *Теперь построим дерево на обучающей выборке и оценим ошибку на тестовой.
# ядро генератора случайных чисел по номеру варианта
my.seed <- 9
set.seed(my.seed)
# обучающая выборка 50%
train <- sample(1:nrow(Auto), 200) #nrow(Auto)*0.5 - даёт слишком мало узлов
# тестовая выборка
auto.test <- Auto[-train,]
high.mpg.test <- high.mpg[-train]
# строим дерево на обучающей выборке
tree.auto <- tree(high.mpg ~ (.-name-mpg), Auto, subset = train)
summary(tree.auto)##
## Classification tree:
## tree(formula = high.mpg ~ (. - name - mpg), data = Auto, subset = train)
## Variables actually used in tree construction:
## [1] "weight" "year" "displacement" "horsepower"
## Number of terminal nodes: 8
## Residual mean deviance: 0.1513 = 29.06 / 192
## Misclassification error rate: 0.03 = 6 / 200# делаем прогноз
tree.pred <- predict(tree.auto, auto.test, type = "class")
# матрица неточностей
tbl <- table(tree.pred, high.mpg.test)
tbl## high.mpg.test
## tree.pred yes no
## yes 132 7
## no 8 45# ACC на тестовой
acc.test <- sum(diag(tbl))/sum(tbl)
names(acc.test)[length(acc.test)] <- 'Auto.class.tree.all'
acc.test## Auto.class.tree.all
## 0.921875Обобщённая характеристика точности: доля верных прогнозов: 0,921
Теперь обрезаем дерево, используя в качестве критерия частоту ошибок классификации. Функция cv.tree() проводит кросс-валидацию для выбора лучшего дерева, аргумент prune.misclass означает, что мы минимизируем ошибку классификации.
set.seed(my.seed)
cv.auto <- cv.tree(tree.auto, FUN = prune.misclass)
# имена элементов полученного объекта
names(cv.auto)## [1] "size" "dev" "k" "method"# сам объект
cv.auto## $size
## [1] 8 6 2 1
##
## $dev
## [1] 24 23 25 51
##
## $k
## [1] -Inf 0.0 2.5 35.0
##
## $method
## [1] "misclass"
##
## attr(,"class")
## [1] "prune" "tree.sequence"Графики изменения параметров метода по ходу обрезки дерева
# 1. ошибка с кросс-валидацией в зависимости от числа узлов
par(mfrow = c(1, 2))
plot(cv.auto$size, cv.auto$dev, type = "b",
ylab = 'Частота ошибок с кросс-вал. (dev)',
xlab = 'Число узлов (size)')
# размер дерева с минимальной ошибкой
opt.size <- cv.auto$size[cv.auto$dev == min(cv.auto$dev)]
abline(v = opt.size, col = 'red', 'lwd' = 2) # соотв. вертикальная прямая
mtext(opt.size, at = opt.size, side = 1, col = 'red', line = 1)
# 2. ошибка с кросс-валидацией в зависимости от штрафа на сложность
plot(cv.auto$k, cv.auto$dev, type = "b",
ylab = 'Частота ошибок с кросс-вал. (dev)',
xlab = 'Штраф за сложность (k)')
Как видно на графике слева, минимум частоты ошибок достигается при числе узлов 6.
Оценим точность дерева с 6 узлами.
# дерево с 6 узлами
prune.auto <- prune.misclass(tree.auto, best = 6)
# визуализация
plot(prune.auto)
text(prune.auto, pretty = 0)
# прогноз на тестовую выборку
tree.pred <- predict(prune.auto, auto.test, type = "class")
# матрица неточностей
tbl <- table(tree.pred, high.mpg.test)
tbl## high.mpg.test
## tree.pred yes no
## yes 132 7
## no 8 45# ACC на тестовой
acc.test <- c(acc.test, sum(diag(tbl))/sum(tbl))
names(acc.test)[length(acc.test)] <- 'Auto.class.tree.6'
acc.test## Auto.class.tree.all Auto.class.tree.6
## 0.921875 0.921875Точность этой модели не изменилась и составляет 0,921.
Увеличив количество узлов, получим точно такое же дерево:
# дерево с 7 узлами
prune.auto <- prune.misclass(tree.auto, best = 7)
# визуализация
plot(prune.auto)
text(prune.auto, pretty = 0)
# прогноз на тестовую выборку
tree.pred <- predict(prune.auto, auto.test, type = "class")
# матрица неточностей
tbl <- table(tree.pred, high.mpg.test)
tbl## high.mpg.test
## tree.pred yes no
## yes 132 7
## no 8 45# ACC на тестовой
acc.test <- c(acc.test, sum(diag(tbl))/sum(tbl))
names(acc.test)[length(acc.test)] <- 'Carseats.class.tree.7'
acc.test## Auto.class.tree.all Auto.class.tree.6 Carseats.class.tree.7
## 0.921875 0.921875 0.921875# сбрасываем графические параметры
par(mfrow = c(1, 1))Регрессионные деревья
Продолжим использовать набор данных Auto, загрузив его заново
data("Auto")
# ?Auto
str(Auto)## 'data.frame': 392 obs. of 9 variables:
## $ mpg : num 18 15 18 16 17 15 14 14 14 15 ...
## $ cylinders : num 8 8 8 8 8 8 8 8 8 8 ...
## $ displacement: num 307 350 318 304 302 429 454 440 455 390 ...
## $ horsepower : num 130 165 150 150 140 198 220 215 225 190 ...
## $ weight : num 3504 3693 3436 3433 3449 ...
## $ acceleration: num 12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
## $ year : num 70 70 70 70 70 70 70 70 70 70 ...
## $ origin : num 1 1 1 1 1 1 1 1 1 1 ...
## $ name : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...head(Auto)## mpg cylinders displacement horsepower weight acceleration year origin
## 1 18 8 307 130 3504 12.0 70 1
## 2 15 8 350 165 3693 11.5 70 1
## 3 18 8 318 150 3436 11.0 70 1
## 4 16 8 304 150 3433 12.0 70 1
## 5 17 8 302 140 3449 10.5 70 1
## 6 15 8 429 198 4341 10.0 70 1
## name
## 1 chevrolet chevelle malibu
## 2 buick skylark 320
## 3 plymouth satellite
## 4 amc rebel sst
## 5 ford torino
## 6 ford galaxie 500# матричные графики разброса переменных
p <- ggpairs(Auto[, c(1, 2:5)])
suppressMessages(print(p))
p <- ggpairs(Auto[, c(1, 6:8)])
suppressMessages(print(p))
# обучающая выборка
set.seed(my.seed)
train <- sample(1:nrow(Auto), nrow(Auto)/2) # обучающая выборка -- 50%Построим дерево регрессии для зависимой переменной mpg: расход топлива, миль на галлон (miles per gallon)
# обучаем модель
tree.auto <- tree(mpg ~ (.-name), data=Auto, subset = train)
summary(tree.auto)##
## Regression tree:
## tree(formula = mpg ~ (. - name), data = Auto, subset = train)
## Variables actually used in tree construction:
## [1] "cylinders" "weight" "year" "displacement"
## Number of terminal nodes: 8
## Residual mean deviance: 7.624 = 1433 / 188
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -7.0000 -1.4530 -0.2411 0.0000 1.5470 10.9000# визуализация
plot(tree.auto)
text(tree.auto, pretty = 0)
Снова сделаем обрезку дерева в целях улучшения качества прогноза.
# обрезка дерева
cv.auto <- cv.tree(tree.auto)
# размер дерева с минимальной ошибкой
plot(cv.auto$size, cv.auto$dev, type = 'b')
opt.size <- cv.auto$size[cv.auto$dev == min(cv.auto$dev)]
abline(v = opt.size, col = 'red', 'lwd' = 2) # соотв. вертикальная прямая
mtext(opt.size, at = opt.size, side = 1, col = 'red', line = 1)
abline(v = 6, col = 'red', 'lwd' = 2, lty=2) # соотв. вертикальная прямая
mtext(6, at = 6, side = 1, col = 'red', line = 1)
В данном случаем минимум ошибки соответствует самому сложному дереву, с 8 узлами. Покажем, как при желании можно обрезать дерево до 6 узлов (ошибка ненамного выше, чем минимальная).
# дерево с 6 узлами
prune.auto = prune.tree(tree.auto, best = 6)
# визуализация
plot(prune.auto)
text(prune.auto, pretty = 0)
Прогноз сначала сделаем по необрезанному дереву, т.к. там ошибка, оцененная по методу перекрёстной проверки, минимальна.
# прогноз по лучшей модели (8 узлов)
yhat <- predict(tree.auto, newdata = Auto[-train, ])
auto.test <- Auto[-train, "mpg"]
# график "прогноз -- реализация"
plot(yhat, auto.test)
# линия идеального прогноза
abline(0, 1)
# MSE на тестовой выборке
mse.test <- mean((yhat - auto.test)^2)
names(mse.test)[length(mse.test)] <- 'Auto.regr.tree.8'
mse.test## Auto.regr.tree.8
## 15.56498MSE на тестовой выборке равна 15.56.
Теперь сделаем прогноз по обрезанному дереву.
# прогноз по модели из 6 узлов
yhat.prune <- predict(prune.auto, newdata = Auto[-train, ])
auto.test <- Auto[-train, "mpg"]
# график "прогноз -- реализация"
plot(yhat.prune, auto.test)
# линия идеального прогноза
abline(0, 1)
# MSE обрезанного дерева на тестовой выборке
mse.prune.test <- mean((yhat.prune - auto.test)^2)
names(mse.prune.test)[length(mse.prune.test)] <- 'Auto.regr.tree.6'
mse.prune.test## Auto.regr.tree.6
## 17.12248MSE обрезанного дерева на тестовой выборке равна 17.12
Источники
- Джеймс Г., Уиттон Д., Хасти Т., Тибширани Р. Введение в статистическое обучение с примерами на языке R / пер. с англ. С.Э. Мастицкого. – М.: ДМК Пресс, 2016 – 450 с. Репозиторий с примерами к книге на русском языке: https://github.com/ranalytics/islr-ru