<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="http://www.r-bloggers.com/visualizing-bootrapped-stepwise-regression-in-r-using-plotly/"></div></div><div><a data-url="http://www.r-bloggers.com/visualizing-bootrapped-stepwise-regression-in-r-using-plotly/" data-counturl="http://www.r-bloggers.com/visualizing-bootrapped-stepwise-regression-in-r-using-plotly/" data-count="horizontal" data-text="Visualizing Bootrapped Stepwise Regression in R using Plotly" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://moderndata.plot.ly/visualizing-bootrapped-stepwise-regression-in-r-using-plotly/"> R – Modern Data</a></strong>, and kindly contributed to <a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>We all have used stepwise regression at some point. Stepwise regression is known to be sensitive to initial inputs. One way to mitigate this sensitivity is to repeatedly run stepwise regression on bootstrap samples.</p>
<p>R has a nice package called <strong>bootStepAIC()</strong> which (from its description) “<em>Implements a Bootstrap procedure to investigate the variability of model selection under the stepAIC() stepwise algorithm of package MASS.</em>”</p>
<p>It provides a lot of information as an output and sometimes it can get challenging to keep track of all of this information especially if there are a lot of covariates. In this post we’ll try to come up with a simple visualization aimed at summarizing the output from the function <code>boot.stepAIC()</code>.</p>
<div>
<h3>Running <code>boot.stepAIC()</code></h3>
<p>Using the <code>boot.stepAIC()</code> is fairly simple. Just input an already fitted <code>lm/glm</code> model and th associated dataset.</p>
<p>We’ll use the <code>BostonHousing</code> dataset from the <code>mlbench</code> package. More details <a target="_blank" href="https://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.names">here</a></p>

<pre>library(bootStepAIC)
library(plotly)
library(mlbench)

# Load Boston housing dataset
data(&quot;BostonHousing&quot;)

# Fit Linear regression model
fit &lt;- lm(crim ~ ., data = BostonHousing)

# Run bootstrapped stepwise regression
fit.boot &lt;- boot.stepAIC(fit, data = BostonHousing, B = 100) # That's it !</pre>

</div>
<div>
<h3>Collecting required information</h3>
<p>The output from <code>boot.stepAIC()</code> contains the following. Note that each output is shown as a percentage (based on the total number of bootstrapped samples)</p>
<ul>
<li>No of times a covariate was featured in the final model from <code>stepAIC()</code></li>
<li>No of times a covariate’s coefficient sign was positive / negative</li>
<li>No of times a covariate was statistically significant (default at alpha = 5%)</li>
</ul>
<p>We’ll collect all of this information first and create data frames so as to make charting easier later on.</p>
<p>Note that in this particualr example there is a variable by the name <strong>chas</strong> which is a factor with levels 0 and 1. R renames the variable as <strong>chas1</strong> by default.</p>

<pre># Extract data
nBoot &lt;- summary(fit.boot)[8,1]
origModel &lt;- paste(names(coef(fit.boot$OrigModel)), collapse = &quot; + &quot;)
stepModel &lt;- paste(names(coef(fit.boot$OrigStepAIC)), collapse = &quot; + &quot;)

# Names of covariates
covariates &lt;- rownames(fit.boot$Covariates)
nCovariates &lt;- length(covariates)

# Matrix of number of times each covariate was picked
coef.pick &lt;- fit.boot$Covariates

# Matrix for the consistency of sign on each covariate
coef.sign &lt;- fit.boot$Sign

# Change name for &quot;chas&quot; since it is a factor
rownames(coef.sign)[7] &lt;- &quot;chas&quot;
coef.sign &lt;- coef.sign[match(rownames(coef.pick), rownames(coef.sign)),]

# Matrix for statistical significance
coef.stat &lt;- fit.boot$Significance

# Change name for &quot;chas&quot; since it is a factor
rownames(coef.stat)[11] &lt;- &quot;chas&quot;
coef.stat &lt;- coef.stat[match(rownames(coef.pick), rownames(coef.stat)),]

# Make into long form for charting later
coef.stat.long &lt;- data.frame()

for(i in 1:length(coef.stat)){
  n &lt;- round(coef.stat[i],0)
  vec &lt;- seq(0, n, by = 2)
  mat &lt;- data.frame(rep(names(coef.stat)[i], length(vec)), vec, paste(&quot;% Sig&quot;, n))
  names(mat) &lt;- c(&quot;variable&quot;, &quot;sig&quot;, &quot;text&quot;)
  
  # We'll use mode = &quot;line&quot;. NA helps separate line segments
  coef.stat.long &lt;- rbind(coef.stat.long, mat, c(NA, NA))
}

# Convert to dataframes
coef.pick &lt;- as.data.frame(coef.pick)
coef.stat &lt;- as.data.frame(coef.stat)
coef.sign &lt;- as.data.frame(coef.sign)

names(coef.pick) &lt;- &quot;pick&quot;
names(coef.sign) &lt;- c(&quot;pos&quot;, &quot;neg&quot;)
names(coef.stat) &lt;- &quot;stat&quot;</pre>

</div>
<div>
<h3>Plot</h3>
<p>Now that we have all the information we need, we just need to plot. The plot is arranged as such:</p>
<ul>
<li>One layer for the number of times a variable was picked up by <code>stepAIC()</code> (barplot)</li>
<li>One layer for the positive and negative coefficients (scatter plot using triangles)</li>
<li>One layer for the number of times a variable was significant (vertical line chart)</li>
<li>Annotation for some other information</li>
</ul>
<pre># Base plot for number of times a variable was picked by stepAIC
plot_ly(coef.pick, x = rownames(coef.pick), y = pick,
        type = &quot;bar&quot;, opacity = 0.75, name = &quot;Times picked (%)&quot;,
        hoverinfo = &quot;text&quot;, text = pick.text,
        marker = list(color = &quot;#00994d&quot;, line = list(width = 2))) %>% 
  
  # Layer for number of times a variable was statistically significant at 5%
  add_trace(data = coef.stat.long, x = variable, y = sig, 
            type = &quot;scatter&quot;, mode = &quot;markers + line&quot;, name = &quot;Stat. Sig (%)&quot;,
            line = list(color = &quot;#ffdb4d&quot;, width = 15),
            hoverinfo = &quot;text&quot;, text = text) %>% 
  
  # Layer for number of times a variable's coefficient was positive
  add_trace(data = coef.sign, x = rownames(coef.pick), y = rep(-5, nCovariates), 
            type = &quot;scatter&quot;, mode = &quot;markers&quot;, name = &quot;Coef Sign(% pos)&quot;,
            marker = list(symbol = &quot;triangle-up&quot;, size = pos/scale, color = &quot;#4da6ff&quot;,
                          line = list(color = &quot;black&quot;, width = 2)),
            hoverinfo = &quot;text&quot;, text = sign.text.up) %>% 
  
  # Layer for number of times a variable's coefficient was negative
  add_trace(data = coef.sign, x = rownames(coef.pick), y = rep(-10, nCovariates), 
            type = &quot;scatter&quot;, mode = &quot;markers&quot;, name = &quot;Coef Sign(% neg)&quot;,
            marker = list(symbol = &quot;triangle-down&quot;, size = neg/scale, color = &quot;#ff704d&quot;,
                          line = list(color = &quot;black&quot;, width = 2)),
            hoverinfo = &quot;text&quot;, text = sign.text.down) %>% 
  
  # Layout, annotations, axis options etc
  layout(xaxis = list(title = &quot;&lt;b>Covariates&lt;/b>&quot;),
         yaxis = list(title = &quot;&lt;b>Percentage(%)&lt;/b>&quot;,
                      tickmode = &quot;array&quot;, 
                      tickvals = round(seq(0, 100, length.out = 10), 0),
                      domain = c(0.2, 1)),
         plot_bgcolor = &quot;#e1efc3&quot;,
         paper_bgcolor = &quot;#e1efc3&quot;,
         
         annotations = list(
           list(x = 0.1, y = 1, 
                xref = &quot;paper&quot;, yref = &quot;paper&quot;, 
                xanchor = &quot;left&quot;, yanchor = &quot;top&quot;,
                ax = 0, ay = 0,
                text = &quot;Visualizing &lt;em>boot.stepAIC()&lt;/em>&quot;,
                font = list(family = &quot;serif&quot;, size = 30)),
           
           list(x = 0.3, y = 0.1, 
                xref = &quot;paper&quot;, yref = &quot;paper&quot;, 
                xanchor = &quot;left&quot;, yanchor = &quot;top&quot;,
                ax = 0, ay = 0,
                text = paste(&quot;&lt;em>Original Model:&lt;/em>&quot;, origModel),
                font = list(family = &quot;PT Sans Narrow&quot;, size = 15)),
           
           list(x = 0.21, y = 0.05, 
                xref = &quot;paper&quot;, yref = &quot;paper&quot;, 
                xanchor = &quot;left&quot;, yanchor = &quot;top&quot;, align = &quot;left&quot;,
                ax = 0, ay = 0,
                text = paste(&quot;&lt;em>Stepwise Model:&lt;/em>&quot;, stepModel),
                font = list(family = &quot;PT Sans Narrow&quot;, size = 15)),
           
           list(x = 0.8, y = 0.90, 
                xref = &quot;paper&quot;, yref = &quot;paper&quot;, 
                xanchor = &quot;left&quot;, yanchor = &quot;top&quot;, align = &quot;left&quot;,
                ax = 0, ay = 0,
                text = paste0(&quot;&lt;em>No. of Covariates:&lt;/em>&quot;, nCovariates, &quot;&lt;br>&quot;,
                              &quot;&lt;em>No. of bootstrap samples:&lt;/em>&quot;, nBoot, &quot;&lt;br>&quot;),
                font = list(family = &quot;PT Sans Narrow&quot;, size = 15))
         ))</pre>
<p><iframe height="600" width="1024" src="https://plot.ly/~riddhiman/292.embed"></iframe></p>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="http://www.r-bloggers.com/visualizing-bootrapped-stepwise-regression-in-r-using-plotly/"></div></div><div><a data-url="http://www.r-bloggers.com/visualizing-bootrapped-stepwise-regression-in-r-using-plotly/" data-counturl="http://www.r-bloggers.com/visualizing-bootrapped-stepwise-regression-in-r-using-plotly/" data-count="vertical" data-text="Visualizing Bootrapped Stepwise Regression in R using Plotly" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://moderndata.plot.ly/visualizing-bootrapped-stepwise-regression-in-r-using-plotly/"> R – Modern Data</a></strong>.</div>
<hr>
<a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="http://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="http://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="http://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="http://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="http://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="http://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="http://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="http://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="http://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="http://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="http://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="http://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="http://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="http://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="http://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="http://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/pXv6e_MmE9k"></a></div></div>