Математическое моделирование

Модели на основе деревьев

Необходимо построить две модели для прогноза на основе дерева решений:

  • для непрерывной зависимой переменной;
  • для категориальной зависимой переменной.

Данные и переменные указаны в таблице с вариантами.
Ядро генератора случайных чисел – номер варианта.

Задания

Для каждой модели:

  1. Указать настроечные параметры метода из своего варианта (например: количество узлов, количество предикторов, скорость обучения).
  2. Подогнать модель на обучающей выборке (50% наблюдений). Рассчитать MSE на тестовой выборке.
  3. Перестроить модель с помощью метода, указанного в варианте.
  4. Сделать прогноз по модели с подобранными в п.3 параметрами на тестовой выборке, оценить его точность и построить график «прогноз-реализация».

Как сдавать: прислать на почту преподавателя ссылки: * на html-отчёт с видимыми блоками кода (блоки кода с параметром echo = T), размещённый на rpubs.com.
* на код, генерирующий отчёт, в репозитории на github.com. В текст отчёта включить постановку задачи и ответы на вопросы задания.

Вариант - 13

Модели: дерево с обрезкой ветвей (настроечный параметр: количество узлов).
Данные: `Boston {MASS}’.

Деревья решений

# Загрузка пакетов
library('tree')              # деревья tree()
library('GGally')            # матричный график разброса ggpairs()
library('MASS')              # набор данных Boston

# загрузка данных Boston
data('Boston')

# название столбцов переменных
names(Boston)
##  [1] "crim"    "zn"      "indus"   "chas"    "nox"     "rm"      "age"    
##  [8] "dis"     "rad"     "tax"     "ptratio" "black"   "lstat"   "medv"
# размерность данных
dim(Boston)
## [1] 506  14
# ядро генератора случайных чисел
my.seed <- 13

Модель 1 (для непрерывной зависимой переменной crim)

# ?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.7
# матричные графики разброса переменных
p <- ggpairs(Boston[, c(1, 2:5)])
suppressMessages(print(p))

p <- ggpairs(Boston[, c(1, 6:9)])
suppressMessages(print(p))

p <- ggpairs(Boston[, c(1, 10:14)])
suppressMessages(print(p))

# обучающая выборка
set.seed(my.seed)
train <- sample(1:nrow(Boston), nrow(Boston)/2) # обучающая выборка -- 50%

Построим дерево регрессии для зависимой переменной crim: уровень преступности на душу населения.

# обучаем модель
tree.boston <- tree(crim ~ ., Boston, subset = train)
summary(tree.boston)
## 
## Regression tree:
## tree(formula = crim ~ ., data = Boston, subset = train)
## Variables actually used in tree construction:
## [1] "medv"  "dis"   "lstat" "rad"  
## Number of terminal nodes:  7 
## Residual mean deviance:  25.51 = 6276 / 246 
## Distribution of residuals:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -25.5200  -0.3585  -0.2628   0.0000   0.1355  46.6400
# визуализация
plot(tree.boston)
text(tree.boston, pretty = 0)

tree.boston                    # посмотреть всё дерево в консоли
## node), split, n, deviance, yval
##       * denotes terminal node
## 
##  1) root 253 23830.00  4.2350  
##    2) medv < 10.9 23  8880.00 25.3900  
##      4) dis < 1.7291 15  6619.00 32.1300  
##        8) dis < 1.45335 5  4575.00 42.3300 *
##        9) dis > 1.45335 10  1264.00 27.0300  
##         18) lstat < 27.075 5   267.90 19.5400 *
##         19) lstat > 27.075 5   434.50 34.5200 *
##      5) dis > 1.7291 8   297.80 12.7400 *
##    3) medv > 10.9 230  3632.00  2.1200  
##      6) rad < 16 182    70.00  0.3986 *
##      7) rad > 16 48   979.30  8.6450  
##       14) dis < 2.25885 26   562.50 11.1300 *
##       15) dis > 2.25885 22    67.76  5.7130 *
# прогноз по модели 
yhat <- predict(tree.boston, newdata = Boston[-train, ])
boston.test <- Boston[-train, "crim"]

# MSE на тестовой выборке
mse.test <- mean((yhat - boston.test)^2)
names(mse.test)[length(mse.test)] <- 'Boston.regr.tree.all'
mse.test
## Boston.regr.tree.all 
##             33.30986
# точность прогноза на тестовой выборке
acc.test <- sum(abs(yhat-boston.test))/sum(boston.test)
names(acc.test)[length(acc.test)] <- 'Boston.regr.tree.all'
acc.test
## Boston.regr.tree.all 
##             0.617652

Дерево с обрезкой ветвей (модель 1)

Cделаем обрезку дерева в целях улучшения качества прогноза.

# обрезка дерева
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)

Как видно на графике, минимум частоты ошибок достигается при числе узлов 3. Оценим точность дерева с 3 узлами.

# дерево с 3 узлами
prune.boston = prune.tree(tree.boston, best = 3)

# визуализация
plot(prune.boston)
text(prune.boston, pretty = 0)

# прогноз по лучшей модели (3 узла)
yhat <- predict(prune.boston, newdata = Boston[-train, ])
boston.test <- Boston[-train, "crim"]

# MSE на тестовой выборке (3 узла)
mse.test <- c(mse.test, mean((yhat - boston.test)^2))
names(mse.test)[length(mse.test)] <- 'Boston.regr.tree.3'
mse.test
## Boston.regr.tree.all   Boston.regr.tree.3 
##             33.30986             31.85253
# точность прогноза на тестовой выборке (3 узла)
acc.test <- c(acc.test, sum(abs(yhat-boston.test))/sum(boston.test))
names(acc.test)[length(acc.test)] <- 'Boston.regr.tree.3'
acc.test
## Boston.regr.tree.all   Boston.regr.tree.3 
##            0.6176520            0.6839895
# график "прогноз -- реализация"
plot(yhat, boston.test)
# линия идеального прогноза
abline(0, 1)

MSE модели (3 узла) на тестовой выборке равна 31.85, точность прогноза составила 0.68

Модель 2 (для категориальной зависимой переменной high.crim)

Загрузим таблицу с данными по стоимости жилья в пригороде Бостона и добавим к ней переменную high.crim – высокий уровень преступности на душу населения со значениями:

  • 1, если продажи больше 3.5;
  • 0 - в противном случае.
# новая переменная
high.crim <- ifelse(Boston$crim <= 3.5, '0', '1')

# присоединяем к таблице данных
Boston<- cbind(Boston, high.crim)

# название столбцов переменных
names(Boston)
##  [1] "crim"      "zn"        "indus"     "chas"      "nox"       "rm"       
##  [7] "age"       "dis"       "rad"       "tax"       "ptratio"   "black"    
## [13] "lstat"     "medv"      "high.crim"
# размерность данных
dim(Boston)
## [1] 506  15
# матричные графики разброса переменных
p <- ggpairs(Boston[, c(15, 1:5)], aes(color = high.crim))
suppressMessages(print(p))

p <- ggpairs(Boston[, c(15, 6:10)], aes(color = high.crim))
suppressMessages(print(p))

p <- ggpairs(Boston[, c(15, 11:14)], aes(color = high.crim))
suppressMessages(print(p))

Судя по графикам, класс 0 превосходит по размеру класс1 по переменной high.crim приблизительно в 3 раза. Классы на графиках разброса объясняющих переменных сильно смешаны, поэтому модели с непрерывной разрешающей границей вряд ли работают хорошо. Построим дерево для категориального отклика high.crim, отбросив непрерывный отклик crim (мы оставили его на первом графике, чтобы проверить, как сработало разделение по значению crim = 3.5).

# модель бинарного  дерева
tree.boston <- tree(high.crim ~ . -crim, Boston)
summary(tree.boston)
## 
## Classification tree:
## tree(formula = high.crim ~ . - crim, data = Boston)
## Variables actually used in tree construction:
## [1] "rad" "nox" "age" "rm" 
## Number of terminal nodes:  5 
## Residual mean deviance:  0.04749 = 23.79 / 501 
## Misclassification error rate: 0.009881 = 5 / 506
# график результата
plot(tree.boston)              # ветви
text(tree.boston, pretty = 0)  # подписи

tree.boston                    # посмотреть всё дерево в консоли
## node), split, n, deviance, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 506 576.600 0 ( 0.743083 0.256917 )  
##    2) rad < 16 374  24.910 0 ( 0.994652 0.005348 )  
##      4) nox < 0.759 358   0.000 0 ( 1.000000 0.000000 ) *
##      5) nox > 0.759 16  12.060 0 ( 0.875000 0.125000 ) *
##    3) rad > 16 132  35.850 1 ( 0.030303 0.969697 )  
##      6) age < 54.95 5   6.730 0 ( 0.600000 0.400000 ) *
##      7) age > 54.95 127  11.680 1 ( 0.007874 0.992126 )  
##       14) rm < 6.998 122   0.000 1 ( 0.000000 1.000000 ) *
##       15) rm > 6.998 5   5.004 1 ( 0.200000 0.800000 ) *

Теперь построим дерево на обучающей выборке и оценим ошибку на тестовой.

# тестовая выборка
Boston.test <- Boston[-train,]
high.crim.test <- high.crim[-train]

# строим дерево на обучающей выборке
tree.boston <- tree(high.crim ~ . -crim, Boston, subset = train)

# делаем прогноз
tree.pred <- predict(tree.boston, Boston.test, type = "class")

# матрица неточностей
tbl <- table(tree.pred, high.crim.test)
tbl
##          high.crim.test
## tree.pred   0   1
##         0 189   2
##         1   3  59
# 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.9802372

Обобщённая характеристика точности: доля верных прогнозов: 0.98.

Дерево с обрезкой ветвей (модель 2)

Теперь обрезаем дерево, используя в качестве критерия частоту ошибок классификации. Функция 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] 3 2 1
## 
## $dev
## [1]  1  1 69
## 
## $k
## [1] -Inf    0   68
## 
## $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)     # соотв. вертикальная прямая


# 2. ошибка с кросс-валидацией в зависимости от штрафа на сложность
plot(cv.boston$k, cv.boston$dev, type = "b",
     ylab = 'Частота ошибок с кросс-вал. (dev)',
     xlab = 'Штраф за сложность (k)')

Как видно на графике слева, минимум частоты ошибок достигается при числе узлов 2 и 3. Оценим точность дерева с 2 узлами.

# дерево с 2 узлами
prune.boston = prune.tree(tree.boston, best = 2)

# визуализация
plot(prune.boston)
text(prune.boston, pretty = 0)

# прогноз на тестовую выборку
tree.pred <- predict(prune.boston, Boston.test, type = "class")

# матрица неточностей
tbl <- table(tree.pred, high.crim.test)
tbl
##          high.crim.test
## tree.pred   0   1
##         0 189   2
##         1   3  59
# ACC на тестовой
acc.test <- c(acc.test, sum(diag(tbl))/sum(tbl))
names(acc.test)[length(acc.test)] <- 'Boston.class.tree.2'
acc.test
## Boston.class.tree.all   Boston.class.tree.2 
##             0.9802372             0.9802372
# график "прогноз -- реализация"
plot(tree.pred, Boston$high.crim[-train])

Точности моделей на тестовой выборке (при двух узлах и максимальном числе узлов) совпадают и равны 0.98.