Упражнение 7
Нестерова А.И.
11 04 2020
Математическое моделирование
Нелинейные модели
Необходимо построить две модели:
- зависимости непрерывного отклика от одного непрерывного предиктора;
- зависимости вероятности (логит) от одного непрерывного предиктора.
Для каждой модели:
1 Указать смысл переменных модели, метод оценки и настроечный параметр (степень полинома, гиперпараметр \(λ\), ширина окна \(s\), число узлов – в зависимости от метода).
2 Подогнать модель на всех наблюдениях, меняя значение настроечного параметра.
3 Обосновать оптимальное значение настроечного параметра подходящим методом (кросс-валидация, ANOVA).
4 Сделать прогноз на обучающую выборку по лучшей модели: модельные значения и ошибки прогноза.
5 Построить график с фактическими наблюдениями, модельной кривой и 95% доверительными интервалами прогноза.
В таблице ниже указаны набор данных, столбцы с переменными для модели и метод подгонки.
Как сдавать: прислать на почту преподавателя ссылки:
- на html-отчёт с видимыми блоками кода (блоки кода с параметром echo = T), размещённый на rpubs.com.
- на код, генерирующий отчёт, в репозитории на github.com. В текст отчёта включить постановку задачи и ответы на вопросы задания.
Вариант - 13
Модели: локальная регрессия.
Данные: `Boston {MASS}’.
library('MASS') # набор данных Boston
library('splines') # сплайны
library('gam') # обобщённые аддитивные модели
library('akima') # график двумерной плоскости
library('ggplot2') # красивые графики
# загрузка данных Boston
data('Boston')
# ядро
my.seed <- 1Работаем с набором данных по стоимости жилья в пригороде Бостона. Присоединяем его к пространству имён функцией attach(), и дальше обращаемся напрямую к столбцам таблицы.
attach(Boston)Работаем со столбцами:
* indus – доля нерозничных торговых площадей на город;
* dis – средневзвешенное расстояние до пяти бостонских центров занятости.
Локальная регрессия
Судя по графику ниже, взаимосвязь доли торговых площадей и средневзвешенного расстояния до центров занятости нелинейна. Наблюдается также группа наблюдений с высоким значением indus, граница проходит примерно на уровне 18.
gp <- ggplot(data = Boston, aes(x = dis, y = indus))
gp <- gp + geom_point() + geom_abline(slope = 0, intercept = 16.5, col = 'red')
gp
Зависимость доли торговых площадей от средневзвешенного расстояния до центров занятости (модель 1)
# границы изменения переменной dis
dislims <- range(dis)
# значения dis, для которых делаем прогноз (от min до max с шагом 1)
dis.grid <- seq(from = dislims[1], to = dislims[2])
# график
plot(dis, indus, xlim = dislims, cex = 0.5, col = 'darkgrey')
title('Локальная регрессия')
# подгоняем модель c окном 0.2
fit1 <- loess(indus ~ dis, cv = 0.2, data = Boston)
# подгоняем модель c окном 0.7
fit2 <- loess(indus ~ dis, span = 0.7, data = Boston)
# рисум модели
lines(dis.grid, predict(fit1, data.frame(dis = dis.grid)),
col = 'red', lwd = 2)
lines(dis.grid, predict(fit2, data.frame(dis = dis.grid)),
col = 'blue', lwd = 2)
# легенда
legend('topright',
c('s = 0.2', 's = 0.7'),
col = c('red', 'blue'), lty = 1, lwd = 2, cex = 0.8)
Определение оптимального настроечного параметра (модель 1)
Подбор настроечного параметра (ширина окна \(s\)) осуществим с помощью кросс-валидации. Для того, чтобы параметр был оптимальным, необходимо минимизировать ошибку СV. Чем меньше окажется ширина окна, тем более локальной и извилистой будет модедь; в то же время очень большое значение параметра \(s\) приведёт к глобальной модели, построенной по всей обучающей выборке.
# подгоняем модель с подбором ширины окна с помощью перекрёстной проверки
fit <- loess(indus ~ dis, cv = T, data = Boston)
# параметр span
fit$pars## $span
## [1] 0.75
##
## $degree
## [1] 2
##
## $normalize
## [1] TRUE
##
## $parametric
## [1] FALSE
##
## $drop.square
## [1] FALSE
##
## $surface
## [1] "interpolate"
##
## $cell
## [1] 0.2
##
## $family
## [1] "gaussian"
##
## $trace.hat
## [1] "exact"
##
## $iterations
## [1] 1# s=0.75
# лучшая модель (с окном 0.75)
best.fit1 <- loess(indus ~ dis, span = 0.75, data = Boston)Прогноз на обучающую выборку по лучшей модели (модель 1)
# прогноз по лучшей модели
preds <- predict(best.fit1, data.frame(dis = dis.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 20.07 23.15
## 2 15.10 16.26
## 3 8.95 10.49
## 4 6.69 8.28
## 5 5.83 7.42
## 6 4.69 6.27# стандартные ошибки
round(preds$se.fit, 2)## 1 2 3 4 5 6 7 8 9 10 11
## 0.77 0.29 0.39 0.40 0.40 0.40 0.42 0.50 0.71 1.11 1.71График с фактическими наблюдениями, модельной кривой и 95% доверительными интервалами прогноза (модель 1)
# сетка для графика
plot(dis, indus, xlim = dislims, type = 'n', ylab = 'P(Indus | Dis)')
# фактические наблюдения показываем засечками
points(jitter(dis), indus, cex = 0.5, pch = '|', col = 'darkgrey')
# модель
pfit <- preds$fit
lines(dis.grid, pfit, lwd = 2, col = 'blue')
# доверительные интервалы
matlines(dis.grid, se.bands, lwd = 1, col = 'blue', lty = 3)
# заголовок
title('Локальная регрессия (модель 1)')
Зависимость вероятности доли торговых площадей > 16.5 от средневзвешенного расстояния до центров занятости (модель 2)
# график
plot(dis, I(indus>16.5), xlim = dislims, cex = 0.5, col = 'darkgrey')
title('Локальная регрессия')
# подгоняем модель c окном 0.2
fit3 <- loess(I(indus>16.5) ~ dis, cv = 0.2, data = Boston)
# подгоняем модель c окном 0.7
fit4 <- loess(I(indus>16.5) ~ dis, span = 0.7, data = Boston)
# рисум модели
lines(dis.grid, predict(fit3, data.frame(dis = dis.grid)),
col = 'red', lwd = 2)
lines(dis.grid, predict(fit4, data.frame(dis = dis.grid)),
col = 'blue', lwd = 2)
# легенда
legend('topright',
c('s = 0.2', 's = 0.7'),
col = c('red', 'blue'), lty = 1, lwd = 2, cex = 0.8)
Определение оптимального настроечного параметра (модель 2)
# подгоняем модель с подбором ширины окна с помощью перекрёстной проверки
fit <- loess(I(indus>16.5) ~ dis, cv = T, data = Boston)
# параметр span
fit$pars## $span
## [1] 0.75
##
## $degree
## [1] 2
##
## $normalize
## [1] TRUE
##
## $parametric
## [1] FALSE
##
## $drop.square
## [1] FALSE
##
## $surface
## [1] "interpolate"
##
## $cell
## [1] 0.2
##
## $family
## [1] "gaussian"
##
## $trace.hat
## [1] "exact"
##
## $iterations
## [1] 1best.fit2 <- loess(I(indus>16.5) ~ dis, span = 0.75, data = Boston)Прогноз на обучающую выборку по лучшей модели (модель 2)
# прогнозы
preds <- predict(best.fit2, data.frame(dis = dis.grid), se = T)
# пересчитываем доверительные интервалы и прогнозы в исходные ЕИ
pfit <- exp(preds$fit) / (1 + exp(preds$fit))
se.bands.loess <- cbind(lower.bound = preds$fit - 2*preds$se.fit,
upper.bound = preds$fit + 2*preds$se.fit)
se.bands <- exp(se.bands.loess)/(1 + exp(se.bands.loess))
# результат - доверительный интервал для вероятности события "доля торговой площади больше 16.5".
round(head(se.bands), 3)## lower.bound upper.bound
## 1 0.755 0.791
## 2 0.660 0.677
## 3 0.544 0.569
## 4 0.491 0.518
## 5 0.482 0.509
## 6 0.481 0.507# стандартные ошибки
round(preds$se.fit, 2)## 1 2 3 4 5 6 7 8 9 10 11
## 0.05 0.02 0.03 0.03 0.03 0.03 0.03 0.03 0.05 0.07 0.11График с фактическими наблюдениями, модельной кривой и 95% доверительными интервалами прогноза (модель 2)
# сетка для графика (изображаем вероятности, поэтому интервал изменения y мал)
plot(dis, I(indus > 16.5), xlim = dislims, type = 'n', ylim = c(0, 1),
ylab = 'P(Indus > 16.5 | Dis)')
# фактические наблюдения показываем засечками
points(jitter(dis), I((indus > 16.5) ), cex = 0.5, pch = '|', col = 'darkgrey')
# модель
lines(dis.grid, pfit, lwd = 2, col = 'blue')
# доверительные интервалы
matlines(dis.grid, se.bands, lwd = 1, col = 'blue', lty = 3)
# заголовок
title('Локальная регрессия (модель 2)')