본문 바로가기

R-PAGO 노트

능형(Ridge) 회귀, Lasso 이론 및 예제

검정MSE가 최소가 되는 모델을 찾아라!

그렇다. 우리의 목적은 언제나 검정오차(또는 검정MSE)가 최소인 모델을 찾는 것이다. 능형회귀와 Lasso는 이 목적때문에 탄생했다. 일종의 최적값(Optimization) 찾기이다. 편향(bias)과 분산(variance)사이의 최적점을 찾아 검정MSE가 최소가 되는 모델을 찾아주기 때문이다. 그러고 보면 지금까지의 다중회귀모델의 최적 모델은 검정MSE 관점에서는 좀 더 개선의 여지가 있었다는 말이 된다.


어떻게 만들까?

능형회귀 식

Ridge 식 

Lasso 식Lasso 식

능형회귀와 Lasso는 다중회귀모델이 사용하는 최소제곱(RSS)에 추가 항을 더해 식을 완성한다. 각각의 계수 추정치 는 상기식 각각을 최소로 하는 값이 된다. 우변의 두번째 항은 수축패널티(shrinkage penalty) 를 포함한다. 수축패널티는 >=0 이다.


왜 작동하는가?

ISLR Figure 6.5 결론부터 말하면 수축패널티를 0에서 무한까지 조절하면 (그래서 조율 파라미터(tuning parameter)라고 한다) 편향과 분산값 사이의 절충(trade-off)이 발생한다. 결국 이 둘의 합인 검정MSE가 최소인 지점을 찾게 되고 이 때의 계수 추정치를 구하게 된다. 상기 왼쪽 그래프는 능형회귀 그래프로 제곱편차(검정색), 분산(녹색), 검정MSE(보라색)와 가능한 최소 MSE(파선)사이의 관계를 나타낸다.(Lasso도 이와 유사하다.) 제곱편향과 분산이 절충되는 X 지점이 검정MSE가 최소인 곳이다.(가 약 30일때)


조금만 더 깊이: 편향-분산 절충이란?

앞에서 편향과 분산이란 용어가 등장한다. 나는 이전글 회귀모델의 변수 선택 방법에서 변수 선택과 편향, 분산과의 관계를 살짝 언급한 바 있다. 능형회귀와 Lasso가 왜 작동하는지를 이해하기 위해 좀더 자세히 설명하고자 한다.

이전 여러 모델들을 설명할 때 참값(Y)에 대한 추정값(Y_hat)의 정확성은 검정MSE를 통해 확인했다. 이 검정MSE 오차는 축소가능 오차와 축소불가능 오차로 나뉜다. 이 둘은 다음과 같이 정의된다.

축소가능 오차(reducible error):
아무리 좋은 모델이라고 해도 f에 대한 추정 f_hat은 f를 완벽하게 추정하지 못한다. 그러나 어느 정도까지는(좋은 모델) 오차를 최소화하여 f와 f_hat의 차이가 크지 않게 할 수 있다. 축소가능 오차는 이 두 사이를 최대한 줄일 수 있는 값(오차)를 의미한다.

축소불가능 오차(irreducible error):
설령 f를 완벽히 추정하여 Y에 대한 예측 결과인 Y_hat이 f의 형태가 된다 하더라도 예측 값은 여전히 오차를 가지게 된다. 이것은 Y는 e의 함수인데 설명변수 X로는 e를 예측할 수 없기 때문이다. 아무리 좋은 모델도 e와 관련된 변동성은 줄일 수 없는데 이것을 축소 불가능 오차라 한다.

상기 설명을 식으로 나타내면 다음와 같다. 

ISLR 식 2.2 ISLR 식 2.3

상기 식처럼 검정 MSE는 편향의 제곱(축소가능 오차)오차의 분산(축소불가능한 오차)로 구성된다. 다시 강조하지만 능형회귀와 Lasso는 검정MSE가 최소이기 위한 두 오차사이의 최적점을 찾는 것이다.


결과

능형회귀 결과

ISLR Figure 6.4: 능형회귀 

Lasso 결과 그래프ISLR Figure 6.6: Lasso

상기 능형회귀와 Lasso의 식으로 돌아가보자. 수축패널티가 0 일 경우, 두 식 모두 최소제곱 추정치와 같다. 반면 무한대일 경우, 두 회귀식의 계수 추정치는 0에 근접하게 될 것이다.(Lasso의 경우 일부 계수는 0이다.) 계수 추정치는 일반적으로 수축패널티가 증가하면 그것의 절대값은 감소한다. 그러나 어느 부근에서는 간혹 증가할 수 있다.(상기 두 그림의 왼쪽 그림 참조, y값은 각 표준화된 설명변수의 계수이다.)

이 모델을 마루(Ridge)와 올가미(Lasso)라 부르는 것은 이 두 모델의 모양 때문이라 추측된다.(상기 두 그림의 왼쪽 참조)

조율파라미터 대신 norm값으로 만든을 기준으로 표현하기도 한다. 분모의 beta는
최소제곱 계수 추정치들의 벡터로 로 정의된다. 능형회귀일 경우 분모는 0에서 최소제곱계수 추정치 벡터의 거리를 나타낸다.(l2 norm) 반면, Lasso는 최소제곱 계수 추정치의 절대값 합이다.(l1 norm) 따라서 상기 오른쪽 그림의 x축 값은 최대값이 1, 최소값은 0이 된다. 이때 수축패널티는 반대로 0과 무한대가 된다.


유의점

상기 결과를 만들 때 다음 사항을 유의하자. 

1. 수축페널티는 절편(b0)에는 적용되지 않는다. 절편을 수축할 이유가 없기 때문이다. 

2. 능형회귀와 Lasso는 다중회귀처럼 변수의 unit 변환시, 예를 들어 100배일때 단순히 계수는 1/100이 되지 않는다. 오히려 계수는 현저하게 바뀔 수 있다.(식의 우변의 두번째 항 때문에) 만약 scale에 의존적이지 않은 계수를 구하고 싶다면 표준편차가 1인 설명 변수의 표준화를 통해 계수 추정치를 나타내곤 한다.


능형회귀와 Lasso의 비교

  • 상기 식과 같이 우변의 두번째 항이 서로 다르다. 능형회귀는 계수의 제곱, Lasso는 절대값을 사용한다.
  • 능형회귀는 모든 설명 변수 p를 모두 포함한다. 이것은 수축페널티가 크다고 해도 각각의 설명변수의 계수가 0에 가까워지지 0은 아니기 때문이다.(수축페널티가 무한대일 경우는 제외)(상기 Figure6.4 참조)
  • 반면 Lasso는 일부 변수만 포함하고 나머지 변수의 계수는 0이 된다.(상기 Lasso 그래프 참조) 이는 수축패널티가 클수록 일부 변수이 계수는 0이 되기 때문이다.


능형회귀와 Lasso의 장단점

1. 능형회귀는 변수가 많고 계수의 크기가 거의 동일한 크기일 때 성능이 좋다.

 - 능형회귀는 높은 분산을 가지는 상황에서 가장 잘 작동하기 때문이다.

 - 변수가 많을 수록 f에 근접한 모델이 되어 편향은 줄어드나 오차의 폭은 커져 분산은 커지기 때문이다. 능형회귀는 약간의 편향 증가대비 분산은 크게 감소할 것이다.

2. Lasso는 적은 수의 설명변수가 상당히 큰 계수를 가질때 잘 작동한다. 

- 일부의 설명변수만 포함하므로 단순하고 해석력 높은 모델을 만든다. - 반응변수와 관련있는 설명변수는 신호이고 변수는 잡음이 된다.

ISLR Figure 6.8 ISLR Figure 6.9

위 두 그래프의 왼쪽은 Lasso에 대한 제곱편향(검은색), 분산(녹색), 검정MSE(보라색)이며 오른쪽 그래프는 각각의 값에 대한 Lasso(실선)와 능형회귀(파선)를 비교한 것이다. 첫번째 그림은 45개 변수일 때이며 아래 두번째는 설명변수가 2개일 때이다.변수가 많을 때는 능형회귀가 약간 좋은 성능인 반면 변수가 적을 때는 Lasso가 월등히 뛰어남을 알 수 있다.


교차검증: 수축페널티 및 변수선택

순서는 다음과 같다.

1. 훈련셋을 사용하여 조율파라미터에 따른 교차검증 오차를 계산한다. 

2. 오차가 가장 작을 때의 조율파라미터를 선택한다. 

3. 검정셋을 통해 검정오차를 확인한다. 

4. 모든 데이터와 선택된 조율파라미터를 사용해 적합된 최종 모델(변수의 계수)을 얻는다.


능형회귀 실전 예제

Lasso와 Ridge는 glmnet 패키지의 glmnet()을 사용한다. glmnet사용시 질적변수는 인식 못하므로 가변수로 바꿔준다. 또한 기존 회귀함수의 (y~x)형태가 아닌 x는 행렬, y는 벡터 인자로 바꿔야한다. 이때 절편 회귀계수는 수축할 필요가 없으므로 x행렬은 절편의 행렬 값 1을 제외한다. 다중회귀와의 비교를 위해 이전글 다중선형회귀 모델의 변수 선택을 위한 최종 Guideline과 같이 19개 설명 변수를 가지는 야구 타자 선수 데이터를 살펴보자.

library(ISLR)
names(Hitters)
##  [1] "AtBat"     "Hits"      "HmRun"     "Runs"      "RBI"      
##  [6] "Walks"     "Years"     "CAtBat"    "CHits"     "CHmRun"   
## [11] "CRuns"     "CRBI"      "CWalks"    "League"    "Division" 
## [16] "PutOuts"   "Assists"   "Errors"    "Salary"    "NewLeague"
dim(Hitters)
## [1] 322  20
Hitters<-na.omit(Hitters)
dim(Hitters)
## [1] 263  20
x<-model.matrix(Salary~.,Hitters)[,-1] #-1은 절편 1제거
y<-Hitters$Salary

다중선형회귀와 능형회귀훈련셋과 검정셋을 glmnet함수에서 능형회귀는 alpha=0, Lasso는 alpha=1이다. 조율파라미터는 0.01에서 큰수까지 설정해 최소제곱적합에서 영(수축)모델 모두가 검토될 수 있게 한다. default가 변수를 표준화하므로 비표준화를 할 경우 standardize=FALSE로 하면 된다. ridge 가 20 x 100인 것은 100개의 조율파라미터별 능형회귀의 회귀계수를 나타낸다.

library(glmnet)
sh<-10^seq(10,-2,length=100)
ridge<-glmnet(x,y,alpha=0, lambda=sh)
dim(coef(ridge))
## [1]  20 100
plot(ridge)

먼저 다중선형회귀(최소제곱)으로 할 경우의 검정MSE를 구해보자. 수축파라미터별 계수를 구할 때 어느정도 값의 변화가 없을 때 계수가 선택 되는데 thresh는 이 수렴 선택의 기준값을 의미한다.default는 1e-7이나 여기선 더 작은 1e-12로 하자. 구하고 싶은 조율파라미터 값을 입력하면 해당 값의 회귀계수를 얻을 수 있다. 이것은 predict함수로 가능하다. 함수내 s는 조율파라미터이다. s=0이면 최소제곱임을 상기하자. 인자 exact=T를 추가하면 더 정확한 값을 얻게 된다. 다중선형회귀, exact=T 포함 여부에 따라 계수가 얼마나 차이가 있는지 확인해 보면 거의 차이가 없으나 exact=T가 더 근접함을 알 수 있다. 다중선형회귀를 할 경우엔 glmnet보다는 lm을 사용하자

set.seed(1)
train<-sample(1:nrow(x),nrow(x)/2) #비복원 추출
test<-(-train)
y.test<-y[test]
ridge.tr<-glmnet(x[train,],y[train],alpha=0,lambda=sh,thresh=1e-12)
lm(y~x,subset=train)
## 
## Call:
## lm(formula = y ~ x, subset = train)
## 
## Coefficients:
## (Intercept)       xAtBat        xHits       xHmRun        xRuns  
##   299.42849     -2.54027      8.36682     11.64512     -9.09923  
##        xRBI       xWalks       xYears      xCAtBat       xCHits  
##     2.44105      9.23440    -22.93673     -0.18154     -0.11598  
##     xCHmRun       xCRuns        xCRBI      xCWalks     xLeagueN  
##    -1.33888      3.32838      0.07536     -1.07841     59.76065  
##  xDivisionW     xPutOuts     xAssists      xErrors  xNewLeagueN  
##   -98.86233      0.34087      0.34165     -0.64207     -0.67442
predict(ridge.tr,s=0,newx=x[test,],exact=T,type="coefficients")[1:20,]
##  (Intercept)        AtBat         Hits        HmRun         Runs 
## 299.42883596  -2.54014665   8.36611719  11.64400720  -9.09877719 
##          RBI        Walks        Years       CAtBat        CHits 
##   2.44152119   9.23403909 -22.93584442  -0.18160843  -0.11561496 
##       CHmRun        CRuns         CRBI       CWalks      LeagueN 
##  -1.33836534   3.32817777   0.07511771  -1.07828647  59.76529059 
##    DivisionW      PutOuts      Assists       Errors   NewLeagueN 
## -98.85996590   0.34086400   0.34165605  -0.64205839  -0.67606314
predict(ridge.tr,s=0,newx=x[test,],type="coefficients")[1:20,]
##  (Intercept)        AtBat         Hits        HmRun         Runs 
## 299.44467220  -2.53538355   8.33585019  11.59830815  -9.05971371 
##          RBI        Walks        Years       CAtBat        CHits 
##   2.45326546   9.21776006 -22.98239583  -0.18191651  -0.10565688 
##       CHmRun        CRuns         CRBI       CWalks      LeagueN 
##  -1.31721358   3.31152519   0.06590689  -1.07244477  59.75587273 
##    DivisionW      PutOuts      Assists       Errors   NewLeagueN 
## -98.94393005   0.34083276   0.34155445  -0.65312471  -0.65882930

다중선형회귀의 검정MSE는 114,783이다.

multi.linear<-predict(ridge.tr,s=0,newx=x[test,],exact=T)
mean((multi.linear-y.test)^2)
## [1] 114783.1

교차검증을 통해 최소검정 MSE가 예상하는 조율파라미터를 구해 보자. 알다시피 교차검증 자체로는 검정MSE를 알 수 없다. 교차검증은 검정MSE를 최소하게 만들 것으로 추정되는

어떤 변수를 선택하는 하는 것이 목적이다. 이 과정은 이전 다중회귀의 변수선택다중선형회귀 모델의 변수 선택을 위한 최종 Guideline이나 회귀트리회귀트리(regression tree)와 분류트리(classification tree) 예측와도 같다. 각각은 교차검증을 통해 변수 갯수와 트리 갯수를 선택했다.

이전 사용한 cv.glm와 비슷한 cv.glmnet()을 사용한다. 이전 nfolds의 default는 10-fold이나 nfolds를 사용하여 인자 변경이 가능하다. 263개는 1~10사이 약수가 없으므로 10개로 나누는 10folds로 하자. 조율파라미터가 211.7일 때 교차검증 오차가 가장 작다.

set.seed(1)
cv<-cv.glmnet(x[train,],y[train],alpha=0)
plot(cv)

bestlam<-cv$lambda.min
bestlam
## [1] 211.7416

검정MSE는 96,015로 다중회귀값 114,783보다 향상이 많이 되었음을 확인할 수 있다.

best.ridge.tr<-predict(ridge.tr,s=bestlam,newx=x[test,])
mean((best.ridge.tr-y.test)^2)
## [1] 96015.51

참고로 l2 norm값을 구하는 법도 알아보자. 능형회귀 식 우변의 두번째 항은 절편은 포함안되는 것을 상기하자. 조율파라미터 211가 가까운 값 64번째를 예시로 이용해 보자. l2 norm은 약 8이다.

ridge.tr$lambda[64]
## [1] 231.013
sqrt(sum(coef(ridge.tr)[-1,50]^2))
## [1] 8.012802

이제 위에서 구한 조율파라미터를 가지고 전체 데이터를 통해 최종 모델을 구하자. 능형회귀의 계수는 추가되는 각 회귀계수들의 제곱합 때문에 하나의 설명변수의 데이터 값(x)이 나눠지거나 곱해지질 경우 계수가 이에 비례해서 변하지 는다. 또한 심지어 다른 설명변수의 데이터 값이 변해도 영향을 받게 된다. 따라서 능형회귀는 항상 표준화 모델을 사용해야한다. 표준화 모델의 회귀 계수 해석은 예측 변수 1단위 증가가 아닌 예측변수 표준편차 1 증가했을 경우의 종속변수의 표준편차 변화로 해석한다. 이점은 다중선형회귀분성과 비교해 단점이다. 그러나 표준화된 능형회귀의 회귀 계수는 크기 자체가 종속변수의 영향을 나타내는 크기라는 점에서 장점도 가진다. 

아래 coding은 비표준화모델(standardize=FALSE)도 함께 구현해 보왔다.

best.ridge<-glmnet(x,y,alpha=0)
best.ridge.nonst<-glmnet(x,y,alpha=0,standardize=FALSE)
predict(best.ridge,s=bestlam,type="coefficients")
## 20 x 1 sparse Matrix of class "dgCMatrix"
##                        1
## (Intercept)   9.88487157
## AtBat         0.03143991
## Hits          1.00882875
## HmRun         0.13927624
## Runs          1.11320781
## RBI           0.87318990
## Walks         1.80410229
## Years         0.13074381
## CAtBat        0.01113978
## CHits         0.06489843
## CHmRun        0.45158546
## CRuns         0.12900049
## CRBI          0.13737712
## CWalks        0.02908572
## LeagueN      27.18227535
## DivisionW   -91.63411299
## PutOuts       0.19149252
## Assists       0.04254536
## Errors       -1.81244470
## NewLeagueN    7.21208390
predict(best.ridge.nonst,s=bestlam,type="coefficients")
## 20 x 1 sparse Matrix of class "dgCMatrix"
##                       1
## (Intercept) 82.26513333
## AtBat       -0.87493383
## Hits         2.90353827
## HmRun       -0.08722841
## Runs         0.93003807
## RBI          0.48449901
## Walks        2.78420201
## Years        0.01730279
## CAtBat      -0.32903975
## CHits        0.79338815
## CHmRun       0.22685420
## CRuns        0.99091253
## CRBI         0.62159540
## CWalks      -0.33052995
## LeagueN      0.08906050
## DivisionW   -0.24222149
## PutOuts      0.29053986
## Assists      0.28533890
## Errors      -0.64561559
## NewLeagueN   0.06345912


개선점 및 의문점

1. 30개의 랜덤 훈련셋과 검정셋 Case를 통해 검정MSE가 1표준편차내의 작은 모델을 선택하자.
2. 교차검증의 대전제는 선형모델로 예측할 경우 실제 f값도 선형이어야 한다는 점이다. 만약 실제 f가 선형이 아니라면 교차검증 자체가 무의미하다. 3. 능형회귀 랜덤포리스트, 회귀트리,부스팅과 같이 선택된 변수의 p값 기준 유의미성을 살펴볼 필요가 있다. 만약 계수가 0에 가깞다면 이 자체가 계수 제거와 유사하다고 할 수 있을까?(p-value>0.05)


Lasso 실전 예제

능형회귀를 하니 Lasso는 덤으로 할 수 있겠다. alpha=1이고 나머지는 같다.

lasso.tr<-glmnet(x[train,],y[train],alpha=1,lambda=sh)
plot(lasso.tr)

교차검증 통해 조율파라미터 값과 검정MSE를 구하자. 조율파라미터가 16.8일때 교차검증 오차가 최소임을 알 수 있다.

set.seed(1)
cv.lasso<-cv.glmnet(x[train,],y[train],alpha=1)
plot(cv.lasso)

bestlam.lasso<-cv.lasso$lambda.min
bestlam.lasso
## [1] 16.78016

검정 MSE는 100,743이다. 비록 능형회귀보다 높은 값이지만 계수를 확인해 보면 7개 변수외에 다른 변수는 0이다. lambda별로 구할 필요가 없을 때는 glmnet함수내 lambda를 표시 하지 않아도 된다.

lasso.tr.pred<-predict(lasso.tr,s=bestlam.lasso,newx=x[test,])
mean((lasso.tr.pred-y.test)^2)
## [1] 100743.4
best.lasso<-glmnet(x,y,alpha=1,lambda=sh)
best.lasso.nonst<-glmnet(x,y,alpha=0,standardize=FALSE)
predict(best.lasso,s=bestlam.lasso,type="coefficients")
## 20 x 1 sparse Matrix of class "dgCMatrix"
##                        1
## (Intercept)   18.5394844
## AtBat          .        
## Hits           1.8735390
## HmRun          .        
## Runs           .        
## RBI            .        
## Walks          2.2178444
## Years          .        
## CAtBat         .        
## CHits          .        
## CHmRun         .        
## CRuns          0.2071252
## CRBI           0.4130132
## CWalks         .        
## LeagueN        3.2666677
## DivisionW   -103.4845458
## PutOuts        0.2204284
## Assists        .        
## Errors         .        
## NewLeagueN     .
predict(best.ridge.nonst,s=bestlam.lasso,type="coefficients")
## 20 x 1 sparse Matrix of class "dgCMatrix"
##                       1
## (Intercept) 82.26513333
## AtBat       -0.87493383
## Hits         2.90353827
## HmRun       -0.08722841
## Runs         0.93003807
## RBI          0.48449901
## Walks        2.78420201
## Years        0.01730279
## CAtBat      -0.32903975
## CHits        0.79338815
## CHmRun       0.22685420
## CRuns        0.99091253
## CRBI         0.62159540
## CWalks      -0.33052995
## LeagueN      0.08906050
## DivisionW   -0.24222149
## PutOuts      0.29053986
## Assists      0.28533890
## Errors      -0.64561559
## NewLeagueN   0.06345912


의문점

glmnet함수내에서 lambda값을 지정했을 경우와 안지정했을 경우 회귀계수의 결과값은 같으나 절편계수는 다르다. 왜 그럴까?