In supervised learning, standardization is a pre-processing step where we rescale input variables to have zero mean and unit variance. There are several reasons why standardizing your variables before fitting your supervised learning model is a good idea:
Easier comparison of coefficients: For linear models, model coefficients for standardized features can be directly compared, since each coefficient represents the importance of a feature in terms of standard deviations. If the features are on different scales, comparing their coefficients does not make sense.
Improved numerical stability: Many machine learning algorithms tend to perform better if the features are on the same scale.
Prevents some features from dominating the model: In regularized regression models where we penalize the size of the coefficients, the scale of the input variables can greatly affect how much each coefficient gets penalized, which in turn affects the model fit. If some features have vastly different scales, they can overwhelmingly influence the model fit, which is typically undesirable behavior.
Are there situations where you might not want to rescale the input variables? If all the input variables are measured in the same units, then we typically do not want to standardize the variables. This is often the case with “omics” data such as gene expression measurements which are all measured on the same scale. Those genes exhibiting high standard deviation (fold change) are often the ones associated with the response, and by standardizing we effectively penalize them more.
5.1 Standardization in glmnet
The glmnet and cv.glmnet functions give the user the option of standardizing the input matrix x before fitting the regularized linear model via the standardize function argument. By default standardize = TRUE: all variables are rescaled to have mean 0 and variance 1 before model-fitting happens. To fit the model with x as is, set standardize = FALSE.
Note that whether standardize is set to TRUE or FALSE, the coefficients are always returned on the original scale. That means that if we want predictions for new observations, we do not have to standardize the input variables for these observations before feeding them to the prediction model. Also, when we interpret our coefficients, we should reason about them in their original units.
Let’s walk through an example of this to be clear. First, we create an x matrix where the first variable is on a much smaller scale than the other variables.
library(knitr)# create x matrix of 5 variablesset.seed(1)n <-100p <-5x <-matrix(rnorm(n * p), ncol = p)x[, 1] <- x[, 1] *0.01# x1 is on a much smaller scale than the otherskable(head(x))
-0.0062645
-0.6203667
0.4094018
0.8936737
1.0744410
0.0018364
0.0421159
1.6888733
-1.0472981
1.8956548
-0.0083563
-0.9109216
1.5865884
1.9713374
-0.6029973
0.0159528
0.1580288
-0.3309078
-0.3836321
-0.3908678
0.0032951
-0.6545846
-2.2852355
1.6541453
-0.4162220
-0.0082047
1.7672873
2.4976616
1.5122127
-0.3756574
The response variable is generated by \(y = 100x_1 + x_2 + x_3 + \text{noise}\):
We split the data into training and test sets, then fit a lasso model on the training set:
# split into training and test setsn_train <-95x_train <- x[1:n_train, ]y_train <- y[1:n_train]x_test <- x[(n_train +1):n, ]set.seed(1)cvfit <-cv.glmnet(x_train, y_train, standardize =TRUE)# lambda with smallest CV error and the coefficientscvfit$lambda.min#> [1] 0.08059172coef(cvfit, s ="lambda.min")#> 6 x 1 sparse Matrix of class "dgCMatrix"#> lambda.min#> (Intercept) -0.01254783#> V1 77.91103098#> V2 0.90374795#> V3 0.89177199#> V4 . #> V5 .
Based on the output of coef, the lasso model at the value of \(\lambda\) giving the smallest CV error is
\[y = -0.01 + 77.91x_1 + 0.90x_2 + 0.89x_3,\]
where the \(x\)’s are on the original scale, as in x_train. (We have rounded the coefficients to 2 decimal places for clarity.) Note that the returned coefficient for \(x_1\), 77.91, is on the same scale as the true coefficient for the original variable (100), even though \(x_1\) was scaled internally in glmnet for model-fitting.
When obtaining predictions on the test set, we don’t have to do any rescaling to the input matrix before passing it to the predict method:
Let’s dig deeper into the details of the previous example. When glmnet is called with standardize = TRUE, it standardizes the features in the x matrix to have mean 0 and variance 1. If feature \(x_i\) has estimated mean \(\mu_i\) and estimated standard deviation \(\sigma_i\), then the transformed feature is
\[z_i = \frac{x_i - \mu_i}{\sigma_i}.\]
Here is the equivalent transformation but in code. Note that the variance is estimated using the 1/n rather than the 1/(n-1) formula.
# get feature meansx_train_means <-colMeans(x_train)x_train_means#> [1] 0.001460855 -0.033568759 0.015458435 0.045226101 -0.015574477# get feature standard deviationsx_train_sds <-apply(x_train, 2, function(x) sd(x) *sqrt((n_train-1) / n_train))x_train_sds#> [1] 0.008887228 0.953384536 1.048073791 0.988229286 1.179332729# get transformed x matrixz_train <-scale(x_train, center = x_train_means, scale = x_train_sds)
If we are calling glmnet with family = "gaussian" (the default), then the response variable \(y\) is also standardized to have mean 0 and variance 1:
How can we compare unstd_fit and cvfit$glmnet.fit (the object returned from glmnet(x_train, y_train, standardize = TRUE))? First, let’s compare their lambda sequences. Since
we expect the lambda sequence for cvfit$glmnet.fit to be exactly the lambda sequence for unstd_fit multiplied by \(\sigma_y\). This is what we see (up to machine precision):
If we assign weights to our observations, the formulas for standardization are a bit different. In this case, glmnet will standardize the input features such that the weighted mean is 0 and the weighted variance is 1. (If we are calling glmnet with family = "gaussian" (the default), then the response variable \(y\) is also standardized to have weighted mean 0 and variance 1.) If observation \(i\) has weight \(w_i\), then for input feature \(x_j\) it is standardized using
The code below shows how we can compare the standardized and unstandardized fits when observation weights are involved. First, we generate some observation weights and fit a glmnet model with standardize = TRUE.
weights <-rep(1, n)weights[1:5] <-10w_std_fit <-glmnet(x, y, weights = weights)
Next, we standardize the input features and the response variable according to the weighted formulas, and fit a glmnet model with these and standardize = FALSE.
# get weighted feature means & standard deviationsx_w_means <-apply(x, 2, function(x) sum(weights * x) /sum(weights))x_w_vars <-sapply(1:p, function(j) sum(weights * (x[, j] - x_w_means[j])^2) /sum(weights))x_w_sds <-sqrt(x_w_vars)# get transformed x matrixz_w <-scale(x, center = x_w_means, scale = x_w_sds)# transform yy_w_mean <-sum(weights * y) /sum(weights)y_w_sd <-sqrt(sum(weights * (y - y_w_mean)^2) /sum(weights))y_prime <- (y - y_w_mean) / y_w_sd# fit model with standardize = FALSEw_unstd_fit <-glmnet(z_w, y_prime, weights = weights, standardize =FALSE)
The formulas for comparing the standardized and unstandardized fits in this setting are exactly the same as in Section 5.2, namely Equation 5.1 and Equation 5.2. First, let’s check that the lambda sequences match (up to machine precision):