Ex10 Shrinkage
Samuel Hudec
11 decembra 2018
Cvičenia na výber modelu načrtli neistú tématiku hľadania vhodného modelu. Štruktúra Lineárneho modelu nám ponúka variabilitu ako modelovať. Pri robustnej verzii sme využili medián, kvantil a podobne.
Naviažme na posledné cvičenie, kde sme spomínali Veľa a málo regresorov. Vtedy keď regresorov je toľko, že si každý môžeme poriadne prezrieť,sme trénovali posledných 9 cvičení. V prípade veľa regresorov sa situácia mení, nemôžeme si prezrieť regresory, ale môžeme využiť “triky” rôzne miery na redukciu.
Redukovať počet regresorov môžete rôzne. Ešte pred tým než sa pustíte do popisovania odozvy lineárnym modelom si vezmete prediktory a urobíte ortogonálnu transformáciu. Proces sa nazýva anaýza hlavných komponentov (PCA) je to jeden z využívaných a najelegantnejších postupov na redukciu regresorov a zároveň získate nové regresory, ktoré majú dobré vlastnosti.
Regresné modely obsahujúce penáltu (cenu) za počet prediktorov sú tiež veľmi účinné a aj keď by ste ich nepoužili, tak vám vedia dobre poradiť alebo prezradiť budúcu štruktúru lineárneho modelu.
Analýza hlavných komponentov (budúcich prediktorov)
Metóda sa hojne využíva v rôznych odvetviach štatistiky, nie len pre lineárnou regresiou a má aj niekoľko verzií. Na dátach hitters z knižnice ISLR si ukážeme redukciu parametrov na základe hlavných komponentov. Ako prvé načítajte dáta a využite doterajšie znalosti na ich obhliadku.
## AtBat Hits HmRun Runs
## Min. : 16.0 Min. : 1 Min. : 0.00 Min. : 0.00
## 1st Qu.:255.2 1st Qu.: 64 1st Qu.: 4.00 1st Qu.: 30.25
## Median :379.5 Median : 96 Median : 8.00 Median : 48.00
## Mean :380.9 Mean :101 Mean :10.77 Mean : 50.91
## 3rd Qu.:512.0 3rd Qu.:137 3rd Qu.:16.00 3rd Qu.: 69.00
## Max. :687.0 Max. :238 Max. :40.00 Max. :130.00
##
## RBI Walks Years CAtBat
## Min. : 0.00 Min. : 0.00 Min. : 1.000 Min. : 19.0
## 1st Qu.: 28.00 1st Qu.: 22.00 1st Qu.: 4.000 1st Qu.: 816.8
## Median : 44.00 Median : 35.00 Median : 6.000 Median : 1928.0
## Mean : 48.03 Mean : 38.74 Mean : 7.444 Mean : 2648.7
## 3rd Qu.: 64.75 3rd Qu.: 53.00 3rd Qu.:11.000 3rd Qu.: 3924.2
## Max. :121.00 Max. :105.00 Max. :24.000 Max. :14053.0
##
## CHits CHmRun CRuns CRBI
## Min. : 4.0 Min. : 0.00 Min. : 1.0 Min. : 0.00
## 1st Qu.: 209.0 1st Qu.: 14.00 1st Qu.: 100.2 1st Qu.: 88.75
## Median : 508.0 Median : 37.50 Median : 247.0 Median : 220.50
## Mean : 717.6 Mean : 69.49 Mean : 358.8 Mean : 330.12
## 3rd Qu.:1059.2 3rd Qu.: 90.00 3rd Qu.: 526.2 3rd Qu.: 426.25
## Max. :4256.0 Max. :548.00 Max. :2165.0 Max. :1659.00
##
## CWalks League Division PutOuts Assists
## Min. : 0.00 A:175 E:157 Min. : 0.0 Min. : 0.0
## 1st Qu.: 67.25 N:147 W:165 1st Qu.: 109.2 1st Qu.: 7.0
## Median : 170.50 Median : 212.0 Median : 39.5
## Mean : 260.24 Mean : 288.9 Mean :106.9
## 3rd Qu.: 339.25 3rd Qu.: 325.0 3rd Qu.:166.0
## Max. :1566.00 Max. :1378.0 Max. :492.0
##
## Errors Salary NewLeague
## Min. : 0.00 Min. : 67.5 A:176
## 1st Qu.: 3.00 1st Qu.: 190.0 N:146
## Median : 6.00 Median : 425.0
## Mean : 8.04 Mean : 535.9
## 3rd Qu.:11.00 3rd Qu.: 750.0
## Max. :32.00 Max. :2460.0
## NA's :59- Podarilo sa vám dáta zobraziť?
Použitie PCA má aj svoje nevýhody, vyžaduje plnú maticu bez NA a regresory môžu byť len v podobe scale. Na druhú stranu, faktorové regresory, ktoré poprípade nezahrniete môžete naspať vrátiť do matice po použití metódy. My vracať nič nebudeme a využijeme len nové špecialne regresory získané metódou.
- vyberte regresory na PCA
- očistite tieto regresory od
NA
Náhodne rozdelime riadky matice na trénovacie a testovacie v pomere 3:1 za použitia knižnice dplyr.
set.seed(1112)
Hitters = Hitters %>% mutate(name = rownames(Hitters))
train = Hitters %>% sample_frac(3/4, replace=FALSE)
test = anti_join(Hitters, train, by="name")
train = train %>% select(-name)
test = test %>% select(-name)- Prečítajte kód.
Na začiatok, zoberte trénovacie dáta a fitnite do nich čo najlepší lineárny fit so Salary ako odozvou. Spôsob akým to spravíte, je na vás. Ja som našiel napríklad takýto
##
## Call:
## lm(formula = Salary ~ AtBat + Hits + Walks + CRuns + CRBI + CWalks +
## PutOuts, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -797.54 -161.92 -21.13 121.79 1207.38
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.1377 67.3885 -0.091 0.927526
## AtBat -1.7744 0.5702 -3.112 0.002146 **
## Hits 6.5407 1.7624 3.711 0.000271 ***
## Walks 5.9301 1.5624 3.796 0.000198 ***
## CRuns 0.6134 0.2464 2.489 0.013669 *
## CRBI 0.4446 0.1955 2.274 0.024107 *
## CWalks -0.6777 0.2606 -2.600 0.010046 *
## PutOuts 0.3250 0.0769 4.226 3.7e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 287.6 on 189 degrees of freedom
## Multiple R-squared: 0.5212, Adjusted R-squared: 0.5035
## F-statistic: 29.39 on 7 and 189 DF, p-value: < 2.2e-16- Koľko percent variancie Váš alebo môj model opisuje?
Pokračujem využítím PCA, ale len na budúce prediktory. Na metódu využite funkciu prcomp(). Pozor, budeme robiť transformáciu a na to je vhodné aby matica prediktorov bola preškálovaná.
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.6149 2.1043 1.3601 0.94311 0.80350 0.7408 0.48868
## Proportion of Variance 0.4274 0.2767 0.1156 0.05559 0.04035 0.0343 0.01493
## Cumulative Proportion 0.4274 0.7041 0.8197 0.87533 0.91568 0.9500 0.96490
## PC8 PC9 PC10 PC11 PC12 PC13
## Standard deviation 0.41611 0.34297 0.31447 0.25533 0.24100 0.16440
## Proportion of Variance 0.01082 0.00735 0.00618 0.00407 0.00363 0.00169
## Cumulative Proportion 0.97572 0.98307 0.98925 0.99333 0.99696 0.99865
## PC14 PC15 PC16
## Standard deviation 0.12408 0.07029 0.03609
## Proportion of Variance 0.00096 0.00031 0.00008
## Cumulative Proportion 0.99961 0.99992 1.00000- Čo by sa stalo, keby sme maticu nepreškálovali?
- Podľa znalosti z prednážok interpretujte výstup.
Vypíšme si podiel variancií prvých dvoch komponentov pomocou $rot ,aby sme získali ako takú predstavu o nich.
## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits
## -0.14 -0.14 -0.17 -0.15 -0.19 -0.18 -0.31 -0.35 -0.35
## CHmRun CRuns CRBI CWalks PutOuts Assists Errors
## -0.34 -0.36 -0.36 -0.34 -0.05 0.02 0.04## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits
## 0.41 0.40 0.28 0.41 0.36 0.26 -0.22 -0.15 -0.14
## CHmRun CRuns CRBI CWalks PutOuts Assists Errors
## -0.08 -0.13 -0.12 -0.15 0.18 0.14 0.16- Vysvetlite čo znamenajú čísla pri jednotlivých regresoroch.
- Ako by ste interpretovali tieto nové regresory?
Na voľbu ideálneho počtu komponentov (regresorov) použime napríklad obrázok.
- Vyberme zvolený počet komoponentov (ako na to zistite z help)
Od teraz naše nové regresory reprezentujú poctane niečo iné ako tie pôvodné. Priložme odozvu k regresorom a znazornime si pairs graf.
- Sme na tom teraz lepšie?
Určite by sa zišlo urobiť poriadnu diagnostiku a dáta odmodelovať. Pre krátkosť času tento krok preskočíme a len slepo fitneme plný regresný model.
##
## Call:
## lm(formula = Salary ~ ., data = pca_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -866.02 -172.01 -33.83 123.48 1217.69
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 520.452 21.583 24.114 < 2e-16 ***
## PC1 -95.927 8.275 -11.593 < 2e-16 ***
## PC2 43.747 10.283 4.254 3.28e-05 ***
## PC3 3.145 15.909 0.198 0.843505
## PC4 -76.713 22.943 -3.344 0.000995 ***
## PC5 -28.906 26.929 -1.073 0.284443
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 302.9 on 191 degrees of freedom
## Multiple R-squared: 0.4633, Adjusted R-squared: 0.4492
## F-statistic: 32.97 on 5 and 191 DF, p-value: < 2.2e-16- Porovnajte tento model s klasickým lineárnym zo začiatku cvičení.
Teraz nemusite nič robiť. Ideme sa pozrieť na predikčnú silu modelu na nových testovacích dátach.
Najprv porovnajme root mean square error medzi fitom lineárneho modelu a odozvou. Potom na ukážku porovnajme predikčnú silu modelu na testovacích dátach.
## [,1]
## RMSE Modelu 281.7247
## RMSE testovacich dat 407.3190Tak isto aj pre lineárny model, kde prediktory boli zložené z hlavných komponentov.
## [,1]
## RMSE Modelu 298.2790
## RMSE testovacich dat 436.9174Plastickou redukciou parametrov sme dosiahli skoro rovnaký model z hladiska RMSE. Pripomenme, že pri druhom modeli sme nepoužili žiadny model seleciton, nastroj. Poznámka: funkciou pcr() fitnete pca regresiu jedním ťahom.
Scvrkávacie regresné metódy
Prejdime k metódam nevyžadujúcim žiadne transformácie regresorov, prvá z nich je Ridge regression (hrebeňová regresia). V tomto prípade vieme odvodiť aj normálne rovnice: \[\hat\beta = (X'X + \lambda I)^{-1} X'y\]
Lineárny model počas odhadovania parametrov kladie penáltu (cenu) na prediktory. Niektoré možno “pošle” k nule a urobí za vás model selection.
Na ukážku metódy vezmime rovnaké dáta a fitnime do trénovacich dát hrebeňovú regresiu. Na model použite funkciu lm.ridge z knižnice MASS. Hyper parameter \(\lambda\) nech je tentokrát vektor seq(0, 100, len=21), aby sme mohli sledovať vývoj odhadov prarametrov pri meniacom sa hyper parametri.
Zobrazme si odhady parametrov pri jednotlivých lambdách pomocou príkazu:
matplot(rmo$lambda, coef(rmo), type="l", ylim=c(-10,10),
xlab=expression(lambda),ylab=expression(hat(beta)),col=1)
Tá čiara idúca dole je intercept.
- Chápete prečo sa volá hrebeňová?
- Prečo je parameter \(\lambda\) tak neumerne zvýšený?
Na zistenie optimalného paramera \(\lambda\) použijme zovšeobecnenú “krížovúvalidáciu” (detailnejšie bude v naväzujúcom kurze) a to za pomoci $GCV. Na základe toho fitnime do dát model a zmerajme rovnaké hodnoty RMSE.
## [,1]
## RMSE Modelu 288.6547
## RMSE testovacich dat 428.5591- Ako hodnotite výsledok?
Hodnoty sú vlastne rovnaké aj bez toho, aby sme použili niekedy zdĺhavé modelovanie.
Pokračujme s ďalšou scvrkávacou metódou v skratke LASSO (Least absolute shrinkage and selection operator). Tu už nevieme spočítať normálne rovnice, ale parametre sa dajú počítať inak. Algoritmus na odhad prametrov má veľkú výhodu, robí v jednom kroku odhadovnie prametrov a model selection!
Funkcií implementovaných v Rku je viacej my použime lars z rovnomennej knižnice. Znova fitnite model, pozor sme v oblasti machine learning (ML), tu sa už model zadáva inak ako pri linearnej regresii a jej podobných. Zároveň sa ťažko počíta inferencia.
Na zistenie vhodného parametra modelu znova využime krížovú validáciu s funkciou cv.lars(). Funkcia priamo ponúkne obrázok a stačí použiť oko.
Znova urobme porovnanie predikčnej sily modelu na zaklade RMSE.
## [,1]
## RMSE Modelu 294.9951
## RMSE testovacich dat 435.2773Dostávame v hrubom zmysle rovnaké výsledky. Dokonca sa javý lineárny model ako najlepší z hladiska predikčnej sily meranej pomocou RMSE. Predstavme si prípad , kde regresorov bude ďaleko viacej ako máme teraz. Následne by voľba najvhodnejšieho modelu a zároveň odhadu bola zložitejšia s klasickým lineárnym modelom. Všetko má svoje pre aj proti, LASSO navyše nevyžaduje také prísne predpoklady.
S penáltov kladenou na parametre sa dá ďalej “hrať”, preto vznikli rôzne kombinácie a verzie týchto metód.
Samostatná práca
Vezmite dáta z predchádzajúcich cvičení seatpos, rozdelte dataset na trénovaciu a testovaciu vzorku a fitnite do dát podľa vás optimálny model s hipcenter ako odozvou. Potom využite nasledovné:
- PCA regresiu
- Ridge regresiu
- LASSO
Na nájdenie najvhodnejšieho modelu na základe RMSE a porovnajte s lineárnym modelom.