TDM 30200: Project 08 - Feature Selection (Regularization)

We’ve previously utilized the Batting table from the Lahman dataset, which contains batting statistics. Specifically, we examined the relationship between Homeruns (HR) and Runs Batted In (RBI), noting how HR influences RBI. However, the Batting table offers numerous other variables that can enhance our models. To explore these variables, execute the following R code first:

To gain detailed information about the variables in the Batting table, use the ?Batting command in R to access the documentation.

People is another table including player names, date of births and demographics for baseball players. It provides more details about players and it shares unique 'playerID' which can be used to merge Batting and People data sets. People contains 26 variables, but many of them are expected to have a large number of missing (NA) values. For instance, since the data includes records up to 2023, the death- columns will have many empty entries, as most players in recent years are alive. Additionally, these variables are not meaningfully related to RBI, so they do not need to be considered in the modeling process.

Some variables may provide interesting insights about the data but be less important for modeling. For example, while the age at which players make their debut can be an intriguing variable, their actual age might contribute more significantly to the model. Additionally, the most common birthplaces highlight Major League Baseball’s international reach, but its contribution to the model may be less. For curious readers, below are chart and tables related to debut ages and the most common (first ten) birthplaces. (Reproducing the visual and table is not required for the project.)

Prior to modeling, let’s examine a heatmap of the numeric variables intended for use in our model.

Data visualization offers essential clues to understanding your data before initiating the modeling process. For instance: At Bats (AB), Hits (H), and Doubles (X2B) are highly correlated because they follow a natural hierarchy in baseball statistics. AB represents a player’s batting opportunities, H is a subset of AB, counting successful hits, and X2B is a further subset, representing only doubles. Since more at-bats generally lead to more hits, and more hits increase the likelihood of doubles, these variables are inherently linked, resulting in strong correlations. In most cases, H (Hits) may be the best choice because it captures a player’s ability to reach base successfully, encompassing both singles and extra-base hits while avoiding redundancy.

The source of the picture is apnews accessed at 03/08/2025.

The following linear model uses as much variables as possible to explain the changes in target, RBI variable (The age variable was generated in Question 1.3).

Although a linear model provides P-values as evidence for variable selection, shrinkage or regularization (e.g., LASSO, Ridge) methods are used for variable selection because P-values can sometimes be misleading, especially when the sample size is large, or when variables are highly correlated with one another. In such cases, a variable with a high P-value might still be relevant for the model, but its contribution isn’t significant enough to pass the threshold. In some cases, especially when dealing with a large number of variables, instead of focusing on interpreting the output, you may want the model to perform both parameter estimation and variable selection simultaneously. Furthermore, many machine learning algorithms, due to their lack of interpretability, do not provide evidence (such as P-values) of which variables contribute to the model.

Variable selection (feature engineering) in statistical models is crucial for improving both the model’s performance and interpretability. By choosing only the most relevant variables, we can simplify the model, making it easier to understand and interpret, which is particularly important in fields where the relationships between variables are essential. Additionally, variable selection helps prevent overfitting, a common issue when too many irrelevant variables are included, which can lead to a model that fits the noise in the data rather than the underlying patterns. A more focused model with fewer predictors tends to generalize better to new data, leading to improved prediction accuracy. By selecting the right variables, we can also reduce computational costs, as fewer predictors mean less memory and processing power are required, which is especially important in large datasets.

While shrinkage methods continue to evolve and new ones (such as Elastic net) are introduced, understanding Ridge and LASSO provides a strong foundation for grasping more advanced techniques. Now, let’s recall the linear model we discussed earlier and review its mathematical representation with a single predictor.

$RBI_{i} = \beta_0 + \beta_1 HR_{i} + \epsilon_{i}$,

where:

We can replace the variable names in the model with their symbolic representations:

$y_{i} = \beta_0 + \beta_1 x_{i} + \epsilon_{i}$.

Generalize it for $p$ variables:

$y_{ij} = \beta_0 + \sum_{j = 1}^{p}\beta_j x_{ij} + \epsilon_{ij}$.

Since we are looking for the best estimates for the unknown $\beta$ parameters in this model, we want to minimize $\epsilon_{ij}$ as much as possible. The method that provides the best estimates for these parameters is the least squares approach, which finds $\beta$ values that minimize the residual sum of squares (RSS):

$RSS = \sum_{i = 1}^{n}(y_i - \beta_0 - \sum_{j = 1}^{p}\beta_j x_{ij})^2$.

Even without prior knowledge of this method, if we were to consider the best way to predict $y_i$ based on the given model, our natural goal would be to make both sides of the equation as close as possible. The variables in the model are known, and the parameters are estimated from the data, while $\epsilon$ represents an uncontrollable term arising from inherent uncertainty. Since this term cannot be eliminated, our objective is to minimize its impact. Rather than minimizing $\epsilon$ directly, we minimize its squared value to ensure that positive and negative errors do not cancel each other out.

Shrinkage/regularization methods introduce penalties when minimizing the Residual Sum of Squares (RSS), also known as the least squares loss function or objective function. These methods incorporate a penalty term into the RSS minimization process when estimating parameter values. Adjusting the magnitude of this penalty term influences the estimated parameters, helping to control model complexity and prevent overfitting.

Imagine you’re packing a suitcase for a trip. You want to bring everything you might need, but if you pack too much, your suitcase becomes heavy and difficult to carry. Shrinkage methods like Ridge and LASSO work similarly in a statistical model. Without any penalty, the model can include as many variables as possible, making it complex and potentially overfitting the data. However, by adding a penalty term (like an airline imposing a weight limit on luggage), the model is forced to prioritize important variables while reducing the impact of less significant ones. Ridge regression acts like a soft weight limit—allowing you to bring all your items but compressing them slightly to make the suitcase more manageable. LASSO, on the other hand, is stricter, forcing you to completely remove unnecessary items to meet the weight limit. This way, shrinkage methods prevent your model from becoming too complex while ensuring it still performs well.

Ridge regression adds a penalty term to the least squares objective function. This penalty term (L-2 penalty) is proportional to the squared magnitude of the coefficients, shrinking them towards zero:

$\sum_{i = 1}^{n}(y_i - \beta_0 - \sum_{j = 1}^{p}\beta_j x_{ij})^2 + \lambda \sum_{j = 1}^{p} \beta_j^2$

where $\lambda$ is a hyperparameter or it is called as tuning parameter.

Before moving on to the next paragraph, I recommend taking a minute to pause and consider what a hyperparameter is, how it differs from a parameter, and how it can be found in a general model setting.

As illustrated in the figure above, we can control the magnitude of the $\beta$ coefficients by adjusting the value of the tuning parameter, $\lambda$. During the minimization process, the Residual Sum of Squares (RSS) is minimized while simultaneously penalizing the size of the coefficients. This leads to the following relationships:

The programming languages such as R and Python estimates the parameter for different values of $\lambda$. glmnet is one of the libraries in R which can be used to run a regularization method. glmnet function can run several types of regression models with a grid of values for the regularization parameter, $\lambda$. Here is the example code:

The plot below illustrates how coefficient values vary with changes in log-lambda in Ridge regression.

While Ridge Regression effectively shrinks coefficients towards zero, it rarely eliminates them entirely. This means that even with Ridge, you might end up with a model containing many features, some of which may be redundant. To address this, we can use LASSO (Least Absolute Shrinkage and Selection Operator). LASSO, like Ridge, uses a penalty term to regularize the model, but it employs the L1 norm (absolute value) instead of the L2 norm (square). This crucial difference allows LASSO to perform feature selection by driving some coefficients to exactly zero, effectively removing those features from the model. This time, the penalized RSS looks as follows:

$\sum_{i = 1}^{n}(y_i - \beta_0 - \sum_{j = 1}^{p}\beta_j x_{ij})^2 + \lambda \sum_{j = 1}^{p} |\beta_j|$

Similar to ridge, we have a tuning parameter, $\lambda$ and changing its value changes the magnitude of parameter coefficients. R code is also exactly the same as the ridge regression, the only difference is that using alpha = 1 instead of 0 in the glmnet function.

The following plot shows how coefficient values are changing with the change in Log Lambda in LASSO regression.

Up to this point, we have run models but have not yet determined the optimized model or final parameter estimates. Why is that?

So far, we have not performed any feature selection; we have only applied regularization methods across different $\lambda$ values and reported the results. To identify the final model, we must use cross-validation to evaluate each $\lambda$ on unseen test data. The optimal $\lambda$ is then selected based on model performance. Only after this step can we determine which explanatory variables contribute most effectively to explaining the variation in the dependent variable.

With the coding and cross-validation knowledge we’ve gained so far, we have the tools to implement cross-validation for regularization in our models. However, R provides the cv.glmnet function, which streamlines this process with a concise set of code. If you feel confident in your understanding of shrinkage/regularization and cross-validation, you can use this function to perform feature selection. But if you think you need more practice with these concepts, I recommend implementing cross-validation manually before using cv.glmnet.

Let’s proceed with feature selection with cross-validation for LASSO:

R Code for Correlation Hep Map