R Notebook

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

Практика 7

Нелинейные модели

В практических примерах ниже показано как:

оценивать полиномиальную регрессию;
аппроксимировать нелинейные модели ступенчатыми функциями;
строить сплайны;
работать с локальной регрессией;
строить обобщённые линейные модели (GAM).

Модели: полиномиальная регрессия, полиномиальная логистическая регрессия, ступенчатая модель, обобщённая линейная модель.
Данные: Wage {ISLR}

Подробные комментарии к коду лабораторных см. в [1], глава 7.

library('ISLR')              # набор данных Auto
library('splines')           # сплайны
library('gam')               # обобщённые аддитивные модели
library('akima')             # график двумерной плоскости
library('ggplot2')           # красивые графики

my.seed <- 1

Работаем с набором данных по зарплатам 3000 работников-мужчин среднеатлантического региона Wage. Присоединяем его к пространству имён функцией attach(), и дальше обращаемся напрямую к столбцам таблицы.

attach(Wage)

Работаем со столбцами:
* wage – заработная плата работника до уплаты налогов;
* age – возраст работника в годах.

Полиномиальная регрессия

Зависимость зарплаты от возраста

Судя по графику ниже, ззаимосвязь заработной платы и возраста нелинейна. Наблюдается также группа наблюдений с высоким значением wage, граница проходит примерно на уровне 250.

gp <- ggplot(data = Wage, aes(x = age, y = wage))
gp <- gp + geom_point() + geom_abline(slope = 0, intercept = 250, col = 'red')
gp

Подгоняем полином четвёртой степени для зависимости заработной платы от возраста.

fit <- lm(wage ~ poly(age, 4), data = Wage)
round(coef(summary(fit)), 2)

##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)     111.70       0.73  153.28     0.00
## poly(age, 4)1   447.07      39.91   11.20     0.00
## poly(age, 4)2  -478.32      39.91  -11.98     0.00
## poly(age, 4)3   125.52      39.91    3.14     0.00
## poly(age, 4)4   -77.91      39.91   -1.95     0.05

Функция poly(age, 4) создаёт таблицу с базисом ортогональных полиномов: линейные комбинации значений переменной age в степенях от 1 до 4.

round(head(poly(age, 4)), 3)

##           1      2      3      4
## [1,] -0.039  0.056 -0.072  0.087
## [2,] -0.029  0.026 -0.015 -0.003
## [3,]  0.004 -0.015  0.000  0.014
## [4,]  0.001 -0.015  0.005  0.013
## [5,]  0.012 -0.010 -0.011  0.010
## [6,]  0.018 -0.002 -0.017 -0.001

# можно получить сами значения age в заданных степенях
round(head(poly(age, 4, raw = T)), 3)

##       1    2      3       4
## [1,] 18  324   5832  104976
## [2,] 24  576  13824  331776
## [3,] 45 2025  91125 4100625
## [4,] 43 1849  79507 3418801
## [5,] 50 2500 125000 6250000
## [6,] 54 2916 157464 8503056

# на прогноз не повлияет, но оценки параметров изменяются
fit.2 <- lm(wage ~ poly(age, 4, raw = T), data = Wage)
round(coef(summary(fit.2)), 2)

##                        Estimate Std. Error t value Pr(>|t|)
## (Intercept)             -184.15      60.04   -3.07     0.00
## poly(age, 4, raw = T)1    21.25       5.89    3.61     0.00
## poly(age, 4, raw = T)2    -0.56       0.21   -2.74     0.01
## poly(age, 4, raw = T)3     0.01       0.00    2.22     0.03
## poly(age, 4, raw = T)4     0.00       0.00   -1.95     0.05

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

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

# рассчитать прогнозы и их стандартные ошибки
preds <- predict(fit, newdata = list(age = age.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       41.33       62.53
## 2       49.76       67.24
## 3       57.39       71.76
## 4       64.27       76.09
## 5       70.44       80.27
## 6       75.94       84.28

Рисуем левую панель графика со слайда 4 презентации (рис. 7.1 книги). Функция matlines() рисует грфик столбцов одной матрицы против столбцов другой.

# наблюдения
plot(age, wage, xlim = agelims, cex = 0.5, col = 'darkgrey')

# заголовок
title('Полином четвёртой степени')

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

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

Убедимся, что прогнозы по моделям с различными вызовами poly() совпадают.

# прогнозы по второму вызову модели
preds2 <- predict(fit.2, newdata = list(age = age.grid), se = T)

# максимальное расхождение между прогнозами по двум вариантам вызова модели
max(abs(preds$fit - preds2$fit))

## [1] 7.81597e-11

Теперь подбираем степень полинома, сравнивая модели со степенями от 1 до 5 с помощью дисперсионного анализа (ANOVA).

fit.1 <- lm(wage ~ age, data = Wage)
fit.2 <- lm(wage ~ poly(age, 2), data = Wage)
fit.3 <- lm(wage ~ poly(age, 3), data = Wage)
fit.4 <- lm(wage ~ poly(age, 4), data = Wage)
fit.5 <- lm(wage ~ poly(age, 5), data = Wage)

round(anova(fit.1, fit.2, fit.3, fit.4, fit.5), 2)

Рассматриваются пять моделей, в которых степени полинома от age идут по возрастанию. В крайнем правом столбце таблице приводятся p-значения для проверки нулевой гипотезы: текущая модель не даёт статистически значимого сокращения RSS по сравнению с предыдущей моделью. Можно сделать вывод, что степени 3 достаточно, дальнейшее увеличение степени не даёт значимого улучшения качества модели.

Зависимость вероятности получать зарплату > 250 от возраста

Теперь вернёмся к группе наблюдений с высоким wage. Рассмотрим зависимость вероятности того, что величина зарплаты больше 250, от возраста.
Подгоняем логистическую регрессию и делаем прогнозы, для этого используем функцию для оценки обобщённой линейной модели glm() и указываем тип модели binomial:

fit <- glm(I(wage > 250) ~ poly(age, 4), data = Wage, family = 'binomial')

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

# результат - доверительный интервал для вероятности события 
#   "Заработная плата выше 250".   
round(head(se.bands), 3)

##   lower.bound upper.bound
## 1           0       0.002
## 2           0       0.003
## 3           0       0.004
## 4           0       0.005
## 5           0       0.006
## 6           0       0.007

Достраиваем график с 4 слайда презентации (рис. 7.1 книги). Рисуем правую панель.

# сетка для графика (изображаем вероятности, поэтому интервал изменения y мал)
plot(age, I(wage > 250), xlim = agelims, type = 'n', ylim = c(0, 0.2),
     ylab = 'P(Wage > 250 | Age)')

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

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

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

# заголовок
title('Полином четвёртой степени')

Ступенчатые функции

Для начала определим несколько интервалов, на каждом из которых будем моделировать зависимость wage от age своим средним уровнем.

# нарезаем предиктор age на 4 равных интервала
table(cut(age, 4))

## 
## (17.9,33.5]   (33.5,49]   (49,64.5] (64.5,80.1] 
##         750        1399         779          72

# подгоняем линейную модель на интервалах
fit <- lm(wage ~ cut(age, 4), data = Wage)
round(coef(summary(fit)), 2)

##                        Estimate Std. Error t value Pr(>|t|)
## (Intercept)               94.16       1.48   63.79     0.00
## cut(age, 4)(33.5,49]      24.05       1.83   13.15     0.00
## cut(age, 4)(49,64.5]      23.66       2.07   11.44     0.00
## cut(age, 4)(64.5,80.1]     7.64       4.99    1.53     0.13

# прогноз -- это средние по `wage` на каждом интервале
preds.cut <- predict(fit, newdata = list(age = age.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(age, wage, xlim = agelims, cex = 0.5, col = 'darkgrey')

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

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

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

Правая часть графика, для вероятности того, что зарплата выше 250.

fit <- glm(I(wage > 250) ~ cut(age, 4), data = Wage, family = 'binomial')

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

# результат - доверительный интервал для вероятности события 
#   "Заработная плата выше 250".   
round(head(se.bands), 3)

##   lower.bound upper.bound
## 1       0.003       0.016
## 2       0.003       0.016
## 3       0.003       0.016
## 4       0.003       0.016
## 5       0.003       0.016
## 6       0.003       0.016

# сетка для графика (изображаем вероятности, поэтому интервал изменения y мал)
plot(age, I(wage > 250), xlim = agelims, type = 'n', ylim = c(0, 0.2),
     ylab = 'P(Wage > 250 | Age)')

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

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

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

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

Сплайны

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

# кубический сплайн с тремя узлами
fit <- lm(wage ~ bs(age, knots = c(25, 40, 60)), data = Wage)
# прогноз
preds.spl <- predict(fit, newdata = list(age = age.grid), se = T)

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

# 3 узла -- 6 степеней свободы (столбцы матрицы)
dim(bs(age, knots = c(25, 40, 60)))

## [1] 3000    6

# если не указываем узлы явно...
dim(bs(age, df = 6))

## [1] 3000    6

#  они привязываются к квартилям
attr(bs(age, df = 6), 'knots')

##   25%   50%   75% 
## 33.75 42.00 51.00

# натуральный сплайн
fit2 <- lm(wage ~ ns(age, df = 4), data = Wage)
preds.spl2 <- predict(fit2, newdata = list(age = age.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(age, wage, col = 'grey')

# модель кубического сплайна
lines(age.grid, preds.spl$fit, lwd = 2)

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

# натуральный сплайн
lines(age.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(age, wage, xlim = agelims, cex = 0.5, col = 'darkgrey')

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

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

# подгоняем модель с подбором лямбды с помощью перекрёстной проверки
fit2 <- smooth.spline(age, wage, cv = T)

## Warning in smooth.spline(age, wage, cv = T): cross-validation with non-unique
## 'x' values seems doubtful

fit2$df

## [1] 6.794596

# рисуем модель
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)

Локальная регрессия

Строим график со слайда 24 (рис. 7.10).

plot(age, wage, xlim = agelims, cex = 0.5, col = 'darkgrey')

title('Локальная регрессия')

# подгоняем модель c окном 0.2
fit <- loess(wage ~ age, span = 0.2, data = Wage)

# подгоняем модель c окном 0.5
fit2 <- loess(wage ~ age, span = 0.5, data = Wage)

# рисум модели
lines(age.grid, predict(fit, data.frame(age = age.grid)),
      col = 'red', lwd = 2)
lines(age.grid, predict(fit2, data.frame(age = age.grid)),
      col = 'blue', lwd = 2)

# легенда
legend('topright', 
       c('s = 0.2', 's = 0.5'),
       col = c('red', 'blue'), lty = 1, lwd = 2, cex = 0.8)

Обобщённые аддитивные модели (GAM) с непрерывным откликом

Построим GAM на натуральных сплайнах степеней 4 (year), 5 (age) с категориальным предиктором edication.

# GAM на натуральных сплайнах
gam.ns <- gam(wage ~ ns(year, 4) + ns(age, 5) + education, data = Wage)

Также построим модель на сглаживающих сплайнах.

# GAM на сглаживающих сплайнах
gam.m3 <- gam(wage ~ s(year, 4) + s(age, 5) + education, data = Wage)

График со слайда 28 (рис. 7.12).

par(mfrow = c(1, 3))
plot(gam.m3, se = T, col = 'blue')

par(mfrow = c(1, 1))

График со слайда 27 (рис. 7.11).

par(mfrow = c(1, 3))
plot(gam.ns, se = T, col = 'red')

par(mfrow = c(1, 1))

График функции от year похож на прямую. Сделаем ANOVA, чтобы понять, какая степень для year лучше.

gam.m1 <- gam(wage ~ s(age, 5) + education, data = Wage)          # без year
gam.m2 <- gam(wage ~ year + s(age, 5) + education, data = Wage)   # year^1

anova(gam.m1, gam.m2, gam.m3, test = 'F')

Третья модель статистически не лучше второй. Кроме того, один из параметров этой модели незначим.

# сводка по модели gam.m3
summary(gam.m3)

## 
## Call: gam(formula = wage ~ s(year, 4) + s(age, 5) + education, data = Wage)
## Deviance Residuals:
##     Min      1Q  Median      3Q     Max 
## -119.43  -19.70   -3.33   14.17  213.48 
## 
## (Dispersion Parameter for gaussian family taken to be 1235.69)
## 
##     Null Deviance: 5222086 on 2999 degrees of freedom
## Residual Deviance: 3689770 on 2986 degrees of freedom
## AIC: 29887.75 
## 
## Number of Local Scoring Iterations: NA 
## 
## Anova for Parametric Effects
##              Df  Sum Sq Mean Sq F value    Pr(>F)    
## s(year, 4)    1   27162   27162  21.981 2.877e-06 ***
## s(age, 5)     1  195338  195338 158.081 < 2.2e-16 ***
## education     4 1069726  267432 216.423 < 2.2e-16 ***
## Residuals  2986 3689770    1236                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Anova for Nonparametric Effects
##             Npar Df Npar F  Pr(F)    
## (Intercept)                          
## s(year, 4)        3  1.086 0.3537    
## s(age, 5)         4 32.380 <2e-16 ***
## education                            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Работаем с моделью gam.m2.

# прогноз по обучающей выборке
preds <- predict(gam.m2, newdata = Wage)

Также можно использовать в GAM локальные регрессии.

# GAM на локальных регрессиях
gam.lo <- gam(wage ~ s(year, df = 4) + lo(age, span = 0.7) + education, 
              data = Wage)

par(mfrow = c(1, 3))
plot(gam.lo, se = T, col = 'green')

par(mfrow = c(1, 1))

# модель со взаимодействием регрессоров year и age
gam.lo.i <- gam(wage ~ lo(year, age, span = 0.5) + education, data = Wage)
par(mfrow = c(1, 2))
plot(gam.lo.i)

par(mfrow = c(1, 1))

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

Построим логистическую GAM для всероятности того, что wage превышает 250.

gam.lr <- gam(I(wage > 250) ~ year + s(age, df = 5) + education, 
              family = 'binomial', data = Wage)
par(mfrow = c(1, 3))
plot(gam.lr, se = T, col = 'green')

# уровни образования по группам разного достатка
table(education, I(wage > 250))

##                     
## education            FALSE TRUE
##   1. < HS Grad         268    0
##   2. HS Grad           966    5
##   3. Some College      643    7
##   4. College Grad      663   22
##   5. Advanced Degree   381   45

В категории с самым низким уровнем образования нет wage > 250, поэтому убираем её.

gam.lr.s <- gam(I(wage > 250) ~ year + s(age, df = 5) + education,
                family = 'binomial', data = Wage, 
                subset = (education != "1. < HS Grad"))
# график
par(mfrow = c(1, 3))
plot(gam.lr.s, se = T, col = 'green')

detach(Wage)

Упражнение 7

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

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

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

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

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

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

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

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

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

Варианты

Номер варианта	Данные	Зависимая переменная	Объясняющая переменная	Вероятность для второй модели	Метод подгонки моделей
1	`Boston {MASS}`	`crim`	`nox`	\(P(crim>30)\)	Полиномиальный сплайн
2	`Boston {MASS}`	`indus`	`dis`	\(P(indus>16.5)\)	Сглаживающие сплайны
3	`Boston {MASS}`	`crim`	`nox`	\(P(crim>30)\)	Локальная регрессия
4	`Boston {MASS}`	`indus`	`dis`	\(P(indus>16.5)\)	Ступенчатая функция
5	`Boston {MASS}`	`nox`	`dis`	\(P(nox>0.8)\)	Натуральный кубический сплайн
6	`Boston {MASS}`	`indus`	`nox`	\(P(indus>16.5)\)	Полиномиальный сплайн
7	`Auto {ISLR}`	`displacement`	`acceleration`	\(P(displacement>280)\)	Сглаживающие сплайны
8	`Auto {ISLR}`	`displacement`	`weight`	\(P(displacement>280)\)	Локальная регрессия
9	`Auto {ISLR}`	`displacement`	`horsepowers`	\(P(displacement>280)\)	Ступенчатая функция
10	`Auto {ISLR}`	`displacement`	`acceleration`	\(P(displacement>280)\)	Натуральный кубический сплайн
11	`Auto {ISLR}`	`displacement`	`weight`	\(P(displacement>280)\)	Полиномиальный сплайн
12	`Auto {ISLR}`	`displacement`	`horsepowers`	\(P(displacement>280)\)	Сглаживающие сплайны
13	`Boston {MASS}`	`indus`	`dis`	\(P(indus>16.5)\)	Локальная регрессия
14	`Auto {ISLR}`	`displacement`	`weight`	\(P(displacement>280)\)	Ступенчатая функция
15	`Boston {MASS}`	`crim`	`nox`	\(P(crim>30)\)	Натуральный кубический сплайн
16	`Boston {MASS}`	`indus`	`dis`	\(P(indus>16.5)\)	Полиномиальный сплайн
17	`Boston {MASS}`	`crim`	`nox`	\(P(crim>30)\)	Сглаживающие сплайны
18	`Auto {ISLR}`	`displacement`	`horsepowers`	\(P(displacement>280)\)	Локальная регрессия

Источники

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