<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/choosing-hyper-parameters-in-penalized-regression/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://privefl.github.io/blog/choosing-hyper-parameters-in-penalized-regression/"> Florian Privé</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<section>
<p>In this post, I’m evaluating some ways of choosing hyper-parameters (<span>\(\alpha\)</span> and <span>\(\lambda\)</span>) in penalized linear regression. The same principles can be applied to other types of penalized regresions (e.g. logistic).</p>
<div>
<h2>Model</h2>
<p>In penalized linear regression, we find regression coefficients <span>\(\hat{\beta}_0\)</span> and <span>\(\hat{\beta}\)</span> that minimize the following regularized loss function <span>\[L(\lambda, \alpha) = \underbrace{ \frac{1}{2n} \sum_{i=1}^n \left( y_i – \hat{y}_i \right)^2 }_\text{Loss function} + \underbrace{ \lambda \left((1-\alpha)\frac{1}{2}\|\hat{\beta}\|_2^2 + \alpha \|\hat{\beta}\|_1\right) }_\text{Penalization} ~,\]</span> where <span>\(\hat{y}_i=\hat{\beta}_0 + x_i^T\hat{\beta}\)</span>, <span>\(0 \le \alpha \le 1\)</span> and <span>\(\lambda > 0\)</span>.</p>
<p>This regularization is called elastic-net and has two particular cases, namely LASSO (<span>\(\alpha = 1\)</span>) and ridge (<span>\(\alpha = 0\)</span>). So, in elastic-net regularization, hyper-parameter <span>\(\alpha\)</span> accounts for the relative importance of the L1 (LASSO) and L2 (ridge) regularizations. There is another hyper-parameter, <span>\(\lambda\)</span>, that accounts for the amount of regularization used in the model.</p>
<p>I won’t discuss the benefits of using regularization here.</p>
</div>
<div>
<h2>Data</h2>
<p>In this blog post, I’m using the <a rel="nofollow" href="http://myweb.uiowa.edu/pbreheny/data/bcTCGA.html" target="_blank">Breast cancer gene expression data from The Cancer Genome Atlas</a>.</p>
<div>
<pre><code><span>dir.create</span>(dir &lt;-<span> &quot;/tmp/data&quot;</span>, <span>showWarnings =</span> <span>FALSE</span>)
file &lt;-<span> </span><span>file.path</span>(dir, <span>&quot;bcTCGA.rds&quot;</span>)
<span>if</span> (<span>!</span><span>file.exists</span>(file)) {
  <span>download.file</span>(<span>&quot;https://s3.amazonaws.com/pbreheny-data-sets/bcTCGA.rds&quot;</span>,
                <span>destfile =</span> file)
}
bcTCGA &lt;-<span> </span><span>readRDS</span>(file)</code></pre>
</div>
<div>
<pre><code>X &lt;-<span> </span>bcTCGA<span>$</span>X
<span>dim</span>(X)</code></pre>
</div>
<pre><code>## [1]   536 17322</code></pre>
<div>
<pre><code>X[<span>1</span><span>:</span><span>5</span>, <span>1</span><span>:</span><span>5</span>]</code></pre>
</div>
<pre><code>##        15E1.2     2'-PDE     7A5       A1BG      A2BP1
## [1,] -1.44775  0.0153750  2.4105  0.9493333  0.5391667
## [2,] -2.29950  0.4669375  0.3635  0.2798333 -1.5346667
## [3,] -1.94400 -0.0295625  1.8550  0.7486667 -1.3048333
## [4,] -2.09800  0.7919375  1.4080  0.7500000 -1.3578333
## [5,] -1.28500 -0.1911250 -1.3005 -0.5471667 -2.1845000</code></pre>
<div>
<pre><code>y &lt;-<span> </span>bcTCGA<span>$</span>y
<span>hist</span>(y)</code></pre>
</div>
<p><img src="https://i1.wp.com/privefl.github.io/blog/knitr_files/post-glmnet_files/figure-html/unnamed-chunk-2-1.png?w=456&ssl=1" data-recalc-dims="1"></p>
</div>
<div>
<h2>Methods compared</h2>
<p>I’m comparing 6 different estimations of <span>\(y\)</span> with the following code. I’m explaining what are these estimations after the code.</p>
<div>
<pre><code>timing &lt;-<span> </span><span>function</span>(expr) <span>system.time</span>(expr)[<span>3</span>]

RMSE &lt;-<span> </span><span>function</span>(pred, y.test) {
  <span>stopifnot</span>(<span>length</span>(pred) <span>==</span><span> </span><span>length</span>(y.test))
  <span>mean</span>((pred <span>-</span><span> </span>y.test)<span>^</span><span>2</span>)
}

<span>set.seed</span>(<span>1</span>)
<span># !! 5h of computations on my computer !! </span>
res_all &lt;-<span> </span><span>replicate</span>(<span>200</span>, <span>simplify =</span> <span>FALSE</span>, {

  ind.train &lt;-<span> </span><span>sample</span>(<span>nrow</span>(X), <span>400</span>)
  ind.test  &lt;-<span> </span><span>setdiff</span>(<span>1</span><span>:</span><span>nrow</span>(X), ind.train)

  <span>library</span>(glmnet)
  t_all &lt;-<span> </span><span>timing</span>(mod_all &lt;-<span> </span><span>glmnet</span>(X[ind.train, ], y[ind.train]))
  preds_all &lt;-<span> </span><span>predict</span>(mod_all, X[ind.test, ])
  rmse_all &lt;-<span> </span><span>apply</span>(preds_all, <span>2</span>, RMSE, y[ind.test])

  t_cv &lt;-<span> </span><span>timing</span>(mod_cv &lt;-<span> </span><span>cv.glmnet</span>(X[ind.train, ], y[ind.train]))
  preds_1se &lt;-<span> </span><span>predict</span>(mod_cv, X[ind.test, ], <span>s =</span> <span>&quot;lambda.1se&quot;</span>)
  rmse_1se &lt;-<span> </span><span>RMSE</span>(preds_1se, y[ind.test])
  preds_min &lt;-<span> </span><span>predict</span>(mod_cv, X[ind.test, ], <span>s =</span> <span>&quot;lambda.min&quot;</span>)
  rmse_min &lt;-<span> </span><span>RMSE</span>(preds_min, y[ind.test])

  <span>library</span>(bigstatsr)
  t_CMSA &lt;-<span> </span><span>timing</span>({
    X2 &lt;-<span> </span><span>as_FBM</span>(X)
    mod_CMSA &lt;-<span> </span><span>big_spLinReg</span>(X2, y[ind.train], ind.train)
  })
  preds_CMSA &lt;-<span> </span><span>predict</span>(mod_CMSA, X2, ind.test)
  rmse_CMSA &lt;-<span> </span><span>RMSE</span>(preds_CMSA, y[ind.test])

  <span>library</span>(glmnetUtils)
  ALPHA &lt;-<span> </span><span>c</span>(<span>1</span>, <span>0.5</span>, <span>0.1</span>)
  t_cva &lt;-<span> </span><span>timing</span>(mod_cva &lt;-<span> </span><span>cva.glmnet</span>(X[ind.train, ], y[ind.train], <span>alpha =</span> ALPHA))
  alpha &lt;-<span> </span>ALPHA[<span>which.min</span>(<span>sapply</span>(mod_cva<span>$</span>modlist, <span>function</span>(mod) <span>min</span>(mod<span>$</span>cvm)))]
  rmse_cva &lt;-<span> </span><span>RMSE</span>(<span>predict</span>(mod_cva, X[ind.test, ], <span>alpha =</span> alpha), y[ind.test])

  t_CMSA2 &lt;-<span> </span><span>timing</span>({
    X2 &lt;-<span> </span><span>as_FBM</span>(X)
    mod_CMSA2 &lt;-<span> </span><span>big_spLinReg</span>(X2, y[ind.train], ind.train, <span>alphas =</span> ALPHA)
  })
  preds_CMSA2 &lt;-<span> </span><span>predict</span>(mod_CMSA2, X2, ind.test)
  rmse_CMSA2 &lt;-<span> </span><span>RMSE</span>(preds_CMSA2, y[ind.test])

  tibble<span>::</span><span>tribble</span>(
    <span>~</span>method,        <span>~</span>timing,  <span>~</span>rmse,
    <span>&quot;glmnet_best&quot;</span>,  t_all,    <span>min</span>(rmse_all),
    <span>&quot;glmnet_min&quot;</span>,   t_cv,     rmse_1se,
    <span>&quot;glmnet_1se&quot;</span>,   t_cv,     rmse_min,
    <span>&quot;CMSA&quot;</span>,         t_CMSA,   rmse_CMSA,
    <span>&quot;glmnet_cva&quot;</span>,   t_cva,    rmse_cva,
    <span>&quot;CMSA2&quot;</span>,        t_CMSA2,  rmse_CMSA2
  )

})

res &lt;-<span> </span><span>do.call</span>(rbind, res_all)
res<span>$</span>run_number &lt;-<span> </span><span>rep</span>(<span>seq_along</span>(res_all), <span>each =</span> <span>6</span>)</code></pre>
</div>
<p>The methods compared are:</p>
<ul>
<li><strong>glmnet_best</strong>: the best prediction (on the test set) among all 100 different <span>\(\lambda\)</span> values (using <span>\(\alpha = 1\)</span>). This is an upper bound of predictive performance using package {glmnet} with L1 regularization.</li>
<li><strong>glmnet_min</strong>: the prediction corresponding to the <span>\(\lambda\)</span> that minimizes the error of cross-validation.</li>
<li><strong>glmnet_1se</strong>: similar to <strong>glmnet_min</strong>, but a little more regularized (the maximum regularization whose error of cross-validation is only one standard error away from the minimum).</li>
<li><strong>CMSA</strong>: using cross-model selection and averaging (CMSA, explained in the next session) to “choose” <span>\(\lambda\)</span>.</li>
<li><strong>glmnet_cva</strong>: using cross-validation of glmnet on both <span>\(\lambda\)</span> AND <span>\(\alpha\)</span>. Using <code>&quot;lambda.1se&quot;</code> (the default) and <span>\(\alpha\)</span> that minimizes the error of cross-validation.</li>
<li><strong>CMSA2</strong>: using cross-model selection and averaging (CMSA, explained in the next session) to “choose” <span>\(\lambda\)</span> AND <span>\(\alpha\)</span>.</li>
</ul>
</div>
<div>
<h2>Cross-Model Selection and Averaging (CMSA)</h2>
<p><strong>Source:</strong> <a rel="nofollow" href="https://doi.org/10.1101/403337" target="_blank">https://doi.org/10.1101/403337</a></p>
<p><img src="https://i1.wp.com/privefl.github.io/blog/images/simple-CMSA.png?w=80%25&ssl=1" data-recalc-dims="1"></p>
<p>Illustration of one turn of the Cross-Model Selection and Averaging (CMSA) procedure. First, this procedure separates the training set in <span>\(K\)</span> folds (e.g. 10 folds). Secondly, in turn, each fold is considered as an inner validation set (red) and the other (<span>\(K – 1\)</span>) folds form an inner training set (blue). A “regularization path” of models is trained on the inner training set and the corresponding predictions (scores) for the inner validation set are computed. The model that minimizes the loss on the inner validation set is selected. Finally, the <span>\(K\)</span> resulting models are averaged. We also use this procedure to derive an early stopping criterion so that the algorithm does not need to evaluate the whole regularization paths, making this procedure much faster.</p>
</div>
<div>
<h2>Results</h2>
<div>
<pre><code><span>library</span>(dplyr, <span>warn.conflicts =</span> <span>FALSE</span>)

res <span>%>%</span>
<span> </span><span>group_by</span>(method) <span>%>%</span>
<span> </span><span>summarise_at</span>(<span>c</span>(<span>&quot;rmse&quot;</span>, <span>&quot;timing&quot;</span>), mean) <span>%>%</span>
<span> </span><span>arrange</span>(rmse)</code></pre>
</div>
<pre><code>## # A tibble: 6 x 3
##   method       rmse timing
##   </code></pre></div></section>