# Hitters data - lasso regression

library(ISLR)
names(Hitters)
dim(Hitters)
sum(is.na(Hitters$Salary))
Hitters=na.omit(Hitters)
dim(Hitters)
sum(is.na(Hitters))


Hitters2 = Hitters
#Scale numerical variables
for(j in 1:20)
{
 if((names(Hitters)[j]!="Salary") & (class(Hitters[,j])!="factor"))
    Hitters2[,j] = scale(Hitters2[,j])
}
fit.lm = lm(Salary~.,data=Hitters2)

#Make design matrix
x=model.matrix(Salary~.,Hitters2)[,-1]
y=Hitters2$Salary

library(glmnet)
#Lambda values (s in glmnet)
grid=10^seq(5,-2,length=100)
#Fit with lasso regression, alpha=1
lasso.mod=glmnet(x,y,alpha=1,lambda=grid)
#Plotting estimates for different values of lambda/s
plot(lasso.mod,xvar="lambda")
dim(coef(lasso.mod))
lasso.mod$lambda[50]
coef(lasso.mod)[,50]
sqrt(sum(coef(lasso.mod)[-1,50]^2))
predict(lasso.mod,s=grid[50],type="coefficients")[1:20,]

#Dividing into training/test set
set.seed(1)
train=sample(1:nrow(x), nrow(x)/2)
test=(-train)
y.test=y[test]
lasso.mod=glmnet(x[train,],y[train],alpha=1,lambda=grid, thresh=1e-12)

#Prediction error for specific value of s
lasso.pred=predict(lasso.mod,s=4,newx=x[test,])
sqrt(mean((lasso.pred-y.test)^2))
#Compare with just using the mean.
sqrt(mean((mean(y[train])-y.test)^2))
#Calculate RSS for each value of lambda
rss = rep(NA,length(grid))
for(k in 1:length(grid))
{
 lasso.pred=predict(lasso.mod,s=grid[k],newx=x[test,])
 rss[k] = sqrt(mean((lasso.pred-y.test)^2))
}
plot(log(grid),rss,type="l",xlab="log(lambda)")

#Cross-validation
set.seed(1)
cv.out=cv.glmnet(x[train,],y[train],alpha=1)
pdf("../fig/Hitters_lasso_cv.pdf",height=5,width=8)
plot(cv.out)
dev.off()
cv.out


plot(lasso.mod,xvar="lambda")
abline(v=log(cv.out$lambda.min))
abline(v=log(cv.out$lambda.1se),col=2)

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

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