Реузроковање

Методе реузорковања су један од основних алата у модерној статистици. Рецимо, на примјер, да имамо проблем линеарне регресије и да желимо да оцијенимо регресионе коефицијенте. Сходно том циљу, подијелили смо наше податке на тренинг и на тест скуп, како бисмо на тренинг скупу добили саму оцјену, а на тест скупу провјерили грешку.

Уколико радимо са подацима великог обима, ово је углавном све што нам треба. Међутим, уколико обим нашег узорка није задовољавајући, или уколико из било ког разлога немамо приступ тест скупу, мора се прибјећи провјери квалитета на самом тренинг скупу.

Такође, у таквим ситуацијама може се десити да је тренинг скуп нерепрезентативан, то јест да смо, случајно, подијелили податке на тренинг и тест скуп на тај начин да су оцјене коефицијената лоше. Стога се прибјегава подјели самог тренинг скупа на два дијела, где ћемо на једном правити модел, а на другом, који ћемо звати валидациони скуп, провјеравати квалитет модела. О чему причамо размотрићемо на примјеру.

Примјер 1. Посматрајмо \(Auto\) базу из пакета \(ISLR2\). База садржи 392 опсервације (аутомобила) за 9 промјенљивих. Рецимо да желимо да моделујемо број миља за галон \(mpg\) на основу броја коњских снага \(horsepower\). Стандардно, прва ствар која пада на памет јесте да правимо линеаран модел, али ћемо прије тога подијелити податке на тренинг и тест скуп, тј. издвојићемо само тренинг скуп јер нам је он овде до интереса. Претвараћемо се да тест скуп немамо, већ да њега чува налогодавац истраживања, како би за себе провјерио да ли наш модел стварно ваља.

set.seed(10)
library(ISLR2)

n <- length(Auto$mpg)
# Dijelimo 70:30
obim <- floor(0.7*n) 

train_idx <- sample(1:n, size = obim)
Auto_train <- Auto[train_idx, ]
# Pravimo model
model1 <- lm(mpg ~ horsepower, data = Auto_train)

Сада ћемо да извршимо предвиђање на основу овог модел.

predvidjanje1 <- predict(model1)

Сада ћемо да видимо колика нам је грешка. Користимо средњеквадратну грешку као показатељ.

mean((Auto_train$mpg - predvidjanje1)^2)
## [1] 23.60343

Хајде сада да исцртамо ову зависност.

plot(mpg ~ horsepower, data = Auto_train)

Овако са слике, није да је веза баш линеарна. Стога се намеће додавање члана \(mpg^2\) у регресију.

model2 <- lm(mpg ~ horsepower + I(horsepower^2), data = Auto_train)

Вршимо предвиђање на основу другог модела и рачунамо грешку.

predvidjanje2 <- predict(model2)
mean((Auto_train$mpg - predvidjanje2)^2)
## [1] 19.75487

У овом случају добијамо значајно мању средњеквадратну грешку. Отуда закључујемо да је бољи модел онај код којег смо додали и квадратни члан.

У реду, рећи ће свако ко је догурао довде. У чему је фора? Ово већ знамо! Е па фора је у томе што радимо са тренинг скупом. Можда се у ових наших 70 посто са којима радимо заломило баш оних 70 посто који дају обрнут резултат. Можда додавање квадратног члана чак и квари прецизност предвиђања, само су се нама подаци “потрефили” да нас наведу на погрешан закључак. Можда је грешка на нама недоступном тестном скупу другачија.

Ту на сцену ступа реузорковање. Дијелимо доступне нам податке на “нове тренинг” податке (зовимо их под-тренинг, да не дође до забуне) и на валидациони скуп. Модел правимо искључиво на под-тренинг подацима, а затим грешку провјеравамо на валидационом. То понављамо довољан број пута, те гледамо мијења ли се шта у оцјенама наших коефицијената од случаја до случаја. То се може урадити на много начина, а ми ћемо кренути од најнаивнијег до оних тежих.

Напомена: Може доћи до забуне у именима. Наиме, ми смо овдје рекли да је неки имагинарни наручилац пројекта задржао дио података за себе (тест), да би на њима провјерио ваља ли оно што ми радимо, а нама послао остатак (тренинг). Ми смо то симулирали издвајањем \(Auto\_train\). Затим из тог тренинга издвајамо валидациони скуп, углавном више пута. Неко податке које смо добили зове просто подаци, а оно што ми издвајамо зове тренинг и тест, те је онда тест скуп само специјалан случај валидационог, када дијелимо само једном. То је само језички проблем, а ми ћемо се овдје држати номенклатуре: [тренинг + валидациони] + тест.

Чисто да видимо шта се дешава са грешкама, одрадићемо мини демонстрацију. Подијелићемо на пола тренинг скуп, на једној половини правити модел, на другој валидирати, те гледати какве су грешке на обе половине, све поновивши 10 пута.

# Obim treninga i validacionog (pola)
n_tr <- length(Auto_train$mpg)

# Ovdje cemo da upisemo greske
MSE1 <- c()
MSE2 <- c()

for (i in 1:10) {
  # Biramo indekse za pod-trening
  tr_idx_tmp <- sample(1:n_tr, size = floor(n_tr / 2))
  
  # Izdvajamo validacioni
  Auto_tr_tmp <- Auto_train[tr_idx_tmp, ]
  Auto_vld_tmp <- Auto_train[-tr_idx_tmp, ]

  # Pravimo oba modela na Auto_tr_tmp
  model1_tmp <- lm(mpg ~ horsepower, data = Auto_tr_tmp)
  model2_tmp <- lm(mpg ~ horsepower + I(horsepower^2), data = Auto_tr_tmp)

  # Vrsimo predvidjanje na validacionom skupu za oba modela
  pred1_tmp <- predict(model1_tmp, newdata = Auto_vld_tmp)
  pred2_tmp <- predict(model2_tmp, newdata = Auto_vld_tmp)

  # Sada pravimo vektore gresaka za modele
  MSE1 <- c(MSE1, mean((Auto_vld_tmp$mpg - pred1_tmp)^2))
  MSE2 <- c(MSE2, mean((Auto_vld_tmp$mpg - pred2_tmp)^2))
}

Дакле, 10 пута смо половили тренинг скуп и на једној половини правили модел, а на другој валидирали. Сада ћемо грешке претворити у разлике: посматраћемо разлику грешака за модел \(2\) и модел \(1\).

MSE_diff <- MSE1 - MSE2
barplot(MSE_diff)

Видимо да се у једном случају, седмом, добија негативна разлика, то јест недодавање квадратног члана даје бољи резултат! Да смо податке које смо добили од послодавца подијелили баш на тај начин, извукли бисмо погрешан закључак! Тек након већег броја реузорковања, видимо да је у већини случајева додавање квадратног члана ипак бољи избор.

Напомена: Како у валидационом приступу правимо модел на мањем скупу, он ће да буде лошији (доказ: ЛСМ), те валидациони приступ генерално грешке приказује већим него што јесу.

Кросфалидација

Претходна прича о валидацији била је јасна и праволинијска, али ипак у неком смислу наивна. Види се да ту има простора за побољшање. Једно од таквих побољшања приче о валидацији јесте унакрсна валидација, позната још и као кросвалидација (енг. Cross Validation).

LOOCV

Први тип кросвалидације јесте валидација са избацивањем једног елемента (енг. Leave One Out Cross Validation). Малоприје смо имали једноставан приступ: скуп опсервација смо дијелили на два једнака дијела. Овај пут, ипак, урадићемо нешто сасвим другачије. Издвојићемо само једну опсервацију, на остатку направити модел, а на њој валидирати.

Нека смо издвојили опсервацију \((x_1, y_1)\), на остатку направили модел и срачунали грешку \(MSE_1 = (y_1-\hat{y}_1)^2\). Затим избацимо \((x_2, y_2)\), направимо модел на остатку и срачунамо \(MSE_2\). Настављамо даље све док не срачунамо и \(MSE_n\). \(LOOCV\) оцјена за грешку на цијелом скупу сада је дата са

\[CV_{(n)}=\frac{1}{n}\sum_{i=1}^n MSE_i\]

Овај приступ, за разлику од претходног (половљење), обично не прецјењује грешку. Такође, тамо смо вршили много мање подјела од могућег броја подјела, док смо овде сваку подјелу типа \(1 + (n − 1)\) узели у обзир. Сходно томе, овај приступ је за разлику од претходног слободан од случајности.

Дјелује да смо на добром путу. Међутим, појављује се један велики проблем. Модел се прави чак \(n\) пута! Ово је рачунарски веома захтјеван процес. Срећом математичари су ухватили пречицу. Ако дефинишемо

\[h_i = \frac{1}{n} + \frac{(x_i-\bar{x})^2}{\sum_{i=1}^n (x_i-\bar{x})^2},\]

за просту линеарну регресију важи да је:

\[CV_{(n)}=\frac{1}{n}\sum_{i=1}^n \left(\frac{y_i-\hat{y}_i}{1-h_i}\right)^2,\]

гдје су \(y_i\) оригиналне фитоване вриједности на цијелом тренинг скупу. Формула важи и за линеарну регресију која није проста, само се тада тежине \(h_i\) дефинишу мало другачије (требало би да буде на ЛСМ).

\(LOOCV\) може да буде коришћен не само са регресијом, већ и са разним типовима класификације. Нажалост, горња формула важи само у случају регресије, те у свим другим случајевима фитовање мора бити извршено \(n\) пута.

\(k-fold\) кросвалидација

Једна од алтернатива за \(LOOCV\) јесте \(k-fold\) кросвалидација. Идеја је да се скуп опсервација подијели на \(k\) дисјунктних скупова (\(folds\)), приближно исте величине. Затим се прва група третира као валидациони скуп, а модел се прави на осталих \(k − 1\) узетих заједно. Затим се на том првом рачуна грешка \(MSE_1\). Слично се понавља са осталима и добија се низ грешака \(MSE_1\), \(MSE_2\),…, \(MSE_k\), од којих се прави оцјена

\[CV_{(k)} = \frac{1}{k} \sum_{i=1}^k MSE_i.\]

Овдје се поставља ново питање: које \(k\) користити? Не постоји добар одговор. У суштини веће \(k\) даје боље резултате јер су “међумодели” бољи, али може да буде веома захтјевно за рачунар. Резултати за мало \(k\) су грубљи, али лакше израчунљиви. Свако треба да бира своје \(k\) спрам модела који користи, спрам обима узорка који има и спрам рачунара на којем ради.

Пристрасност vs дисперзија за \(k-fold\)

Када смо помињали однос \(LOOCV\) алгоритма и приступа пола-пола, поменули смо да \(LOOCV\) даје непристраснију оцјену средњеквадратне грешке на тестном скупу неголи метод пола-пола, јер потоњи модел прави само на половини података. Дакле, спрам непристрасности, \(LOOCV\) је свакако пожељнији метод. Из истог разлога \(LOOCV\) је пожељнији у односу на \(k-fold\).

Међутим, није пристрасност оцјене једина ствар која треба да нас брине. Наиме, испоставља се да \(LOOCV\) процедура има значајно већу дисперзију у односу на \(k-fold\). Шта то значи и откуда то? Наиме, када вршимо \(LOOCV\), ми упросјечавамо излазе \(n\) модела који су прављени на скоро идентичним подацима. Стога су ти излази међусобно јако позитивно корелисани. Са друге стране, ако примјењујемо \(𝑘-fold\), ми упросјечавамо \(𝑘 < 𝑛\) много мање корелисаних модела, јер су преклапања мања. Стога је оцјена коју даје \(𝑘-fold\) са мањом дисперзијом од \(LOOCV\).

Да сумирамо, \(LOOCV\) има мању пристрасност оцјене, али већу дисперзију, док је за \(𝑘-fold\) обрнуто.

Примјер 2. Сада ћемо се вратити примјеру са базом \(Auto\), и примјенити кросвалидациони приступ. Користићемо модел са квадратним чланом.

За почетак, нацртајмо грешке (не разлике - већ баш грешке) за приступ пола/пола:

barplot(MSE2)

Види се значајна варијабилност у грешкама, за различите подјеле. Одрадићемо и \(10-fold\), те упоредити варијабилност оцјене за \(MSE\).

У \(R\)-у постоји библиотека \(caret\) (енг. Classification and Regression Training), која је једна од најпознатијих библиотека за машинско учење које постоје за \(R\).

library(caret)

Сада ћемо да направимо линеарни модел у овом пакету и да извршимо кросвалидацију.

train_control_10f <- trainControl(method = "cv", number = 10)

Ово је објекат (листа), која заправо служи као упакован начин да се наредној функцији каже шта да ради.

model_10fold <- train(mpg ~ horsepower + I(horsepower^2), 
                      data = Auto_train, 
                      trControl = train_control_10f, 
                      method = "lm")

У повратној вриједности \(model\)resample\(RMSE\) чувају се коријени средњеквадратних грешака за сваки фолд. Одавде може да се добије кросвалидациона оцјена као

sum((model_10fold$resample$RMSE)^2)/10
## [1] 20.21792

Сада ћемо поновити поступак десет пута, и видјети како варирају оцјене.

cv_10_err <- c()

for (i in 1:10) {
  model_i <- train(mpg ~ horsepower + I(horsepower^2), 
                   data = Auto_train,
                   trControl = train_control_10f, 
                   method = "lm")
  cv_i <- sum((model_i$resample$RMSE)^2) / 10
  cv_10_err <- c(cv_10_err, cv_i)
}

barplot(cv_10_err)

Видимо много мању варијабилност оцјене \(MSE\) у односу на пола-пола приступ. Наравно, нема баш смисла радити ово исто за \(LOOCV\), јер бисмо свих десет пута добили исту вриједност, јер немамо случајан одабир фолдова.

Кросвалидација и класификациони проблеми

До сада смо кросвалидацију демонстрирали на проблему регресије гдје је зависна промјенљива квантитативно обиљежје. Међутим, кросвалидација је подједнако корисна и за проблеме класификације, гдје је зависна промјенљива категоричког типа.

На примјер, за \(LOOCV\) при класификацији имамо оцјену:

\[CV_{(n)} = \frac{1}{n}\sum_{i=1}^n ERR_i,\]

гдје је \(ERR_i = 𝐼\{y_i\neq \hat{y}_i\}\). Овдје, чак, \(CV_{(𝑛)}\) има више смисла него код регресије: он представља проценат погрешно класификованих опсервација.

Кросвалидациона оцена за \(𝑘-fold\) се дефинише аналогно.

Бутстреп

Код бутстреп приступа прескочићемо теоријски увод, те се одмах бацити на примјер - тако ће нам бити јасније. Рецимо да желимо да уложимо фиксну своту новца у два асета на берзи, који, редом, имају приносе \(𝑋\) и \(𝑌\), гдје су \(𝑋\) и \(𝑌\) случајне величине. Желимо да уложимо фракцију \(\alpha \in (0, 1)\) нашег новца у први асет, а \(1 − \alpha\) у други. Овакво улагање желимо да направимо тако да дисперзија нашег портфолија буде што мања. Другим речима, желимо да минимизујемо

\[D(\alpha X+(1-\alpha)Y).\]

Може се показати (није тема овог курса), да је оно \(\alpha\) које минимизује горњу дисперзију дато са:

\[\alpha = \frac{\sigma_Y^2-\sigma_{XY}}{\sigma_X^2+\sigma_Y^2-2\sigma_{XY}},\]

гдје је \(\sigma_X^2 = DX\), \(\sigma_Y^2 = DY\) и \(\sigma_{XY}=Cov(X,Y)\).

Нама (инвеститору) је циљ да се домогне овог \(\alpha\). Наизглед је све ту. Међутим, проблем је тај што су величине \(\sigma_X^2\), \(\sigma_Y^2\) и \(\sigma_{XY}\) у начелу непознате, те мора да се користи оцјена

\[\hat{\alpha} = \frac{\hat{\sigma}_Y^2-\hat{\sigma}_{XY}}{\hat{\sigma}_X^2+\hat{\sigma}_Y^2-2\hat{\sigma}_{XY}}.\]

Хајде да претпоставимо да је позната расподјела за \(𝑋\) и \(𝑌\), као и \((𝑋,𝑌)\). Тада можемо да урадимо сљедеће:

  1. Генеришемо узорак обима \(𝑛\) за обе расподјеле;
  2. Из тих узорака оцијенимо све неопходне параметре;
  3. Поновимо кораке 1. и 2. \(𝑁\) пута (велики број);
  4. За сваки од три параметра од интереса, назовимо га привремено \(𝑝\), добили смо низ оцјена \(𝑝_1\), \(𝑝_2\),…, \(𝑝_𝑁\). За ове оцјене цртамо хистограм, те узимамо просјек (или медијану, моду итд) за оцјену параметра.

Ове класичне Монте Карло симулације дају изванредне резултате у случајевима када знамо (макар приближно) расподјелу обиљежја од интереса. Оно што ми у пракси имамо јесте низ прошлих вриједности за \(𝑋\) и \(𝑌\). Стога ми морамо на основу доступних узорака, који су коначног обима, без могућности генерисања нових, да оцијенимо параметре.

Идеја бутстреп алгоритма је у томе да се из доступних узорака ваде узорци. И то не било какви, већ са понављањем, истог обима као почетни узорак, како бисмо добили што већи варијабилитет извучених узорака. Уколико имамо скуп података \(Z\), примјеном бутстреп узорковања \(R\) пута добићемо \(R\) узорака са понављањем истог обима као и \(Z\). Нека су то \(Z_1^*\), \(Z_2^*\),…, \(Z_R^*\). На основу њих правимо низ оцјена за \(\alpha\): \(\hat{\alpha}_1^*\), \(\hat{\alpha}_2^*\),…, \(\hat{\alpha}_R^*\). Оцјену за \(\alpha\) стандардно рачунамо као просјек свих ових оцјена, а оцјену дисперзије као узорачку дисперзију:

\[\widehat{D\hat{\alpha}} = \frac{1}{R-1}\sum_{i=1}^R\left(\hat{\alpha}_i^* - \frac{1}{R}\sum_{j=1}^R \hat{\alpha}_j^* \right)^2.\]

Да видимо како то сада да испрограмирамо у \(R\)-у.

Процес спровођења бутстрепа у \(R\)-у одвија се у два корака. Прво дефинишемо функцију која рачуна параметар (оцјену) од интереса. Затим користимо функцију \(boot\) из истоименог пакета која врши узроковање са понављањем и комбинује оцјене.

Примјер 3. У пакету \(ISLR2\) постоји база података \(Portfolio\) у којој су вјештачки генерисани подаци обима \(100\) који представљају принос два асета које смо малоприје поменули. За почетак правимо функцију која рачуна оцјену параметра. Њој прослијеђујемо базу, као и вектор индекса на основу којих се врши оцјењивање.

alfa_fja <- function(data, index) {
  X <- data$X[index]
  Y <- data$Y[index]
  (var(Y) - cov(X, Y))/(var(X) + var(Y) - 2*cov(X,Y))
}

Функција \(boot\) као аргумент прима базу података (\(data.frame\)), функцију за рачунање оцјене, као и аргумент \(R\) који говори о томе колико пута се врши узорковање.

library(boot)
boot(Portfolio, alfa_fja, R = 1000)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Portfolio, statistic = alfa_fja, R = 1000)
## 
## 
## Bootstrap Statistics :
##      original     bias    std. error
## t1* 0.5758321 0.00444577  0.09394233

Задаци за самосталан рад

Задатак 1. У бази \(BreastCancer\) (пакет \(mlbench\)) налазе се подаци о резултатима тумора који је класификован као бенигни или малигни. Направити класификатор по избору који погађа да ли је тумор бенигни или малигни користећи приступ унакрсне валидације.

Задатак 2. У пакету \(ISLR\) пронаћи базу \(Weekly\).

а) Направити логистички модел за промјенљиву \(Direction\) и предикторе \(Lag1\), \(Lag2\), \(Lag3\), \(Lag4\), \(Lag5\) и \(Volume\). Да ли су неки од тих предиктора значајни? Провјерити колико је добар тај модел примјењен на скупу за обучавање.

б) Издвојити скуп од 1990. до 2007. године као скуп за обучавање и спровести модел поново. Пронаћи најбољи могући логистички модел за овај случај (у смислу да најмање гријеши). Које предикторе садржи тај модел?

в) Поновити дио под (б) за разне комбинације предиктора и њихових степена само користећи приступ унакрсне валидације за \(k = 10\).

Задатак 3. Ручно написати функцију која ради \(k-fold\) валидацију.