Ridge_Lasso_ElasticNet.Rmd

---
title: "Second Model Build"
output:
  word_document: default
  html_document: default
---

In this second build, we will first do some subset selection.  First is a best subset analysis, followed by a forward stepwise analysis; both are run against the reduced dataset, “DE1_0_2008_Beneficiary_Summary_File_Sample_2_Reduced.csv”, which has been cleaned to remove the Beneficiary Identification field, Birth and Death Date fields (only the Age at either Death in 2008 or 12/31/2008 is included in the reduced dataset), State and County indicators, the number of months of Part A, Part B, Part D and HMO coverage for beneficiary, Beneficiary Responsibility Costs, and Primary Payer Reimbursement Costs.

Then Ridge Regression, Lasso Regression, and Elastic-Net Regression models will be created and tested.


```{r}
#First read in the reduced data file 
library(leaps)
synpuf <- read.csv("C:/Users/marze/OneDrive/Documents/UAlbany/CINF 624 FA21/SynPUF Data Files/DE1_0_2008_Beneficiary_Summary_File_Sample_2_Reduced.csv")
```

First we will do a Best Subset Analysis
credit for code to: https://www.statlearning.com/online-course
```{r}
#Do a preliminary subset analysis, first a full best subset regression
#default is max size 8, but we specify 16 to get all the variables included

regfit.full = regsubsets(MEDREIMB_TOTAL ~.,data=synpuf, nvmax=16)
summary(regfit.full)
reg.summary = summary(regfit.full)
names(reg.summary)
#Cp is an estimate of prediction error; we plot this estimate against the number of variables in the subset
#we want to choose the model with the lowest Cp
plot(reg.summary$cp,xlab="Number of Variables",ylab="Cp",main="Full Best Subset Regression")
min_Cp = which.min(reg.summary$cp)
points(min_Cp,reg.summary$cp[min_Cp],pch=20,col="red")  #make the minimum Cp a red dot
legend("topright",c("Minimum Cp"),pch=c(20),col="red") #make a legend on the plot for the red dot
#This is a plot of the Cp statistic (smallest is best) vs. the variables.  black squares are included in model
#bad Cps correspond to models with all the variables in them or fewer variables included
plot(regfit.full,scale="Cp",main="Best Subset Cp Plot")
#print out the best subset (one with the best Cp)
which.min(reg.summary$cp)
#print out the coefficients for the full subset regression
coef(regfit.full,min_Cp)
```

Forward Stepwise Selection
credit for code to: https://www.statlearning.com/online-course

```{r}


#now a forward-stepwise subset regression
regfit.fwd=regsubsets(MEDREIMB_TOTAL ~., nvmax=16,data=synpuf,nbest = 1, method="forward")
summary(regfit.fwd)
plot(regfit.fwd,scale="Cp",main="Forward Step Subset Cp Plot")
reg_fwd_summary = summary(regfit.fwd)
names(reg_fwd_summary)
plot(reg_fwd_summary$cp,xlab="Number of Variables",ylab="Cp",main="Forward Step Subset Regression")




```
Create a training and a validation set to assist in choosing a good subset model.  The response variable is MEDREIMB_TOTAL (total medicare reimbursements) modeled against all the 16 predictor variables
credit for code to: https://www.statlearning.com/online-course

```{r}
dim(synpuf)
set.seed(1)
train=sample(seq(37099),25970,replace=FALSE)
regfit.fwd=regsubsets(MEDREIMB_TOTAL ~., nvmax=16,data=synpuf[train,],method="forward")
#now make predictions on the observations not used for training. We have 16 models, so setup some 
#vectors to record the errors
val.errors=rep(NA,16)  #this creates the vector
x.test=model.matrix(MEDREIMB_TOTAL~.,data=synpuf[-train,]) #using the sample not in train
for(i in 1:16) {
  coefi=coef(regfit.fwd,id=i)
  pred=x.test[,names(coefi)]%*%coefi #indexing the names in x.test by using the coefi column names to get to      #the variables then multiplied by the coefficient vector (from lab in ch.6)
  val.errors[i]=mean((synpuf$MEDREIMB_TOTAL[-train]-pred)^2) #computing the mean squared error
}
plot(sqrt(val.errors),ylab="Root MSE",ylim=c(11400,11900),pch=19,type="b",main="Root MSE vs. Index")
#plot the MSE
points(sqrt(regfit.fwd$rss[-1]/25970),col="blue",pch=19,type="b") #rss on the training data on same graph
legend("topright",legend=c("Training","Validation"),col=c("blue","black"),pch=19)

plot(regfit.fwd$rss, xlab="Number of Variables", ylab="RSS", type="l")

#summary(regfit.fwd)
#plot(regfit.fwd,scale="Cp")
```


Model Selection by Cross-validation using 10-fold cross-validation.
credit for code to: https://www.statlearning.com/online-course

```{r}
library(leaps)
synpuf <- read.csv("C:/Users/marze/OneDrive/Documents/UAlbany/CINF 624 FA21/SynPUF Data Files/DE1_0_2008_Beneficiary_Summary_File_Sample_2_Reduced.csv")
set.seed(11)
folds=sample(rep(1:10,length=nrow(synpuf)))  #sampling by assigning each row in dataset a random # 1-10
#folds
table(folds)
cv.errors=matrix(NA,10,16)
for(k in 1:10){
  best.fit=regsubsets(MEDREIMB_TOTAL~.,data=synpuf[folds!=k,],nvmax=16,method="forward")
  for(i in 1:16) {
    pred=predict(best.fit,synpuf[folds==k,],id=i)
    cv.errors[k,i]=mean( (synpuf$MEDREIMB_TOTAL[folds==k]-pred)^2)
  }
}
rmse.cv=sqrt(apply(cv.errors,2,mean))
plot(rmse.cv,pch=16,type="b",main="10-fold Cross-validation Selection")
```

Ridge, Lasso, and Elastic-Net Regression

credit for code to:
https://www.pluralsight.com/guides/linear-lasso-and-ridge-regression-with-r


First split the data into training and testing data sets

```{r}
library(glmnet)
library(caret)
library(dplyr)

synpuf <- read.csv("C:/Users/marze/OneDrive/Documents/UAlbany/CINF 624 FA21/SynPUF Data Files/DE1_0_2008_Beneficiary_Summary_File_Sample_2_Diabetes_numeric_cleaned.csv")

#glmnet does not work with data frames, so we have to build a matrix for the training features, and
#a vector of target values to use for the response variable
library(glmnet)
synpuf <- read.csv("C:/Users/marze/OneDrive/Documents/UAlbany/CINF 624 FA21/SynPUF Data Files/DE1_0_2008_Beneficiary_Summary_File_Sample_2_Reduced.csv")
#train <- sample(seq(37099),25963,replace=FALSE)
index = sample(1:nrow(synpuf), 0.7*nrow(synpuf),replace=FALSE)
train = synpuf[index,]
test = synpuf[-index,]
x.train=model.matrix(MEDREIMB_TOTAL~.,data=train) #create a matrix for the training data
y.train=synpuf$MEDREIMB_TOTAL[index]   #create the response variables for the training data
x.test=model.matrix(MEDREIMB_TOTAL~.,data=test) #create a matrix for the testing data
y.test=synpuf$MEDREIMB_TOTAL[-index]   #create the response variables for the testing data




```

Ridge Regression

```{r}
#ridge

#the following command gives a list of lambdas, you have to pick the best one
ridge_reg = glmnet(x.train, y.train, nlambda=25, alpha=0, family = "gaussian")
plot(ridge_reg,main="Ridge Regression")
print(ridge_reg)
#the cv.glmnet command will automate finding the optimal lambda value
cv.ridge=cv.glmnet(x.train,y.train,alpha=0)
summary(cv.ridge)
plot(cv.ridge,main="Ridge Regression MSE vs Log Lambda",label = TRUE)  
#plot MSE vs. log(lambda)
plot(cv.ridge$glmnet.fit, "lambda", label=FALSE, main="Ridge coefficents vs log lamgda") 
#plot coefficients vs. log lambda
print(cv.ridge)

optimal_lambda <- cv.ridge$lambda.min #select the minimum lambda
optimal_lambda


#now that we have the optimal lambda, we can use that to re-build the ridge regression model

best_ridge <- glmnet(x.train, y.train, alpha=0, lambda = optimal_lambda)
coef(best_ridge)
summary(best_ridge)


#first create a function for model performance metrics, used for all the following regression models
eval_results <- function(true, predicted, df) {
  SSE <- sum((predicted - true)^2)
  SST <- sum((true - mean(true))^2)
  SSE
  SST
  R_square <- 1 - SSE / SST
  RMSE = sqrt(SSE/nrow(df))

  data.frame(RMSE = RMSE, Rsquare = R_square)

}


#prediction and evaluation on training data set
predictions_train <- predict(best_ridge, s = optimal_lambda, newx = x.train)
eval_results(y.train, predictions_train, train)


#prediction and evaluation on test data set
predictions_test <- predict(best_ridge, s = optimal_lambda, newx = x.test)
eval_results(y.test, predictions_test, test)







```



Now we try Lasso Regression

```{r}
#lasso
x.train=model.matrix(MEDREIMB_TOTAL~.,data=train) #create a matrix for the training data
y.train=synpuf$MEDREIMB_TOTAL[index]   #create the response variables for the training data
x.test=model.matrix(MEDREIMB_TOTAL~.,data=test) #create a matrix for the testing data
y.test=synpuf$MEDREIMB_TOTAL[-index]   #create the response variables for the testing data

#Use the earlier train/validation sets to select the 'lambda' for the lasso
lasso_reg=cv.glmnet(x.train,y.train, alpha=1, standardize=TRUE,nfolds=5)
plot(lasso_reg)
coef(lasso_reg)

#choose the best lambda
lambda_best <- lasso_reg$lambda.min
lambda_best

#train the lasso model with the chosen lambda
lasso.model=glmnet(x.train,y.train, alpha=1, lambda = lambda_best, standardize=TRUE)
pred_l_train <- predict(lasso.model, s= lambda_best, newx=x.train)  #use the non-training data as test data
eval_results(y.train, pred_l_train, train)

pred_l_test <- predict(lasso.model, s = lambda_best, newx = x.test)
eval_results(y.test, pred_l_test, test)


```

Elastic-Net Regression

```{r}
#elastic-net regression

x.train=model.matrix(MEDREIMB_TOTAL~.,data=train) #create a matrix for the training data
y.train=synpuf$MEDREIMB_TOTAL[index]   #create the response variables for the training data
x.test=model.matrix(MEDREIMB_TOTAL~.,data=test) #create a matrix for the testing data
y.test=synpuf$MEDREIMB_TOTAL[-index]   #create the response variables for the testing data

library(caret)
library(dplyr)
#trying an alternative method to generate train and test matrices
dummy_dataset<- dummyVars(MEDREIMB_TOTAL~., data = synpuf)
train_dummies = predict(dummy_dataset, newdata = train)
test_dummies = predict(dummy_dataset, newdata = test)

elastic_x = as.matrix(train_dummies)  #example's x
elastic_y_train = train$MEDREIMB_TOTAL #example's y_train
elastic_x_test = as.matrix(test_dummies) #examples x_test
elastic_y_test = test$MEDREIMB_TOTAL #example's y_test

#setup the training control
train_cont <- trainControl(method = "repeatedcv",
                           number = 10,
                           repeats = 5,
                           search = "random",
                           verboseIter = TRUE)

#train the model
elastic_reg <- train(MEDREIMB_TOTAL~.,
                     data = train,
                     method = "glmnet",
                     preProcess = c("center", "scale"),
                     tuneLength = 10,
                     trControl = train_cont)

# Best tuning parameter
elastic_reg$bestTune

#Make predictions on training set
predictions_ElNet_train <- predict(elastic_reg, newdata = train)
eval_results(elastic_y_train, predictions_ElNet_train, train_dummies)
#Make predictions on testing set
predictions_ElNet_test <- predict(elastic_reg, newdata = test)
eval_results(elastic_y_test, predictions_ElNet_test, test_dummies)

plot(elastic_reg)

```