Tutorials › Exam PA - Predictive Analytics › Fitting Ridge, Lasso, And Elastic Net In R

Exam PA · Tuning and Regularized Model Use · Tutorial 461 of 966

Fitting Ridge, Lasso, And Elastic Net In R

Fit penalized regression models in R, select the family and alpha deliberately, and interpret fitted values and basic diagnostics without confusing link-scale output with response-scale predictions.

Advanced 12 min read

What You'll Learn

  • Prepare a numeric predictor matrix and response for glmnet
  • Specify Gaussian or binomial family and ridge, lasso, or elastic-net alpha
  • Explain the lambda path and key contents of a glmnet fitted object
  • Request predictions on the link, response, or class scale as appropriate
  • Calculate RMSE, Brier score, and threshold-based classification accuracy

From a Prepared Design Matrix to a Fitted Model

The earlier material on preprocessing and model matrices addresses how to create predictors without leaking information from validation data. Once the training predictors and outcome are ready, the next task is to fit a penalized regression model correctly. In R, the glmnet package fits ridge, lasso, and elastic-net models using a common interface. The important choices are the outcome family, the mixing parameter \(\alpha\), and the penalty value \(\lambda\).

This tutorial focuses on fitting Gaussian models for quantitative outcomes and binomial models for binary outcomes. It also explains what the fitted object contains, how to obtain predictions in the desired form, and how to calculate a few basic diagnostics. A fitted model is not automatically a validated or well-calibrated model: these calculations are summaries of predictions on specified data, not evidence by themselves that a model will perform well on new cases.

Definition: In glmnet, the predictor input \(x\) is a numeric matrix with one row per observation and one column per predictor. The response \(y\) is a vector with one entry per row. The family specifies the outcome model: use "gaussian" for a quantitative response and "binomial" for a two-class response. The parameter \(\alpha\) sets the penalty mixture: \(\alpha=0\) is ridge, \(\alpha=1\) is lasso, and \(0<\alpha<1\) is elastic net.

The Model and Its Penalty

For a Gaussian model, the fitted conditional mean is linear in the predictors. With the usual independent, constant-variance Gaussian error model, glmnet minimizes a squared-error criterion plus a penalty. For a binomial model with a binary response coded 0 or 1, it models the probability of outcome 1 using a logistic link and minimizes binomial negative log-likelihood plus a penalty. These likelihood descriptions specify the model; the penalty changes how the coefficients are estimated.

In the following objectives, \(N\) is the number of observations, \(\beta_0\) is the intercept, and \(\boldsymbol{\beta}\) is the vector of predictor coefficients. The predictor matrix is understood to be the matrix used in the optimization; with glmnet's default standardization, predictors are standardized internally. The intercept is not penalized.

$$ \text{Gaussian:}\qquad \frac{1}{2N}\sum_{i=1}^{N}(y_i-\beta_0-\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta})^2 +\lambda\left[\frac{1-\alpha}{2}\|\boldsymbol{\beta}\|_2^2 +\alpha\|\boldsymbol{\beta}\|_1\right] $$
$$ \text{Binomial:}\qquad -\frac{1}{N}\sum_{i=1}^{N} \left[y_i\log(p_i)+(1-y_i)\log(1-p_i)\right] +\lambda\left[\frac{1-\alpha}{2}\|\boldsymbol{\beta}\|_2^2 +\alpha\|\boldsymbol{\beta}\|_1\right], \qquad \log\left(\frac{p_i}{1-p_i}\right)=\beta_0+\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta} $$

The objective shows why \(\alpha\) and \(\lambda\) play different roles. \(\alpha\) sets the relative contribution of the \(L_1\) and squared \(L_2\) penalties; \(\lambda\) controls the overall penalty strength. At \(\alpha=0\), the \(L_1\) term disappears, giving ridge regression. At \(\alpha=1\), the squared \(L_2\) term disappears, giving lasso. Intermediate values combine them. Larger \(\lambda\) generally shrinks coefficients more strongly. Coefficients and their selected predictors should not be interpreted as causal effects merely because they appear in a fitted model.

Fitting a Model and Reading the Fitted Object

A glmnet() call ordinarily fits a sequence of models across a decreasing sequence of penalty values, rather than just one model. If a particular lambda sequence is supplied, the returned path corresponds to that sequence. The fitted object is commonly of class "glmnet". It stores the lambda values and coefficient solutions along the path, along with information about the fit. The coefficient matrix is often stored sparsely; coef(fit, s=...) requests coefficients at a specified lambda. A requested value between values in the stored path may be interpolated.

The intercept is handled by the model fit, so do not add a column of ones to \(x\) when using the default intercept=TRUE. With the default standardize=TRUE, glmnet standardizes predictor columns for fitting and returns coefficients on the original predictor scale. For a prediction call, newx must have the same predictor columns, in the same order, as the matrix used for fitting. Any factor encoding or other preprocessing must also match the training design.

Worked Example: Fitting a Gaussian Elastic-Net Model

A hypothetical analyst has six training records, each with two numeric predictors: a device's monthly operating hours and its number of prior service alerts. The quantitative outcome is repair cost in hundreds of dollars. The analyst wants to fit an elastic-net Gaussian model at \(\lambda=0.20\), then predict for a device with 20 operating hours and one alert. Identify the fitting choices and show how to use the fitted coefficients.

Method. Arrange the two predictors as columns of a numeric matrix with six rows and put the six repair costs in a response vector \(y\). The response is quantitative, so use family="gaussian". Elastic net requires a value strictly between 0 and 1; choose \(\alpha=0.5\). The R call is fit <- glmnet(x, y, family="gaussian", alpha=0.5, lambda=0.20). The matrix has \(N=6\) rows and \(p=2\) predictor columns, and the specified fit uses \(\lambda=0.20\).

Suppose that, for illustration, extracting the coefficients at this lambda gives an intercept of 4.2, an operating-hours coefficient of 0.30, and an alerts coefficient of \(-0.50\). These illustrative values show how to use the returned coefficients; they are not asserted to be the numerical result of fitting any particular six observations. For the new device, the linear predictor is \(4.2+(0.30)(20)+(-0.50)(1)\). The contributions are \(4.2\), \(6.0\), and \(-0.5\), so the predicted outcome is \(4.2+6.0-0.5=9.7\) hundreds of dollars, or \(\$970\).

Check. The prediction is on the original outcome scale because this is a Gaussian fit. A matching R prediction call uses predict(fit, newx=newx, s=0.20, type="response"), where newx has one row and the same two columns in the same order. A common mistake is to treat \(\alpha\) as the strength of shrinkage; here \(\alpha=0.5\) chooses the mixture, while \(\lambda=0.20\) sets its overall strength.

A short justification explains the link-scale and probability-scale predictions for a binomial fit. The model states that \(\eta=\log(p/(1-p))\). Exponentiating gives \(e^\eta=p/(1-p)\), so \(e^\eta(1-p)=p\), or \(e^\eta=p(1+e^\eta)\). Solving for \(p\) gives \(p=e^\eta/(1+e^\eta)=1/(1+e^{-\eta})\). Thus a link prediction is a log-odds value, while a response prediction is a probability between 0 and 1.

Choosing the Prediction Type

The predict() method can return different quantities. For a Gaussian fit, type="response" returns predicted outcome values; in this setting, the linear predictor and response prediction are on the same scale. For a binomial fit, type="link" returns log odds and type="response" returns estimated probabilities. For binomial classification, type="class" returns classes using a threshold of 0.5. When a different decision threshold is required, request probabilities and apply that threshold explicitly. The threshold is a classification rule, not a change to the fitted probability model.

Worked Example: From Binomial Coefficients to a Probability

A hypothetical service team models whether a device requires emergency review in the next month. The response is coded 1 for review and 0 otherwise. There are two numeric predictors: days since last inspection and a binary indicator for a reported warning. The analyst wants a lasso fit at \(\lambda=0.05\). For a worked prediction, suppose the illustrative fitted coefficients at that lambda are intercept \(-2.0\), days coefficient \(0.04\), and warning coefficient \(0.70\). Predict for a device with 30 days since inspection and a warning.

Fit choices. The response is binary, so use family="binomial". Lasso uses \(\alpha=1\). With a two-column numeric matrix \(x\) and a 0/1 response vector \(y\), the call is fit <- glmnet(x, y, family="binomial", alpha=1, lambda=0.05). The fitted coefficient values in this example are stipulated for illustrating prediction, not claimed as the numerical output from an unspecified training data set.

Link prediction. Substitute the predictor values into the linear predictor: \(\eta=-2.0+(0.04)(30)+(0.70)(1)\). The days contribution is \(1.20\), the warning contribution is \(0.70\), and therefore \(\eta=-2.0+1.20+0.70=-0.10\). The link-scale prediction is \(-0.10\) log odds.

Probability prediction. Convert the link value using the inverse-logit formula: \(p=1/(1+e^{-(-0.10)})=1/(1+e^{0.10})\). Since \(e^{0.10}\approx1.1052\), the denominator is \(2.1052\), giving \(p\approx0.4750\) to four decimal places. At a 0.5 threshold the predicted class is 0, because \(0.4750<0.5\).

Check. The probability is below 0.5 and the class is 0, even though the log odds are only slightly negative. The common error is to interpret \(-0.10\), the link prediction, as a probability. A prediction call with type="response" returns approximately 0.4750; type="link" returns \(-0.10\).

Basic Diagnostics on Predictions

A diagnostic calculation requires observed outcomes and predictions on a compatible scale. For quantitative outcomes, the root mean squared error (RMSE) is the square root of the average squared prediction error. For binary outcomes, the Brier score is the average squared difference between the observed 0/1 outcome and the predicted probability. Both scores are smaller when predictions are closer to their observed values, but they summarize different outcome types. Classification accuracy instead compares predicted classes with observed classes at a stated threshold.

$$ \operatorname{RMSE} =\sqrt{\frac{1}{N}\sum_{i=1}^{N}(y_i-\hat{y}_i)^2}, \qquad \operatorname{Brier} =\frac{1}{N}\sum_{i=1}^{N}(y_i-\hat{p}_i)^2 $$

These metrics are meaningful only when the predictions are evaluated on data appropriate to the question, such as a held-out set for a predictive assessment. Do not calculate a binary Brier score from predicted classes: it uses probabilities. Likewise, do not calculate RMSE from binomial link predictions when the task is to assess probabilities or classes. A confusion matrix tabulates observed and predicted classes, so its counts depend on the chosen threshold.

Worked Example: Checking Binary Predictions and Accuracy

A hypothetical model produces probabilities \(0.8\), \(0.3\), \(0.4\), and \(0.6\) for four validation records. The observed outcomes in the same order are \(1\), \(0\), \(1\), and \(0\). Calculate the Brier score and the number and proportion classified correctly at threshold 0.5.

Brier score. The squared errors are \((1-0.8)^2=0.2^2=0.04\), \((0-0.3)^2=(-0.3)^2=0.09\), \((1-0.4)^2=0.6^2=0.36\), and \((0-0.6)^2=(-0.6)^2=0.36\). Their sum is \(0.04+0.09+0.36+0.36=0.85\). Divide by \(N=4\): the Brier score is \(0.85/4=0.2125\).

Classes and accuracy. With threshold 0.5, the probabilities give predicted classes \(1,0,0,1\), respectively. Comparing with observed classes \(1,0,1,0\), the first two are correct and the last two are incorrect. The count correct is 2 of 4, and the accuracy proportion is \(2/4=0.5000\).

Check. The Brier score averages four squared probability errors, while accuracy counts correct thresholded classes. A common mistake is to report the count 2 as an accuracy proportion; the correct proportion is 0.5000. Changing the threshold could change the classes and accuracy without changing the model's four predicted probabilities or their Brier score.

Common Fitting and Prediction Traps

  • Using the wrong family. A continuous cost outcome calls for a Gaussian model; a 0/1 outcome calls for a binomial model. The family is not selected by whether predictors are numeric or categorical.
  • Reversing alpha and lambda. Alpha specifies the ridge–lasso mixture, with 0 for ridge and 1 for lasso. Lambda controls the strength of the penalty.
  • Sending incompatible new data. The prediction matrix must have the same columns, encoding, and order as the training matrix. A matrix with the right number of columns but a different order can produce incorrect predictions without an obvious error message.
  • Confusing log odds and probability. For binomial models, use type="response" when probabilities are required. A link-scale value can be negative or greater than 1; it is not a probability.
  • Mixing up a classification count and a proportion. If 2 of 3 records are classified correctly, the count is 2 and the accuracy proportion is \(2/3\approx0.6667\). State which quantity is being reported.
  • Assuming the fitted path selects the best lambda. A call to glmnet() fits coefficient solutions along a lambda path. The existence of a path does not establish which lambda performs best for a predictive objective.
Key takeaway: Prepare a consistent numeric matrix, match the family to the response, set \(\alpha\) for the desired penalty mixture, and use a clearly specified \(\lambda\). For predictions, distinguish the linear predictor, response scale, and classification threshold. These fitting and interpretation practices are aligned with the SOA Exam PA (Predictive Analytics) syllabus.

Check Your Understanding

For each problem, identify the relevant fitting or prediction choice and check the arithmetic before comparing with the answer.

  1. A training set has 90 observations, five numeric predictor columns, and a quantitative response. State the dimensions of \(x\), the family, and the alpha value for a ridge fit. (Answer: \(x\) is \(90\times5\); use family="gaussian" and \(\alpha=0\).)
  2. A binomial fit gives \(\eta=\log(3)\) for a record. Convert to a predicted probability and classify at threshold 0.5. (Answer: \(p=3/(1+3)=0.7500\), class 1.)
  3. A team wants an elastic-net binomial model with \(\lambda=0.08\), and its design matrix includes four predictor columns. State the appropriate family, alpha range, and the matrix shape if there are 240 training records. (Answer: use family="binomial", choose \(0<\alpha<1\), and use a \(240\times4\) matrix.)
  4. Observed binary outcomes are \(1,0,1\), and predicted probabilities are \(0.8,0.3,0.4\). Calculate Brier score and accuracy proportion at threshold 0.5. (Answer: squared errors are \(0.04,0.09,0.36\), so Brier score is \(0.49/3\approx0.1633\); predicted classes are \(1,0,0\), so accuracy is \(2/3\approx0.6667\).)
  5. For three quantitative outcomes, observed values are \(10,12,15\) and predictions are \(9,14,14\). Calculate RMSE. (Answer: squared errors are \(1,4,1\), their mean is \(6/3=2\), and RMSE is \(\sqrt{2}\approx1.4142\).)
Clinical Biostats