Задание:

Необходимо построить две модели:
* зависимости непрерывного отклика от одного непрерывного предиктора;
* зависимости вероятности (логит) от одного непрерывного предиктора.

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

  1. Указать смысл переменных модели, метод оценки и настроечный параметр (степень полинома, гиперпараметр \(λ\), ширина окна \(s\), число узлов – в зависимости от метода).

  2. Подогнать модель на всех наблюдениях, меняя значение настроечного параметра.

  3. Обосновать оптимальное значение настроечного параметра подходящим методом (кросс-валидация, \(ANOVA\)).

  4. Сделать прогноз на обучающую выборку по лучшей модели: модельные значения и ошибки прогноза.

  5. Построить график с фактическими наблюдениями, модельной кривой и 95% доверительными интервалами прогноза.

В таблице ниже указаны набор данных, столбцы с переменными для модели и метод подгонки.

Вариант 5

Номер варианта Данные Зависимая переменная Объясняющая переменная Вероятность для второй модели Метод подгонки моделей
5 Boston {MASS} nox dis \(P(nox>0.8)\) Натуральный кубический сплайн

Как сдавать: прислать на почту преподавателя ссылки:

* на html-отчёт с видимыми блоками кода (блоки кода с параметром echo = T), размещённый на rpubs.com.

* на код, генерирующий отчёт, в репозитории на github.com. В текст отчёта включить постановку задачи и ответы на вопросы задания.

Решение

Исходные данные — таблица Boston с данными по стоимости жилья в пригороде Бостона
Зависимая переменная nox (концентрация оксидов азота (частей на 10 млн).)
Объясняющая переменная dis (средневзвешенное расстояние до пяти центров занятости г. Бостона).
Вероятность для второй модели \(P(nox>0.8)\)
Подключаем набор данных Boston

library('splines', quietly = T)           # сплайны
library('gam', quietly = T, verbose = F)               # обобщённые аддитивные модели
## Loaded gam 1.20
library(MASS)
attach(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.7
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 ...
# ?Boston # описание набора данных

Построим кубический сплайн с тремя узлами.

# границы изменения переменной dis
dis_lims <- range(dis)

# значения age, для которых делаем прогноз (от min до max с шагом 1)
dis_grid <- seq(from = dis_lims[1], to = dis_lims[2])

# кубический сплайн с тремя узлами
fit <- lm(nox ~ bs(dis, knots = c(25, 40, 60)), data = Boston)
# прогноз
preds_spl <- predict(fit, newdata = list(dis = dis_grid), se = T)
## Warning in predict.lm(fit, newdata = list(dis = dis_grid), se = T): prediction
## from a rank-deficient fit may be misleading

Теперь построим натуральный по трём узлам. Три узла это 6 степеней свободы. Если функции bs(), которая создаёт матрицу с базисом для полиномиального сплайна, передать только степени свободы, она распределит узлы равномерно. В данном случае это квартили распределения dis.

# 3 узла -- 6 степеней свободы (столбцы матрицы)
dim(bs(dis, knots = c(25, 40, 60)))
## [1] 506   6
# если не указываем узлы явно...
dim(bs(dis, df = 6))
## [1] 506   6
#  они привязываются к квартилям
attr(bs(dis, df = 6), 'knots')
##      25%      50%      75% 
## 2.100175 3.207450 5.188425
# натуральный сплайн
fit2 <- lm(nox ~ ns(dis, df = 4), data = Boston)
preds_spl2 <- predict(fit2, newdata = list(dis = dis_grid), se = T)

График сравнения кубического и натурального сплайнов.

par(mfrow = c(1, 1), mar = c(4.5, 4.5, 1, 8.5), oma = c(0, 0, 0, 0), xpd = T)

# наблюдения
plot(dis, nox, col = 'grey')

# модель кубического сплайна
lines(dis_grid, preds_spl$fit, lwd = 2)

# доверительный интервал
lines(dis_grid, preds_spl$fit + 2*preds_spl$se, lty = 'dashed')
lines(dis_grid, preds_spl$fit - 2*preds_spl$se, lty = 'dashed')

# натуральный сплайн
lines(dis_grid, preds_spl2$fit, col = 'red', lwd = 2)

# легенда
legend("topright", inset = c(-0.7, 0),
       c('Кубический \n с 3 узлами', 'Натуральный'),
       lwd = rep(2, 2), col = c('black', 'red'))

# заголовок
title("Сплайны")

Построим график со слайда 20 (рисунок 7.8 книги).

par(mfrow = c(1, 1), mar = c(4.5, 4.5, 1, 1), oma = c(0, 0, 4, 0))

# наблюдения
plot(dis, nox, xlim = dis_lims, cex = 0.5, col = 'darkgrey')

# заголовок
title('Сглаживающий сплайн')

# подгоняем модель с 16 степенями свободы
fit <- smooth.spline(dis, nox, df = 16)

# подгоняем модель с подбором лямбды с помощью перекрёстной проверки
fit2 <- smooth.spline(dis, nox, cv = T)
## Warning in smooth.spline(dis, nox, cv = T): cross-validation with non-unique 'x'
## values seems doubtful
fit2$df
## [1] 15.42984
# рисуем модель
lines(fit, col = 'red', lwd = 2)
lines(fit2, col = 'blue', lwd = 2)
legend('topright', 
       c('16 df', '6.8 df'),
       col = c('red', 'blue'), lty = 1, lwd = 2, cex = 0.8)

Логистическая GAM

Построим логистическую GAM на натуральном кубическом сплайне ns для всероятности того, что nox превышает 0.8

# knots <- quantile(nox, p = c(0.25, 0.5, 0.75))
# knots <- quantile(dis, p = c(0.25, 0.5, 0.75))
gam_lr <- gam(I(nox > 0.8) ~ ns(dis, knots = c(2.1, 3.2, 5.1)), family = 'binomial', data = Boston)
plot.Gam(gam_lr, se = T, col = 'green')

Источники

  1. Джеймс Г., Уиттон Д., Хасти Т., Тибширани Р. Введение в статистическое обучение с примерами на языке R / пер. с англ. С.Э. Мастицкого. – М.: ДМК Пресс, 2016 – 450 с. Репозиторий с примерами к книге на русском языке: https://github.com/ranalytics/islr-ru