<p>As part of <a href="https://tjmahr.github.io/rstanarm-tutorial-slides/">my tutorial talk on RStanARM</a>, I
presented some examples of how to visualize the uncertainty in Bayesian linear
regression models. This post is an expanded demonstration of the approaches I
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 href="https://vincentarelbundock.github.io/Rdatasets/doc/ggplot2/msleep.html">mammal sleep dataset from ggplot2</a>. This dataset
contains the number of hours spent sleeping per day for 83 different species of
mammals along with each species’ brain mass (kg) and body mass (kg), among other
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></p>
<p>Hmmm, not very helpful! We should put our measures on a log-10 scale. Also, 27
of the species don’t have brain mass data, so we’ll exclude those rows for the
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
the points for some example critters :cat: so that we can get some intuition
about the data in this scaling. (Plus, I wanted to try out the <a href="https://tjmahr.github.io/r4ds.had.co.nz/graphics-for-communication.html#annotations">annotation
tips</a> from the R4DS 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></p>
<p>As a child growing up on a dairy farm :cow:, it was remarkable to me how little I saw cows sleeping, compared to dogs or cats. Were they okay? Are they constantly tired and groggy? Maybe they are asleep when I’m asleep? Here, it
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
measure so that the regression line doesn’t imply negative sleep (even though
brains never get that 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) of brain mass sleeps 100.74 = 5.5 hours per
day. A mammal with a tenth of that brain mass (-1 log-kg) sleeps 100.74 + 0.13 = 7.4 hours.</p>
<p>We illustrate the regression results to show the predicted mean of y and
its 95% confidence interval. This task is readily accomplished in ggplot2 using <code>stat_smooth()</code>. This function fits a model and plots the mean and CI for each aesthetic grouping of data<a href="https://tjmahr.github.io/#fn:1">1</a> 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></p>
<p>This interval conveys some uncertainty in the estimate of the mean, but this interval has a <a href="https://stats.stackexchange.com/questions/2272/whats-the-difference-between-a-confidence-interval-and-a-credible-interval">frequentist interpretation</a> which can be
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
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
line of best fit that satisfies a least-squares or maximum-likelihood objective. Our Bayesian model estimates an entire distribution of plausible
regression lines. The first way to visualize our uncertainty is to plot our
own line of best fit along with a sample of other lines from the posterior
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 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
the classical model’s intercept and slope. The median line serves as the “point
estimate” for our model: If we had to summarize the modeled relationship using
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
plus a sample of the other credible lines from our model</strong>. First, we create a 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, 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></p>
<p>Each of these light lines represents a credible prediction of the mean across the values of x. As these line pile up on top of each other, they create an uncertainty band around our line of best fit. More plausible lines are more likely to be sampled, so these lines overlap and create a uniform color around
the median line. As we move left or right, getting farther away from the mean of
x, the lines start to fan out and we see very faint individual lines for some
of the more extreme (yet still plausible) lines.</p>
<p>The advantage of this plot is that it is a direct visualization of posterior samples—one line per sample. It provides an estimate for the central tendency
in the data but it also converys uncertainty around that estimate.</p>
<p>This approach has limitations, however. Lines for subgroups require a little more effort to undo interactions. Also, the regression lines span the whole x axis which is not appropriate when subgroups only use a portion of the x-axis. (This limitation is solvable though.) Finally, I haven’t found good defaults for the aesthetic options: The number of samples, the colors to use, and the transparency level. One can lose of lots and lots and lots time fiddling with
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
best fit and the 95% uncertainty interval around it.</p>
<p>To limit the amount of the x axis used by the lines, we’re going to create a 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
of new data. I say means because the function computes 80 predicted means for
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 point. I do some tidying to get the data into a long format (one row per fitted mean per posterior sample), and then do a table-join with the <code>observation</code> column included in <code>new_data</code>. I store these steps in a function because I
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
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></p>
<p>This plot is just like the <code>stat_smooth()</code> plot, except the interval here is interpreted in terms of post-data probabilities: We’re 95% certain—given the data, model and our prior information—that the “true” average sleep duration is contained in this interval. I put “true” in quotes because this is truth in
the “small world” of the model, to quote <a href="http://xcelab.net/rm/statistical-rethinking/">Statistical Rethinking</a>, not necessarily the real world.</p>
<p>Although the interpretation of the interval changes (compared to a classical confidence interval), its location barely changes at all. If we overlay a <code>stat_smooth()</code> layer onto this plot, we can see that two sets of intervals are virtually identical. With this much data and for this simple of a model, both
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></p>
<p>The previous plot illustrates one limitation of this approach: Pragmatically
speaking, <code>stat_smooth()</code> basically does the same thing, and we’re
not taking advantage of the affordances provided by our model. This is why RStanARM, in a kind of amusing way, disowns <code>posterior_linpred()</code> in its documentation:</p>
<blockquote>
<p>This function is occasionally convenient, but it should be used sparingly. Inference and model checking should generally be carried out using the posterior predictive distribution (see <code>posterior_predict</code>).</p>
</blockquote>
<p>Occasionally convenient. :open_mouth: And elsewhere:</p>
<blockquote>
<p>See also: <code>posterior_predict</code> to draw from the posterior predictive 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 information from our model, namely the error term <code>sigma</code>. <code>poseterior_linpred()</code> predicts averages; <code>posterior_predict()</code> predicts new observations. This posterior predictive checking helps us confirm whether our model—a story of how the data could have been generated—can produce new data
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></p>
<p>First, we can appreciate that this interval is <strong>much wider</strong>. That’s because the
interval doesn’t summarize a particular statistic (like an average) but all of
the observations that can generated by our model. Okay, not all of the
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 simulation randomness</strong>. Each prediction is a random number draw, and at each
value of x, we have 4000 such random draws. We computed a median and 95% interval at each x, but due to randomness from simulating new data, these
medians do not smoothly connect together in the plot. That’s okay, because these
fluctuations are relatively small.</p>
<p>Finally, we can see that there are only two <strong>points outside of the interval</strong>. These appear to be the restless roe deer and the ever-sleepy giant armadillo.
These two represent the main outliers for our model because they fall slight
outside of the 95% prediction interval. In this way, the posterior predictive
interval can help us discover which data points are relative outliers for our
model.</p>
<p>(Maybe outliers isn’t the right word. It makes perfect sense that 2/56 =
3.6% of the observations fall outside of the 95%
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
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></p>
<p>One faulty consequence of how our model was specified is that it predicts that
some mammals sleep more than 24 hours per day—oh, what a life to live :sleeping:.</p>
<h2>Wrap up</h2>
<p>In the post, I covered three different ways to plot the results of an RStanARM model, while demonstrating some of the key functions for working with RStanARM models. Time well spent, I think.</p>
<p>As for future directions, I learned about the under-development (as of November
2016) R package <a href="https://github.com/stan-dev/bayesplot">bayesplot</a> by the Stan
team. The README package shows off a lot of different ways to visualize
posterior samples from a model. I’ll be sure to demo it on this data-set once it
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 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 href="https://tjmahr.github.io/#fnref:1">↩</a></p>
</li>
</ol></div>