Оценим стандартную ошибку модели для линейных регрессионных моделей а) со всеми объясняющими переменными; б) только с непрерывными объясняющими переменными Будем использовать методы:
методом проверочной выборки с долей обучающей 50%; методом LOOCV; k-кратной кросс-валидацией с k=5 и k=10 Загрузим данные и проведем все необходимые вычисления:
library(MASS)
data(Boston)
model1 <- lm(crim ~ indus + age + indus:chas + age:chas,
data = Boston)
summary(model1)##
## Call:
## lm(formula = crim ~ indus + age + indus:chas + age:chas, data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.798 -2.866 -0.598 1.087 81.222
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.99849 0.91212 -4.384 1.42e-05 ***
## indus 0.39465 0.06754 5.843 9.22e-09 ***
## age 0.05015 0.01627 3.082 0.00217 **
## indus:chas -0.15806 0.31836 -0.496 0.61977
## age:chas -0.01554 0.05544 -0.280 0.77931
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.775 on 501 degrees of freedom
## Multiple R-squared: 0.1894, Adjusted R-squared: 0.183
## F-statistic: 29.27 on 4 and 501 DF, p-value: < 2.2e-16model2 <- lm(crim ~ indus + age + indus:chas,
data = Boston)
summary(model2)##
## Call:
## lm(formula = crim ~ indus + age + indus:chas, data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.832 -2.848 -0.602 1.084 81.220
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.00331 0.91112 -4.394 1.36e-05 ***
## indus 0.39838 0.06616 6.022 3.34e-09 ***
## age 0.04949 0.01609 3.076 0.00221 **
## indus:chas -0.24299 0.09782 -2.484 0.01331 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.768 on 502 degrees of freedom
## Multiple R-squared: 0.1893, Adjusted R-squared: 0.1845
## F-statistic: 39.08 on 3 and 502 DF, p-value: < 2.2e-16#общее число наблюдений
m <- nrow(Boston)
# доля обучающей выборки
train.percent <- 0.5
# выбрать наблюдения в обучающую выборку
inTrain <- sample(m, m * train.percent)
inTrain## [1] 260 105 340 351 460 321 211 120 320 358 421 365 116 56 230 385 152
## [18] 199 145 143 136 188 190 296 503 295 307 115 185 79 488 294 231 240
## [35] 279 158 432 389 80 147 450 350 99 248 44 416 20 462 400 370 444
## [52] 8 371 221 43 467 395 191 87 454 472 213 179 392 388 19 61 285
## [69] 97 407 71 149 2 355 438 24 409 363 499 173 268 272 463 23 204
## [86] 360 202 119 429 411 487 84 357 393 446 457 70 34 76 428 368 102
## [103] 14 171 282 238 18 280 154 243 469 311 451 401 491 439 433 445 352
## [120] 384 100 237 9 121 324 124 142 46 48 31 144 313 41 210 11 156
## [137] 3 110 465 484 453 497 1 394 259 455 299 234 281 339 269 346 375
## [154] 246 67 386 198 319 94 335 53 65 315 424 175 62 264 274 332 49
## [171] 500 498 255 180 382 309 336 390 222 220 205 101 196 325 4 177 440
## [188] 442 329 414 232 16 397 98 502 405 75 262 112 104 423 95 38 174
## [205] 477 108 113 431 349 459 45 458 425 217 33 157 447 373 195 316 68
## [222] 241 417 448 165 310 235 475 347 148 208 129 441 160 367 489 245 224
## [239] 138 109 378 59 93 306 478 300 141 123 293 369 81 212 139# Линейная модель(a) ##############################################################
# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(Boston)
# подгонка линейной модели на обучающей выборке
fit.lm.1 <- lm(crim ~ indus + age + indus:chas,
subset = inTrain)
# считаем MSE на тестовой выборке
mean.a.1 <- mean((crim[-inTrain] - predict(fit.lm.1,
Boston[-inTrain, ]))^2)
# отсоединить таблицу с данными
detach(Boston)
# Квадратичная модель ##########################################################
# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(Boston)
# подгонка линейной модели на обучающей выборке
fit.lm.2 <- lm(crim ~ poly(age,2) + poly(indus, 2) + indus:chas,
subset = inTrain)
summary(fit.lm.2)##
## Call:
## lm(formula = crim ~ poly(age, 2) + poly(indus, 2) + indus:chas,
## subset = inTrain)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.226 -3.246 -0.605 1.619 42.692
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.5391 0.3819 9.268 < 2e-16 ***
## poly(age, 2)1 23.3584 11.3953 2.050 0.0414 *
## poly(age, 2)2 15.9602 9.3043 1.715 0.0875 .
## poly(indus, 2)1 59.7007 11.4077 5.233 3.56e-07 ***
## poly(indus, 2)2 -17.8190 8.7368 -2.040 0.0425 *
## indus:chas -0.1880 0.1004 -1.872 0.0624 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.848 on 247 degrees of freedom
## Multiple R-squared: 0.283, Adjusted R-squared: 0.2684
## F-statistic: 19.49 on 5 and 247 DF, p-value: 2.35e-16fit.lm.2.1 <- lm(crim ~ poly(age,2) + poly(indus, 2),
subset = inTrain)
summary(fit.lm.2.1)##
## Call:
## lm(formula = crim ~ poly(age, 2) + poly(indus, 2), subset = inTrain)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.984 -3.156 -0.558 1.622 43.016
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.3671 0.3725 9.039 < 2e-16 ***
## poly(age, 2)1 21.9676 11.4284 1.922 0.0557 .
## poly(age, 2)2 16.0476 9.3511 1.716 0.0874 .
## poly(indus, 2)1 58.2974 11.4404 5.096 6.89e-07 ***
## poly(indus, 2)2 -16.1585 8.7355 -1.850 0.0655 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.878 on 248 degrees of freedom
## Multiple R-squared: 0.2728, Adjusted R-squared: 0.2611
## F-statistic: 23.26 on 4 and 248 DF, p-value: 2.445e-16# считаем MSE на тестовой выборке
mean.a.2 <- mean((crim[-inTrain] - predict(fit.lm.2.1,
Boston[-inTrain, ]))^2)
# отсоединить таблицу с данными
detach(Boston)
# Кубическая модель ############################################################
# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(Boston)
# подгонка линейной модели на обучающей выборке
fit.lm.3 <- lm(crim ~ poly(age,3) + poly(indus, 3) + indus:chas,
subset = inTrain)
summary(fit.lm.3)##
## Call:
## lm(formula = crim ~ poly(age, 3) + poly(indus, 3) + indus:chas,
## subset = inTrain)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.544 -1.903 -0.326 1.024 41.179
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.45196 0.35492 9.726 < 2e-16 ***
## poly(age, 3)1 23.84287 10.58623 2.252 0.02519 *
## poly(age, 3)2 8.92303 8.76261 1.018 0.30954
## poly(age, 3)3 4.96348 8.20359 0.605 0.54571
## poly(indus, 3)1 59.11759 10.59758 5.578 6.41e-08 ***
## poly(indus, 3)2 -22.01279 8.32821 -2.643 0.00874 **
## poly(indus, 3)3 -50.50186 8.15377 -6.194 2.46e-09 ***
## indus:chas -0.21440 0.09334 -2.297 0.02246 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.431 on 245 degrees of freedom
## Multiple R-squared: 0.3866, Adjusted R-squared: 0.3691
## F-statistic: 22.06 on 7 and 245 DF, p-value: < 2.2e-16# считаем MSE на тестовой выборке
mean.a.3 <- mean((crim[-inTrain] - predict(fit.lm.3,
Boston[-inTrain, ]))^2)
# отсоединить таблицу с данными
detach(Boston)
# k-кратная перекрёстная проверка ==============================================
# оценим точность полиномиальных моделей, меняя степень
# вектор с ошибками по 10-кратной кросс-валидации
cv.err.k.fold <- rep(0, 5)
names(cv.err.k.fold) <- 1:5
# цикл по степеням полиномов
for (i in 1:5){
fit.glm <- glm(crim ~ poly(age,3) + poly(indus, 3) + indus: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
## 52.64665 53.25052 52.83313 53.01211 53.11844Наименьшая ошибка у кубической модели. Она равна 54.13. У к-кратной выборки наименьшая ошибка у i=3.
Проведем аналогичную работу для нашей модели, исключив из нее фактор, оставим только дискретные величины.
# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(Boston)
# подгонка линейной модели на обучающей выборке
fit.lm.1b <- lm(crim ~ indus + age,
subset = inTrain)
summary(fit.lm.1b)##
## Call:
## lm(formula = crim ~ indus + age, subset = inTrain)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.055 -2.767 -0.621 1.173 43.777
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.67566 1.06496 -3.451 0.000654 ***
## indus 0.41231 0.07134 5.780 2.22e-08 ***
## age 0.03572 0.01801 1.984 0.048400 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.923 on 250 degrees of freedom
## Multiple R-squared: 0.2555, Adjusted R-squared: 0.2495
## F-statistic: 42.9 on 2 and 250 DF, p-value: < 2.2e-16# считаем MSE на тестовой выборке
mean.b.1 <- mean((crim[-inTrain] - predict(fit.lm.1b,
Boston[-inTrain, ]))^2)
# отсоединить таблицу с данными
detach(Boston)
# Квадратичная модель ##########################################################
# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(Boston)
# подгонка линейной модели на обучающей выборке
fit.lm.2b <- lm(crim ~ poly(age,2) + poly(indus, 2),
subset = inTrain)
summary(fit.lm.2b)##
## Call:
## lm(formula = crim ~ poly(age, 2) + poly(indus, 2), subset = inTrain)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.984 -3.156 -0.558 1.622 43.016
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.3671 0.3725 9.039 < 2e-16 ***
## poly(age, 2)1 21.9676 11.4284 1.922 0.0557 .
## poly(age, 2)2 16.0476 9.3511 1.716 0.0874 .
## poly(indus, 2)1 58.2974 11.4404 5.096 6.89e-07 ***
## poly(indus, 2)2 -16.1585 8.7355 -1.850 0.0655 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.878 on 248 degrees of freedom
## Multiple R-squared: 0.2728, Adjusted R-squared: 0.2611
## F-statistic: 23.26 on 4 and 248 DF, p-value: 2.445e-16# считаем MSE на тестовой выборке
mean.b.2 <- mean((crim[-inTrain] - predict(fit.lm.2b,
Boston[-inTrain, ]))^2)
# отсоединить таблицу с данными
detach(Boston)
# Кубическая модель ############################################################
# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(Boston)
# подгонка линейной модели на обучающей выборке
fit.lm.3b <- lm(crim ~ poly(age,3) + poly(indus, 3),
subset = inTrain)
summary(fit.lm.3b)##
## Call:
## lm(formula = crim ~ poly(age, 3) + poly(indus, 3), subset = inTrain)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.189 -1.932 -0.237 1.068 41.564
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.2574 0.3477 9.370 < 2e-16 ***
## poly(age, 3)1 22.2582 10.6551 2.089 0.0377 *
## poly(age, 3)2 9.1111 8.8381 1.031 0.3036
## poly(age, 3)3 5.0719 8.2745 0.613 0.5405
## poly(indus, 3)1 57.5353 10.6667 5.394 1.62e-07 ***
## poly(indus, 3)2 -20.0971 8.3581 -2.405 0.0169 *
## poly(indus, 3)3 -49.6581 8.2160 -6.044 5.53e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.478 on 246 degrees of freedom
## Multiple R-squared: 0.3734, Adjusted R-squared: 0.3581
## F-statistic: 24.43 on 6 and 246 DF, p-value: < 2.2e-16# считаем MSE на тестовой выборке
mean.b.3 <- mean((crim[-inTrain] - predict(fit.lm.3b,
Boston[-inTrain, ]))^2)
# отсоединить таблицу с данными
detach(Boston)
# k-кратная перекрёстная проверка ==============================================
# оценим точность полиномиальных моделей, меняя степень
# вектор с ошибками по 10-кратной кросс-валидации
cv.err.k.fold.b <- rep(0, 5)
names(cv.err.k.fold.b) <- 1:5
# цикл по степеням полиномов
for (i in 1:5){
fit.glm.b <- glm(crim ~ poly(age,3) + poly(indus, 3), data = Boston)
cv.err.k.fold.b[i] <- cv.glm(Boston, fit.glm,
K = 10)$delta[1]
}
# результат
cv.err.k.fold.b## 1 2 3 4 5
## 52.75442 52.59038 53.19932 52.72476 53.05409в этом случае наименьшая оштбка также у кубической модели. Она равна 54,52. У к-кратной выборки наименьшая ошибка у i=4.
Применим бутстреп к нашей модели, сравним ее показатели с оценками, полученными с помощью МНК.
attach(Boston)
boot.fn <- function(data, index){
coef(lm(crim ~ poly(age,3) + poly(indus, 3) + indus:chas, subset = index))
}
boot.fn(Boston, 1:506)## (Intercept) poly(age, 3)1 poly(age, 3)2 poly(age, 3)3
## 3.8738819 30.1766409 20.2827291 15.6602755
## poly(indus, 3)1 poly(indus, 3)2 poly(indus, 3)3 indus:chas
## 57.2257259 -28.2091081 -49.6572388 -0.2959349# применяем функцию 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.8738819 -0.0104830183 0.35556340
## t2* 30.1766409 -0.0302793177 6.71623769
## t3* 20.2827291 -0.0780726953 6.58510902
## t4* 15.6602755 0.1160958063 7.20312225
## t5* 57.2257259 -0.4994779639 7.33903987
## t6* -28.2091081 -0.9736492584 5.77745850
## t7* -49.6572388 -0.8332395312 6.44087038
## t8* -0.2959349 0.0007056421 0.05144689# сравним с МНК
attach(Boston)
summary(lm(crim ~ poly(age,3) + poly(indus, 3) + indus:chas))$coef## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.8738819 0.33243436 11.653073 6.449355e-28
## poly(age, 3)1 30.1766409 9.81390385 3.074887 2.221466e-03
## poly(age, 3)2 20.2827291 7.76149057 2.613252 9.239761e-03
## poly(age, 3)3 15.6602755 7.48825414 2.091312 3.700621e-02
## poly(indus, 3)1 57.2257259 10.16909934 5.627413 3.056997e-08
## poly(indus, 3)2 -28.2091081 7.42224589 -3.800616 1.622093e-04
## poly(indus, 3)3 -49.6572388 7.44514574 -6.669747 6.829902e-11
## indus:chas -0.2959349 0.09188861 -3.220583 1.362901e-03Оценки получились в обоих случаях практически одинаковыми, различаются они лишь ошибками, так что нельзя сказать какая модель лучше.