<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/practical-machine-learning-with-r-and-python-part-3/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/"> R – Giga thoughts …</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>In this post ‘Practical Machine Learning with R and Python – Part 3’,  I discuss ‘Feature Selection’ methods. This post is a continuation of my 2 earlier posts</p>
<ol>
<li><a rel="nofollow" href="https://gigadom.wordpress.com/2017/10/06/practical-machine-learning-with-r-and-python-part-1/" target="_blank">Practical Machine Learning with R and Python – Part 1</a></li>
<li><a rel="nofollow" href="https://gigadom.wordpress.com/2017/10/13/practical-machine-learning-with-r-and-python-part-2/" target="_blank">Practical Machine Learning with R and Python – Part 2</a></li>
</ol>
<p>While applying Machine Learning techniques, the data set will usually include a large number of predictors for a target variable. It is quite likely, that not all the predictors or feature variables will have an impact on the output. Hence it is becomes necessary to choose only those features which influence the output variable thus simplifying  to a reduced feature set on which to train the ML model on. The techniques that are used are the following</p>
<ul>
<li>Best fit</li>
<li>Forward fit</li>
<li>Backward fit</li>
<li>Ridge Regression or L2 regularization</li>
<li>Lasso or L1 regularization</li>
</ul>
<p>This post includes the equivalent ML code in R and Python.</p>
<p>All these methods remove those features which do not sufficiently influence the output. As in my previous 2 posts on “Practical Machine Learning with R and Python’, this post is largely based on the topics in the following 2 MOOC courses<br>
1. <a rel="nofollow" href="https://lagunita.stanford.edu/courses/HumanitiesandScience/StatLearning/Winter2015/about" target="_blank">Statistical Learning</a>, Prof Trevor Hastie &amp; Prof Robert Tibesherani, Online Stanford<br>
2. <a rel="nofollow" href="https://www.coursera.org/learn/python-machine-learning" target="_blank">Applied Machine Learning in Python</a> Prof Kevyn-Collin Thomson, University Of Michigan, Coursera</p>
<p>You can download this R Markdown file and the associated data from Github – <a rel="nofollow" href="https://github.com/tvganesh/MachineLearning-RandPython-Part3" target="_blank">Machine Learning-RandPython-Part3. </a></p>
<div>
<h1>1.1 Best Fit</h1>
<p>For a dataset with features f1,f2,f3…fn, the ‘Best fit’ approach, chooses all possible combinations of features and creates separate ML models for each of the different combinations. The best fit algotithm then uses some filtering criteria based on Adj Rsquared, Cp, BIC or AIC to pick out the best model among all models.</p>
<p>Since the Best Fit approach searches the entire solution space it is computationally infeasible. The number of models that have to be searched increase exponentially as the number of predictors increase. For ‘p’ predictors a total of <img alt="2^{p}" title="2^{p}" src="https://s0.wp.com/latex.php?latex=2%5E%7Bp%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> ML models have to be searched. This can be shown as follows</p>
<p>There are <img alt="C_{1}" title="C_{1}" src="https://s0.wp.com/latex.php?latex=C_%7B1%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> ways to choose single feature ML models among ‘n’ features, <img alt="C_{2}" title="C_{2}" src="https://s0.wp.com/latex.php?latex=C_%7B2%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> ways to choose 2 feature models among ‘n’ models and so on, or<br>
<img alt="1+C_{1} + C_{2} +... + C_{n}" title="1+C_{1} + C_{2} +... + C_{n}" src="https://s0.wp.com/latex.php?latex=1%2BC_%7B1%7D+%2B+C_%7B2%7D+%2B...+%2B+C_%7Bn%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"><br>
= Total number of models in Best Fit.  Since from Binomial theorem we have<br>
<img alt="(1+x)^{n} = 1+C_{1}x + C_{2}x^{2} +... + C_{n}x^{n}" title="(1+x)^{n} = 1+C_{1}x + C_{2}x^{2} +... + C_{n}x^{n}" src="https://s0.wp.com/latex.php?latex=%281%2Bx%29%5E%7Bn%7D+%3D+1%2BC_%7B1%7Dx+%2B+C_%7B2%7Dx%5E%7B2%7D+%2B...+%2B+C_%7Bn%7Dx%5E%7Bn%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"><br>
When x=1 in the equation (1) above, this becomes<br>
<img alt="2^{n} = 1+C_{1} + C_{2} +... + C_{n}" title="2^{n} = 1+C_{1} + C_{2} +... + C_{n}" src="https://s0.wp.com/latex.php?latex=2%5E%7Bn%7D+%3D+1%2BC_%7B1%7D+%2B+C_%7B2%7D+%2B...+%2B+C_%7Bn%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p>Hence there are <img alt="2^{n}" title="2^{n}" src="https://s0.wp.com/latex.php?latex=2%5E%7Bn%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> models to search amongst in Best Fit. For 10 features this is <img alt="2^{10}" title="2^{10}" src="https://s0.wp.com/latex.php?latex=2%5E%7B10%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> or ~1000 models and for 40 features this becomes <img alt="2^{40} " title="2^{40} " src="https://s0.wp.com/latex.php?latex=2%5E%7B40%7D+&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> which almost 1 trillion. Usually there are datasets with 1000 or maybe even 100000 features and Best fit becomes computationally infeasible.</p>
<p>Anyways I have included the Best Fit approach as I use the Boston crime datasets which is available both the MASS package in R and Sklearn in Python and it has 13 features. Even this small feature set takes a bit of time since the Best fit needs to search among ~<img alt="2^{13}= 8192" title="2^{13}= 8192" src="https://s0.wp.com/latex.php?latex=2%5E%7B13%7D%3D+8192&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">  models</p>
<p>Initially I perform a simple Linear Regression Fit to estimate the features that are statistically insignificant. By looking at the p-values of the features it can be seen that ‘indus’ and ‘age’ features have high p-values and are not significant</p>
</div>
<div>
<h1>1.1a Linear Regression – R code</h1>
<pre><code><span>source</span><span>(</span><span>'RFunctions-1.R'</span><span>)</span></code><code></code><code></code><code></code><code></code><code></code></pre>
<pre><code><span>#Read the Boston crime data</span>
<span>df</span><span>=</span><span>read.csv</span><span>(</span><span>&quot;Boston.csv&quot;</span>,<span>stringsAsFactors</span> <span>=</span> <span>FALSE</span><span>)</span> <span># Data from MASS - SL</span>
<span># Rename the columns</span>
<span>names</span><span>(</span><span>df</span><span>)</span> <span>&lt;-</span><span>c</span><span>(</span><span>&quot;no&quot;</span>,<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;charles&quot;</span>,<span>&quot;nox&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;highways&quot;</span>,<span>&quot;tax&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span><span>)</span>
<span># Select specific columns</span>
<span>df1</span> <span>&lt;-</span> <span>df</span> <span>%>%</span> <span>dplyr</span><span>:</span><span>:</span><span>select</span><span>(</span><span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;charles&quot;</span>,<span>&quot;nox&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
                            <span>&quot;distances&quot;</span>,<span>&quot;highways&quot;</span>,<span>&quot;tax&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span><span>)</span>
<span>dim</span><span>(</span><span>df1</span><span>)</span></code></pre>
<pre><code>## [1] 506  14</code></pre>
<pre><code><span># Linear Regression fit</span>
<span>fit</span> <span>&lt;-</span> <span>lm</span><span>(</span><span>cost</span><span>~</span>. ,<span>data</span><span>=</span><span>df1</span><span>)</span>
<span>summary</span><span>(</span><span>fit</span><span>)</span></code></pre>
<pre><code>## 
## Call:
## lm(formula = cost ~ ., data = df1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.595  -2.730  -0.518   1.777  26.199 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.646e+01  5.103e+00   7.144 3.28e-12 ***
## crimeRate    -1.080e-01  3.286e-02  -3.287 0.001087 ** 
## zone          4.642e-02  1.373e-02   3.382 0.000778 ***
## indus         2.056e-02  6.150e-02   0.334 0.738288    
## charles       2.687e+00  8.616e-01   3.118 0.001925 ** 
## nox          -1.777e+01  3.820e+00  -4.651 4.25e-06 ***
## rooms         3.810e+00  4.179e-01   9.116  &lt; 2e-16 ***
## age           6.922e-04  1.321e-02   0.052 0.958229    
## distances    -1.476e+00  1.995e-01  -7.398 6.01e-13 ***
## highways      3.060e-01  6.635e-02   4.613 5.07e-06 ***
## tax          -1.233e-02  3.760e-03  -3.280 0.001112 ** 
## teacherRatio -9.527e-01  1.308e-01  -7.283 1.31e-12 ***
## color         9.312e-03  2.686e-03   3.467 0.000573 ***
## status       -5.248e-01  5.072e-02 -10.347  &lt; 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.745 on 492 degrees of freedom
## Multiple R-squared:  0.7406, Adjusted R-squared:  0.7338 
## F-statistic: 108.1 on 13 and 492 DF,  p-value: &lt; 2.2e-16</code></pre>
<p>Next we apply the different feature selection models to automatically remove features that are not significant below</p>
</div>
<div>
<h1>1.1a Best Fit – R code</h1>
<p>The Best Fit requires the ‘leaps’ R package</p>
<pre><code><span>library</span><span>(</span><span>leaps</span><span>)</span></code></pre>
<pre><code><span>source</span><span>(</span><span>'RFunctions-1.R'</span><span>)</span></code><code></code><code></code><code></code><code></code></pre>
<pre><code><span>#Read the Boston crime data</span>
<span>df</span><span>=</span><span>read.csv</span><span>(</span><span>&quot;Boston.csv&quot;</span>,<span>stringsAsFactors</span> <span>=</span> <span>FALSE</span><span>)</span> <span># Data from MASS - SL</span>
<span># Rename the columns</span>
<span>names</span><span>(</span><span>df</span><span>)</span> <span>&lt;-</span><span>c</span><span>(</span><span>&quot;no&quot;</span>,<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;charles&quot;</span>,<span>&quot;nox&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;highways&quot;</span>,<span>&quot;tax&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span><span>)</span>
<span># Select specific columns</span>
<span>df1</span> <span>&lt;-</span> <span>df</span> <span>%>%</span> <span>dplyr</span><span>:</span><span>:</span><span>select</span><span>(</span><span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;charles&quot;</span>,<span>&quot;nox&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
                            <span>&quot;distances&quot;</span>,<span>&quot;highways&quot;</span>,<span>&quot;tax&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span><span>)</span>

<span># Perform a best fit</span>
<span>bestFit</span><span>=</span><span>regsubsets</span><span>(</span><span>cost</span><span>~</span>.,<span>df1</span>,<span>nvmax</span><span>=</span><span>13</span><span>)</span>

<span># Generate a summary of the fit</span>
<span>bfSummary</span><span>=</span><span>summary</span><span>(</span><span>bestFit</span><span>)</span>

<span># Plot the Residual Sum of Squares vs number of variables </span>
<span>plot</span><span>(</span><span>bfSummary</span><span>$</span><span>rss</span>,<span>xlab</span><span>=</span><span>&quot;Number of Variables&quot;</span>,<span>ylab</span><span>=</span><span>&quot;RSS&quot;</span>,<span>type</span><span>=</span><span>&quot;l&quot;</span>,<span>main</span><span>=</span><span>&quot;Best fit RSS vs No of features&quot;</span><span>)</span>
<span># Get the index of the minimum value</span>
<span>a</span><span>=</span><span>which.min</span><span>(</span><span>bfSummary</span><span>$</span><span>rss</span><span>)</span>
<span># Mark this in red</span>
<span>points</span><span>(</span><span>a</span>,<span>bfSummary</span><span>$</span><span>rss</span><span>[</span><span>a</span><span>]</span>,<span>col</span><span>=</span><span>&quot;red&quot;</span>,<span>cex</span><span>=</span><span>2</span>,<span>pch</span><span>=</span><span>20</span><span>)</span></code></pre>
<p><img>The plot below shows that the Best fit occurs with all 13 features included. Notice that there is no significant change in RSS from 11 features onward.<img data-orig-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-13.png" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-13.png?w=300" data-attachment-id="5671" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-13.png?w=676&h=483" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/unnamed-chunk-1-1-4/" data-large-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-13.png?w=676&amp;h=483" alt="" data-image-title="unnamed-chunk-1-1" sizes="(max-width: 676px) 100vw, 676px" width="450" data-comments-opened="1" class="wp-image-5671"></p>
<pre><code><span># Plot the CP statistic vs Number of variables</span>
<span>plot</span><span>(</span><span>bfSummary</span><span>$</span><span>cp</span>,<span>xlab</span><span>=</span><span>&quot;Number of Variables&quot;</span>,<span>ylab</span><span>=</span><span>&quot;Cp&quot;</span>,<span>type</span><span>=</span><span>'l'</span>,<span>main</span><span>=</span><span>&quot;Best fit Cp vs No of features&quot;</span><span>)</span>
<span># Find the lowest CP value</span>
<span>b</span><span>=</span><span>which.min</span><span>(</span><span>bfSummary</span><span>$</span><span>cp</span><span>)</span>
<span># Mark this in red</span>
<span>points</span><span>(</span><span>b</span>,<span>bfSummary</span><span>$</span><span>cp</span><span>[</span><span>b</span><span>]</span>,<span>col</span><span>=</span><span>&quot;red&quot;</span>,<span>cex</span><span>=</span><span>2</span>,<span>pch</span><span>=</span><span>20</span><span>)</span></code></pre>
<p><img>Based on Cp metric the best fit occurs at 11 features as seen below. The values of the coefficients are also included below<img data-orig-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-23.png" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-23.png?w=300" data-attachment-id="5672" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-23.png?w=676&h=483" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/unnamed-chunk-1-2-4/" data-large-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-23.png?w=676&amp;h=483" alt="" data-image-title="unnamed-chunk-1-2" sizes="(max-width: 676px) 100vw, 676px" width="450" data-comments-opened="1" class="wp-image-5672"></p>
<pre><code><span># Display the set of features which provide the best fit</span>
<span>coef</span><span>(</span><span>bestFit</span>,<span>b</span><span>)</span></code></pre>
<pre><code>##   (Intercept)     crimeRate          zone       charles           nox 
##  36.341145004  -0.108413345   0.045844929   2.718716303 -17.376023429 
##         rooms     distances      highways           tax  teacherRatio 
##   3.801578840  -1.492711460   0.299608454  -0.011777973  -0.946524570 
##         color        status 
##   0.009290845  -0.522553457</code></pre>
<pre><code><span># Plot the BIC value</span>
<span>plot</span><span>(</span><span>bfSummary</span><span>$</span><span>bic</span>,<span>xlab</span><span>=</span><span>&quot;Number of Variables&quot;</span>,<span>ylab</span><span>=</span><span>&quot;BIC&quot;</span>,<span>type</span><span>=</span><span>'l'</span>,<span>main</span><span>=</span><span>&quot;Best fit BIC vs No of Features&quot;</span><span>)</span>
<span># Find and mark the min value</span>
<span>c</span><span>=</span><span>which.min</span><span>(</span><span>bfSummary</span><span>$</span><span>bic</span><span>)</span>
<span>points</span><span>(</span><span>c</span>,<span>bfSummary</span><span>$</span><span>bic</span><span>[</span><span>c</span><span>]</span>,<span>col</span><span>=</span><span>&quot;red&quot;</span>,<span>cex</span><span>=</span><span>2</span>,<span>pch</span><span>=</span><span>20</span><span>)</span></code></pre>
<p><img><img data-orig-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-32.png" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-32.png?w=300" data-attachment-id="5673" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-32.png?w=676&h=483" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/unnamed-chunk-1-3-3/" data-large-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-32.png?w=676&amp;h=483" alt="" data-image-title="unnamed-chunk-1-3" sizes="(max-width: 676px) 100vw, 676px" width="450" data-comments-opened="1" class="wp-image-5673"></p>
<pre><code><span># R has some other good plots for best fit</span>
<span>plot</span><span>(</span><span>bestFit</span>,<span>scale</span><span>=</span><span>&quot;r2&quot;</span>,<span>main</span><span>=</span><span>&quot;Rsquared vs No Features&quot;</span><span>)</span></code></pre>
<p><img>R has the following set of really nice visualizations. The plot below shows the Rsquared for a set of predictor variables. It can be seen when Rsquared starts at 0.74- indus, charles and age have not been included. <img data-orig-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-42.png" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-42.png?w=300" data-attachment-id="5674" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-42.png?w=676&h=483" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/unnamed-chunk-1-4-3/" data-large-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-42.png?w=676&amp;h=483" alt="" data-image-title="unnamed-chunk-1-4" sizes="(max-width: 676px) 100vw, 676px" width="450" data-comments-opened="1" class="wp-image-5674"></p>
<pre><code><span>plot</span><span>(</span><span>bestFit</span>,<span>scale</span><span>=</span><span>&quot;Cp&quot;</span>,<span>main</span><span>=</span><span>&quot;Cp vs NoFeatures&quot;</span><span>)</span></code></pre>
<p><img>The Cp plot below for value shows indus, charles and age as not included in the Best fit</p>
<p><img data-orig-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-5.png" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-5.png?w=300" data-attachment-id="5675" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-5.png?w=676&h=483" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/unnamed-chunk-1-5/" data-large-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-5.png?w=676&amp;h=483" alt="" data-image-title="unnamed-chunk-1-5" sizes="(max-width: 676px) 100vw, 676px" width="450" data-comments-opened="1" class="wp-image-5675"></p>
<pre><code><span>plot</span><span>(</span><span>bestFit</span>,<span>scale</span><span>=</span><span>&quot;bic&quot;</span>,<span>main</span><span>=</span><span>&quot;BIC vs Features&quot;</span><span>)</span></code></pre>
<p><img><img data-orig-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-6.png" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-6.png?w=300" data-attachment-id="5676" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-6.png?w=676&h=483" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/unnamed-chunk-1-6/" data-large-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-1-6.png?w=676&amp;h=483" alt="" data-image-title="unnamed-chunk-1-6" sizes="(max-width: 676px) 100vw, 676px" width="450" data-comments-opened="1" class="wp-image-5676"></p>
<div>
<h2>1.1b Best fit (Exhaustive Search ) – Python code</h2>
<p>The Python package for performing a Best Fit is the Exhaustive Feature Selector <a rel="nofollow" href="http://rasbt.github.io/mlxtend/user_guide/feature_selection/ExhaustiveFeatureSelector/" target="_blank">EFS</a>.</p>
<pre><code><span>import</span> numpy <span>as</span> np
<span>import</span> pandas <span>as</span> pd
<span>import</span> os
<span>import</span> matplotlib.pyplot <span>as</span> plt
<span>from</span> sklearn.model_selection <span>import</span> train_test_split
<span>from</span> sklearn.linear_model <span>import</span> LinearRegression
<span>from</span> mlxtend.feature_selection <span>import</span> ExhaustiveFeatureSelector <span>as</span> EFS

<span># Read the Boston crime data</span>
df = pd.read_csv(<span>&quot;Boston.csv&quot;</span>,encoding = <span>&quot;ISO-8859-1&quot;</span>)

<span>#Rename the columns</span>
df.columns=[<span>&quot;no&quot;</span>,<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;chasRiver&quot;</span>,<span>&quot;NO2&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;idxHighways&quot;</span>,<span>&quot;taxRate&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span>]
<span># Set X and y </span>
X=df[[<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;chasRiver&quot;</span>,<span>&quot;NO2&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;idxHighways&quot;</span>,<span>&quot;taxRate&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>]]
y=df[<span>'cost'</span>]

<span># Perform an Exhaustive Search. The EFS and SFS packages use 'neg_mean_squared_error'. The 'mean_squared_error' seems to have been deprecated. I think this is just the MSE with the a negative sign.</span>
lr = LinearRegression()
efs1 = EFS(lr, 
           min_features=<span>1</span>,
           max_features=<span>13</span>,
           scoring=<span>'neg_mean_squared_error'</span>,
           print_progress=<span>True</span>,
           cv=<span>5</span>)


<span># Create a efs fit</span>
efs1 = efs1.fit(X.as_matrix(), y.as_matrix())

<span>print</span>(<span>'Best negtive mean squared error: %.2f'</span> % efs1.best_score_)
<span>## Print the IDX of the best features </span>
<span>print</span>(<span>'Best subset:'</span>, efs1.best_idx_)</code><code>
</code></pre>
</div>
</div>
<div>
<pre><code>Features: 8191/8191Best negtive mean squared error: -28.92
## ('Best subset:', (0, 1, 4, 6, 7, 8, 9, 10, 11, 12))</code></pre>
<p>The indices for the best subset are shown above.</p>
<h1>1.2 Forward fit</h1>
<p>Forward fit is a greedy algorithm that tries to optimize the feature selected, by minimizing the selection criteria (adj Rsqaured, Cp, AIC or BIC) at every step. For a dataset with features f1,f2,f3…fn, the forward fit starts with the NULL set. It then pick the ML model with a single feature from n features which has the highest adj Rsquared, or minimum Cp, BIC or some such criteria. After picking the 1 feature from n which satisfies the criteria the most, the next feature from the remaining n-1 features is chosen. When the 2 feature model which satisfies the selection criteria the best is chosen, another feature from the remaining n-2 features are added and so on. The forward fit is a sub-optimal algorithm. There is no guarantee that the final list of features chosen will be the best among the lot. The computation required for this is of  <img alt="n + n-1 + n -2 + .. 1 = n(n+1)/2" title="n + n-1 + n -2 + .. 1 = n(n+1)/2" src="https://s0.wp.com/latex.php?latex=n+%2B+n-1+%2B+n+-2+%2B+..+1+%3D+n%28n%2B1%29%2F2&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> which is of the order of <img alt="n^{2}" title="n^{2}" src="https://s0.wp.com/latex.php?latex=n%5E%7B2%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">. Though forward fit is a sub optimal solution it is far more computationally efficient than best fit</p>
<div>
<h2>1.2a Forward fit – R code</h2>
<p>Forward fit in R determines that 11 features are required for the best fit. The features are shown below</p>
<pre><code><span>library</span><span>(</span><span>leaps</span><span>)</span>
<span># Read the data</span>
<span>df</span><span>=</span><span>read.csv</span><span>(</span><span>&quot;Boston.csv&quot;</span>,<span>stringsAsFactors</span> <span>=</span> <span>FALSE</span><span>)</span> <span># Data from MASS - SL</span>
<span># Rename the columns</span>
<span>names</span><span>(</span><span>df</span><span>)</span> <span>&lt;-</span><span>c</span><span>(</span><span>&quot;no&quot;</span>,<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;charles&quot;</span>,<span>&quot;nox&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;highways&quot;</span>,<span>&quot;tax&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span><span>)</span>

<span># Select columns</span>
<span>df1</span> <span>&lt;-</span> <span>df</span> <span>%>%</span> <span>dplyr</span><span>:</span><span>:</span><span>select</span><span>(</span><span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;charles&quot;</span>,<span>&quot;nox&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
                     <span>&quot;distances&quot;</span>,<span>&quot;highways&quot;</span>,<span>&quot;tax&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span><span>)</span>

<span>#Split as training and test </span>
<span>train_idx</span> <span>&lt;-</span> <span>trainTestSplit</span><span>(</span><span>df1</span>,<span>trainPercent</span><span>=</span><span>75</span>,<span>seed</span><span>=</span><span>5</span><span>)</span>
<span>train</span> <span>&lt;-</span> <span>df1</span><span>[</span><span>train_idx</span>, <span>]</span>
<span>test</span> <span>&lt;-</span> <span>df1</span><span>[</span><span>-</span><span>train_idx</span>, <span>]</span>

<span># Find the best forward fit</span>
<span>fitFwd</span><span>=</span><span>regsubsets</span><span>(</span><span>cost</span><span>~</span>.,<span>data</span><span>=</span><span>train</span>,<span>nvmax</span><span>=</span><span>13</span>,<span>method</span><span>=</span><span>&quot;forward&quot;</span><span>)</span>

<span># Compute the MSE</span>
<span>valErrors</span><span>=</span><span>rep</span><span>(</span><span>NA</span>,<span>13</span><span>)</span>
<span>test.mat</span><span>=</span><span>model.matrix</span><span>(</span><span>cost</span><span>~</span>.,<span>data</span><span>=</span><span>test</span><span>)</span>
<span>for</span><span>(</span><span>i</span> <span>in</span> <span>1</span><span>:</span><span>13</span><span>)</span><span>{</span>
    <span>coefi</span><span>=</span><span>coef</span><span>(</span><span>fitFwd</span>,<span>id</span><span>=</span><span>i</span><span>)</span>
    <span>pred</span><span>=</span><span>test.mat</span><span>[</span>,<span>names</span><span>(</span><span>coefi</span><span>)</span><span>]</span><span>%*%</span><span>coefi</span>
    <span>valErrors</span><span>[</span><span>i</span><span>]</span><span>=</span><span>mean</span><span>(</span><span>(</span><span>test</span><span>$</span><span>cost</span><span>-</span><span>pred</span><span>)</span><span>^</span><span>2</span><span>)</span>
<span>}</span>

<span># Plot the Residual Sum of Squares</span>
<span>plot</span><span>(</span><span>valErrors</span>,<span>xlab</span><span>=</span><span>&quot;Number of Variables&quot;</span>,<span>ylab</span><span>=</span><span>&quot;Validation Error&quot;</span>,<span>type</span><span>=</span><span>&quot;l&quot;</span>,<span>main</span><span>=</span><span>&quot;Forward fit RSS vs No of features&quot;</span><span>)</span>
<span># Gives the index of the minimum value</span>
<span>a</span><span>&lt;-</span><span>which.min</span><span>(</span><span>valErrors</span><span>)</span>
<span>print</span><span>(</span><span>a</span><span>)</span></code></pre>
<pre><code>## [1] 11</code></pre>
<pre><code><span># Highlight the smallest value</span>
<span>points</span><span>(</span><span>c</span>,<span>valErrors</span><span>[</span><span>a</span><span>]</span>,<span>col</span><span>=</span><span>&quot;blue&quot;</span>,<span>cex</span><span>=</span><span>2</span>,<span>pch</span><span>=</span><span>20</span><span>)</span></code></pre>
<p><img>Forward fit R selects 11 predictors as the best ML model to predict the ‘cost’ output variable. The values for these 11 predictors are included below<img data-orig-file="https://gigadom.files.wordpress.com/2017/10/fig2-1.png" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/10/fig2-1.png?w=300" data-attachment-id="5678" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/fig2-1.png?w=676&h=483" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/fig2-1-2/" data-large-file="https://gigadom.files.wordpress.com/2017/10/fig2-1.png?w=676&amp;h=483" alt="" data-image-title="fig2-1" sizes="(max-width: 676px) 100vw, 676px" width="450" data-comments-opened="1" class="wp-image-5678"></p>
<pre><code><span>#Print the 11 ccoefficients</span>
<span>coefi</span><span>=</span><span>coef</span><span>(</span><span>fitFwd</span>,<span>id</span><span>=</span><span>i</span><span>)</span>
<span>coefi</span></code></pre>
<pre><code>##   (Intercept)     crimeRate          zone         indus       charles 
##  2.397179e+01 -1.026463e-01  3.118923e-02  1.154235e-04  3.512922e+00 
##           nox         rooms           age     distances      highways 
## -1.511123e+01  4.945078e+00 -1.513220e-02 -1.307017e+00  2.712534e-01 
##           tax  teacherRatio         color        status 
## -1.330709e-02 -8.182683e-01  1.143835e-02 -3.750928e-01</code></pre>
</div>
<div>
<h2>1.2b Forward fit with Cross Validation – R code</h2>
<p>The Python package SFS includes N Fold Cross Validation errors for forward and backward fit so I decided to add this code to R. This is not available in the ‘leaps’ R package, however the implementation is quite simple. Another implementation is also available at <a rel="nofollow" href="https://lagunita.stanford.edu/courses/HumanitiesandScience/StatLearning/Winter2015/about" target="_blank">Statistical Learning</a>, Prof Trevor Hastie &amp; Prof Robert Tibesherani, Online Stanford 2.</p>
<pre><code><span>library</span><span>(</span><span>dplyr</span><span>)</span>
<span>df</span><span>=</span><span>read.csv</span><span>(</span><span>&quot;Boston.csv&quot;</span>,<span>stringsAsFactors</span> <span>=</span> <span>FALSE</span><span>)</span> <span># Data from MASS - SL</span>
<span>names</span><span>(</span><span>df</span><span>)</span> <span>&lt;-</span><span>c</span><span>(</span><span>&quot;no&quot;</span>,<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;charles&quot;</span>,<span>&quot;nox&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;highways&quot;</span>,<span>&quot;tax&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span><span>)</span>

<span># Select columns</span>
<span>df1</span> <span>&lt;-</span> <span>df</span> <span>%>%</span> <span>dplyr</span><span>:</span><span>:</span><span>select</span><span>(</span><span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;charles&quot;</span>,<span>&quot;nox&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
                     <span>&quot;distances&quot;</span>,<span>&quot;highways&quot;</span>,<span>&quot;tax&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span><span>)</span>

<span>set.seed</span><span>(</span><span>6</span><span>)</span>
<span># Set max number of features</span>
<span>nvmax</span><span>&lt;-</span><span>13</span>
<span>cvError</span> <span>&lt;-</span> <span>NULL</span>
<span># Loop through each features</span>
<span>for</span><span>(</span><span>i</span> <span>in</span> <span>1</span><span>:</span><span>nvmax</span><span>)</span><span>{</span>
    <span># Set no of folds</span>
    <span>noFolds</span><span>=</span><span>5</span>
    <span># Create the rows which fall into different folds from 1..noFolds</span>
    <span>folds</span> <span>=</span> <span>sample</span><span>(</span><span>1</span><span>:</span><span>noFolds</span>, <span>nrow</span><span>(</span><span>df1</span><span>)</span>, <span>replace</span><span>=</span><span>TRUE</span><span>)</span> 
    <span>cv</span><span>&lt;-</span><span>0</span>
    <span># Loop through the folds</span>
    <span>for</span><span>(</span><span>j</span> <span>in</span> <span>1</span><span>:</span><span>noFolds</span><span>)</span><span>{</span>
        <span># The training is all rows for which the row is != j (k-1 folds -> training)</span>
        <span>train</span> <span>&lt;-</span> <span>df1</span><span>[</span><span>folds</span><span>!=</span><span>j</span>,<span>]</span>
        <span># The rows which have j as the index become the test set</span>
        <span>test</span> <span>&lt;-</span> <span>df1</span><span>[</span><span>folds</span><span>==</span><span>j</span>,<span>]</span>
        <span># Create a forward fitting model for this</span>
        <span>fitFwd</span><span>=</span><span>regsubsets</span><span>(</span><span>cost</span><span>~</span>.,<span>data</span><span>=</span><span>train</span>,<span>nvmax</span><span>=</span><span>13</span>,<span>method</span><span>=</span><span>&quot;forward&quot;</span><span>)</span>
        <span># Select the number of features and get the feature coefficients</span>
        <span>coefi</span><span>=</span><span>coef</span><span>(</span><span>fitFwd</span>,<span>id</span><span>=</span><span>i</span><span>)</span>
        <span>#Get the value of the test data</span>
        <span>test.mat</span><span>=</span><span>model.matrix</span><span>(</span><span>cost</span><span>~</span>.,<span>data</span><span>=</span><span>test</span><span>)</span>
        <span># Multiply the tes data with teh fitted coefficients to get the predicted value</span>
        <span># pred = b0 + b1x1+b2x2... b13x13</span>
        <span>pred</span><span>=</span><span>test.mat</span><span>[</span>,<span>names</span><span>(</span><span>coefi</span><span>)</span><span>]</span><span>%*%</span><span>coefi</span>
        <span># Compute mean squared error</span>
        <span>rss</span><span>=</span><span>mean</span><span>(</span><span>(</span><span>test</span><span>$</span><span>cost</span> <span>-</span> <span>pred</span><span>)</span><span>^</span><span>2</span><span>)</span>
        <span># Add all the Cross Validation errors</span>
        <span>cv</span><span>=</span><span>cv</span><span>+</span><span>rss</span>
    <span>}</span>
    <span># Compute the average of MSE for K folds for number of features 'i'</span>
    <span>cvError</span><span>[</span><span>i</span><span>]</span><span>=</span><span>cv</span><span>/</span><span>noFolds</span>
<span>}</span>
<span>a</span> <span>&lt;-</span> <span>seq</span><span>(</span><span>1</span>,<span>13</span><span>)</span>
<span>d</span> <span>&lt;-</span> <span>as.data.frame</span><span>(</span><span>t</span><span>(</span><span>rbind</span><span>(</span><span>a</span>,<span>cvError</span><span>)</span><span>)</span><span>)</span>
<span>names</span><span>(</span><span>d</span><span>)</span> <span>&lt;-</span> <span>c</span><span>(</span><span>&quot;Features&quot;</span>,<span>&quot;CVError&quot;</span><span>)
</span>#Plot the CV Error vs No of Features
<span>ggplot</span><span>(</span><span>d</span>,<span>aes</span><span>(</span><span>x</span><span>=</span><span>Features</span>,<span>y</span><span>=</span><span>CVError</span><span>)</span>,<span>color</span><span>=</span><span>&quot;blue&quot;</span><span>)</span> <span>+</span> <span>geom_point</span><span>(</span><span>)</span> <span>+</span> <span>geom_line</span><span>(</span><span>color</span><span>=</span><span>&quot;blue&quot;</span><span>)</span> <span>+</span>
    <span>xlab</span><span>(</span><span>&quot;No of features&quot;</span><span>)</span> <span>+</span> <span>ylab</span><span>(</span><span>&quot;Cross Validation Error&quot;</span><span>)</span> <span>+</span>
    <span>ggtitle</span><span>(</span><span>&quot;Forward Selection - Cross Valdation Error vs No of Features&quot;</span><span>)</span></code></pre>
<p><img>Forward fit with 5 fold cross validation indicates that all 13 features are required<img data-orig-file="https://gigadom.files.wordpress.com/2017/10/fig3-1.png" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/10/fig3-1.png?w=300" data-attachment-id="5679" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/fig3-1.png?w=676&h=483" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/fig3-1-2/" data-large-file="https://gigadom.files.wordpress.com/2017/10/fig3-1.png?w=676&amp;h=483" alt="" data-image-title="fig3-1" sizes="(max-width: 676px) 100vw, 676px" width="450" data-comments-opened="1" class="wp-image-5679"></p>
<pre><code><span># This gives the index of the minimum value</span>
<span>a</span><span>=</span><span>which.min</span><span>(</span><span>cvError</span><span>)</span>
<span>print</span><span>(</span><span>a</span><span>)</span></code></pre>
<pre><code>## [1] 13</code></pre>
<pre><code><span>#Print the 13 coefficients of these features</span>
<span>coefi</span><span>=</span><span>coef</span><span>(</span><span>fitFwd</span>,<span>id</span><span>=</span><span>a</span><span>)</span>
<span>coefi</span></code></pre>
<pre><code>##   (Intercept)     crimeRate          zone         indus       charles 
##  36.650645380  -0.107980979   0.056237669   0.027016678   4.270631466 
##           nox         rooms           age     distances      highways 
## -19.000715500   3.714720418   0.019952654  -1.472533973   0.326758004 
##           tax  teacherRatio         color        status 
##  -0.011380750  -0.972862622   0.009549938  -0.582159093</code></pre>
</div>
<div>
<h2>1.2c Forward fit – Python code</h2>
<p>The Backward Fit in Python uses the Sequential feature selection (SFS) package (SFS)(<a rel="nofollow" href="https://rasbt.github.io/mlxtend/user_guide/feature_selection/SequentialFeatureSelector/" target="_blank">https://rasbt.github.io/mlxtend/user_guide/feature_selection/SequentialFeatureSelector/</a>)</p>
<p><strong>Note:</strong> The Cross validation error for SFS in Sklearn is negative, possibly because it computes the ‘neg_mean_squared_error’. The earlier ‘mean_squared_error’ in the package seems to have been deprecated. I have taken the -ve of this neg_mean_squared_error. I think this would give mean_squared_error.</p>
<pre><code><span>import</span> numpy <span>as</span> np
<span>import</span> pandas <span>as</span> pd
<span>import</span> os
<span>import</span> matplotlib.pyplot <span>as</span> plt
<span>from</span> sklearn.model_selection <span>import</span> train_test_split
<span>from</span> sklearn.linear_model <span>import</span> LinearRegression
<span>from</span> sklearn.datasets <span>import</span> load_boston
<span>from</span> mlxtend.plotting <span>import</span> plot_sequential_feature_selection <span>as</span> plot_sfs
<span>import</span> matplotlib.pyplot <span>as</span> plt
<span>from</span> mlxtend.feature_selection <span>import</span> SequentialFeatureSelector <span>as</span> SFS
<span>from</span> sklearn.linear_model <span>import</span> LinearRegression


df = pd.read_csv(<span>&quot;Boston.csv&quot;</span>,encoding = <span>&quot;ISO-8859-1&quot;</span>)
<span>#Rename the columns</span>
df.columns=[<span>&quot;no&quot;</span>,<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;chasRiver&quot;</span>,<span>&quot;NO2&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;idxHighways&quot;</span>,<span>&quot;taxRate&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span>]
X=df[[<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;chasRiver&quot;</span>,<span>&quot;NO2&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;idxHighways&quot;</span>,<span>&quot;taxRate&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>]]
y=df[<span>'cost'</span>]
lr = LinearRegression()
<span># Create a forward fit model</span>
sfs = SFS(lr, 
          k_features=(<span>1</span>,<span>13</span>), 
          forward=<span>True</span>, <span># Forward fit</span>
          floating=<span>False</span>, 
          scoring=<span>'neg_mean_squared_error'</span>,
          cv=<span>5</span>)

<span># Fit this on the data</span>
sfs = sfs.fit(X.as_matrix(), y.as_matrix())
<span># Get all the details of the forward fits</span>
a=sfs.get_metric_dict()
n=[]
o=[]

<span># Compute the mean cross validation scores</span>
<span>for</span> i <span>in</span> np.arange(<span>1</span>,<span>13</span>):
    n.append(-np.mean(a[i][<span>'cv_scores'</span>]))  
m=np.arange(<span>1</span>,<span>13</span>)
<span># Get the index of the minimum CV score</span>

<span># Plot the CV scores vs the number of features</span>
fig1=plt.plot(m,n)
fig1=plt.title(<span>'Mean CV Scores vs No of features'</span>)
fig1.figure.savefig(<span>'fig1.png'</span>, bbox_inches=<span>'tight'</span>)

<span>print</span>(pd.DataFrame.from_dict(sfs.get_metric_dict(confidence_interval=<span>0.90</span>)).T)

idx = np.argmin(n)
<span>print</span> <span>&quot;No of features=&quot;</span>,idx
<span>#Get the features indices for the best forward fit and convert to list</span>
b=list(a[idx][<span>'feature_idx'</span>])
<span>print</span>(b)

<span># Index the column names. </span>
<span># Features from forward fit</span>
<span>print</span>(<span>&quot;Features selected in forward fit&quot;</span>)
<span>print</span>(X.columns[b])</code></pre>
<pre><code>##    avg_score ci_bound                                          cv_scores  \
## 1   -42.6185  19.0465  [-23.5582499971, -41.8215743748, -73.993608929...   
## 2   -36.0651  16.3184  [-18.002498199, -40.1507894517, -56.5286659068...   
## 3   -34.1001    20.87  [-9.43012884381, -25.9584955394, -36.184188174...   
## 4   -33.7681  20.1638  [-8.86076528781, -28.650217633, -35.7246353855...   
## 5   -33.6392  20.5271  [-8.90807628524, -28.0684679108, -35.827463022...   
## 6   -33.6276  19.0859  [-9.549485942, -30.9724602876, -32.6689523347,...   
## 7   -32.4082  19.1455  [-10.0177149635, -28.3780298492, -30.926917231...   
## 8   -32.3697   18.533  [-11.1431684243, -27.5765510172, -31.168994094...   
## 9   -32.4016  21.5561  [-10.8972555995, -25.739780653, -30.1837430353...   
## 10  -32.8504  22.6508  [-12.3909282079, -22.1533250755, -33.385407342...   
## 11  -34.1065  24.7019  [-12.6429253721, -22.1676650245, -33.956999528...   
## 12  -35.5814   25.693  [-12.7303397453, -25.0145323483, -34.211898373...   
## 13  -37.1318  23.2657  [-12.4603005692, -26.0486211062, -33.074137979...   
## 
##                                    feature_idx  std_dev  std_err  
## 1                                        (12,)  18.9042  9.45212  
## 2                                     (10, 12)  16.1965  8.09826  
## 3                                  (10, 12, 5)  20.7142  10.3571  
## 4                               (10, 3, 12, 5)  20.0132  10.0066  
## 5                            (0, 10, 3, 12, 5)  20.3738  10.1869  
## 6                         (0, 3, 5, 7, 10, 12)  18.9433  9.47167  
## 7                      (0, 2, 3, 5, 7, 10, 12)  19.0026  9.50128  
## 8                   (0, 1, 2, 3, 5, 7, 10, 12)  18.3946  9.19731  
## 9               (0, 1, 2, 3, 5, 7, 10, 11, 12)  21.3952  10.6976  
## 10           (0, 1, 2, 3, 4, 5, 7, 10, 11, 12)  22.4816  11.2408  
## 11        (0, 1, 2, 3, 4, 5, 6, 7, 10, 11, 12)  24.5175  12.2587  
## 12     (0, 1, 2, 3, 4, 5, 6, 7, 9, 10, 11, 12)  25.5012  12.7506  
## 13  (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12)  23.0919   11.546  
## No of features= 7
## [0, 2, 3, 5, 7, 10, 12]
## #################################################################################
## Features selected in forward fit
## Index([u'crimeRate', u'indus', u'chasRiver', u'rooms', u'distances',
##        u'teacherRatio', u'status'],
##       dtype='object')</code></pre>
<div>
<p><img>The table above shows the average score, 10 fold CV errors, the features included at every step, std. deviation and std. error<img data-orig-file="https://gigadom.files.wordpress.com/2017/10/fig11.png?w=676" data-orig-size="546,435" data-medium-file="https://gigadom.files.wordpress.com/2017/10/fig11.png?w=676?w=300" data-attachment-id="5681" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/fig11.png?w=676" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/fig1-5/" data-large-file="https://gigadom.files.wordpress.com/2017/10/fig11.png?w=676?w=546" alt="" data-image-title="fig1" sizes="(max-width: 546px) 100vw, 546px" data-comments-opened="1" class="wp-image-5681"></p>
<p>The above plot indicates that 8 features provide the lowest Mean CV error</p>
</div>
</div>
<div>
<h2>1.3 Backward Fit</h2>
<p>Backward fit belongs to the class of greedy algorithms which tries to optimize the feature set, by dropping a feature at every stage which results in the worst performance for a given criteria of Adj RSquared, Cp, BIC or AIC. For a dataset with features f1,f2,f3…fn, the backward fit starts with the all the features f1,f2.. fn to begin with. It then pick the ML model with a n-1 features by dropping the feature,<img alt="f_{j}" title="f_{j}" src="https://s0.wp.com/latex.php?latex=f_%7Bj%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, for e.g., the inclusion of which results in the worst performance in adj Rsquared, or minimum Cp, BIC or some such criteria. At every step 1 feature is dopped. There is no guarantee that the final list of features chosen will be the best among the lot. The computation required for this is of <img alt="n + n-1 + n -2 + .. 1 = n(n+1)/2" title="n + n-1 + n -2 + .. 1 = n(n+1)/2" src="https://s0.wp.com/latex.php?latex=n+%2B+n-1+%2B+n+-2+%2B+..+1+%3D+n%28n%2B1%29%2F2&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> which is of the order of <img alt="n^{2}" title="n^{2}" src="https://s0.wp.com/latex.php?latex=n%5E%7B2%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">. Though backward fit is a sub optimal solution it is far more computationally efficient than best fit</p>
</div>
<div>
<h2>1.3a Backward fit – R code</h2>
<pre><code><span>library</span><span>(</span><span>dplyr</span><span>)</span>
<span># Read the data</span>
<span>df</span><span>=</span><span>read.csv</span><span>(</span><span>&quot;Boston.csv&quot;</span>,<span>stringsAsFactors</span> <span>=</span> <span>FALSE</span><span>)</span> <span># Data from MASS - SL</span>
<span># Rename the columns</span>
<span>names</span><span>(</span><span>df</span><span>)</span> <span>&lt;-</span><span>c</span><span>(</span><span>&quot;no&quot;</span>,<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;charles&quot;</span>,<span>&quot;nox&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;highways&quot;</span>,<span>&quot;tax&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span><span>)</span>

<span># Select columns</span>
<span>df1</span> <span>&lt;-</span> <span>df</span> <span>%>%</span> <span>dplyr</span><span>:</span><span>:</span><span>select</span><span>(</span><span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;charles&quot;</span>,<span>&quot;nox&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
                     <span>&quot;distances&quot;</span>,<span>&quot;highways&quot;</span>,<span>&quot;tax&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span><span>)</span>

<span>set.seed</span><span>(</span><span>6</span><span>)</span>
<span># Set max number of features</span>
<span>nvmax</span><span>&lt;-</span><span>13</span>
<span>cvError</span> <span>&lt;-</span> <span>NULL</span>
<span># Loop through each features</span>
<span>for</span><span>(</span><span>i</span> <span>in</span> <span>1</span><span>:</span><span>nvmax</span><span>)</span><span>{</span>
    <span># Set no of folds</span>
    <span>noFolds</span><span>=</span><span>5</span>
    <span># Create the rows which fall into different folds from 1..noFolds</span>
    <span>folds</span> <span>=</span> <span>sample</span><span>(</span><span>1</span><span>:</span><span>noFolds</span>, <span>nrow</span><span>(</span><span>df1</span><span>)</span>, <span>replace</span><span>=</span><span>TRUE</span><span>)</span> 
    <span>cv</span><span>&lt;-</span><span>0</span>
    <span>for</span><span>(</span><span>j</span> <span>in</span> <span>1</span><span>:</span><span>noFolds</span><span>)</span><span>{</span>
        <span># The training is all rows for which the row is != j </span>
        <span>train</span> <span>&lt;-</span> <span>df1</span><span>[</span><span>folds</span><span>!=</span><span>j</span>,<span>]</span>
        <span># The rows which have j as the index become the test set</span>
        <span>test</span> <span>&lt;-</span> <span>df1</span><span>[</span><span>folds</span><span>==</span><span>j</span>,<span>]</span>
        <span># Create a backward fitting model for this</span>
        <span>fitFwd</span><span>=</span><span>regsubsets</span><span>(</span><span>cost</span><span>~</span>.,<span>data</span><span>=</span><span>train</span>,<span>nvmax</span><span>=</span><span>13</span>,<span>method</span><span>=</span><span>&quot;backward&quot;</span><span>)</span>
        <span># Select the number of features and get the feature coefficients</span>
        <span>coefi</span><span>=</span><span>coef</span><span>(</span><span>fitFwd</span>,<span>id</span><span>=</span><span>i</span><span>)</span>
        <span>#Get the value of the test data</span>
        <span>test.mat</span><span>=</span><span>model.matrix</span><span>(</span><span>cost</span><span>~</span>.,<span>data</span><span>=</span><span>test</span><span>)</span>
        <span># Multiply the tes data with teh fitted coefficients to get the predicted value</span>
        <span># pred = b0 + b1x1+b2x2... b13x13</span>
        <span>pred</span><span>=</span><span>test.mat</span><span>[</span>,<span>names</span><span>(</span><span>coefi</span><span>)</span><span>]</span><span>%*%</span><span>coefi</span>
        <span># Compute mean squared error</span>
        <span>rss</span><span>=</span><span>mean</span><span>(</span><span>(</span><span>test</span><span>$</span><span>cost</span> <span>-</span> <span>pred</span><span>)</span><span>^</span><span>2</span><span>)</span>
        <span># Add the Residual sum of square</span>
        <span>cv</span><span>=</span><span>cv</span><span>+</span><span>rss</span>
    <span>}</span>
    <span># Compute the average of MSE for K folds for number of features 'i'</span>
    <span>cvError</span><span>[</span><span>i</span><span>]</span><span>=</span><span>cv</span><span>/</span><span>noFolds</span>
<span>}</span>
<span>a</span> <span>&lt;-</span> <span>seq</span><span>(</span><span>1</span>,<span>13</span><span>)</span>
<span>d</span> <span>&lt;-</span> <span>as.data.frame</span><span>(</span><span>t</span><span>(</span><span>rbind</span><span>(</span><span>a</span>,<span>cvError</span><span>)</span><span>)</span><span>)</span>
<span>names</span><span>(</span><span>d</span><span>)</span> <span>&lt;-</span> <span>c</span><span>(</span><span>&quot;Features&quot;</span>,<span>&quot;CVError&quot;</span><span>)</span>
<span># Plot the Cross Validation Error vs Number of features</span>
<span>ggplot</span><span>(</span><span>d</span>,<span>aes</span><span>(</span><span>x</span><span>=</span><span>Features</span>,<span>y</span><span>=</span><span>CVError</span><span>)</span>,<span>color</span><span>=</span><span>&quot;blue&quot;</span><span>)</span> <span>+</span> <span>geom_point</span><span>(</span><span>)</span> <span>+</span> <span>geom_line</span><span>(</span><span>color</span><span>=</span><span>&quot;blue&quot;</span><span>)</span> <span>+</span>
    <span>xlab</span><span>(</span><span>&quot;No of features&quot;</span><span>)</span> <span>+</span> <span>ylab</span><span>(</span><span>&quot;Cross Validation Error&quot;</span><span>)</span> <span>+</span>
    <span>ggtitle</span><span>(</span><span>&quot;Backward Selection - Cross Valdation Error vs No of Features&quot;</span><span>)</span></code></pre>
<p><img><img data-orig-file="https://gigadom.files.wordpress.com/2017/10/fig4-1.png" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/10/fig4-1.png?w=300" data-attachment-id="5682" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/fig4-1.png?w=676&h=483" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/fig4-1-2/" data-large-file="https://gigadom.files.wordpress.com/2017/10/fig4-1.png?w=676&amp;h=483" alt="" data-image-title="fig4-1" sizes="(max-width: 676px) 100vw, 676px" width="450" data-comments-opened="1" class="wp-image-5682"></p>
<pre><code><span># This gives the index of the minimum value</span>
<span>a</span><span>=</span><span>which.min</span><span>(</span><span>cvError</span><span>)</span>
<span>print</span><span>(</span><span>a</span><span>)</span></code></pre>
<pre><code>## [1] 13</code></pre>
<pre><code><span>#Print the 13 coefficients of these features</span>
<span>coefi</span><span>=</span><span>coef</span><span>(</span><span>fitFwd</span>,<span>id</span><span>=</span><span>a</span><span>)</span>
<span>coefi</span></code></pre>
<pre><code>##   (Intercept)     crimeRate          zone         indus       charles 
##  36.650645380  -0.107980979   0.056237669   0.027016678   4.270631466 
##           nox         rooms           age     distances      highways 
## -19.000715500   3.714720418   0.019952654  -1.472533973   0.326758004 
##           tax  teacherRatio         color        status 
##  -0.011380750  -0.972862622   0.009549938  -0.582159093</code></pre>
</div>
<div>
<p>Backward selection in R also indicates the 13 features and the corresponding coefficients as providing the best fit</p>
<h2>1.3b Backward fit – Python code</h2>
<p>The Backward Fit in Python uses the Sequential feature selection (SFS) package (SFS)(<a rel="nofollow" href="https://rasbt.github.io/mlxtend/user_guide/feature_selection/SequentialFeatureSelector/" target="_blank">https://rasbt.github.io/mlxtend/user_guide/feature_selection/SequentialFeatureSelector/</a>)</p>
<pre><code><span>import</span> numpy <span>as</span> np
<span>import</span> pandas <span>as</span> pd
<span>import</span> os
<span>import</span> matplotlib.pyplot <span>as</span> plt
<span>from</span> sklearn.model_selection <span>import</span> train_test_split
<span>from</span> sklearn.linear_model <span>import</span> LinearRegression
<span>from</span> mlxtend.plotting <span>import</span> plot_sequential_feature_selection <span>as</span> plot_sfs
<span>import</span> matplotlib.pyplot <span>as</span> plt
<span>from</span> mlxtend.feature_selection <span>import</span> SequentialFeatureSelector <span>as</span> SFS
<span>from</span> sklearn.linear_model <span>import</span> LinearRegression

<span># Read the data</span>
df = pd.read_csv(<span>&quot;Boston.csv&quot;</span>,encoding = <span>&quot;ISO-8859-1&quot;</span>)
<span>#Rename the columns</span>
df.columns=[<span>&quot;no&quot;</span>,<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;chasRiver&quot;</span>,<span>&quot;NO2&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;idxHighways&quot;</span>,<span>&quot;taxRate&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span>]
X=df[[<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;chasRiver&quot;</span>,<span>&quot;NO2&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;idxHighways&quot;</span>,<span>&quot;taxRate&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>]]
y=df[<span>'cost'</span>]
lr = LinearRegression()

<span># Create the SFS model</span>
sfs = SFS(lr, 
          k_features=(<span>1</span>,<span>13</span>), 
          forward=<span>False</span>, <span># Backward</span>
          floating=<span>False</span>, 
          scoring=<span>'neg_mean_squared_error'</span>,
          cv=<span>5</span>)

<span># Fit the model</span>
sfs = sfs.fit(X.as_matrix(), y.as_matrix())
a=sfs.get_metric_dict()
n=[]
o=[]

<span># Compute the mean of the validation scores</span>
<span>for</span> i <span>in</span> np.arange(<span>1</span>,<span>13</span>):
    n.append(-np.mean(a[i][<span>'cv_scores'</span>])) 
m=np.arange(<span>1</span>,<span>13</span>)

<span># Plot the Validation scores vs number of features</span>
fig2=plt.plot(m,n)
fig2=plt.title(<span>'Mean CV Scores vs No of features'</span>)
fig2.figure.savefig(<span>'fig2.png'</span>, bbox_inches=<span>'tight'</span>)

<span>print</span>(pd.DataFrame.from_dict(sfs.get_metric_dict(confidence_interval=<span>0.90</span>)).T)

<span># Get the index of minimum cross validation error</span>
idx = np.argmin(n)
<span>print</span> <span>&quot;No of features=&quot;</span>,idx
<span>#Get the features indices for the best forward fit and convert to list</span>
b=list(a[idx][<span>'feature_idx'</span>])
<span># Index the column names. </span>
<span># Features from backward fit</span>
<span>print</span>(<span>&quot;Features selected in bacward fit&quot;</span>)
<span>print</span>(X.columns[b])
</code></pre>
<pre><code>##    avg_score ci_bound                                          cv_scores  \
## 1   -42.6185  19.0465  [-23.5582499971, -41.8215743748, -73.993608929...   
## 2   -36.0651  16.3184  [-18.002498199, -40.1507894517, -56.5286659068...   
## 3   -35.4992  13.9619  [-17.2329292677, -44.4178648308, -51.633177846...   
## 4    -33.463  12.4081  [-20.6415333292, -37.3247852146, -47.479302977...   
## 5   -33.1038  10.6156  [-20.2872309863, -34.6367078466, -45.931870352...   
## 6   -32.0638  10.0933  [-19.4463829372, -33.460638577, -42.726257249,...   
## 7   -30.7133  9.23881  [-19.4425181917, -31.1742902259, -40.531266671...   
## 8   -29.7432  9.84468  [-19.445277268, -30.0641187173, -40.2561247122...   
## 9   -29.0878  9.45027  [-19.3545569877, -30.094768669, -39.7506036377...   
## 10  -28.9225  9.39697  [-18.562171585, -29.968504938, -39.9586835965,...   
## 11  -29.4301  10.8831  [-18.3346152225, -30.3312847532, -45.065432793...   
## 12  -30.4589  11.1486  [-18.493389527, -35.0290639374, -45.1558231765...   
## 13  -37.1318  23.2657  [-12.4603005692, -26.0486211062, -33.074137979...   
## 
##                                    feature_idx  std_dev  std_err  
## 1                                        (12,)  18.9042  9.45212  
## 2                                     (10, 12)  16.1965  8.09826  
## 3                                  (10, 12, 7)  13.8576  6.92881  
## 4                               (12, 10, 4, 7)  12.3154  6.15772  
## 5                            (4, 7, 8, 10, 12)  10.5363  5.26816  
## 6                         (4, 7, 8, 9, 10, 12)  10.0179  5.00896  
## 7                      (1, 4, 7, 8, 9, 10, 12)  9.16981  4.58491  
## 8                  (1, 4, 7, 8, 9, 10, 11, 12)  9.77116  4.88558  
## 9               (0, 1, 4, 7, 8, 9, 10, 11, 12)  9.37969  4.68985  
## 10           (0, 1, 4, 6, 7, 8, 9, 10, 11, 12)   9.3268   4.6634  
## 11        (0, 1, 3, 4, 6, 7, 8, 9, 10, 11, 12)  10.8018  5.40092  
## 12     (0, 1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12)  11.0653  5.53265  
## 13  (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12)  23.0919   11.546  
## No of features= 9
## Features selected in bacward fit
## Index([u'crimeRate', u'zone', u'NO2', u'distances', u'idxHighways', u'taxRate',
##        u'teacherRatio', u'color', u'status'],
##       dtype='object')</code></pre>
<div>
<p><img>The table above shows the average score, 10 fold CV errors, the features included at every step, std. deviation and std. error<img data-orig-file="https://gigadom.files.wordpress.com/2017/10/fig2.png?w=676" data-orig-size="546,435" data-medium-file="https://gigadom.files.wordpress.com/2017/10/fig2.png?w=676?w=300" data-attachment-id="5683" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/fig2.png?w=676" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/fig2-3/" data-large-file="https://gigadom.files.wordpress.com/2017/10/fig2.png?w=676?w=546" alt="" data-image-title="fig2" sizes="(max-width: 546px) 100vw, 546px" data-comments-opened="1" class="wp-image-5683"></p>
<p>Backward fit in Python indicate that 10 features provide the best fit</p>
</div>
</div>
<div>
<h2>1.3c Sequential Floating Forward Selection (SFFS) – Python code</h2>
<p>The Sequential Feature search also includes ‘floating’ variants which include or exclude features conditionally, once they were excluded or included. The SFFS can conditionally include features which were excluded from the previous step, if it results in a better fit. This option will tend to a better solution, than plain simple SFS. These variants are included below</p>
<pre><code><span>import</span> numpy <span>as</span> np
<span>import</span> pandas <span>as</span> pd
<span>import</span> os
<span>import</span> matplotlib.pyplot <span>as</span> plt
<span>from</span> sklearn.model_selection <span>import</span> train_test_split
<span>from</span> sklearn.linear_model <span>import</span> LinearRegression
<span>from</span> sklearn.datasets <span>import</span> load_boston
<span>from</span> mlxtend.plotting <span>import</span> plot_sequential_feature_selection <span>as</span> plot_sfs
<span>import</span> matplotlib.pyplot <span>as</span> plt
<span>from</span> mlxtend.feature_selection <span>import</span> SequentialFeatureSelector <span>as</span> SFS
<span>from</span> sklearn.linear_model <span>import</span> LinearRegression


df = pd.read_csv(<span>&quot;Boston.csv&quot;</span>,encoding = <span>&quot;ISO-8859-1&quot;</span>)
<span>#Rename the columns</span>
df.columns=[<span>&quot;no&quot;</span>,<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;chasRiver&quot;</span>,<span>&quot;NO2&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;idxHighways&quot;</span>,<span>&quot;taxRate&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span>]
X=df[[<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;chasRiver&quot;</span>,<span>&quot;NO2&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;idxHighways&quot;</span>,<span>&quot;taxRate&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>]]
y=df[<span>'cost'</span>]
lr = LinearRegression()

<span># Create the floating forward search</span>
sffs = SFS(lr, 
          k_features=(<span>1</span>,<span>13</span>), 
          forward=<span>True</span>,  <span># Forward</span>
          floating=<span>True</span>,  <span>#Floating</span>
          scoring=<span>'neg_mean_squared_error'</span>,
          cv=<span>5</span>)

<span># Fit a model</span>
sffs = sffs.fit(X.as_matrix(), y.as_matrix())
a=sffs.get_metric_dict()
n=[]
o=[]
<span># Compute mean validation scores</span>
<span>for</span> i <span>in</span> np.arange(<span>1</span>,<span>13</span>):
    n.append(-np.mean(a[i][<span>'cv_scores'</span>])) 
   
    
    
m=np.arange(<span>1</span>,<span>13</span>)


<span># Plot the cross validation score vs number of features</span>
fig3=plt.plot(m,n)
fig3=plt.title(<span>'SFFS:Mean CV Scores vs No of features'</span>)
fig3.figure.savefig(<span>'fig3.png'</span>, bbox_inches=<span>'tight'</span>)

<span>print</span>(pd.DataFrame.from_dict(sffs.get_metric_dict(confidence_interval=<span>0.90</span>)).T)
<span># Get the index of the minimum CV score</span>
idx = np.argmin(n)
<span>print</span> <span>&quot;No of features=&quot;</span>,idx
<span>#Get the features indices for the best forward floating fit and convert to list</span>
b=list(a[idx][<span>'feature_idx'</span>])
<span>print</span>(b)

<span>print</span>(<span>&quot;#################################################################################&quot;</span>)
<span># Index the column names. </span>
<span># Features from forward fit</span>
<span>print</span>(<span>&quot;Features selected in forward fit&quot;</span>)
<span>print</span>(X.columns[b])</code></pre>
<pre><code>##    avg_score ci_bound                                          cv_scores  \
## 1   -42.6185  19.0465  [-23.5582499971, -41.8215743748, -73.993608929...   
## 2   -36.0651  16.3184  [-18.002498199, -40.1507894517, -56.5286659068...   
## 3   -34.1001    20.87  [-9.43012884381, -25.9584955394, -36.184188174...   
## 4   -33.7681  20.1638  [-8.86076528781, -28.650217633, -35.7246353855...   
## 5   -33.6392  20.5271  [-8.90807628524, -28.0684679108, -35.827463022...   
## 6   -33.6276  19.0859  [-9.549485942, -30.9724602876, -32.6689523347,...   
## 7   -32.1834  12.1001  [-17.9491036167, -39.6479234651, -45.470227740...   
## 8   -32.0908  11.8179  [-17.4389015788, -41.2453629843, -44.247557798...   
## 9   -31.0671  10.1581  [-17.2689542913, -37.4379370429, -41.366372300...   
## 10  -28.9225  9.39697  [-18.562171585, -29.968504938, -39.9586835965,...   
## 11  -29.4301  10.8831  [-18.3346152225, -30.3312847532, -45.065432793...   
## 12  -30.4589  11.1486  [-18.493389527, -35.0290639374, -45.1558231765...   
## 13  -37.1318  23.2657  [-12.4603005692, -26.0486211062, -33.074137979...   
## 
##                                    feature_idx  std_dev  std_err  
## 1                                        (12,)  18.9042  9.45212  
## 2                                     (10, 12)  16.1965  8.09826  
## 3                                  (10, 12, 5)  20.7142  10.3571  
## 4                               (10, 3, 12, 5)  20.0132  10.0066  
## 5                            (0, 10, 3, 12, 5)  20.3738  10.1869  
## 6                         (0, 3, 5, 7, 10, 12)  18.9433  9.47167  
## 7                      (0, 1, 2, 3, 7, 10, 12)  12.0097  6.00487  
## 8                   (0, 1, 2, 3, 7, 8, 10, 12)  11.7297  5.86484  
## 9                (0, 1, 2, 3, 7, 8, 9, 10, 12)  10.0822  5.04111  
## 10           (0, 1, 4, 6, 7, 8, 9, 10, 11, 12)   9.3268   4.6634  
## 11        (0, 1, 3, 4, 6, 7, 8, 9, 10, 11, 12)  10.8018  5.40092  
## 12     (0, 1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12)  11.0653  5.53265  
## 13  (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12)  23.0919   11.546  
## No of features= 9
## [0, 1, 2, 3, 7, 8, 9, 10, 12]
## #################################################################################
## Features selected in forward fit
## Index([u'crimeRate', u'zone', u'indus', u'chasRiver', u'distances',
##        u'idxHighways', u'taxRate', u'teacherRatio', u'status'],
##       dtype='object')</code></pre>
<div>
<p><img>The table above shows the average score, 10 fold CV errors, the features included at every step, std. deviation and std. error<img data-orig-file="https://gigadom.files.wordpress.com/2017/10/fig3.png?w=676" data-orig-size="546,435" data-medium-file="https://gigadom.files.wordpress.com/2017/10/fig3.png?w=676?w=300" data-attachment-id="5685" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/fig3.png?w=676" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/fig3-2/" data-large-file="https://gigadom.files.wordpress.com/2017/10/fig3.png?w=676?w=546" alt="" data-image-title="fig3" sizes="(max-width: 546px) 100vw, 546px" data-comments-opened="1" class="wp-image-5685"></p>
<p>SFFS provides the best fit with 10 predictors</p>
</div>
</div>
<div>
<h2>1.3d Sequential Floating Backward Selection (SFBS) – Python code</h2>
<p>The SFBS is an extension of the SBS. Here features that are excluded at any stage can be conditionally included if the resulting feature set gives a better fit.</p>
<pre><code><span>import</span> numpy <span>as</span> np
<span>import</span> pandas <span>as</span> pd
<span>import</span> os
<span>import</span> matplotlib.pyplot <span>as</span> plt
<span>from</span> sklearn.model_selection <span>import</span> train_test_split
<span>from</span> sklearn.linear_model <span>import</span> LinearRegression
<span>from</span> sklearn.datasets <span>import</span> load_boston
<span>from</span> mlxtend.plotting <span>import</span> plot_sequential_feature_selection <span>as</span> plot_sfs
<span>import</span> matplotlib.pyplot <span>as</span> plt
<span>from</span> mlxtend.feature_selection <span>import</span> SequentialFeatureSelector <span>as</span> SFS
<span>from</span> sklearn.linear_model <span>import</span> LinearRegression


df = pd.read_csv(<span>&quot;Boston.csv&quot;</span>,encoding = <span>&quot;ISO-8859-1&quot;</span>)
<span>#Rename the columns</span>
df.columns=[<span>&quot;no&quot;</span>,<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;chasRiver&quot;</span>,<span>&quot;NO2&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;idxHighways&quot;</span>,<span>&quot;taxRate&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span>]
X=df[[<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;chasRiver&quot;</span>,<span>&quot;NO2&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;idxHighways&quot;</span>,<span>&quot;taxRate&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>]]
y=df[<span>'cost'</span>]
lr = LinearRegression()

sffs = SFS(lr, 
          k_features=(<span>1</span>,<span>13</span>), 
          forward=<span>False</span>, <span># Backward</span>
          floating=<span>True</span>, <span># Floating</span>
          scoring=<span>'neg_mean_squared_error'</span>,
          cv=<span>5</span>)

sffs = sffs.fit(X.as_matrix(), y.as_matrix())
a=sffs.get_metric_dict()
n=[]
o=[]
<span># Compute the mean cross validation score</span>
<span>for</span> i <span>in</span> np.arange(<span>1</span>,<span>13</span>):
    n.append(-np.mean(a[i][<span>'cv_scores'</span>]))  
    
m=np.arange(<span>1</span>,<span>13</span>)

fig4=plt.plot(m,n)
fig4=plt.title(<span>'SFBS: Mean CV Scores vs No of features'</span>)
fig4.figure.savefig(<span>'fig4.png'</span>, bbox_inches=<span>'tight'</span>)

<span>print</span>(pd.DataFrame.from_dict(sffs.get_metric_dict(confidence_interval=<span>0.90</span>)).T)

<span># Get the index of the minimum CV score</span>
idx = np.argmin(n)
<span>print</span> <span>&quot;No of features=&quot;</span>,idx
<span>#Get the features indices for the best backward floating fit and convert to list</span>
b=list(a[idx][<span>'feature_idx'</span>])
<span>print</span>(b)

<span>print</span>(<span>&quot;#################################################################################&quot;</span>)
<span># Index the column names. </span>
<span># Features from forward fit</span>
<span>print</span>(<span>&quot;Features selected in backward floating fit&quot;</span>)
<span>print</span>(X.columns[b])</code></pre>
<pre><code>##    avg_score ci_bound                                          cv_scores  \
## 1   -42.6185  19.0465  [-23.5582499971, -41.8215743748, -73.993608929...   
## 2   -36.0651  16.3184  [-18.002498199, -40.1507894517, -56.5286659068...   
## 3   -34.1001    20.87  [-9.43012884381, -25.9584955394, -36.184188174...   
## 4    -33.463  12.4081  [-20.6415333292, -37.3247852146, -47.479302977...   
## 5   -32.3699  11.2725  [-20.8771078371, -34.9825657934, -45.813447203...   
## 6   -31.6742  11.2458  [-20.3082500364, -33.2288990522, -45.535507868...   
## 7   -30.7133  9.23881  [-19.4425181917, -31.1742902259, -40.531266671...   
## 8   -29.7432  9.84468  [-19.445277268, -30.0641187173, -40.2561247122...   
## 9   -29.0878  9.45027  [-19.3545569877, -30.094768669, -39.7506036377...   
## 10  -28.9225  9.39697  [-18.562171585, -29.968504938, -39.9586835965,...   
## 11  -29.4301  10.8831  [-18.3346152225, -30.3312847532, -45.065432793...   
## 12  -30.4589  11.1486  [-18.493389527, -35.0290639374, -45.1558231765...   
## 13  -37.1318  23.2657  [-12.4603005692, -26.0486211062, -33.074137979...   
## 
##                                    feature_idx  std_dev  std_err  
## 1                                        (12,)  18.9042  9.45212  
## 2                                     (10, 12)  16.1965  8.09826  
## 3                                  (10, 12, 5)  20.7142  10.3571  
## 4                               (4, 10, 7, 12)  12.3154  6.15772  
## 5                            (12, 10, 4, 1, 7)  11.1883  5.59417  
## 6                        (4, 7, 8, 10, 11, 12)  11.1618  5.58088  
## 7                      (1, 4, 7, 8, 9, 10, 12)  9.16981  4.58491  
## 8                  (1, 4, 7, 8, 9, 10, 11, 12)  9.77116  4.88558  
## 9               (0, 1, 4, 7, 8, 9, 10, 11, 12)  9.37969  4.68985  
## 10           (0, 1, 4, 6, 7, 8, 9, 10, 11, 12)   9.3268   4.6634  
## 11        (0, 1, 3, 4, 6, 7, 8, 9, 10, 11, 12)  10.8018  5.40092  
## 12     (0, 1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12)  11.0653  5.53265  
## 13  (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12)  23.0919   11.546  
## No of features= 9
## [0, 1, 4, 7, 8, 9, 10, 11, 12]
## #################################################################################
## Features selected in backward floating fit
## Index([u'crimeRate', u'zone', u'NO2', u'distances', u'idxHighways', u'taxRate',
##        u'teacherRatio', u'color', u'status'],
##       dtype='object')</code></pre>
<div>
<p><img>The table above shows the average score, 10 fold CV errors, the features included at every step, std. deviation and std. error<img data-orig-file="https://gigadom.files.wordpress.com/2017/10/fig4.png?w=676" data-orig-size="546,435" data-medium-file="https://gigadom.files.wordpress.com/2017/10/fig4.png?w=676?w=300" data-attachment-id="5686" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/fig4.png?w=676" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/fig4-2/" data-large-file="https://gigadom.files.wordpress.com/2017/10/fig4.png?w=676?w=546" alt="" data-image-title="fig4" sizes="(max-width: 546px) 100vw, 546px" data-comments-opened="1" class="wp-image-5686"></p>
<p>SFBS indicates that 10 features are needed for the best fit</p>
</div>
</div>
</div>
<div>
<h1>1.4 Ridge regression</h1>
<p>In Linear Regression the Residual Sum of Squares (RSS) is given as</p>
<p><img alt="RSS = sum_{i=1}^{n} (y_{i} - beta_{0} - sum_{j=1}^{p}beta_jx_{ij})^{2}" title="RSS = sum_{i=1}^{n} (y_{i} - beta_{0} - sum_{j=1}^{p}beta_jx_{ij})^{2}" src="https://s0.wp.com/latex.php?latex=RSS+%3D+%5Csum_%7Bi%3D1%7D%5E%7Bn%7D+%28y_%7Bi%7D+-+%5Cbeta_%7B0%7D+-+%5Csum_%7Bj%3D1%7D%5E%7Bp%7D%5Cbeta_jx_%7Bij%7D%29%5E%7B2%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"><br>
Ridge regularization =<img alt="sum_{i=1}^{n} (y_{i} - beta_{0} - sum_{j=1}^{p}beta_jx_{ij})^{2} + lambda sum_{j=1}^{p}beta^{2}" title="sum_{i=1}^{n} (y_{i} - beta_{0} - sum_{j=1}^{p}beta_jx_{ij})^{2} + lambda sum_{j=1}^{p}beta^{2}" src="https://s0.wp.com/latex.php?latex=%5Csum_%7Bi%3D1%7D%5E%7Bn%7D+%28y_%7Bi%7D+-+%5Cbeta_%7B0%7D+-+%5Csum_%7Bj%3D1%7D%5E%7Bp%7D%5Cbeta_jx_%7Bij%7D%29%5E%7B2%7D+%2B+%5Clambda+%5Csum_%7Bj%3D1%7D%5E%7Bp%7D%5Cbeta%5E%7B2%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p>where is the regularization or tuning parameter. Increasing increases the penalty on the coefficients thus shrinking them. However in Ridge Regression features that do not influence the target variable will shrink closer to zero but never become zero except for very large values of</p>
<p>Ridge regression in R requires the ‘glmnet’ package</p>
<div>
<h2>1.4a Ridge Regression – R code</h2>
<pre><code><span>library</span><span>(</span><span>glmnet</span><span>)</span></code></pre>
<pre><code><span>library</span><span>(</span><span>dplyr</span><span>)</span>
<span># Read the data</span>
<span>df</span><span>=</span><span>read.csv</span><span>(</span><span>&quot;Boston.csv&quot;</span>,<span>stringsAsFactors</span> <span>=</span> <span>FALSE</span><span>)</span> <span># Data from MASS - SL</span>
<span>#Rename the columns</span>
<span>names</span><span>(</span><span>df</span><span>)</span> <span>&lt;-</span><span>c</span><span>(</span><span>&quot;no&quot;</span>,<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;charles&quot;</span>,<span>&quot;nox&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;highways&quot;</span>,<span>&quot;tax&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span><span>)</span>
<span># Select specific columns</span>
<span>df1</span> <span>&lt;-</span> <span>df</span> <span>%>%</span> <span>dplyr</span><span>:</span><span>:</span><span>select</span><span>(</span><span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;charles&quot;</span>,<span>&quot;nox&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
                            <span>&quot;distances&quot;</span>,<span>&quot;highways&quot;</span>,<span>&quot;tax&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span><span>)</span>

<span># Set X and y as matrices</span>
<span>X</span><span>=</span><span>as.matrix</span><span>(</span><span>df1</span><span>[</span>,<span>1</span><span>:</span><span>13</span><span>]</span><span>)</span>
<span>y</span><span>=</span><span>df1</span><span>$</span><span>cost</span>

<span># Fit a Ridge model</span>
<span>fitRidge</span> <span>&lt;-</span><span>glmnet</span><span>(</span><span>X</span>,<span>y</span>,<span>alpha</span><span>=</span><span>0</span><span>)</span>

<span>#Plot the model where the coefficient shrinkage is plotted vs log lambda</span>
<span>plot</span><span>(</span><span>fitRidge</span>,<span>xvar</span><span>=</span><span>&quot;lambda&quot;</span>,<span>label</span><span>=</span><span>TRUE</span>,<span>main</span><span>=</span> <span>&quot;Ridge regression coefficient shrikage vs log lambda&quot;</span><span>)</span></code></pre>
<p><img>The plot below shows how the 13 coefficients for the 13 predictors vary when lambda is increased. The x-axis includes log (lambda). We can see that increasing lambda from <img alt="10^{2}" title="10^{2}" src="https://s0.wp.com/latex.php?latex=10%5E%7B2%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> to <img alt="10^{6}" title="10^{6}" src="https://s0.wp.com/latex.php?latex=10%5E%7B6%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> significantly shrinks the coefficients. We can draw a vertical line from the x-axis and read the values of the 13 coefficients. Some of them will be close to zero<img data-orig-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-7-11.png" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-7-11.png?w=300" data-attachment-id="5688" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-7-11.png?w=676&h=483" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/unnamed-chunk-7-1-2/" data-large-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-7-11.png?w=676&amp;h=483" alt="" data-image-title="unnamed-chunk-7-1" sizes="(max-width: 676px) 100vw, 676px" width="450" data-comments-opened="1" class="wp-image-5688"></p>
<pre><code><span># Compute the cross validation error</span>
<span>cvRidge</span><span>=</span><span>cv.glmnet</span><span>(</span><span>X</span>,<span>y</span>,<span>alpha</span><span>=</span><span>0</span><span>)</span>

<span>#Plot the cross validation error</span>
<span>plot</span><span>(</span><span>cvRidge</span>, <span>main</span><span>=</span><span>&quot;Ridge regression Cross Validation Error (10 fold)&quot;</span><span>)</span></code></pre>
<p><img>This gives the 10 fold Cross Validation  Error with respect to log (lambda) As lambda increase the MSE increases<img data-orig-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-7-2.png" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-7-2.png?w=300" data-attachment-id="5689" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-7-2.png?w=676&h=483" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/unnamed-chunk-7-2/" data-large-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-7-2.png?w=676&amp;h=483" alt="" data-image-title="unnamed-chunk-7-2" sizes="(max-width: 676px) 100vw, 676px" width="450" data-comments-opened="1" class="wp-image-5689"></p>
</div>
<div>
<h2>1.4a Ridge Regression – Python code</h2>
<p>The coefficient shrinkage for Python can be plotted like R using Least Angle Regression model a.k.a. LARS package. This is included below</p>
<pre><code><span>import</span> numpy <span>as</span> np
<span>import</span> pandas <span>as</span> pd
<span>import</span> os
<span>import</span> matplotlib.pyplot <span>as</span> plt
<span>from</span> sklearn.model_selection <span>import</span> train_test_split


df = pd.read_csv(<span>&quot;Boston.csv&quot;</span>,encoding = <span>&quot;ISO-8859-1&quot;</span>)
<span>#Rename the columns</span>
df.columns=[<span>&quot;no&quot;</span>,<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;chasRiver&quot;</span>,<span>&quot;NO2&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;idxHighways&quot;</span>,<span>&quot;taxRate&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span>]
X=df[[<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;chasRiver&quot;</span>,<span>&quot;NO2&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;idxHighways&quot;</span>,<span>&quot;taxRate&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>]]
y=df[<span>'cost'</span>]
<span>from</span> sklearn.preprocessing <span>import</span> MinMaxScaler
scaler = MinMaxScaler()

<span>from</span> sklearn.linear_model <span>import</span> Ridge
X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                   random_state = <span>0</span>)

<span># Scale the X_train and X_test</span>
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

<span># Fit a ridge regression with alpha=20</span>
linridge = Ridge(alpha=<span>20.0</span>).fit(X_train_scaled, y_train)

<span># Print the training R squared</span>
<span>print</span>(<span>'R-squared score (training): {:.3f}'</span>
     .format(linridge.score(X_train_scaled, y_train)))
<span># Print the test Rsquared</span>
<span>print</span>(<span>'R-squared score (test): {:.3f}'</span>
     .format(linridge.score(X_test_scaled, y_test)))
<span>print</span>(<span>'Number of non-zero features: {}'</span>
     .format(np.sum(linridge.coef_ != <span>0</span>)))

trainingRsquared=[]
testRsquared=[]
<span># Plot the effect of alpha on the test Rsquared</span>
<span>print</span>(<span>'Ridge regression: effect of alpha regularization parameter\n'</span>)
<span># Choose a list of alpha values</span>
<span>for</span> this_alpha <span>in</span> [<span>0.001</span>,.<span>01</span>,.<span>1</span>,<span>0</span>, <span>1</span>, <span>10</span>, <span>20</span>, <span>50</span>, <span>100</span>, <span>1000</span>]:
    linridge = Ridge(alpha = this_alpha).fit(X_train_scaled, y_train)
    <span># Compute training rsquared</span>
    r2_train = linridge.score(X_train_scaled, y_train)
    <span># Compute test rsqaured</span>
    r2_test = linridge.score(X_test_scaled, y_test)
    num_coeff_bigger = np.sum(abs(linridge.coef_) > <span>1.0</span>)
    trainingRsquared.append(r2_train)
    testRsquared.append(r2_test)
    
<span># Create a dataframe</span>
alpha=[<span>0.001</span>,.<span>01</span>,.<span>1</span>,<span>0</span>, <span>1</span>, <span>10</span>, <span>20</span>, <span>50</span>, <span>100</span>, <span>1000</span>]    
trainingRsquared=pd.DataFrame(trainingRsquared,index=alpha)
testRsquared=pd.DataFrame(testRsquared,index=alpha)

<span># Plot training and test R squared as a function of alpha</span>
df3=pd.concat([trainingRsquared,testRsquared],axis=<span>1</span>)
df3.columns=[<span>'trainingRsquared'</span>,<span>'testRsquared'</span>]

fig5=df3.plot()
fig5=plt.title(<span>'Ridge training and test squared error vs Alpha'</span>)
fig5.figure.savefig(<span>'fig5.png'</span>, bbox_inches=<span>'tight'</span>)

<span># Plot the coefficient shrinage using the LARS package</span>

<span>from</span> sklearn <span>import</span> linear_model
<span># #############################################################################</span>
<span># Compute paths</span>

n_alphas = <span>200</span>
alphas = np.logspace(<span>0</span>, <span>8</span>, n_alphas)

coefs = []
<span>for</span> a <span>in</span> alphas:
    ridge = linear_model.Ridge(alpha=a, fit_intercept=<span>False</span>)
    ridge.fit(X_train_scaled, y_train)
    coefs.append(ridge.coef_)

<span># #############################################################################</span>
<span># Display results</span>

ax = plt.gca()

fig6=ax.plot(alphas, coefs)
fig6=ax.set_xscale(<span>'log'</span>)
fig6=ax.set_xlim(ax.get_xlim()[::-<span>1</span>])  <span># reverse axis</span>
fig6=plt.xlabel(<span>'alpha'</span>)
fig6=plt.ylabel(<span>'weights'</span>)
fig6=plt.title(<span>'Ridge coefficients as a function of the regularization'</span>)
fig6=plt.axis(<span>'tight'</span>)
plt.savefig(<span>'fig6.png'</span>, bbox_inches=<span>'tight'</span>)
</code></pre>
<pre><code>## R-squared score (training): 0.620
## R-squared score (test): 0.438
## Number of non-zero features: 13
## Ridge regression: effect of alpha regularization parameter</code></pre>
<p><img><img>The plot below shows the training and test error when increasing the tuning or regularization parameter ‘alpha’<img data-orig-file="https://gigadom.files.wordpress.com/2017/10/fig5.png?w=676" data-orig-size="565,435" data-medium-file="https://gigadom.files.wordpress.com/2017/10/fig5.png?w=676?w=300" data-attachment-id="5705" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/fig5.png?w=676" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/fig5-2/" data-large-file="https://gigadom.files.wordpress.com/2017/10/fig5.png?w=676?w=565" alt="" data-image-title="fig5" sizes="(max-width: 565px) 100vw, 565px" data-comments-opened="1" class="wp-image-5705"></p>
<p>For Python the coefficient shrinkage with LARS must be viewed from right to left, where you have increasing alpha. As alpha increases the coefficients shrink to 0.<img data-orig-file="https://gigadom.files.wordpress.com/2017/10/fig6.png?w=676" data-orig-size="578,459" data-medium-file="https://gigadom.files.wordpress.com/2017/10/fig6.png?w=676?w=300" data-attachment-id="5704" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/fig6.png?w=676" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/fig6-2/" data-large-file="https://gigadom.files.wordpress.com/2017/10/fig6.png?w=676?w=578" alt="" data-image-title="fig6" sizes="(max-width: 578px) 100vw, 578px" data-comments-opened="1" class="wp-image-5704"></p>
</div>
<div>
<h2>1.5 Lasso regularization</h2>
<p>The Lasso is another form of regularization, also known as L1 regularization. Unlike the Ridge Regression where the coefficients of features which do not influence the target tend to zero, in the lasso regualrization the coefficients become 0. The general form of Lasso is as follows</p>
<p><img alt="sum_{i=1}^{n} (y_{i} - beta_{0} - sum_{j=1}^{p}beta_jx_{ij})^{2} + lambda sum_{j=1}^{p}|beta|" title="sum_{i=1}^{n} (y_{i} - beta_{0} - sum_{j=1}^{p}beta_jx_{ij})^{2} + lambda sum_{j=1}^{p}|beta|" src="https://s0.wp.com/latex.php?latex=%5Csum_%7Bi%3D1%7D%5E%7Bn%7D+%28y_%7Bi%7D+-+%5Cbeta_%7B0%7D+-+%5Csum_%7Bj%3D1%7D%5E%7Bp%7D%5Cbeta_jx_%7Bij%7D%29%5E%7B2%7D+%2B+%5Clambda+%5Csum_%7Bj%3D1%7D%5E%7Bp%7D%7C%5Cbeta%7C&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
</div>
<div>
<h2>1.5a Lasso regularization – R code</h2>
<pre><code><span>library</span><span>(</span><span>glmnet</span><span>)</span>
<span>library</span><span>(</span><span>dplyr</span><span>)</span>
<span>df</span><span>=</span><span>read.csv</span><span>(</span><span>&quot;Boston.csv&quot;</span>,<span>stringsAsFactors</span> <span>=</span> <span>FALSE</span><span>)</span> <span># Data from MASS - SL</span>
<span>names</span><span>(</span><span>df</span><span>)</span> <span>&lt;-</span><span>c</span><span>(</span><span>&quot;no&quot;</span>,<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;charles&quot;</span>,<span>&quot;nox&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;highways&quot;</span>,<span>&quot;tax&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span><span>)</span>
<span>df1</span> <span>&lt;-</span> <span>df</span> <span>%>%</span> <span>dplyr</span><span>:</span><span>:</span><span>select</span><span>(</span><span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;charles&quot;</span>,<span>&quot;nox&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
                            <span>&quot;distances&quot;</span>,<span>&quot;highways&quot;</span>,<span>&quot;tax&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span><span>)</span>

<span># Set X and y as matrices</span>
<span>X</span><span>=</span><span>as.matrix</span><span>(</span><span>df1</span><span>[</span>,<span>1</span><span>:</span><span>13</span><span>]</span><span>)</span>
<span>y</span><span>=</span><span>df1</span><span>$</span><span>cost</span>

<span># Fit the lasso model</span>
<span>fitLasso</span> <span>&lt;-</span> <span>glmnet</span><span>(</span><span>X</span>,<span>y</span><span>)</span>
<span># Plot the coefficient shrinkage as a function of log(lambda)</span>
<span>plot</span><span>(</span><span>fitLasso</span>,<span>xvar</span><span>=</span><span>&quot;lambda&quot;</span>,<span>label</span><span>=</span><span>TRUE</span>,<span>main</span><span>=</span><span>&quot;Lasso regularization - Coefficient shrinkage vs log lambda&quot;</span><span>)</span></code></pre>
<p><img>The plot below shows that in L1 regularization the coefficients actually become zero with increasing lambda<img data-orig-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-9-12.png" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-9-12.png?w=300" data-attachment-id="5706" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-9-12.png?w=676&h=483" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/unnamed-chunk-9-1-3/" data-large-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-9-12.png?w=676&amp;h=483" alt="" data-image-title="unnamed-chunk-9-1" sizes="(max-width: 676px) 100vw, 676px" width="450" data-comments-opened="1" class="wp-image-5706"></p>
<pre><code><span># Compute the cross validation error (10 fold)</span>
<span>cvLasso</span><span>=</span><span>cv.glmnet</span><span>(</span><span>X</span>,<span>y</span>,<span>alpha</span><span>=</span><span>0</span><span>)</span>
<span># Plot the cross validation error</span>
<span>plot</span><span>(</span><span>cvLasso</span><span>)</span></code></pre>
<p><img>This gives the MSE for the lasso model<img data-orig-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-9-21.png" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-9-21.png?w=300" data-attachment-id="5707" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-9-21.png?w=676&h=483" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/unnamed-chunk-9-2-2/" data-large-file="https://gigadom.files.wordpress.com/2017/10/unnamed-chunk-9-21.png?w=676&amp;h=483" alt="" data-image-title="unnamed-chunk-9-2" sizes="(max-width: 676px) 100vw, 676px" width="450" data-comments-opened="1" class="wp-image-5707"></p>
</div>
<div>
<h2>1.5 b Lasso regularization – Python code</h2>
<pre><code>import numpy as np
import pandas as pd
import os
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import Lasso
from sklearn.preprocessing import MinMaxScaler
from sklearn import linear_model

scaler = MinMaxScaler()
df = pd.read_csv(&quot;Boston.csv&quot;,encoding = &quot;ISO-8859-1&quot;)
#Rename the columns
df.columns=[&quot;no&quot;,&quot;crimeRate&quot;,&quot;zone&quot;,&quot;indus&quot;,&quot;chasRiver&quot;,&quot;NO2&quot;,&quot;rooms&quot;,&quot;age&quot;,
              &quot;distances&quot;,&quot;idxHighways&quot;,&quot;taxRate&quot;,&quot;teacherRatio&quot;,&quot;color&quot;,&quot;status&quot;,&quot;cost&quot;]
X=df[[&quot;crimeRate&quot;,&quot;zone&quot;,&quot;indus&quot;,&quot;chasRiver&quot;,&quot;NO2&quot;,&quot;rooms&quot;,&quot;age&quot;,
              &quot;distances&quot;,&quot;idxHighways&quot;,&quot;taxRate&quot;,&quot;teacherRatio&quot;,&quot;color&quot;,&quot;status&quot;]]
y=df['cost']
X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                   random_state = 0)

X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

linlasso = Lasso(alpha=0.1, max_iter = 10).fit(X_train_scaled, y_train)

print('Non-zero features: {}'
     .format(np.sum(linlasso.coef_ != 0)))
print('R-squared score (training): {:.3f}'
     .format(linlasso.score(X_train_scaled, y_train)))
print('R-squared score (test): {:.3f}\n'
     .format(linlasso.score(X_test_scaled, y_test)))
print('Features with non-zero weight (sorted by absolute magnitude):')

for e in sorted (list(zip(list(X), linlasso.coef_)),
                key = lambda e: -abs(e[1])):
    if e[1] != 0:
        print('\t{}, {:.3f}'.format(e[0], e[1]))
        

print('Lasso regression: effect of alpha regularization\n\
parameter on number of features kept in final model\n')

trainingRsquared=[]
testRsquared=[]
#for alpha in [0.01,0.05,0.1, 1, 2, 3, 5, 10, 20, 50]:
for alpha in [0.01,0.07,0.05, 0.1, 1,2, 3, 5, 10]:
    linlasso = Lasso(alpha, max_iter = 10000).fit(X_train_scaled, y_train)
    r2_train = linlasso.score(X_train_scaled, y_train)
    r2_test = linlasso.score(X_test_scaled, y_test)
    trainingRsquared.append(r2_train)
    testRsquared.append(r2_test)
    
alpha=[0.01,0.07,0.05, 0.1, 1,2, 3, 5, 10]    
#alpha=[0.01,0.05,0.1, 1, 2, 3, 5, 10, 20, 50]
trainingRsquared=pd.DataFrame(trainingRsquared,index=alpha)
testRsquared=pd.DataFrame(testRsquared,index=alpha)

df3=pd.concat([trainingRsquared,testRsquared],axis=1)
df3.columns=['trainingRsquared','testRsquared']

fig7=df3.plot()
fig7=plt.title('LASSO training and test squared error vs Alpha')
fig7.figure.savefig('fig7.png', bbox_inches='tight')



</code></pre>
<pre><code>## Non-zero features: 7
## R-squared score (training): 0.726
## R-squared score (test): 0.561
## 
## Features with non-zero weight (sorted by absolute magnitude):
##  status, -18.361
##  rooms, 18.232
##  teacherRatio, -8.628
##  taxRate, -2.045
##  color, 1.888
##  chasRiver, 1.670
##  distances, -0.529
## Lasso regression: effect of alpha regularization
## parameter on number of features kept in final model
## 
## Computing regularization path using the LARS ...
## .C:\Users\Ganesh\ANACON~1\lib\site-packages\sklearn\linear_model\coordinate_descent.py:484: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Fitting data with very small alpha may cause precision problems.
##   ConvergenceWarning)</code></pre>
<div>
<p><img>The plot below gives the training and test R squared error<img data-orig-file="https://gigadom.files.wordpress.com/2017/10/fig7.png?w=676" data-orig-size="556,435" data-medium-file="https://gigadom.files.wordpress.com/2017/10/fig7.png?w=676?w=300" data-attachment-id="5708" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/fig7.png?w=676" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/fig7/" data-large-file="https://gigadom.files.wordpress.com/2017/10/fig7.png?w=676?w=556" alt="" data-image-title="fig7" sizes="(max-width: 556px) 100vw, 556px" data-comments-opened="1" class="wp-image-5708"></p>
<p>
</div>
</div>
<div>
<h2>1.5c Lasso coefficient shrinkage – Python code</h2>
<p>To plot the coefficient shrinkage for Lasso the Least Angle Regression model a.k.a. LARS package. This is shown below</p>
<pre><code><span>import</span> numpy <span>as</span> np
<span>import</span> pandas <span>as</span> pd
<span>import</span> os
<span>import</span> matplotlib.pyplot <span>as</span> plt
<span>from</span> sklearn.model_selection <span>import</span> train_test_split
<span>from</span> sklearn.linear_model <span>import</span> Lasso
<span>from</span> sklearn.preprocessing <span>import</span> MinMaxScaler
<span>from</span> sklearn <span>import</span> linear_model
scaler = MinMaxScaler()
df = pd.read_csv(<span>&quot;Boston.csv&quot;</span>,encoding = <span>&quot;ISO-8859-1&quot;</span>)
<span>#Rename the columns</span>
df.columns=[<span>&quot;no&quot;</span>,<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;chasRiver&quot;</span>,<span>&quot;NO2&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;idxHighways&quot;</span>,<span>&quot;taxRate&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>,<span>&quot;cost&quot;</span>]
X=df[[<span>&quot;crimeRate&quot;</span>,<span>&quot;zone&quot;</span>,<span>&quot;indus&quot;</span>,<span>&quot;chasRiver&quot;</span>,<span>&quot;NO2&quot;</span>,<span>&quot;rooms&quot;</span>,<span>&quot;age&quot;</span>,
              <span>&quot;distances&quot;</span>,<span>&quot;idxHighways&quot;</span>,<span>&quot;taxRate&quot;</span>,<span>&quot;teacherRatio&quot;</span>,<span>&quot;color&quot;</span>,<span>&quot;status&quot;</span>]]
y=df[<span>'cost'</span>]
X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                   random_state = <span>0</span>)

X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)


<span>print</span>(<span>&quot;Computing regularization path using the LARS ...&quot;</span>)
alphas, _, coefs = linear_model.lars_path(X_train_scaled, y_train, method=<span>'lasso'</span>, verbose=<span>True</span>)

xx = np.sum(np.abs(coefs.T), axis=<span>1</span>)
xx /= xx[-<span>1</span>]

fig8=plt.plot(xx, coefs.T)

ymin, ymax = plt.ylim()
fig8=plt.vlines(xx, ymin, ymax, linestyle=<span>'dashed'</span>)
fig8=plt.xlabel(<span>'|coef| / max|coef|'</span>)
fig8=plt.ylabel(<span>'Coefficients'</span>)
fig8=plt.title(<span>'LASSO Path - Coefficient Shrinkage vs L1'</span>)
fig8=plt.axis(<span>'tight'</span>)
plt.savefig(<span>'fig8.png'</span>, bbox_inches=<span>'tight'</span>)
</code></pre>
<div><img>This plot show the coefficient shrinkage for lasso.<img data-orig-file="https://gigadom.files.wordpress.com/2017/10/fig8.png?w=676" data-orig-size="578,455" data-medium-file="https://gigadom.files.wordpress.com/2017/10/fig8.png?w=676?w=300" data-attachment-id="5709" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/10/fig8.png?w=676" data-permalink="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/fig8-2/" data-large-file="https://gigadom.files.wordpress.com/2017/10/fig8.png?w=676?w=578" alt="" data-image-title="fig8" sizes="(max-width: 578px) 100vw, 578px" data-comments-opened="1" class="wp-image-5709"></div>
<div>This 3rd part of the series covers the main ‘feature selection’ methods. I hope these posts serve as a quick and useful reference to ML code both for R and Python!</div>
<div></div>
</div>
<div>Stay tuned for further updates to this series!</div>
<div>Watch this space!</div>
<div></div>
<p> </p>
<div>You may also like</div>
<p>1. <a rel="nofollow" href="https://gigadom.wordpress.com/2015/10/02/natural-language-processing-what-would-shakespeare-say/" target="_blank">Natural language processing: What would Shakespeare say?</a><br>
2. <a rel="nofollow" href="https://gigadom.wordpress.com/2016/06/23/introducing-qcsimulator-a-5-qubit-quantum-computing-simulator-in-r/" target="_blank">Introducing QCSimulator: A 5-qubit quantum computing simulator in R</a><br>
3. <a rel="nofollow" href="https://gigadom.wordpress.com/2017/01/05/googlyplus-yorkr-analyzes-ipl-players-teams-matches-with-plots-and-tables/" target="_blank">GooglyPlus: yorkr analyzes IPL players, teams, matches with plots and tables</a><br>
4. <a rel="nofollow" href="https://gigadom.wordpress.com/2017/09/11/my-travels-through-the-realms-of-data-science-machine-learning-deep-learning-and-ai/" target="_blank">My travels through the realms of Data Science, Machine Learning, Deep Learning and (AI)</a><br>
5. <a rel="nofollow" href="https://gigadom.wordpress.com/2011/11/09/experiments-with-deblurring-using-opencv/" target="_blank">Experiments with deblurring using OpenCV</a><br>
6. <a rel="nofollow" href="https://gigadom.wordpress.com/2017/05/22/r-vs-python-different-similarities-and-similar-differences/" target="_blank">R vs Python: Different similarities and similar differences</a></p>
<p>To see all posts see <a rel="nofollow" href="https://gigadom.wordpress.com/aa-2/" target="_blank">Index of posts</a></p>
<div></div>
</div>
<p> <a rel="nofollow" href="http://feeds.wordpress.com/1.0/gocomments/gigadom.wordpress.com/5653/" target="_blank"><img border="0" alt="" src="http://feeds.wordpress.com/1.0/comments/gigadom.wordpress.com/5653/"></a> <img border="0" alt="" width="1" src="https://pixel.wp.com/b.gif?host=gigadom.wordpress.com&blog=14647342&%23038;post=5653&%23038;subd=gigadom&%23038;ref=&%23038;feed=1" height="1"></p>
	<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://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/"> R – Giga thoughts …</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 alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/A14oYFirZe0" height="1"></a></div>