#####################################################################################################
##
##                                    Ridge and Lasso Regression
##
## E. Flachaire, 30/11/2021
##
#####################################################################################################

library(ISLR)
Hitters=na.omit(Hitters)

## Ridge and Lasso 

library(glmnet)
x=model.matrix(Salary~.,Hitters)[,-1] # matrix of covariates (transform qualitative variable into dummies)
y=Hitters$Salary

ridge.model=glmnet(x,y,alpha=0)
lasso.model=glmnet(x,y,alpha=1)

par(mfrow=c(1,2))
plot(ridge.model, main="Ridge")
plot(lasso.model, main="Lasso")

## Cross-Validation

cv.ridge=cv.glmnet(x,y,alpha=0)
cv.lasso=cv.glmnet(x,y,alpha=1)

plot(cv.ridge)
plot(cv.lasso)

## Adaptive Lasso and Adaptive Elastic-net

ols=lm(Salary~.,Hitters)
w=1/abs(coef(ols))
cv.adalasso <- cv.glmnet(x,y,alpha=1,penalty.factor=w)
cv.adaelast <- cv.glmnet(x,y,alpha=.5,penalty.factor=w)
plot(cv.adalasso, main="Adaptive Lasso")
plot(cv.adaelast, main="Adaptive Elastic-net")

## Compare the coefficients

coef1=coef(ols)
coef2=coef(cv.ridge, s="lambda.min")
coef3=coef(cv.lasso, s="lambda.min")
coef4=coef(cv.lasso, s="lambda.1se")
coef5=coef(cv.adalasso, s="lambda.min")
coef6=coef(cv.adaelast, s="lambda.min")
options(scipen = 999) # disable scientific notation
coeff=cbind(coef1,coef2,coef3,coef4, coef5, coef6)
colnames(coeff) <- c("ols","ridge","lasso","lasso.1se","adaLasso", "adaElastic")
coeff





























######################################################################################################

summary(lm(Salary~., Hitters))
cbind(car::vif(ols))
###
set.seed(1)
train=sample(1:NROW(x),NROW(x)/2)
set.seed(1)
cv.out=cv.glmnet(x[train,],y[train],alpha=0)
plot(cv.out)
cv.out