<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/which-linear-model-is-best/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://realdataweb.wordpress.com/2017/05/10/which-linear-model-is-best/"> R – Real Data</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><img sizes="(max-width: 626px) 100vw, 626px" data-image-description="" data-orig-file="https://realdataweb.files.wordpress.com/2017/05/trains-railroad_00389817.jpg?w=456" data-orig-size="626,358" data-image-title="trains-railroad_00389817" data-comments-opened="1" class="wp-image-552" data-recalc-dims="1" data-attachment-id="552" data-permalink="https://realdataweb.wordpress.com/2017/05/10/which-linear-model-is-best/trains-railroad_00389817/" alt="trains-railroad_00389817" data-medium-file="https://realdataweb.files.wordpress.com/2017/05/trains-railroad_00389817.jpg?w=456?w=300" data-large-file="https://realdataweb.files.wordpress.com/2017/05/trains-railroad_00389817.jpg?w=456?w=626" data-image-meta="{" src="https://realdataweb.files.wordpress.com/2017/05/trains-railroad_00389817.jpg?w=456"></p>
<p>Recently I have been working on <a target="_blank" rel="nofollow" href="https://www.kaggle.com/c/sberbank-russian-housing-market">a Kaggle competition</a> where participants are tasked with predicting Russian housing prices. In developing a model for the challenge, I came across a few methods for selecting the best regression<br>
model for a given dataset. Let’s load up some data and take a look.</p>
<pre>
library(ISLR)
set.seed(123)
swiss &lt;- data.frame(datasets::swiss)
dim(swiss)
</pre>
<pre><code>## [1] 47  6
</code></pre>
<pre>
head(swiss)
</pre>
<pre><code>##              Fertility Agriculture Examination Education Catholic
## Courtelary        80.2        17.0          15        12     9.96
## Delemont          83.1        45.1           6         9    84.84
## Franches-Mnt      92.5        39.7           5         5    93.40
## Moutier           85.8        36.5          12         7    33.77
## Neuveville        76.9        43.5          17        15     5.16
## Porrentruy        76.1        35.3           9         7    90.57
##              Infant.Mortality
## Courtelary               22.2
## Delemont                 22.2
## Franches-Mnt             20.2
## Moutier                  20.3
## Neuveville               20.6
## Porrentruy               26.6
</code></pre>
<p>This data provides data on Swiss fertility and some socioeconomic<br>
indicators. Suppose we want to predict fertility rate. Each observation<br>
is as a percentage. In order to assess our prediction ability, create<br>
test and training data sets.</p>
<pre>
test &lt;- swiss[sample(nrow(swiss),5),]
test &lt;- data.frame(test)

train &lt;- swiss[-(sample(nrow(swiss), 5)),]
train &lt;- data.frame(train)
</pre>
<p>Next, have a quick look at the data.</p>
<pre>
par(mfrow=c(2,2))
plot(train$Education, train$Fertility)
plot(train$Catholic, train$Fertility)
plot(train$Infant.Mortality, train$Fertility)
hist(train$Fertility)
</pre>
<p><img sizes="(max-width: 600px) 100vw, 600px" data-image-description="" data-orig-file="https://realdataweb.files.wordpress.com/2017/05/swissplots.png?w=456" data-orig-size="600,600" data-image-title="swissplots" data-comments-opened="1" class="wp-image-563" data-recalc-dims="1" data-attachment-id="563" data-permalink="https://realdataweb.wordpress.com/2017/05/10/which-linear-model-is-best/swissplots/" alt="swissplots" data-medium-file="https://realdataweb.files.wordpress.com/2017/05/swissplots.png?w=456?w=300" data-large-file="https://realdataweb.files.wordpress.com/2017/05/swissplots.png?w=456?w=600" data-image-meta="{" src="https://realdataweb.files.wordpress.com/2017/05/swissplots.png?w=456"></p>
<p>Looks like there could be some interesting relationships here. The<br>
following block of code will take a model formula (or model matrix) and<br>
search for the best combination of predictors.</p>
<pre>
library(leaps)
leap &lt;- regsubsets(Fertility~., data = train, nbest = 10)
leapsum &lt;- summary(leap)
plot(leap, scale = 'adjr2')
</pre>
<p><img alt="leap plot" src="https://realdataweb.files.wordpress.com/2017/05/leap3.png?w=456" data-recalc-dims="1"></p>
<p>According to the adjusted R-squared value (larger is better), the best<br>
two models are: those with all predictors and all predictors less<br>
Examination. Both have adjusted R-squared values around .64 – a decent fit. Fit these models so we can evaluate them further.</p>
<pre>
swisslm &lt;- lm(Fertility~., data = train)
swisslm2 &lt;- lm(Fertility~.-Examination, data = train)
#use summary() for a more detailed look.
</pre>
<p>First we'll compute the mean squared error (MSE).</p>
<pre>
mean(swisslm$residuals^2)
</pre>
<pre><code>## [1] 45.44987
</code></pre>
<pre>
mean(swisslm2$residuals^2)
</pre>
<pre><code>## [1] 47.16062
</code></pre>
<p>Looks like the model with all predictors is marginally better. We're<br>
looking for a low value of MSE here. We can also look at the Akaike<br>
information criterion (AIC) for more information. Lower is better here<br>
as well.</p>
<pre>
library(MASS)
step1 &lt;- stepAIC(swisslm, direction = &quot;both&quot;)
</pre>
<pre>
step1$anova
</pre>
<pre><code>## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## Fertility ~ Agriculture + Examination + Education + Catholic + 
##     Infant.Mortality
## 
## Final Model:
## Fertility ~ Education + Catholic + Infant.Mortality
## 
## 
##            Step Df Deviance Resid. Df Resid. Dev      AIC
## 1                                  36   1908.895 172.2976
## 2 - Examination  1 71.85140        37   1980.746 171.8495
## 3 - Agriculture  1 78.11025        38   2058.856 171.4739
</code></pre>
<p>Here, the model without <code>Examination</code> scores lower than the full<br>
model. It seems that both models are evenly matched, though I might be<br>
inclined to use the full model since the MSE is lower.</p>
<p>Since we sampled our original data to make the train and test datasets,<br>
the difference in these tests is subject to change based on the training<br>
data used. I encourage anyone who wants to test themselves to change the<br>
<code>set.seed</code> at the beginning and evaluate the above results again. Even<br>
better: use your own data!</p>
<p>If you learned something or have questions feel free to leave a comment<br>
or shoot me an email.</p>
<p>Happy modeling,</p>
<p>Kiefer</p>
<p> <a target="_blank" rel="nofollow" href="http://feeds.wordpress.com/1.0/gocomments/realdataweb.wordpress.com/550/"><img alt="" src="http://feeds.wordpress.com/1.0/comments/realdataweb.wordpress.com/550/" border="0"></a> <img alt="" src="https://i1.wp.com/pixel.wp.com/b.gif?resize=1%2C1&ssl=1" border="0" 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://realdataweb.wordpress.com/2017/05/10/which-linear-model-is-best/"> R – Real Data</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/NdtFSDPSjpc"></a></div>