laba6
Semin A.S.
29 марта 2018 г
РЕГУЛЯРИЗАЦИЯ ЛИНЕЙНЫХ МОДЕЛЕЙ
В практических примерах ниже показано:
отобрать предикторы методами пошагового включения
как построить ридж-регрессию
как применять эти методы в сочетании с кросс-валидацией
Модели: линейная регрессия, ридж.
Данные: Carseats {ISLR}
## [1] "Sales" "CompPrice" "Income" "Advertising" "Population"
## [6] "Price" "ShelveLoc" "Age" "Education" "Urban"
## [11] "US"Пошаговое включение
regfit.fwd <- regsubsets(Sales ~ ., data = Carseats,
nvmax = 10, method = 'forward')
summary(regfit.fwd)## Subset selection object
## Call: regsubsets.formula(Sales ~ ., data = Carseats, nvmax = 10, method = "forward")
## 11 Variables (and intercept)
## Forced in Forced out
## CompPrice FALSE FALSE
## Income FALSE FALSE
## Advertising FALSE FALSE
## Population FALSE FALSE
## Price FALSE FALSE
## ShelveLocGood FALSE FALSE
## ShelveLocMedium FALSE FALSE
## Age FALSE FALSE
## Education FALSE FALSE
## UrbanYes FALSE FALSE
## USYes FALSE FALSE
## 1 subsets of each size up to 10
## Selection Algorithm: forward
## CompPrice Income Advertising Population Price ShelveLocGood
## 1 ( 1 ) " " " " " " " " " " "*"
## 2 ( 1 ) " " " " " " " " "*" "*"
## 3 ( 1 ) "*" " " " " " " "*" "*"
## 4 ( 1 ) "*" " " "*" " " "*" "*"
## 5 ( 1 ) "*" " " "*" " " "*" "*"
## 6 ( 1 ) "*" " " "*" " " "*" "*"
## 7 ( 1 ) "*" "*" "*" " " "*" "*"
## 8 ( 1 ) "*" "*" "*" " " "*" "*"
## 9 ( 1 ) "*" "*" "*" " " "*" "*"
## 10 ( 1 ) "*" "*" "*" " " "*" "*"
## ShelveLocMedium Age Education UrbanYes USYes
## 1 ( 1 ) " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " "
## 3 ( 1 ) " " " " " " " " " "
## 4 ( 1 ) " " " " " " " " " "
## 5 ( 1 ) "*" " " " " " " " "
## 6 ( 1 ) "*" "*" " " " " " "
## 7 ( 1 ) "*" "*" " " " " " "
## 8 ( 1 ) "*" "*" " " " " "*"
## 9 ( 1 ) "*" "*" "*" " " "*"
## 10 ( 1 ) "*" "*" "*" "*" "*"round(coef(regfit.fwd, 10), 3)## (Intercept) CompPrice Income Advertising
## 5.762 0.093 0.016 0.125
## Price ShelveLocGood ShelveLocMedium Age
## -0.095 4.847 1.952 -0.046
## Education UrbanYes USYes
## -0.022 0.119 -0.199k-кратная кросс-валидация
# отбираем 10 блоков наблюдений
k <- 10
set.seed(my.seed)
folds <- sample(1:k, nrow(Carseats), replace = T)
# заготовка под матрицу с ошибками
cv.errors <- matrix(NA, k, 10, dimnames = list(NULL, paste(1:10)))
predict.regsubsets = function(object, newdata, id, ...) {
form = as.formula(object$call[[2]])
mat = model.matrix(form, newdata)
coefi = coef(object, id = id)
mat[, names(coefi)] %*% coefi
}
# заполняем матрицу в цикле по блокам данных
for (j in 1:k){
best.fit <- regsubsets(Sales ~ ., data = Carseats[folds != j, ],
nvmax = 10)
# теперь цикл по количеству объясняющих переменных
for (i in 1:10){
# модельные значения Sales
pred <- predict(best.fit, Carseats[folds == j, ], id = i)
# вписываем ошибку в матрицу
cv.errors[j, i] <- mean((Carseats$Sales[folds == j] - pred)^2)
}
}
# усредняем матрицу по каждому столбцу (т.е. по блокам наблюдений),
# чтобы получить оценку MSE для каждой модели с фиксированным
# количеством объясняющих переменных
mean.cv.errors <- apply(cv.errors, 2, mean)
round(mean.cv.errors, 2)## 1 2 3 4 5 6 7 8 9 10
## 6.12 4.28 3.03 2.48 2.00 1.27 1.09 1.12 1.12 1.11# на графике
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)
# перестраиваем модель с 7 объясняющими переменными на всём наборе данных
reg.best <- regsubsets(Sales ~ ., data = Carseats, nvmax = 10)
round(coef(reg.best, 7), 3)## (Intercept) CompPrice Income Advertising
## 5.475 0.093 0.016 0.116
## Price ShelveLocGood ShelveLocMedium Age
## -0.095 4.836 1.952 -0.046K-кратная кросс-валидация показала, что наименьшая ошибка выходит, если предикторов 7 штук, поэтому на модели, не считая константы, мы наблюдаем 7 регрессов вместо 10.
MSE модели на тестовой выборке
## [1] 5.88 4.15 3.05 2.22 1.61 1.13 0.95 0.94 0.94 0.95Гребневая регрессия
# из-за синтаксиса glmnet() формируем явно матрицу объясняющих...
x <- model.matrix(Sales ~ ., Carseats)[, -1]
# и вектор значений зависимой переменной
y <- Carseats$Sales
# вектор значений гиперпараметра лямбда
grid <- 10^seq(10, -2, length = 100)
# подгоняем серию моделей ридж-регрессии
ridge.mod <- glmnet(x, y, alpha = 0, lambda = grid)
# размерность матрицы коэффициентов моделей
dim(coef(ridge.mod))## [1] 12 100# значение лямбда под номером 50
round(ridge.mod$lambda[50], 2)## [1] 11497.57# коэффициенты соответствующей модели
round(coef(ridge.mod)[, 50], 3)## (Intercept) CompPrice Income Advertising
## 7.497 0.000 0.000 0.000
## Population Price ShelveLocGood ShelveLocMedium
## 0.000 0.000 0.001 0.000
## Age Education UrbanYes USYes
## 0.000 0.000 0.000 0.000# норма эль-два
round(sqrt(sum(coef(ridge.mod)[-1, 50]^2)), 2)## [1] 0# всё то же для лямбды под номером 60
# значение лямбда под номером 60
round(ridge.mod$lambda[60], 2)## [1] 705.48# коэффициенты соответствующей модели
round(coef(ridge.mod)[, 60], 3)## (Intercept) CompPrice Income Advertising
## 7.514 0.000 0.000 0.000
## Population Price ShelveLocGood ShelveLocMedium
## 0.000 0.000 0.014 -0.002
## Age Education UrbanYes USYes
## 0.000 0.000 0.000 0.004# норма эль-два
round(sqrt(sum(coef(ridge.mod)[-1, 60]^2)), 2)## [1] 0.01# мы можем получить значения коэффициентов для новой лямбды
round(predict(ridge.mod, s = 50, type = 'coefficients')[1:12, ], 3)## (Intercept) CompPrice Income Advertising
## 7.725 0.001 0.001 0.006
## Population Price ShelveLocGood ShelveLocMedium
## 0.000 -0.003 0.187 -0.017
## Age Education UrbanYes USYes
## -0.002 -0.003 -0.005 0.053set.seed(my.seed)
train <- sample(1:nrow(x), nrow(x)/2)
test <- -train
y.test <- y[test]
# подгоняем ридж-модели с большей точностью (thresh ниже значения по умолчанию)
ridge.mod <- glmnet(x[train, ], y[train], alpha = 0, lambda = grid,
thresh = 1e-12)
plot(ridge.mod)
# прогнозы для модели с лямбда = 4
ridge.pred <- predict(ridge.mod, s = 4, newx = x[test, ])
round(mean((ridge.pred - y.test)^2), 2)## [1] 4.95# сравним с MSE для нулевой модели (прогноз = среднее)
round(mean((mean(y[train]) - y.test)^2), 2)## [1] 8.33# насколько модель с лямбда = 4 отличается от обычной ПЛР
ridge.pred <- predict(ridge.mod, s = 0, newx = x[test, ], exact = T,
x = x[train, ], y = y[train])
round(mean((ridge.pred - y.test)^2), 2)## [1] 1.04# predict с лямбдой (s) = 0 даёт модель ПЛР
lm(y ~ x, subset = train)##
## Call:
## lm(formula = y ~ x, subset = train)
##
## Coefficients:
## (Intercept) xCompPrice xIncome xAdvertising
## 6.4680663 0.0877798 0.0144526 0.1047663
## xPopulation xPrice xShelveLocGood xShelveLocMedium
## 0.0009013 -0.0969584 4.8680820 2.0575944
## xAge xEducation xUrbanYes xUSYes
## -0.0485719 -0.0320330 0.3270469 -0.0180527round(predict(ridge.mod, s = 0, exact = T, type = 'coefficients',
x = x[train, ], y = y[train])[1:12, ], 3)## (Intercept) CompPrice Income Advertising
## 6.468 0.088 0.014 0.105
## Population Price ShelveLocGood ShelveLocMedium
## 0.001 -0.097 4.868 2.058
## Age Education UrbanYes USYes
## -0.049 -0.032 0.327 -0.018# k-кратная кросс-валидация
set.seed(my.seed)
# оценка ошибки
cv.out <- cv.glmnet(x[train, ], y[train], alpha = 0)
plot(cv.out)
# значение лямбда, обеспечивающее минимальную ошибку перекрёстной проверки
bestlam <- cv.out$lambda.min
round(bestlam, 2)## [1] 0.14# MSE на тестовой для этого значения лямбды
ridge.pred <- predict(ridge.mod, s = bestlam, newx = x[test, ])
val.lam <- mean((ridge.pred - y.test)^2)
round(val.lam, 2)## [1] 1.13# наконец, подгоняем модель для оптимальной лямбды,
# найденной по перекрёстной проверке
out <- glmnet(x, y, alpha = 0)
round(predict(out, type = 'coefficients', s = bestlam)[1:12, ], 3)## (Intercept) CompPrice Income Advertising
## 6.313 0.081 0.015 0.112
## Population Price ShelveLocGood ShelveLocMedium
## 0.000 -0.086 4.416 1.670
## Age Education UrbanYes USYes
## -0.043 -0.021 0.097 -0.074Подведем итог: MSE на тестовой выборке методом пошагового включения равняется 0.95, в то время как, метод гребневой регрессии дал нам результат 1.13.
Поэтому метод пошагового включения в этом случае является наиболее оптимальным.