R Notebook

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

Практика 5

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

В практических примерах ниже показано:

• как оценить точность модели методом перекрёстной выборки;
  
• методом проверочной выборки;
  
• методом перекрёстной проверки по отдельным наблюдениям (LOOCV);
  
• методом k-кратной перекрёстной проверки;
  
• как применять бутстреп для оценки точности статистического параметра и оценок параметров модели

Модели: линейная регрессия, kNN.
Данные: Auto {ISLR}, Portfolio {ISLR}

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

my.seed <- 1

Метод перекрёстной проверки

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

# Пример на данных по автомобилям: Auto {ISLR}
DF.auto <- Auto

head(DF.auto)

summary(DF.auto)

##       mpg          cylinders      displacement     horsepower        weight    
##  Min.   : 9.00   Min.   :3.000   Min.   : 68.0   Min.   : 46.0   Min.   :1613  
##  1st Qu.:17.00   1st Qu.:4.000   1st Qu.:105.0   1st Qu.: 75.0   1st Qu.:2225  
##  Median :22.75   Median :4.000   Median :151.0   Median : 93.5   Median :2804  
##  Mean   :23.45   Mean   :5.472   Mean   :194.4   Mean   :104.5   Mean   :2978  
##  3rd Qu.:29.00   3rd Qu.:8.000   3rd Qu.:275.8   3rd Qu.:126.0   3rd Qu.:3615  
##  Max.   :46.60   Max.   :8.000   Max.   :455.0   Max.   :230.0   Max.   :5140  
##                                                                                
##   acceleration        year           origin                      name    
##  Min.   : 8.00   Min.   :70.00   Min.   :1.000   amc matador       :  5  
##  1st Qu.:13.78   1st Qu.:73.00   1st Qu.:1.000   ford pinto        :  5  
##  Median :15.50   Median :76.00   Median :1.000   toyota corolla    :  5  
##  Mean   :15.54   Mean   :75.98   Mean   :1.577   amc gremlin       :  4  
##  3rd Qu.:17.02   3rd Qu.:79.00   3rd Qu.:2.000   amc hornet        :  4  
##  Max.   :24.80   Max.   :82.00   Max.   :3.000   chevrolet chevette:  4  
##                                                  (Other)           :365

В таблице данных 392 наблюдений и 9 переменных, среди которых есть непрерывные количественные, дискретные количественные и одна номинальная (name, название модели автомобиля, сохранено как фактор). В данном случае по функции summary() сложно определить реальные типы переменных, помогает table() от отдельных столбцов таблицы: если уникальных значений немного, перед нами фактор.

# количество цилиндров
table(DF.auto$cylinders)

## 
##   3   4   5   6   8 
##   4 199   3  83 103

# происхождение автомобиля: 1 - США, 2 - Европа, 3 - Япония
table(DF.auto$origin)

## 
##   1   2   3 
## 245  68  79

Построим графики разброса, показав факторы cylinders (число цилиндров) и origin (регион происхождения автомобиля) цветом. Зависимой переменной модели будет mpg, её покажем в первой строке / столбце матричного графика. Во вторую строку / столбец поставим фактор.

# переводим дискретные количественные переменные в факторы
DF.auto$cylinders <- as.factor(DF.auto$cylinders)
DF.auto$origin <- factor(DF.auto$origin, levels = c(1, 2, 3), 
                            labels = c('Am', 'Eur', 'Jap'))

# графики разброса, цвет -- количество цилиндров
ggpairs(DF.auto[, 1:5], ggplot2::aes(color = cylinders))

ggpairs(DF.auto[, c(1:2, 6:7)], ggplot2::aes(color = cylinders))

# графики разброса, цвет -- регион происхождения
ggpairs(DF.auto[, c(1, 8, 3:5)], ggplot2::aes(color = origin))

ggpairs(DF.auto[, c(1, 8, 6:7)], ggplot2::aes(color = origin))

# только mpg ~ horsepower
plot(DF.auto$horsepower, DF.auto$mpg,
     xlab = 'horsepower', ylab = 'mpg', pch = 21,
     col = rgb(0, 0, 1, alpha = 0.4), bg = rgb(0, 0, 1, alpha = 0.4))

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

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

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

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

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

# фактические значения Y на тестовой выборке
y.test.fact <- DF.auto$mpg[-inTrain]

# рисуем разными цветами обучающую и тестовую
plot(DF.auto$horsepower[inTrain], DF.auto$mpg[inTrain],
     xlab = 'horsepower', ylab = 'mpg', pch = 21,
     col = rgb(0, 0, 1, alpha = 0.4), bg = rgb(0, 0, 1, alpha = 0.4))
points(DF.auto$horsepower[-inTrain], DF.auto$mpg[-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'))

Построим модели для проверки точности. Вид моделей:
\[ \hat{mpg} = f(horsepower) \] Линейная модель: \(\hat{mpg} = \hat{\beta}_0 + \hat{\beta}_1 \cdot horsepower\).

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

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

# подгонка линейной модели на обучающей выборке
fit.lm.1 <- lm(mpg ~ horsepower, 
               subset = inTrain)
# прогноз на тестовую
y.test.lm.1 <- predict(fit.lm.1, DF.auto[-inTrain, ])

# считаем MSE на тестовой выборке
MSE.lm.1 <- mean((y.test.fact - y.test.lm.1)^2)

# отсоединить таблицу с данными
detach(DF.auto)

# смотрим ошибку
MSE.lm.1

## [1] 23.26601

Строим квадратичную модель: \(\hat{mpg} = \hat{\beta}_0 + \hat{\beta}_1 \cdot horsepower + \hat{\beta}_2 \cdot horsepower^2\).

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

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

# прогноз на тестовую
y.test.lm.2 <- predict(fit.lm.2, DF.auto[-inTrain, ])

# считаем MSE на тестовой выборке
MSE.lm.2 <- round(mean((y.test.fact - y.test.lm.2)^2), 2)

# отсоединить таблицу с данными

detach(DF.auto)

# смотрим ошибку
MSE.lm.2

## [1] 18.72

Строим кубическую модель: \(\hat{mpg} = \hat{\beta}_0 + \hat{\beta}_1 \cdot horsepower + \hat{\beta}_2 \cdot horsepower^2 + \hat{\beta}_3 \cdot horsepower^3\).

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

# подгонка модели на обучающей выборке
fit.lm.3 <- lm(mpg ~ poly(horsepower, 3), 
               subset = inTrain)

# прогноз на тестовую
y.test.lm.3 <- predict(fit.lm.3, DF.auto[-inTrain, ])

# считаем MSE на тестовой выборке
MSE.lm.3 <- round(mean((y.test.fact - y.test.lm.3)^2), 2)

# отсоединить таблицу с данными
detach(DF.auto)

# смотрим ошибку
MSE.lm.3

## [1] 18.79

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

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

# подгонка линейной модели на обучающей выборке
fit.glm <- glm(mpg ~ horsepower, data = DF.auto)

# считаем LOOCV-ошибку
cv.err <- cv.glm(DF.auto, fit.glm)

# результат: первое число -- по формуле LOOCV-ошибки,
#  второе -- с поправкой на смещение
cv.err$delta[1]

## [1] 24.23151

Теперь оценим точность полиномиальных моделей, меняя степень, в которой стоит регрессор.

# вектор с LOOCV-ошибками
cv.err.loocv <- rep(0, 5)
# имена элементов вектора
names(cv.err.loocv) <- 1:5

# цикл по степеням полиномов
for (i in 1:5) {
  # оценка модели
  fit.glm <- glm(mpg ~ poly(horsepower, i), data = DF.auto)
  # расчёт ошибки
  cv.err.loocv[i] <- cv.glm(DF.auto, fit.glm)$delta[1]
}

# результат
cv.err.loocv

##        1        2        3        4        5 
## 24.23151 19.24821 19.33498 19.42443 19.03321

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

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

# оценим точность полиномиальных моделей, меняя степень
# вектор с ошибками по 10-кратной кросс-валидации
cv.err.k.fold <- rep(0, 5)
# имена элементов вектора
names(cv.err.k.fold) <- 1:5

# цикл по степеням полиномов
for (i in 1:5) {
  # оценка модели
  fit.glm <- glm(mpg ~ poly(horsepower, i), data = DF.auto)
  # расчёт ошибки
  cv.err.k.fold[i] <- cv.glm(DF.auto, fit.glm, K = 10)$delta[1]
}

# результат
cv.err.k.fold

##        1        2        3        4        5 
## 24.10821 19.13160 19.33517 19.53766 18.97701

Объединим все ошибки в одну таблицу и отсортируем её по возрастанию MSE:

# записываем все ошибки в таблицу
df.MSE <- data.frame(Модель = c('Линейная', 'Полином 2 степени',
                                'Полином 3 степени', 
                                rep(paste('Полином', 1:5, 'степени'), 2)), 
                     Проверка.точности = c(rep('Проверочная выборка 50%', 3),
                                           rep('LOOCV', 5), 
                                           rep('Кросс-валидация, k = 10', 5)),
                     MSE = round(c(MSE.lm.1, MSE.lm.2, MSE.lm.3, 
                                   cv.err.loocv, cv.err.k.fold), 2))

# все модели по возрастанию ошибки
df.MSE[order(df.MSE$MSE), ]

Опираясь на результаты расчётов с кросс-валидацией, можно заключить, что на самом деле ошибка вне выборки у линейной модели выше, чем показывала MSE на тестовой выборке. В целом, ошибка методом проверочной выборки размером 50% от числа наблюдений занижает MSE и, следовательно, завышает точность моделей.

Бутстреп

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

Пример с инвестиционным портфелем из двух активов: Portfolio {ISLR}. В наборе данных две переменных: * X – доход от актива X,
* Y – доход от актива Y.
Портфель инвестиций состоит из активов \(X\) и \(Y\), долю актива \(X\) обозначим как \(\alpha\). Минимум дисперсии доходности:

\[ \mathrm{Var}(\alpha X + (1 - \alpha) Y) \rightarrow \mathrm{min} \]

– достигается при значении параметра:
\[ \alpha = \frac{\sigma_Y^2 - \sigma_{XY}}{\sigma_X^2 + \sigma_Y^2 - 2\sigma_{XY}} \]

Взглянем на таблицу Portfolio:

head(Portfolio)

str(Portfolio)

## 'data.frame':    100 obs. of  2 variables:
##  $ X: num  -0.895 -1.562 -0.417 1.044 -0.316 ...
##  $ Y: num  -0.235 -0.885 0.272 -0.734 0.842 ...

Данных для оценки \(\hat{\sigma_X^2}\), \(\hat{\sigma_Y^2}\) и \(\hat{\sigma_{XY}}\) немного, всего 100 наблюдений, поэтому применим бутстреп.

# функция для вычисления искомого параметра
alpha.fn <- function(data, index) {
  # читаем переменные
  X = data$X[index]
  Y = data$Y[index]
  # считаем оценку параметра alpha по формуле
  (var(Y) - cov(X, Y)) / (var(X) + var(Y) - 2*cov(X, Y))
}

# рассчитать alpha по всем 100 наблюдениям
alpha.100.obs <- alpha.fn(Portfolio, 1:100)
alpha.100.obs

## [1] 0.5758321

# создать бутстреп-выборку и повторно вычислить alpha
set.seed(my.seed)
alpha.boot.1 <- alpha.fn(Portfolio, sample(100, 100, replace = T))
alpha.boot.1

## [1] 0.7368375

# теперь -- многократное повторение предыдущей операции
alpha.boot.1000 <- boot(Portfolio, alpha.fn, R = 1000)
alpha.boot.1000

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Portfolio, statistic = alpha.fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##      original       bias    std. error
## t1* 0.5758321 -0.001695873  0.09366347

Бутстреп повторяет расчёт параметра много раз, делая повторные выборки из наших 100 наблюдений. В итоге этим методом можно вычислить стандартную ошибку параметра, не опираясь на допущения о законе его распределения. В нашем случае \(\alpha = 0.576\) со стандартной ошибкой \(s_{\hat{\alpha}} = 0.094\).

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

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

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

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

# функция для расчёта коэффициентов ПЛР по выборке из данных
boot.fn <- function(data, index){
  coef(lm(mpg ~ horsepower, data = data, subset = index))
}
boot.fn(DF.auto, 1:n)

## (Intercept)  horsepower 
##  39.9358610  -0.1578447

# пример применения функции к бутстреп-выборке
set.seed(my.seed)
boot.fn(DF.auto, sample(n, n, replace = T))

## (Intercept)  horsepower 
##  40.3404517  -0.1634868

# применяем функцию boot для вычисления стандартных ошибок параметров
#  (1000 выборок с повторами)
boot(DF.auto, boot.fn, 1000)

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = DF.auto, statistic = boot.fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##       original        bias    std. error
## t1* 39.9358610  0.0549915227 0.841925746
## t2* -0.1578447 -0.0006210818 0.007348956

# сравним с ошибками параметров по МНК
summary(fit.lm.1)$coef

##               Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 41.2835482 1.04435157  39.53032 9.226000e-95
## horsepower  -0.1696587 0.00955597 -17.75421 1.613925e-42

# график остатков модели
plot(fit.lm.1, 3)

# оценки отличаются из-за того, что МНК -- параметрический метод с допущениями,
#  и, судя по графику остатков, допущения не выполняются

# вычислим оценки параметров квадратичной модели регрессии
boot.fn.2 <- function(data, index){
  coef(lm(mpg ~ poly(horsepower, 2), data = data, subset = index))
}
# применим функцию к 1000 бутсреп-выборкам
set.seed(my.seed)
boot(DF.auto, boot.fn.2, 1000)

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = DF.auto, statistic = boot.fn.2, R = 1000)
## 
## 
## Bootstrap Statistics :
##       original       bias    std. error
## t1*   23.44592 -0.003660358   0.2195369
## t2* -120.13774  0.002769239   3.6138046
## t3*   44.08953  0.101767465   4.1998076

# сравним с ошибками параметров по МНК
summary(fit.lm.2)$coef

##                        Estimate Std. Error    t value      Pr(>|t|)
## (Intercept)            23.54955  0.3174554  74.182235 8.123709e-144
## poly(horsepower, 2)1 -123.58813  6.4587373 -19.135029  1.858719e-46
## poly(horsepower, 2)2   47.71889  6.3612782   7.501463  2.247784e-12

# график остатков модели
plot(fit.lm.2, 3)

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

Упражнение 5

1 Оценить стандартную ошибку модели для линейных регрессионных моделей из упражнения 4 (варианты ниже): а) со всеми объясняющими переменными; б) только с непрерывными объясняющими переменными:

• методом проверочной выборки с долей обучающей 50%;

• методом LOOCV;

• k-кратной кросс-валидацией с \(k = 5\) и \(k = 10\).

Выбрать лучшую модель по минимуму ошибки. Все ли методы кросс-валидации сходятся на одной и той же модели?

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

Как сдавать: прислать на почту преподавателя ссылки:

• на html-отчёт с видимыми блоками кода (блоки кода с параметром echo = T), размещённый на rpubs.com.

• на код, генерирующий отчёт, в репозитории на github.com.

В текст отчёта включить постановку задачи и ответы на вопросы задания.

Варианты

См. варианты из упражнения 4.

Источники

1. James G., Witten D., Hastie T. and Tibshirani R. An Introduction to Statistical Learning with Applications in R. URL: http://www-bcf.usc.edu/~gareth/ISL/ISLR%20First%20Printing.pdf