<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/elastic-net-regularization-of-a-model-of-burned-calories/"></div></div><div><a data-url="https://www.r-bloggers.com/elastic-net-regularization-of-a-model-of-burned-calories/" data-counturl="https://www.r-bloggers.com/elastic-net-regularization-of-a-model-of-burned-calories/" data-count="horizontal" data-text="Elastic net regularization of a model of burned calories" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://ellisp.github.io/blog/2016/08/13/fitbit-lasso"> Peter's stats stuff - R</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h2>Deal with feature selection and collinearity</h2>
<p>Recently I’ve been making more use of <a target="_blank" href="https://en.wikipedia.org/wiki/Elastic_net_regularization">elastic net regularization</a> as a way of fitting linear models to data when I have more candidate explanatory variables than I know what to do with and some of them are collinear ie their information doubles up on what is in other variables.</p>
<p>Elastic net regularization is a more general form of the ridge regression and lasso methods. Ridge regression was primarily a method to deal with the instability of estimates of the coefficients of linear models when they are collinear; the lasso a method of feature selection that forces coefficients for some/many explanatory variables to be zero and provides good estimates of the coefficients of the remaining features. Both are preferable to stepwise selection which provides estimates of effect sizes biased upwards in absolute magnitude. I deliberately use the term “was” for ridge regression because I don’t see any reason for using either of these methods any more; they are superseded by the more general elastic net. A hyperparameter <code>alpha</code> means the method is equivalent to lasso when <code>alpha == 1</code> and to ridge regression when <code>alpha == 0</code>. We can use cross-validation assessing the fit of the best models with different values of <code>alpha</code> to work out the best value of <code>alpha</code>.</p>
<p>The method works as an extension of ordinary least squares (OLS), the most common and basic (ie Stats 101) method of estimating the coefficients in a linear model. Instead of minimising the sum of the squares of the differences between predicted and actual values of the response variables (the OLS method), the elastic net regularization estimation process minimises that sum of those squares <em>and</em> a penalty based on the size of the estimated coefficients (on a standardised scale) in the model. I like to think of the impact of this as to create an increased burden of proof on the data for including a variable’s effect in the model; the estimation process drags coefficients towards zero in a spirit of conservatism. The <code>alpha</code> hyperparameter controls the <em>balance</em> between the different penalty methods of ridge and lasso. A second hyperparameter, <code>lambda</code> sets the total <em>weight</em> of the penalty, and can also be set at the most effective level by trying different values with cross validation.</p>
<p>While ridge regression has been around for decades and the <a target="_blank" href="http://www.jstor.org/stable/2346178">lasso since the 1990s</a>, the generalisation to elastic net regularations was only <a target="_blank" href="https://web.stanford.edu/~hastie/Papers/B67.2%20(2005)%20301-320%20Zou%20&amp;%20Hastie.pdf">properly elaborated in 2005</a> and software to implement it is fairly new. The <a target="_blank" href="https://cran.r-project.org/web/packages/glmnet/index.html"><code>glmnet</code> package in R</a> gives very fast estimation and good cross-validation methods, and scales up well to large and wide datasets.</p>
<h2>A modest sized dataset with collinearity</h2>
<p>To build more familiarity with this tool I grabbed a small dataset of 30 observations from my fitbit tracker. I originally downloaded the data from the web browser interface in accordance with the <a target="_blank" href="https://help.fitbit.com/articles/en_US/Help_article/1133">fitbit instructions</a>. For reproducibility I’ve made it <a target="_blank" href="https://raw.githubusercontent.com/ellisp/ellisp.github.io/source/data/fitbit_export_20160806.csv">available to others</a> on the web.</p>
<p>My aim with this data is to reverse engineer the way fitbit estimates calories burned in a day. I know from <a target="_blank" href="https://help.fitbit.com/articles/en_US/Help_article/1381">their website</a> that this is done based on an estimate of my basal metabolic rate and the activities my tracker records (plus any activities manually logged by me, which is none in my case).</p>
<p>Basal metabolic rate is dependent on my gender, age, height and weight, none of which changed in this thirty day period, so I would expect that component of the calories burnt to be constant over time. As well as the calories burned per day, the information the fitbit provides on my activities is available in these variables:</p>
<ul>
<li>steps</li>
<li>distance (estimated from steps)</li>
<li>floors (ie of stairs and stair-like equivalents – some of this period I walked up hills in Western Australia’s Stirling Ranges)</li>
<li>sedentary minutes</li>
<li>lightly active minutes</li>
<li>fairly active minutes</li>
<li>very active minutes</li>
</ul>
<p>Here’s how those seven candidate explanatory variables relate to eachother:</p>
<p><img alt="pairs" src="http://ellisp.github.io/img/0050-pairs.svg"></p>
<p>From this chart we see a few things of interest:</p>
<ul>
<li><code>Distance</code> can be calculated almost exactly directly from <code>Steps</code></li>
<li>The number of sedentary minutes in a day is inversely related to the other minutes measured – lightly active, fairly active, and very active</li>
</ul>
<p>I wondered if the four minutes measures would be perfectly collinear – there’s a set number of minutes each day after all, so if I know how long I’ve been lightly active, fairly active and very active couldn’t I precisely say that the remaining minutes were sedentary? A bit of exploration shows that the four “minutes” measures tend to add up to 15 or 16 hours a day, suggesting that sedentary minutes actually means “sedentary but awake minutes”. So they are not completely collinear because the amount I sleep varies from day to day.</p>
<p>For those who want to follow along in R, here’s how I did that chart with the correlation coefficients in the lower triangle of panels. This adapts an example from the helpfile for the <code>pairs()</code> function in base graphics.</p>
<figure>
<pre><code data-lang="r"><span>library</span><span>(</span><span>dplyr</span><span>)</span><span>
</span><span>library</span><span>(</span><span>tidyr</span><span>)</span><span>
</span><span>library</span><span>(</span><span>ggplot2</span><span>)</span><span>
</span><span>library</span><span>(</span><span>caret</span><span>)</span><span>
</span><span>library</span><span>(</span><span>glmnet</span><span>)</span><span>
</span><span>library</span><span>(</span><span>boot</span><span>)</span><span>
</span><span>library</span><span>(</span><span>RColorBrewer</span><span>)</span><span>
</span><span>fitbit</span><span> </span><span>&lt;-</span><span> </span><span>read.csv</span><span>(</span><span>&quot;https://raw.githubusercontent.com/ellisp/ellisp.github.io/source/data/fitbit_export_20160806.csv&quot;</span><span>,</span><span> </span><span>skip</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span> </span><span>stringsAsFactors</span><span> </span><span>=</span><span> </span><span>FALSE</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>mutate</span><span>(</span><span>
</span><span>Calories.Burned</span><span> </span><span>=</span><span> </span><span>as.numeric</span><span>(</span><span>gsub</span><span>(</span><span>&quot;,&quot;</span><span>,</span><span> </span><span>&quot;&quot;</span><span>,</span><span> </span><span>Calories.Burned</span><span>)),</span><span>
</span><span>Steps</span><span> </span><span>=</span><span> </span><span>as.numeric</span><span>(</span><span>gsub</span><span>(</span><span>&quot;,&quot;</span><span>,</span><span> </span><span>&quot;&quot;</span><span>,</span><span> </span><span>Steps</span><span>)),</span><span>
</span><span>Activity.Calories</span><span> </span><span>=</span><span> </span><span>as.numeric</span><span>(</span><span>gsub</span><span>(</span><span>&quot;,&quot;</span><span>,</span><span> </span><span>&quot;&quot;</span><span>,</span><span> </span><span>Activity.Calories</span><span>)),</span><span>
</span><span>Date</span><span> </span><span>=</span><span> </span><span>as.Date</span><span>(</span><span>Date</span><span>,</span><span> </span><span>format</span><span> </span><span>=</span><span> </span><span>&quot;%d/%m/%Y&quot;</span><span>)</span><span>
</span><span>)</span><span>
</span><span>fitbit2</span><span> </span><span>&lt;-</span><span> </span><span>fitbit</span><span>
</span><span>fitbit2</span><span>$</span><span>Activity.Calories</span><span> </span><span>&lt;-</span><span> </span><span>NULL</span><span>
</span><span>fitbit2</span><span>$</span><span>Date</span><span> </span><span>&lt;-</span><span> </span><span>NULL</span><span>
</span><span>fitbit2</span><span>$</span><span>Steps</span><span> </span><span>&lt;-</span><span> </span><span>fitbit2</span><span>$</span><span>Steps</span><span> </span><span>/</span><span> </span><span>1000</span><span> </span><span># so coefficients will be calories per thousand steps, easier to describe
</span><span>
</span><span>par</span><span>(</span><span>family</span><span> </span><span>=</span><span> </span><span>&quot;myfont&quot;</span><span>)</span><span>
</span><span># let's look at all the candidate variables
</span><span>panel.cor</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>x</span><span>,</span><span> </span><span>y</span><span>,</span><span> </span><span>digits</span><span> </span><span>=</span><span> </span><span>2</span><span>,</span><span> </span><span>prefix</span><span> </span><span>=</span><span> </span><span>&quot;&quot;</span><span>,</span><span> </span><span>cex.cor</span><span>,</span><span> </span><span>...</span><span>){</span><span>
</span><span># function for calculating correlation coefficients and using
</span><span> </span><span># as panels in pairs() function. Taken from ?pairs
</span><span> </span><span>usr</span><span> </span><span>&lt;-</span><span> </span><span>par</span><span>(</span><span>&quot;usr&quot;</span><span>);</span><span> </span><span>on.exit</span><span>(</span><span>par</span><span>(</span><span>usr</span><span>))</span><span>
</span><span>par</span><span>(</span><span>usr</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>0</span><span>,</span><span> </span><span>1</span><span>,</span><span> </span><span>0</span><span>,</span><span> </span><span>1</span><span>))</span><span>
</span><span>r</span><span> </span><span>&lt;-</span><span> </span><span>abs</span><span>(</span><span>cor</span><span>(</span><span>x</span><span>,</span><span> </span><span>y</span><span>))</span><span>
</span><span>txt</span><span> </span><span>&lt;-</span><span> </span><span>format</span><span>(</span><span>c</span><span>(</span><span>r</span><span>,</span><span> </span><span>0.123456789</span><span>),</span><span> </span><span>digits</span><span> </span><span>=</span><span> </span><span>digits</span><span>)[</span><span>1</span><span>]</span><span>
</span><span>txt</span><span> </span><span>&lt;-</span><span> </span><span>paste0</span><span>(</span><span>prefix</span><span>,</span><span> </span><span>txt</span><span>)</span><span>
</span><span>if</span><span>(</span><span>missing</span><span>(</span><span>cex.cor</span><span>))</span><span> </span><span>cex.cor</span><span> </span><span>&lt;-</span><span> </span><span>0.8</span><span>/</span><span>strwidth</span><span>(</span><span>txt</span><span>)</span><span>
</span><span>text</span><span>(</span><span>0.5</span><span>,</span><span> </span><span>0.5</span><span>,</span><span> </span><span>txt</span><span>,</span><span> </span><span>cex</span><span> </span><span>=</span><span> </span><span>cex.cor</span><span> </span><span>*</span><span> </span><span>r</span><span>)</span><span>
</span><span>}</span><span>
</span><span>pairs</span><span>(</span><span>fitbit2</span><span>[</span><span> </span><span>,</span><span> </span><span>-1</span><span>],</span><span> </span><span>lower.panel</span><span> </span><span>=</span><span> </span><span>panel.cor</span><span>,</span><span> </span><span>main</span><span> </span><span>=</span><span> </span><span>&quot;Pairwise relationships of Fitbit's measured activities&quot;</span><span>)</span></code></pre>
</figure>
<p>The close relationship between steps and distance (correlation coefficient of 1.00, which means that knowing one you can almost exactly predict the other) isn’t surprising to me because fitbit’s website makes clear it estimates distance based on steps and stride length. It applies one stride length to convert from steps to distance for walking, and one for running (which it can detect, presumably, from the pace and violence of motion), but I know that I did virtually no running in this period so it’s only converting from steps to walking. Here’s the distribution of fitbit’s estimated distance divided by steps:</p>
<p><img alt="stride" src="http://ellisp.github.io/img/0050-stridelenth.svg"></p>
<figure>
<pre><code data-lang="r"><span>ggplot</span><span>(</span><span>fitbit2</span><span>,</span><span> </span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>Distance</span><span> </span><span>/</span><span> </span><span>Steps</span><span>))</span><span> </span><span>+</span><span> </span><span>geom_rug</span><span>()</span><span> </span><span>+</span><span> </span><span>geom_density</span><span>()</span><span> </span><span>+</span><span>
</span><span>ggtitle</span><span>(</span><span>&quot;Stride length reverse-engineered from Fitbit data&quot;</span><span>,</span><span>
</span><span>subtitle</span><span> </span><span>=</span><span> </span><span>&quot;Not all strides identical, due to rounding or other jitter&quot;</span><span>)</span></code></pre>
</figure>
<h2>69 calories per 1,000 steps on average</h2>
<p>To get a sense of the data I was interested in a super-simple model with a single explanatory variable. “What’s the relationship between the number of steps walked and calories burned, ignoring height gains and intensity of heart rate”, I asked myself? Judging from various sites like <a target="_blank" href="http://www.livestrong.com/article/320124-how-many-calories-does-the-average-person-use-per-step/">this</a> and <a target="_blank" href="http://www.livestrong.com/article/238020-how-to-convert-pedometer-steps-to-calories/">this</a> the answer might be about 40-50 calories per thousand steps for walking. Fitbit would use a more sophisticated model than this and take into account the fact that some of my walking was on tougher terrain and/or up stairs or slopes which would reflect in increased heart rate, so it will be interesting to see its estimates of calories per thousand steps compared to the simple website estimates.</p>
<p>The way to do this is just an ordinary least squares estimation of a model with steps as the sole explanatory variable and calories burned as the response variable. This returns the estimate of my base daily rate of 1,926 calories per day, and 69 additional calories burned per thousand steps. This figure of 69 is more than the “moderate walking” rates on various websites, presumably because of the terrain and slope issues mentioned above.</p>
<figure>
<pre><code data-lang="r"><span>mod0</span><span> </span><span>&lt;-</span><span> </span><span>lm</span><span>(</span><span>Calories.Burned</span><span> </span><span>~</span><span> </span><span>Steps</span><span>,</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>fitbit2</span><span>)</span><span>
</span><span>round</span><span>(</span><span>coef</span><span>(</span><span>mod0</span><span>))</span><span>
</span><span>## (Intercept) Steps </span><span>##</span><span> </span><span>1926</span><span> </span><span>69</span><span> </span></code></pre>
</figure>
<p>That figure of 1,926 per day as my “no steps” base calory burn is also somewhat more than the 1,629 that <a target="_blank" href="http://www.bodybuilding.com/fun/bmr_calculator.htm">this website</a> says should be my basal metabolic rate based on the Harris-Benedict equation. I’m not sure what’s going on here – either fitbit uses a different calculation, or perhaps it’s judging various activities from me from my heart rate that don’t feature in steps, or something else. Fitbit doesn’t know my real calory burn rates, so it can’t be something mysterious to do with my metabolism as it has no way of knowing, it just estimates based on steps, heart rate, weight, height, etc. Anyway, the simple model above isn’t too bad an estimate of how many calories fitbit <em>says</em> I use; for the record, here’s a comparison of the predicted calories compared to residuals from that model:</p>
<p><img alt="residuals" src="http://ellisp.github.io/img/0050-1var-resid.svg"></p>
<p>And just to put one issue to rest, here’s the partial autocorrelation function of the residuals, showing that there’s no particular pattern over time not captured in the model. This is not a particularly interesting chart but I’m including it because it should be standard practice to check for such patterns (which would be shown by one of the vertical bars crossing the horizontal blue lines) when working with data that was captured over a time period.</p>
<p><img alt="pacf" src="http://ellisp.github.io/img/0050-1var-resid-pacf.svg"></p>
<figure>
<pre><code data-lang="r"><span># check residuals v. fit
</span><span>plot</span><span>(</span><span>mod0</span><span>,</span><span> </span><span>which</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span> </span><span>bty</span><span> </span><span>=</span><span> </span><span>&quot;l&quot;</span><span>);</span><span> </span><span>grid</span><span>()</span><span>
</span><span># check autocorrelation:
</span><span>pacf</span><span>(</span><span>resid</span><span>(</span><span>mod0</span><span>),</span><span> </span><span>main</span><span> </span><span>=</span><span> </span><span>&quot;Partial autocorrelation of residuals from single variable regression&quot;</span><span>)</span></code></pre>
</figure>
<h2>Choose alpha by cross-validation over a grid of candidate values</h2>
<p>Now I’m ready to move to my full model, predicting daily calories based on all the seven explanatory variables available.</p>
<p>The first step is choosing an appropriate alpha – the balance between the extremes of ridge regression and lasso estimation. I choose repeated cross-validation to do this. Basically this means fitting the model to slightly different re-samples many times at different values of alpha, using the fit model to predict the “out of bag” points from the original sample that weren’t in the re-sample. The results show only a weak impact of the choice of alpha:</p>
<p><img alt="alphas" src="http://ellisp.github.io/img/0050-cv-alpha.svg"></p>
<p>Here’s the code that did this.</p>
<figure>
<pre><code data-lang="r"><span># create X matrix and Y vector for use with glmnet, which doesn't take Rmodel formulae
# First column of fitbit2 is the response variable &quot;Calories burned&quot; so we exclude from X:
</span><span>X</span><span> </span><span>&lt;-</span><span> </span><span>as.matrix</span><span>(</span><span>fitbit2</span><span>[</span><span> </span><span>,</span><span> </span><span>-1</span><span>])</span><span> </span><span># standardisation happens as part of glmnet
</span><span>Y</span><span> </span><span>&lt;-</span><span> </span><span>fitbit2</span><span>$</span><span>Calories.Burned</span><span>
</span><span>set.seed</span><span>(</span><span>123</span><span>)</span><span>
</span><span>alphas</span><span> </span><span>&lt;-</span><span> </span><span>seq</span><span>(</span><span>from</span><span> </span><span>=</span><span> </span><span>0</span><span>,</span><span> </span><span>to</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span> </span><span>length.out</span><span> </span><span>=</span><span> </span><span>10</span><span>)</span><span>
</span><span>res</span><span> </span><span>&lt;-</span><span> </span><span>matrix</span><span>(</span><span>0</span><span>,</span><span> </span><span>nrow</span><span> </span><span>=</span><span> </span><span>length</span><span>(</span><span>alphas</span><span>),</span><span> </span><span>ncol</span><span> </span><span>=</span><span> </span><span>6</span><span>)</span><span> </span><span># five columns for results - five repeats of each CV run
</span><span>for</span><span>(</span><span>i</span><span> </span><span>in</span><span> </span><span>1</span><span>:</span><span>length</span><span>(</span><span>alphas</span><span>)){</span><span>
</span><span>for</span><span>(</span><span>j</span><span> </span><span>in</span><span> </span><span>2</span><span>:</span><span>6</span><span>){</span><span>
</span><span>cvmod</span><span> </span><span>&lt;-</span><span> </span><span>cv.glmnet</span><span>(</span><span>X</span><span>,</span><span> </span><span>Y</span><span>,</span><span> </span><span>alpha</span><span> </span><span>=</span><span> </span><span>alphas</span><span>[</span><span>i</span><span>])</span><span>
</span><span>res</span><span>[</span><span>i</span><span>,</span><span> </span><span>c</span><span>(</span><span>1</span><span>,</span><span> </span><span>j</span><span>)]</span><span> </span><span>&lt;-</span><span> </span><span>c</span><span>(</span><span>alphas</span><span>[</span><span>i</span><span>],</span><span> </span><span>sqrt</span><span>(</span><span>min</span><span>(</span><span>cvmod</span><span>$</span><span>cvm</span><span>)))</span><span>
</span><span>}</span><span>
</span><span>}</span><span>
</span><span>res</span><span> </span><span>&lt;-</span><span> </span><span>data.frame</span><span>(</span><span>res</span><span>)</span><span>
</span><span>res</span><span>$</span><span>average_rmse</span><span> </span><span>&lt;-</span><span> </span><span>apply</span><span>(</span><span>res</span><span>[</span><span> </span><span>,</span><span> </span><span>2</span><span>:</span><span>6</span><span>],</span><span> </span><span>1</span><span>,</span><span> </span><span>mean</span><span>)</span><span>
</span><span>res</span><span> </span><span>&lt;-</span><span> </span><span>res</span><span>[</span><span>order</span><span>(</span><span>res</span><span>$</span><span>average_rmse</span><span>),</span><span> </span><span>]</span><span>
</span><span>names</span><span>(</span><span>res</span><span>)[</span><span>1</span><span>]</span><span> </span><span>&lt;-</span><span> </span><span>&quot;alpha&quot;</span><span>
</span><span>res</span><span> </span><span>%>%</span><span>
</span><span>select</span><span>(</span><span>-</span><span>average_rmse</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>gather</span><span>(</span><span>trial</span><span>,</span><span> </span><span>rmse</span><span>,</span><span> </span><span>-</span><span>alpha</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>ggplot</span><span>(</span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>alpha</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>rmse</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_point</span><span>()</span><span> </span><span>+</span><span>
</span><span>geom_smooth</span><span>(</span><span>se</span><span> </span><span>=</span><span> </span><span>FALSE</span><span>)</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;Root mean square error&quot;</span><span>)</span><span> </span><span>+</span><span>
</span><span>ggtitle</span><span>(</span><span>&quot;Cross validation best RMSE for differing values of alpha&quot;</span><span>)</span><span>
</span><span># best alpha varies according to the random seed set earlier but with seed 123 it is 0.667
</span><span>bestalpha</span><span> </span><span>&lt;-</span><span> </span><span>res</span><span>[</span><span>1</span><span>,</span><span> </span><span>1</span><span>]</span></code></pre>
</figure>
<h2>Elastic net gives more sensible results in presence of collinearity</h2>
<p>Now that I have a preferred value of alpha I can use the elastic net to get estimated values of the eight coefficients (seven explanatory variables plus an intercept), and compare them to the ordinary least squares equivalents. Remember that there’s a second hyper-parameter, <code>lambda</code> to estimate too – again I choose this from cross-validation with the chosen value of <code>alpha</code>, and following the implied suggesting on the <code>glmnet</code> package I try two versions of <code>lambda</code> – the level of shrinkage that gives the smallest errors in cross validation, and the higher level of shrinkage that gets acceptable close (within one standard deviation in the cross validation process) to that smallest error. This practice is predicated on the idea that a model with a higher degree of shrinkage may be cheaper (in data management or interpretation) than a lower, so it is worth taking a small hit of accuracy for it.</p>
<p>Here are the results from the three estimation methods:</p>
<figure>
<pre><code data-lang="r"><span> </span><span>original</span><span> </span><span>shrunk</span><span> </span><span>very.shrunk</span><span>
</span><span>(</span><span>Intercept</span><span>)</span><span> </span><span>1941.189</span><span> </span><span>1971.199</span><span> </span><span>2157.172</span><span>
</span><span>Steps</span><span> </span><span>-66.827</span><span> </span><span>9.176</span><span> </span><span>15.703</span><span>
</span><span>Distance</span><span> </span><span>116.077</span><span> </span><span>15.045</span><span> </span><span>21.835</span><span>
</span><span>Floors</span><span> </span><span>0.027</span><span> </span><span>0.000</span><span> </span><span>0.000</span><span>
</span><span>Minutes.Sedentary</span><span> </span><span>-0.246</span><span> </span><span>-0.236</span><span> </span><span>-0.187</span><span>
</span><span>Minutes.Lightly.Active</span><span> </span><span>2.247</span><span> </span><span>1.985</span><span> </span><span>0.888</span><span>
</span><span>Minutes.Fairly.Active</span><span> </span><span>4.402</span><span> </span><span>4.384</span><span> </span><span>3.842</span><span>
</span><span>Minutes.Very.Active</span><span> </span><span>3.895</span><span> </span><span>3.295</span><span> </span><span>1.019</span></code></pre>
</figure>
<p>As experienced modellers would have expected, the “original” (ie ordinary least squares) method gives non-sensible estimates of the coefficients of the two highly collinear variabls, <code>Steps</code> and <code>Distance</code>. it clearly makes no sense to say that my calories consumed go <em>down</em> by 67 for each additional 1,000 steps but <em>up</em> by 116 for each kilometre. The elastic net method also struggles with these two, but returns values that are better for interpretation.</p>
<p>Here’s the code for obtaining those coefficient estimates:</p>
<figure>
<pre><code data-lang="r"><span># Cross-validated version to determine lambda at best value of alpha, for use later
</span><span>cvmod</span><span> </span><span>&lt;-</span><span> </span><span>cv.glmnet</span><span>(</span><span>X</span><span>,</span><span> </span><span>Y</span><span>,</span><span> </span><span>alpha</span><span> </span><span>=</span><span> </span><span>bestalpha</span><span>)</span><span>
</span><span># the model itself
</span><span>mod1</span><span> </span><span>&lt;-</span><span> </span><span>glmnet</span><span>(</span><span>X</span><span>,</span><span> </span><span>Y</span><span>,</span><span> </span><span>alpha</span><span> </span><span>=</span><span> </span><span>bestalpha</span><span>)</span><span>
</span><span># the OLS model
</span><span>mod2</span><span> </span><span>&lt;-</span><span> </span><span>lm</span><span>(</span><span>Calories.Burned</span><span> </span><span>~</span><span> </span><span>.</span><span>,</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>fitbit2</span><span>)</span><span>
</span><span>coefs</span><span> </span><span>&lt;-</span><span> </span><span>data.frame</span><span>(</span><span>original</span><span> </span><span>=</span><span> </span><span>coef</span><span>(</span><span>mod2</span><span>),</span><span> </span><span>shrunk</span><span> </span><span>=</span><span> </span><span>as.vector</span><span>(</span><span>coef</span><span>(</span><span>mod1</span><span>,</span><span> </span><span>s</span><span> </span><span>=</span><span> </span><span>cvmod</span><span>$</span><span>lambda.min</span><span>)),</span><span>
</span><span>very.shrunk</span><span> </span><span>=</span><span> </span><span>as.vector</span><span>(</span><span>coef</span><span>(</span><span>mod1</span><span>,</span><span> </span><span>s</span><span> </span><span>=</span><span> </span><span>cvmod</span><span>$</span><span>lambda.1se</span><span>)))</span><span>
</span><span>round</span><span>(</span><span>coefs</span><span>,</span><span> </span><span>3</span><span>)</span></code></pre>
</figure>
<h2>Even with <code>lambda=0</code>, elastic net gives different result to OLS due to numerical calculation issues</h2>
<p>I thought I would check what happens if lambda is set to zero so there is no penalty for the total size of coefficients. In theory, the result should be the same as the ordinary least squares estimate, but in the presence of severe collinearity like we have here (between <code>Steps</code> and <code>Distance</code>) this turns out not to be the case.</p>
<figure>
<pre><code data-lang="r"><span>></span><span> </span><span>mod3</span><span> </span><span>&lt;-</span><span> </span><span>glmnet</span><span>(</span><span>X</span><span>,</span><span> </span><span>Y</span><span>,</span><span> </span><span>lambda</span><span> </span><span>=</span><span> </span><span>0</span><span>)</span><span>
</span><span>></span><span> </span><span>round</span><span>(</span><span>data.frame</span><span>(</span><span>&quot;elastic, lambda = 0&quot;</span><span> </span><span>=</span><span> </span><span>as.vector</span><span>(</span><span>coef</span><span>(</span><span>mod3</span><span>)),</span><span> </span><span>&quot;lm&quot;</span><span> </span><span>=</span><span> </span><span>coef</span><span>(</span><span>mod2</span><span>),</span><span> </span><span>check.names</span><span> </span><span>=</span><span> </span><span>FALSE</span><span>),</span><span> </span><span>3</span><span>)</span><span>
</span><span>elastic</span><span>,</span><span> </span><span>lambda</span><span> </span><span>=</span><span> </span><span>0</span><span> </span><span>lm</span><span>
</span><span>(</span><span>Intercept</span><span>)</span><span> </span><span>1937.924</span><span> </span><span>1941.189</span><span>
</span><span>Steps</span><span> </span><span>15.455</span><span> </span><span>-66.827</span><span>
</span><span>Distance</span><span> </span><span>0.653</span><span> </span><span>116.077</span><span>
</span><span>Floors</span><span> </span><span>0.011</span><span> </span><span>0.027</span><span>
</span><span>Minutes.Sedentary</span><span> </span><span>-0.241</span><span> </span><span>-0.246</span><span>
</span><span>Minutes.Lightly.Active</span><span> </span><span>2.236</span><span> </span><span>2.247</span><span>
</span><span>Minutes.Fairly.Active</span><span> </span><span>4.415</span><span> </span><span>4.402</span><span>
</span><span>Minutes.Very.Active</span><span> </span><span>3.894</span><span> </span><span>3.895</span></code></pre>
</figure>
<p>This comes about from numerical issues to do with the estimation method under the hood in <code>glmnet</code> and I won’t try to understand or explain exactly why. If we remove the severe collinearity (which, in a real life situation, I would have done as the very first step when I realised <code>Distance</code> contains no real information additional to <code>Steps</code>) the estimated coefficients are almost identical, as they should be:</p>
<figure>
<pre><code data-lang="r"><span>></span><span> </span><span># what about a simpler case with less numerical instability - </span><span>></span><span> </span><span># drop &quot;distance&quot; (which is the second column)
</span><span>></span><span> </span><span>mod4</span><span> </span><span>&lt;-</span><span> </span><span>glmnet</span><span>(</span><span>X</span><span>[</span><span> </span><span>,</span><span> </span><span>-2</span><span>],</span><span> </span><span>Y</span><span>,</span><span> </span><span>lambda</span><span> </span><span>=</span><span> </span><span>0</span><span>)</span><span>
</span><span>></span><span> </span><span>mod5</span><span> </span><span>&lt;-</span><span> </span><span>lm</span><span>(</span><span>Y</span><span> </span><span>~</span><span> </span><span>X</span><span>[</span><span> </span><span>,</span><span> </span><span>-2</span><span>])</span><span>
</span><span>></span><span> </span><span>round</span><span>(</span><span>data.frame</span><span>(</span><span>&quot;elastic, lambda = 0&quot;</span><span> </span><span>=</span><span> </span><span>as.vector</span><span>(</span><span>coef</span><span>(</span><span>mod4</span><span>)),</span><span> </span><span>&quot;lm&quot;</span><span> </span><span>=</span><span> </span><span>coef</span><span>(</span><span>mod5</span><span>),</span><span> </span><span>check.names</span><span> </span><span>=</span><span> </span><span>FALSE</span><span>),</span><span> </span><span>3</span><span>)</span><span>
</span><span>elastic</span><span>,</span><span> </span><span>lambda</span><span> </span><span>=</span><span> </span><span>0</span><span> </span><span>lm</span><span>
</span><span>(</span><span>Intercept</span><span>)</span><span> </span><span>1938.103</span><span> </span><span>1938.129</span><span>
</span><span>Steps</span><span> </span><span>15.885</span><span> </span><span>15.798</span><span>
</span><span>Floors</span><span> </span><span>0.011</span><span> </span><span>0.010</span><span>
</span><span>Minutes.Sedentary</span><span> </span><span>-0.241</span><span> </span><span>-0.241</span><span>
</span><span>Minutes.Lightly.Active</span><span> </span><span>2.236</span><span> </span><span>2.239</span><span>
</span><span>Minutes.Fairly.Active</span><span> </span><span>4.415</span><span> </span><span>4.413</span><span>
</span><span>Minutes.Very.Active</span><span> </span><span>3.897</span><span> </span><span>3.906</span></code></pre>
</figure>
<h2>In this particular dataset, predictive power is similar in OLS and after elastic net regularisation</h2>
<p>To compare the predictive strength of these different models I used the simple bootstrap, where the modelling approach is applied to bootstrap resamples of the data and the estimated model then used to predict the full original dataset. The elastic net regularized model slightly out-performs the ordinary least squares model with all variables; both significantly out-perform the super-simple one variable model with only <code>Steps</code> in the explanatory side of the formula. Here are the root mean square errors of the three different estimation methods:</p>
<figure>
<pre><code data-lang="r"><span> </span><span>elastic</span><span> </span><span>lm_allvar</span><span> </span><span>lm_1var</span><span> </span><span>95.7</span><span> </span><span>99.5</span><span> </span><span>165.0</span><span> </span></code></pre>
</figure>
<p>Code for the validation process:</p>
<figure>
<pre><code data-lang="r"><span>#----------------validation--------
# function to feed to boot that does elastic modelling
</span><span>mod1f</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>data</span><span>,</span><span> </span><span>i</span><span>){</span><span>
</span><span>X</span><span> </span><span>&lt;-</span><span> </span><span>as.matrix</span><span>(</span><span>data</span><span>[</span><span>i</span><span> </span><span>,</span><span> </span><span>-1</span><span>])</span><span>
</span><span>Y</span><span> </span><span>&lt;-</span><span> </span><span>data</span><span>[</span><span>i</span><span> </span><span>,</span><span> </span><span>1</span><span>]</span><span>
</span><span>cvmod</span><span> </span><span>&lt;-</span><span> </span><span>cv.glmnet</span><span>(</span><span>X</span><span>,</span><span> </span><span>Y</span><span>,</span><span> </span><span>alpha</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span> </span><span>nfolds</span><span> </span><span>=</span><span> </span><span>30</span><span>)</span><span>
</span><span>mod1</span><span> </span><span>&lt;-</span><span> </span><span>glmnet</span><span>(</span><span>X</span><span>,</span><span> </span><span>Y</span><span>,</span><span> </span><span>alpha</span><span> </span><span>=</span><span> </span><span>1</span><span>)</span><span>
</span><span>RMSE</span><span>(</span><span>predict</span><span>(</span><span>mod1</span><span>,</span><span> </span><span>newx</span><span> </span><span>=</span><span> </span><span>as.matrix</span><span>(</span><span>data</span><span>[</span><span> </span><span>,</span><span> </span><span>-1</span><span>]),</span><span> </span><span>s</span><span> </span><span>=</span><span> </span><span>cvmod</span><span>$</span><span>lambda.min</span><span>),</span><span> </span><span>data</span><span>[</span><span> </span><span>,</span><span> </span><span>1</span><span>])</span><span>
</span><span>}</span><span>
</span><span>elastic_boot</span><span> </span><span>&lt;-</span><span> </span><span>boot</span><span>(</span><span>fitbit2</span><span>,</span><span> </span><span>statistic</span><span> </span><span>=</span><span> </span><span>mod1f</span><span>,</span><span> </span><span>R</span><span> </span><span>=</span><span> </span><span>99</span><span>)</span><span>
</span><span># function to feed to boot that does OLS modelling
</span><span>mod2f</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>data</span><span>,</span><span> </span><span>i</span><span>){</span><span>
</span><span>mod2</span><span> </span><span>&lt;-</span><span> </span><span>lm</span><span>(</span><span>Calories.Burned</span><span> </span><span>~</span><span> </span><span>.</span><span>,</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>data</span><span>[</span><span>i</span><span>,</span><span> </span><span>])</span><span>
</span><span>RMSE</span><span>(</span><span>predict</span><span>(</span><span>mod2</span><span>,</span><span> </span><span>newdata</span><span> </span><span>=</span><span> </span><span>data</span><span>),</span><span> </span><span>data</span><span>[</span><span> </span><span>,</span><span> </span><span>1</span><span>])</span><span>
</span><span>}</span><span>
</span><span>lm_boot</span><span> </span><span>&lt;-</span><span> </span><span>boot</span><span>(</span><span>fitbit2</span><span>,</span><span> </span><span>statistic</span><span> </span><span>=</span><span> </span><span>mod2f</span><span>,</span><span> </span><span>R</span><span> </span><span>=</span><span> </span><span>99</span><span>)</span><span>
</span><span># for boot with OLS modelling, only one explanatory variable
</span><span>mod0f</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>data</span><span>,</span><span> </span><span>i</span><span>){</span><span>
</span><span>mod0</span><span> </span><span>&lt;-</span><span> </span><span>lm</span><span>(</span><span>Calories.Burned</span><span> </span><span>~</span><span> </span><span>Steps</span><span>,</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>data</span><span>[</span><span>i</span><span>,</span><span> </span><span>])</span><span>
</span><span>RMSE</span><span>(</span><span>predict</span><span>(</span><span>mod0</span><span>,</span><span> </span><span>newdata</span><span> </span><span>=</span><span> </span><span>data</span><span>),</span><span> </span><span>data</span><span>[</span><span> </span><span>,</span><span> </span><span>1</span><span>])</span><span>
</span><span>}</span><span>
</span><span>lm0_boot</span><span> </span><span>&lt;-</span><span> </span><span>boot</span><span>(</span><span>fitbit2</span><span>,</span><span> </span><span>statistic</span><span> </span><span>=</span><span> </span><span>mod0f</span><span>,</span><span> </span><span>R</span><span> </span><span>=</span><span> </span><span>99</span><span>)</span><span>
</span><span>round</span><span>(</span><span>c</span><span>(</span><span>&quot;elastic&quot;</span><span> </span><span>=</span><span> </span><span>mean</span><span>(</span><span>elastic_boot</span><span>$</span><span>t</span><span>),</span><span> </span><span>&quot;lm_allvar&quot;</span><span> </span><span>=</span><span> </span><span>mean</span><span>(</span><span>lm_boot</span><span>$</span><span>t</span><span>),</span><span>
</span><span>&quot;lm_1var&quot;</span><span> </span><span>=</span><span> </span><span>mean</span><span>(</span><span>lm0_boot</span><span>$</span><span>t</span><span>)),</span><span> </span><span>1</span><span>)</span></code></pre>
</figure>
<h2>“Floors” not needed for a good prediction of calories</h2>
<p>One substantive issue of interest to me that I haven’t mentioned was the relationship of “floors climbed” to total calories. This is of interest because it’s the sort of thing that it’s easy to set a personal target for (“steps” is similar). However, none of the estimation methods showed significant evidence of <code>Floors</code> as an explanatory variable. Basically, the “minutes active” variables (at various levels of intensity) contains all the information needed by fitbit to make its estimates of calories. Walking up stairs contributes to more active minutes as shown by heart rate but other than through that is not directly informative itself.</p>
<p>This can be illustrated by what has become one of the standard plots for this sort of estimation method – a graphic illustration of how the size of coefficients for different variables increase as the penalty against coefficients is relaxed:</p>
<p><img alt="variables" src="http://ellisp.github.io/img/0050-elastic-coefs.svg"></p>
<p>We see from this that even the fullest model leaves a zero coefficient for variable #3, <code>Floors</code> – the horizontal green line. <code>Minutes.Very.Active</code> is the second last variable to be allowed to increase from zero – it contains relatively little information after all the other “minutes” variables are included.</p>
<p>Code for that graphic:</p>
<figure>
<pre><code data-lang="r"><span># refit the model with scaled variables just for the graphic.
# working with scaled variables just for visual distinguishing
# shouldn't make a difference as under the hood glmnet does it anyway
# in the models fit earlier.
</span><span>mod1s</span><span> </span><span>&lt;-</span><span> </span><span>glmnet</span><span>(</span><span>scale</span><span>(</span><span>X</span><span>),</span><span> </span><span>Y</span><span>,</span><span> </span><span>alpha</span><span> </span><span>=</span><span> </span><span>bestalpha</span><span>)</span><span>
</span><span># set a palette
</span><span>thepal</span><span> </span><span>&lt;-</span><span> </span><span>brewer.pal</span><span>(</span><span>7</span><span>,</span><span> </span><span>&quot;Set1&quot;</span><span>)</span><span>
</span><span># set the ordering for ease of reading on graphic
</span><span>ordering</span><span> </span><span>&lt;-</span><span> </span><span>c</span><span>(</span><span>7</span><span>,</span><span>5</span><span>,</span><span>6</span><span>,</span><span>2</span><span>,</span><span>1</span><span>,</span><span>3</span><span>,</span><span>4</span><span>)</span><span>
</span><span>par</span><span>(</span><span>mar</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>5.1</span><span>,</span><span> </span><span>4.1</span><span>,</span><span> </span><span>6.5</span><span>,</span><span> </span><span>1</span><span>),</span><span> </span><span>bg</span><span> </span><span>=</span><span> </span><span>&quot;grey90&quot;</span><span>)</span><span>
</span><span>plot</span><span>(</span><span>mod1s</span><span>,</span><span> </span><span>xvar</span><span> </span><span>=</span><span> </span><span>&quot;dev&quot;</span><span>,</span><span> </span><span>label</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>,</span><span> </span><span>col</span><span> </span><span>=</span><span> </span><span>thepal</span><span>,</span><span> </span><span>lwd</span><span> </span><span>=</span><span> </span><span>2</span><span>,</span><span> </span><span>main</span><span> </span><span>=</span><span> </span><span>&quot;Increasing contribution of different explanatory variables\nas penalty for including them is relaxed&quot;</span><span>)</span><span>
</span><span>legend</span><span>(</span><span>&quot;topleft&quot;</span><span>,</span><span> </span><span>legend</span><span> </span><span>=</span><span> </span><span>colnames</span><span>(</span><span>X</span><span>)[</span><span>ordering</span><span>],</span><span> </span><span>text.col</span><span> </span><span>=</span><span> </span><span>thepal</span><span>[</span><span>ordering</span><span>],</span><span> </span><span>lwd</span><span> </span><span>=</span><span> </span><span>2</span><span>,</span><span> </span><span>bty</span><span> </span><span>=</span><span> </span><span>&quot;n&quot;</span><span>,</span><span> </span><span>col</span><span> </span><span>=</span><span> </span><span>thepal</span><span>[</span><span>ordering</span><span>])</span></code></pre>
</figure>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="https://www.r-bloggers.com/elastic-net-regularization-of-a-model-of-burned-calories/"></div></div><div><a data-url="https://www.r-bloggers.com/elastic-net-regularization-of-a-model-of-burned-calories/" data-counturl="https://www.r-bloggers.com/elastic-net-regularization-of-a-model-of-burned-calories/" data-count="vertical" data-text="Elastic net regularization of a model of burned calories" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://ellisp.github.io/blog/2016/08/13/fitbit-lasso"> Peter's stats stuff - R</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 title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/FMGatgTqnio"></a></div>