Упражнение 6
Нестерова А.И.
28 03 2020
Математическое моделирование
Упражнение 6
1 Примените указанный в варианте метод к набору данных по своему варианту (см. таблицу ниже). Не забудьте предварительно сделать из категориальных переменных факторы. Выберите оптимальную модель с помощью кросс-валидации. Выведите её коэффициенты с помощью функции coef(). Рассчитайте MSE модели на тестовой выборке.
2 Примените указанный в варианте метод к набору данных по своему варианту (см. таблицу ниже). Для модели:
Подогнать модель на всей выборке и вычислить ошибку (MSE) с кросс-валидацией. По наименьшей MSE подобрать оптимальное значение настроечного параметра метода (гиперпараметр λ или число главных компонент M).
Подогнать модель с оптимальным значением параметра на обучающей выборке, посчитать MSE на тестовой.
Подогнать модель с оптимальным значением параметра на всех данных, вывести характеристики модели функцией summary().
3 Сравните оптимальные модели, полученные в заданиях 1 и 2 по MSE на тестовой выборке. Какой метод дал лучший результат? Доля тестовой выборки: 50%.
Как сдавать: прислать на почту преподавателя ссылки:
- на html-отчёт с видимыми блоками кода (блоки кода с параметром echo = T), размещённый на rpubs.com.
- на код, генерирующий отчёт, в репозитории на github.com.
В текст отчёта включить постановку задачи и ответы на вопросы задания. В начале отчёта нужно описать рассматриваемый набор данных: тип и содержание переменных, количество наблюдений.
Вариант - 13
Методы: отбор оптимального подмножества, гребневая регрессия. Данные: `Prestige {car}’.
Набор данных Prestige содержит переменные:
education- среднее образование профессиональных сотрудников в 1971 г. (в годах);income– средний доход должностных лиц в 1971 г. (в долларах);
women– процент работающих женщин;prestige– баллы престижа по профессии из социального опроса, проведённого в середине 1960-х годов;census- канадская перепись профессиональных кодов;type- род занятия. Фактор с уровнями:
bc – синие воротнички; prof – профессионалы, менеджеры и техники; wc - белые воротнички.
library('knitr') # пакет для генерации отчёта
library('car') # набор данных Prestige
library('leaps') # функция regsubset() -- отбор оптимального подмножества переменных
library('glmnet') # функция glmnet() -- лассо
my.seed <- 1
# загрузка данных Prestige
data('Prestige')
# переводим дискретные количественные переменные в факторы
Prestige$type <- as.factor(Prestige$type)Набор данных по престижу канадских профессий Prestige.
# название столбцов переменных
names(Prestige)## [1] "education" "income" "women" "prestige" "census" "type"# размерность данных
dim(Prestige)## [1] 102 6Считаем число пропусков в данных и убираем их.
# считаем пропуски
sum(is.na(Prestige))## [1] 4# убираем пропуски
Prestige <- na.omit(Prestige)
# проверяем результат
dim(Prestige)## [1] 98 6sum(is.na(Prestige))## [1] 0Отбор оптимального подмножества
# подгоняем модели с сочетаниями предикторов до 6 (максимум в данных)
regfit.full <- regsubsets(prestige ~ ., Prestige)
reg.summary <- summary(regfit.full)
reg.summary## Subset selection object
## Call: regsubsets.formula(prestige ~ ., Prestige)
## 6 Variables (and intercept)
## Forced in Forced out
## education FALSE FALSE
## income FALSE FALSE
## women FALSE FALSE
## census FALSE FALSE
## typeprof FALSE FALSE
## typewc FALSE FALSE
## 1 subsets of each size up to 6
## Selection Algorithm: exhaustive
## education income women census typeprof typewc
## 1 ( 1 ) "*" " " " " " " " " " "
## 2 ( 1 ) "*" "*" " " " " " " " "
## 3 ( 1 ) "*" "*" " " " " "*" " "
## 4 ( 1 ) "*" "*" " " "*" "*" " "
## 5 ( 1 ) "*" "*" "*" "*" "*" " "
## 6 ( 1 ) "*" "*" "*" "*" "*" "*"# структура отчёта по модели (ищем характеристики качества)
names(reg.summary)## [1] "which" "rsq" "rss" "adjr2" "cp" "bic" "outmat" "obj"# R^2 и скорректированный R^2
round(reg.summary$rsq, 3)## [1] 0.751 0.814 0.833 0.841 0.841 0.841# на графике
plot(1:6, reg.summary$rsq, type = 'b',
xlab = 'Количество предикторов', ylab = 'R-квадрат')
# сюда же добавим скорректированный R-квадрат
points(1:6, reg.summary$adjr2, col = 'red')
# модель с максимальным скорректированным R-квадратом
which.max(reg.summary$adjr2)## [1] 4### 4
points(which.max(reg.summary$adjr2),
reg.summary$adjr2[which.max(reg.summary$adjr2)],
col = 'red', cex = 2, pch = 20)
legend('bottomright', legend = c('R^2', 'R^2_adg'),
col = c('black', 'red'), lty = c(1, NA),
pch = c(1, 1))
# C_p
reg.summary$cp## [1] 48.672200 14.481794 5.747918 3.221668 5.008399 7.000000# число предикторов у оптимального значения критерия
which.min(reg.summary$cp)## [1] 4### 4
# график
plot(reg.summary$cp, xlab = 'Число предикторов',
ylab = 'C_p', type = 'b')
points(which.min(reg.summary$cp),
reg.summary$cp[which.min(reg.summary$cp)],
col = 'red', cex = 2, pch = 20)
# BIC
reg.summary$bic## [1] -126.9960 -151.0834 -156.9115 -157.0723 -152.7167 -148.1408# число предикторов у оптимального значения критерия
which.min(reg.summary$bic)## [1] 4### 4
# график
plot(reg.summary$bic, xlab = 'Число предикторов',
ylab = 'BIC', type = 'b')
points(which.min(reg.summary$bic),
reg.summary$bic[which.min(reg.summary$bic)],
col = 'red', cex = 2, pch = 20)
# метод plot для визуализации результатов
plot(regfit.full, scale = 'r2')
plot(regfit.full, scale = 'adjr2')
plot(regfit.full, scale = 'Cp')
plot(regfit.full, scale = 'bic')
# коэффициенты модели с наименьшим BIC
round(coef(regfit.full, 4), 3)## (Intercept) education income census typeprof
## -11.238 3.996 0.001 0.001 10.226Нахождение оптимальной модели при помощи метода перекрёстной проверки
k-кратная кросс-валидация
# отбираем 10 блоков наблюдений
k <- 10
set.seed(my.seed)
folds <- sample(1:k, nrow(Prestige), replace = T)
# заготовка под матрицу с ошибками
cv.errors <- matrix(NA, k, 6, dimnames = list(NULL, paste(1:6)))
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(prestige ~ ., data = Prestige[folds != j, ],
nvmax = 6)
# теперь цикл по количеству объясняющих переменных
for (i in 1:6){
# модельные значения prestige
pred <- predict(best.fit, Prestige[folds == j, ], id = i)
# вписываем ошибку в матрицу
cv.errors[j, i] <- mean((Prestige$prestige[folds == j] - pred)^2)
}
}
# усредняем матрицу по каждому столбцу (т.е. по блокам наблюдений),
# чтобы получить оценку MSE для каждой модели с фиксированным
# количеством объясняющих переменных
mean.cv.errors <- apply(cv.errors, 2, mean)
round(mean.cv.errors, 0)## 1 2 3 4 5 6
## 77 58 53 48 50 50# на графике
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)
# перестраиваем модель с 4 объясняющими переменными на всём наборе данных
reg.best <- regsubsets(prestige ~ ., data = Prestige, nvmax = 7)
round(coef(reg.best, 4), 3)## (Intercept) education income census typeprof
## -11.238 3.996 0.001 0.001 10.226Гребневая регрессия
# из-за синтаксиса glmnet() формируем явно матрицу объясняющих...
x <- model.matrix(prestige ~ ., Prestige)[, -1]
# и вектор значений зависимой переменной
y <- Prestige$prestige# вектор значений гиперпараметра лямбда
grid <- 10^seq(10, -2, length = 100)
# подгоняем ридж-модели с большей точностью (thresh ниже значения по умолчанию) на всех данных
ridge.mod <- glmnet(x, y, alpha = 0, lambda = grid,
thresh = 1e-12)
plot(ridge.mod)
Подбор оптимального значения лямбда с помощью перекрёстной проверки
# k-кратная кросс-валидация
set.seed(my.seed)
# оценка ошибки
cv.out <- cv.glmnet(x, y, alpha = 0)
plot(cv.out)
# значение лямбда, обеспечивающее минимальную ошибку перекрёстной проверки
bestlam <- cv.out$lambda.min
round(bestlam, 0)## [1] 1# наконец, подгоняем модель для оптимальной лямбды,
# найденной по перекрёстной проверке
out <- glmnet(x, y, alpha = 0)
round(predict(out, type = 'coefficients', s = bestlam)[1:7, ], 3)## (Intercept) education income women census typeprof
## 3.347 3.015 0.001 0.006 0.000 9.767
## typewc
## -0.754set.seed(my.seed)
train <- sample(1:nrow(x), nrow(x)/2)
test <- -train
y.test <- y[test]
# модель с оптимальным значением параметра на обучающей выборке
ridge.train <- glmnet(x[train, ], y[train], alpha = 0, lambda = grid,
thresh = 1e-12)
round(predict(ridge.train, type = 'coefficients', s = bestlam)[1:7, ], 3)## (Intercept) education income women census typeprof
## 2.886 3.234 0.001 -0.019 0.000 9.988
## typewc
## -1.605# summary по всей выборке
ridge.mod1 <- glmnet(x, y, alpha = 0, lambda = bestlam)
ridge.sum <- summary(ridge.mod1)
ridge.sum## Length Class Mode
## a0 1 -none- numeric
## beta 6 dgCMatrix S4
## df 1 -none- numeric
## dim 2 -none- numeric
## lambda 1 -none- numeric
## dev.ratio 1 -none- numeric
## nulldev 1 -none- numeric
## npasses 1 -none- numeric
## jerr 1 -none- numeric
## offset 1 -none- logical
## call 5 -none- call
## nobs 1 -none- numericy_predicted <- predict(ridge.mod1, newx = x)
sst <- sum((y-mean(y))^2)
sse <- sum((y_predicted-y)^2)
# R^2
rsq <- 1-sse/sst
rsq## [1] 0.8340552Сравнение оптимальных моделей, полученных в заданиях 1 и 2 по MSE на тестовой выборке (объём выобрки - 50%)
# MSE на тестовой выборке с 4 объясняющими переменными (отбор оптимального подмножества)
opt.test <- predict(best.fit, Prestige[test, ], id = 4)
opt.mse.test <- round(mean((opt.test - y.test)^2), 0)
# MSE на тестовой выборке (гребневая регрессия)
ridge.test <- predict(ridge.mod, s = bestlam, newx = x[test, ])
ridge.mse.test <- round(mean((ridge.test - y.test)^2), 0)
MSE.test <- rbind(opt.mse.test, ridge.mse.test)
row.names(MSE.test) <- c('MSE (отбор оптимального подмножества)', 'MSE (гребневая регрессия)')
kable(MSE.test)| MSE (отбор оптимального подмножества) | 54 |
| MSE (гребневая регрессия) | 55 |
Сравнивая результаты расчётов MSE на тестовой выборке для двух оптимальных моделей, полученных в заданиях 1 и 2, можно заключить, что стандартная ошибка MSE модели №1 (отбор оптимального подмножества) оказалась меньше, чем MSE модели №2. Таким образом, модель №1 (отбор оптимального подмножества) оказалась лучшей.