Home work 6
Romanova ER
9 апреля 2018 г
Упражнение 6
Регуляризация линейных моделей
В домашней работе необходимо:
- Применить отбор путём пошагового включения к набору данных Auto {ISLR}. Выбрать оптимальную модель с помощью кросс-валидации. Вывести её коэффициенты с помощью функции coef(). Рассчитать MSE модели на тестовой выборке.
2 Применить частный метод наименьших квадратов к набору данных Auto {ISLR}. Для модели:
Подогнать модель на всей выборке и вычислить ошибку (MSE) с кросс-валидацией. По наименьшей MSE подобрать оптимальное значение настроечного параметра метода.
Подогнать модель с оптимальным значением параметра на обучающей выборке, посчитать MSE на тестовой.
Подогнать модель с оптимальным значением параметра на всех данных, вывести характеристики модели функцией summary().
3 Сравнить оптимальные модели, полученные в заданиях 1 и 2 по MSE на тестовой выборке. Какой метод дал лучший результат? Доля тестовой выборки: 50%.
library('ISLR') # набор данных Auto## Warning: package 'ISLR' was built under R version 3.4.3library('leaps') # функция regsubset() -- отбор оптимального ## Warning: package 'leaps' was built under R version 3.4.4 # подмножества переменных
library('glmnet') # функция glmnet() -- лассо## Warning: package 'glmnet' was built under R version 3.4.4## Loading required package: Matrix## Loading required package: foreach## Warning: package 'foreach' was built under R version 3.4.4## Loaded glmnet 2.0-16library('pls') # регрессия на главные компоненты -- pcr()## Warning: package 'pls' was built under R version 3.4.4##
## Attaching package: 'pls'## The following object is masked from 'package:stats':
##
## loadings # и частный МНК -- plsr()
my.seed <- 1
str(Auto)## 'data.frame': 392 obs. of 9 variables:
## $ mpg : num 18 15 18 16 17 15 14 14 14 15 ...
## $ cylinders : num 8 8 8 8 8 8 8 8 8 8 ...
## $ displacement: num 307 350 318 304 302 429 454 440 455 390 ...
## $ horsepower : num 130 165 150 150 140 198 220 215 225 190 ...
## $ weight : num 3504 3693 3436 3433 3449 ...
## $ acceleration: num 12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
## $ year : num 70 70 70 70 70 70 70 70 70 70 ...
## $ origin : num 1 1 1 1 1 1 1 1 1 1 ...
## $ name : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...Auto <- Auto[,-9] #исключаем name
# считаем пропуски
sum(is.na(Auto$origin))## [1] 0Отбор путём пошагового включения переменных
#пошаговое включение
regfit.fwd <- regsubsets(mpg ~ ., data = Auto,
nvmax = 8, method = 'forward')
summary(regfit.fwd)## Subset selection object
## Call: regsubsets.formula(mpg ~ ., data = Auto, nvmax = 8, method = "forward")
## 7 Variables (and intercept)
## Forced in Forced out
## cylinders FALSE FALSE
## displacement FALSE FALSE
## horsepower FALSE FALSE
## weight FALSE FALSE
## acceleration FALSE FALSE
## year FALSE FALSE
## origin FALSE FALSE
## 1 subsets of each size up to 7
## Selection Algorithm: forward
## cylinders displacement horsepower weight acceleration year origin
## 1 ( 1 ) " " " " " " "*" " " " " " "
## 2 ( 1 ) " " " " " " "*" " " "*" " "
## 3 ( 1 ) " " " " " " "*" " " "*" "*"
## 4 ( 1 ) " " "*" " " "*" " " "*" "*"
## 5 ( 1 ) " " "*" "*" "*" " " "*" "*"
## 6 ( 1 ) "*" "*" "*" "*" " " "*" "*"
## 7 ( 1 ) "*" "*" "*" "*" "*" "*" "*"round(coef(regfit.fwd, 6), 3)## (Intercept) cylinders displacement horsepower weight
## -15.563 -0.507 0.019 -0.024 -0.006
## year origin
## 0.748 1.428Проверка кросс-валидацией
# функция для прогноза для функции regsubset()
predict.regsubsets <- function(object, newdata, id, ...){
form <- as.formula(object$call[[2]])
mat <- model.matrix(form, newdata)
coefi <- coef(object, id = id)
xvars <- names(coefi)
mat[, xvars] %*% coefi
}
# отбираем 10 блоков наблюдений
k <- 10
set.seed(my.seed)
folds <- sample(1:k, nrow(Auto), replace = T)
# заготовка под матрицу с ошибками
cv.errors <- matrix(NA, k, 7, dimnames = list(NULL, paste(1:7)))
# заполняем матрицу в цикле по блокам данных
for (j in 1:k){
best.fit <- regsubsets(mpg ~ ., data = Auto[folds != j, ],
nvmax = 7)
# теперь цикл по количеству объясняющих переменных
for (i in 1:7){
# модельные значения mpg
pred <- predict(best.fit, Auto[folds == j, ], id = i)
# вписываем ошибку в матрицу
cv.errors[j, i] <- mean((Auto$mpg[folds == j] - pred)^2)
}
}
# усредняем матрицу по каждому столбцу (т.е. по блокам наблюдений),
# чтобы получить оценку MSE для каждой модели с фиксированным
# количеством объясняющих переменных
mean.cv.errors <- apply(cv.errors, 2,
mean)
round(mean.cv.errors, 3)## 1 2 3 4 5 6 7
## 18.909 11.947 11.419 11.619 11.561 11.452 11.389# на графике
plot(mean.cv.errors, type = 'b')
points(which.min(mean.cv.errors), mean.cv.errors[which.min(mean.cv.errors)],
col = 'red', pch = 20, cex = 2)
K-кратная кросс-валидация показала, что наименьшая ошибка выходит, если все предикторы включены в модель
Регрессия по методу частных наименьших квадратов
set.seed(my.seed)
train <- sample(c(T, F), nrow(Auto), rep = T)
test <- !train
# из-за синтаксиса glmnet() формируем явно матрицу объясняющих...
x <- model.matrix(mpg ~ ., Auto)[,-1]
# и вектор значений зависимой переменной
y <- Auto$mpg
y.test <- y[test]
set.seed(my.seed)
pls.fit <- plsr(mpg ~ ., data = Auto, subset = train, scale = T,
validation = 'CV')
summary(pls.fit)## Data: X dimension: 212 7
## Y dimension: 212 1
## Fit method: kernelpls
## Number of components considered: 7
##
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
## (Intercept) 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps
## CV 7.634 3.679 3.316 3.199 3.126 3.116 3.119
## adjCV 7.634 3.676 3.308 3.179 3.123 3.109 3.108
## 7 comps
## CV 3.125
## adjCV 3.114
##
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps
## X 64.86 75.33 80.22 89.35 97.67 99.36 100.00
## mpg 77.46 82.56 84.08 84.73 84.98 85.18 85.21# теперь подгоняем модель для найденного оптимального M = 5
# и оцениваем MSE на тестовой
pls.pred <- predict(pls.fit, x[test, ], ncomp = 5)
round(mean((pls.pred - y.test)^2), 3)## [1] 14.275# подгоняем модель на всей выборке
pls.fit <- plsr(mpg ~ ., data = Auto, scale = T, ncomp = 5)
summary(pls.fit)## Data: X dimension: 392 7
## Y dimension: 392 1
## Fit method: kernelpls
## Number of components considered: 5
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps 4 comps 5 comps
## X 65.75 76.46 82.97 91.01 97.97
## mpg 74.44 79.88 80.93 81.60 81.86MSE на тестовой выборке методом пошагового включения является 11,389, в то время как MSE модели, полученной методом частных наименьших квадратов равна 14,275. Поэтому лучшей моделью является модель, построенная методом пошагового включения предикторов