Fancy Plot (with Posterior Samples) for Bayesian Regressions

[This article was first published on Dominique Makowski, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

As Bayesian models usually generate a lot of samples (iterations), one could want to plot them as well, instead (or along) the posterior “summary” (with indices like the 90% HDI). This can be done quite easily by extracting all the iterations in get_predicted from the psycho package.

The Model

<span class="c1"># devtools::install_github("neuropsychology/psycho.R")  # Install the last psycho version if needed</span><span class="w">

</span><span class="c1"># Load packages</span><span class="w">
</span><span class="n">library</span><span class="p">(</span><span class="n">tidyverse</span><span class="p">)</span><span class="w">
</span><span class="n">library</span><span class="p">(</span><span class="n">psycho</span><span class="p">)</span><span class="w">

</span><span class="c1"># Import data</span><span class="w">
</span><span class="n">df</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">psycho</span><span class="o">::</span><span class="n">affective</span><span class="w">

</span><span class="c1"># Fit a logistic regression model</span><span class="w">
</span><span class="n">fit</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">rstanarm</span><span class="o">::</span><span class="n">stan_glm</span><span class="p">(</span><span class="n">Sex</span><span class="w"> </span><span class="o">~</span><span class="w"> </span><span class="n">Adjusting</span><span class="p">,</span><span class="w"> </span><span class="n">data</span><span class="o">=</span><span class="n">df</span><span class="p">,</span><span class="w"> </span><span class="n">family</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="s2">"binomial"</span><span class="p">)</span><span class="w">
</span>

We fitted a Bayesian logistic regression to predict the sex (W / M) with one’s ability to flexibly adjust to his/her emotional reaction.

Plot

To visualize the model, the most neat way is to extract a “reference grid” (i.e., a theorethical dataframe with balanced data). Our refgrid is made of equally spaced predictor values. With it, we can make predictions using the previously fitted model. This will compute the median of the posterior prediction, as well as the 90% credible interval. However, we’re interested in keeping all the prediction samples (iterations). Note that get_predicted automatically transformed log odds ratios (the values in which the model is expressed) to probabilities, easier to apprehend.

<span class="c1"># Generate a new refgrid</span><span class="w">
</span><span class="n">refgrid</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">df</span><span class="w"> </span><span class="o">%>%</span><span class="w"> 
  </span><span class="n">dplyr</span><span class="o">::</span><span class="n">select</span><span class="p">(</span><span class="n">Adjusting</span><span class="p">)</span><span class="w"> </span><span class="o">%>%</span><span class="w"> 
  </span><span class="n">psycho</span><span class="o">::</span><span class="n">refdata</span><span class="p">(</span><span class="n">length.out</span><span class="o">=</span><span class="m">10</span><span class="p">)</span><span class="w">

</span><span class="c1"># Get predictions and keep iterations</span><span class="w">
</span><span class="n">predicted</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">psycho</span><span class="o">::</span><span class="n">get_predicted</span><span class="p">(</span><span class="n">fit</span><span class="p">,</span><span class="w"> </span><span class="n">newdata</span><span class="o">=</span><span class="n">refgrid</span><span class="p">,</span><span class="w"> </span><span class="n">keep_iterations</span><span class="o">=</span><span class="kc">TRUE</span><span class="p">)</span><span class="w">

</span><span class="c1"># Reshape this dataframe to have iterations as factor</span><span class="w">
</span><span class="n">predicted</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">predicted</span><span class="w"> </span><span class="o">%>%</span><span class="w"> 
  </span><span class="n">tidyr</span><span class="o">::</span><span class="n">gather</span><span class="p">(</span><span class="n">Iteration</span><span class="p">,</span><span class="w"> </span><span class="n">Iteration_Value</span><span class="p">,</span><span class="w"> </span><span class="n">starts_with</span><span class="p">(</span><span class="s2">"iter"</span><span class="p">))</span><span class="w">

</span><span class="c1"># Plot all iterations with the median prediction</span><span class="w">
</span><span class="n">ggplot</span><span class="p">(</span><span class="n">predicted</span><span class="p">,</span><span class="w"> </span><span class="n">aes</span><span class="p">(</span><span class="n">x</span><span class="o">=</span><span class="n">Adjusting</span><span class="p">))</span><span class="w"> </span><span class="o">+</span><span class="w">
  </span><span class="n">geom_line</span><span class="p">(</span><span class="n">aes</span><span class="p">(</span><span class="n">y</span><span class="o">=</span><span class="n">Iteration_Value</span><span class="p">,</span><span class="w"> </span><span class="n">group</span><span class="o">=</span><span class="n">Iteration</span><span class="p">),</span><span class="w"> </span><span class="n">size</span><span class="o">=</span><span class="m">0.3</span><span class="p">,</span><span class="w"> </span><span class="n">alpha</span><span class="o">=</span><span class="m">0.01</span><span class="p">)</span><span class="w"> </span><span class="o">+</span><span class="w">
  </span><span class="n">geom_line</span><span class="p">(</span><span class="n">aes</span><span class="p">(</span><span class="n">y</span><span class="o">=</span><span class="n">Sex_Median</span><span class="p">),</span><span class="w"> </span><span class="n">size</span><span class="o">=</span><span class="m">1</span><span class="p">)</span><span class="w"> </span><span class="o">+</span><span class="w"> 
  </span><span class="n">ylab</span><span class="p">(</span><span class="s2">"Probability of being a man\n"</span><span class="p">)</span><span class="w"> </span><span class="o">+</span><span class="w">
  </span><span class="n">theme_classic</span><span class="p">()</span><span class="w">
</span>

Credits

This package helped you? Don’t forget to cite the various packages you used 🙂

You can cite psycho as follows:

  • Makowski, (2018). The psycho Package: an Efficient and Publishing-Oriented Workflow for Psychological Science. Journal of Open Source Software, 3(22), 470. https://doi.org/10.21105/joss.00470

To leave a comment for the author, please follow the link and comment on their blog: Dominique Makowski.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)