<div>
<p>The Beta-Binomial model is the “hello world” of Bayesian statistics. That is, it’s the first model you get to run, often before you even know what you are doing. There are many reasons for this: </p>
<ul>
<li>It only has one parameter, the underlying proportion of success, so it’s easy to visualize and reason about.</li>
<li>It’s easy to come up with a scenario where it can be used, for example: “What is the proportion of patients that will be cured by this drug?”</li>
<li>The model can be computed analytically (no need for any messy MCMC).</li>
<li>It’s relatively easy to come up with an informative prior for the underlying proportion.</li>
<li>Most importantly: It’s fun to see some results before diving into the theory! <img alt="😁" src="https://s.w.org/images/core/emoji/2.4/72x72/1f601.png" class="wp-smiley"></li>
</ul>
<p>That’s why I also introduced the Beta-Binomial model as the first model in my DataCamp course <a rel="nofollow" href="http://www.sumsar.net/blog/2018/12/my-introductory-course-on-bayesian-statistics/" target="_blank">Fundamentals of Bayesian Data Analysis in R</a> and quite a lot of people have asked me for the code I used to visualize the Beta-Binomial. Scroll to the bottom of this post if that’s what you want, otherwise, here is how I visualized the Beta-Binomial in my course given two successes and four failures:</p>
<p><center>
<p>
<img src="https://i1.wp.com/www.sumsar.net/images/posts/2018-12-13-visualizing-the-beta-binomial/prop_model_animation.gif?w=450" data-recalc-dims="1">
</p>
<p></p>
<p></p>
<p>The function that produces these plots is called <code>prop_model</code> (<code>prop</code> as in proportion) and takes a vector of <code>TRUE</code>s and <code>FALSE</code>s representing successes and failures. The visualization is created using the excellent <a rel="nofollow" href="https://cran.r-project.org/package=ggridges" target="_blank"><code>ggridges</code> package</a> (<a rel="nofollow" href="https://serialmentor.com/blog/2017/9/15/goodbye-joyplots" target="_blank">previously called joyplot</a>). Here’s how you would use <code>prop_model</code> to produce the last plot in the animation above:</p>
<pre><code><div><pre>data <span>&lt;-</span> c(<span>FALSE</span>, <span>TRUE</span>, <span>FALSE</span>, <span>FALSE</span>, <span>FALSE</span>, <span>TRUE</span>)
prop_model(data)
</pre>
</div>
<p></p>
<p><img alt="" src="https://i0.wp.com/www.sumsar.net/figures/2018-12-13-visualizing-the-beta-binomial/unnamed-chunk-2-1.png?w=456" data-recalc-dims="1"></p>
<p>The result is, I think, a quite nice visualization of how the model’s knowledge about the parameter changes as data arrives. At <code>n=0</code> the model doesn’t know anything and — as the default prior states that it’s equally likely the proportion of success is anything from 0.0 to 1.0 — the result is a big, blue, and uniform square. As more data arrives the probability distribution becomes more concentrated, with the final posterior distribution at <code>n=6</code>. </p>
<p>Some added features of <code>prop_model</code> is that it also plots larger data somewhat gracefully and that it returns a random sample from the posterior that can be further explored. For example: </p>
<pre><code><div><pre>big_data <span>&lt;-</span> sample(c(<span>TRUE</span>, <span>FALSE</span>), prob <span>=</span> c(<span>0.75</span>, <span>0.25</span>),
                   size <span>=</span> <span>100</span>, replace <span>=</span> <span>TRUE</span>)
posterior <span>&lt;-</span> prop_model(big_data)
</pre>
</div>
<p></p>
<p><img alt="" src="https://i2.wp.com/www.sumsar.net/figures/2018-12-13-visualizing-the-beta-binomial/unnamed-chunk-3-1.png?w=456" data-recalc-dims="1"></p>
<pre><code><div><pre>quantile(posterior, c(<span>0.025</span>, <span>0.5</span>, <span>0.975</span>))
</pre>
</div>
<p></p>
<pre><code><div><pre><span>## 2.5% 50% 98% </span>
<span>## 0.68 0.77 0.84</span>
</pre>
</div>
<p></p>
<p>So here we calculated that the underlying proportion of success is most likely 0.77 with a 95% CI of [0.68, 0.84] (which nicely includes the correct value of 0.75 which we used to simulate <code>big_data</code>).</p>
<p>To be clear, <code>prop_model</code> is not intended as anything serious, it’s just meant as a nice way of exploring the Beta-Binomial model when learning Bayesian statistics, maybe as part of a workshop exercise.</p>
<h2>The <code>prop_model</code> function</h2>
<pre><code><div><pre><span># This function takes a number of successes and failuers coded as a TRUE/FALSE</span>
<span># or 0/1 vector. This should be given as the data argument.</span>
<span># The result is a visualization of the how a Beta-Binomial</span>
<span># model gradualy learns the underlying proportion of successes </span>
<span># using this data. The function also returns a sample from the</span>
<span># posterior distribution that can be further manipulated and inspected.</span>
<span># The default prior is a Beta(1,1) distribution, but this can be set using the</span>
<span># prior_prop argument.</span>

<span># Make sure the packages tidyverse and ggridges are installed, otherwise run:</span>
<span># install.packages(c(&quot;tidyverse&quot;, &quot;ggridges&quot;))</span>

<span># Example usage:</span>
<span># data &lt;-c(TRUE, FALSE, TRUE, TRUE, FALSE, TRUE, TRUE)</span>
<span># prop_model(data)</span>
prop_model <span>&lt;-</span> <span>function</span>(data <span>=</span> c(), prior_prop <span>=</span> c(<span>1</span>, <span>1</span>), n_draws <span>=</span> <span>10000</span>) {
  library(tidyverse)
  data <span>&lt;-</span> as.logical(data)
  <span># data_indices decides what densities to plot between the prior and the posterior</span>
  <span># For 20 datapoints and less we're plotting all of them.</span>
  data_indices <span>&lt;-</span> round(seq(<span>0</span>, length(data), length.out <span>=</span> min(length(data) <span>+</span> <span>1</span>, <span>20</span>)))

  <span># dens_curves will be a data frame with the x &amp; y coordinates for the </span>
  <span># denities to plot where x = proportion_success and y = probability</span>
  proportion_success <span>&lt;-</span> c(<span>0</span>, seq(<span>0</span>, <span>1</span>, length.out <span>=</span> <span>100</span>), <span>1</span>)
  dens_curves <span>&lt;-</span> map_dfr(data_indices, <span>function</span>(i) {
    value <span>&lt;-</span> ifelse(i <span>==</span> <span>0</span>, <span>&quot;Prior&quot;</span>, ifelse(data[i], <span>&quot;Success&quot;</span>, <span>&quot;Failure&quot;</span>))
    label <span>&lt;-</span> paste0(<span>&quot;n=&quot;</span>, i)
    probability <span>&lt;-</span> dbeta(proportion_success,
                         prior_prop[<span>1</span>] <span>+</span> sum(data[seq_len(i)]),
                         prior_prop[<span>2</span>] <span>+</span> sum(<span>!</span>data[seq_len(i)]))
    probability <span>&lt;-</span> probability <span>/</span> max(probability)
    data_frame(value, label, proportion_success, probability)
  })
  <span># Turning label and value into factors with the right ordering for the plot</span>
  dens_curves<span>$</span>label <span>&lt;-</span> fct_rev(factor(dens_curves<span>$</span>label, levels <span>=</span>  paste0(<span>&quot;n=&quot;</span>, data_indices )))
  dens_curves<span>$</span>value <span>&lt;-</span> factor(dens_curves<span>$</span>value, levels <span>=</span> c(<span>&quot;Prior&quot;</span>, <span>&quot;Success&quot;</span>, <span>&quot;Failure&quot;</span>))

  p <span>&lt;-</span> ggplot(dens_curves, aes(x <span>=</span> proportion_success, y <span>=</span> label,
                               height <span>=</span> probability, fill <span>=</span> value)) <span>+</span>
    ggridges<span>::</span>geom_density_ridges(stat<span>=</span><span>&quot;identity&quot;</span>, color <span>=</span> <span>&quot;white&quot;</span>, alpha <span>=</span> <span>0.8</span>,
                                  panel_scaling <span>=</span> <span>TRUE</span>, size <span>=</span> <span>1</span>) <span>+</span>
    scale_y_discrete(<span>&quot;&quot;</span>, expand <span>=</span> c(<span>0.01</span>, <span>0</span>)) <span>+</span>
    scale_x_continuous(<span>&quot;Underlying proportion of success&quot;</span>) <span>+</span>
    scale_fill_manual(values <span>=</span> hcl(<span>120</span> <span>*</span> <span>2:0</span> <span>+</span> <span>15</span>, <span>100</span>, <span>65</span>), name <span>=</span> <span>&quot;&quot;</span>, drop <span>=</span> <span>FALSE</span>,
                      labels <span>=</span>  c(<span>&quot;Prior &quot;</span>, <span>&quot;Success &quot;</span>, <span>&quot;Failure &quot;</span>)) <span>+</span>
    ggtitle(paste0(
      <span>&quot;Binomial model - Data: &quot;</span>, sum(data),  <span>&quot; successes, &quot;</span> , sum(<span>!</span>data), <span>&quot; failures&quot;</span>)) <span>+</span>
    theme_light() <span>+</span>
    theme(legend.position <span>=</span> <span>&quot;top&quot;</span>)
  print(p)

  <span># Returning a sample from the posterior distribution that can be further </span>
  <span># manipulated and inspected</span>
  posterior_sample <span>&lt;-</span> rbeta(n_draws, prior_prop[<span>1</span>] <span>+</span> sum(data), prior_prop[<span>2</span>] <span>+</span> sum(<span>!</span>data))
  invisible(posterior_sample)
}
</pre>
</div>
<p></p>
<p></p>

				
		
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://www.sumsar.net/blog/2018/12/visualizing-the-beta-binomial/"> Publishable Stuff</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 rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/_gRzIJVItGU" height="1"></a></div></code></pre></code></pre></code></pre></code></pre></code></pre></center></div>