Задание:

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

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

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

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

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

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

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

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

Вариант 9

Номер варианта Данные Зависимая переменная Объясняющая переменная Вероятность для второй модели Метод подгонки моделей
9 Auto {ISLR} displacement horsepowers \(P(displacement>280)\) Ступенчатая функция

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

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

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

Решение

Подключаем набор данных Auto

library('ISLR')
attach(Auto)
data(Auto)

Зависимая переменная displacement (объём двигателя в кубических дюймах) Объясняющая переменная horsepowers Вероятность для второй модели P(displacement>280)

# нарезаем предиктор horsepower на 4 равных интервала
table(cut(horsepower, 4))
## 
## (45.8,92]  (92,138] (138,184] (184,230] 
##       195       111        69        17
# подгоняем линейную модель на интервалах
fit <- lm(displacement ~ cut(horsepower, 4), data = Auto)
round(coef(summary(fit)), 2)
##                             Estimate Std. Error t value Pr(>|t|)
## (Intercept)                   120.89       3.65   33.15        0
## cut(horsepower, 4)(92,138]     79.82       6.05   13.18        0
## cut(horsepower, 4)(138,184]   219.55       7.13   30.78        0
## cut(horsepower, 4)(184,230]   282.93      12.88   21.97        0
horsepowerlims = range(horsepower)
horsepower.grid = seq(from = horsepowerlims[1], to = horsepowerlims[2])

# прогноз -- это средние по `displacement` на каждом интервале
preds.cut <-predict(fit,newdata = list(horsepower = horsepower.grid),se = T)

# интервальный прогноз
se.bands.cut <-cbind(
    lower.bound = preds.cut$fit - 2 * preds.cut$se.fit,
    upper.bound = preds.cut$fit + 2 * preds.cut$se.fit
  )

Воспроизведём график со слайда 7 презентации (рис. 7.2 книги).

# наблюдения
plot(horsepower,displacement,xlim = horsepowerlims,cex = 0.5,col = 'darkgrey')

# модель
lines(horsepower.grid,preds.cut$fit,lwd = 2,col = 'darkgreen')

# доверительные интервалы прогноза
matlines(x = horsepower.grid,y = se.bands.cut,lwd = 1,col = 'darkgreen',lty = 3)

# заголовок
title('Ступенчатая функция')

Правая часть графика, для вероятности того, что объём двигателя выше 280.

fit <-glm(I(displacement > 280) ~ cut(horsepower, 4),data = Auto,family = 'binomial')

# прогнозы
preds <-predict(fit,newdata = list(horsepower = horsepower.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))

# результат - доверительный интервал для вероятности события "Объём двигателя выше 280".
round(head(se.bands), 3)
##   lower.bound upper.bound
## 1           0         NaN
## 2           0         NaN
## 3           0         NaN
## 4           0         NaN
## 5           0         NaN
## 6           0         NaN
# сетка для графика (изображаем вероятности, поэтому интервал изменения y мал)
plot(horsepower,I(displacement > 280),xlim = horsepowerlims,type = 'n',ylab = 'P(displacement > 280 | horsepower)')

# фактические наблюдения показываем засечками
points(jitter(horsepower),I((displacement > 280) / 5),cex = 0.5,pch = '|',col = 'darkgrey')

# модель
lines(horsepower.grid, pfit, lwd = 2, col = 'darkgreen')

# доверительные интервалы
matlines(horsepower.grid,se.bands,lwd = 1,col = 'darkgreen',lty = 3)

# заголовок
title('Ступенчатая функция')

Источники

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