Ridge and Lasso
AN JUNGCHAN
2021 9 30
DATA set
데이터 불러와서 보기
library(ISLR)
library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.4 v dplyr 1.0.7
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 2.0.1 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()data(Hitters)
names(Hitters)## [1] "AtBat" "Hits" "HmRun" "Runs" "RBI" "Walks"
## [7] "Years" "CAtBat" "CHits" "CHmRun" "CRuns" "CRBI"
## [13] "CWalks" "League" "Division" "PutOuts" "Assists" "Errors"
## [19] "Salary" "NewLeague"head(Hitters)## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun
## -Andy Allanson 293 66 1 30 29 14 1 293 66 1
## -Alan Ashby 315 81 7 24 38 39 14 3449 835 69
## -Alvin Davis 479 130 18 66 72 76 3 1624 457 63
## -Andre Dawson 496 141 20 65 78 37 11 5628 1575 225
## -Andres Galarraga 321 87 10 39 42 30 2 396 101 12
## -Alfredo Griffin 594 169 4 74 51 35 11 4408 1133 19
## CRuns CRBI CWalks League Division PutOuts Assists Errors
## -Andy Allanson 30 29 14 A E 446 33 20
## -Alan Ashby 321 414 375 N W 632 43 10
## -Alvin Davis 224 266 263 A W 880 82 14
## -Andre Dawson 828 838 354 N E 200 11 3
## -Andres Galarraga 48 46 33 N E 805 40 4
## -Alfredo Griffin 501 336 194 A W 282 421 25
## Salary NewLeague
## -Andy Allanson NA A
## -Alan Ashby 475.0 N
## -Alvin Davis 480.0 A
## -Andre Dawson 500.0 N
## -Andres Galarraga 91.5 N
## -Alfredo Griffin 750.0 ANULL 값 관리
Hitters %>%
is.na() %>%
sum()## [1] 59NA값의 갯수를 확인한다.
Hitters데이터 중에서 na를 가지고 있는 관측치를 TRUE, FALSE로 식별해준 후,
다 더해서 NA값을 확인한다.
Hitters %>%
na.omit() -> HittersHitters데이터에서 NA값을 없애버린다.
그리고 그러한 데이터 셋을 Hitters에 다시 저장한다.
Ridge
개념
- RSS를 최소화 시키는데, L2-norm(제곱의 절댓값)을 제약조건으로 둔 regression model 이다.
- lamda가 작으면, b에 대한 패널티가 감소하여, b값이 커질 수가 있다.
- lamda가 크면, b에 대한 패널티가 증가하여, b의 값이 제약을 받으면서 축소된다.
- lamda값을 어떤 것으로 정하느냐에 따라서 회귀 모델이 달라지므로, lamda를 tunning하는 것이 중요하다.
lamda가 크면, 축소를 많이 시킨 것이고, lamda가 작으면 축소를 안시킨 것이다.
패키지
#install.packages("glmnet")
library(glmnet)## 필요한 패키지를 로딩중입니다: Matrix##
## 다음의 패키지를 부착합니다: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## Loaded glmnet 4.1-2분석
- 데이터 나누기
우선, 회귀분석은 x를 가지고 y를 추측하는 것이다. 우리는 타자들의 연봉을 추청하고 싶은 것이므로, Salary를 제외한 나머지의 값을 x에 저장하고, 우리가 추측하고 싶은 변수는 y에 저장해야 할 것이다. 우선 데이터의 형태를 보자.
str(Hitters)## 'data.frame': 263 obs. of 20 variables:
## $ AtBat : int 315 479 496 321 594 185 298 323 401 574 ...
## $ Hits : int 81 130 141 87 169 37 73 81 92 159 ...
## $ HmRun : int 7 18 20 10 4 1 0 6 17 21 ...
## $ Runs : int 24 66 65 39 74 23 24 26 49 107 ...
## $ RBI : int 38 72 78 42 51 8 24 32 66 75 ...
## $ Walks : int 39 76 37 30 35 21 7 8 65 59 ...
## $ Years : int 14 3 11 2 11 2 3 2 13 10 ...
## $ CAtBat : int 3449 1624 5628 396 4408 214 509 341 5206 4631 ...
## $ CHits : int 835 457 1575 101 1133 42 108 86 1332 1300 ...
## $ CHmRun : int 69 63 225 12 19 1 0 6 253 90 ...
## $ CRuns : int 321 224 828 48 501 30 41 32 784 702 ...
## $ CRBI : int 414 266 838 46 336 9 37 34 890 504 ...
## $ CWalks : int 375 263 354 33 194 24 12 8 866 488 ...
## $ League : Factor w/ 2 levels "A","N": 2 1 2 2 1 2 1 2 1 1 ...
## $ Division : Factor w/ 2 levels "E","W": 2 2 1 1 2 1 2 2 1 1 ...
## $ PutOuts : int 632 880 200 805 282 76 121 143 0 238 ...
## $ Assists : int 43 82 11 40 421 127 283 290 0 445 ...
## $ Errors : int 10 14 3 4 25 7 9 19 0 22 ...
## $ Salary : num 475 480 500 91.5 750 ...
## $ NewLeague: Factor w/ 2 levels "A","N": 2 1 2 2 1 1 1 2 1 1 ...
## - attr(*, "na.action")= 'omit' Named int [1:59] 1 16 19 23 31 33 37 39 40 42 ...
## ..- attr(*, "names")= chr [1:59] "-Andy Allanson" "-Billy Beane" "-Bruce Bochte" "-Bob Boone" ...데이터를 보았을 때, League와 Division, NewLeague는 Factor형 변수를 가진다는 것을 볼 수 있다.
우리는 이러한 factor형 변수들을 수치형 데이터로 바꿔서 문제를 해결해야 한다.
이럴때 사용하는 함수가 바로 model.matrix()이다.
model.matirx(예측하고자 하는 변수~.,데이터셋)[,-1]
여기서 [,-1]은 intercept값은 제외시키겠다는 이야기이다.
x<-model.matrix(Salary~.,Hitters)[,-1]
x %>% head()## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun
## -Alan Ashby 315 81 7 24 38 39 14 3449 835 69
## -Alvin Davis 479 130 18 66 72 76 3 1624 457 63
## -Andre Dawson 496 141 20 65 78 37 11 5628 1575 225
## -Andres Galarraga 321 87 10 39 42 30 2 396 101 12
## -Alfredo Griffin 594 169 4 74 51 35 11 4408 1133 19
## -Al Newman 185 37 1 23 8 21 2 214 42 1
## CRuns CRBI CWalks LeagueN DivisionW PutOuts Assists Errors
## -Alan Ashby 321 414 375 1 1 632 43 10
## -Alvin Davis 224 266 263 0 1 880 82 14
## -Andre Dawson 828 838 354 1 0 200 11 3
## -Andres Galarraga 48 46 33 1 0 805 40 4
## -Alfredo Griffin 501 336 194 0 1 282 421 25
## -Al Newman 30 9 24 1 0 76 127 7
## NewLeagueN
## -Alan Ashby 1
## -Alvin Davis 0
## -Andre Dawson 1
## -Andres Galarraga 1
## -Alfredo Griffin 0
## -Al Newman 0데이터가 수치형으로 잘 변환되었다는 것을 알 수 있다.
y에는 당연히 우리가 예측하려고 하는 변수인 y를 넣는다.
y<-Hitters$Salary
y %>% head()## [1] 475.0 480.0 500.0 91.5 750.0 70.02.glmnet()함수
ridge와 Lasso를 fit해주는 함수이다.
glmnet(관측할 값, 예측할 값, alpha=, lamda=)
alpha=0 : ridge 적합
alpha=1 : Lasso 적합
lamda=람다 값
lamda에 아무것도 쓰지 않으면 그냥 알아서 lamda를 설정해준다.
grid<-10^seq(10,-2,length=100)
grid## [1] 1.000000e+10 7.564633e+09 5.722368e+09 4.328761e+09 3.274549e+09
## [6] 2.477076e+09 1.873817e+09 1.417474e+09 1.072267e+09 8.111308e+08
## [11] 6.135907e+08 4.641589e+08 3.511192e+08 2.656088e+08 2.009233e+08
## [16] 1.519911e+08 1.149757e+08 8.697490e+07 6.579332e+07 4.977024e+07
## [21] 3.764936e+07 2.848036e+07 2.154435e+07 1.629751e+07 1.232847e+07
## [26] 9.326033e+06 7.054802e+06 5.336699e+06 4.037017e+06 3.053856e+06
## [31] 2.310130e+06 1.747528e+06 1.321941e+06 1.000000e+06 7.564633e+05
## [36] 5.722368e+05 4.328761e+05 3.274549e+05 2.477076e+05 1.873817e+05
## [41] 1.417474e+05 1.072267e+05 8.111308e+04 6.135907e+04 4.641589e+04
## [46] 3.511192e+04 2.656088e+04 2.009233e+04 1.519911e+04 1.149757e+04
## [51] 8.697490e+03 6.579332e+03 4.977024e+03 3.764936e+03 2.848036e+03
## [56] 2.154435e+03 1.629751e+03 1.232847e+03 9.326033e+02 7.054802e+02
## [61] 5.336699e+02 4.037017e+02 3.053856e+02 2.310130e+02 1.747528e+02
## [66] 1.321941e+02 1.000000e+02 7.564633e+01 5.722368e+01 4.328761e+01
## [71] 3.274549e+01 2.477076e+01 1.873817e+01 1.417474e+01 1.072267e+01
## [76] 8.111308e+00 6.135907e+00 4.641589e+00 3.511192e+00 2.656088e+00
## [81] 2.009233e+00 1.519911e+00 1.149757e+00 8.697490e-01 6.579332e-01
## [86] 4.977024e-01 3.764936e-01 2.848036e-01 2.154435e-01 1.629751e-01
## [91] 1.232847e-01 9.326033e-02 7.054802e-02 5.336699e-02 4.037017e-02
## [96] 3.053856e-02 2.310130e-02 1.747528e-02 1.321941e-02 1.000000e-02우선 그리드를 자동으로 찾아보지 않고 내가 지정한다. 이 때, 그리드는 \(10^10\) 부터 \(10^-2\) 까지의 값을 100개의 간격으로 자른 수이다.
ridge.mod<-glmnet(x,y,alpha=0,lamda=grid)그 후, 능형회귀에 lamda를 grid로 설정해주고 fitting해주었다. grid가 100개니까 100개의 모델이 fit되었을 것이다.
names(ridge.mod)## [1] "a0" "beta" "df" "dim" "lambda" "dev.ratio"
## [7] "nulldev" "npasses" "jerr" "offset" "call" "nobs"이 때, ridge.mod가 가진 names들을 보니까, 위와같이 있다는 것을 알았다.
ridge.mod %>%
coef() %>%
dim()## [1] 20 100이 때, ridge.mod의 회귀계수를 보고, dim을 봐서, 몇개의 model이 fitting 되었는지 보자.
값은 행이 20개고 열이 100개인 차원이 나온다.
즉, 값은 행에 각각의 변수들이 저장되고, 열에는 lamda의 값들이 저장되는 것을 확인할 수 있다.
100개의 lamda에 대한 각각의 20개의 회귀변수가 나왔다고 해석하면 될 것이다.
ridge.mod %>%
coef ->ridge.mod.coef
ridge.mod.coef[,50]## (Intercept) AtBat Hits HmRun Runs
## 213.066444060 0.090095728 0.371252755 1.180126954 0.596298285
## RBI Walks Years CAtBat CHits
## 0.594502389 0.772525465 2.473494235 0.007597952 0.029272172
## CHmRun CRuns CRBI CWalks LeagueN
## 0.217335715 0.058705097 0.060722036 0.058698830 3.276567808
## DivisionW PutOuts Assists Errors NewLeagueN
## -21.889942546 0.052667119 0.007463678 -0.145121335 2.972759111ridge.mod$lambda[50]## [1] 2674.375이 때, 50번째 lamda에 대한 회귀계수를 한번 보자.
ridge.mod.coef[,10]## (Intercept) AtBat Hits HmRun Runs
## 5.193354e+02 4.787004e-03 1.742001e-02 6.974621e-02 2.941856e-02
## RBI Walks Years CAtBat CHits
## 3.102036e-02 3.661001e-02 1.486101e-01 4.102827e-04 1.511717e-03
## CHmRun CRuns CRBI CWalks LeagueN
## 1.139588e-02 3.032875e-03 3.130119e-03 3.304540e-03 -4.416363e-02
## DivisionW PutOuts Assists Errors NewLeagueN
## -6.992187e-01 1.940563e-03 3.162091e-04 -1.571121e-03 -4.618317e-03ridge.mod$lambda[10]## [1] 110505.5이 때, 10번째 lamda에 대한 회귀계수도 같이 비교해보자.
보면 50번째 lamda가 더 lamda의 값이 작고, 회귀계수가 축소가 덜 된 것을 확인할 수 있다. 즉, 람다의 값이 작을수록 회귀계수의 축소가 덜 일어난다는 것을 알 수 있다.
반면, 10번째 람다에 대해서는 람다의 값이 훨씬 더 크게 형성되었다는 것을 알 수 있고, 이때 회귀계수돌은 매우 많이 축소되었다는 것을 알 수 있다.
L2-norm과 같이 비교해보자.
ridge.mod.coef[-1,50]^2 %>%
sum() %>%
sqrt->model_50
model_50## [1] 22.53415ridge.mod.coef[-1,10]^2 %>%
sum() %>%
sqrt->model_10
model_10## [1] 0.7221433이 때, model10이 L2-norm의 값이 더 작게 나오고 있다. 즉, 람다의 값이 더 크면, b가 축소되어, b와 관련된 L2-norm의 값이 작아진다. L2-norm의 값이 작을수록, 축소가 많이 일어났다는 것을 말해준다.
plot(ridge.mod, "lambda")
이 때, 람다에 대한 각 회귀계수의 값을 그래프로 그려보았는데, 람다가 작아질수록 회귀계수들이 수축이 안된다는 사실을 알 수 있다.
LASSO
개념
- ridge regression의 단점을 없앤 기법
- ridge regression은 변수를 제거할 수 없다. 절대로 b가 0으로 축소되지 않기 때문이다.(L2-norm이기 때문에!)
- lasso regression은 L1-norm을 사용하기 때문에, b가 0으로 축소될 수 있어서, ridge regression의 단점을 없애버렸다.
- 그렇기 때문에, 자유도가 항상 변수의 개수와 일치하는 ridge regression과 다르게, lasso는 자유도가 축소될 수 있다. (아예 b를 0으로 보낼 수 있기 때문이다!)
lamda가 크면, 축소를 많이 시킨 것이고, lamda가 작으면 축소를 안시킨 것이다.
분석
lasso.mod<-glmnet(x,y,alpha=1)alpha에 1을 주게 되면, lasso를 실시한다. 람다에 아무값도 주지 않는다면, 스스로 람다를 정해준다.
names(lasso.mod)## [1] "a0" "beta" "df" "dim" "lambda" "dev.ratio"
## [7] "nulldev" "npasses" "jerr" "offset" "call" "nobs"lasso.mod의 구성요소들을 볼 수 있다.
lasso.mod$lambda %>%
length()## [1] 80우선 임의로 람다를 80개를 뽑았다는 것을 알 수 있다.
lasso.mod$lambda## [1] 255.2820965 232.6035387 211.9396814 193.1115442 175.9560469 160.3245967
## [7] 146.0818014 133.1042968 121.2796779 110.5055256 100.6885193 91.7436287
## [13] 83.5933776 76.1671724 69.4006907 63.2353246 57.6176727 52.4990774
## [19] 47.8352041 43.5856564 39.7136268 36.1855777 32.9709507 30.0419023
## [25] 27.3730625 24.9413151 22.7255974 20.7067180 18.8671902 17.1910810
## [31] 15.6638728 14.2723375 13.0044224 11.8491453 10.7964999 9.8373686
## [37] 8.9634439 8.1671562 7.4416086 6.7805166 6.1781542 5.6293040
## [43] 5.1292122 4.6735471 4.2583620 3.8800609 3.5353670 3.2212947
## [49] 2.9351238 2.6743755 2.4367913 2.2203135 2.0230670 1.8433433
## [55] 1.6795857 1.5303760 1.3944216 1.2705450 1.1576733 1.0548288
## [61] 0.9611207 0.8757374 0.7979393 0.7270526 0.6624632 0.6036118
## [67] 0.5499886 0.5011291 0.4566102 0.4160462 0.3790858 0.3454089
## [73] 0.3147237 0.2867645 0.2612891 0.2380769 0.2169268 0.1976557
## [79] 0.1800965 0.1640972선택한 람다들을 볼 수 있다.
lasso.mod %>%
coef() %>%
dim()## [1] 20 80이 때, 람다 80개 각각에 대해서, \(b_0\)부터 \(b_19\)까지 총 20개의 회귀계수가 추정되었다는 사실을 알 수 있다.
lasso_lambda_df<-cbind(lasso.mod$lambda,lasso.mod$df)
colnames(lasso_lambda_df)<-c("lambda","df")
head(lasso_lambda_df,5)## lambda df
## [1,] 255.2821 0
## [2,] 232.6035 1
## [3,] 211.9397 2
## [4,] 193.1115 2
## [5,] 175.9560 3tail(lasso_lambda_df,5)## lambda df
## [76,] 0.2380769 18
## [77,] 0.2169268 18
## [78,] 0.1976557 18
## [79,] 0.1800965 18
## [80,] 0.1640972 19이 때, lambda값이 작아질수록 df의 값이 커진다는 사실을 알 수 있다. 즉, 람다의 값이 작으면 축소가 작게된 것이고, 람다의 값이 크면, 축소가 크게되었다는 사실을 알 수 있다.
dim(lasso.mod$beta)## [1] 19 80apply(lasso.mod$beta, 2, function(t) sum(t!=0))## s0 s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18 s19
## 0 1 2 2 3 4 4 4 4 4 5 5 5 5 6 6 6 6 6 6
## s20 s21 s22 s23 s24 s25 s26 s27 s28 s29 s30 s31 s32 s33 s34 s35 s36 s37 s38 s39
## 6 6 6 6 6 6 6 6 6 7 7 7 9 9 9 9 9 11 11 12
## s40 s41 s42 s43 s44 s45 s46 s47 s48 s49 s50 s51 s52 s53 s54 s55 s56 s57 s58 s59
## 12 13 13 13 13 13 13 13 13 13 13 14 15 15 17 17 17 17 17 17
## s60 s61 s62 s63 s64 s65 s66 s67 s68 s69 s70 s71 s72 s73 s74 s75 s76 s77 s78 s79
## 17 17 17 17 18 18 18 17 18 18 18 18 18 18 18 18 18 18 18 19다음으로는, 0이 아닌 beta의 갯수를 살펴볼 것이다. beta의 dim을 보니까, 행이 19개고, 열이 80개라고 나와있다. 즉, 열은 lambda를 행은 beta를 나타낸다고 할 수 있다.
즉, 열 방향으로 0이 아닌것을 모두 더하면, 0이 아닌 beta의 갯수가 나올 것이다.
이를 apply함수를 사용하여 계산하였다.
이 때, t!=0는 논리연산자이기 때문에 TRUE, FALSE로 계산되어 갯수를 출력할 수 있다.
apply(데이터, 옵션, function(t) 함수(t))
옵션 = 1 : 행의 방향으로 모두 더한다
옵션 = 2 : 열의 방향으로 모두 더한다
apply(lasso.mod$beta, 2, function(t) sum(abs(t)))## s0 s1 s2 s3 s4 s5
## 0.00000000 0.07026614 0.13476614 0.19402153 0.32607217 0.66902681
## s6 s7 s8 s9 s10 s11
## 1.16062643 1.60917353 2.01773750 2.38998186 2.72609874 2.97402205
## s12 s13 s14 s15 s16 s17
## 3.20123049 3.40846344 11.64817218 22.98691627 33.32443344 42.73991182
## s18 s19 s20 s21 s22 s23
## 51.32307688 59.13982854 66.26049080 72.75317515 78.66537137 84.05609839
## s24 s25 s26 s27 s28 s29
## 88.96411117 93.43957427 97.51365487 101.22906065 104.61046950 108.52993494
## s30 s31 s32 s33 s34 s35
## 114.71476794 120.36084436 125.89147670 131.25050408 136.15014891 140.62777705
## s36 s37 s38 s39 s40 s41
## 144.71425308 148.18864814 151.59045005 154.77018794 157.86174328 160.79016820
## s42 s43 s44 s45 s46 s47
## 163.63677945 166.20598165 168.54928183 170.69298424 172.63820733 174.40982707
## s48 s49 s50 s51 s52 s53
## 176.02382743 177.49419703 178.84082943 180.15787370 181.56430483 182.27355302
## s54 s55 s56 s57 s58 s59
## 183.66972594 189.17522968 192.91656387 196.32739300 199.43820868 202.27621917
## s60 s61 s62 s63 s64 s65
## 204.85652984 207.20743988 209.35199945 211.29939556 213.11247799 215.55494811
## s66 s67 s68 s69 s70 s71
## 217.47724806 219.02184399 220.60118640 222.07917527 223.51628146 224.81215797
## s72 s73 s74 s75 s76 s77
## 226.03163842 227.15990943 228.20005790 229.15914906 230.03587822 230.83954274
## s78 s79
## 231.57502970 231.92115738다음은 L1-norm을 구한 것을 나타내었다. t를 절댓값으로 묶어준 다음에, apply함수를 사용하여 합을 계산하였다.
lasso.L1.norm<-apply(lasso.mod$beta, 2, function(t) sum(abs(t)))
names(lasso.L1.norm)<-log(lasso.mod$lambda)
lasso.L1.norm## 5.54236919451432 5.44933545346989 5.35630171242546 5.26326797138104
## 0.00000000 0.07026614 0.13476614 0.19402153
## 5.17023423033661 5.07720048929218 4.98416674824775 4.89113300720332
## 0.32607217 0.66902681 1.16062643 1.60917353
## 4.7980992661589 4.70506552511447 4.61203178407004 4.51899804302561
## 2.01773750 2.38998186 2.72609874 2.97402205
## 4.42596430198118 4.33293056093676 4.23989681989233 4.1468630788479
## 3.20123049 3.40846344 11.64817218 22.98691627
## 4.05382933780347 3.96079559675904 3.86776185571462 3.77472811467019
## 33.32443344 42.73991182 51.32307688 59.13982854
## 3.68169437362576 3.58866063258133 3.4956268915369 3.40259315049247
## 66.26049080 72.75317515 78.66537137 84.05609839
## 3.30955940944805 3.21652566840362 3.12349192735919 3.03045818631476
## 88.96411117 93.43957427 97.51365487 101.22906065
## 2.93742444527033 2.84439070422591 2.75135696318148 2.65832322213705
## 104.61046950 108.52993494 114.71476794 120.36084436
## 2.56528948109262 2.47225574004819 2.37922199900377 2.28618825795934
## 125.89147670 131.25050408 136.15014891 140.62777705
## 2.19315451691491 2.10012077587048 2.00708703482605 1.91405329378162
## 144.71425308 148.18864814 151.59045005 154.77018794
## 1.8210195527372 1.72798581169277 1.63495207064834 1.54191832960391
## 157.86174328 160.79016820 163.63677945 166.20598165
## 1.44888458855948 1.35585084751506 1.26281710647063 1.1697833654262
## 168.54928183 170.69298424 172.63820733 174.40982707
## 1.07674962438177 0.983715883337344 0.890682142292916 0.797648401248488
## 176.02382743 177.49419703 178.84082943 180.15787370
## 0.704614660204059 0.611580919159631 0.518547178115203 0.425513437070775
## 181.56430483 182.27355302 183.66972594 189.17522968
## 0.332479696026347 0.239445954981919 0.146412213937491 0.0533784728930625
## 192.91656387 196.32739300 199.43820868 202.27621917
## -0.0396552681513655 -0.132689009195794 -0.225722750240222 -0.31875649128465
## 204.85652984 207.20743988 209.35199945 211.29939556
## -0.411790232329078 -0.504823973373506 -0.597857714417934 -0.690891455462362
## 213.11247799 215.55494811 217.47724806 219.02184399
## -0.78392519650679 -0.876958937551219 -0.969992678595647 -1.06302641964007
## 220.60118640 222.07917527 223.51628146 224.81215797
## -1.1560601606845 -1.24909390172893 -1.34212764277336 -1.43516138381779
## 226.03163842 227.15990943 228.20005790 229.15914906
## -1.52819512486222 -1.62122886590664 -1.71426260695107 -1.8072963479955
## 230.03587822 230.83954274 231.57502970 231.92115738람다값과 비교해보면, 람다값이 클수록 L1 norm이 아주 크게 축소되었고, 람다가 작게 나올수록 축소가 별로 되지 않은 것을 확인할 수 있다.
plot(lasso.mod,"lambda", label= TRUE)
마찬가지로, lambda에 대한 회귀계수의 각각의 값들을 plot한 결과, 람다가 커질수록 회귀계수가 축소되고, 람다가 작을수록 회귀계수가 축소가 잘 안된다는 사실을 알 수 있다.
plot(lasso.mod, "norm", label=TRUE)
이번엔 L1 norm을 기준으로 본 결과, L1 norm이 커지면 커질수록 회귀모형이 더 커진다는 사실을 알 수 있는데 이는 당연하다. 왜냐하면, L1 norm은 회귀계수의 절댓값을 모두 더한 값이기 때문이다.
파라미터 튜닝
개념
- 어떤 람다가 가장 좋은 람다인가를 파악하는 기법이다.
- 가장 좋은 모델에 대한 모수를 찾는 과정을 파라미터 튜닝이라고 한다.
- Ridge와 Lasso에는 sigle split과 cross validation을 이용한 파라미터 튜닝이 있다.
Ridge
single split
set.seed(123)
train<-sample(1:nrow(x), nrow(x)/2)
test<- (-train)sigle split은 전체 데이터에 대해서, train set과 test set을 한번만 나누는 것을 말한다.
이 때, test set은 오직 검증만을 위해 사용되고, train set은 파라미터 튜닝을 위해 사용된다.
즉 test set안에 있는 x값은 전혀 사용되지 않고, 오직 y값만 사용된다
이 떄, train set과 test set은 각각 반반으로 지정하였다.
y.test<-y[test]이 때, test set은 y값만 사용할 것이므로, test set안에있는 y값은 미리 추출하여 따로 저장해주자.
1. train set에 대해서 모델 fit 해주기
grid<-10^seq(10,-2,length=100)
fit.tran<-glmnet(x[train,],y[train],alpha=0,lambda=grid)
dim(coef(fit.tran))## [1] 20 100모델을 피팅한 결과, 100개에 람다에 대해서 20개의 회귀계수가 추출된 것을 알 수 있다.
이제 우리는 best lambda를 구하고 회귀계수를 추출해야 한다.
2. 모델 명 설정해주기
tran.model.name<-0:(length(fit.tran$lambda)-1)
tran.model.name## [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
## [26] 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
## [51] 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74
## [76] 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99그 후, 모델의 이름을 설정해준다. 모델의 이름을 0부터 99까지 설정한 이유는, glmnet()의 고유한 특성 때문이다.
colnames(fit.tran$beta)## [1] "s0" "s1" "s2" "s3" "s4" "s5" "s6" "s7" "s8" "s9" "s10" "s11"
## [13] "s12" "s13" "s14" "s15" "s16" "s17" "s18" "s19" "s20" "s21" "s22" "s23"
## [25] "s24" "s25" "s26" "s27" "s28" "s29" "s30" "s31" "s32" "s33" "s34" "s35"
## [37] "s36" "s37" "s38" "s39" "s40" "s41" "s42" "s43" "s44" "s45" "s46" "s47"
## [49] "s48" "s49" "s50" "s51" "s52" "s53" "s54" "s55" "s56" "s57" "s58" "s59"
## [61] "s60" "s61" "s62" "s63" "s64" "s65" "s66" "s67" "s68" "s69" "s70" "s71"
## [73] "s72" "s73" "s74" "s75" "s76" "s77" "s78" "s79" "s80" "s81" "s82" "s83"
## [85] "s84" "s85" "s86" "s87" "s88" "s89" "s90" "s91" "s92" "s93" "s94" "s95"
## [97] "s96" "s97" "s98" "s99"glmnet으로 fit시킨 beta들의 값들은 s00의 값으로 0부터 저장되기 때문에, 이러한 과정을 거쳐야한다. 즉, 임의로 glmnet이 이름을 지은 것을 그대로 따라간다고 보면 된다.
Error<-NULL
for (i in 1:length(fit.tran$lambda)){
fit.tran.pred<-predict(fit.tran, tran.model.name=tran.model.name[i], newx=x[test,])
Error[i]<-mean((fit.tran.pred[,i] - y.test)^2)
}
Error## [1] 214949.5 214949.5 214949.5 214949.5 214949.4 214949.4 214949.3 214949.2
## [9] 214949.1 214948.9 214948.7 214948.5 214948.1 214947.6 214947.0 214946.1
## [17] 214945.0 214943.5 214941.6 214939.0 214935.6 214931.0 214925.1 214917.2
## [25] 214906.7 214892.9 214874.7 214850.6 214818.8 214776.8 214721.2 214647.9
## [33] 214551.1 214423.5 214255.2 214033.5 213741.8 213358.5 212856.1 212203.8
## [41] 211350.9 210245.7 208822.6 207004.7 204706.9 201841.1 198327.6 194111.7
## [49] 189186.6 183616.8 177558.9 171262.5 165047.6 159250.3 154159.5 149960.1
## [57] 146709.7 144345.4 142745.1 141762.3 141247.9 141090.6 141205.5 141524.3
## [65] 141990.6 142559.1 143187.1 143837.9 144479.6 145086.8 145657.9 146153.9
## [73] 146571.1 146915.8 147189.8 147410.8 147600.0 147788.0 148007.6 148282.7
## [81] 148631.8 149052.7 149513.9 149992.4 150470.4 150929.1 151349.3 151721.0
## [89] 152037.5 152301.8 152515.5 152687.3 152822.9 152928.8 153009.4 153070.3
## [97] 153117.8 153152.2 153181.1 153200.9각각 하나씩 설명하도록 하겠다
Error<-NULL #우선 오차값을 저장할 Error라는 공간을 만든다.
for(i in 1:length(fit.tran$lambda)){ #tran에 fit한 모델들의 람다의 갯수만큼 반복한다.
tran.fit.pred<-predict(fit.tran, tran.model.names=tran.model.names[i], xnew=x[test,])
tran으로 fit한 모델의 예측값을 구한다. 이 때,
predict
fit.tran <- fitting한 모델
tran.model.names=tran.model.name[i]
fitting한 모델이 여러가지 인데,(우리 모델로 치면 람다가 100개라서 모델이 100개), 이러한 모델들 중 어느 것을 fit 할 것인가?
newx<-예측을 하려고 하는 test set은 무엇인가? 즉, test set의 x관측값들을 넣으면 된다.
which.min(Error)->loc.lambda
which.min(Error)## [1] 62fit.tran$lambda[which.min(Error)]->opt.lambda
fit.tran$lambda[which.min(Error)]## [1] 403.7017이 때, tran의 fit한 모델에서 가장 작은 오차를 같는 lambda를 찾는다.
full.model<-glmnet(x,y,alpha=0,lambda=grid)
full.model$beta[,loc.lambda]## AtBat Hits HmRun Runs RBI Walks
## 0.09238344 0.79744297 0.80505438 1.03411394 0.87728221 1.53862798
## Years CAtBat CHits CHmRun CRuns CRBI
## 1.83395501 0.01129748 0.05406625 0.38509943 0.10815217 0.11427651
## CWalks LeagueN DivisionW PutOuts Assists Errors
## 0.06137171 19.67455371 -72.32155123 0.15293830 0.02448004 -1.15688766
## NewLeagueN
## 9.45368575predict(full.model,type="coefficients", s=opt.lambda)## 20 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) 22.07163454
## AtBat 0.09238344
## Hits 0.79744297
## HmRun 0.80505438
## Runs 1.03411394
## RBI 0.87728221
## Walks 1.53862798
## Years 1.83395501
## CAtBat 0.01129748
## CHits 0.05406625
## CHmRun 0.38509943
## CRuns 0.10815217
## CRBI 0.11427651
## CWalks 0.06137171
## LeagueN 19.67455371
## DivisionW -72.32155123
## PutOuts 0.15293830
## Assists 0.02448004
## Errors -1.15688766
## NewLeagueN 9.45368575이 때, 다시 모든 관측치들로 최상의 lambda를 넣어서 fitting을 시켜주면 된다.
K-fold
1.glmnet는 cv.glmnet()으로 cross validation 옵션을 제공한다 2.cv.glmnet(x,y,alpha=옵션,nfolds=folding의 갯수)
set.seed(1)
cv.model<-cv.glmnet(x,y,alpha=0,nfolds=10)
names(cv.model)## [1] "lambda" "cvm" "cvsd" "cvup" "cvlo"
## [6] "nzero" "call" "name" "glmnet.fit" "lambda.min"
## [11] "lambda.1se" "index"우선, cv.glmnet로 cv된 모델을 만들었다. 자기 혼자서 cross validation을 진행 한 것이다.
이 때 names로 옵션을 살펴보자.
cvlo : cross validation lower bound : 로어바운드
cvm : cross validation median : 중앙값
cvup : cross validation upper bound : 어퍼바운드
이 세개는 one-standard-error-rule에 쓰인다
range<-cbind(cv.model$cvlo,cv.model$cvm,cv.model$cvup)
colnames(range)<-c("lower", "middle", "upper")
plot(cv.model)
이 때, 각각의 로어바운드, 어퍼바운드, 중앙을 나타내주었고, fit한 모델을 플랏에 표시해준 결과,
one-standard-error-Rule에 따라서 수축시킨 모델은 람다 값이 큰 쪽에서 그어져 있는 것을 볼 수 있고,
반대로, 에러가 제일 작은 것은 람다가 제일 작은 곳에 선이 그어져 있는 것을 확인할 수 있다.
log(cv.model$lambda.min)## [1] 3.239784log(cv.model$lambda.1se)## [1] 7.519336위와 같은 수식을 써서 구할수도 있다. lambda.min의 경우 오차가 가장 작을 때를 나타내고,
lambda.1se의 경우에는 최대한으로 축소시킬 수 있는 모델을 나타낸다.
which(cv.model$lambda==cv.model$lambda.min)## [1] 100which(cv.model$lambda==cv.model$lambda.1se)## [1] 54이 때, 100번째 람다에서는 오차를 가장 최소화하는 람다가 뽑한 것을 볼 수 있고, 54번째 람다에서는 축소를 최대로 시킨 모델임을 알 수 있다.
b.min=predict(cv.model, type="coefficients", s=cv.model$lambda.min)
b.1se=predict(cv.model, type="coefficients", s=cv.model$lambda.1se)
b.min## 20 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) 8.112693e+01
## AtBat -6.815959e-01
## Hits 2.772312e+00
## HmRun -1.365680e+00
## Runs 1.014826e+00
## RBI 7.130225e-01
## Walks 3.378558e+00
## Years -9.066800e+00
## CAtBat -1.199478e-03
## CHits 1.361029e-01
## CHmRun 6.979958e-01
## CRuns 2.958896e-01
## CRBI 2.570711e-01
## CWalks -2.789666e-01
## LeagueN 5.321272e+01
## DivisionW -1.228345e+02
## PutOuts 2.638876e-01
## Assists 1.698796e-01
## Errors -3.685645e+00
## NewLeagueN -1.810510e+01b.1se## 20 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) 159.796625704
## AtBat 0.102483884
## Hits 0.446840518
## HmRun 1.289060569
## Runs 0.702915317
## RBI 0.686866068
## Walks 0.925962427
## Years 2.714623467
## CAtBat 0.008746278
## CHits 0.034359576
## CHmRun 0.253594870
## CRuns 0.068874010
## CRBI 0.071334608
## CWalks 0.066114944
## LeagueN 5.396487429
## DivisionW -29.096663728
## PutOuts 0.067805862
## Assists 0.009201998
## Errors -0.235989097
## NewLeagueN 4.457548059이를 predict함수로 최종적인 모델을 추출할 수 있다.
predict(적합시킨 모델, type=“coefficient”, s=람다값)
(b.min[-1,])^2 %>% sum() ->L2.b.min
(b.1se[-1,])^2 %>% sum() ->L2.b.1se
cbind(L2.b.min,L2.b.1se)## L2.b.min L2.b.1se
## [1,] 18367.29 906.8121L2 norm을 보았을 때, 확실히 1se가 모델이 더욱 축소된 모습을 보여주고 있다.
Lasso
single split
set.seed(12345)
train<-sample(1:nrow(x), nrow(x)/2)
test<- (-train)sigle split은 전체 데이터에 대해서, train set과 test set을 한번만 나누는 것을 말한다.
이 때, test set은 오직 검증만을 위해 사용되고, train set은 파라미터 튜닝을 위해 사용된다.
즉 test set안에 있는 x값은 전혀 사용되지 않고, 오직 y값만 사용된다
이 떄, train set과 test set은 각각 반반으로 지정하였다.
y.test<-y[test]이 때, test set은 y값만 사용할 것이므로, test set안에있는 y값은 미리 추출하여 따로 저장해주자.
1. train set에 대해서 모델 fit 해주기
grid<-10^seq(10,-2,length=100)
fit.tran<-glmnet(x[train,],y[train],alpha=1,lambda=grid)
dim(coef(fit.tran))## [1] 20 100모델을 피팅한 결과, 100개에 람다에 대해서 20개의 회귀계수가 추출된 것을 알 수 있다.
이제 우리는 best lambda를 구하고 회귀계수를 추출해야 한다.
2. 모델 명 설정해주기
tran.model.name<-0:(length(fit.tran$lambda)-1)
tran.model.name## [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
## [26] 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
## [51] 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74
## [76] 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99그 후, 모델의 이름을 설정해준다. 모델의 이름을 0부터 99까지 설정한 이유는, glmnet()의 고유한 특성 때문이다.
colnames(fit.tran$beta)## [1] "s0" "s1" "s2" "s3" "s4" "s5" "s6" "s7" "s8" "s9" "s10" "s11"
## [13] "s12" "s13" "s14" "s15" "s16" "s17" "s18" "s19" "s20" "s21" "s22" "s23"
## [25] "s24" "s25" "s26" "s27" "s28" "s29" "s30" "s31" "s32" "s33" "s34" "s35"
## [37] "s36" "s37" "s38" "s39" "s40" "s41" "s42" "s43" "s44" "s45" "s46" "s47"
## [49] "s48" "s49" "s50" "s51" "s52" "s53" "s54" "s55" "s56" "s57" "s58" "s59"
## [61] "s60" "s61" "s62" "s63" "s64" "s65" "s66" "s67" "s68" "s69" "s70" "s71"
## [73] "s72" "s73" "s74" "s75" "s76" "s77" "s78" "s79" "s80" "s81" "s82" "s83"
## [85] "s84" "s85" "s86" "s87" "s88" "s89" "s90" "s91" "s92" "s93" "s94" "s95"
## [97] "s96" "s97" "s98" "s99"glmnet으로 fit시킨 beta들의 값들은 s00의 값으로 0부터 저장되기 때문에, 이러한 과정을 거쳐야한다. 즉, 임의로 glmnet이 이름을 지은 것을 그대로 따라간다고 보면 된다.
Error<-NULL
for (i in 1:length(fit.tran$lambda)){
fit.tran.pred<-predict(fit.tran, tran.model.name=tran.model.name[i], newx=x[test,])
Error[i]<-mean((fit.tran.pred[,i] - y.test)^2)
}
Error## [1] 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5
## [9] 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5
## [17] 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5
## [25] 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5
## [33] 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5
## [41] 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5
## [49] 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5
## [57] 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5 199378.5 182432.9
## [65] 156370.4 140265.6 130781.3 124290.5 118509.1 114790.8 112574.2 111248.3
## [73] 110445.9 110745.7 111025.7 110958.8 109744.6 108692.9 108219.6 108104.2
## [81] 108505.5 108893.4 108933.1 109221.9 110511.2 111767.3 112847.0 113766.4
## [89] 114462.8 115014.6 115452.3 115449.6 116117.6 116116.1 116115.7 116115.9
## [97] 116709.2 116707.6 116707.0 116706.6각각 하나씩 설명하도록 하겠다
Error<-NULL #우선 오차값을 저장할 Error라는 공간을 만든다.
for(i in 1:length(fit.tran$lambda)){ #tran에 fit한 모델들의 람다의 갯수만큼 반복한다.
tran.fit.pred<-predict(fit.tran, tran.model.names=tran.model.names[i], xnew=x[test,])
tran으로 fit한 모델의 예측값을 구한다. 이 때,
predict
fit.tran <- fitting한 모델
tran.model.names=tran.model.name[i]
fitting한 모델이 여러가지 인데,(우리 모델로 치면 람다가 100개라서 모델이 100개), 이러한 모델들 중 어느 것을 fit 할 것인가?
newx<-예측을 하려고 하는 test set은 무엇인가? 즉, test set의 x관측값들을 넣으면 된다.
which.min(Error)->loc.lambda
which.min(Error)## [1] 80fit.tran$lambda[which.min(Error)]->opt.lambda
fit.tran$lambda[which.min(Error)]## [1] 2.656088이 때, tran의 fit한 모델에서 가장 작은 오차를 같는 lambda를 찾는다.
full.model<-glmnet(x,y,alpha=1,lambda=grid)
full.model$beta[,loc.lambda]## AtBat Hits HmRun Runs RBI Walks
## -1.5600984 5.6931685 0.0000000 0.0000000 0.0000000 4.7505395
## Years CAtBat CHits CHmRun CRuns CRBI
## -9.5180241 0.0000000 0.0000000 0.5191611 0.6604074 0.3915415
## CWalks LeagueN DivisionW PutOuts Assists Errors
## -0.5326868 32.1125493 -119.2583540 0.2726207 0.1748164 -2.0567432
## NewLeagueN
## 0.0000000predict(full.model,type="coefficients", s=opt.lambda)## 20 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) 124.0894873
## AtBat -1.5600984
## Hits 5.6931685
## HmRun .
## Runs .
## RBI .
## Walks 4.7505395
## Years -9.5180241
## CAtBat .
## CHits .
## CHmRun 0.5191611
## CRuns 0.6604074
## CRBI 0.3915415
## CWalks -0.5326868
## LeagueN 32.1125493
## DivisionW -119.2583540
## PutOuts 0.2726207
## Assists 0.1748164
## Errors -2.0567432
## NewLeagueN .이 때, 다시 모든 관측치들로 최상의 lambda를 넣어서 fitting을 시켜주면 된다.
K-fold
1.glmnet는 cv.glmnet()으로 cross validation 옵션을 제공한다 2.cv.glmnet(x,y,alpha=옵션,nfolds=folding의 갯수)
set.seed(2)
cv.model<-cv.glmnet(x,y,alpha=1,nfolds=10)
names(cv.model)## [1] "lambda" "cvm" "cvsd" "cvup" "cvlo"
## [6] "nzero" "call" "name" "glmnet.fit" "lambda.min"
## [11] "lambda.1se" "index"우선, cv.glmnet로 cv된 모델을 만들었다. 자기 혼자서 cross validation을 진행 한 것이다.
이 때 names로 옵션을 살펴보자.
cvlo : cross validation lower bound : 로어바운드
cvm : cross validation median : 중앙값
cvup : cross validation upper bound : 어퍼바운드
이 세개는 one-standard-error-rule에 쓰인다
range<-cbind(cv.model$cvlo,cv.model$cvm,cv.model$cvup)
colnames(range)<-c("lower", "middle", "upper")
plot(cv.model)
이 때, 각각의 로어바운드, 어퍼바운드, 중앙을 나타내주었고, fit한 모델을 플랏에 표시해준 결과,
one-standard-error-Rule에 따라서 수축시킨 모델은 람다 값이 큰 쪽에서 그어져 있는 것을 볼 수 있고,
반대로, 에러가 제일 작은 것은 람다가 제일 작은 곳에 선이 그어져 있는 것을 확인할 수 있다.
log(cv.model$lambda.min)## [1] 0.8906821log(cv.model$lambda.1se)## [1] 4.239897위와 같은 수식을 써서 구할수도 있다. lambda.min의 경우 오차가 가장 작을 때를 나타내고,
lambda.1se의 경우에는 최대한으로 축소시킬 수 있는 모델을 나타낸다.
which(cv.model$lambda==cv.model$lambda.min)## [1] 51which(cv.model$lambda==cv.model$lambda.1se)## [1] 15이 때, 51번째 람다에서는 오차를 가장 최소화하는 람다가 뽑한 것을 볼 수 있고, 15번째 람다에서는 축소를 최대로 시킨 모델임을 알 수 있다.
b.min=predict(cv.model, type="coefficients", s=cv.model$lambda.min)
b.1se=predict(cv.model, type="coefficients", s=cv.model$lambda.1se)
cbind(b.min,b.1se)## 20 x 2 sparse Matrix of class "dgCMatrix"
## s1 s1
## (Intercept) 129.4155569 127.95694771
## AtBat -1.6130155 .
## Hits 5.8058915 1.42342566
## HmRun . .
## Runs . .
## RBI . .
## Walks 4.8469340 1.58214110
## Years -9.9724045 .
## CAtBat . .
## CHits . .
## CHmRun 0.5374550 .
## CRuns 0.6811938 0.16027975
## CRBI 0.3903563 0.33667715
## CWalks -0.5560143 .
## LeagueN 32.4646094 .
## DivisionW -119.3480842 -8.06171247
## PutOuts 0.2741895 0.08393604
## Assists 0.1855978 .
## Errors -2.1650837 .
## NewLeagueN . .이를 predict함수로 최종적인 모델을 추출할 수 있다.
predict(적합시킨 모델, type=“coefficient”, s=람다값)
(b.min[-1,])^2 %>% sum() ->L2.b.min
(b.1se[-1,])^2 %>% sum() ->L2.b.1se
cbind(L2.b.min,L2.b.1se)## L2.b.min L2.b.1se
## [1,] 15463.18 69.66661L2 norm을 보았을 때, 확실히 1se가 모델이 더욱 축소된 모습을 보여주고 있다.