24 glmnet’s limitations

While we’ve spent several chapters extolling the benefits of the glmnet R package, we are well aware that glmnet might not meet all of your modeling needs. In this chapter, we discuss some of glmnet’s limitations and recommend alternatives.

24.1 Problem size

The main limitation faced by glmnet users in R is that the x and y inputs have to be loaded in R’s memory before glmnet can use them. This means that the amount of training data we can give glmnet is bottle- necked by the size of R’s memory. [TODO: Add some rough limits we have seen here.]

Here are some possible ways to work around this issue:

Use pre-screening as a way to decrease the number of input features for the model. (Just make sure to do pre-screening correctly: see Section 18.2 for a discussion.)

24.2 Elastic net penalty

glmnet fits elastic net models really well, but that’s all it does. The elastic net gives sparsity of input features for the regression model, but there might be other types of sparsity that we want from our regularized regression. In Chapter 26 we discuss some of these other regression models and recommend implementations for them.

24.3 Linear relationship with input features

In the family = "gaussian" context, if the input features are \(x_1, \dots, x_p\), then the model prediction is forced to be of the form

\[ y = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p. \]

This forces each input feature to have a linear relationship with the response, and there are no interactions between any two input features.

One way to improve the expressiveness of the model is to add new features which are derived from the original \(p\) features. For example, if we want to make the relationship between \(y\) and \(x_1\) more flexible, we could add \(x_{p+1} = x_1^2\) and \(x_{p+2} = x_1^3\) to capture quadratic and cubic relationships between \(y\) and \(x_1\).

Another common approach is to expand each feature into multiple features using a spline basis. The code below demonstrates how you can expand every feature in an input matrix with a B-spline basis.

# simulated data example
data("QuickStartExample")
train_x <- QuickStartExample$x[1:90,]
train_y <- QuickStartExample$y[1:90]
test_x <- QuickStartExample$x[91:100,]
test_y <- QuickStartExample$y[91:100]

# original train x matrix has 20 features
p <- ncol(train_x)
p
#> [1] 20

# create transformed input matrix
# df is the number of features each original features expands into
library(splines)
df <- 3
spline_basis_list <- lapply(1:p, function(j) {
  bs(train_x[, j], df = df)
})
transformed_train_x <- do.call(cbind, spline_basis_list)

# transformed train x has 20 * 3 = 60 features
ncol(transformed_train_x)
#> [1] 60

# fit model on the transformed input matrix
set.seed(1)
cvfit <- cv.glmnet(transformed_train_x, train_y)

Note that we have to save the spline objects from the training data in order to transform the input matrix for test data correctly:

# get transformed test x
transformed_test_x <- do.call(cbind, lapply(1:p, function(j) {
  predict(spline_basis_list[[j]], newx = test_x[, j])
}))

# transformed test x has 60 features as well
ncol(transformed_test_x)
#> [1] 60

# get predictions at lambda.min
test_y_pred <- predict(cvfit, newx = transformed_test_x, s = "lambda.min")
test_y_pred
#>       lambda.min
#>  [1,] -2.4270384
#>  [2,] -0.8443221
#>  [3,] -0.4033902
#>  [4,] -1.9719815
#>  [5,]  1.3264066
#>  [6,] -3.9959630
#>  [7,] -1.5715346
#>  [8,] -2.9605005
#>  [9,]  5.0397046
#> [10,] -0.5040410

There are many other possible methods for creating these new features: see for example Chapter 5 of Hastie et al. (2001).

To add interactions between the input features, we could add pairwise products \(x_i x_j\) (with \(i \neq j\)) as new features. This will greatly increase the total number of features from \(p\) to \(p^2 - p\). This may still be computationally feasible for small to moderate \(p\). For large \(p\) we will need more sophisticated methods, some of which we outline in Chapter 26.