<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/rstanarm-basics-visualizing-uncertainty-in-linear-regression/"></div></div><div><a data-url="https://www.r-bloggers.com/rstanarm-basics-visualizing-uncertainty-in-linear-regression/" data-counturl="https://www.r-bloggers.com/rstanarm-basics-visualizing-uncertainty-in-linear-regression/" data-count="horizontal" data-text="RStanARM basics: visualizing uncertainty in linear regression" 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="https://tjmahr.github.io/visualizing-uncertainty-rstanarm/"> Higher Order Functions</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>As part of <a target="_blank" rel="nofollow" href="https://tjmahr.github.io/rstanarm-tutorial-slides/">my tutorial talk on RStanARM</a>, I<br>
presented some examples of how to visualize the uncertainty in Bayesian linear<br>
regression models. This post is an expanded demonstration of the approaches I<br>
presented in that tutorial.</p>
<h2>Data: Does brain mass predict how much mammals sleep in a day?</h2>
<p>Let’s use the <a target="_blank" rel="nofollow" href="https://vincentarelbundock.github.io/Rdatasets/doc/ggplot2/msleep.html">mammal sleep dataset from ggplot2</a>. This dataset<br>
contains the number of hours spent sleeping per day for 83 different species of<br>
mammals along with each species’ brain mass (kg) and body mass (kg), among other<br>
measures. Here’s a first look at the data.</p>
<div>
<pre><code><span>library</span><span>(</span><span>dplyr</span><span>,</span><span> </span><span>warn.conflicts</span><span> </span><span>=</span><span> </span><span>FALSE</span><span>)</span><span>
</span><span>library</span><span>(</span><span>ggplot2</span><span>)</span><span>
</span><span># Preview sorted by brain/body ratio. I chose this sorting so that humans would
# show up in the preview.
</span><span>msleep</span><span> </span><span>%>%</span><span> </span><span>select</span><span>(</span><span>name</span><span>,</span><span> </span><span>sleep_total</span><span>,</span><span> </span><span>brainwt</span><span>,</span><span> </span><span>bodywt</span><span>,</span><span> </span><span>everything</span><span>())</span><span> </span><span>%>%</span><span> </span><span>arrange</span><span>(</span><span>desc</span><span>(</span><span>brainwt</span><span> </span><span>/</span><span> </span><span>bodywt</span><span>))</span><span>
</span><span>#> # A tibble: 83 × 11
#> name sleep_total brainwt bodywt genus
#> &lt;chr> &lt;dbl> &lt;dbl> &lt;dbl> &lt;chr>
#> 1 Thirteen-lined ground squirrel 13.8 0.00400 0.101 Spermophilus
#> 2 Owl monkey 17.0 0.01550 0.480 Aotus
#> 3 Lesser short-tailed shrew 9.1 0.00014 0.005 Cryptotis
#> 4 Squirrel monkey 9.6 0.02000 0.743 Saimiri
#> 5 Macaque 10.1 0.17900 6.800 Macaca
#> 6 Little brown bat 19.9 0.00025 0.010 Myotis
#> 7 Galago 9.8 0.00500 0.200 Galago
#> 8 Mole rat 10.6 0.00300 0.122 Spalax
#> 9 Tree shrew 8.9 0.00250 0.104 Tupaia
#> 10 Human 8.0 1.32000 62.000 Homo
#> # ... with 73 more rows, and 6 more variables: vore &lt;chr>, order &lt;chr>,
#> # conservation &lt;chr>, sleep_rem &lt;dbl>, sleep_cycle &lt;dbl>, awake &lt;dbl>
</span><span>
</span><span>ggplot</span><span>(</span><span>msleep</span><span>)</span><span> </span><span>+</span><span> </span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>brainwt</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>sleep_total</span><span>)</span><span> </span><span>+</span><span> </span><span>geom_point</span><span>()</span><span>
</span><span>#> Warning: Removed 27 rows containing missing values (geom_point).
</span></code></pre>
</div>
<p><img title="Brain mass by sleep hours. This plot looks terrible because the masses span many orders of magnitudes." alt="Brain mass by sleep hours. This plot looks terrible because the masses span many orders of magnitudes." src="https://i1.wp.com/tjmahr.github.io/figs//2016-11-18-visualizing-uncertainty-rstanarm/brain-sleep-1.png?w=80%25&ssl=1" data-recalc-dims="1"></p>
<p>Hmmm, not very helpful! We should put our measures on a log-10 scale. Also, 27<br>
of the species don’t have brain mass data, so we’ll exclude those rows for the<br>
rest of this tutorial.</p>
<div>
<pre><code><span>msleep</span><span> </span><span>&lt;-</span><span> </span><span>msleep</span><span> </span><span>%>%</span><span> </span><span>filter</span><span>(</span><span>!</span><span>is.na</span><span>(</span><span>brainwt</span><span>))</span><span> </span><span>%>%</span><span> </span><span>mutate</span><span>(</span><span>log_brainwt</span><span> </span><span>=</span><span> </span><span>log10</span><span>(</span><span>brainwt</span><span>),</span><span> </span><span>log_bodywt</span><span> </span><span>=</span><span> </span><span>log10</span><span>(</span><span>bodywt</span><span>),</span><span> </span><span>log_sleep_total</span><span> </span><span>=</span><span> </span><span>log10</span><span>(</span><span>sleep_total</span><span>))</span><span>
</span></code></pre>
</div>
<p>Now, plot the log-transformed data. But let’s also get a little fancy and label<br>
the points for some example critters :cat: so that we can get some intuition<br>
about the data in this scaling. (Plus, I wanted to try out the <a target="_blank" rel="nofollow" href="https://tjmahr.github.io/visualizing-uncertainty-rstanarm/r4ds.had.co.nz/graphics-for-communication.html#annotations">annotation<br>
tips</a> from the <em>R4DS</em> book.)</p>
<div>
<pre><code><span># Create a separate data-frame of species to highlight
</span><span>ex_mammals</span><span> </span><span>&lt;-</span><span> </span><span>c</span><span>(</span><span>&quot;Domestic cat&quot;</span><span>,</span><span> </span><span>&quot;Human&quot;</span><span>,</span><span> </span><span>&quot;Dog&quot;</span><span>,</span><span> </span><span>&quot;Cow&quot;</span><span>,</span><span> </span><span>&quot;Rabbit&quot;</span><span>,</span><span>
</span><span>&quot;Big brown bat&quot;</span><span>,</span><span> </span><span>&quot;House mouse&quot;</span><span>,</span><span> </span><span>&quot;Horse&quot;</span><span>,</span><span> </span><span>&quot;Golden hamster&quot;</span><span>)</span><span>
</span><span># We will give some familiar species shorter names
</span><span>renaming_rules</span><span> </span><span>&lt;-</span><span> </span><span>c</span><span>(</span><span>
</span><span>&quot;Domestic cat&quot;</span><span> </span><span>=</span><span> </span><span>&quot;Cat&quot;</span><span>,</span><span> </span><span>&quot;Golden hamster&quot;</span><span> </span><span>=</span><span> </span><span>&quot;Hamster&quot;</span><span>,</span><span> </span><span>&quot;House mouse&quot;</span><span> </span><span>=</span><span> </span><span>&quot;Mouse&quot;</span><span>)</span><span>
</span><span>ex_points</span><span> </span><span>&lt;-</span><span> </span><span>msleep</span><span> </span><span>%>%</span><span> </span><span>filter</span><span>(</span><span>name</span><span> </span><span>%in%</span><span> </span><span>ex_mammals</span><span>)</span><span> </span><span>%>%</span><span> </span><span>mutate</span><span>(</span><span>name</span><span> </span><span>=</span><span> </span><span>stringr</span><span>::</span><span>str_replace_all</span><span>(</span><span>name</span><span>,</span><span> </span><span>renaming_rules</span><span>))</span><span>
</span><span># Define these labels only once for all the plots
</span><span>lab_lines</span><span> </span><span>&lt;-</span><span> </span><span>list</span><span>(</span><span>
</span><span>brain_log</span><span> </span><span>=</span><span> </span><span>&quot;Brain mass (kg., log-scaled)&quot;</span><span>,</span><span> </span><span>sleep_raw</span><span> </span><span>=</span><span> </span><span>&quot;Sleep per day (hours)&quot;</span><span>,</span><span>
</span><span>sleep_log</span><span> </span><span>=</span><span> </span><span>&quot;Sleep per day (log-hours)&quot;</span><span>
</span><span>)</span><span>
</span><span>ggplot</span><span>(</span><span>msleep</span><span>)</span><span> </span><span>+</span><span> </span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>brainwt</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>sleep_total</span><span>)</span><span> </span><span>+</span><span> </span><span>geom_point</span><span>(</span><span>color</span><span> </span><span>=</span><span> </span><span>&quot;grey40&quot;</span><span>)</span><span> </span><span>+</span><span>
</span><span># Circles around highlighted points + labels
</span><span> </span><span>geom_point</span><span>(</span><span>size</span><span> </span><span>=</span><span> </span><span>3</span><span>,</span><span> </span><span>shape</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span> </span><span>color</span><span> </span><span>=</span><span> </span><span>&quot;grey40&quot;</span><span>,</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>ex_points</span><span>)</span><span> </span><span>+</span><span>
</span><span>ggrepel</span><span>::</span><span>geom_text_repel</span><span>(</span><span>aes</span><span>(</span><span>label</span><span> </span><span>=</span><span> </span><span>name</span><span>),</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>ex_points</span><span>)</span><span> </span><span>+</span><span> </span><span># Use log scaling on x-axis
</span><span> </span><span>scale_x_log10</span><span>(</span><span>breaks</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>.001</span><span>,</span><span> </span><span>.01</span><span>,</span><span> </span><span>.1</span><span>,</span><span> </span><span>1</span><span>))</span><span> </span><span>+</span><span> </span><span>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>lab_lines</span><span>$</span><span>brain_log</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>lab_lines</span><span>$</span><span>sleep_raw</span><span>)</span><span>
</span></code></pre>
</div>
<p><img title="Brain mass by sleep hours, now with both on a log-10 scale. Some species have their data highlighted." alt="Brain mass by sleep hours, now with both on a log-10 scale. Some species have their data highlighted." src="https://i2.wp.com/tjmahr.github.io/figs//2016-11-18-visualizing-uncertainty-rstanarm/log-brain-sleep-1.png?w=80%25&ssl=1" data-recalc-dims="1"></p>
<p>As a child growing up on a dairy farm :cow:, it was remarkable to me how little<br>
I saw cows sleeping, compared to dogs or cats. Were they okay? Are they<br>
constantly tired and groggy? Maybe they are asleep when I’m asleep? Here, it<br>
looks like they just don’t need very much sleep.</p>
<p>Next, let’s fit a classical regression model. We will use a log-scaled sleep<br>
measure so that the regression line doesn’t imply negative sleep (even though<br>
brains never get <em>that</em> large).</p>
<div>
<pre><code><span>m</span><span>1</span><span>_</span><span>classical</span><span> </span><span>&lt;-</span><span> </span><span>lm</span><span>(</span><span>log_sleep_total</span><span> </span><span>~</span><span> </span><span>log_brainwt</span><span>,</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>msleep</span><span>)</span><span> </span><span>arm</span><span>::</span><span>display</span><span>(</span><span>m</span><span>1</span><span>_</span><span>classical</span><span>)</span><span>
</span><span>#> lm(formula = log_sleep_total ~ log_brainwt, data = msleep)
#> coef.est coef.se
#> (Intercept) 0.74 0.04 #> log_brainwt -0.13 0.02 #> ---
#> n = 56, k = 2
#> residual sd = 0.17, R-Squared = 0.40
</span></code></pre>
</div>
<p>We can interpret the model in the usual way: A mammal with 1 kg (0 log-kg)<br>
of brain mass sleeps 10<sup>0.74</sup> = 5.5 hours per<br>
day. A mammal with a tenth of that brain mass (-1 log-kg) sleeps<br>
10<sup>0.74 + 0.13</sup> = 7.4 hours.</p>
<p>We illustrate the regression results to show the predicted mean of <em>y</em> and<br>
its 95% confidence interval. This task is readily accomplished in ggplot2 using<br>
<code>stat_smooth()</code>. This function fits a model and plots the mean and CI for each<br>
aesthetic grouping of data<sup><a target="_blank" rel="nofollow" href="https://tjmahr.github.io/visualizing-uncertainty-rstanarm/#fn:1">1</a></sup> in a plot.</p>
<div>
<pre><code><span>ggplot</span><span>(</span><span>msleep</span><span>)</span><span> </span><span>+</span><span> </span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>log_brainwt</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>log_sleep_total</span><span>)</span><span> </span><span>+</span><span> </span><span>geom_point</span><span>()</span><span> </span><span>+</span><span>
</span><span>stat_smooth</span><span>(</span><span>method</span><span> </span><span>=</span><span> </span><span>&quot;lm&quot;</span><span>,</span><span> </span><span>level</span><span> </span><span>=</span><span> </span><span>.95</span><span>)</span><span> </span><span>+</span><span> </span><span>scale_x_continuous</span><span>(</span><span>labels</span><span> </span><span>=</span><span> </span><span>function</span><span>(</span><span>x</span><span>)</span><span> </span><span>10</span><span> </span><span>^</span><span> </span><span>x</span><span>)</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>lab_lines</span><span>$</span><span>brain_log</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>lab_lines</span><span>$</span><span>sleep_log</span><span>)</span><span>
</span></code></pre>
</div>
<p><img title="Brain mass by sleep hours, log-10 scale, plus the predicted mean and 95% CI from a linear regression." alt="Brain mass by sleep hours, log-10 scale, plus the predicted mean and 95% CI from a linear regression." src="https://i2.wp.com/tjmahr.github.io/figs//2016-11-18-visualizing-uncertainty-rstanarm/log-brain-sleep-lm-fit-1.png?w=80%25&ssl=1" data-recalc-dims="1"></p>
<p>This interval conveys some uncertainty in the estimate of the mean, but this<br>
interval has a <a target="_blank" rel="nofollow" href="https://stats.stackexchange.com/questions/2272/whats-the-difference-between-a-confidence-interval-and-a-credible-interval">frequentist interpretation</a> which can be<br>
unintuitive for this sort of data.</p>
<p>Now, for the point of this post: <strong>What’s the Bayesian version of this kind of<br>
visualization</strong>? Specifically, we want to illustrate:</p>
<ul>
<li>Predictions from a regression model</li>
<li>Some uncertainty about those predictions</li>
<li>Raw data used to train the model</li>
</ul>
<h2>Option 1: The pile-of-lines plot</h2>
<p>The regression line in the classical plot is just one particular line. It’s the<br>
line of best fit that satisfies a least-squares or maximum-likelihood objective.<br>
Our Bayesian model estimates an entire distribution of plausible<br>
regression lines. The first way to visualize our uncertainty is to plot our<br>
own line of best fit along with a sample of other lines from the posterior<br>
distribution of the model.</p>
<p>First, we fit a model RStanARM using weakly informative priors.</p>
<div>
<pre><code><span>library</span><span>(</span><span>&quot;rstanarm&quot;</span><span>)</span><span>
</span><span>m</span><span>1</span><span> </span><span>&lt;-</span><span> </span><span>stan_glm</span><span>(</span><span>
</span><span>log_sleep_total</span><span> </span><span>~</span><span> </span><span>log_brainwt</span><span>,</span><span> </span><span>family</span><span> </span><span>=</span><span> </span><span>gaussian</span><span>(),</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>msleep</span><span>,</span><span> </span><span>prior</span><span> </span><span>=</span><span> </span><span>normal</span><span>(</span><span>0</span><span>,</span><span> </span><span>3</span><span>),</span><span>
</span><span>prior_intercept</span><span> </span><span>=</span><span> </span><span>normal</span><span>(</span><span>0</span><span>,</span><span> </span><span>3</span><span>))</span><span>
</span></code></pre>
</div>
<p>We now have 4,000 credible regressions lines for our data.</p>
<div>
<pre><code><span>summary</span><span>(</span><span>m</span><span>1</span><span>)</span><span>
</span><span>#> stan_glm(formula = log_sleep_total ~ log_brainwt, family = gaussian(), #> data = msleep, prior = normal(0, 3), prior_intercept = normal(0, #> 3))
#> #> Family: gaussian (identity)
#> Algorithm: sampling
#> Posterior sample size: 4000
#> Observations: 56
#> #> Estimates:
#> mean sd 2.5% 25% 50% 75% 97.5%
#> (Intercept) 0.7 0.0 0.6 0.7 0.7 0.8 0.8 #> log_brainwt -0.1 0.0 -0.2 -0.1 -0.1 -0.1 -0.1 #> sigma 0.2 0.0 0.1 0.2 0.2 0.2 0.2 #> mean_PPD 1.0 0.0 0.9 0.9 1.0 1.0 1.0 #> log-posterior 12.0 1.2 9.0 11.5 12.3 12.9 13.4 #> #> Diagnostics:
#> mcse Rhat n_eff
#> (Intercept) 0.0 1.0 3040 #> log_brainwt 0.0 1.0 3046 #> sigma 0.0 1.0 2862 #> mean_PPD 0.0 1.0 3671 #> log-posterior 0.0 1.0 2159 #> #> For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
</span></code></pre>
</div>
<p>For models fit by RStanARM, the generic coefficient function <code>coef()</code> returns<br>
the median parameter values.</p>
<div>
<pre><code><span>coef</span><span>(</span><span>m</span><span>1</span><span>)</span><span>
</span><span>#> (Intercept) log_brainwt #> 0.7354829 -0.1263922
</span><span>coef</span><span>(</span><span>m</span><span>1</span><span>_</span><span>classical</span><span>)</span><span>
</span><span>#> (Intercept) log_brainwt #> 0.7363492 -0.1264049
</span></code></pre>
</div>
<p>We can see that the intercept and slope of the median line is pretty close to<br>
the classical model’s intercept and slope. The median line serves as the “point<br>
estimate” for our model: If we had to summarize the modeled relationship using<br>
just a single number for each parameter, we can use the medians.</p>
<p>One way to visualize our model therefore is to <strong>plot our point-estimate line<br>
plus a sample of the other credible lines from our model</strong>. First, we create a<br>
data-frame with all 4,000 regression lines.</p>
<div>
<pre><code><span># Coercing a model to a data-frame returns data-frame of posterior samples. # One row per sample.
</span><span>fits</span><span> </span><span>&lt;-</span><span> </span><span>m</span><span>1</span><span> </span><span>%>%</span><span> </span><span>as_data_frame</span><span> </span><span>%>%</span><span> </span><span>rename</span><span>(</span><span>intercept</span><span> </span><span>=</span><span> </span><span>`(Intercept)`</span><span>)</span><span> </span><span>%>%</span><span> </span><span>select</span><span>(</span><span>-</span><span>sigma</span><span>)</span><span>
</span><span>fits</span><span>
</span><span>#> # A tibble: 4,000 × 2
#> intercept log_brainwt
#> &lt;dbl> &lt;dbl>
#> 1 0.7529824 -0.1369554
#> 2 0.7243708 -0.1266290
#> 3 0.7575502 -0.1171410
#> 4 0.7855554 -0.1031353
#> 5 0.6327073 -0.1795992
#> 6 0.6474521 -0.1714347
#> 7 0.7512467 -0.1155559
#> 8 0.7363273 -0.1162038
#> 9 0.7490401 -0.1276618
#> 10 0.7238091 -0.1305896
#> # ... with 3,990 more rows
</span></code></pre>
</div>
<p>We now plot the 500 randomly sampled lines from our model with light,<br>
semi-transparent lines.</p>
<div>
<pre><code><span># aesthetic controllers
</span><span>n_draws</span><span> </span><span>&lt;-</span><span> </span><span>500</span><span>
</span><span>alpha_level</span><span> </span><span>&lt;-</span><span> </span><span>.15</span><span>
</span><span>col_draw</span><span> </span><span>&lt;-</span><span> </span><span>&quot;grey60&quot;</span><span>
</span><span>col_median</span><span> </span><span>&lt;-</span><span> </span><span>&quot;#3366FF&quot;</span><span>
</span><span>ggplot</span><span>(</span><span>msleep</span><span>)</span><span> </span><span>+</span><span> </span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>log_brainwt</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>log_sleep_total</span><span>)</span><span> </span><span>+</span><span> </span><span># Plot a random sample of rows as gray semi-transparent lines
</span><span> </span><span>geom_abline</span><span>(</span><span>aes</span><span>(</span><span>intercept</span><span> </span><span>=</span><span> </span><span>intercept</span><span>,</span><span> </span><span>slope</span><span> </span><span>=</span><span> </span><span>log_brainwt</span><span>),</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>sample_n</span><span>(</span><span>fits</span><span>,</span><span> </span><span>n_draws</span><span>),</span><span> </span><span>color</span><span> </span><span>=</span><span> </span><span>col_draw</span><span>,</span><span> </span><span>alpha</span><span> </span><span>=</span><span> </span><span>alpha_level</span><span>)</span><span> </span><span>+</span><span> </span><span># Plot the median values in blue
</span><span> </span><span>geom_abline</span><span>(</span><span>intercept</span><span> </span><span>=</span><span> </span><span>median</span><span>(</span><span>fits</span><span>$</span><span>intercept</span><span>),</span><span> </span><span>slope</span><span> </span><span>=</span><span> </span><span>median</span><span>(</span><span>fits</span><span>$</span><span>log_brainwt</span><span>),</span><span> </span><span>size</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span> </span><span>color</span><span> </span><span>=</span><span> </span><span>col_median</span><span>)</span><span> </span><span>+</span><span>
</span><span>geom_point</span><span>()</span><span> </span><span>+</span><span> </span><span>scale_x_continuous</span><span>(</span><span>labels</span><span> </span><span>=</span><span> </span><span>function</span><span>(</span><span>x</span><span>)</span><span> </span><span>10</span><span> </span><span>^</span><span> </span><span>x</span><span>)</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>lab_lines</span><span>$</span><span>brain_log</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>lab_lines</span><span>$</span><span>sleep_log</span><span>)</span><span>
</span></code></pre>
</div>
<p><img title="Brain mass by sleep hours, log-10 scale, plus the median regression line and 500 random regressions lines sampled from the posterior." alt="Brain mass by sleep hours, log-10 scale, plus the median regression line and 500 random regressions lines sampled from the posterior." src="https://i2.wp.com/tjmahr.github.io/figs//2016-11-18-visualizing-uncertainty-rstanarm/pile-of-lines-plot-1.png?w=80%25&ssl=1" data-recalc-dims="1"></p>
<p>Each of these light lines represents a credible prediction of the mean across<br>
the values of <em>x</em>. As these line pile up on top of each other, they create an<br>
uncertainty band around our line of best fit. More plausible lines are more<br>
likely to be sampled, so these lines overlap and create a uniform color around<br>
the median line. As we move left or right, getting farther away from the mean of<br>
<em>x</em>, the lines start to fan out and we see very faint individual lines for some<br>
of the more extreme (yet still plausible) lines.</p>
<p>The advantage of this plot is that it is a direct visualization of posterior<br>
samples—one line per sample. It provides an estimate for the central tendency<br>
in the data but it also converys uncertainty around that estimate.</p>
<p>This approach has limitations, however. Lines for subgroups require a little<br>
more effort to undo interactions. Also, the regression lines span the whole <em>x</em><br>
axis which is not appropriate when subgroups only use a portion of the x-axis.<br>
(This limitation is solvable though.) Finally, I haven’t found good defaults<br>
for the aesthetic options: The number of samples, the colors to use, and the<br>
transparency level. One can lose of lots and lots and lots time fiddling with<br>
those knobs!</p>
<h2>Option 2: Mean and its 95% interval</h2>
<p>Another option is a direct port of the <code>stat_smooth()</code> plot: Draw a line of<br>
best fit and the 95% uncertainty interval around it.</p>
<p>To limit the amount of the <em>x</em> axis used by the lines, we’re going to create a<br>
sequence of 80 points along the range of the data.</p>
<div>
<pre><code><span>x_rng</span><span> </span><span>&lt;-</span><span> </span><span>range</span><span>(</span><span>msleep</span><span>$</span><span>log_brainwt</span><span>)</span><span> </span><span>x_steps</span><span> </span><span>&lt;-</span><span> </span><span>seq</span><span>(</span><span>x_rng</span><span>[</span><span>1</span><span>],</span><span> </span><span>x_rng</span><span>[</span><span>2</span><span>],</span><span> </span><span>length.out</span><span> </span><span>=</span><span> </span><span>80</span><span>)</span><span>
</span><span>new_data</span><span> </span><span>&lt;-</span><span> </span><span>data_frame</span><span>(</span><span>
</span><span>observation</span><span> </span><span>=</span><span> </span><span>seq_along</span><span>(</span><span>x_steps</span><span>),</span><span> </span><span>log_brainwt</span><span> </span><span>=</span><span> </span><span>x_steps</span><span>)</span><span>
</span><span>new_data</span><span>
</span><span>#> # A tibble: 80 × 2
#> observation log_brainwt
#> &lt;int> &lt;dbl>
#> 1 1 -3.853872
#> 2 2 -3.795509
#> 3 3 -3.737146
#> 4 4 -3.678784
#> 5 5 -3.620421
#> 6 6 -3.562058
#> 7 7 -3.503695
#> 8 8 -3.445332
#> 9 9 -3.386970
#> 10 10 -3.328607
#> # ... with 70 more rows
</span></code></pre>
</div>
<p>The function <code>posterior_linpred()</code> returns the model-fitted means for a data-frame<br>
of new data. I say <em>means</em> because the function computes 80 predicted means for<br>
each sample from the posterior. The result is 4000 x 80 matrix of fitted means.</p>
<div>
<pre><code><span>pred_lin</span><span> </span><span>&lt;-</span><span> </span><span>posterior_linpred</span><span>(</span><span>m</span><span>1</span><span>,</span><span> </span><span>newdata</span><span> </span><span>=</span><span> </span><span>new_data</span><span>)</span><span>
</span><span>dim</span><span>(</span><span>pred_lin</span><span>)</span><span>
</span><span>#> [1] 4000 80
</span></code></pre>
</div>
<p>We are going to reduce this down to just a median and 95% interval around each<br>
point. I do some tidying to get the data into a long format (one row per fitted<br>
mean per posterior sample), and then do a table-join with the <code>observation</code><br>
column included in <code>new_data</code>. I store these steps in a function because I<br>
have to do them again later in this post.</p>
<div>
<pre><code><span>tidy_predictions</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>mat_pred</span><span>,</span><span> </span><span>df_data</span><span>,</span><span> </span><span>obs_name</span><span> </span><span>=</span><span> </span><span>&quot;observation&quot;</span><span>,</span><span>
</span><span>prob_lwr</span><span> </span><span>=</span><span> </span><span>.025</span><span>,</span><span> </span><span>prob_upr</span><span> </span><span>=</span><span> </span><span>.975</span><span>)</span><span> </span><span>{</span><span>
</span><span># Get data-frame with one row per fitted value per posterior sample
</span><span> </span><span>df_pred</span><span> </span><span>&lt;-</span><span> </span><span>mat_pred</span><span> </span><span>%>%</span><span> </span><span>as_data_frame</span><span> </span><span>%>%</span><span> </span><span>setNames</span><span>(</span><span>seq_len</span><span>(</span><span>ncol</span><span>(</span><span>.</span><span>)))</span><span> </span><span>%>%</span><span> </span><span>tibble</span><span>::</span><span>rownames_to_column</span><span>(</span><span>&quot;posterior_sample&quot;</span><span>)</span><span> </span><span>%>%</span><span> </span><span>tidyr</span><span>::</span><span>gather_</span><span>(</span><span>obs_name</span><span>,</span><span> </span><span>&quot;fitted&quot;</span><span>,</span><span> </span><span>setdiff</span><span>(</span><span>names</span><span>(</span><span>.</span><span>),</span><span> </span><span>&quot;posterior_sample&quot;</span><span>))</span><span>
</span><span>df_pred</span><span>
</span><span># Helps with joining later
</span><span> </span><span>class</span><span>(</span><span>df_pred</span><span>[[</span><span>obs_name</span><span>]])</span><span> </span><span>&lt;-</span><span> </span><span>class</span><span>(</span><span>df_data</span><span>[[</span><span>obs_name</span><span>]])</span><span>
</span><span># Summarise prediction interval for each observation
</span><span> </span><span>df_pred</span><span> </span><span>%>%</span><span> </span><span>group_by_</span><span>(</span><span>obs_name</span><span>)</span><span> </span><span>%>%</span><span> </span><span>summarise</span><span>(</span><span>median</span><span> </span><span>=</span><span> </span><span>median</span><span>(</span><span>fitted</span><span>),</span><span>
</span><span>lower</span><span> </span><span>=</span><span> </span><span>quantile</span><span>(</span><span>fitted</span><span>,</span><span> </span><span>prob_lwr</span><span>),</span><span> </span><span>upper</span><span> </span><span>=</span><span> </span><span>quantile</span><span>(</span><span>fitted</span><span>,</span><span> </span><span>prob_upr</span><span>))</span><span> </span><span>%>%</span><span> </span><span>left_join</span><span>(</span><span>df_data</span><span>,</span><span> </span><span>by</span><span> </span><span>=</span><span> </span><span>obs_name</span><span>)</span><span>
</span><span>}</span><span>
</span><span>df_pred_lin</span><span> </span><span>&lt;-</span><span> </span><span>tidy_predictions</span><span>(</span><span>pred_lin</span><span>,</span><span> </span><span>new_data</span><span>)</span><span>
</span><span>df_pred_lin</span><span>
</span><span>#> # A tibble: 80 × 5
#> observation median lower upper log_brainwt
#> &lt;int> &lt;dbl> &lt;dbl> &lt;dbl> &lt;dbl>
#> 1 1 1.223770 1.128224 1.320591 -3.853872
#> 2 2 1.216516 1.122147 1.311214 -3.795509
#> 3 3 1.209222 1.117190 1.301462 -3.737146
#> 4 4 1.201831 1.112268 1.291821 -3.678784
#> 5 5 1.194506 1.107512 1.282047 -3.620421
#> 6 6 1.187240 1.102580 1.272930 -3.562058
#> 7 7 1.179955 1.096945 1.263415 -3.503695
#> 8 8 1.172608 1.091237 1.254113 -3.445332
#> 9 9 1.165268 1.085800 1.244733 -3.386970
#> 10 10 1.157932 1.080823 1.235356 -3.328607
#> # ... with 70 more rows
</span></code></pre>
</div>
<p>We can do the line-plus-interval plot using <code>geom_ribbon()</code> for the uncertainty<br>
band.</p>
<div>
<pre><code><span>p_linpread</span><span> </span><span>&lt;-</span><span> </span><span>ggplot</span><span>(</span><span>msleep</span><span>)</span><span> </span><span>+</span><span> </span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>log_brainwt</span><span>)</span><span> </span><span>+</span><span> </span><span>geom_ribbon</span><span>(</span><span>aes</span><span>(</span><span>ymin</span><span> </span><span>=</span><span> </span><span>lower</span><span>,</span><span> </span><span>ymax</span><span> </span><span>=</span><span> </span><span>upper</span><span>),</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>df_pred_lin</span><span>,</span><span> </span><span>alpha</span><span> </span><span>=</span><span> </span><span>0.4</span><span>,</span><span> </span><span>fill</span><span> </span><span>=</span><span> </span><span>&quot;grey60&quot;</span><span>)</span><span> </span><span>+</span><span> </span><span>geom_line</span><span>(</span><span>aes</span><span>(</span><span>y</span><span> </span><span>=</span><span> </span><span>median</span><span>),</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>df_pred_lin</span><span>,</span><span> </span><span>colour</span><span> </span><span>=</span><span> </span><span>&quot;#3366FF&quot;</span><span>,</span><span> </span><span>size</span><span> </span><span>=</span><span> </span><span>1</span><span>)</span><span> </span><span>+</span><span> </span><span>geom_point</span><span>(</span><span>aes</span><span>(</span><span>y</span><span> </span><span>=</span><span> </span><span>log_sleep_total</span><span>))</span><span> </span><span>+</span><span> </span><span>scale_x_continuous</span><span>(</span><span>labels</span><span> </span><span>=</span><span> </span><span>function</span><span>(</span><span>x</span><span>)</span><span> </span><span>10</span><span> </span><span>^</span><span> </span><span>x</span><span>)</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>lab_lines</span><span>$</span><span>brain_log</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>lab_lines</span><span>$</span><span>sleep_log</span><span>)</span><span>
</span><span>p_linpread</span><span>
</span></code></pre>
</div>
<p><img title="Brain mass by sleep hours, log-10 scale, plus the median and 95% uncertainty interval for the model-predicted mean." alt="Brain mass by sleep hours, log-10 scale, plus the median and 95% uncertainty interval for the model-predicted mean." src="https://i1.wp.com/tjmahr.github.io/figs//2016-11-18-visualizing-uncertainty-rstanarm/posterior-linpred-plot-1.png?w=80%25&ssl=1" data-recalc-dims="1"></p>
<p>This plot is just like the <code>stat_smooth()</code> plot, except the interval here is<br>
interpreted in terms of post-data probabilities: We’re 95% certain—given the<br>
data, model and our prior information—that the “true” average sleep duration<br>
is contained in this interval. I put “true” in quotes because this is truth in<br>
the “small world” of the model, to quote<br>
<a target="_blank" rel="nofollow" href="http://xcelab.net/rm/statistical-rethinking/"><em>Statistical Rethinking</em></a>, not<br>
necessarily the real world.</p>
<p>Although the interpretation of the interval changes (compared to a classical<br>
confidence interval), its location barely changes at all. If we overlay a<br>
<code>stat_smooth()</code> layer onto this plot, we can see that two sets of intervals are<br>
virtually identical. With this much data and for this simple of a model, both<br>
types of models can make very similar estimates.</p>
<div>
<pre><code><span>p_linpread</span><span> </span><span>+</span><span> </span><span>stat_smooth</span><span>(</span><span>aes</span><span>(</span><span>y</span><span> </span><span>=</span><span> </span><span>log_sleep_total</span><span>),</span><span> </span><span>method</span><span> </span><span>=</span><span> </span><span>&quot;lm&quot;</span><span>)</span><span>
</span></code></pre>
</div>
<p><img title="Previous line-plus-interval plot with the classical regression line and confidence interval overlaid." alt="Previous line-plus-interval plot with the classical regression line and confidence interval overlaid." src="https://i2.wp.com/tjmahr.github.io/figs//2016-11-18-visualizing-uncertainty-rstanarm/posterior-linpred-plot-and-smooth-1.png?w=80%25&ssl=1" data-recalc-dims="1"></p>
<p>The previous plot illustrates one limitation of this approach: Pragmatically<br>
speaking, <code>stat_smooth()</code> basically does the same thing, and we’re<br>
not taking advantage of the affordances provided by our model. This is why<br>
RStanARM, in a kind of amusing way, disowns <code>posterior_linpred()</code> in its<br>
documentation:</p>
<blockquote>
<p>This function is occasionally convenient, but it should be used sparingly.<br>
Inference and model checking should generally be carried out using the<br>
posterior predictive distribution (see <code>posterior_predict</code>).</p>
</blockquote>
<p><em>Occasionally convenient.</em> :open_mouth: And elsewhere:</p>
<blockquote>
<p>See also: <code>posterior_predict</code> to draw from the posterior predictive<br>
distribution of the outcome, which is almost always preferable.</p>
</blockquote>
<h2>Option 3: Mean and 95% interval for model-generated data</h2>
<p>The reason why <code>posterior_predict()</code> is preferable is that it uses more<br>
information from our model, namely the error term <code>sigma</code>.<br>
<code>poseterior_linpred()</code> predicts averages; <code>posterior_predict()</code> predicts new<br>
observations. This <em>posterior predictive checking</em> helps us confirm whether our<br>
model—a story of how the data could have been generated—can produce new data<br>
that resembles our data.</p>
<p>Here, we can use the function we defined earlier to get prediction intervals.</p>
<div>
<pre><code><span># Still a matrix with one row per posterior draw and one column per observation
</span><span>pred_post</span><span> </span><span>&lt;-</span><span> </span><span>posterior_predict</span><span>(</span><span>m</span><span>1</span><span>,</span><span> </span><span>newdata</span><span> </span><span>=</span><span> </span><span>new_data</span><span>)</span><span>
</span><span>dim</span><span>(</span><span>pred_post</span><span>)</span><span>
</span><span>#> [1] 4000 80
</span><span>
</span><span>df_pred_post</span><span> </span><span>&lt;-</span><span> </span><span>tidy_predictions</span><span>(</span><span>pred_post</span><span>,</span><span> </span><span>new_data</span><span>)</span><span>
</span><span>df_pred_post</span><span>
</span><span>#> # A tibble: 80 × 5
#> observation median lower upper log_brainwt
#> &lt;int> &lt;dbl> &lt;dbl> &lt;dbl> &lt;dbl>
#> 1 1 1.224866 0.8685090 1.577798 -3.853872
#> 2 2 1.207392 0.8395285 1.560691 -3.795509
#> 3 3 1.209352 0.8499785 1.569175 -3.737146
#> 4 4 1.203873 0.8333415 1.563349 -3.678784
#> 5 5 1.204020 0.8537000 1.554171 -3.620421
#> 6 6 1.183633 0.8284588 1.552674 -3.562058
#> 7 7 1.182420 0.8234048 1.549418 -3.503695
#> 8 8 1.177556 0.8111187 1.543201 -3.445332
#> 9 9 1.164234 0.8238208 1.524496 -3.386970
#> 10 10 1.161509 0.8130019 1.526353 -3.328607
#> # ... with 70 more rows
</span></code></pre>
</div>
<p>And we can plot the interval in the same way.</p>
<div>
<pre><code><span>ggplot</span><span>(</span><span>msleep</span><span>)</span><span> </span><span>+</span><span> </span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>log_brainwt</span><span>)</span><span> </span><span>+</span><span> </span><span>geom_ribbon</span><span>(</span><span>aes</span><span>(</span><span>ymin</span><span> </span><span>=</span><span> </span><span>lower</span><span>,</span><span> </span><span>ymax</span><span> </span><span>=</span><span> </span><span>upper</span><span>),</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>df_pred_post</span><span>,</span><span> </span><span>alpha</span><span> </span><span>=</span><span> </span><span>0.4</span><span>,</span><span> </span><span>fill</span><span> </span><span>=</span><span> </span><span>&quot;grey60&quot;</span><span>)</span><span> </span><span>+</span><span> </span><span>geom_line</span><span>(</span><span>aes</span><span>(</span><span>y</span><span> </span><span>=</span><span> </span><span>median</span><span>),</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>df_pred_post</span><span>,</span><span> </span><span>colour</span><span> </span><span>=</span><span> </span><span>&quot;#3366FF&quot;</span><span>,</span><span> </span><span>size</span><span> </span><span>=</span><span> </span><span>1</span><span>)</span><span> </span><span>+</span><span> </span><span>geom_point</span><span>(</span><span>aes</span><span>(</span><span>y</span><span> </span><span>=</span><span> </span><span>log_sleep_total</span><span>))</span><span> </span><span>+</span><span> </span><span>scale_x_continuous</span><span>(</span><span>labels</span><span> </span><span>=</span><span> </span><span>function</span><span>(</span><span>x</span><span>)</span><span> </span><span>10</span><span> </span><span>^</span><span> </span><span>x</span><span>)</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>lab_lines</span><span>$</span><span>brain_log</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>lab_lines</span><span>$</span><span>sleep_log</span><span>)</span><span>
</span></code></pre>
</div>
<p><img title="Brain mass by sleep hours, log-10 scale, plus the median and 95% interval for posterior predicted observations." alt="Brain mass by sleep hours, log-10 scale, plus the median and 95% interval for posterior predicted observations." src="https://i2.wp.com/tjmahr.github.io/figs//2016-11-18-visualizing-uncertainty-rstanarm/posterior-predict-1.png?w=80%25&ssl=1" data-recalc-dims="1"></p>
<p>First, we can appreciate that this interval is <strong>much wider</strong>. That’s because the<br>
interval doesn’t summarize a particular statistic (like an average) but all of<br>
the observations that can generated by our model. Okay, not <em>all</em> of the<br>
observations—just the 95% most probable observations.</p>
<p>Next, we can also appreciate that the line and the ribbon are <strong>jagged due to<br>
simulation randomness</strong>. Each prediction is a random number draw, and at each<br>
value of <em>x</em>, we have 4000 such random draws. We computed a median and 95%<br>
interval at each <em>x</em>, but due to randomness from simulating new data, these<br>
medians do not smoothly connect together in the plot. That’s okay, because these<br>
fluctuations are relatively small.</p>
<p>Finally, we can see that there are only two <strong>points outside of the interval</strong>.<br>
These appear to be the restless roe deer and the ever-sleepy giant armadillo.<br>
These two represent the main outliers for our model because they fall slight<br>
outside of the 95% prediction interval. In this way, the posterior predictive<br>
interval can help us discover which data points are relative outliers for our<br>
model.</p>
<p>(Maybe <em>outliers</em> isn’t the right word. It makes perfect sense that 2/56 =<br>
3.6% of the observations fall outside of the 95%<br>
interval.)</p>
<div>
<pre><code>#> # A tibble: 5 × 4
#>              name sleep_total log_sleep_total log_brainwt
#>             &lt;chr>       &lt;dbl>           &lt;dbl>       &lt;dbl>
#> 1        Roe deer         3.0       0.4771213  -1.0078885
#> 2            Goat         5.3       0.7242759  -0.9393022
#> 3             Dog        10.1       1.0043214  -1.1549020
#> 4    Patas monkey        10.9       1.0374265  -0.9393022
#> 5 Giant armadillo        18.1       1.2576786  -1.0915150
</code></pre>
</div>
<p>This posterior prediction plot does reveal a shortcoming of our model, when<br>
plotted in a different manner.</p>
<div>
<pre><code><span>last_plot</span><span>()</span><span> </span><span>+</span><span> </span><span>geom_hline</span><span>(</span><span>yintercept</span><span> </span><span>=</span><span> </span><span>log10</span><span>(</span><span>24</span><span>),</span><span> </span><span>color</span><span> </span><span>=</span><span> </span><span>&quot;grey50&quot;</span><span>)</span><span> </span><span>+</span><span> </span><span>geom_label</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>0</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>log10</span><span>(</span><span>24</span><span>),</span><span> </span><span>label</span><span> </span><span>=</span><span> </span><span>&quot;24 hours&quot;</span><span>)</span><span>
</span></code></pre>
</div>
<p><img title="Previous plot updated to include a line indicating 24 hours. Same of the 95% interval goes above the line." alt="Previous plot updated to include a line indicating 24 hours. Same of the 95% interval goes above the line." src="https://i0.wp.com/tjmahr.github.io/figs//2016-11-18-visualizing-uncertainty-rstanarm/posterior-predict-24-hours-1.png?w=80%25&ssl=1" data-recalc-dims="1"></p>
<p>One faulty consequence of how our model was specified is that it predicts that<br>
some mammals sleep more than 24 hours per day—oh, what a life to live<br>
:sleeping:.</p>
<h2>Wrap up</h2>
<p>In the post, I covered three different ways to plot the results of an RStanARM<br>
model, while demonstrating some of the key functions for working with RStanARM<br>
models. Time well spent, I think.</p>
<p>As for future directions, I learned about the under-development (as of November<br>
2016) R package <a target="_blank" rel="nofollow" href="https://github.com/stan-dev/bayesplot">bayesplot</a> by the Stan<br>
team. The README package shows off a lot of different ways to visualize<br>
posterior samples from a model. I’ll be sure to demo it on this data-set once it<br>
goes live.</p>
<div>
<ol>
<li>
<p>That is, if we map the plot’s color aesthetic to a categorical variable in the data, <code>stat_smooth()</code> will fit a separate model for each color/category. I figured this out when I tried to write my own function <code>stat_smooth_stan()</code> based on <a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/ggplot2/vignettes/extending-ggplot2.html">ggplot2’s extensions vignette</a> and noticed that RStanARM was printing out MCMC sampling information for each color/category of the data. <a target="_blank" rel="nofollow" href="https://tjmahr.github.io/visualizing-uncertainty-rstanarm/#fnref:1">↩</a></p>
</li>
</ol>
</div>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="https://www.r-bloggers.com/rstanarm-basics-visualizing-uncertainty-in-linear-regression/"></div></div><div><a data-url="https://www.r-bloggers.com/rstanarm-basics-visualizing-uncertainty-in-linear-regression/" data-counturl="https://www.r-bloggers.com/rstanarm-basics-visualizing-uncertainty-in-linear-regression/" data-count="vertical" data-text="RStanARM basics: visualizing uncertainty in linear regression" 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="https://tjmahr.github.io/visualizing-uncertainty-rstanarm/"> Higher Order Functions</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/1lddd9YFQ5Q"></a></div>