Математическое моделирование

Кросс-валидация и бутстреп

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

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

Как сдавать: прислать на почту преподавателя ссылки:
- на html-отчёт с видимыми блоками кода (блоки кода с параметром echo = T), размещённый на rpubs.com. - на код, генерирующий отчёт, в репозитории на github.com.
В текст отчёта включить постановку задачи и ответы на вопросы задания.

Вариант - 13

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

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

  • Sales - объем продаж в каждом месте (в тысячах);
  • Price – взимаемая плата за автокресла на каждом участке;
  • Population – численность населения в регионе (в тысячах);
  • ShelveLoc – коэффициент, указывающий на качество места для размещения автомобильных кресел на каждом участке:

1 – плохое (Bad); 2 – хорошее (Good); 3 - среднее (Medium).

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

my.seed <- 1

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

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

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

# просмотр первых записей
head(DF.carseats)
##   Sales Price Population ShelveLoc
## 1  9.50   120        276       Bad
## 2 11.22    83        260      Good
## 3 10.06    80        269    Medium
## 4  7.40    97        466    Medium
## 5  4.15   128        340       Bad
## 6 10.81    72        501       Bad
# описательные статистики
summary(DF.carseats)
##      Sales            Price         Population     ShelveLoc  
##  Min.   : 0.000   Min.   : 24.0   Min.   : 10.0   Bad   : 96  
##  1st Qu.: 5.390   1st Qu.:100.0   1st Qu.:139.0   Good  : 85  
##  Median : 7.490   Median :117.0   Median :272.0   Medium:219  
##  Mean   : 7.496   Mean   :115.8   Mean   :264.8               
##  3rd Qu.: 9.320   3rd Qu.:131.0   3rd Qu.:398.5               
##  Max.   :16.270   Max.   :191.0   Max.   :509.0

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

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

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

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

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

# общее число наблюдений
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'))

# Population
plot(DF.carseats$Population[inTrain], DF.carseats$Sales[inTrain],
     xlab = 'Population', ylab = 'Sales', pch = 21,
     col = rgb(0, 0, 1, alpha = 0.4), bg = rgb(0, 0, 1, alpha = 0.4))
points(DF.carseats$Population[-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, Population, ShelveLoc); \] б) Только с непрерывными объясняющими переменными \[ \hat{Sales} = f(Price, Population). \]

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

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

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

# считаем MSE на тестовой выборке
mean((Sales[-inTrain] - predict(fit.lm.1, DF.carseats[-inTrain, ]))^2)
## [1] 3.279182
# отсоединить таблицу с данными
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 Population\).

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

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

# считаем MSE на тестовой выборке
mean((Sales[-inTrain] - predict(fit.lm.2, DF.carseats[-inTrain, ]))^2)
## [1] 6.55673
# отсоединить таблицу с данными
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 + Population + ShelveLoc, 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 + Population, 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 
## 3.708368 6.459368

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+ Population + ShelveLoc, data = DF.carseats)
# расчёт ошибки
cv.err.k.fold5[1] <- cv.glm(DF.carseats, fit.glm, K = 5)$delta[1]

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

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

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

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

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

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

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

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

Для определения лучшей модели по стандартной ошибке 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)
Модель а) Модель б)
Проверочная выборка 3.279182 6.556730
LOOCV 3.708368 6.459368
5-кратная кросс-валидация 3.710963 6.467952
10-кратная кросс-валидация 3.704528 6.514330

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

Бутстреп

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

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

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

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

# функция для расчёта коэффициентов ЛР по выборке из данных
boot.fn <- function(data, index){
  coef(lm(Sales ~ Price + Population + ShelveLoc, data = data, subset = index))
}
boot.fn(DF.carseats, 1:n)
##     (Intercept)           Price      Population   ShelveLocGood ShelveLocMedium 
##    11.715906257    -0.056625562     0.001008605     4.903903216     1.877949406
# пример применения функции к бутстреп-выборке
set.seed(my.seed)
boot.fn(DF.carseats, sample(n, n, replace = T))
##     (Intercept)           Price      Population   ShelveLocGood ShelveLocMedium 
##    11.781567367    -0.058350146     0.001080649     4.879445141     2.036233247
# применяем функцию 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* 11.715906257  8.551657e-03 0.5604663005
## t2* -0.056625562 -1.271621e-05 0.0039743701
## t3*  0.001008605 -1.835743e-05 0.0006328749
## t4*  4.903903216 -9.411253e-03 0.2908960551
## t5*  1.877949406 -1.952405e-03 0.2394661598
# сравним с МНК
attach(DF.carseats)
summary(lm(Sales ~ Price + Population + ShelveLoc))$coef
##                     Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)     11.715906257 0.535347674  21.884668 4.154355e-70
## Price           -0.056625562 0.004052425 -13.973254 2.451630e-36
## Population       0.001008605 0.000650822   1.549741 1.220046e-01
## ShelveLocGood    4.903903216 0.285463175  17.178759 8.471642e-50
## ShelveLocMedium  1.877949406 0.234558505   8.006316 1.334818e-14
detach(DF.carseats)

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