Математическое моделирование. Практика 7
Данила Боковой
11 мая, 2021
Задание:
Необходимо построить две модели:
* зависимости непрерывного отклика от одного непрерывного предиктора;
* зависимости вероятности (логит) от одного непрерывного предиктора.
Для каждой модели:
Указать смысл переменных модели, метод оценки и настроечный параметр (степень полинома, гиперпараметр \(λ\), ширина окна \(s\), число узлов – в зависимости от метода).
Подогнать модель на всех наблюдениях, меняя значение настроечного параметра.
Обосновать оптимальное значение настроечного параметра подходящим методом (кросс-валидация, \(ANOVA\)).
Сделать прогноз на обучающую выборку по лучшей модели: модельные значения и ошибки прогноза.
Построить график с фактическими наблюдениями, модельной кривой и 95% доверительными интервалами прогноза.
В таблице ниже указаны набор данных, столбцы с переменными для модели и метод подгонки.
Вариант 6
| Номер варианта | Данные | Зависимая переменная | Объясняющая переменная | Вероятность для второй модели | Метод подгонки моделей |
|---|---|---|---|---|---|
| 6 | Boston {MASS} |
indus |
nox |
\(P(indus>16.5)\) | Полиномиальный сплайн |
Как сдавать: прислать на почту преподавателя ссылки:
* на html-отчёт с видимыми блоками кода (блоки кода с параметром echo = T), размещённый на rpubs.com.
* на код, генерирующий отчёт, в репозитории на github.com. В текст отчёта включить постановку задачи и ответы на вопросы задания.
Решение
Исходные данные — таблица Boston с данными по стоимости жилья в пригороде Бостона
Зависимая переменная indus (доля акров, не относящихся к розничной торговле)
Объясняющая переменная nox (концентрация оксидов азота (частей на 10 млн).).
Вероятность для второй модели \(P(indus>16.5)\)
Подключаем набор данных Boston
# Графики
library('ggplot2')
# сплайны
library('splines', quietly = T)
# обобщённые аддитивные модели
library('gam', quietly = T, verbose = F)## Loaded gam 1.20# набор данных Boston
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.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Зависимость доли акров от концентрации азота
График ниже показывает, что взаимосвязь переменных nox и indus нелинейна. Наблюдается также группа наблюдений с высоким значением indus, граница проходит примерно на уровне 16.5.
gp <- ggplot(data = Boston, aes(x = nox, y = indus))
gp <- gp + geom_point() + geom_abline(slope = 0, intercept = 16.5, col = 'red')
gp
Подгоняем полином четвёртой степени:.
fit <- lm(indus ~ poly(nox, 4), data = Boston)
round(coef(summary(fit)), 2)## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.14 0.18 62.61 0
## poly(nox, 4)1 117.73 4.00 29.42 0
## poly(nox, 4)2 -32.47 4.00 -8.12 0
## poly(nox, 4)3 -14.12 4.00 -3.53 0
## poly(nox, 4)4 25.14 4.00 6.28 0Функция poly(nox, 4) создаёт таблицу с базисом ортогональных полиномов: линейные комбинации значений переменной nox в степенях от 1 до 4.
round(head(poly(nox, 4)), 3)## 1 2 3 4
## [1,] -0.006 -0.033 0.040 0.006
## [2,] -0.033 0.003 0.028 -0.049
## [3,] -0.033 0.003 0.028 -0.049
## [4,] -0.037 0.011 0.019 -0.047
## [5,] -0.037 0.011 0.019 -0.047
## [6,] -0.037 0.011 0.019 -0.047# можно получить сами значения nox в заданных степенях
round(head(poly(nox, 4, raw = T)), 3)## 1 2 3 4
## [1,] 0.538 0.289 0.156 0.084
## [2,] 0.469 0.220 0.103 0.048
## [3,] 0.469 0.220 0.103 0.048
## [4,] 0.458 0.210 0.096 0.044
## [5,] 0.458 0.210 0.096 0.044
## [6,] 0.458 0.210 0.096 0.044# на прогноз не повлияет, но оценки параметров изменяются
fit_2 <- lm(indus ~ poly(nox, 4, raw = T), data = Boston)
round(coef(summary(fit_2)), 2)## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 682.67 107.26 6.36 0
## poly(nox, 4, raw = T)1 -4726.47 730.65 -6.47 0
## poly(nox, 4, raw = T)2 11889.17 1827.01 6.51 0
## poly(nox, 4, raw = T)3 -12762.64 1987.03 -6.42 0
## poly(nox, 4, raw = T)4 4981.94 793.08 6.28 0# границы изменения переменной nox
nox_lims <- range(nox)
# значения nox, для которых делаем прогноз (от min до max с шагом 0.05)
nox_grid <- seq(from = nox_lims[1], to = nox_lims[2], by=0.05)
# рассчитать прогнозы и их стандартные ошибки
preds <- predict(fit, newdata = list(nox = nox_grid), se = T)
# границы доверительного интервала
se_bands <- cbind(lower_bound = preds$fit - 2*preds$se.fit,
upper_bound = preds$fit + 2*preds$se.fit)
# смотрим результат
round(head(se_bands), 2)## lower_bound upper_bound
## 1 4.57 8.21
## 2 3.64 4.85
## 3 5.98 7.23
## 4 10.24 11.34
## 5 14.19 15.48
## 6 16.86 18.23Рисуем левую панель графика со слайда 4 презентации (рис. 7.1 книги). Функция matlines() рисует грфик столбцов одной матрицы против столбцов другой.
# наблюдения
plot(nox, indus, xlim = nox_lims, cex = 0.5, col = 'darkgrey')
# заголовок
title('Полином четвёртой степени')
# модель
lines(nox_grid, preds$fit, lwd = 2, col = 'blue')
# доверительные интервалы прогноза
matlines(nox_grid, se_bands, lwd = 1, col = 'blue', lty = 3)
Убедимся, что прогнозы по моделям с различными вызовами poly() совпадают.
# прогнозы по второму вызову модели
preds2 <- predict(fit_2, newdata = list(nox = nox_grid), se = T)
# максимальное расхождение между прогнозами по двум вариантам вызова модели
max(abs(preds$fit - preds2$fit))## [1] 4.835243e-12Теперь подбираем степень полинома, сравнивая модели со степенями от 1 до 5 с помощью дисперсионного анализа (ANOVA).
fit1 <- lm(indus ~ poly(nox, 1), data = Boston)
fit2 <- lm(indus ~ poly(nox, 2), data = Boston)
fit3 <- lm(indus ~ poly(nox, 3), data = Boston)
fit4 <- lm(indus ~ poly(nox, 4), data = Boston)
fit5 <- lm(indus ~ poly(nox, 5), data = Boston)
round(anova(fit1, fit2, fit3, fit4, fit5), 2)## Analysis of Variance Table
##
## Model 1: indus ~ poly(nox, 1)
## Model 2: indus ~ poly(nox, 2)
## Model 3: indus ~ poly(nox, 3)
## Model 4: indus ~ poly(nox, 4)
## Model 5: indus ~ poly(nox, 5)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 504 9907.2
## 2 503 8852.8 1 1054.42 66.64 <2e-16 ***
## 3 502 8653.5 1 199.25 12.59 <2e-16 ***
## 4 501 8021.7 1 631.82 39.93 <2e-16 ***
## 5 500 7911.2 1 110.49 6.98 0.01 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anovas <- anova(fit1, fit2, fit3, fit4, fit5)
#Библиотека для оформления таблиц
library(kableExtra)
kable(anovas, digits = 2, format = "simple", caption = "Табл. 1 — Дисперсионный анализ (ANOVA)") %>% kable_styling()## Warning in kable_styling(.): Please specify format in kable. kableExtra can
## customize either HTML or LaTeX outputs. See https://haozhu233.github.io/
## kableExtra/ for details.| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 504 | 9907.18 | NA | NA | NA | NA |
| 503 | 8852.76 | 1 | 1054.42 | 66.64 | 0.00 |
| 502 | 8653.51 | 1 | 199.25 | 12.59 | 0.00 |
| 501 | 8021.69 | 1 | 631.82 | 39.93 | 0.00 |
| 500 | 7911.19 | 1 | 110.49 | 6.98 | 0.01 |
#Оценим критерий Акаике у пяти моделей и построим график
degree <- 5
aics <- numeric(degree)
for (i in 1:degree) aics[i] <- AIC(lm(indus ~ poly(nox, i), data = Boston))
plot(1:degree, aics, xlab = 'Степень полинома', ylab = 'Информационный критерий Акаике (AIC)', type = 'l')
Рассматриваются пять моделей, в которых степени полинома от nox идут по возрастанию.
В крайнем правом столбце таблице приводятся p-значения для проверки нулевой гипотезы:
текущая модель не даёт статистически значимого сокращения RSS по сравнению с предыдущей моделью.
Можно сделать вывод, что степени 4 достаточно, дальнейшее увеличение степени не даёт значимого улучшения качества модели.
Зависимость вероятности переменной indus > 16.5 от nox
Рассмотрим зависимость вероятности того, что доля акров, не относящихся к розничной торговле больше 16,5, от концентрации оксидов азота.
Подгоняем логистическую регрессию и делаем прогнозы, для этого используем функцию для оценки обобщённой линейной модели glm() и указываем тип модели binomial:
fit_logit <- glm(I(indus > 16.5) ~ poly(nox, 4), data = Boston, family = 'binomial')## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred# прогнозы
preds <- predict(fit_logit, newdata = list(nox = nox_grid), se = T)
# пересчитываем доверительные интервалы и прогнозы в исходные ЕИ
pfit <- exp(preds$fit) / (1 + exp(preds$fit))
se_bands_logit <- cbind(lower_bound = preds$fit - 2*preds$se.fit,
upper_bound = preds$fit + 2*preds$se.fit)
se_bands <- exp(se_bands_logit)/(1 + exp(se_bands_logit))
# результат - доверительный интервал для вероятности события
# "indus выше 16.5".
round(head(se_bands), 3)## lower_bound upper_bound
## 1 0.000 0.999
## 2 0.000 0.000
## 3 0.000 0.001
## 4 0.013 0.118
## 5 0.622 0.843
## 6 0.645 0.896Достраиваем график с 4 слайда презентации (рис. 7.1 книги). Рисуем правую панель.
# сетка для графика (изображаем вероятности, поэтому интервал изменения y мал)
plot(nox, I(indus > 16.5), xlim = nox_lims, type = 'n', ylab = 'P(indus > 16.5 | nox)')
# фактические наблюдения показываем засечками
points(nox, I((indus > 16.5) / 1), cex = 0.5, pch = '|', col = 'darkgrey')
# модель
lines(nox_grid, pfit, lwd = 2, col = 'blue')
# доверительные интервалы
matlines(nox_grid, se_bands, lwd = 1, col = 'blue', lty = 3)
# заголовок
title('Полином четвёртой степени')
Источники
- Джеймс Г., Уиттон Д., Хасти Т., Тибширани Р. Введение в статистическое обучение с примерами на языке R / пер. с англ. С.Э. Мастицкого. – М.: ДМК Пресс, 2016 – 450 с. Репозиторий с примерами к книге на русском языке: https://github.com/ranalytics/islr-ru