본문 바로가기
머신러닝, 딥러닝

[머신러닝] 5. 선형모델 선택 및 Regularization- Ridge, Lasso regression, PCR, PLS (feat. R Code)

ღyeonyღ 2021. 4. 25. 01:06

 

선형 회귀 모델 구축에 사용되는 최소제곱법이 아닌 다른 적합절차를 사용하는 이유는?

  • 예측 정확도: n이 p보다 아주 크지 않으면 최소제곱적합에 많은 변동이 존재할 수 있어 과적합을 초래할 수 있다. p>n가 되면 분산이 무한대가 되어 최소제곱 방법은 전혀 사용할 수 없게 된다.
  • 모델 해석력: 최소제곱방법으로 반응변수와 관련이 없는 변수들도 정확하게 0인 계수 추정치를 얻게 될 가능성은 거의 없다. 관련이 없는 변수들은 모델을 필요없이 복잡하게 만들기 때문에, 관련없는 변수들을 제외하는 기법들이 필요하다.

최소제곱법 대신 사용할 대안은?

  • 서브셋(부분집합) 선택
  • 수축(shrinkage)
  • 차원축소(Dimension Reduction)

부분집합 선택- 단계적 선택

1. 전진 단계적 선택

  • 전진 단계적 선택은 설명변수가 하나도 포함되지 않은 모델에서 시작하여 모든 설명변수가 모델에 포함될 때까지 한번에 하나씩 설명변수를 추가한다.
  • 최상의 부분집합 선택에 p개 설명변수들의 서브셋을 포함하는 2^p개의 가능한 모든 모델들을 고려하는 반면에 전진 단계적 선택은 훨씬 적은 수의 모델을 고려한다는 계산적 장점을 가지고 있다.

2. 후진 단계적 선택

  • p개의 설명변수 모두를 포함하는 완전 최소제곱 모델을 가지고 시작하고, 한번에 하나씩 반복적으로 유용성이 가장 적은 설명변수를 제외한다.
  • 후진 단계적 선택법은 표본의 수 n이 설명변수의 수 p보다 커야 사용할 수 있다.

3. 하이브리드 방식

  • 변수들이 모델에 순차적으로 추가된다는 점에서 전진 단계적 선택과 비슷하다. 하지만 새로운 변수를 추가한 후에 모델 적합을 더 이상 향상시키지 않는 변수가 있으면 제거할 수도 있다.

Shrinkage 방법

  • 계수 추정치를 제한(constrains)하거나 규칙화(regularizes)하는 기법을 사용하여 p개의 설명변수 모두를 포함하는 모델 적합할 수 있다.
  • 계수 추정치들을 수축하는 것은 추정치들의 분산을 상당히 줄일 수 있는 것으로 밝혀져 있다.
  • 회귀계수들을 영으로 수축하기 위한 두 가지 가장 잘 알려진 기법은 능형회귀(ridge regression)과 라쏘회귀(lasso regression)이다.

1. 능형회귀

능형회귀 계수 추정치는 다음 식을 최소로 하는 값이다.

  • 약간 다른 수량을 최소화하여 계수들을 추정한다는 점을 제외하면 최소제곱과 아주 유사하다.
  • lambda=0일 때 페널티 항의 영향이 없어 능형회귀는 최소제곱 추정치를 제공한다. 하지만 lamda가 무한대에 가까워짐에 따라 수축 페널티의 영향이 커져 능형회귀 계수 추정치는 0으로 근접한다.
  • 단 하나의 계수 추정치들의 집합을 생성하는 최소제곱과 달리, 능형회귀는 lambda의 값 각각에 대해 다른 집합의 계수 추정치를 생성한다.

 

능형회귀가 최소제곱보다 나은 이유

  • 최소제곱에 대한 능형회귀의 장점은 편향-분산 절충에 있다. lambda가 증가하면 능형회귀 적합의 유연성이 감소하게 되어 분산은 감소하지만 편향은 증가한다.
  • 변수의 수 p가 관측치의 수 n만큼 클 때 최소제곱 추정치들은 변동이 아주 클 것이다. 또한 p>n이면 최소제곱 추정치는 심지어 유일한 해도 가지지 않지만 능형회귀는 약간의 편향 증가로 분산을 크게 감소하도록 절충하여 잘 동작할 수 있다. 따라서 능형회귀는 최소제곱 추정치가 높은 분산을 가지는 상항에서 가장 잘 동작한다.
  • 능형회귀는 2^p개의 모델을 검색해야 하는 부분집합 선택에 비해 상당한 계산상의 장점이 있다.

하지만 한 가지 분명한 단점이 있다

  • 일반적으로 변수들의 서브셋만 포함하는 모델들을 선택하는 최상의 부분집합 선택과 단계적 선택과 달리, 능형회귀는 최종 모델에 p개 설명변수 모두를 포함할 것이다. 모든 계수들은 0을 향해 수축할 뿐 0이 되지 않을 것이기 때문이다.
  • 이것을 예측 정확도에 있어서는 문제가 되지 않을 수 있지만 변수의 수 p가 상당히 큰 설정에서 모델을 해석하는 데 어려움을 초래할 수 있다.

2. 라쏘회귀

  • 능형회귀에서와 같이 계수 추정치들을 0으로 수축하지만, lasso에서는 L1 페널티는 lambda가 충분히 클 경우 계수 추정치들의 일부를 정확히 0이 되게 하는 효과를 갖는다.

왼쪽이 라쏘회귀, 오른쪽이 능형회귀에 대한 오차의 등고선과 제한함수. 파란색 영역은 제한영역이고 붉은색 타원들은 RSS의 등고선이다.

차원축소 방법

1. 주성분회귀

  • PCA는 변동이 크게 되는 새로운 축을 주성분으로 하여 데이터의 차원을 줄이는 기법이다.
  • PCR기법은 설명변수 X를 가장 잘 나타내는 선형결합 또는 방향을 찾아내는 것이 관련된다. 이러한 방향들은 비지도 방식으로 식별된다.

2. 부분최소제곱(PLS; partial least squares)

  • PCR처럼 차원축소 방법이며, 원래 변수들의 선형결합인 새로운 변수들의 집합을 먼저 찾은 다음 이 M개 새로운 변수들을 이용한 최소제곱을 통해 선형모델을 적합한다.
  • 그러나 PCR과는 달리 PLS는 이러한 새로운 변수들을 지도식 방식으로 찾는다. 즉, PLS는 반응변수 Y를 이용하여 원래의 변수들을 잘 근사할 뿐만 아니라 반응변수와 관련이 있는 새로운 변수들을 식별한다.
  • PLS는 디지털화된 분광분석 신호들로부터 많은 변수들이 발생되는 계량분석화학 분야에서 널리 사용된다.
  • 실질적으로는 능형회귀나 PCR과 비슷한 성능을 보인다. PLS의 지도식 차원축소는 편향을 줄일 수 있지만 반면에 분산을 증가시킬 가능성이 있다.

고차원의 고려

1. 고차원에서 무엇이 문제인가?

  • 변수의 수 p가 관측치 수 n만큼 크거나 또는 n보다 더 클 때, 최소제곱법을 사용할 수 없다.
  • 설명변수들과 반응변수 사이에 정말로 상관관계가 있는지 여부에 관계없이 최소제곱은 데이터에 완벽하게 적합되어 잔차들이 0인 계수 추정치들의 집합을 제공할 것이다.

2. 고차원에서의 회귀

  • 변수의 수 p를 이용해 lasso가 수행되었고, p개의 변수 중 20개는 반응변수와 관련이 있다고 가정한다. boxplot은 조율 파라미터 lambda가 3가지 다른 값을 가질 때의 검정 MSE를 보여준다.
  • p=20일 때, 가장 낮은 검정 MSE는 조절 정도가 가장 작을 때 얻어진다. p=50일 때, 가장 낮은 검정 MSE는 조절 정도가 상당히 클 때 얻어진다. p=2000일 때, lasso는 조절 정도에 관계없이 결과가 좋지 않은데, 그 이유는 2,000개의 변수 중 20개만이 실제로 결과와 관련이 있기 때문이다.
  • 모델을 적합하는 데 사용된 변수의 수가 증가함에 따라 적합된 모델의 질도 높아지지만은 않는다는 것을 알 수 있다.

 

 

R Code

최상의 서브셋 선택

library(ISLR)
head(Hitters)
dim(Hitters)
##[1] 322  20
sum(is.na(Hitters$Salary))
##[1] 59

Hitters=na.omit(Hitters)
dim(Hitters)
##[1] 263  20
sum(is.na(Hitters$Salary))
##[1] 0
  • Salary 변수에 누락된 선수 59명의 데이터 모두 제거
library(leaps)
regfit.full=regsubsets(Salary~., Hitters)
summary(regfit.full)
##Subset selection object
##Call: regsubsets.formula(Salary ~ ., Hitters)
##19 Variables  (and intercept)
##           Forced in Forced out
##AtBat          FALSE      FALSE
##Hits           FALSE      FALSE
##HmRun          FALSE      FALSE
##Runs           FALSE      FALSE
##RBI            FALSE      FALSE
##Walks          FALSE      FALSE
##Years          FALSE      FALSE
##CAtBat         FALSE      FALSE
##CHits          FALSE      FALSE
##CHmRun         FALSE      FALSE
##CRuns          FALSE      FALSE
##CRBI           FALSE      FALSE
##CWalks         FALSE      FALSE
##LeagueN        FALSE      FALSE
##DivisionW      FALSE      FALSE
##PutOuts        FALSE      FALSE
##Assists        FALSE      FALSE
##Errors         FALSE      FALSE
##NewLeagueN     FALSE      FALSE
##1 subsets of each size up to 8
##Selection Algorithm: exhaustive
##         AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns
##1  ( 1 ) " "   " "  " "   " "  " " " "   " "   " "    " "   " "    " "  
##2  ( 1 ) " "   "*"  " "   " "  " " " "   " "   " "    " "   " "    " "  
##3  ( 1 ) " "   "*"  " "   " "  " " " "   " "   " "    " "   " "    " "  
##4  ( 1 ) " "   "*"  " "   " "  " " " "   " "   " "    " "   " "    " "  
##5  ( 1 ) "*"   "*"  " "   " "  " " " "   " "   " "    " "   " "    " "  
##6  ( 1 ) "*"   "*"  " "   " "  " " "*"   " "   " "    " "   " "    " "  
##7  ( 1 ) " "   "*"  " "   " "  " " "*"   " "   "*"    "*"   "*"    " "  
##8  ( 1 ) "*"   "*"  " "   " "  " " "*"   " "   " "    " "   "*"    "*"  
##         CRBI CWalks LeagueN DivisionW PutOuts Assists Errors NewLeagueN
##1  ( 1 ) "*"  " "    " "     " "       " "     " "     " "    " "       
##2  ( 1 ) "*"  " "    " "     " "       " "     " "     " "    " "       
##3  ( 1 ) "*"  " "    " "     " "       "*"     " "     " "    " "       
##4  ( 1 ) "*"  " "    " "     "*"       "*"     " "     " "    " "       
##5  ( 1 ) "*"  " "    " "     "*"       "*"     " "     " "    " "       
##6  ( 1 ) "*"  " "    " "     "*"       "*"     " "     " "    " "       
##7  ( 1 ) " "  " "    " "     "*"       "*"     " "     " "    " "       
##8  ( 1 ) " "  "*"    " "     "*"       "*"     " "     " "    " "
  • leaps라이브러리에 포함된 regsubsets함수는 주어진 수의 설명변수를 포함하는 최고의 모델을 식별함으로써 최상의 서브셋을 선택한다. 여기서 최상 또는 최고는 RSS를 사용하여 수량화한다.
regfit.full=regsubsets(Salary~., data=Hitters, nvmax=19)
reg.summary=summary(regfit.full)

par(mfrow=c(2, 2))
plot(reg.summary$rss, xlab='Number of Variables', ylab='RSS', type='l')
plot(reg.summary$adjr2, xlab='Number of Variables', ylab='Adjusted RSq', type='l')
points(which.max(reg.summary$adjr2), reg.summary$adjr2[11], col='red', cex=2, pch=20)
plot(reg.summary$cp, xlab='Number of Variables', ylab='Cp', type='l')
points(which.min(reg.summary$cp), reg.summary$cp[10], col='red', cex=2, pch=20)
plot(reg.summary$bic, xlab='Number of Variables', ylab='BIC', type='l')
points(which.min(reg.summary$bic), reg.summary$bic[6], col='red', cex=2, pch=20)

  • regsubsets에 nvmax옵션을 사용해 원하는 변수 수만큼 결과를 얻을 수 있다.
  • 모델에 대한 RSS, 조정된 R^2, Cp, BIC를 한꺼번에 그려 어느 모델을 선택할지 결정한다.
plot(regfit.full, scale='r2')
plot(regfit.full, scale='adjr2')
plot(regfit.full, scale='Cp')
plot(regfit.full, scale='bic')

  • regsubsets함수는 내장된 plot함수를 가지며, 이것은 주어진 수의 설명변수를 갖는 최상의 모델에 포함되는 변수들을 나타내는 데 사용될 수 있다.
coef(regfit.full, 6)
## (Intercept)        AtBat         Hits        Walks         CRBI 
##  91.5117981   -1.8685892    7.6043976    3.6976468    0.6430169 
##   DivisionW      PutOuts 
##-122.9515338    0.2643076
  • 가장 낮은 BIC를 가지는 모델은 AtBat, Hits, Walk, CRBI, DivisionW, PutOuts만 포함하는 6-변수 모델이다. coef함수를 사용하여 이 모델과 관련된 계수 추정치를 확인할 수 있다.

전진 및 후진 단계적 선택

regfit.fwd=regsubsets(Salary~., data=Hitters, nvmax=19, method='forward')
summary(regfit.fwd)
regfit.bwd=regsubsets(Salary~., data=Hitters, nvmax=19, method='backward')
summary(regfit.bwd)

coef(regfit.full, 7)
## (Intercept)        AtBat         Hits        Walks         CRBI 
## 109.7873062   -1.9588851    7.4498772    4.9131401    0.8537622 
##      CWalks    DivisionW      PutOuts 
##  -0.3053070 -127.1223928    0.2533404 
coef(regfit.fwd, 7)
## (Intercept)        AtBat         Hits        Walks        CRuns 
## 105.6487488   -1.9762838    6.7574914    6.0558691    1.1293095 
##      CWalks    DivisionW      PutOuts 
##  -0.7163346 -116.1692169    0.3028847 
coef(regfit.bwd, 7)
## (Intercept)        AtBat         Hits        Walks        CRuns 
## 105.6487488   -1.9762838    6.7574914    6.0558691    1.1293095 
##      CWalks    DivisionW      PutOuts 
##  -0.7163346 -116.1692169    0.3028847
  • regsubsets함수는 method='forward', method='backward'를 사용하여 전진 또는 후진 단계적 선택을 수행하는 데도 사용할 수 있다.
  • 전진 단계적 선택, 후진 단계적 선택, 최상의 서브셋 선택에 의해 식별된 최상의 7-변수 모델들은 서로 다른 것을 확인할 수 있다.

검증셋 기법과 교차검증을 사용한 모델 선택

set.seed(1)
train=sample(c(TRUE, FALSE), nrow(Hitters), rep=TRUE)
test=(!train)

regfit.best=regsubsets(Salary~., data=Hitters[train,], nvmax=19)
test.mat=model.matrix(Salary~., data=Hitters[test,])

val.errors=rep(NA, 19)
for (i in 1:19){
	coefi=coef(regfit.best, id=i)
	pred=test.mat[,names(coefi)]%*%coefi
	val.errors[i]=mean((Hitters$Salary[test]-pred)^2)}
val.errors
## [1] 164377.3 144405.5 152175.7 145198.4 137902.1 139175.7 126849.0 136191.4
## [9] 132889.6 135434.9 136963.3 140694.9 140690.9 141951.2 141508.2 142164.4
## [17] 141767.4 142339.6 142238.2

coef(regfit.best, which.min(val.errors))
##  (Intercept)        AtBat         Hits        Walks        CRuns 
##  67.1085369   -2.1462987    7.0149547    8.0716640    1.2425113 
##      CWalks    DivisionW      PutOuts 
##  -0.8337844 -118.4364998    0.2526925
  • model.matrix함수는 데이터를 모형행렬로 만들어주는 함수다.
  • for문을 실행하여 모델에 대한 계수들을 regfit.best에서 추출하여 coefi에 저장하고, 이 계수들을 검정모델 행렬에 곱하여 예측값을 구하고 검정 MSE를 계산한다.
predict.regsubsets=function(object, newdata, id, ...){
	form=as.formula(object$call[[2]])
	mat=model.matrix(form, newdata)
	coefi=coef(object, id=id)
	xvars=names(coefi)
	mat[,xvars]%*%coefi}
  • 이렇게 함수를 사용할 수도 있다. 실제로 적용할 때는 object에 regfit.best와 같은 모델 명을 입력하면 된다.
k=10
set.seed(1)
folds=sample(1:k, nrow(Hitters), replace=TRUE)
folds
##  [1]  9  4  7  1  2  7  2  3  1  5  5 10  6 10  7  9  5  5  9  9  5  5  2 10
## [25]  9  1  4  3  6 10 10  6  4  4 10  9  7  6  9  8  9  7  8  6 10  7  3 10
## [49]  6  8  2  2  6  6  1  3  3  8  6  7  6  8  7  1  4  8  9  9  7  4  7  6
## [73]  1  5  6  1  9  7  7  3  6  2 10 10  7  3  2 10  1 10 10  8 10  5  7  8
## [97]  5  6  8  1  3 10  3  1  6  6  4  9  5  1  3  6  3  7  3  3  1  9  2  8
##[121]  6  1  2  7  7  4  9  8  3  5  3  4  2  1  7  9 10 10  2  2  3  1  2  3
##[145]  3  3  8  9  2 10  8 10  4  5  9  5  7  5  6  4  2  1  3  8  9  6  1  4
##[169]  5  9  5  8  4  1  9  5  1  5  4 10 10  9  8  5  5  6  6  2  2  8  4 10
##[193]  8  5  5  8  8  7  4  4  1 10  4  9  9  9  9  6  6  4  3  3  9  9  7  9
##[217]  5  7  4  4 10  8  1 10  2 10  1  1  4  5  5  6  9  8  5  1  2  1  8  5
##[241]  8 10  7  7  2  9  4  2  5  2  4  3  6  9  7  5  5  1  1 10  1  3 10

cv.errors=matrix(NA, k, 19, dimnames=list(NULL, paste(1:19)))
cv.errors
##       1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19
## [1,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [2,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [3,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [4,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [5,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [6,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [7,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [8,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [9,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
##[10,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
  • k-fold 교차검증을 사용하여 최상의 서브셋을 선택하기 위해, 10-fold 중 하나의 그룹에 할당하는 벡터와 그 결과를 저장할 행렬을 생성한다.
for (j in 1:k){
	best.fit=regsubsets(Salary~., data=Hitters[folds!=j,], nvmax=19)
	for (i in 1:19){
		pred=predict.regsubsets(best.fit, Hitters[folds==j, ], id=i)
		cv.errors[j, i]=mean((Hitters$Salary[folds==j]-pred)^2)
		}
	}
    
mean.cv.errors=apply(cv.errors, 2, mean)
mean.cv.errors
##       1        2        3        4        5        6        7        8 
##149821.1 130922.0 139127.0 131028.8 131050.2 119538.6 124286.1 113580.0 
##       9       10       11       12       13       14       15       16 
##115556.5 112216.7 113251.2 115755.9 117820.8 119481.2 120121.6 120074.3 
##      17       18       19 
##120084.8 120085.8 120403.5
  • for문을 통해 10*19행렬이 얻어지며, 원소(i, j)는 i번째 교차검증 fold와 최고의 j-변수 모델에 대한 검정MSE다.
  • regsubsets함수에는 predict가 적용되지 않으므로 위에서 새롭게 만들어준 predict.regsubsets함수로 예측값을 구했다.
  • apply함수를 이용해 이 행렬의 열별로 평균을 구하면 벡터가 얻어지는데, 이 벡터의 j번째 원소는 j-변수 모델에 대한 교차검증 오차다.
plot(mean.cv.errors, type='b')

  • 교차검증은 11-변수 모델을 선택한다.

릿지(Ridge)회귀

x=model.matrix(Salary~., Hitters)[,-1]
y=Hitters$Salary

library(glmnet)
grid=10^seq(10, -2, length=100)
ridge.mod=glmnet(x, y, alpha=0, lambda=grid)
  • glmnet라이브러리의 glmnet모델을 사용해 릿지회귀와 라쏘회귀를 적합할 수 있다.
  • alpha인자가 0이면 릿지회귀를 적합하고, 1이라면 라쏘모델을 적합한다.
  • glmnet함수는 기본적으로 변수들을 표준화하여 scale이 동일하게 한다. 기본 설정을 끄려면 standardize=FALSE 옵션을 사용하면 된다.
dim(coef(ridge.mod))
##[1]  20 100
  • coef의 행렬의 경우, 20행(각 설명변수와 절편에 대해 하나씩)과 100열(각 lambda값에 대해 하나씩)을 가진다.
ridge.mod$lambda[50]
##[1] 11497.57
coef(ridge.mod)[,50]
##  (Intercept)         AtBat          Hits         HmRun          Runs 
##407.356050200   0.036957182   0.138180344   0.524629976   0.230701523 
##          RBI         Walks         Years        CAtBat         CHits 
##  0.239841459   0.289618741   1.107702929   0.003131815   0.011653637 
##       CHmRun         CRuns          CRBI        CWalks       LeagueN 
##  0.087545670   0.023379882   0.024138320   0.025015421   0.085028114 
##    DivisionW       PutOuts       Assists        Errors    NewLeagueN 
## -6.215440973   0.016482577   0.002612988  -0.020502690   0.301433531 
sqrt(sum(coef(ridge.mod)[-1,50]^2))
##[1] 6.360612
  • lambda가 11,498일 때의 계수들과 이들의 L2 norm을 보여준다.
predict(ridge.mod, s=50, type='coefficients')[1:20,]
##  (Intercept)         AtBat          Hits         HmRun          Runs 
## 4.876610e+01 -3.580999e-01  1.969359e+00 -1.278248e+00  1.145892e+00 
##          RBI         Walks         Years        CAtBat         CHits 
## 8.038292e-01  2.716186e+00 -6.218319e+00  5.447837e-03  1.064895e-01 
##       CHmRun         CRuns          CRBI        CWalks       LeagueN 
## 6.244860e-01  2.214985e-01  2.186914e-01 -1.500245e-01  4.592589e+01 
##    DivisionW       PutOuts       Assists        Errors    NewLeagueN 
##-1.182011e+02  2.502322e-01  1.215665e-01 -3.278600e+00 -9.496680e+00
  • predict함수와 s=50옵션을 설정하여 다음과 같이 새로운 lambda값 50에 대한 릿지회귀계수를 얻을 수 있다.
set.seed(1)
train=sample(1:nrow(x), nrow(x)/2)
test=(-train)
y.test=y[test]

ridge.mod=glmnet(x[train,], y[train], alpha=0, lambda=grid, thresh=1e-12)
ridge.pred=predict(ridge.mod, s=4, newx=x[test,])
mean((ridge.pred-y.test)^2)
##[1] 142199.2
mean((mean(y[train])-y.test)^2) #절편만 가진 모델의 MSE
##[1] 224669.9
  • thresh옵션은 수렴의 기준값이 된다. default는 1e-7이지만, 여기에선 1e-12로 더 작은 값을 설정했다.
  • lambda가 4인 릿지회귀를 적합하고 검정셋으로 MSE를 계산한 결과다. lambda가 4인 릿지회귀모델을 적합하는 것이 절편만 가진 모델을 적합하는 것보다 훨씬 낮은 MSE를 초래한다.
set.seed(1)
cv.out=cv.glmnet(x[train,], y[train], alpha=0)
plot(cv.out)
bestlam=cv.out$lambda.min
bestlam
##[1] 326.0828

ridge.pred=predict(ridge.mod, s=bestlam, newx=x[test,])
mean((ridge.pred-y.test)^2)
##[1] 139856.6

out=glmnet(x, y, alpha=0)
predict(out, type='coefficients', s=bestlam)[1:20,]
## (Intercept)        AtBat         Hits        HmRun         Runs 
## 15.44383135   0.07715547   0.85911581   0.60103107   1.06369007 
##         RBI        Walks        Years       CAtBat        CHits 
##  0.87936105   1.62444616   1.35254780   0.01134999   0.05746654 
##      CHmRun        CRuns         CRBI       CWalks      LeagueN 
##  0.40680157   0.11456224   0.12116504   0.05299202  22.09143189 
##   DivisionW      PutOuts      Assists       Errors   NewLeagueN 
##-79.04032637   0.16619903   0.02941950  -1.36092945   9.12487767

  • cv.glmnet함수를 사용하여 조율 파라미터를 선택할 수 있다. 기본적으로 10-fold교차검증을 수행하는데, nfolds인자를 사용하여 변경이 가능하다.
  • cv.glmnet의 lambda.min을 통해 교차검증 오차가 가장 작은 lambda값을 구할 수 있다.
  • lambda가 4일 때 얻은 검정 MSE보다 더 나은 결과인 것을 확인할 수 있다.
  • 마지막으로 교차검증에 의해 선택된 lambda값을 사용하여 릿지회귀를 적합하고, predict에 type='coefficients'옵션을 추가하여 계수 추정치를 출력했다.

라쏘(LASSO)회귀

lasso.mod=glmnet(x[train,], y[train], alpha=1, lambda=grid)
plot(lasso.mod)

  • 계수 그래프를 통해 조율 파라미터의 선택에 따라 몇몇 계수는 그 값이 정확하게 0이 됨을 알 수 있다.
set.seed(1)
cv.out=cv.glmnet(x[train,], y[train], alpha=1)
plot(cv.out)

bestlam=cv.out$lambda.min
lasso.pred=predict(lasso.mod, s=bestlam, newx=x[test,])
mean((lasso.pred-y.test)^2)
##[1] 143673.6

  • 라쏘회귀 적합 결과, 교차검증으로 선택된 lambda를 사용한 능형회귀의 검정 MSE와 거의 같다.
out=glmnet(x, y, alpha=1, lambda=grid)
lasso.coef=predict(out, type='coefficients', s=bestlam)[1:20, ]
lasso.coef
##  (Intercept)         AtBat          Hits         HmRun          Runs 
##   1.27479059   -0.05497143    2.18034583    0.00000000    0.00000000 
##          RBI         Walks         Years        CAtBat         CHits 
##   0.00000000    2.29192406   -0.33806109    0.00000000    0.00000000 
##       CHmRun         CRuns          CRBI        CWalks       LeagueN 
##   0.02825013    0.21628385    0.41712537    0.00000000   20.28615023 
##    DivisionW       PutOuts       Assists        Errors    NewLeagueN 
##-116.16755870    0.23752385    0.00000000   -0.85629148    0.00000000 

lasso.coef[lasso.coef!=0]
##  (Intercept)         AtBat          Hits         Walks         Years 
##   1.27479059   -0.05497143    2.18034583    2.29192406   -0.33806109 
##       CHmRun         CRuns          CRBI       LeagueN     DivisionW 
##   0.02825013    0.21628385    0.41712537   20.28615023 -116.16755870 
##      PutOuts        Errors 
##   0.23752385   -0.85629148
  • 이 예에서 19개의 계수 추정치 중 12개가 0이다. 따라서, 교차검증으로 선택된 lambda값의 lasso모델은 단지 7개의 변수를 포함한다.

주성분회귀(PCR)

library(pls)
set.seed(2)
pcr.fit=pcr(Salary~., data=Hitters, scale=TRUE, validation='CV')
summary(pcr.fit)
##Data:   X dimension: 263 19 
##        Y dimension: 263 1
##Fit method: svdpc
##Number of components considered: 19
##
##VALIDATION: RMSEP
##Cross-validated using 10 random segments.
##       (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
##CV             452    351.9    353.2    355.0    352.8    348.4    343.6
##adjCV          452    351.6    352.7    354.4    352.1    347.6    342.7
##       7 comps  8 comps  9 comps  10 comps  11 comps  12 comps  13 comps
##CV       345.5    347.7    349.6     351.4     352.1     353.5     358.2
##adjCV    344.7    346.7    348.5     350.1     350.7     352.0     356.5
##       14 comps  15 comps  16 comps  17 comps  18 comps  19 comps
##CV        349.7     349.4     339.9     341.6     339.2     339.6
##adjCV     348.0     347.7     338.2     339.7     337.2     337.6
##
##TRAINING: % variance explained
##        1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps
##X         38.31    60.16    70.84    79.03    84.29    88.63    92.26
##Salary    40.63    41.58    42.17    43.22    44.90    46.48    46.69
##        8 comps  9 comps  10 comps  11 comps  12 comps  13 comps  14 comps
##X         94.96    96.28     97.26     97.98     98.65     99.15     99.47
##Salary    46.75    46.86     47.76     47.82     47.85     48.10     50.40
##        15 comps  16 comps  17 comps  18 comps  19 comps
##X          99.75     99.89     99.97     99.99    100.00
##Salary     50.55     53.01     53.85     54.61     54.61

validationplot(pcr.fit, val.type="MSEP")
  • 주성분회귀(PCR)은 pls라이브러리의 pcr함수를 사용하여 수행한다.
  • validation='CV'를 설정하면 pcr은 사용된 주성분의 수 M에 대한 10-fold 교차검증 오차를 계산한다.
  • 교차검증 결과는 validationplot 함수를 사용하여 그래프로 나타낼 수도 있다. val.type='MSEP'를 사용하면 교차검증 MSE가 그래프로 표현될 것이다.

  • M=16일 때 교차검증 오차가 가장 작다. 이것은 M=19일 때와 거의 차이가 없다. M=19이면 단순히 최소제곱을 수행하는 것이된다.
  • 하지만 그래프를 살펴보면 교차검증 오차는 하나의 성분만 포함하는 모델과 거의 같다는 것을 알 수 있다. 이것은 작은 수의 성분을 사용하는 모델이면 충분할 수 있다는 것을 시사한다.
set.seed(1)
pcr.fit=pcr(Salary~., data=Hitters, subsets=train, scale=TRUE, validation="CV")
validationplot(pcr.fit, val.type="MSEP")

  • 훈련셋에 대해 PCR을 수행하고 검정셋으로 성능을 평가한다. subset=train옵션으로 훈련셋의 인덱스를 지정해준다.
  • 교차검증 오차가 가장 낮은 것은 M=7개의 주성분이 사용된 경우이다.
pcr.pred=predict(pcr.fit, x[test,], ncomp=7)
mean((pcr.pred-y.test)^2)
##[1] 119516.5
  • 검정 MSE는 다음과 같이 계산한다.
  • PCR은 변수 선택을 수행하지 않고 심지어 직접적으로 계수 추정치도 제공하지 않기 때문에 모델 해석이 어렵다.
pcr.fit=pcr(y~x, xcale=TRUE, ncomp=7)
summary(pcr.fit)
##Data:   X dimension: 263 19 
##        Y dimension: 263 1
##Fit method: svdpc
##Number of components considered: 7
##TRAINING: % variance explained
##   1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps
##X    97.44    98.78    99.27    99.55    99.77    99.90    99.97
##y    28.15    37.22    38.23    43.61    43.70    43.77    48.07
  • 교차검증에 의해 선택된 주성분의 수 M=7을 사용하여 PCR을 전체 자료에 적합한다.

부분최소제곱(PLS)

set.seed(1)
pls.fit=plsr(Salary~., data=Hitters, subset=train, scale=TRUE, validation="CV")
summary(pls.fit)
##Data:   X dimension: 131 19 
##        Y dimension: 131 1
##Fit method: kernelpls
##Number of components considered: 19
##
##VALIDATION: RMSEP
##Cross-validated using 10 random segments.
##       (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
##CV           428.3    325.5    329.9    328.8    339.0    338.9    340.1
##adjCV        428.3    325.0    328.2    327.2    336.6    336.1    336.6
##       7 comps  8 comps  9 comps  10 comps  11 comps  12 comps  13 comps
##CV       339.0    347.1    346.4     343.4     341.5     345.4     356.4
##adjCV    336.2    343.4    342.8     340.2     338.3     341.8     351.1
##       14 comps  15 comps  16 comps  17 comps  18 comps  19 comps
##CV        348.4     349.1     350.0     344.2     344.5     345.0
##adjCV     344.2     345.0     345.9     340.4     340.6     341.1
##
##TRAINING: % variance explained
##        1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps
##X         39.13    48.80    60.09    75.07    78.58    81.12    88.21
##Salary    46.36    50.72    52.23    53.03    54.07    54.77    55.05
##        8 comps  9 comps  10 comps  11 comps  12 comps  13 comps  14 comps
##X         90.71    93.17     96.05     97.08     97.61     97.97     98.70
##Salary    55.66    55.95     56.12     56.47     56.68     57.37     57.76
##        15 comps  16 comps  17 comps  18 comps  19 comps
##X          99.12     99.61     99.70     99.95    100.00
##Salary     58.08     58.17     58.49     58.56     58.62

validationplot(pls.fit, val.type="MSEP")

  • 문법은 pcr과 동일하다.
  • 가장 낮은 교차검증 오차는 M=2개의 부분최소제곱 방향이 사용된 경우에 발생한다.
pls.pred=predict(pls.fit, x[test,], ncomp=2)
mean((pls.pred-y.test)^2)
##[1] 145367.7
  • 위의 결과에 대응하는 검정 MSE다.
pls.fit=plsr(Salary~., data=Hitters, scale=TRUE, ncomp=2)
summary(pls.fit)
##Data:   X dimension: 263 19 
##        Y dimension: 263 1
##Fit method: kernelpls
##Number of components considered: 2
##TRAINING: % variance explained
##        1 comps  2 comps
##X         38.08    51.03
##Salary    43.05    46.40
  • 부분최소제곱을 전체 자료에 적합했다.

'머신러닝, 딥러닝' 카테고리의 다른 글

[R] 사용되지 않는 level 제거하기 (feat. drop.levels함수)  (0) 2021.05.10
[머신러닝] 6. 선형성을 넘어서(비선형모델)- 다항식회귀, 계단함수, 스플라인, GAMs (feat. R Code)  (0) 2021.04.26
[머신러닝] 3. 분류- 로지스틱 회귀, LDA, QDA (feat. R Code)  (0) 2021.04.12
[머신러닝] 2. 선형회귀분석- 고려해야 할 요소, 잠재적 문제(feat. R Code)  (0) 2021.04.04
[머신러닝] 1. 통계학습  (0) 2021.04.02

'머신러닝, 딥러닝' Related Articles

티스토리툴바