R Notebook
Математическое моделирование
Кросс-валидация и бутстреп
В данной задаче выполняются следующие пункты:
- оценка точности модели методом перекрёстной выборки;
- методом проверочной выборки;
- методом перекрёстной проверки по отдельным наблюдениям (LOOCV);
- методом k-кратной перекрёстной проверки;
- применение бутстреп для оценки точности статистического параметра и оценок параметров модели
Модели: линейная регрессия, kNN.
Данные: статистика стоимости жилья в пригороде Бостона.
library('GGally') # графики совместного разброса переменных## Warning: package 'GGally' was built under R version 3.5.3## Loading required package: ggplot2library('lmtest') # тесты остатков регрессионных моделей## Warning: package 'lmtest' was built under R version 3.5.3## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary('FNN') # алгоритм kNN## Warning: package 'FNN' was built under R version 3.5.3library('mlbench')## Warning: package 'mlbench' was built under R version 3.5.3library('MASS')## Warning: package 'MASS' was built under R version 3.5.3library('ISLR')## Warning: package 'ISLR' was built under R version 3.5.3library('boot')
my.seed <- 1Метод перекрёстной проверки
# Пример на данных по районам: Boston
data(Boston) # открываем данные
?Boston## starting httpd help server ... doneBoston <- Boston[,-c(2,5,6,8,9,10,11,12,13,14)]
Boston$chas <- as.factor(Boston$chas)
# графики разброса
ggpairs(Boston)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Метод проверочной выборки
Он состоит в том, что мы отбираем одну тестовую выборку и будем считать на ней ошибку модели.
n <- nrow(Boston)
# доля обучающей выборки
train.percent <- 0.5
# выбрать наблюдения в обучающую выборку
set.seed(my.seed)
inTrain <- sample(n, n*train.percent)
inTrain## [1] 135 188 289 457 102 451 473 330 314 31 103 88 340 190 379 245 352
## [18] 486 186 492 455 496 316 61 129 488 7 184 416 163 230 285 234 89
## [35] 391 315 374 51 339 193 383 301 364 257 491 464 11 220 336 317 218
## [52] 392 199 111 32 45 143 233 297 182 408 131 204 148 288 114 211 337
## [69] 37 466 443 366 151 145 206 385 372 168 333 411 481 303 170 138 320
## [86] 86 299 469 485 60 100 25 266 362 321 328 187 429 331 247 414 144
## [103] 110 401 255 421 52 191 368 238 387 290 141 424 59 6 280 41 174
## [120] 248 437 399 426 67 504 173 195 79 87 225 217 29 14 240 346 222
## [137] 208 380 363 410 250 459 378 94 265 164 64 269 38 309 219 198 117
## [154] 160 177 360 384 27 97 74 99 357 154 268 302 394 22 441 462 468
## [171] 212 282 286 445 127 448 213 244 454 296 96 63 442 477 283 483 243
## [188] 232 465 388 377 123 452 291 500 479 35 260 98 241 82 386 158 306
## [205] 55 157 169 39 77 214 475 30 484 278 373 90 189 276 275 308 76
## [222] 48 92 433 261 405 73 13 354 237 326 438 418 381 107 367 271 150
## [239] 116 121 425 344 307 358 109 56 112 310 119 422 196 397 120Построим модели для проверки точности. Вид моделей:
\[
\hat{crim} = f(age,indus, chas)
\] Линейная модель: \(\hat{crim} = \hat{\beta}_0 + \hat{\beta}_1 \cdot age + \hat{\beta}_2 \cdot indus + \hat{\beta}_3 \cdot chas\).
# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(Boston)
# подгонка линейной модели на обучающей выборке
fit.lm.1 <- lm(crim ~ age + indus + chas, subset = inTrain)
# считаем MSE на тестовой выборке
mean((crim[-inTrain] - predict(fit.lm.1, Boston[-inTrain,]))^2)## [1] 58.21777# отсоединить таблицу с данными
detach(Boston)Строим первую квадратичную модель: \(\hat{crim} = \hat{\beta}_0 + \hat{\beta}_1 \cdot age^2 + \hat{\beta}_2 \cdot indus+ \hat{\beta}_3 \cdot chas\).
# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(Boston)
# подгонка линейной модели на обучающей выборке
fit.lm.2 <- lm(crim ~ poly(age,2)+ indus + chas ,subset = inTrain)
# считаем MSE на тестовой выборке
mean((crim[-inTrain] - predict(fit.lm.2, Boston[-inTrain,]))^2)## [1] 56.99055# отсоединить таблицу с данными
detach(Boston)Строим вторую квадратичную модель: \(\hat{crim} = \hat{\beta}_0 + \hat{\beta}_1 \cdot indus^2 + \hat{\beta}_2 \cdot age+ \hat{\beta}_3 \cdot chas\).
# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(Boston)
# подгонка линейной модели на обучающей выборке
fit.lm.3 <- lm(crim ~ age + poly(indus,2) + chas ,subset = inTrain)
# считаем MSE на тестовой выборке
mean((crim[-inTrain] - predict(fit.lm.3, Boston[-inTrain,]))^2)## [1] 56.69466# отсоединить таблицу с данными
detach(Boston)Строим первую кубическую модель: \(\hat{crim} = \hat{\beta}_0 + \hat{\beta}_1 \cdot age^3 + \hat{\beta}_2 \cdot indus + \hat{\beta}_3 \cdot chas\).
# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(Boston)
# подгонка линейной модели на обучающей выборке
fit.lm.5 <- lm(crim ~ poly(age,3)+ indus + chas ,subset = inTrain)
# считаем MSE на тестовой выборке
mean((crim[-inTrain] - predict(fit.lm.5, Boston[-inTrain,]))^2)## [1] 55.85914# отсоединить таблицу с данными
detach(Boston)Строим вторую кубическую модель: \(\hat{crim} = \hat{\beta}_0 + \hat{\beta}_1 \cdot indus^3 + \hat{\beta}_2 \cdot age + \hat{\beta}_3 \cdot chas\).
# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(Boston)
# подгонка линейной модели на обучающей выборке
fit.lm.4 <- lm(crim ~ age + poly(indus,3) + chas ,subset = inTrain)
# считаем MSE на тестовой выборке
mean((crim[-inTrain] - predict(fit.lm.4, Boston[-inTrain,]))^2)## [1] 50.83844# отсоединить таблицу с данными
detach(Boston)Ошибка последней кубической модели оказалась наименьшей из всех построенных, следовательно она и будет являться наиболее пригодной для прогнозирования.
Перекрёстная проверка по отдельным наблюдениям (LOOCV)
Это самый затратный в вычислительном плане метод, но и самый надёжный в плане оценки ошибки вне выборки. Попробуем применить его к линейной модели.
# подгонка линейной модели на обучающей выборке
fit.glm <- glm(crim ~ age + indus + chas, data = Boston)
# считаем LOOCV-ошибку
cv.err <- cv.glm(Boston, fit.glm)
# результат: первое число -- по формуле LOOCV-ошибки,
# второе -- с поправкой на смещение
cv.err$delta[1]## [1] 60.68151Теперь оценим точность полиномиальных моделей, меняя степень, в которой стоит регрессор.
Вид модели 1: crim ~ poly(age,i) + indus + chas
# вектор с LOOCV-ошибками
cv.err.loocv <- rep(0, 5)
names(cv.err.loocv) <- 1:5
# цикл по степеням полиномов
for (i in 1:5){
fit.glm <- glm(crim ~ poly(age,i) + indus + chas, data = Boston)
cv.err.loocv[i] <- cv.glm(Boston, fit.glm)$delta[1]
}
# результат
cv.err.loocv## 1 2 3 4 5
## 60.68151 59.81116 59.22108 59.10194 58.78232Вид модели 2: crim ~ age + poly(indus,i) + chas
# вектор с LOOCV-ошибками
cv.err.loocv <- rep(0, 5)
names(cv.err.loocv) <- 1:5
# цикл по степеням полиномов
for (i in 1:5){
fit.glm <- glm(crim ~ age + poly(indus,i) + chas, data = Boston)
cv.err.loocv[i] <- cv.glm(Boston, fit.glm)$delta[1]
}
# результат
cv.err.loocv## 1 2 3 4 5
## 60.68151 59.62857 54.14634 54.08104 50.54975k-кратная перекрёстная проверка
K-кратная кросс-валидация – компромисс между методом проверочной выборки и LOOCV. Оценка ошибки вне выборки ближе к правде, по сравнению с проверочной выборкой, а объём вычислений меньше, чем при LOOCV. Проведём 5-кратную и 10-кратную кросс-валидацию моделей разных степеней вида модели 2 (“crim ~ age + poly(indus,i) + chas”), так как предыдущий пункт показал, что ошибки в данном случае оказались ниже.
Проведём 5-кратную кросс-валидацию:
# оценим точность полиномиальных моделей, меняя степень
# вектор с ошибками по 5-кратной кросс-валидации
cv.err.k.fold <- rep(0, 5)
names(cv.err.k.fold) <- 1:5
# цикл по степеням полиномов
for (i in 1:5){
fit.glm <- glm(crim ~ age + poly(indus,i) + chas, data = Boston)
cv.err.k.fold[i] <- cv.glm(Boston, fit.glm,
K = 5)$delta[1]
}
# результат
cv.err.k.fold## 1 2 3 4 5
## 60.91788 59.34413 53.87487 54.89812 51.08626Проведём 10-кратную кросс-валидацию:
# оценим точность полиномиальных моделей, меняя степень
# вектор с ошибками по 10-кратной кросс-валидации
cv.err.k.fold <- rep(0, 10)
names(cv.err.k.fold) <- 1:10
# цикл по степеням полиномов
for (i in 1:10){
fit.glm <- glm(crim ~ age + poly(indus,i) + chas, data = Boston)
cv.err.k.fold[i] <- cv.glm(Boston, fit.glm,
K = 10)$delta[1]
}
# результат
cv.err.k.fold## 1 2 3 4 5 6 7 8
## 60.80360 59.74739 54.50540 54.08287 50.59756 50.32679 47.81233 47.90084
## 9 10
## 47.26915 46.40347Опираясь на результаты расчётов с кросс-валидацией, можно заключить, что на самом деле ошибка вне выборки у линейной модели выше, чем показывала MSE на тестовой выборке. Методы кросс-валидации сходятся на одной и той же модели: на первом этапе (первая степень) модель показывает наименьшую ошибку.
Бутстреп
Точность оценки параметра регрессии
При построении модели регрессии проблемы в остатках приводят к неверной оценке ошибок параметров. Обойти эту проблему можно, применив для расчёта этих ошибок бутстреп.
# Оценивание точности линейной регрессионной модели ----------------------------
# оценить стандартные ошибки параметров модели
# mpg = beta_0 + beta_1 * horsepower с помощью бутстрепа,
# сравнить с оценками ошибок по МНК
# функция для расчёта коэффициентов ПЛР по выборке из данных
boot.fn <- function(data, index){
coef(lm(crim ~ age + indus + chas, data = data, subset = index))
}
boot.fn(Boston, 1:n)## (Intercept) age indus chas1
## -3.86979127 0.04951804 0.38581635 -3.02343988# пример применения функции к бутстреп-выборке
set.seed(my.seed)
boot.fn(Boston, sample(n, n, replace = T))## (Intercept) age indus chas1
## -3.45370810 0.04865728 0.34109612 -3.51748129# применяем функцию boot для вычисления стандартных ошибок параметров
# (1000 выборок с повторами)
boot(Boston, boot.fn, 1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston, statistic = boot.fn, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* -3.86979127 -0.029351624 0.495122773
## t2* 0.04951804 0.000224052 0.008702608
## t3* 0.38581635 0.002985588 0.051688604
## t4* -3.02343988 -0.019640138 0.621956134# сравним с МНК
attach(Boston)
summary(lm(crim ~ age + indus + chas))$coef## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.86979127 0.91184502 -4.243913 2.616812e-05
## age 0.04951804 0.01611478 3.072834 2.235573e-03
## indus 0.38581635 0.06600382 5.845364 9.107782e-09
## chas1 -3.02343988 1.36781343 -2.210418 2.752638e-02detach(Boston)
#оценим наилучшую найденную модель
boot.fn.2 <- function(data, index){
coef(lm(crim ~ age + poly(indus,3) + chas, data = data, subset = index))
}
# применим функцию к 1000 бутсреп-выборкам
set.seed(my.seed)
boot(Boston, boot.fn, 1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston, statistic = boot.fn, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* -3.86979127 -0.0287147264 0.494620951
## t2* 0.04951804 0.0002217635 0.008703169
## t3* 0.38581635 0.0029202760 0.051693557
## t4* -3.02343988 -0.0194853173 0.622807107В модели регрессии, для которой проводился расчёт, похоже, не нарушаются требования к остаткам, и оценки стандартных ошибок параметров, рассчитанные по МНК, очень близки к ошибкам этих же параметров, полученных бутстрепом.