Regularized regression

In this lab, we will see how to apply what we have learned about the ridge regression, lasso and elastic net. We will mainly use the glmnet package so please install it first.

1. Just to be sure to understand the math

In this section we will check if the maths related to regularized regression are properly understood. As always, we will fit a ridge regression, lasso and elastic net models on simulated data (so that the true model is perfectly known) from scratch and using the dedicated builtin R function.

a) Generate some artificial data

Let’s first simulate some data

set.seed(1)##fix the random seed so that we all have the same results
p <- 50
n <- 51
X <- matrix(rnorm(n * p), n)
beta.theo <- rep(0, p)
beta.theo[1:6] <- 2
Y.theo <- X %*% beta.theo
Y <- Y.theo + rnorm(n)

b) Analysis “from scratch”

1. Try to understand every lines of the code below (which computes the penalized least-square)

penalizedleastsquare <- function(par, features, response, lambda, alpha = 0){
  ## par is the vector of parameters, i.e., beta_0, beta_1
  ## features is the data.frame storing features
  ## response is the observed response (binary response here)
  ## lambda is the overall penalty coefficient
  ## alpha is the weight to balance between ridge and lasso penalties (default -> ridge)
  
  return(0.5 * mean((response - features %*% par)^2) + lambda * (alpha * sum(abs(par)) + 0.5 * (1 - alpha) * sum(par^2)))
}

2. Minimize the above function to get the ridge estimate when \(\lambda = 1\)

## Fill in this part

3. Redo the previous step with different values for \(\lambda\) and plot the evolution of the parameter estimate as \(\lambda\) increases.

## Fill in this part

4. Do the same with the lasso.

## Fill in this part

5. Do the same with an elastic net with \(\alpha = 0.5\).

## Fill in this part

c) Analysis using the dedicated R function

1. Read the help file of the glmnet function.
2. Redo the previous fits with this function and compare your results

## Fill in this part

4. Run the following code and comment

plot(fit)##fit is the output of the previous glm call, I'm here to explain each plot ;-)
plot(fit, xvar = "lambda", label = TRUE)

5. What does the following pieces of code

coef(fit, s = 1)
cvfit <- cv.glmnet(X, Y)
plot(cvfit)
coef(cvfit, s = "lambda.min")

6. Read (again) the documentation of the cv.glmnet function and do a 5-fold cross validation for tuning \(\lambda\).

2. Analyzing the Hitter dataset

The Hitter dataset stores 322 observation of major league Baseball players from the 1986 and 1987 season. The aim of this study is to be able to predict the Salary given some features. Here is a quick recap of the available feature and response.

• AtBat Number of times at bat in 1986
• Hits Number of hits in 1986
• HmRun Number of home runs in 1986
• Runs Number of runs in 1986
• RBI Number of runs batted in 1986
• Walks Number of walks in 1986
• Years Number of years in the major leagues
• CAtBat Number of times at bat during his career
• CHits Number of hits during his career
• CHmRun Number of home runs during his career
• CRuns Number of runs during his career
• CRBI Number of runs batted in during his career
• CWalks Number of walks during his career
• League A factor with levels A and N indicating player’s league at the end of 1986 Division A factor with levels E and W indicating player’s division at the end of 1986
• PutOuts Number of put outs in 1986
• Assists Number of assists in 1986
• Errors Number of errors in 1986
• Salary 1987 annual salary on opening day in thousands of dollars
• NewLeague A factor with levels A and N indicating player’s league at the beginning of 1987

You can get the data from here.

1. Import the data and familiarize your self with this dataset.
2. Fit a ridge regression model, plot the “Ridge path” and comment.
3. Select the best model according to a 5-fold cross validation procedure.
4. Redo the previous two steps with the lasso.
5. According to your model which features appear to be significant in predicting the salary?

3. Going a bit further

Actually the glmnet function can do much more than ridge, lasso and elastic net regression. For instance we can do a lasso logistic regression by passing the argument family = binomial to the glmnet function. Reanalyze the Default dataset we previously model with a penalized logistic regression.