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

Практика 7

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

Модели: полиномиальная регрессия, полиномиальная логистическая регрессия, обобщённая линейная модель.

Данные: Boston {MASS}- статистика стоимости жилья в пригороде Бостона.

library('MASS')              # набор данных Auto
library('splines')           # сплайны
library('gam')               # обобщённые аддитивные модели
library('akima')             # график двумерной плоскости
library('ggplot2')           # красивые графики
my.seed <- 1
attach(Boston)

Работаем со столбцами:
* crim – уровень преступности на душу населения по городам;
* nox – концентрация оксидов азота (частей на 10 миллионов).

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

Подгоняем полином четвёртой степени для зависимости уровень преступности на душу населения по городам от концентрации оксидов азота.

fit <- lm(crim ~ poly(nox, 4), data = Boston)
coef(summary(fit))
##                 Estimate Std. Error    t value     Pr(>|t|)
## (Intercept)     3.613524  0.3218873 11.2260528 3.065306e-26
## poly(nox, 4)1  81.372015  7.2406752 11.2381806 2.746929e-26
## poly(nox, 4)2 -28.828594  7.2406752 -3.9814787 7.862697e-05
## poly(nox, 4)3 -60.361894  7.2406752 -8.3365008 7.426202e-16
## poly(nox, 4)4  -1.027048  7.2406752 -0.1418442 8.872601e-01

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

head(poly(nox,4))
##                 1            2          3            4
## [1,] -0.006411247 -0.032904358 0.03992435  0.006469961
## [2,] -0.032908669  0.003246481 0.02812126 -0.048609604
## [3,] -0.032908669  0.003246481 0.02812126 -0.048609604
## [4,] -0.037132896  0.011480916 0.01888395 -0.046792235
## [5,] -0.037132896  0.011480916 0.01888395 -0.046792235
## [6,] -0.037132896  0.011480916 0.01888395 -0.046792235
# можно получить сами значения age в заданных степенях
head(poly(nox, 4, raw = T))
##          1        2          3          4
## [1,] 0.538 0.289444 0.15572087 0.08377783
## [2,] 0.469 0.219961 0.10316171 0.04838284
## [3,] 0.469 0.219961 0.10316171 0.04838284
## [4,] 0.458 0.209764 0.09607191 0.04400094
## [5,] 0.458 0.209764 0.09607191 0.04400094
## [6,] 0.458 0.209764 0.09607191 0.04400094
# на прогноз не повлияет, но оценки параметров изменяются
fit.2 <- lm(crim ~ poly(nox, 4, raw = T))
coef(summary(fit.2))
##                          Estimate Std. Error    t value  Pr(>|t|)
## (Intercept)              205.9748   194.0815  1.0612803 0.2890739
## poly(nox, 4, raw = T)1 -1093.4021  1322.1291 -0.8270010 0.4086302
## poly(nox, 4, raw = T)2  1781.2912  3306.0214  0.5388021 0.5902628
## poly(nox, 4, raw = T)3  -736.1309  3595.5839 -0.2047319 0.8378647
## poly(nox, 4, raw = T)4  -203.5596  1435.0931 -0.1418442 0.8872601
# границы изменения переменной age
noxlims <- range(nox)
# значения age, для которых делаем прогноз (от min до max с шагом 1)
nox.grid <- seq(noxlims[1], noxlims[2], length = 200) 
# рассчитать прогнозы и их стандартные ошибки
preds <- predict(fit, newdata = list(nox = nox.grid), se = T)
# границы доверительного интервала
se.bands <- cbind(preds$fit + 2*preds$se.fit,
                  preds$fit - 2*preds$se.fit)
par(mfrow = c(1, 2), mar = c(4.5, 4.5, 1, 1), oma = c(0, 0, 4, 0))
# наблюдения
plot(nox, crim, xlim = noxlims, 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)
## integer(0)

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

# прогнозы по второму вызову модели
preds2 <- predict(fit.2, newdata = list(nox = nox.grid), se = T)
# максимальное расхождение между прогнозами по двум вариантам вызова модели
max(abs(preds$fit - preds2$fit))
## [1] 1.910028e-12

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

fit.1 <- lm(crim ~ nox, data = Boston)
fit.2 <- lm(crim ~ poly(nox, 2), data = Boston)
fit.3 <- lm(crim ~ poly(nox, 3), data = Boston)
fit.4 <- lm(crim ~ poly(nox, 4), data = Boston)
fit.5 <- lm(crim ~ poly(nox, 5), data = Boston)
round(anova(fit.1, fit.2, fit.3, fit.4, fit.5), 2)

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

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

noxlims <- range(nox)
nox.grid <- seq(noxlims[1], noxlims[2], length = 200)
fit <- loess(crim ~ nox, span = 0.2, data = Boston)
fit2 <- loess(crim ~ nox, span = 0.5, data = Boston)
plot(nox, crim, xlim = noxlims, cex = 0.5, col = 'darkgrey')+
title('Локальная регрессия')+
lines(nox.grid, predict(fit, data.frame(nox = nox.grid)),col = 'red', lwd = 2)+
lines(nox.grid, predict(fit2, data.frame(nox = nox.grid)),col = 'blue', lwd = 2)   
## integer(0)
legend('topright',c('s = 0.2', 's = 0.5'))

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

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

# GAM на натуральных сплайнах
gam.ns <- gam(crim ~ ns(indus, 4) + ns(nox, 5) + chas, data = Boston)

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

# GAM на сглаживающих сплайнах
gam.m3 <- gam(crim ~ s(indus, 4) + s(nox, 5) + chas, data = Boston)
par(mfrow = c(1, 3))
plot(gam.m3, se = T, col = 'blue')

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

par(mfrow = c(1, 1))

Сделаем ANOVA, чтобы понять, какая степень для chas лучше.

gam.m1 <- gam(crim ~ s(nox, 5) + chas, data = Boston)          
gam.m2 <- gam(crim ~ indus + s(nox, 5) + chas, data = Boston)   
anova(gam.m1, gam.m2, gam.m3, test = 'F')

Третья модель статистически лучше второй.

# сводка по модели gam.m3
summary(gam.m3)
## 
## Call: gam(formula = crim ~ s(indus, 4) + s(nox, 5) + chas, data = Boston)
## Deviance Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.4891 -1.9341  0.1456  0.3604 75.1531 
## 
## (Dispersion Parameter for gaussian family taken to be 45.027)
## 
##     Null Deviance: 37363.22 on 505 degrees of freedom
## Residual Deviance: 22288.37 on 495.0002 degrees of freedom
## AIC: 3375.319 
## 
## Number of Local Scoring Iterations: 2 
## 
## Anova for Parametric Effects
##              Df  Sum Sq Mean Sq  F value    Pr(>F)    
## s(indus, 4)   1  5144.7  5144.7 114.2576 < 2.2e-16 ***
## s(nox, 5)     1   604.5   604.5  13.4251 0.0002751 ***
## chas          1   118.1   118.1   2.6234 0.1059351    
## Residuals   495 22288.4    45.0                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Anova for Nonparametric Effects
##             Npar Df Npar F     Pr(F)    
## (Intercept)                             
## s(indus, 4)       3 17.358 9.903e-11 ***
## s(nox, 5)         4 21.311 3.331e-16 ***
## chas                                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

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

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

# GAM на локальных регрессиях
gam.lo <- gam(crim ~ s(indus, df = 4) + lo(nox, span = 0.7) + chas, 
              data = Boston)
par(mfrow = c(1, 3))
plot(gam.lo, se = T, col = 'green')

par(mfrow = c(1, 1))
# модель со взаимодействием регрессоров indus и nox
gam.lo.i <- gam(crim ~ lo(indus, nox, span = 0.5) + chas, data = Boston)
par(mfrow = c(1, 2))
plot(gam.lo.i)

par(mfrow = c(1, 1))
detach(Boston)