<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-5/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/"> R – Giga thoughts …</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>This is the 5th and probably penultimate part of my series on ‘Practical Machine Learning with R and Python’. The earlier parts of this series included</p>
<p>1. <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> In this initial post, I touch upon univariate, multivariate, polynomial regression and KNN regression in R and Python<br>
2.<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> In this post, I discuss Logistic Regression, KNN classification and cross validation error for both LOOCV and K-Fold in both R and Python<br>
3.<a rel="nofollow" href="https://gigadom.wordpress.com/2017/10/20/practical-machine-learning-with-r-and-python-part-3/" target="_blank">Practical Machine Learning with R and Python – Part 3</a> This post covered ‘feature selection’ in Machine Learning. Specifically I touch best fit, forward fit, backward fit, ridge(L2 regularization) &amp; lasso (L1 regularization). The post includes equivalent code in R and Python.<br>
4.<a rel="nofollow" href="https://gigadom.wordpress.com/2017/10/29/practical-machine-learning-with-r-and-python-part-4/" target="_blank">Practical Machine Learning with R and Python – Part 4</a> In this part I discussed SVMs, Decision Trees, validation, precision recall, and roc curves</p>
<p>This post ‘Practical Machine Learning with R and Python – Part 5’ discusses regression with B-splines, natural splines, smoothing splines, generalized additive models (GAMS), bagging, random forest and boosting</p>
<p>As with my previous posts in this series, this post is largely based on the following 2 MOOC courses</p>
<p>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 associated data files from Github at <a rel="nofollow" href="https://github.com/tvganesh/MachineLearning-RandPython-Part5" target="_blank">MachineLearning-RandPython-Part5</a></p>
<p>For this part I have used the data sets from <a rel="nofollow" href="https://archive.ics.uci.edu/ml/datasets.html" target="_blank">UCI Machine Learning repository</a>(<a rel="nofollow" href="https://archive.ics.uci.edu/ml/datasets/Communities+and+Crime" target="_blank">Communities and Crime</a> and <a rel="nofollow" href="https://archive.ics.uci.edu/ml/datasets/Auto+MPG" target="_blank">Auto MPG</a>)</p>
<div>
<h2>1. Splines</h2>
<p>When performing regression (continuous or logistic) between a target variable and a feature (or a set of features), a single polynomial for the entire range of the data set usually does not perform a good fit.Rather we would need to provide we could fit<br>
regression curves for different section of the data set.</p>
<p>There are several techniques which do this for e.g. piecewise-constant functions, piecewise-linear functions, piecewise-quadratic/cubic/4th order polynomial functions etc. One such set of functions are the cubic splines which fit cubic polynomials to successive sections of the dataset. The points where the cubic splines join, are called ‘knots’.</p>
<p>Since each section has a different cubic spline, there could be discontinuities (or breaks) at these knots. To prevent these discontinuities ‘natural splines’ and ‘smoothing splines’ ensure that the seperate cubic functions have 2nd order continuity at these knots with the adjacent splines. 2nd order continuity implies that the value, 1st order derivative and 2nd order derivative at these knots are equal.</p>
<p>A cubic spline with knots <img alt="alpha_{k}" title="alpha_{k}" src="https://s0.wp.com/latex.php?latex=%5Calpha_%7Bk%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> , k=1,2,3,..K is a piece-wise cubic polynomial with continuous derivative up to order 2 at each knot. We can write <img alt="y_{i} = beta_{0} +beta_{1}b_{1}(x_{i}) +beta_{2}b_{2}(x_{i}) + .. + beta_{K+3}b_{K+3}(x_{i}) + epsilon_{i}" title="y_{i} = beta_{0} +beta_{1}b_{1}(x_{i}) +beta_{2}b_{2}(x_{i}) + .. + beta_{K+3}b_{K+3}(x_{i}) + epsilon_{i}" src="https://s0.wp.com/latex.php?latex=y_%7Bi%7D+%3D+%5Cbeta_%7B0%7D+%2B%5Cbeta_%7B1%7Db_%7B1%7D%28x_%7Bi%7D%29+%2B%5Cbeta_%7B2%7Db_%7B2%7D%28x_%7Bi%7D%29+%2B+..+%2B+%5Cbeta_%7BK%2B3%7Db_%7BK%2B3%7D%28x_%7Bi%7D%29+%2B+%5Cepsilon_%7Bi%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">.<br>
For each (<img alt="x{i},y{i}" title="x{i},y{i}" src="https://s0.wp.com/latex.php?latex=x%7Bi%7D%2Cy%7Bi%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">), <img alt="b_{i}" title="b_{i}" src="https://s0.wp.com/latex.php?latex=b_%7Bi%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> are called ‘basis’ functions, where  <img alt="b_{1}(x_{i})=x_{i}" title="b_{1}(x_{i})=x_{i}" src="https://s0.wp.com/latex.php?latex=b_%7B1%7D%28x_%7Bi%7D%29%3Dx_%7Bi%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">,  <img alt="b_{2}(x_{i})=x_{i}^2" title="b_{2}(x_{i})=x_{i}^2" src="https://s0.wp.com/latex.php?latex=b_%7B2%7D%28x_%7Bi%7D%29%3Dx_%7Bi%7D%5E2&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, <img alt="b_{3}(x_{i})=x_{i}^3" title="b_{3}(x_{i})=x_{i}^3" src="https://s0.wp.com/latex.php?latex=b_%7B3%7D%28x_%7Bi%7D%29%3Dx_%7Bi%7D%5E3&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, <img alt="b_{k+3}(x_{i})=(x_{i} -alpha_{k})^3" title="b_{k+3}(x_{i})=(x_{i} -alpha_{k})^3" src="https://s0.wp.com/latex.php?latex=b_%7Bk%2B3%7D%28x_%7Bi%7D%29%3D%28x_%7Bi%7D+-%5Calpha_%7Bk%7D%29%5E3&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> where k=1,2,3… K The 1st and 2nd derivatives of cubic splines are continuous at the knots. Hence splines provide a smooth continuous fit to the data by fitting different splines to different sections of the data</p>
</div>
<div>
<h2>1.1a Fit a 4th degree polynomial – R code</h2>
<p>In the code below a non-linear function (a 4th order polynomial) is used to fit the data. Usually when we fit a single polynomial to the entire data set the tails of the fit tend to vary a lot particularly if there are fewer points at the ends. Splines help in reducing this variation at the extremities</p>
<pre><code><span>library</span><span>(</span><span>dplyr</span><span>)</span>
<span>library</span><span>(</span><span>ggplot2</span><span>)</span>
<span>source</span><span>(</span><span>'RFunctions-1.R'</span><span>)</span>
<span># Read the data</span>
<span>df</span><span>=</span><span>read.csv</span><span>(</span><span>&quot;auto_mpg.csv&quot;</span>,<span>stringsAsFactors</span> <span>=</span> <span>FALSE</span><span>)</span> <span># Data from UCI</span>
<span>df1</span> <span>&lt;-</span> <span>as.data.frame</span><span>(</span><span>sapply</span><span>(</span><span>df</span>,<span>as.numeric</span><span>)</span><span>)</span>
<span>#Select specific columns</span>
<span>df2</span> <span>&lt;-</span> <span>df1</span> <span>%>%</span> <span>dplyr</span><span>:</span><span>:</span><span>select</span><span>(</span><span>cylinder</span>,<span>displacement</span>, <span>horsepower</span>,<span>weight</span>, <span>acceleration</span>, <span>year</span>,<span>mpg</span><span>)</span>
<span>auto</span> <span>&lt;-</span> <span>df2</span><span>[</span><span>complete.cases</span><span>(</span><span>df2</span><span>)</span>,<span>]</span>
<span># Fit a 4th degree polynomial</span>
<span>fit</span><span>=</span><span>lm</span><span>(</span><span>mpg</span><span>~</span><span>poly</span><span>(</span><span>horsepower</span>,<span>4</span><span>)</span>,<span>data</span><span>=</span><span>auto</span><span>)</span>
<span>#Display a summary of fit</span>
<span>summary</span><span>(</span><span>fit</span><span>)</span></code></pre>
<pre><code>## 
## Call:
## lm(formula = mpg ~ poly(horsepower, 4), data = auto)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.8820  -2.5802  -0.1682   2.2100  16.1434 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            23.4459     0.2209 106.161   &lt;2e-16 ***
## poly(horsepower, 4)1 -120.1377     4.3727 -27.475   &lt;2e-16 ***
## poly(horsepower, 4)2   44.0895     4.3727  10.083   &lt;2e-16 ***
## poly(horsepower, 4)3   -3.9488     4.3727  -0.903    0.367    
## poly(horsepower, 4)4   -5.1878     4.3727  -1.186    0.236    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.373 on 387 degrees of freedom
## Multiple R-squared:  0.6893, Adjusted R-squared:  0.6861 
## F-statistic: 214.7 on 4 and 387 DF,  p-value: &lt; 2.2e-16</code></pre>
<pre><code><span>#Get the range of horsepower</span>
<span>hp</span> <span>&lt;-</span> <span>range</span><span>(</span><span>auto</span><span>$</span><span>horsepower</span><span>)</span>
<span>#Create a sequence to be used for plotting</span>
<span>hpGrid</span> <span>&lt;-</span> <span>seq</span><span>(</span><span>hp</span><span>[</span><span>1</span><span>]</span>,<span>hp</span><span>[</span><span>2</span><span>]</span>,<span>by</span><span>=</span><span>10</span><span>)</span>
<span>#Predict for these values of horsepower. Set Standard error as TRUE</span>
<span>pred</span><span>=</span><span>predict</span><span>(</span><span>fit</span>,<span>newdata</span><span>=</span><span>list</span><span>(</span><span>horsepower</span><span>=</span><span>hpGrid</span><span>)</span>,<span>se</span><span>=</span><span>TRUE</span><span>)</span>
<span>#Compute bands on either side that is 2xSE</span>
<span>seBands</span><span>=</span><span>cbind</span><span>(</span><span>pred</span><span>$</span><span>fit</span><span>+</span><span>2</span><span>*</span><span>pred</span><span>$</span><span>se.fit</span>,<span>pred</span><span>$</span><span>fit</span><span>-</span><span>2</span><span>*</span><span>pred</span><span>$</span><span>se.fit</span><span>)</span>
<span>#Plot the fit with Standard Error bands</span>
<span>plot</span><span>(</span><span>auto</span><span>$</span><span>horsepower</span>,<span>auto</span><span>$</span><span>mpg</span>,<span>xlim</span><span>=</span><span>hp</span>,<span>cex</span><span>=</span><span>.5</span>,<span>col</span><span>=</span><span>&quot;black&quot;</span>,<span>xlab</span><span>=</span><span>&quot;Horsepower&quot;</span>,
     <span>ylab</span><span>=</span><span>&quot;MPG&quot;</span>, <span>main</span><span>=</span><span>&quot;Polynomial of degree 4&quot;</span><span>)</span>
<span>lines</span><span>(</span><span>hpGrid</span>,<span>pred</span><span>$</span><span>fit</span>,<span>lwd</span><span>=</span><span>2</span>,<span>col</span><span>=</span><span>&quot;blue&quot;</span><span>)</span>
<span>matlines</span><span>(</span><span>hpGrid</span>,<span>seBands</span>,<span>lwd</span><span>=</span><span>2</span>,<span>col</span><span>=</span><span>&quot;blue&quot;</span>,<span>lty</span><span>=</span><span>3</span><span>)</span></code></pre>
<p><img><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig1-1.png?w=456" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig1-1.png?w=456?w=300" data-attachment-id="5763" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig1-1.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig1-1-4/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig1-1.png?w=456?w=676" alt="fig1-1" data-image-title="fig1-1" sizes="(max-width: 676px) 100vw, 676px" data-comments-opened="1" class="wp-image-5763" data-recalc-dims="1"></p>
</div>
<div>
<h2>1.1b Fit a 4th degree polynomial – Python code</h2>
<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.preprocessing <span>import</span> PolynomialFeatures
<span>from</span> sklearn.linear_model <span>import</span> LinearRegression
<span>#Read the auto data</span>
autoDF =pd.read_csv(<span>&quot;auto_mpg.csv&quot;</span>,encoding=<span>&quot;ISO-8859-1&quot;</span>)
<span># Select columns</span>
autoDF1=autoDF[[<span>'mpg'</span>,<span>'cylinder'</span>,<span>'displacement'</span>,<span>'horsepower'</span>,<span>'weight'</span>,<span>'acceleration'</span>,<span>'year'</span>]]
<span># Convert all columns to numeric</span>
autoDF2 = autoDF1.apply(pd.to_numeric, errors=<span>'coerce'</span>)

<span>#Drop NAs</span>
autoDF3=autoDF2.dropna()
autoDF3.shape
X=autoDF3[[<span>'horsepower'</span>]]
y=autoDF3[<span>'mpg'</span>]
<span>#Create a polynomial of degree 4</span>
poly = PolynomialFeatures(degree=<span>4</span>)
X_poly = poly.fit_transform(X)

<span># Fit a polynomial regression line</span>
linreg = LinearRegression().fit(X_poly, y)
<span># Create a range of values</span>
hpGrid = np.arange(np.min(X),np.max(X),<span>10</span>)
hp=hpGrid.reshape(-<span>1</span>,<span>1</span>)
<span># Transform to 4th degree</span>
poly = PolynomialFeatures(degree=<span>4</span>)
hp_poly = poly.fit_transform(hp)

<span>#Create a scatter plot</span>
plt.scatter(X,y)
<span># Fit the prediction</span>
ypred=linreg.predict(hp_poly)
plt.title(<span>&quot;Poylnomial of degree 4&quot;</span>)
fig2=plt.xlabel(<span>&quot;Horsepower&quot;</span>)
fig2=plt.ylabel(<span>&quot;MPG&quot;</span>)
<span># Draw the regression curve</span>
plt.plot(hp,ypred,c=<span>&quot;red&quot;</span>)
plt.savefig(<span>'fig1.png'</span>, bbox_inches=<span>'tight'</span>)</code></pre>
</div>
<div>
<h2><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig1.png?w=456" data-orig-size="566,455" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig1.png?w=456?w=300" data-attachment-id="5765" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig1.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig1-7/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig1.png?w=456?w=566" alt="fig1" data-image-title="fig1" sizes="(max-width: 566px) 100vw, 566px" data-comments-opened="1" class="wp-image-5765" data-recalc-dims="1"></h2>
<h2>1.1c Fit a B-Spline – R Code</h2>
<p>In the code below a B- Spline is fit to data. The B-spline requires the manual selection of knots</p>
<pre><code><span>#Splines</span>
<span>library</span><span>(</span><span>splines</span><span>)</span>
<span># Fit a B-spline to the data. Select knots at 60,75,100,150</span>
<span>fit</span><span>=</span><span>lm</span><span>(</span><span>mpg</span><span>~</span><span>bs</span><span>(</span><span>horsepower</span>,<span>df</span><span>=</span><span>6</span>,<span>knots</span><span>=</span><span>c</span><span>(</span><span>60</span>,<span>75</span>,<span>100</span>,<span>150</span><span>)</span><span>)</span>,<span>data</span><span>=</span><span>auto</span><span>)</span>
<span># Use the fitted regresion to predict</span>
<span>pred</span><span>=</span><span>predict</span><span>(</span><span>fit</span>,<span>newdata</span><span>=</span><span>list</span><span>(</span><span>horsepower</span><span>=</span><span>hpGrid</span><span>)</span>,<span>se</span><span>=</span><span>T</span><span>)</span>
<span># Create a scatter plot</span>
<span>plot</span><span>(</span><span>auto</span><span>$</span><span>horsepower</span>,<span>auto</span><span>$</span><span>mpg</span>,<span>xlim</span><span>=</span><span>hp</span>,<span>cex</span><span>=</span><span>.5</span>,<span>col</span><span>=</span><span>&quot;black&quot;</span>,<span>xlab</span><span>=</span><span>&quot;Horsepower&quot;</span>,
     <span>ylab</span><span>=</span><span>&quot;MPG&quot;</span>, <span>main</span><span>=</span><span>&quot;B-Spline with 4 knots&quot;</span><span>)</span>
<span>#Draw lines with 2 Standard Errors on either side</span>
<span>lines</span><span>(</span><span>hpGrid</span>,<span>pred</span><span>$</span><span>fit</span>,<span>lwd</span><span>=</span><span>2</span><span>)</span>
<span>lines</span><span>(</span><span>hpGrid</span>,<span>pred</span><span>$</span><span>fit</span><span>+</span><span>2</span><span>*</span><span>pred</span><span>$</span><span>se</span>,<span>lty</span><span>=</span><span>&quot;dashed&quot;</span><span>)</span>
<span>lines</span><span>(</span><span>hpGrid</span>,<span>pred</span><span>$</span><span>fit</span><span>-</span><span>2</span><span>*</span><span>pred</span><span>$</span><span>se</span>,<span>lty</span><span>=</span><span>&quot;dashed&quot;</span><span>)</span>
<span>abline</span><span>(</span><span>v</span><span>=</span><span>c</span><span>(</span><span>60</span>,<span>75</span>,<span>100</span>,<span>150</span><span>)</span>,<span>lty</span><span>=</span><span>2</span>,<span>col</span><span>=</span><span>&quot;darkgreen&quot;</span><span>)</span></code></pre>
<p><img><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig2-1.png?w=456" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig2-1.png?w=456?w=300" data-attachment-id="5768" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig2-1.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig2-1-4/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig2-1.png?w=456?w=676" alt="fig2-1" data-image-title="fig2-1" sizes="(max-width: 676px) 100vw, 676px" data-comments-opened="1" class="wp-image-5768" data-recalc-dims="1"></p>
</div>
<div>
<h2>1.1d Fit a Natural Spline – R Code</h2>
<p>Here a ‘Natural Spline’ is used to fit .The Natural Spline extrapolates beyond the boundary knots and the ends of the function are much more constrained than a regular spline or a global polynomoial where the ends can wag a lot more. Natural splines do not require the explicit selection of knots</p>
<pre><code><span># There is no need to select the knots here. There is a smoothing parameter which</span>
<span># can be specified by the degrees of freedom 'df' parameter. The natural spline</span>

<span>fit2</span><span>=</span><span>lm</span><span>(</span><span>mpg</span><span>~</span><span>ns</span><span>(</span><span>horsepower</span>,<span>df</span><span>=</span><span>4</span><span>)</span>,<span>data</span><span>=</span><span>auto</span><span>)</span>
<span>pred</span><span>=</span><span>predict</span><span>(</span><span>fit2</span>,<span>newdata</span><span>=</span><span>list</span><span>(</span><span>horsepower</span><span>=</span><span>hpGrid</span><span>)</span>,<span>se</span><span>=</span><span>T</span><span>)</span>
<span>plot</span><span>(</span><span>auto</span><span>$</span><span>horsepower</span>,<span>auto</span><span>$</span><span>mpg</span>,<span>xlim</span><span>=</span><span>hp</span>,<span>cex</span><span>=</span><span>.5</span>,<span>col</span><span>=</span><span>&quot;black&quot;</span>,<span>xlab</span><span>=</span><span>&quot;Horsepower&quot;</span>,
     <span>ylab</span><span>=</span><span>&quot;MPG&quot;</span>, <span>main</span><span>=</span><span>&quot;Natural Splines&quot;</span><span>)</span>
<span>lines</span><span>(</span><span>hpGrid</span>,<span>pred</span><span>$</span><span>fit</span>,<span>lwd</span><span>=</span><span>2</span><span>)</span>
<span>lines</span><span>(</span><span>hpGrid</span>,<span>pred</span><span>$</span><span>fit</span><span>+</span><span>2</span><span>*</span><span>pred</span><span>$</span><span>se</span>,<span>lty</span><span>=</span><span>&quot;dashed&quot;</span><span>)</span>
<span>lines</span><span>(</span><span>hpGrid</span>,<span>pred</span><span>$</span><span>fit</span><span>-</span><span>2</span><span>*</span><span>pred</span><span>$</span><span>se</span>,<span>lty</span><span>=</span><span>&quot;dashed&quot;</span><span>)</span></code></pre>
<p><img><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig3-1.png?w=456" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig3-1.png?w=456?w=300" data-attachment-id="5770" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig3-1.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig3-1-4/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig3-1.png?w=456?w=676" alt="fig3-1" data-image-title="fig3-1" sizes="(max-width: 676px) 100vw, 676px" data-comments-opened="1" class="wp-image-5770" data-recalc-dims="1"></p>
</div>
<div>
<h2>1.1.e Fit a Smoothing Spline – R code</h2>
<p>Here a smoothing spline is used. Smoothing splines also do not require the explicit setting of knots. We can change the ‘degrees of freedom(df)’ paramater to get the best fit</p>
<pre><code><span># Smoothing spline has a smoothing parameter, the degrees of freedom</span>
<span># This is too wiggly</span>
<span>plot</span><span>(</span><span>auto</span><span>$</span><span>horsepower</span>,<span>auto</span><span>$</span><span>mpg</span>,<span>xlim</span><span>=</span><span>hp</span>,<span>cex</span><span>=</span><span>.5</span>,<span>col</span><span>=</span><span>&quot;black&quot;</span>,<span>xlab</span><span>=</span><span>&quot;Horsepower&quot;</span>,
     <span>ylab</span><span>=</span><span>&quot;MPG&quot;</span>, <span>main</span><span>=</span><span>&quot;Smoothing Splines&quot;</span><span>)</span>

<span># Here df is set to 16. This has a lot of variance</span>
<span>fit</span><span>=</span><span>smooth.spline</span><span>(</span><span>auto</span><span>$</span><span>horsepower</span>,<span>auto</span><span>$</span><span>mpg</span>,<span>df</span><span>=</span><span>16</span><span>)</span>
<span>lines</span><span>(</span><span>fit</span>,<span>col</span><span>=</span><span>&quot;red&quot;</span>,<span>lwd</span><span>=</span><span>2</span><span>)</span>

<span># We can use Cross Validation to allow the spline to pick the value of this smpopothing paramter. We do not need to set the degrees of freedom 'df'</span>
<span>fit</span><span>=</span><span>smooth.spline</span><span>(</span><span>auto</span><span>$</span><span>horsepower</span>,<span>auto</span><span>$</span><span>mpg</span>,<span>cv</span><span>=</span><span>TRUE</span><span>)</span>
<span>lines</span><span>(</span><span>fit</span>,<span>col</span><span>=</span><span>&quot;blue&quot;</span>,<span>lwd</span><span>=</span><span>2</span><span>)</span></code></pre>
<p><img><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig4-1.png?w=456" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig4-1.png?w=456?w=300" data-attachment-id="5773" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig4-1.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig4-1-4/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig4-1.png?w=456?w=676" alt="fig4-1" data-image-title="fig4-1" sizes="(max-width: 676px) 100vw, 676px" data-comments-opened="1" class="wp-image-5773" data-recalc-dims="1"></p>
</div>
<div>
<h2>1.1e Splines – Python</h2>
<p>There isn’t as much treatment of splines in Python and SKLearn. I did find the LSQUnivariate, UnivariateSpline spline. The LSQUnivariate spline requires the explcit setting of knots</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> scipy.interpolate <span>import</span> LSQUnivariateSpline
autoDF =pd.read_csv(<span>&quot;auto_mpg.csv&quot;</span>,encoding=<span>&quot;ISO-8859-1&quot;</span>)
autoDF.shape
autoDF.columns
autoDF1=autoDF[[<span>'mpg'</span>,<span>'cylinder'</span>,<span>'displacement'</span>,<span>'horsepower'</span>,<span>'weight'</span>,<span>'acceleration'</span>,<span>'year'</span>]]
autoDF2 = autoDF1.apply(pd.to_numeric, errors=<span>'coerce'</span>)
auto=autoDF2.dropna()
auto=auto[[<span>'horsepower'</span>,<span>'mpg'</span>]].sort_values(<span>'horsepower'</span>)

<span># Set the knots manually</span>
knots=[<span>65</span>,<span>75</span>,<span>100</span>,<span>150</span>]
<span># Create an array for X &amp; y</span>
X=np.array(auto[<span>'horsepower'</span>])
y=np.array(auto[<span>'mpg'</span>])
<span># Fit a LSQunivariate spline</span>
s = LSQUnivariateSpline(X,y,knots)

<span>#Plot the spline</span>
xs = np.linspace(<span>40</span>,<span>230</span>,<span>1000</span>)
ys = s(xs)
plt.scatter(X, y)
plt.plot(xs, ys)
plt.savefig(<span>'fig2.png'</span>, bbox_inches=<span>'tight'</span>)
</code></pre>
</div>
<div>
<h2><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig2.png?w=456" data-orig-size="546,416" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig2.png?w=456?w=300" data-attachment-id="5775" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig2.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig2-5/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig2.png?w=456?w=546" alt="fig2" data-image-title="fig2" sizes="(max-width: 546px) 100vw, 546px" data-comments-opened="1" class="wp-image-5775" data-recalc-dims="1"></h2>
<h2>1.2 Generalized Additiive models (GAMs)</h2>
<p>Generalized Additive Models (GAMs) is a really powerful ML tool.</p>
<div><img alt="y_{i} = beta_{0} + f_{1}(x_{i1}) + f_{2}(x_{i2}) + .. +f_{p}(x_{ip}) + epsilon_{i}" title="y_{i} = beta_{0} + f_{1}(x_{i1}) + f_{2}(x_{i2}) + .. +f_{p}(x_{ip}) + epsilon_{i}" src="https://s0.wp.com/latex.php?latex=y_%7Bi%7D+%3D+%5Cbeta_%7B0%7D+%2B+f_%7B1%7D%28x_%7Bi1%7D%29+%2B+f_%7B2%7D%28x_%7Bi2%7D%29+%2B+..+%2Bf_%7Bp%7D%28x_%7Bip%7D%29+%2B+%5Cepsilon_%7Bi%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></div>
<p>In GAMs we use a different functions for each of the variables. GAMs give a much better fit since we can choose any function for the different sections</p>
</div>
<div>
<h2>1.2a Generalized Additive Models (GAMs) – R Code</h2>
<p>The plot below show the smooth spline that is fit for each of the features horsepower, cylinder, displacement, year and acceleration. We can use any function for example loess, 4rd order polynomial etc.</p>
<pre><code><span>library</span><span>(</span><span>gam</span><span>)</span>
<span># Fit a smoothing spline for horsepower, cyliner, displacement and acceleration</span>
<span>gam</span><span>=</span><span>gam</span><span>(</span><span>mpg</span><span>~</span><span>s</span><span>(</span><span>horsepower</span>,<span>4</span><span>)</span><span>+</span><span>s</span><span>(</span><span>cylinder</span>,<span>5</span><span>)</span><span>+</span><span>s</span><span>(</span><span>displacement</span>,<span>4</span><span>)</span><span>+</span><span>s</span><span>(</span><span>year</span>,<span>4</span><span>)</span><span>+</span><span>s</span><span>(</span><span>acceleration</span>,<span>5</span><span>)</span>,<span>data</span><span>=</span><span>auto</span><span>)</span>
<span># Display the summary of the fit. This give the significance of each of the paramwetr</span>
<span># Also an ANOVA is given for each combination of the features</span>
<span>summary</span><span>(</span><span>gam</span><span>)</span></code></pre>
<pre><code>## 
## Call: gam(formula = mpg ~ s(horsepower, 4) + s(cylinder, 5) + s(displacement, 
##     4) + s(year, 4) + s(acceleration, 5), data = auto)
## Deviance Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.3190 -1.4436 -0.0261  1.2279 12.0873 
## 
## (Dispersion Parameter for gaussian family taken to be 6.9943)
## 
##     Null Deviance: 23818.99 on 391 degrees of freedom
## Residual Deviance: 2587.881 on 370 degrees of freedom
## AIC: 1898.282 
## 
## Number of Local Scoring Iterations: 3 
## 
## Anova for Parametric Effects
##                     Df  Sum Sq Mean Sq  F value    Pr(>F)    
## s(horsepower, 4)     1 15632.8 15632.8 2235.085 &lt; 2.2e-16 ***
## s(cylinder, 5)       1   508.2   508.2   72.666 3.958e-16 ***
## s(displacement, 4)   1   374.3   374.3   53.514 1.606e-12 ***
## s(year, 4)           1  2263.2  2263.2  323.583 &lt; 2.2e-16 ***
## s(acceleration, 5)   1   372.4   372.4   53.246 1.809e-12 ***
## Residuals          370  2587.9     7.0                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Anova for Nonparametric Effects
##                    Npar Df Npar F     Pr(F)    
## (Intercept)                                    
## s(horsepower, 4)         3 13.825 1.453e-08 ***
## s(cylinder, 5)           3 17.668 9.712e-11 ***
## s(displacement, 4)       3 44.573 &lt; 2.2e-16 ***
## s(year, 4)               3 23.364 7.183e-14 ***
## s(acceleration, 5)       4  3.848  0.004453 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</code></pre>
<pre><code><span>par</span><span>(</span><span>mfrow</span><span>=</span><span>c</span><span>(</span><span>2</span>,<span>3</span><span>)</span><span>)</span>
<span>plot</span><span>(</span><span>gam</span>,<span>se</span><span>=</span><span>TRUE</span><span>)</span></code></pre>
<p><img><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig5-1.png?w=456" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig5-1.png?w=456?w=300" data-attachment-id="5781" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig5-1.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig5-1-3/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig5-1.png?w=456?w=676" alt="fig5-1" data-image-title="fig5-1" sizes="(max-width: 676px) 100vw, 676px" data-comments-opened="1" class="wp-image-5781" data-recalc-dims="1"></p>
</div>
<div>
<h2>1.2b Generalized Additive Models (GAMs) – Python Code</h2>
<p>I did not find the equivalent of GAMs in SKlearn in Python. There was an early prototype (2012) in Github. Looks like it is still work in progress or has probably been abandoned.</p>
</div>
<div>
<h2>1.3 Tree based Machine Learning Models</h2>
<p>Tree based Machine Learning are all based on the ‘bootstrapping’ technique. In bootstrapping given a sample of size N, we create datasets of size N by sampling this original dataset with replacement. Machine Learning models are built on the different bootstrapped samples and then averaged.</p>
<p>Decision Trees as seen above have the tendency to overfit. There are several techniques that help to avoid this namely a) Bagging b) Random Forests c) Boosting</p>
<div>
<h3>Bagging, Random Forest and Gradient Boosting</h3>
<p><strong>Bagging</strong>: Bagging, or Bootstrap Aggregation decreases the variance of predictions, by creating separate Decisiion Tree based ML models on the different samples and then averaging these ML models</p>
<p><strong>Random Forests</strong>: Bagging is a greedy algorithm and tries to produce splits based on all variables which try to minimize the error. However the different ML models have a high correlation. Random Forests remove this shortcoming, by using a variable and random set of features to split on. Hence the features chosen and the resulting trees are uncorrelated. When these ML models are averaged the performance is much better.</p>
<p><strong>Boosting</strong>: Gradient Boosted Decision Trees also use an ensemble of trees but they don’t build Machine Learning models with random set of features at each step. Rather small and simple trees are built. Successive trees try to minimize the error from the earlier trees.</p>
<p><strong>Out of Bag (OOB) Error</strong>: In Random Forest and Gradient Boosting for each bootstrap sample taken from the dataset, there will be samples left out. These are known as Out of Bag samples.Classification accuracy carried out on these OOB samples is known as OOB error</p>
</div>
</div>
<div>
<h2>1.31a Decision Trees – R Code</h2>
<p>The code below creates a Decision tree with the cancer training data. The summary of the fit is output. Based on the ML model, the predict function is used on test data and a confusion matrix is output.</p>
<pre><code><span># Read the cancer data</span>
<span>library</span><span>(</span><span>tree</span><span>)</span>
<span>library</span><span>(</span><span>caret</span><span>)</span>
<span>library</span><span>(</span><span>e1071</span><span>)</span>
<span>cancer</span> <span>&lt;-</span> <span>read.csv</span><span>(</span><span>&quot;cancer.csv&quot;</span>,<span>stringsAsFactors</span> <span>=</span> <span>FALSE</span><span>)</span>
<span>cancer</span> <span>&lt;-</span> <span>cancer</span><span>[</span>,<span>2</span><span>:</span><span>32</span><span>]</span>
<span>cancer</span><span>$</span><span>target</span> <span>&lt;-</span> <span>as.factor</span><span>(</span><span>cancer</span><span>$</span><span>target</span><span>)</span>
<span>train_idx</span> <span>&lt;-</span> <span>trainTestSplit</span><span>(</span><span>cancer</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>cancer</span><span>[</span><span>train_idx</span>, <span>]</span>
<span>test</span> <span>&lt;-</span> <span>cancer</span><span>[</span><span>-</span><span>train_idx</span>, <span>]</span>

<span># Create Decision Tree</span>
<span>cancerStatus</span><span>=</span><span>tree</span><span>(</span><span>target</span><span>~</span>.,<span>train</span><span>)</span>
<span>summary</span><span>(</span><span>cancerStatus</span><span>)</span></code></pre>
<pre><code>## 
## Classification tree:
## tree(formula = target ~ ., data = train)
## Variables actually used in tree construction:
## [1] &quot;worst.perimeter&quot;      &quot;worst.concave.points&quot; &quot;area.error&quot;          
## [4] &quot;worst.texture&quot;        &quot;mean.texture&quot;         &quot;mean.concave.points&quot; 
## Number of terminal nodes:  9 
## Residual mean deviance:  0.1218 = 50.8 / 417 
## Misclassification error rate: 0.02347 = 10 / 426</code></pre>
<pre><code><span>pred</span> <span>&lt;-</span> <span>predict</span><span>(</span><span>cancerStatus</span>,<span>newdata</span><span>=</span><span>test</span>,<span>type</span><span>=</span><span>&quot;class&quot;</span><span>)</span>
<span>confusionMatrix</span><span>(</span><span>pred</span>,<span>test</span><span>$</span><span>target</span><span>)</span></code></pre>
<pre><code>## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  0  1
##          0 49  7
##          1  8 78
##                                           
##                Accuracy : 0.8944          
##                  95% CI : (0.8318, 0.9397)
##     No Information Rate : 0.5986          
##     P-Value [Acc > NIR] : 4.641e-15       
##                                           
##                   Kappa : 0.7795          
##  Mcnemar's Test P-Value : 1               
##                                           
##             Sensitivity : 0.8596          
##             Specificity : 0.9176          
##          Pos Pred Value : 0.8750          
##          Neg Pred Value : 0.9070          
##              Prevalence : 0.4014          
##          Detection Rate : 0.3451          
##    Detection Prevalence : 0.3944          
##       Balanced Accuracy : 0.8886          
##                                           
##        'Positive' Class : 0               
## </code></pre>
<pre><code><span># Plot decision tree with labels</span>
<span>plot</span><span>(</span><span>cancerStatus</span><span>)</span>
<span>text</span><span>(</span><span>cancerStatus</span>,<span>pretty</span><span>=</span><span>0</span><span>)</span></code></pre>
<p><img><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig6-1.png?w=456" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig6-1.png?w=456?w=300" data-attachment-id="5783" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig6-1.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig6-1/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig6-1.png?w=456?w=676" alt="fig6-1" data-image-title="fig6-1" sizes="(max-width: 676px) 100vw, 676px" data-comments-opened="1" class="wp-image-5783" data-recalc-dims="1"></p>
</div>
<p> </p>
<div>
<p> </p>
<h2>1.31b Decision Trees – Cross Validation – R Code</h2>
<p>We can also perform a Cross Validation on the data to identify the Decision Tree which will give the minimum deviance.</p>
<pre><code><span>library</span><span>(</span><span>tree</span><span>)</span>
<span>cancer</span> <span>&lt;-</span> <span>read.csv</span><span>(</span><span>&quot;cancer.csv&quot;</span>,<span>stringsAsFactors</span> <span>=</span> <span>FALSE</span><span>)</span>
<span>cancer</span> <span>&lt;-</span> <span>cancer</span><span>[</span>,<span>2</span><span>:</span><span>32</span><span>]</span>
<span>cancer</span><span>$</span><span>target</span> <span>&lt;-</span> <span>as.factor</span><span>(</span><span>cancer</span><span>$</span><span>target</span><span>)</span>
<span>train_idx</span> <span>&lt;-</span> <span>trainTestSplit</span><span>(</span><span>cancer</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>cancer</span><span>[</span><span>train_idx</span>, <span>]</span>
<span>test</span> <span>&lt;-</span> <span>cancer</span><span>[</span><span>-</span><span>train_idx</span>, <span>]</span>

<span># Create Decision Tree</span>
<span>cancerStatus</span><span>=</span><span>tree</span><span>(</span><span>target</span><span>~</span>.,<span>train</span><span>)</span>

<span># Execute 10 fold cross validation</span>
<span>cvCancer</span><span>=</span><span>cv.tree</span><span>(</span><span>cancerStatus</span><span>)</span>
<span>plot</span><span>(</span><span>cvCancer</span><span>)</span></code></pre>
<p><img><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig7-1.png?w=456" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig7-1.png?w=456?w=300" data-attachment-id="5785" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig7-1.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig7-1/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig7-1.png?w=456?w=676" alt="fig7-1" data-image-title="fig7-1" sizes="(max-width: 676px) 100vw, 676px" data-comments-opened="1" class="wp-image-5785" data-recalc-dims="1"></p>
<pre><code><span># Plot the </span>
<span>plot</span><span>(</span><span>cvCancer</span><span>$</span><span>size</span>,<span>cvCancer</span><span>$</span><span>dev</span>,<span>type</span><span>=</span><span>'b'</span><span>)</span></code></pre>
<p><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig13.png?w=456" data-orig-size="882,437" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig13.png?w=456?w=300" data-attachment-id="5850" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig13.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig1-10/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig13.png?w=456?w=676" alt="fig1" data-image-title="fig1" sizes="(max-width: 676px) 100vw, 676px" data-comments-opened="1" class="wp-image-5850" data-recalc-dims="1"></p>
<pre><code><span>prunedCancer</span><span>=</span><span>prune.tree</span><span>(</span><span>cancerStatus</span>,<span>best</span><span>=</span><span>4</span><span>)</span>
<span>plot</span><span>(</span><span>prunedCancer</span><span>)</span>
<span>text</span><span>(</span><span>prunedCancer</span>,<span>pretty</span><span>=</span><span>0</span><span>)</span></code></pre>
<p><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig24.png?w=456" data-orig-size="848,490" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig24.png?w=456?w=300" data-attachment-id="5851" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig24.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig2-9/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig24.png?w=456?w=676" alt="fig2" data-image-title="fig2" sizes="(max-width: 676px) 100vw, 676px" data-comments-opened="1" class="wp-image-5851" data-recalc-dims="1"></p>
<pre><code><span>pred</span> <span>&lt;-</span> <span>predict</span><span>(</span><span>prunedCancer</span>,<span>newdata</span><span>=</span><span>test</span>,<span>type</span><span>=</span><span>&quot;class&quot;</span><span>)</span>
<span>confusionMatrix</span><span>(</span><span>pred</span>,<span>test</span><span>$</span><span>target</span><span>)</span></code></pre>
<pre><code>## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  0  1
##          0 50  7
##          1  7 78
##                                          
##                Accuracy : 0.9014         
##                  95% CI : (0.8401, 0.945)
##     No Information Rate : 0.5986         
##     P-Value [Acc > NIR] : 7.988e-16      
##                                          
##                   Kappa : 0.7948         
##  Mcnemar's Test P-Value : 1              
##                                          
##             Sensitivity : 0.8772         
##             Specificity : 0.9176         
##          Pos Pred Value : 0.8772         
##          Neg Pred Value : 0.9176         
##              Prevalence : 0.4014         
##          Detection Rate : 0.3521         
##    Detection Prevalence : 0.4014         
##       Balanced Accuracy : 0.8974         
##                                          
##        'Positive' Class : 0              
## </code></pre>
</div>
<div>
<h2>1.31c Decision Trees – Python Code</h2>
<p>Below is the Python code for creating Decision Trees. The accuracy, precision, recall and F1 score is computed on the test data set.</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.metrics <span>import</span> confusion_matrix
<span>from</span> sklearn <span>import</span> tree
<span>from</span> sklearn.datasets <span>import</span> load_breast_cancer
<span>from</span> sklearn.model_selection <span>import</span> train_test_split
<span>from</span> sklearn.tree <span>import</span> DecisionTreeClassifier
<span>from</span> sklearn.datasets <span>import</span> make_classification, make_blobs
<span>from</span> sklearn.metrics <span>import</span> accuracy_score, precision_score, recall_score, f1_score
<span>import</span> graphviz 

cancer = load_breast_cancer()
(X_cancer, y_cancer) = load_breast_cancer(return_X_y = <span>True</span>)

X_train, X_test, y_train, y_test = train_test_split(X_cancer, y_cancer,
                                                   random_state = <span>0</span>)
clf = DecisionTreeClassifier().fit(X_train, y_train)

<span>print</span>(<span>'Accuracy of Decision Tree classifier on training set: {:.2f}'</span>
     .format(clf.score(X_train, y_train)))
<span>print</span>(<span>'Accuracy of Decision Tree classifier on test set: {:.2f}'</span>
     .format(clf.score(X_test, y_test)))

y_predicted=clf.predict(X_test)
confusion = confusion_matrix(y_test, y_predicted)
<span>print</span>(<span>'Accuracy: {:.2f}'</span>.format(accuracy_score(y_test, y_predicted)))
<span>print</span>(<span>'Precision: {:.2f}'</span>.format(precision_score(y_test, y_predicted)))
<span>print</span>(<span>'Recall: {:.2f}'</span>.format(recall_score(y_test, y_predicted)))
<span>print</span>(<span>'F1: {:.2f}'</span>.format(f1_score(y_test, y_predicted)))

<span># Plot the Decision Tree</span>
clf = DecisionTreeClassifier(max_depth=<span>2</span>).fit(X_train, y_train)
dot_data = tree.export_graphviz(clf, out_file=<span>None</span>, 
                         feature_names=cancer.feature_names,  
                         class_names=cancer.target_names,  
                         filled=<span>True</span>, rounded=<span>True</span>,  
                         special_characters=<span>True</span>)  
graph = graphviz.Source(dot_data)  
graph</code></pre>
<pre><code>## Accuracy of Decision Tree classifier on training set: 1.00
## Accuracy of Decision Tree classifier on test set: 0.87
## Accuracy: 0.87
## Precision: 0.97
## Recall: 0.82
## F1: 0.89</code></pre>
</div>
<div>
<h2><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/tree.png?w=456" data-orig-size="874,446" data-medium-file="https://gigadom.files.wordpress.com/2017/11/tree.png?w=456?w=300" data-attachment-id="5791" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/tree.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/tree/" data-large-file="https://gigadom.files.wordpress.com/2017/11/tree.png?w=456?w=676" alt="tree" data-image-title="tree" sizes="(max-width: 676px) 100vw, 676px" data-comments-opened="1" class="wp-image-5791" data-recalc-dims="1"></h2>
<h2>1.31d Decision Trees – Cross Validation – Python Code</h2>
<p>In the code below 5-fold cross validation is performed for different depths of the tree and the accuracy is computed. The accuracy on the test set seems to plateau when the depth is 8. But it is seen to increase again from 10 to 12. More analysis needs to be done here</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.datasets <span>import</span> load_breast_cancer
<span>from</span> sklearn.tree <span>import</span> DecisionTreeClassifier
(X_cancer, y_cancer) = load_breast_cancer(return_X_y = <span>True</span>)
<span>from</span> sklearn.cross_validation <span>import</span> train_test_split, KFold
<span><span>def</span> <span>computeCVAccuracy</span><span>(X,y,folds)</span>:</span>
    accuracy=[]
    foldAcc=[]
    depth=[<span>1</span>,<span>2</span>,<span>3</span>,<span>4</span>,<span>5</span>,<span>6</span>,<span>7</span>,<span>8</span>,<span>9</span>,<span>10</span>,<span>11</span>,<span>12</span>]
    nK=len(X)/float(folds)
    xval_err=<span>0</span>
    <span>for</span> i <span>in</span> depth: 
        kf = KFold(len(X),n_folds=folds)
        <span>for</span> train_index, test_index <span>in</span> kf:
            X_train, X_test = X.iloc[train_index], X.iloc[test_index]
            y_train, y_test = y.iloc[train_index], y.iloc[test_index]  
            clf = DecisionTreeClassifier(max_depth = i).fit(X_train, y_train)
            score=clf.score(X_test, y_test)
            accuracy.append(score)     
            
        foldAcc.append(np.mean(accuracy))  
        
    <span>return</span>(foldAcc)
    
    
cvAccuracy=computeCVAccuracy(pd.DataFrame(X_cancer),pd.DataFrame(y_cancer),folds=<span>10</span>)

df1=pd.DataFrame(cvAccuracy)
df1.columns=[<span>'cvAccuracy'</span>]
df=df1.reindex([<span>1</span>,<span>2</span>,<span>3</span>,<span>4</span>,<span>5</span>,<span>6</span>,<span>7</span>,<span>8</span>,<span>9</span>,<span>10</span>,<span>11</span>,<span>12</span>])
df.plot()
plt.title(<span>&quot;Decision Tree - 10-fold Cross Validation Accuracy vs Depth of tree&quot;</span>)
plt.xlabel(<span>&quot;Depth of tree&quot;</span>)
plt.ylabel(<span>&quot;Accuracy&quot;</span>)
plt.savefig(<span>'fig3.png'</span>, bbox_inches=<span>'tight'</span>)</code></pre>
<pre></pre>
<p> </p>
</div>
<p> </p>
<div>
<h2><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig3.png?w=456" data-orig-size="611,455" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig3.png?w=456?w=300" data-attachment-id="5794" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig3.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig3-4/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig3.png?w=456?w=611" alt="fig3" data-image-title="fig3" sizes="(max-width: 611px) 100vw, 611px" data-comments-opened="1" class="wp-image-5794" data-recalc-dims="1"></h2>
<p> </p>
<h2>1.4a Random Forest – R code</h2>
<p>A Random Forest is fit using the Boston data. The summary shows that 4 variables were randomly chosen at each split and the resulting ML model explains 88.72% of the test data. Also the variable importance is plotted. It can be seen that ‘rooms’ and ‘status’ are the most influential features in the model</p>
<pre><code><span>library</span><span>(</span><span>randomForest</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># Select specific columns</span>
<span>Boston</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;medianValue&quot;</span><span>)</span>

<span># Fit a Random Forest on the Boston training data</span>
<span>rfBoston</span><span>=</span><span>randomForest</span><span>(</span><span>medianValue</span><span>~</span>.,<span>data</span><span>=</span><span>Boston</span><span>)</span>
<span># Display the summatu of the fit. It can be seen that the MSE is 10.88 </span>
<span># and the percentage variance explained is 86.14%. About 4 variables were tried at each # #split for a maximum tree of 500.</span>
<span># The MSE and percent variance is on Out of Bag trees</span>
<span>rfBoston</span></code></pre>
<pre><code>## 
## Call:
##  randomForest(formula = medianValue ~ ., data = Boston) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 4
## 
##           Mean of squared residuals: 9.521672
##                     % Var explained: 88.72</code></pre>
<pre><code><span>#List and plot the variable importances</span>
<span>importance</span><span>(</span><span>rfBoston</span><span>)</span></code></pre>
<pre><code>##              IncNodePurity
## crimeRate        2602.1550
## zone              258.8057
## indus            2599.6635
## charles           240.2879
## nox              2748.8485
## rooms           12011.6178
## age              1083.3242
## distances        2432.8962
## highways          393.5599
## tax              1348.6987
## teacherRatio     2841.5151
## color             731.4387
## status          12735.4046</code></pre>
<pre><code><span>varImpPlot</span><span>(</span><span>rfBoston</span><span>)</span></code></pre>
<p><img><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig8-1.png?w=456" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig8-1.png?w=456?w=300" data-attachment-id="5796" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig8-1.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig8-1/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig8-1.png?w=456?w=676" alt="fig8-1" data-image-title="fig8-1" sizes="(max-width: 676px) 100vw, 676px" data-comments-opened="1" class="wp-image-5796" data-recalc-dims="1"></p>
</div>
<div>
<h2>1.4b Random Forest-OOB and Cross Validation Error – R code</h2>
<p>The figure below shows the OOB error and the Cross Validation error vs the ‘mtry’. Here mtry indicates the number of random features that are chosen at each split. The lowest test error occurs when mtry = 8</p>
<pre><code><span>library</span><span>(</span><span>randomForest</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># Select specific columns</span>
<span>Boston</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;medianValue&quot;</span><span>)</span>
<span># Split as training and tst sets</span>
<span>train_idx</span> <span>&lt;-</span> <span>trainTestSplit</span><span>(</span><span>Boston</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>Boston</span><span>[</span><span>train_idx</span>, <span>]</span>
<span>test</span> <span>&lt;-</span> <span>Boston</span><span>[</span><span>-</span><span>train_idx</span>, <span>]</span>

<span>#Initialize OOD and testError</span>
<span>oobError</span> <span>&lt;-</span> <span>NULL</span>
<span>testError</span> <span>&lt;-</span> <span>NULL</span>
<span># In the code below the number of variables to consider at each split is increased</span>
<span># from 1 - 13(max features) and the OOB error and the MSE is computed</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>fitRF</span><span>=</span><span>randomForest</span><span>(</span><span>medianValue</span><span>~</span>.,<span>data</span><span>=</span><span>train</span>,<span>mtry</span><span>=</span><span>i</span>,<span>ntree</span><span>=</span><span>400</span><span>)</span>
    <span>oobError</span><span>[</span><span>i</span><span>]</span> <span>&lt;-</span><span>fitRF</span><span>$</span><span>mse</span><span>[</span><span>400</span><span>]</span>
    <span>pred</span> <span>&lt;-</span> <span>predict</span><span>(</span><span>fitRF</span>,<span>newdata</span><span>=</span><span>test</span><span>)</span>
    <span>testError</span><span>[</span><span>i</span><span>]</span> <span>&lt;-</span> <span>mean</span><span>(</span><span>(</span><span>pred</span><span>-</span><span>test</span><span>$</span><span>medianValue</span><span>)</span><span>^</span><span>2</span><span>)</span>
<span>}</span>

<span># We can see the OOB and Test Error. It can be seen that the Random Forest performs</span>
<span># best with the lowers MSE at mtry=6</span>
<span>matplot</span><span>(</span><span>1</span><span>:</span><span>13</span>,<span>cbind</span><span>(</span><span>testError</span>,<span>oobError</span><span>)</span>,<span>pch</span><span>=</span><span>19</span>,<span>col</span><span>=</span><span>c</span><span>(</span><span>&quot;red&quot;</span>,<span>&quot;blue&quot;</span><span>)</span>,
        <span>type</span><span>=</span><span>&quot;b&quot;</span>,<span>xlab</span><span>=</span><span>&quot;mtry(no of varaibles at each split)&quot;</span>, <span>ylab</span><span>=</span><span>&quot;Mean Squared Error&quot;</span>,
        <span>main</span><span>=</span><span>&quot;Random Forest - OOB and Test Error&quot;</span><span>)</span>
<span>legend</span><span>(</span><span>&quot;topright&quot;</span>,<span>legend</span><span>=</span><span>c</span><span>(</span><span>&quot;OOB&quot;</span>,<span>&quot;Test&quot;</span><span>)</span>,<span>pch</span><span>=</span><span>19</span>,<span>col</span><span>=</span><span>c</span><span>(</span><span>&quot;red&quot;</span>,<span>&quot;blue&quot;</span><span>)</span><span>)</span></code></pre>
<p><img><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig9-1.png?w=456" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig9-1.png?w=456?w=300" data-attachment-id="5798" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig9-1.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig9-1/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig9-1.png?w=456?w=676" alt="fig9-1" data-image-title="fig9-1" sizes="(max-width: 676px) 100vw, 676px" data-comments-opened="1" class="wp-image-5798" data-recalc-dims="1"></p>
</div>
<div>
<h2>1.4c Random Forest – Python code</h2>
<p>The python code for Random Forest Regression is shown below. The training and test score is computed. The variable importance shows that ‘rooms’ and ‘status’ are the most influential of the variables</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.ensemble <span>import</span> RandomForestRegressor
df = pd.read_csv(<span>&quot;Boston.csv&quot;</span>,encoding = <span>&quot;ISO-8859-1&quot;</span>)

X=df[[<span>'crimeRate'</span>,<span>'zone'</span>, <span>'indus'</span>,<span>'charles'</span>,<span>'nox'</span>,<span>'rooms'</span>, <span>'age'</span>,<span>'distances'</span>,<span>'highways'</span>,<span>'tax'</span>,
       <span>'teacherRatio'</span>,<span>'color'</span>,<span>'status'</span>]]
y=df[<span>'medianValue'</span>]

X_train, X_test, y_train, y_test = train_test_split(X, y, random_state = <span>0</span>)

regr = RandomForestRegressor(max_depth=<span>4</span>, random_state=<span>0</span>)
regr.fit(X_train, y_train)

<span>print</span>(<span>'R-squared score (training): {:.3f}'</span>
     .format(regr.score(X_train, y_train)))
<span>print</span>(<span>'R-squared score (test): {:.3f}'</span>
     .format(regr.score(X_test, y_test)))

feature_names=[<span>'crimeRate'</span>,<span>'zone'</span>, <span>'indus'</span>,<span>'charles'</span>,<span>'nox'</span>,<span>'rooms'</span>, <span>'age'</span>,<span>'distances'</span>,<span>'highways'</span>,<span>'tax'</span>,
       <span>'teacherRatio'</span>,<span>'color'</span>,<span>'status'</span>]
<span>print</span>(regr.feature_importances_)
plt.figure(figsize=(<span>10</span>,<span>6</span>),dpi=<span>80</span>)
c_features=X_train.shape[<span>1</span>]
plt.barh(np.arange(c_features),regr.feature_importances_)
plt.xlabel(<span>&quot;Feature importance&quot;</span>)
plt.ylabel(<span>&quot;Feature name&quot;</span>)

plt.yticks(np.arange(c_features), feature_names)
plt.tight_layout()

plt.savefig(<span>'fig4.png'</span>, bbox_inches=<span>'tight'</span>)
</code></pre>
<pre><code>## R-squared score (training): 0.917
## R-squared score (test): 0.734
## [ 0.03437382  0.          0.00580335  0.          0.00731004  0.36461548
##   0.00638577  0.03432173  0.0041244   0.01732328  0.01074148  0.0012638
##   0.51373683]</code></pre>
</div>
<div>
<h2><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig4.png?w=456" data-orig-size="791,472" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig4.png?w=456?w=300" data-attachment-id="5801" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig4.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig4-4/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig4.png?w=456?w=676" alt="fig4" data-image-title="fig4" sizes="(max-width: 676px) 100vw, 676px" data-comments-opened="1" class="wp-image-5801" data-recalc-dims="1"></h2>
<h2>1.4d Random Forest – Cross Validation and OOB Error – Python code</h2>
<p>As with R the ‘max_features’ determines the random number of features the random forest will use at each split. The plot shows that when max_features=8 the MSE is lowest</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.ensemble <span>import</span> RandomForestRegressor
<span>from</span> sklearn.model_selection <span>import</span> cross_val_score
df = pd.read_csv(<span>&quot;Boston.csv&quot;</span>,encoding = <span>&quot;ISO-8859-1&quot;</span>)

X=df[[<span>'crimeRate'</span>,<span>'zone'</span>, <span>'indus'</span>,<span>'charles'</span>,<span>'nox'</span>,<span>'rooms'</span>, <span>'age'</span>,<span>'distances'</span>,<span>'highways'</span>,<span>'tax'</span>,
       <span>'teacherRatio'</span>,<span>'color'</span>,<span>'status'</span>]]
y=df[<span>'medianValue'</span>]

cvError=[]
oobError=[]
oobMSE=[]
<span>for</span> i <span>in</span> range(<span>1</span>,<span>13</span>):
    regr = RandomForestRegressor(max_depth=<span>4</span>, n_estimators=<span>400</span>,max_features=i,oob_score=True,random_state=<span>0</span>)
    mse= np.mean(cross_val_score(regr, X, y, cv=<span>5</span>,scoring = <span>'neg_mean_squared_error'</span>))
    <span># Since this is neg_mean_squared_error I have inverted the sign to get MSE</span>
    cvError.append(-mse)
    # Fit on all data to compute OOB error
    regr.fit(X, y)
    # Record the OOB error for each `max_features=i` setting
    oob = 1 - regr.oob_score_
    oobError.append(oob)
    # Get the Out of Bag prediction
    oobPred=regr.oob_prediction_ 
    # Compute the Mean Squared Error between OOB Prediction and target
    mseOOB=np.mean(np.square(oobPred-y))
    oobMSE.append(mseOOB)

# Plot the CV Error and OOB Error
# Set max_features
maxFeatures=np.arange(1,13) 
cvError=pd.DataFrame(cvError,index=maxFeatures)
oobMSE=pd.DataFrame(oobMSE,index=maxFeatures)
#Plot
fig8=df.plot()
fig8=plt.title('Random forest - CV Error and OOB Error vs max_features')
fig8.figure.savefig('fig8.png', bbox_inches='tight')</code><code>
</code><code>
#Plot the OOB Error vs max_features
plt.plot(range(1,13),oobError)
fig2=plt.title(&quot;Random Forest - OOB Error vs max_features (variable no of features)&quot;)
fig2=plt.xlabel(&quot;max_features (variable no of features)&quot;)
fig2=plt.ylabel(&quot;OOB Error&quot;)
fig2.figure.savefig('fig7.png', bbox_inches='tight')
</code></pre>
<pre><code></code></pre>
</div>
<div></div>
<h2><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig8.png?w=456" data-orig-size="552,435" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig8.png?w=456?w=300" data-attachment-id="5834" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig8.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig8-3/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig8.png?w=456?w=552" alt="fig8" data-image-title="fig8" sizes="(max-width: 552px) 100vw, 552px" data-comments-opened="1" class="wp-image-5834" data-recalc-dims="1"> <img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig73.png?w=456" data-orig-size="625,455" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig73.png?w=456?w=300" data-attachment-id="5831" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig73.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig7-7/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig73.png?w=456?w=625" alt="fig7" data-image-title="fig7" sizes="(max-width: 625px) 100vw, 625px" data-comments-opened="1" class="wp-image-5831" data-recalc-dims="1"></h2>
<div>
<h2>1.5a Boosting – R code</h2>
<p>Here a Gradient Boosted ML Model is built with a n.trees=5000, with a learning rate of 0.01 and depth of 4. The feature importance plot also shows that rooms and status are the 2 most important features. The MSE vs the number of trees plateaus around 2000 trees</p>
<pre><code><span>library</span><span>(</span><span>gbm</span><span>)</span>
<span># Perform gradient boosting on the Boston data set. The distribution is gaussian since we</span>
<span># doing MSE. The interaction depth specifies the number of splits</span>
<span>boostBoston</span><span>=</span><span>gbm</span><span>(</span><span>medianValue</span><span>~</span>.,<span>data</span><span>=</span><span>train</span>,<span>distribution</span><span>=</span><span>&quot;gaussian&quot;</span>,<span>n.trees</span><span>=</span><span>5000</span>,
                <span>shrinkage</span><span>=</span><span>0.01</span>,<span>interaction.depth</span><span>=</span><span>4</span><span>)</span>
<span>#The summary gives the variable importance. The 2 most significant variables are</span>
<span># number of rooms and lower status</span>
<span>summary</span><span>(</span><span>boostBoston</span><span>)</span></code></pre>
<p><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig33.png" data-orig-size="947,542" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig33.png?w=300" data-attachment-id="5853" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig33.png?w=450&h=387&fit=676%2C387" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig3-7/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig33.png?w=450&amp;h=387&amp;fit=676%2C387" alt="" data-image-title="fig3" sizes="(max-width: 676px) 100vw, 676px" data-comments-opened="1" class="wp-image-5853" data-recalc-dims="1"></p>
<pre><code>##                       var    rel.inf
## rooms               rooms 42.2267200
## status             status 27.3024671
## distances       distances  7.9447972
## crimeRate       crimeRate  5.0238827
## nox                   nox  4.0616548
## teacherRatio teacherRatio  3.1991999
## age                   age  2.7909772
## color               color  2.3436295
## tax                   tax  2.1386213
## charles           charles  1.3799109
## highways         highways  0.7644026
## indus               indus  0.7236082
## zone                 zone  0.1001287</code></pre>
<pre><code><span># The plots below show how each variable relates to the median value of the home. As</span>
<span># the number of roomd increase the median value increases and with increase in lower status</span>
<span># the median value decreases</span>
<span>par</span><span>(</span><span>mfrow</span><span>=</span><span>c</span><span>(</span><span>1</span>,<span>2</span><span>)</span><span>)</span>
<span>#Plot the relation between the top 2 features and the target</span>
<span>plot</span><span>(</span><span>boostBoston</span>,<span>i</span><span>=</span><span>&quot;rooms&quot;</span><span>)</span>
<span>plot</span><span>(</span><span>boostBoston</span>,<span>i</span><span>=</span><span>&quot;status&quot;</span><span>)</span></code></pre>
<p><img><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig10-2.png?w=456" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig10-2.png?w=456?w=300" data-attachment-id="5806" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig10-2.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig10-2/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig10-2.png?w=456?w=676" alt="fig10-2" data-image-title="fig10-2" sizes="(max-width: 676px) 100vw, 676px" data-comments-opened="1" class="wp-image-5806" data-recalc-dims="1"></p>
<pre><code><span># Create a sequence of trees between 100-5000 incremented by 50</span>
<span>nTrees</span><span>=</span><span>seq</span><span>(</span><span>100</span>,<span>5000</span>,<span>by</span><span>=</span><span>50</span><span>)</span>
<span># Predict the values for the test data</span>
<span>pred</span> <span>&lt;-</span> <span>predict</span><span>(</span><span>boostBoston</span>,<span>newdata</span><span>=</span><span>test</span>,<span>n.trees</span><span>=</span><span>nTrees</span><span>)</span>
<span># Compute the mean for each of the MSE for each of the number of trees </span>
<span>boostError</span> <span>&lt;-</span> <span>apply</span><span>(</span><span>(</span><span>pred</span><span>-</span><span>test</span><span>$</span><span>medianValue</span><span>)</span><span>^</span><span>2</span>,<span>2</span>,<span>mean</span><span>)</span>
<span>#Plot the MSE vs the number of trees</span>
<span>plot</span><span>(</span><span>nTrees</span>,<span>boostError</span>,<span>pch</span><span>=</span><span>19</span>,<span>col</span><span>=</span><span>&quot;blue&quot;</span>,<span>ylab</span><span>=</span><span>&quot;Mean Squared Error&quot;</span>,
     <span>main</span><span>=</span><span>&quot;Boosting Test Error&quot;</span><span>)</span></code></pre>
<p><img><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig10-3.png?w=456" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig10-3.png?w=456?w=300" data-attachment-id="5810" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig10-3.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig10-3/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig10-3.png?w=456?w=676" alt="fig10-3" data-image-title="fig10-3" sizes="(max-width: 676px) 100vw, 676px" data-comments-opened="1" class="wp-image-5810" data-recalc-dims="1"></p>
</div>
<div>
<h2>1.5b Cross Validation Boosting – R code</h2>
<p>Included below is a cross validation error vs the learning rate. The lowest error is when learning rate = 0.09</p>
<pre><code><span>cvError</span> <span>&lt;-</span> <span>NULL</span>
<span>s</span> <span>&lt;-</span> <span>c</span><span>(</span><span>.001</span>,<span>0.01</span>,<span>0.03</span>,<span>0.05</span>,<span>0.07</span>,<span>0.09</span>,<span>0.1</span><span>)</span>
<span>for</span><span>(</span><span>i</span> <span>in</span> <span>seq_along</span><span>(</span><span>s</span><span>)</span><span>)</span><span>{</span>
    <span>cvBoost</span><span>=</span><span>gbm</span><span>(</span><span>medianValue</span><span>~</span>.,<span>data</span><span>=</span><span>train</span>,<span>distribution</span><span>=</span><span>&quot;gaussian&quot;</span>,<span>n.trees</span><span>=</span><span>5000</span>,
                <span>shrinkage</span><span>=</span><span>s</span><span>[</span><span>i</span><span>]</span>,<span>interaction.depth</span><span>=</span><span>4</span>,<span>cv.folds</span><span>=</span><span>5</span><span>)</span>
    <span>cvError</span><span>[</span><span>i</span><span>]</span> <span>&lt;-</span> <span>mean</span><span>(</span><span>cvBoost</span><span>$</span><span>cv.error</span><span>)</span>
<span>}</span>

<span># Create a data frame for plotting</span>
<span>a</span> <span>&lt;-</span> <span>rbind</span><span>(</span><span>s</span>,<span>cvError</span><span>)</span>
<span>b</span> <span>&lt;-</span> <span>as.data.frame</span><span>(</span><span>t</span><span>(</span><span>a</span><span>)</span><span>)</span>
<span># It can be seen that a shrinkage parameter of 0,05 gives the lowes CV Error</span>
<span>ggplot</span><span>(</span><span>b</span>,<span>aes</span><span>(</span><span>s</span>,<span>cvError</span><span>)</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;Shrinkage&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;Gradient boosted trees - Cross Validation error vs Shrinkage&quot;</span><span>)</span></code></pre>
<p><img><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig11-1.png?w=456" data-orig-size="1344,960" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig11-1.png?w=456?w=300" data-attachment-id="5812" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig11-1.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig11-1/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig11-1.png?w=456?w=676" alt="fig11-1" data-image-title="fig11-1" sizes="(max-width: 676px) 100vw, 676px" data-comments-opened="1" class="wp-image-5812" data-recalc-dims="1"></p>
</div>
<div>
<h2>1.5c Boosting – Python code</h2>
<p>A gradient boost ML model in Python is created below. The Rsquared score is computed on the training and test data.</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.ensemble <span>import</span> GradientBoostingRegressor
df = pd.read_csv(<span>&quot;Boston.csv&quot;</span>,encoding = <span>&quot;ISO-8859-1&quot;</span>)

X=df[[<span>'crimeRate'</span>,<span>'zone'</span>, <span>'indus'</span>,<span>'charles'</span>,<span>'nox'</span>,<span>'rooms'</span>, <span>'age'</span>,<span>'distances'</span>,<span>'highways'</span>,<span>'tax'</span>,
       <span>'teacherRatio'</span>,<span>'color'</span>,<span>'status'</span>]]
y=df[<span>'medianValue'</span>]

X_train, X_test, y_train, y_test = train_test_split(X, y, random_state = <span>0</span>)

regr = GradientBoostingRegressor()
regr.fit(X_train, y_train)

<span>print</span>(<span>'R-squared score (training): {:.3f}'</span>
     .format(regr.score(X_train, y_train)))
<span>print</span>(<span>'R-squared score (test): {:.3f}'</span>
     .format(regr.score(X_test, y_test)))</code></pre>
<pre><code>## R-squared score (training): 0.983
## R-squared score (test): 0.821</code></pre>
</div>
<div>
<h2>1.5c Cross Validation Boosting – Python code</h2>
<p>the cross validation error is computed as the learning rate is varied. The minimum CV eror occurs when lr = 0.04</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.ensemble <span>import</span> RandomForestRegressor
<span>from</span> sklearn.ensemble <span>import</span> GradientBoostingRegressor
<span>from</span> sklearn.model_selection <span>import</span> cross_val_score
df = pd.read_csv(<span>&quot;Boston.csv&quot;</span>,encoding = <span>&quot;ISO-8859-1&quot;</span>)

X=df[[<span>'crimeRate'</span>,<span>'zone'</span>, <span>'indus'</span>,<span>'charles'</span>,<span>'nox'</span>,<span>'rooms'</span>, <span>'age'</span>,<span>'distances'</span>,<span>'highways'</span>,<span>'tax'</span>,
       <span>'teacherRatio'</span>,<span>'color'</span>,<span>'status'</span>]]
y=df[<span>'medianValue'</span>]

cvError=[]
learning_rate =[.<span>001</span>,<span>0.01</span>,<span>0.03</span>,<span>0.05</span>,<span>0.07</span>,<span>0.09</span>,<span>0.1</span>]
<span>for</span> lr <span>in</span> learning_rate:
    regr = GradientBoostingRegressor(max_depth=<span>4</span>, n_estimators=<span>400</span>,learning_rate  =lr,random_state=<span>0</span>)
    mse= np.mean(cross_val_score(regr, X, y, cv=<span>10</span>,scoring = <span>'neg_mean_squared_error'</span>))
    <span># Since this is neg_mean_squared_error I have inverted the sign to get MSE</span>
    cvError.append(-mse)
learning_rate =[.<span>001</span>,<span>0.01</span>,<span>0.03</span>,<span>0.05</span>,<span>0.07</span>,<span>0.09</span>,<span>0.1</span>]
plt.plot(learning_rate,cvError)
plt.title(<span>&quot;Gradient Boosting - 5-fold CV- Mean Squared Error vs max_features (variable no of features)&quot;</span>)
plt.xlabel(<span>&quot;max_features (variable no of features)&quot;</span>)
plt.ylabel(<span>&quot;Mean Squared Error&quot;</span>)
plt.savefig(<span>'fig6.png'</span>, bbox_inches=<span>'tight'</span>)</code></pre>
<p><img data-orig-file="https://gigadom.files.wordpress.com/2017/11/fig6.png?w=456" data-orig-size="788,455" data-medium-file="https://gigadom.files.wordpress.com/2017/11/fig6.png?w=456?w=300" data-attachment-id="5815" data-image-meta="{" data-image-description="" src="https://gigadom.files.wordpress.com/2017/11/fig6.png?w=456" data-permalink="https://gigadom.wordpress.com/2017/11/07/practical-machine-learning-with-r-and-python-part-5/fig6-4/" data-large-file="https://gigadom.files.wordpress.com/2017/11/fig6.png?w=456?w=676" alt="fig6" data-image-title="fig6" sizes="(max-width: 676px) 100vw, 676px" data-comments-opened="1" class="wp-image-5815" data-recalc-dims="1"></p>
<p><strong>Conclusion</strong> This post covered Splines and Tree based ML models like Bagging, Random Forest and Boosting. Stay tuned for further updates.</p>
<p>You may also like</p>
<ol>
<li><a rel="nofollow" href="https://gigadom.wordpress.com/2016/05/14/re-introducing-cricketr-an-r-package-to-analyze-performances-of-cricketers/" target="_blank">Re-introducing cricketr! : An R package to analyze performances of cricketer</a></li>
<li><a rel="nofollow" href="https://gigadom.wordpress.com/2013/02/14/designing-a-social-web-portal/" target="_blank">Designing a Social Web Portal</a></li>
<li><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></li>
<li><a rel="nofollow" href="https://gigadom.wordpress.com/2013/06/20/my-tedx-talk-on-the-internet-of-things/" target="_blank">My TEDx talk on the “Internet of Things”</a></li>
<li><a rel="nofollow" href="https://gigadom.wordpress.com/2011/06/10/singularity/" target="_blank">Singularity</a></li>
<li><a rel="nofollow" href="https://gigadom.wordpress.com/2013/07/08/a-closer-look-at-robot-horse-on-a-canter-in-android/" target="_blank">A closer look at “Robot Horse on a Trot” in Android</a></li>
</ol>
<p>To see all posts see <a rel="nofollow" href="https://gigadom.wordpress.com/aa-2/" target="_blank">Index of posts</a></p>
</div>
<p> <a rel="nofollow" href="https://feeds.wordpress.com/1.0/gocomments/gigadom.wordpress.com/5753/" target="_blank"><img border="0" alt="" src="https://feeds.wordpress.com/1.0/comments/gigadom.wordpress.com/5753/"></a> <img border="0" alt="" src="https://i1.wp.com/pixel.wp.com/b.gif?resize=1%2C1&ssl=1" data-recalc-dims="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/11/07/practical-machine-learning-with-r-and-python-part-5/"> 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/276J2A36UQo" height="1"></a></div>