glmnet

# subject : glmnet
# package 
#  * glmnet 

1. Model Selection ( Statistical Technicals )

# ref : [기본 교제] [ISL ch.6] http://www-bcf.usc.edu/~gareth/ISL/ISLR%20First%20Printing.pdf
# ref : [기본 교제] https://lagunita.stanford.edu/c4x/HumanitiesScience/StatLearning/asset/model_selection.pdf
# ref : [추가 교제] http://www.statground.org/?module=file&act=procFileDownload&file_srl=9145&sid=a0a4e7fd21503e45194d485ccb0c5476&module_srl=137
# ref : [참고]      http://sosal.kr/868
# ref : [참고]      https://rstudio-pubs-static.s3.amazonaws.com/22067_48fad02fb1a944e9a8fb1d56c55119ef.html

# Occam's Razor in Model Selection
#
#  - 통계 모델링에서 간결함의 원리
#     -- 복잡한 설명보다 단순하고 간단한 설명이 좋다.
#     -- 모형에 있어서 모수는 가능한 적게 가지고 있어야 한다.
#     -- 선형 모델이 비선형 모델보다 좋다.
#     -- 적은 수의 가정을 가진 실험은 많은 가정을 가진 실험보다 좋다.
#     -- 데이터로부터 회귀모델을 구할 때는 위해선 최소적합이 될 때까지 축소해야한다. 
#
#  - 단순화 작업 목록
#     -- 유의하지 않은 교호작용 항목을 제거한다.
#     -- 유의하지 않은 2차 또는 비선형 항목을 제거한다.
#     -- 유의하지 않은 설명 변수를 제거한다.
#     -- 유사한 모수값을 갖는 설명 변수를 통합한다.

1.1 Accuracy vs Explainability

# 1.1.1 Prediction Accuracy 
#  - n ≫ p ( n : # of observations, p : # of variables ) => the least squares estimates tend to also have low variance
#  - n > p ( n 이 충분히 크지 않은 경우 ) => can be a lot of variability, OVERFITTING 가능성이 높아짐 
#  - p > n => there is no longer a unique least squares coefficient estimate 
#  
#  - curse of dimensionality !
#
#  - 특정 n 값에 대해서 모델의 안정성을 높이기 위해서는, p 값을 줄이는 것이 효과가 있다 
#    ( 부가적으로 예측 성능을 높인다고 알려져 있다. )
#     -- constraining
#     -- shrinking ( Regularization ) : Lasso , Ridge , ElasticNet
#
# 1.1.2 Model Interpretability
#  - 당연하게도 설명변수의 수가 적으면, Interpretability 은 증가한다. 
#  - coefficient 값을 0으로 셋팅함으로써 linear model 에서 해당 feature 를 제거할 수 있다. ( feature selection )  
#
# ref : [관련 link] http://www.elderresearch.com/company/blog/predictive-model-accuracy-versus-interpretability

1.2 irrelevant feature 를 제거하는 방법

# 1.2.1 Subset Selection 
#  - like a regsubsets function ( estimate model based on AIC , BIC , ADJR2 , CP ..  ) 
#  - Stepwise Selection ( Backward / Forward ) or regsubsets ( http://rpubs.com/monoboy/glm 에서 수행해봤으므로 pass )
#
# 1.2.2 Shrinkage ( Regularization ) 
#  - Lasso , Ridge , ElasticNet
#  - shrinks the coefficient estimates towards zero
#
# 1.2.3 Dimension Reduction 
#  - PCA , ICA , t-SNE ( rpubs 로 정리는 안했지만, 살펴본 내용 pass ) 
# 
# 이번에는 Lasso , Ridge , ElasticNet 만 잘 살펴본다.
#
# 선형회귀모형에서 모형계수 벡터(coefficient estimated value) 에 OLS 보다
# 이것을 약간 축소한 ridge 해를 쓰는 것이 예측 성과가 좋은 것으로 알려져 있다. 
# 또한 모형계수 벡터의 일부요소를 0으로 퇴화시킨 lasso 해를 쓰게 되면 모형이 간결해진다는 잇점이 있다. 

Ridge Regression ( Linear Regression with L2 regularization, α=0 )

# 이 부분은 따로 정리하는 것보다, 그대로 보는 편이 더 정리가 될 것 같다.

# 일단 tunnig parameter (λ) 를 잘 선택해야한다. 
# 변화되는 λ 에 따라 coefficient set 이 각각 결정되며, 
# 우리는 MSE 혹은 error class ratio 등이 최소가되는 시점의 lamda 값을 찾아, 그때의 coef set 을 취한다.  
#
# RSS + Shrinkage Penalty 를 최소화시키는 것이 목표 
# lamda = 0 이면, RSS 를 최소화 시키는 것과 동일
# 이를 glmnet 의 sample data 중 response var. 이 binomial 인 것으로 수행해보자. 
library(glmnet)
data(BinomialExample)
plot(density(y))  # binomial 임.

selected_features <- as.matrix(as.data.frame(x)[c("V1", "V2", "V3", "V4", "V5")]) # 30개의 Feature 중 일부만 취하여 테스트한다.
model.ridge <- glmnet(selected_features, y, family = "binomial", alpha = 0 )      # binomial distribution 이므로 명시해준다. alpha = 0 은 ridge / 1 은 lasso 
# 기본적인 plot 은 위와같이 lamda / norm / dev 에 대해서 볼 수 있고, model 의 attributes 는 다음과 같다. 
attributes(model.ridge)

$names
 [1] "a0"         "beta"       "df"         "dim"        "lambda"     "dev.ratio"  "nulldev"    "npasses"    "jerr"       "offset"    
[11] "classnames" "call"       "nobs"      

$class
[1] "lognet" "glmnet"

# [df] 
#  - coef 가 non-zero 인 feature 의 갯수이다. 
#  - L2 regularization ( == ridge ) 은 coef 값을 조정은 하지만, 이를 0으로 만들지는 않는다. 따라서 5개의 feature 가 모두 살아남는다.
#  - L1 regularization ( == lasso ) 는 coef 값을 0까지 낮춰, 특정 feature 이 소거된다. 
#
# [dev.ratio] 
#  - percent of the null deviance explained 으로 null model 과의 차이를 말한다. 
#  
result <- cbind(df = model.ridge$df, dev.ratio = model.ridge$dev.ratio, lambda = model.ridge$lambda, log_lamda = log(model.ridge$lambda) )
head(result)

     df               dev.ratio   lambda log_lamda
[1,]  5 0.000000000000001618567 240.4762  5.482621
[2,]  5 0.000726975785632607077 219.1130  5.389588
[3,]  5 0.000797692900504678223 199.6476  5.296554
[4,]  5 0.000875254727849297766 181.9115  5.203520
[5,]  5 0.000960335450932397962 165.7509  5.110486
[6,]  5 0.001053659337566821910 151.0261  5.017453

options(encoding = 'utf-8')
par(family="AppleGothic")
Sys.setlocale("LC_ALL", "en_US.UTF-8") 

[1] "en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/C"

par(mfrow=c(2,2))
# model.ridge$beta 에 각 시도별 각 feature 에 대한 coefficient 값이 return 되며,
# lamda 값이 커질수록, 모든 feature 의 coef 값이 0에 수렴함을 알 수 있다. ( 이는 null model , y ~ B0 * intercept 와 동일 )
# lamda 값이 0 으로 가면, 일반 binomial glm 의 coef 값으로 수렴한다. 
# 각 feature 의 coef 값은 lamda 의 증가에 따라 증가 혹은 감소하는 것을 관찰할 수 있다. 
plot(model.ridge, xvar="lambda", label=TRUE, main = "lamda 와 각 feature 의 coef 관계") 
plot(model.ridge, xvar="norm", label=TRUE)

plot(model.ridge, xvar="dev", label=TRUE)

# 최적값을 자동으로 선택해주는 함수도 있다. 
# 최적의 lamda 를 어떤 기준으로 측정할 것인가를 type.measure argument 로 전달하면 된다.
# type.measure 로 선택할 수 있는 지표는 아래와 같다.
#  - “mse”, “deviance”, “class”, “auc”, “mae”
library(doMC)
registerDoMC(cores = 4)
par(mfrow = c(3,2))
plot(cv.glmnet(selected_features, y, alpha = 0, family = "binomial", nfolds = 10, type.measure = "mse"))
plot(cv.glmnet(selected_features, y, alpha = 0, family = "binomial", nfolds = 10, type.measure = "deviance"))

plot(cv.glmnet(selected_features, y, alpha = 0, family = "binomial", nfolds = 10, type.measure = "class"))
plot(cv.glmnet(selected_features, y, alpha = 0, family = "binomial", nfolds = 10, type.measure = "auc"))

plot(cv.glmnet(selected_features, y, alpha = 0, family = "binomial", nfolds = 10, type.measure = "mae"))

# 각 평가 기준에 따라, 우리는 lamda.min 과 lamda.1se 값에 해당하는 coef set 을 조회할 수 있다.
# 위 그래프에서, 두개의 점선이 있는데, 좌측이 min , 우측이 1se 값이다. 
# [lambda.min] minimum mean cross-validated error. 
# [lambda.1se] the most regularized model such that error is within one standard error of the minimum
model.ridge <- cv.glmnet(selected_features, y, alpha = 0, family = "binomial", nfolds = 10, type.measure = "mse")
c(model.ridge$lambda.min, model.ridge$lambda.1se)

[1] 0.03829062 0.27013273

c(coef(model.ridge, s = model.ridge$lambda.min), coef(model.ridge, s = model.ridge$lambda.1se))

[[1]]
6 x 1 sparse Matrix of class "dgCMatrix"
                     1
(Intercept)  0.2804697
V1           0.1538674
V2           0.4038223
V3          -0.2882205
V4          -0.8737910
V5          -0.3991743

[[2]]
6 x 1 sparse Matrix of class "dgCMatrix"
                      1
(Intercept)  0.24509256
V1           0.06364507
V2           0.21778061
V3          -0.14792249
V4          -0.41263055
V5          -0.22006191

# 자 길게 돌아왔다. 이제 predict 를 하자. 
# statistical technical 에 의해 결정된 lambda 값을 사용하여 prediction 후, 
# 서비스에서 살펴보려고 한 주제에 맞는 evaluation 방법으로 평가하자. 
p.min <- predict(model.ridge, newx = selected_features[1:5,], s = "lambda.min")
p.1se <- predict(model.ridge, newx = selected_features[1:5,], s = "lambda.1se")

Lasso Regression

# Ridge ( L2 regression ) 의 경우, feature 의 coef 값을 낮은 값으로 낮추기는 하지만, feature 를 제거하지는 못한다.
# Lasso ( L1 regression ) 의 경우, featrue 의 coef 값을 0 까지 낮춰, linear model 을 단순하게 만들어준다. 

model.lasso <- glmnet(selected_features, y, family = "binomial", alpha = 1 )
result <- cbind(df = model.lasso$df, dev.ratio = model.lasso$dev.ratio, lambda = model.lasso$lambda, log_lamda = log(model.lasso$lambda) )
# df 의 값이 0 까지 떨어지는 것을 볼 수 있다. ( df == 0 : null model )
head(result)

     df               dev.ratio    lambda log_lamda
[1,]  0 0.000000000000001618567 0.2404762 -1.425134
[2,]  1 0.029046784434478587628 0.2191130 -1.518168
[3,]  1 0.053350946639190058307 0.1996476 -1.611201
[4,]  1 0.073880814015632220726 0.1819115 -1.704235
[5,]  1 0.091340484596407103823 0.1657509 -1.797269
[6,]  1 0.106263675457181128170 0.1510261 -1.890303

par(mfrow=c(2,2))
plot(model.lasso, xvar="lambda", label=TRUE); plot(model.lasso, xvar="norm", label=TRUE); plot(model.lasso, xvar="dev", label=TRUE)

model.lasso <- cv.glmnet(selected_features, y, alpha = 1, family = "binomial", nfolds = 10, type.measure = "mse")
c(model.lasso$lambda.min, model.lasso$lambda.1se)

[1] 0.001313548 0.071749625

c(coef(model.lasso, s = model.lasso$lambda.min), coef(model.lasso, s = model.lasso$lambda.1se))

[[1]]
6 x 1 sparse Matrix of class "dgCMatrix"
                     1
(Intercept)  0.3032750
V1           0.2091808
V2           0.5103805
V3          -0.3676812
V4          -1.1793888
V5          -0.4897989

[[2]]
6 x 1 sparse Matrix of class "dgCMatrix"
                      1
(Intercept)  0.22778582
V1           .         
V2           0.13322369
V3          -0.03373276
V4          -0.65808493
V5          -0.20112684

ElasticNet Regression

# ref : https://web.stanford.edu/~hastie/TALKS/enet_talk.pdf ( Zou & Hastie .. EN 제안자들 자료 )
# ref : https://web.stanford.edu/~hastie/Papers/B67.2%20(2005)%20301-320%20Zou%20&%20Hastie.pdf
# ref : https://www.slideshare.net/ShangxuanZhang/ridge-regression-lasso-and-elastic-net
#
# The elastic-net penalty is controlled by α, and bridges the gap between lasso (α=1, the default) and ridge (α=0)
# 
# grouping effect ( related to collinearity )
#  - grouped selection
#    -- if two predictors are highly correlated among themselves, the estimated coef will be similar for them ( on ridge )
#    -- ridge is good for grouped selection                                                 => not good for eliminating trivial features
#    -- lasso select only one among them and remove the others ( shrink the other to zero ) => not good for grouped selection 
#
# R code 를 통한 simulation 은 생략한다. 
#
# feature 들간 다중공선성이 발생하는 경우, 이를 제거하는 방식의 차이가 Ridge / Lasso / EN 이 서로 다름을 인식하고.
# 결국에는 모두 try 를 해보고, 평가를 해봐야 할 것 같다.

기타

foldid <- sample(1:10,size=length(y),replace=TRUE)
cv.1  <- cv.glmnet(x,y,foldid=foldid,alpha=1)
cv.05 <- cv.glmnet(x,y,foldid=foldid,alpha=.5)
cv.0  <- cv.glmnet(x,y,foldid=foldid,alpha=0)

par(mfrow=c(2,2))
plot(cv.1);plot(cv.05);plot(cv.0)
plot(log(cv.1$lambda),cv.1$cvm,pch=19,col="red",xlab="log(Lambda)",ylab=cv.1$name)
points(log(cv.05$lambda),cv.05$cvm,pch=19,col="grey")
points(log(cv.0$lambda),cv.0$cvm,pch=19,col="blue")
legend("topleft",legend=c("alpha= 1","alpha= .5","alpha 0"),pch=19,col=c("red","grey","blue"))

# 주로 다중공선성이 발생되는 경우, redundant 한 feature 를 corrplot 을 보고 임의로 선택하는 경우가 많았는데, 
# Ridge / Lasso / Elastic-Net 을 시도해봄직할 것 같다.
#
# 서비스를 하다보면 의외로 다중공선성을 매우 자주 만난다. 
# 그도 그럴 것이, 로긴을 많이 할수록 글생산량도 높아지고, 글 생산량이 많은 사람이 친구도 많고, 뭐 이런 식이다. 

# 하지만 Random Forest 혹은 SVM 에서는 다중공선성이 Linear Regression 에서처럼 심각한 이슈가 아니다....

# ref : https://onlinecourses.science.psu.edu/stat857/node/155
# ref : https://ncss-wpengine.netdna-ssl.com/wp-content/themes/ncss/pdf/Procedures/NCSS/Ridge_Regression.pdf
# ref : https://lagunita.stanford.edu/c4x/HumanitiesScience/StatLearning/asset/model_selection.pdf
# ref : https://drsimonj.svbtle.com/ridge-regression-with-glmnet
# ref : https://gerardnico.com/wiki/lang/r/ridge_lasso#formula_parameters
# ref : http://r-bong.blogspot.kr/p/blog-page_4.html
# ref : https://web.stanford.edu/~hastie/TALKS/enet_talk.pdf
# 
# Motivation: too many predictors
# Motivation: ill-conditioned X
# Motivation: 다중공분산성 존재 해결