<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/a-deep-dive-into-glmnet-standardize/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://statisticaloddsandends.wordpress.com/2018/11/15/a-deep-dive-into-glmnet-standardize/"> R – Statistical Odds &amp; Ends</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>I’m writing a series of posts on various function options of the <code>glmnet</code> function (from the package of the same name), hoping to give more detail and insight beyond R’s documentation.</p>
<p>In this post, we will focus on the <strong>standardize</strong> option.</p>
<p>For reference, here is the full signature of the <code>glmnet</code> function:</p>
<pre>
glmnet(x, y, family=c(&quot;gaussian&quot;,&quot;binomial&quot;,&quot;poisson&quot;,&quot;multinomial&quot;,&quot;cox&quot;,&quot;mgaussian&quot;),
    weights, offset=NULL, alpha = 1, nlambda = 100,
    lambda.min.ratio = ifelse(nobs&lt;nvars,0.01,0.0001), lambda=NULL,
    standardize = TRUE, intercept=TRUE, thresh = 1e-07,  dfmax = nvars + 1,
    pmax = min(dfmax * 2+20, nvars), exclude, penalty.factor = rep(1, nvars),
    lower.limits=-Inf, upper.limits=Inf, maxit=100000,
    type.gaussian=ifelse(nvars&lt;500,&quot;covariance&quot;,&quot;naive&quot;),
    type.logistic=c(&quot;Newton&quot;,&quot;modified.Newton&quot;),
    standardize.response=FALSE, type.multinomial=c(&quot;ungrouped&quot;,&quot;grouped&quot;))
</pre>
<p>Unless otherwise stated, <img alt="n" title="n" src="https://s0.wp.com/latex.php?latex=n&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> will denote the number of observations, <img alt="p" title="p" src="https://s0.wp.com/latex.php?latex=p&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> will denote the number of features, and <code>fit</code> will denote the output/result of the <code>glmnet</code> call. The data matrix is denoted by <img alt="X in mathbb{R}^{n times p}" title="X in mathbb{R}^{n times p}" src="https://s0.wp.com/latex.php?latex=X+%5Cin+%5Cmathbb%7BR%7D%5E%7Bn+%5Ctimes+p%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> and the response is denoted by <img alt="y in mathbb{R}^n" title="y in mathbb{R}^n" src="https://s0.wp.com/latex.php?latex=y+%5Cin+%5Cmathbb%7BR%7D%5En&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">.</p>
<p><span><strong>standardize</strong></span></p>
<p>When <code>standardize = TRUE</code> (default), columns of the data matrix <code>x</code> are standardized, i.e. each column of <code>x</code> has mean 0 and standard deviation 1. More specifically, we have that for each <img alt="j = 1, dots, p" title="j = 1, dots, p" src="https://s0.wp.com/latex.php?latex=j+%3D+1%2C+%5Cdots%2C+p&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">,</p>
<p><img alt="displaystylesum_{i=1}^n X_{ij} = 0" title="displaystylesum_{i=1}^n X_{ij} = 0" src="https://s0.wp.com/latex.php?latex=%5Cdisplaystyle%5Csum_%7Bi%3D1%7D%5En+X_%7Bij%7D+%3D+0&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">, and <img alt="sqrt{displaystylesum_{i=1}^n frac{X_{ij}^2}{n}} = 1" title="sqrt{displaystylesum_{i=1}^n frac{X_{ij}^2}{n}} = 1" src="https://s0.wp.com/latex.php?latex=%5Csqrt%7B%5Cdisplaystyle%5Csum_%7Bi%3D1%7D%5En+%5Cfrac%7BX_%7Bij%7D%5E2%7D%7Bn%7D%7D+%3D+1&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">.</p>
<p>Why might we want to do this? Standardizing our features before model fitting is common practice in statistical learning. This is because if our features are on vastly different scales, the features with larger scales will tend to dominate the action. (One instance where we might not want to standardize our features is if they are already all measured along the same scale, e.g. meters or kilograms.)</p>
<p>Notice that the standardization here is slightly different from that offered by the <code>scale</code> function: <code>scale(x, center = TRUE, scale = TRUE)</code> gives the standardization</p>
<p><img alt="displaystylesum_{i=1}^n X_{ij} = 0" title="displaystylesum_{i=1}^n X_{ij} = 0" src="https://s0.wp.com/latex.php?latex=%5Cdisplaystyle%5Csum_%7Bi%3D1%7D%5En+X_%7Bij%7D+%3D+0&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">, and <img alt="sqrt{displaystylesum_{i=1}^n frac{X_{ij}^2}{n-1}} = 1" title="sqrt{displaystylesum_{i=1}^n frac{X_{ij}^2}{n-1}} = 1" src="https://s0.wp.com/latex.php?latex=%5Csqrt%7B%5Cdisplaystyle%5Csum_%7Bi%3D1%7D%5En+%5Cfrac%7BX_%7Bij%7D%5E2%7D%7Bn-1%7D%7D+%3D+1&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">.</p>
<p>We verify this with a small data example. Generate data according to the following code:</p>
<pre>
n &lt;- 100; p &lt;- 5; true_p &lt;- 2
set.seed(950)
X &lt;- matrix(rnorm(n * p), nrow = n)
beta &lt;- matrix(c(rep(1, true_p), rep(0, p - true_p)), ncol = 1)
y &lt;- X %*% beta + 3 * rnorm(n)
</pre>
<p>Create a version of the data matrix which has standardized columns:</p>
<pre>
X_centered &lt;- apply(X, 2, function(x) x - mean(x))
Xs &lt;- apply(X_centered, 2, function(x) x / sqrt(sum(x^2) / n))
</pre>
<p>Next, we run <code>glmnet</code> on <code>Xs</code> and <code>y</code> with both possible options for <code>standardize</code>:</p>
<pre>
library(glmnet)
fit &lt;- glmnet(Xs, y, standardize = TRUE)
fit2 &lt;- glmnet(Xs, y, standardize = FALSE)
</pre>
<p>We can check that we get the same fit in both cases (modulo numerical precision):</p>
<pre>
sum(fit$lambda != fit2$lambda)
# 0
max(abs(fit$beta - fit2$beta))
# 6.661338e-16
</pre>
<p>The documentation notes that the coefficients returned are on the original scale. Let’s confirm that with our small data set. Run <code>glmnet</code> with the original data matrix and <code>standardize = TRUE</code>:</p>
<pre>
fit3 &lt;- glmnet(X, y, standardize = TRUE)
</pre>
<p>For each column <img alt="j" title="j" src="https://s0.wp.com/latex.php?latex=j&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">, our standardized variables are <img alt="Z_j = dfrac{X_j - mu_j}{s_j}" title="Z_j = dfrac{X_j - mu_j}{s_j}" src="https://s0.wp.com/latex.php?latex=Z_j+%3D+%5Cdfrac%7BX_j+-+%5Cmu_j%7D%7Bs_j%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">, where <img alt="mu_j" title="mu_j" src="https://s0.wp.com/latex.php?latex=%5Cmu_j&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> and <img alt="s_j" title="s_j" src="https://s0.wp.com/latex.php?latex=s_j&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> are the mean and standard deviation of column <img alt="j" title="j" src="https://s0.wp.com/latex.php?latex=j&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> respectively. If <img alt="beta_j" title="beta_j" src="https://s0.wp.com/latex.php?latex=%5Cbeta_j&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> and <img alt="gamma_j" title="gamma_j" src="https://s0.wp.com/latex.php?latex=%5Cgamma_j&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> represent the model coefficients of <code>fit2</code> and <code>fit3</code> respectively, then we should have</p>
<p><img alt="begin{aligned} beta_0 + sum_{j=1}^p beta_j Z_j &amp;= gamma_0 + sum_{j=1}^p gamma_j X_j, \  beta_0 + sum_{j=1}^p beta_j frac{X_j - mu_j}{s_j} &amp;= gamma_0 + sum_{j=1}^p gamma_j X_j, \  left( beta_0 - sum_{j=1}^p frac{mu_j}{s_j} right) + sum_{j=1}^p frac{beta_j}{s_j} X_j &amp;= gamma_0 + sum_{j=1}^p gamma_j X_j,  end{aligned}" title="begin{aligned} beta_0 + sum_{j=1}^p beta_j Z_j &amp;= gamma_0 + sum_{j=1}^p gamma_j X_j, \  beta_0 + sum_{j=1}^p beta_j frac{X_j - mu_j}{s_j} &amp;= gamma_0 + sum_{j=1}^p gamma_j X_j, \  left( beta_0 - sum_{j=1}^p frac{mu_j}{s_j} right) + sum_{j=1}^p frac{beta_j}{s_j} X_j &amp;= gamma_0 + sum_{j=1}^p gamma_j X_j,  end{aligned}" src="https://s0.wp.com/latex.php?latex=%5Cbegin%7Baligned%7D+%5Cbeta_0+%2B+%5Csum_%7Bj%3D1%7D%5Ep+%5Cbeta_j+Z_j+%26%3D+%5Cgamma_0+%2B+%5Csum_%7Bj%3D1%7D%5Ep+%5Cgamma_j+X_j%2C+%5C%5C++%5Cbeta_0+%2B+%5Csum_%7Bj%3D1%7D%5Ep+%5Cbeta_j+%5Cfrac%7BX_j+-+%5Cmu_j%7D%7Bs_j%7D+%26%3D+%5Cgamma_0+%2B+%5Csum_%7Bj%3D1%7D%5Ep+%5Cgamma_j+X_j%2C+%5C%5C++%5Cleft%28+%5Cbeta_0+-+%5Csum_%7Bj%3D1%7D%5Ep+%5Cfrac%7B%5Cmu_j%7D%7Bs_j%7D+%5Cright%29+%2B+%5Csum_%7Bj%3D1%7D%5Ep+%5Cfrac%7B%5Cbeta_j%7D%7Bs_j%7D+X_j+%26%3D+%5Cgamma_0+%2B+%5Csum_%7Bj%3D1%7D%5Ep+%5Cgamma_j+X_j%2C++%5Cend%7Baligned%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"></p>
<p>i.e. we should have <img alt="gamma_0 = beta_0 - sum_{j=1}^p frac{mu_j}{s_j}" title="gamma_0 = beta_0 - sum_{j=1}^p frac{mu_j}{s_j}" src="https://s0.wp.com/latex.php?latex=%5Cgamma_0+%3D+%5Cbeta_0+-+%5Csum_%7Bj%3D1%7D%5Ep+%5Cfrac%7B%5Cmu_j%7D%7Bs_j%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> and <img alt="gamma_j = frac{beta_j}{s_j}" title="gamma_j = frac{beta_j}{s_j}" src="https://s0.wp.com/latex.php?latex=%5Cgamma_j+%3D+%5Cfrac%7B%5Cbeta_j%7D%7Bs_j%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> for <img alt="j = 1, dots, p" title="j = 1, dots, p" src="https://s0.wp.com/latex.php?latex=j+%3D+1%2C+%5Cdots%2C+p&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">. The code below checks that this is indeed the case (modulo numerical precision):</p>
<pre>
# get column means and SDs
X_mean &lt;- colMeans(X)
X_sd &lt;- apply(X_centered, 2, function(x) sqrt(sum(x^2) / n))

# check difference for intercepts
fit2_int &lt;- coefficients(fit2)[1,]
fit3_int &lt;- coefficients(fit3)[1,]
temp &lt;- fit2_int - colSums(diag(X_mean / X_sd) %*% fit2$beta)
max(abs(temp - fit3_int))
# 1.110223e-16

# check difference for feature coefficients
temp &lt;- diag(1 / X_sd) %*% fit2$beta
max(abs(temp - fit3$beta))
# 1.110223e-15
</pre>
<p>The discussion above has been for the standardization of <code>x</code>. What about standardization for <code>y</code>? The documentation notes that when <code>family = &quot;gaussian&quot;</code>, <code>y</code> is automatically standardized, and the coefficients are unstandardized at the end of the procedure.</p>
<p>More concretely, let the mean and standard deviation of <img alt="y" title="y" src="https://s0.wp.com/latex.php?latex=y&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> be denoted by <img alt="mu_y" title="mu_y" src="https://s0.wp.com/latex.php?latex=%5Cmu_y&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> and <img alt="s_y" title="s_y" src="https://s0.wp.com/latex.php?latex=s_y&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> respectively. If running <code>glmnet</code> on standardized <code>y</code> gives intercept <img alt="beta_0" title="beta_0" src="https://s0.wp.com/latex.php?latex=%5Cbeta_0&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> and coefficients <img alt="beta_1, dots, beta_p" title="beta_1, dots, beta_p" src="https://s0.wp.com/latex.php?latex=%5Cbeta_1%2C+%5Cdots%2C+%5Cbeta_p&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">, then <code>glmnet</code> on unstandardized <code>y</code> will give intercept <img alt="mu_y + s_ybeta_0" title="mu_y + s_ybeta_0" src="https://s0.wp.com/latex.php?latex=%5Cmu_y+%2B+s_y%5Cbeta_0&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> and coefficients <img alt="s_ybeta_1, dots, s_ybeta_p" title="s_ybeta_1, dots, s_ybeta_p" src="https://s0.wp.com/latex.php?latex=s_y%5Cbeta_1%2C+%5Cdots%2C+s_y%5Cbeta_p&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">.</p>
<p>Again, this can be verified empirically:</p>
<pre>
# get mean and SD of y
y_mean &lt;- mean(y)
y_sd &lt;- sqrt(sum((y - y_mean)^2) / n)

# fit model with standardized y
fit4 &lt;- glmnet(X, (y - y_mean) / y_sd, standardize = TRUE)

# check difference for intercepts
fit4_int &lt;- coefficients(fit4)[1,]
temp &lt;- fit4_int * y_sd + y_mean
max(abs(temp - fit3_int))
# 1.110223e-16

# check difference for feature coefficients
max(abs(y_sd * fit4$beta - fit3$beta))
# 8.881784e-16
</pre>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://statisticaloddsandends.wordpress.com/2018/11/15/a-deep-dive-into-glmnet-standardize/"> R – Statistical Odds &amp; Ends</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/fyV2bAjfrEo" height="1"></a></div>