Лабораторная №6

Задачи:

Примените указанный в варианте метод к набору данных по своему варианту (см. таблицу ниже). Не забудьте предварительно сделать из категориальных переменных факторы. Выберите оптимальную модель с помощью кросс-валидации. Выведите её коэффициенты с помощью функции coef(). Рассчитайте MSE модели на тестовой выборке.
Примените указанный в варианте метод к набору данных по своему варианту (см. таблицу ниже). Для модели:

Подогнать модель на всей выборке и вычислить ошибку (MSE) с кросс-валидацией. По наименьшей MSE подобрать оптимальное значение настроечного параметра метода (гиперпараметр λ или число главных компонент M).
Подогнать модель с оптимальным значением параметра на обучающей выборке, посчитать MSE на тестовой.
Подогнать модель с оптимальным значением параметра на всех данных, вывести характеристики модели функцией summary().

Сравните оптимальные модели, полученные в заданиях 1 и 2 по MSE на тестовой выборке. Какой метод дал лучший результат? Доля тестовой выборки: 50%.

Вариант 10

Данные College {ISLR}:

Accept - Количество принятых заявок;
Private - Фактор с уровнями Нет и Да, указывающий частный или государственный университет;
Apps - Количество полученных заявок;
Enroll - Количество новых студентов, зачисленных;
Top10perc - Процент новых студентов из лучших 10% H.S. класс;
Top25perc - Процент новых студентов из лучших 25% H.S. класс;
F.Undergrad - Количество студентов очной формы обучения;
P.Undergrad - Количество заочных студентов;
Outstate - Обучение за пределами штата;
Room.Board - Стоимость проживания и питания;
Book - Ориентировочная стоимость книги;
Personal - Расчетные личные расходы;
PhD - Процент факультета с докторской степенью;
Terminal - Процент факультета с конечной степенью;
S.F.Ratio - Соотношение студентов и преподавателей;
perc.alumni - Процент выпускников, которые жертвуют;
Expend - Расходы на одного учащегося;
Grad.Rate - Выпускной.

Методы:

Для задания 1: отбор оптимального подмножества;
Для задания 2: частный метод наименьших квадратов.

# Загрузка пакетов
library('knitr')             # Пакет для генерации отчёта
library('ISLR')              # Набор данных College
library('leaps')             # Функция regsubset() - отбор оптимального подмножества переменных
library('pls')               # Частный метод наименьших квадратов - pls()

my.seed <- 1

# Загрузка данных College
data('College')
# Переводим дискретные количественные переменные в факторы
College$Peivate <- as.factor(College$Private)

Статистика по большому количеству американских колледжей из выпуска US News and World Report за 1995 год.

# Название столбцов переменных
names(College)

##  [1] "Private"     "Apps"        "Accept"      "Enroll"      "Top10perc"  
##  [6] "Top25perc"   "F.Undergrad" "P.Undergrad" "Outstate"    "Room.Board" 
## [11] "Books"       "Personal"    "PhD"         "Terminal"    "S.F.Ratio"  
## [16] "perc.alumni" "Expend"      "Grad.Rate"   "Peivate"

# Размерность данных
dim(College)

## [1] 777  19

Считаем число пропусков в данных и убираем их.

# Считаем пропуски
sum(is.na(College))

## [1] 0

# Убираем пропуски
College <- na.omit(College)

# Проверяем результат
dim(College)

## [1] 777  19

sum(is.na(College))

## [1] 0

Отбор оптимального подмножества

# Подгоняем модели с сочетаниями предикторов до 8 (максимум в данных)
regfit.full <- regsubsets(Accept ~ ., College)

## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax,
## force.in = force.in, : 1 linear dependencies found

reg.summary <- summary(regfit.full)
reg.summary

## Subset selection object
## Call: regsubsets.formula(Accept ~ ., College)
## 18 Variables  (and intercept)
##             Forced in Forced out
## PrivateYes      FALSE      FALSE
## Apps            FALSE      FALSE
## Enroll          FALSE      FALSE
## Top10perc       FALSE      FALSE
## Top25perc       FALSE      FALSE
## F.Undergrad     FALSE      FALSE
## P.Undergrad     FALSE      FALSE
## Outstate        FALSE      FALSE
## Room.Board      FALSE      FALSE
## Books           FALSE      FALSE
## Personal        FALSE      FALSE
## PhD             FALSE      FALSE
## Terminal        FALSE      FALSE
## S.F.Ratio       FALSE      FALSE
## perc.alumni     FALSE      FALSE
## Expend          FALSE      FALSE
## Grad.Rate       FALSE      FALSE
## PeivateYes      FALSE      FALSE
## 1 subsets of each size up to 8
## Selection Algorithm: exhaustive
##          PrivateYes Apps Enroll Top10perc Top25perc F.Undergrad
## 1  ( 1 ) " "        "*"  " "    " "       " "       " "        
## 2  ( 1 ) " "        "*"  "*"    " "       " "       " "        
## 3  ( 1 ) " "        "*"  "*"    "*"       " "       " "        
## 4  ( 1 ) " "        "*"  "*"    "*"       " "       " "        
## 5  ( 1 ) " "        "*"  "*"    "*"       " "       " "        
## 6  ( 1 ) " "        "*"  "*"    "*"       "*"       " "        
## 7  ( 1 ) " "        "*"  "*"    "*"       "*"       " "        
## 8  ( 1 ) " "        "*"  "*"    "*"       "*"       " "        
##          P.Undergrad Outstate Room.Board Books Personal PhD Terminal
## 1  ( 1 ) " "         " "      " "        " "   " "      " " " "     
## 2  ( 1 ) " "         " "      " "        " "   " "      " " " "     
## 3  ( 1 ) " "         " "      " "        " "   " "      " " " "     
## 4  ( 1 ) " "         "*"      " "        " "   " "      " " " "     
## 5  ( 1 ) " "         "*"      " "        " "   " "      " " " "     
## 6  ( 1 ) " "         "*"      " "        " "   " "      " " " "     
## 7  ( 1 ) "*"         "*"      " "        " "   " "      " " " "     
## 8  ( 1 ) "*"         "*"      " "        " "   " "      "*" " "     
##          S.F.Ratio perc.alumni Expend Grad.Rate PeivateYes
## 1  ( 1 ) " "       " "         " "    " "       " "       
## 2  ( 1 ) " "       " "         " "    " "       " "       
## 3  ( 1 ) " "       " "         " "    " "       " "       
## 4  ( 1 ) " "       " "         " "    " "       " "       
## 5  ( 1 ) " "       " "         "*"    " "       " "       
## 6  ( 1 ) " "       " "         "*"    " "       " "       
## 7  ( 1 ) " "       " "         "*"    " "       " "       
## 8  ( 1 ) " "       " "         "*"    " "       " "

# Структура отчёта по модели (ищем характеристики качества)
names(reg.summary)

## [1] "which"  "rsq"    "rss"    "adjr2"  "cp"     "bic"    "outmat" "obj"

# R^2 и скорректированный R^2
round(reg.summary$rsq, 3)

## [1] 0.890 0.935 0.944 0.948 0.950 0.952 0.952 0.952

# На графике
plot(1:8, reg.summary$rsq, type = 'b',
     xlab = 'Number of predictors', ylab = 'R-square')
# Сюда же добавим скорректированный R-квадрат
points(1:8, reg.summary$adjr2, col = 'red')

# Модель с максимальным скорректированным R-квадратом
which.max(reg.summary$adjr2)

## [1] 8

### 8
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] 1008.13316  282.46065  145.55236   70.49357   37.68531   22.33299
## [7]   18.79450   15.54451

# число предикторов у оптимального значения критерия
which.min(reg.summary$cp)

## [1] 8

### 8

# График
plot(reg.summary$cp, xlab = 'Number of predictors',
     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] -1702.441 -2103.842 -2207.055 -2268.792 -2295.121 -2305.447 -2304.291
## [8] -2302.884

# Число предикторов у оптимального значения критерия
which.min(reg.summary$bic)

## [1] 6

### 6

# График
plot(reg.summary$bic, xlab = 'Number of predictors',
     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)        Apps      Enroll   Top10perc    Outstate 
##     -36.275       0.410       1.068     -20.904       0.054

Нахождение оптимальной модели при помощи метода перекрёстной проверки

k-кратная кросс-валидация

# Отбираем 10 блоков наблюдений
k <- 10
set.seed(my.seed)
folds <- sample(1:k, nrow(Auto), replace = T)

# Заготовка под матрицу с ошибками
cv.errors <- matrix(NA, k, 8, dimnames = list(NULL, paste(1:8)))

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(Accept ~ ., data = College[folds != j, ],
                           nvmax = 8)
    # Теперь цикл по количеству объясняющих переменных
    for (i in 1:8){
        # Модельные значения mpg
        pred <- predict(best.fit, College[folds == j, ], id = i)
        # Вписываем ошибку в матрицу
        cv.errors[j, i] <- mean((College$Accept[folds == j] - pred)^2)
    }
}

## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax,
## force.in = force.in, : 1 linear dependencies found

## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax,
## force.in = force.in, : 1 linear dependencies found

## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax,
## force.in = force.in, : 1 linear dependencies found

## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax,
## force.in = force.in, : 1 linear dependencies found

## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax,
## force.in = force.in, : 1 linear dependencies found

## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax,
## force.in = force.in, : 1 linear dependencies found

## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax,
## force.in = force.in, : 1 linear dependencies found

## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax,
## force.in = force.in, : 1 linear dependencies found

## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax,
## force.in = force.in, : 1 linear dependencies found

## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax,
## force.in = force.in, : 1 linear dependencies found

# Усредняем матрицу по каждому столбцу (т.е. по блокам наблюдений), 
# Чтобы получить оценку MSE для каждой модели с фиксированным 
# Количеством объясняющих переменных
mean.cv.errors <- apply(cv.errors, 2, mean)
round(mean.cv.errors, 0)

##      1      2      3      4      5      6      7      8 
## 662556 411038 352717 330058 325171 314418 320946 321512

# На графике
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)

# Перестраиваем модель с 8 объясняющими переменными на всём наборе данных
reg.best <- regsubsets(Accept ~ ., data = College, nvmax = 8)

## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax,
## force.in = force.in, : 1 linear dependencies found

round(coef(reg.best, 8), 3)

## (Intercept)        Apps      Enroll   Top10perc   Top25perc P.Undergrad 
##    -395.183       0.419       1.050     -27.915       8.558      -0.040 
##    Outstate         PhD      Expend 
##       0.066       3.500      -0.029

Регрессия по методу частных наименьших квадратов

set.seed(my.seed)
x <- model.matrix(Accept ~ ., College)[, -1]
train <- sample(1:nrow(x), nrow(x)/2)
test <- -train
y <- College$Accept
y.test <- y[test]
pls.fit <- plsr(Accept ~ ., data = College, subset = train, scale = T,
                validation = 'CV')
summary(pls.fit)

## Data:    X dimension: 388 18 
##  Y dimension: 388 1
## Fit method: kernelpls
## Number of components considered: 18
## 
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
##        (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
## CV            2619     1291    866.1    796.9    693.9    659.7    609.5
## adjCV         2619     1287    844.9    800.5    685.1    653.6    604.8
##        7 comps  8 comps  9 comps  10 comps  11 comps  12 comps  13 comps
## CV       593.6    590.7    589.5     593.1     596.0     600.0     598.3
## adjCV    590.1    586.5    585.3     588.5     591.2     594.6     593.0
##        14 comps  15 comps  16 comps  17 comps  18 comps
## CV        597.8     597.5     597.5     597.5     597.5
## adjCV     592.5     592.3     592.3     592.3     592.3
## 
## TRAINING: % variance explained
##         1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps
## X         26.67    34.21    63.25    66.31    71.10    74.71    78.89
## Accept    77.35    91.72    92.87    95.05    95.52    95.87    95.94
##         8 comps  9 comps  10 comps  11 comps  12 comps  13 comps  14 comps
## X         80.98    83.53     86.71     89.40     90.89     92.79     94.74
## Accept    96.02    96.07     96.09     96.11     96.14     96.15     96.15
##         15 comps  16 comps  17 comps  18 comps
## X          96.29     98.82    100.00    100.64
## Accept     96.15     96.15     96.15     96.15

# График ошибок
validationplot(pls.fit, val.type = 'MSEP')

# Теперь подгоняем модель для найденного оптимального M = 8
# и оцениваем MSE на тестовой выборке
pls.pred <- predict(pls.fit, x[test, ], ncomp = 8)
round(mean(pls.pred - y.test^2), 0)

## [1] -8843376

# Подгоняем модель на всей выборке
pls.fit <- plsr(Accept ~ ., data = College, scale = T, ncomp = 2)
summary(pls.fit)

## Data:    X dimension: 777 18 
##  Y dimension: 777 1
## Fit method: kernelpls
## Number of components considered: 2
## TRAINING: % variance explained
##         1 comps  2 comps
## X         26.45    38.65
## Accept    78.56    89.81

# MSE на тестовой выборке с 8 объясняющими переменными (отбор оптимального подмножества)
opt.test <- predict(best.fit, College[test, ], id = 8)
opt.mse.test <- round(mean((opt.test - y.test)^2), 0)

# MSE на тестовой выборке (частный метод наименьших квадратов)
sqr.test <- predict(pls.fit, x[test, ], ncomp = 2)
sqr.mse.test <- round(mean((pls.pred - y.test)^2), 0)

MSE.test <- rbind(opt.mse.test, sqr.mse.test)
row.names(MSE.test) <- c('MSE (selection of the optimal subset)', 'MSE (Partial Least Squares)')
kable(MSE.test)

MSE (selection of the optimal subset)	292197
MSE (Partial Least Squares)	315316

Сравнивая результаты расчётов MSE на тестовой выборке для двух оптимальных моделей, полученных в заданиях 1 и 2, можно заключить, что стандартная ошибка MSE модели №1 (отбор оптимального подмножества) оказалась меньше, чем MSE модели №2. Таким образом, модель №1 (отбор оптимального подмножества) оказалась лучшей.