1 Оценить стандартную ошибку модели для линейных регрессионных моделей из упражнения 4 (варианты ниже): а) со всеми объясняющими переменными; б) только с непрерывными объясняющими переменными:
- методом проверочной выборки с долей обучающей 50%;
- методом LOOCV;
- k-кратной кросс-валидацией с \(k = 5\) и \(k = 10\).
Выбрать лучшую модель по минимуму ошибки. Все ли методы кросс-валидации сходятся на одной и той же модели?

2 Оценить стандартные ошибки параметров лучшей модели регрессии методом бутстрепа. Сравнить с оценками стандартных ошибок параметров по МНК.

Вариант 10

Модели: линейная регрессия.
Данные: `Carseats {ISLR}’.

Набор данных Carseats содержит переменные:

# загрузка пакетов
library('knitr')             # пакет для генерации отчёта
library('ISLR')              # набор данных Carseats
library('GGally')            # матричные графики
library('boot')              # расчёт ошибки с кросс-валидацией

my.seed <- 1

# загрузка данных Carseats
data('Carseats')
# отбор необходимых данных для построения моделей
Carseats <- Carseats[,c('Sales', 'Price', 'Advertising', 'US'),drop=FALSE]

Рассмотрим данные с характеристиками автомобилей Carseats из пакета ISLR. Скопируем таблицу во фрейм DF.carseats для дальнейших манипуляций.

# запись данных в фрейм
DF.carseats <- Carseats

# просмотр первых записей
head(DF.carseats)
# описательные статистики
summary(DF.carseats)
##      Sales            Price        Advertising       US     
##  Min.   : 0.000   Min.   : 24.0   Min.   : 0.000   No :142  
##  1st Qu.: 5.390   1st Qu.:100.0   1st Qu.: 0.000   Yes:258  
##  Median : 7.490   Median :117.0   Median : 5.000            
##  Mean   : 7.496   Mean   :115.8   Mean   : 6.635            
##  3rd Qu.: 9.320   3rd Qu.:131.0   3rd Qu.:12.000            
##  Max.   :16.270   Max.   :191.0   Max.   :29.000

В таблице данных 400 наблюдений и 4 переменных, среди которых есть непрерывные количественные и одна дискретная (US, коэффициент, указывающий на то, расположен ли магазин в US. Построим графики разброса, показав фактор US цветом. Зависимой переменной модели будет Sales, её покажем в первой строке / столбце матричного графика.

# переводим дискретные количественные переменные в факторы
Carseats$US <- as.factor(Carseats$US)

# графики разброса, цвет -- US
ggpairs(DF.carseats, ggplot2::aes(color = US))

Метод проверочной выборки

Он состоит в том, что мы отбираем одну тестовую выборку и будем считать на ней ошибку модели.

# общее число наблюдений
n <- nrow(DF.carseats)

# доля обучающей выборки
train.percent <- 0.5

# выбрать наблюдения в обучающую выборку
set.seed(my.seed)
inTrain <- sample(n, n * train.percent)

# рисуем разными цветами обучающую и тестовую (для непрерывных переменных)

# Price
par(mfrow = c(1, 2))
plot(DF.carseats$Price[inTrain], DF.carseats$Sales[inTrain],
     xlab = 'Price', ylab = 'Sales', pch = 21,
     col = rgb(0, 0, 1, alpha = 0.4), bg = rgb(0, 0, 1, alpha = 0.4))
points(DF.carseats$Price[-inTrain], DF.carseats$Sales[-inTrain],
       pch = 21, col = rgb(1, 0, 0, alpha = 0.4), 
       bg = rgb(1, 0, 0, alpha = 0.4))
legend('topright', 
       pch = c(16, 16), col = c('blue', 'red'), legend = c('test', 'train'))

# Advertising
plot(DF.carseats$Advertising[inTrain], DF.carseats$Sales[inTrain],
     xlab = 'Advertising', ylab = 'Sales', pch = 21,
     col = rgb(0, 0, 1, alpha = 0.4), bg = rgb(0, 0, 1, alpha = 0.4))
points(DF.carseats$Advertising[-inTrain], DF.carseats$Sales[-inTrain],
       pch = 21, col = rgb(1, 0, 0, alpha = 0.4), 
       bg = rgb(1, 0, 0, alpha = 0.4))
legend('topright', 
       pch = c(16, 16), col = c('blue', 'red'), legend = c('test', 'train'))

par(mfrow = c(1, 1))

Построим модели для проверки точности. Вид моделей:

а) Со всеми объясняющими переменными \[ \hat{Sales} = f(Price, Advertising, US); \] б) Только с непрерывными объясняющими переменными \[ \hat{Sales} = f(Price, Advertising). \]

Линейная модель a): \(\hat{Sales} = \hat{\beta}_0 + \hat{\beta}_1 \cdot Price + \hat{\beta}_2 \cdot Advertising + \hat{\beta}_3 \cdot US\).

# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(DF.carseats)

# подгонка линейной модели на обучающей выборке
fit.lm.1 <- lm(Sales ~ Price + Advertising + US, subset = inTrain)

# считаем MSE на тестовой выборке
mean((Sales[-inTrain] - predict(fit.lm.1, DF.carseats[-inTrain, ]))^2)
## [1] 5.997982
# отсоединить таблицу с данными
detach(DF.carseats)
# сохраняем ошибку модели (MSE) на проверочной выборке
err.test <- mean((DF.carseats$Sales[-inTrain] - predict(fit.lm.1, 
                                                  DF.carseats[-inTrain, ]))^2)

# сохранять все ошибки будем в один вектор, присваиваем имя первому элементу
#  (имя -- степень объясняющей переменной)
names(err.test) <- 1

Линейная модель б): \(\hat{Sales} = \hat{\beta}_0 + \hat{\beta}_1 \cdot Price + \hat{\beta}_2 \cdot Advertising\).

# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(DF.carseats)

# подгонка линейной модели на обучающей выборке
fit.lm.2 <- lm(Sales ~ Price + Advertising, subset = inTrain)

# считаем MSE на тестовой выборке
mean((Sales[-inTrain] - predict(fit.lm.2, DF.carseats[-inTrain, ]))^2)
## [1] 5.700441
# отсоединить таблицу с данными
detach(DF.carseats)
# сохраняем ошибку модели (MSE) на проверочной выборке
err.test <- c(err.test,
              mean((DF.carseats$Sales[-inTrain] - predict(fit.lm.2,
                                                 DF.carseats[-inTrain, ]))^2))

# имя второго элемента вектора
names(err.test)[length(err.test)] <- 2

Перекрёстная проверка по отдельным наблюдениям (LOOCV)

Это самый затратный в вычислительном плане метод, но и самый надёжный в плане оценки ошибки вне выборки. Попробуем применить его к линейной модели а).

# подгонка линейной модели на обучающей выборке
fit.glm1 <- glm(Sales ~ Price + Advertising + US, data = DF.carseats)

# считаем LOOCV-ошибку
cv.err.loocv <- cv.glm(DF.carseats, fit.glm1)$delta[1]

# сохранять все ошибки будем в один вектор, присваиваем имя первому элементу
#  (имя -- степень объясняющей переменной)
names(cv.err.loocv) <- 1

Теперь оценим точность линейной модели б).

# подгонка линейной модели на обучающей выборке
fit.glm2 <- glm(Sales ~ Price + Advertising, data = DF.carseats)

# считаем LOOCV-ошибку
cv.err.loocv <- c(cv.err.loocv, cv.glm(DF.carseats, fit.glm2)$delta[1])

# сохранять все ошибки будем в один вектор, присваиваем имя второму элементу
names(cv.err.loocv)[length(cv.err.loocv)] <- 2

# результат
cv.err.loocv
##        1        2 
## 5.832970 5.802091

k-кратная перекрёстная проверка

K-кратная кросс-валидация – компромисс между методом проверочной выборки и LOOCV. Оценка ошибки вне выборки ближе к правде, по сравнению с проверочной выборкой, а объём вычислений меньше, чем при LOOCV. Проведём 5-кратную кросс-валидацию моделей а) и б).

# оценим точность линейных моделей а) и б)
# вектор с ошибками по 5-кратной кросс-валидации
cv.err.k.fold5 <- rep(0, 2)

# имена элементов вектора
names(cv.err.k.fold5) <- 1:2

# оценка модели а)
fit.glm <- glm(Sales ~ Price + Advertising + US, data = DF.carseats)
# расчёт ошибки
cv.err.k.fold5[1] <- cv.glm(DF.carseats, fit.glm, K = 5)$delta[1]

# оценка модели б)
fit.glm <- glm(Sales ~ Price + Advertising, data = DF.carseats)
# расчёт ошибки
cv.err.k.fold5[2] <- cv.glm(DF.carseats, fit.glm, K = 5)$delta[1]

# результат
cv.err.k.fold5
##        1        2 
## 5.867203 5.823899

Теперь проведём 5-кратную кросс-валидацию моделей а) и б).

# оценим точность линейных моделей а) и б)
# вектор с ошибками по 10-кратной кросс-валидации
cv.err.k.fold10 <- rep(0, 2)

# имена элементов вектора
names(cv.err.k.fold10) <- 1:2

# оценка модели а)
fit.glm <- glm(Sales ~ Price + Advertising + US, data = DF.carseats)
# расчёт ошибки
cv.err.k.fold10[1] <- cv.glm(DF.carseats, fit.glm, K = 10)$delta[1]

# оценка модели б)
fit.glm <- glm(Sales ~ Price + Advertising, data = DF.carseats)
# расчёт ошибки
cv.err.k.fold10[2] <- cv.glm(DF.carseats, fit.glm, K = 10)$delta[1]

# результат
cv.err.k.fold10
##        1        2 
## 5.849068 5.845285

Для определения лучшей модели по стандартной ошибке MSE объединим все полученные результаты в таблицу.

MSE.tbl <- rbind(err.test, cv.err.loocv, cv.err.k.fold5, cv.err.k.fold10)
colnames(MSE.tbl)<-c('Модель а)', 'Модель б)')
row.names(MSE.tbl) <- c('Проверочная выборка', 'LOOCV', '5-кратная кросс-валидация', '10-кратная кросс-валидация')
kable(MSE.tbl)
<U+041C><U+043E><U+0434><U+0435><U+043B><U+044C> <U+0430>) <U+041C><U+043E><U+0434><U+0435><U+043B><U+044C> <U+0431>)
<U+041F><U+0440><U+043E><U+0432><U+0435><U+0440><U+043E><U+0447><U+043D><U+0430><U+044F> <U+0432><U+044B><U+0431><U+043E><U+0440><U+043A><U+0430> 5.997982 5.700441
LOOCV 5.832970 5.802091
5-<U+043A><U+0440><U+0430><U+0442><U+043D><U+0430><U+044F> <U+043A><U+0440><U+043E><U+0441><U+0441>-<U+0432><U+0430><U+043B><U+0438><U+0434><U+0430><U+0446><U+0438><U+044F> 5.867203 5.823899
10-<U+043A><U+0440><U+0430><U+0442><U+043D><U+0430><U+044F> <U+043A><U+0440><U+043E><U+0441><U+0441>-<U+0432><U+0430><U+043B><U+0438><U+0434><U+0430><U+0446><U+0438><U+044F> 5.849068 5.845285

Опираясь на результаты расчётов с проверочной выборкой, LOOCV и кросс-валидацией (\(k = 5\) и \(k = 10\)), можно заключить, что стандартная ошибка MSE линейной модели б) (только с непрерывными объясняющими переменными) оказалась меньше по всем методам кросс-валидации, чем MSE линейной модели а) (со всеми объясняющими переменными). Таким образом, линейная модель б) оказалась лучшей: \(\hat{Sales} = \hat{\beta}_0 + \hat{\beta}_1 \cdot Price + \hat{\beta}_2 \cdot Advertising\).

Бутстреп

Точность оценки параметра регрессии

При построении модели регрессии проблемы в остатках приводят к неверной оценке ошибок параметров. Обойти эту проблему можно, применив для расчёта этих ошибок бутстреп.

# Оценивание точности лучшей линейной регрессионной модели ----------------------------

#  оценить стандартные ошибки параметров модели 
#  Sales = beta_0 + beta_1 * Price + beta_2 * Advertising с помощью бутстрепа,
#  сравнить с оценками ошибок по МНК

# функция для расчёта коэффициентов ЛР по выборке из данных
boot.fn <- function(data, index){
  coef(lm(Sales ~ Price + Advertising, data = data, subset = index))
}
boot.fn(DF.carseats, 1:n)
## (Intercept)       Price Advertising 
## 13.00342697 -0.05461304  0.12310706
# пример применения функции к бутстреп-выборке
set.seed(my.seed)
boot.fn(DF.carseats, sample(n, n, replace = T))
## (Intercept)       Price Advertising 
##  13.5130453  -0.0570975   0.1041370
# применяем функцию boot для вычисления стандартных ошибок параметров
#  (1000 выборок с повторами)
boot(DF.carseats, boot.fn, 1000)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = DF.carseats, statistic = boot.fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##        original        bias    std. error
## t1* 13.00342697  7.495898e-04 0.638162137
## t2* -0.05461304  2.023930e-06 0.005280631
## t3*  0.12310706 -8.781749e-04 0.019317299
# сравним с МНК
attach(DF.carseats)
summary(lm(Sales ~ Price + Advertising))$coef
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept) 13.00342697 0.606850308  21.427734 3.059318e-68
## Price       -0.05461304 0.005078128 -10.754562 7.596365e-24
## Advertising  0.12310706 0.018079180   6.809327 3.639822e-11
detach(DF.carseats)

В модели регрессии, для которой проводился расчёт, похоже, не нарушаются требования к остаткам, и оценки стандартных ошибок параметров, рассчитанные по МНК, очень близки к ошибкам этих же параметров, полученных бутстрепом.