<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/bayesian-binomial-test-in-r/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://silentspring.org/devblog/2018/01/11/bayesian-binomial-test.html"> Silent Spring Institute Developer Blog</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><strong>Summary</strong>: in this post, I implemenent an R function for computing \( P(\theta_1 > \theta2) \), where \( \theta_1 \) and \( \theta_2 \) are beta-distributed random variables. This is useful for estimating the probability that one binomial proportion is greater than another.</p>
<p>I am working on a project in which I need to compare two binomial proportions to see if one is likely to be greater than the other. I wanted to do this in a Bayesian framework; fortunately, estimating binomial proportions is one of the first subjects taught to introduce Bayesian statistics (and statistics in general, for that matter) so the method is well established. <a rel="nofollow" href="http://www.sumsar.net/blog/2014/01/bayesian-first-aid-binomial-test/" target="_blank">Rasmus Bååth</a> has a very nice article describing the Bayesian binomial test, and an estimation approach using JAGS.</p>
<p>I made the problem even simpler by using the fact that the beta distribution is the conjugate prior to the binomial distribution. That is, if the prior is \( \mathrm{beta}(\alpha, \beta) \) distributed, then the posterior after observing \( k \) successes in \( n \) trials is</p>
<p></p>
<p>If you want to start with a noninformative prior, then you can use the beta distribution with \( \alpha = \beta = 1\), which collapses to the uniform distribution. This makes the posterior</p>
<p></p>
<p>What I am interested is in comparing two binomial proportions, \( \theta_1 \) and \( \theta_2 \). I want to calculate</p>
<p></p>
<p>We observe \( k_1 \) successes out of \( n_1 \) trials and estimate a posterior distribution for \( \theta_1 \), and \( k_2 \) successes out of \( n_2 \) trials for \( \theta_2 \). Then the probability in question is</p>
<p></p>
<p>So far, so textbook. At this point, the simplest way to get an estimate is to draw two sets of random numbers, distributed by the two respective beta distributions, and see how many times the first random numbers is greater than the second:</p>
<figure>
<pre><code data-lang="r"><span>set.seed</span><span>(</span><span>1</span><span>)</span>

k1 <span>&lt;-</span> <span>5</span> <span># Five successes</span>
k2 <span>&lt;-</span> <span>1</span> <span># One success</span>
n1 <span>&lt;-</span> <span>10</span> <span># Ten trials</span>
n2 <span>&lt;-</span> <span>10</span> <span># Ten trials</span>

b1 <span>&lt;-</span> rbeta<span>(</span><span>10000</span><span>,</span> <span>1</span> <span>+</span> k1<span>,</span> <span>1</span> <span>+</span> n1 <span>-</span> k1<span>)</span>
b2 <span>&lt;-</span> rbeta<span>(</span><span>10000</span><span>,</span> <span>1</span> <span>+</span> k2<span>,</span> <span>1</span> <span>+</span> n2 <span>-</span> k2<span>)</span>

<span># Probability estimate of P(theta_1 > theta_2)</span>
<span>sum</span><span>(</span>b1 <span>></span> b2<span>)</span> <span>/</span> <span>10000</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## [1] 0.9688</code></pre>
</figure>
<p>This works fine, but it would be nice to have an exact solution. I found several papers describing different approaches.</p>
<p><span></span></p>
<p><a rel="nofollow" href="http://onlinelibrary.wiley.com/doi/10.1002/sim.2701/abstract" target="_blank">Nadarajah and Kotz (2007)</a> derived a functional form for the distribution of \( P(\theta_1 – \theta_2) \) (and <a rel="nofollow" href="http://onlinelibrary.wiley.com/doi/10.1002/sim.4248/full" target="_blank">Chen and Luo</a> fixed a typo in their work in 2011.) This approach is implemented in the <code>tolerance</code> package via the functions <code>ddiffprop</code>, <code>pdiffprop</code>, <code>qdiffprop</code>, and <code>rdiffprop</code>. The <code>a1</code> and <code>a2</code> parameters let you adjust the priors, where <code>a1 = a2 = 1</code> is a uniform prior for both proportions.</p>
<figure>
<pre><code data-lang="r"><span>library</span><span>(</span>tolerance<span>)</span>

<span># Probability that the first proportion (5 successes out of 10 trials)</span>
<span># is greater than the second (1 success out of 10 trials)</span>
pdiffprop<span>(</span><span>0</span><span>,</span>
          k1 <span>=</span> <span>5</span><span>,</span>
          k2 <span>=</span> <span>1</span><span>,</span>
          n1 <span>=</span> <span>10</span><span>,</span>
          n2 <span>=</span> <span>10</span><span>,</span>
          a1 <span>=</span> <span>1</span><span>,</span>
          a2 <span>=</span> <span>1</span><span>,</span>
          lower.tail <span>=</span> <span>FALSE</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## [1] 0.9682663</code></pre>
</figure>
<p>The functions work well, but I found that they stop working for large values of \( n_1 \) or \( n_2 \). This is because several of the terms in the density function use the gamma function, and the R implementation returns infinity and throws a warning for any input larger than 171.</p>
<p>The second approach I found is described in <a rel="nofollow" href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4019606/" target="_blank">Raineri et al. 2014</a>, in which they derived a recursive formula for \( P(\theta_1 > \theta_2) \). I unrolled the recursion so it could handle large \( n \) without overflowing the stack, and implemented it in R. The first two functions, <code>h</code> and <code>g</code>, come directly from Raineri et al. and <code>pdiff</code> provides a wrapper for estimating differences in binomial proportions.</p>
<figure>
<pre><code data-lang="r">h <span>&lt;-</span> <span>function</span><span>(</span>a<span>,</span> b<span>,</span> <span>c</span><span>,</span> d<span>)</span> <span>{</span>
  <span>exp</span><span>(</span><span>lbeta</span><span>(</span>a <span>+</span> <span>c</span><span>,</span> b <span>+</span> d<span>)</span> <span>-</span> <span>(</span><span>lbeta</span><span>(</span>a<span>,</span> b<span>)</span> <span>+</span> <span>lbeta</span><span>(</span><span>c</span><span>,</span> d<span>)))</span>
<span>}</span>

<span># probability that beta(a, b) > beta(c, d)</span>
g <span>&lt;-</span> <span>function</span><span>(</span>a<span>,</span> b<span>,</span> <span>c</span><span>,</span> d<span>)</span> <span>{</span>
  <span>stopifnot</span><span>(</span>a <span>>=</span> <span>1</span> <span>&amp;&amp;</span> b <span>>=</span> <span>1</span> <span>&amp;&amp;</span> <span>c</span> <span>>=</span> <span>1</span> <span>&amp;&amp;</span> d <span>>=</span> <span>1</span><span>)</span>
  accum <span>&lt;-</span> <span>0.5</span>
  <span>while</span> <span>(</span>a <span>!=</span> <span>c</span> <span>||</span> b <span>!=</span> d<span>)</span> <span>{</span>
    <span>if</span> <span>(</span>a <span>></span> <span>c</span><span>)</span> <span>{</span>
      accum <span>&lt;-</span> accum <span>+</span> h<span>(</span>a <span>-</span> <span>1</span><span>,</span> b<span>,</span> <span>c</span><span>,</span> d<span>)</span> <span>/</span> <span>(</span>a <span>-</span> <span>1</span><span>)</span>
      a <span>&lt;-</span> a <span>-</span> <span>1</span>
    <span>}</span> <span>else</span> <span>if</span> <span>(</span>b <span>></span> d<span>)</span> <span>{</span>
      accum <span>&lt;-</span> accum <span>-</span> h<span>(</span>a<span>,</span> b <span>-</span> <span>1</span><span>,</span> <span>c</span><span>,</span> d<span>)</span> <span>/</span> <span>(</span>b <span>-</span> <span>1</span><span>)</span>
      b <span>&lt;-</span> b <span>-</span> <span>1</span>
    <span>}</span> <span>else</span> <span>if</span> <span>(</span><span>c</span> <span>></span> a<span>)</span> <span>{</span>
      accum <span>&lt;-</span> accum <span>-</span> h<span>(</span>a<span>,</span> b<span>,</span> <span>c</span> <span>-</span> <span>1</span><span>,</span> d<span>)</span> <span>/</span> <span>(</span><span>c</span> <span>-</span> <span>1</span><span>)</span>
      <span>c</span> <span>&lt;-</span> <span>c</span> <span>-</span> <span>1</span>
    <span>}</span> <span>else</span> <span>if</span> <span>(</span>d <span>></span> b<span>)</span> <span>{</span>
      accum <span>&lt;-</span> accum <span>+</span> h<span>(</span>a<span>,</span> b<span>,</span> <span>c</span><span>,</span> d <span>-</span> <span>1</span><span>)</span> <span>/</span> <span>(</span>d <span>-</span> <span>1</span><span>)</span>
      d <span>&lt;-</span> d <span>-</span> <span>1</span>
    <span>}</span>
  <span>}</span>
  accum
<span>}</span>

pdiff <span>&lt;-</span> <span>function</span><span>(</span>k1<span>,</span> k2<span>,</span> n1<span>,</span> n2<span>,</span> a1 <span>=</span> <span>1</span><span>,</span> a2 <span>=</span> <span>1</span><span>)</span> <span>{</span>
  g<span>(</span><span>1</span> <span>+</span> k1<span>,</span> a1 <span>+</span> n1 <span>-</span> k1<span>,</span> <span>1</span> <span>+</span> k2<span>,</span> a2 <span>+</span> n2 <span>-</span> k2<span>)</span>
<span>}</span>

<span># Probability that the first proportion (5 successes out of 10 trials)</span>
<span># is greater than the second (1 success out of 10 trials)</span>
pdiff<span>(</span>k1 <span>=</span> <span>5</span><span>,</span>
      k2 <span>=</span> <span>1</span><span>,</span>
      n1 <span>=</span> <span>10</span><span>,</span>
      n2 <span>=</span> <span>10</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## [1] 0.9682663</code></pre>
</figure>
<p>It can also handle large values for \( n_1 \) and \( n_2 \):</p>
<figure>
<pre><code data-lang="r">pdiff<span>(</span>k1 <span>=</span> <span>420</span><span>,</span>
      k2 <span>=</span> <span>400</span><span>,</span>
      n1 <span>=</span> <span>1000</span><span>,</span>
      n2 <span>=</span> <span>1000</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## [1] 0.8182686</code></pre>
</figure>
<p>The results from the simulation, <code>tolerance</code> package, and this new function all agree, which is a good sign. I ended up using the <code>pdiff</code> function for my project.</p>
<h2>A baseball example</h2>
<p>Let’s apply this to a simple example. Suppose we want to estimate an “ability score” for every baseball team, so we can see which teams are better than others. We model each team’s ability score as a binomial random variable representing the team’s proportion of wins to total games. By feeding in each team’s win and loss record for a season, we’ll estimate a posterior distribution that we can use to rank pairs of teams by how likely they are to have different ability scores.</p>
<p>First, we filter for the team records from 2016:</p>
<figure>
<pre><code data-lang="r"><span>library</span><span>(</span>tidyverse<span>)</span>
<span>library</span><span>(</span>Lahman<span>)</span>
<span>library</span><span>(</span>knitr<span>)</span>

teams <span>&lt;-</span> Teams <span>%>%</span>
  filter<span>(</span>yearID <span>==</span> <span>2016</span><span>)</span> <span>%>%</span>
  select<span>(</span>teamID<span>,</span> W<span>,</span> L<span>)</span> <span>%>%</span>
  mutate<span>(</span>teamID <span>=</span> <span>as.character</span><span>(</span>teamID<span>),</span>
         P <span>=</span> <span>paste0</span><span>(</span><span>round</span><span>(</span>W <span>/</span> <span>(</span>W <span>+</span> L<span>)</span> <span>*</span> <span>100</span><span>),</span> <span>&quot;%&quot;</span><span>))</span></code></pre>
</figure>
<p>And generate a dataframe that has the win and loss numbers from every pair of teams:</p>
<figure>
<pre><code data-lang="r">teams_compared <span>&lt;-</span> <span>expand.grid</span><span>(</span>teams<span>$</span>teamID<span>,</span> teams<span>$</span>teamID<span>,</span> stringsAsFactors <span>=</span> <span>FALSE</span><span>)</span> <span>%>%</span>
  rename<span>(</span>team1 <span>=</span> Var1<span>,</span> team2 <span>=</span> Var2<span>)</span> <span>%>%</span>
  filter<span>(</span>team1 <span>!=</span> team2<span>)</span> <span>%>%</span>
  left_join<span>(</span>teams<span>,</span> by <span>=</span> <span>c</span><span>(</span>team1 <span>=</span> <span>&quot;teamID&quot;</span><span>))</span> <span>%>%</span>
  left_join<span>(</span>teams<span>,</span> by <span>=</span> <span>c</span><span>(</span>team2 <span>=</span> <span>&quot;teamID&quot;</span><span>),</span> suffix <span>=</span> <span>c</span><span>(</span><span>&quot;_1&quot;</span><span>,</span> <span>&quot;_2&quot;</span><span>))</span>
  
<span>head</span><span>(</span>teams_compared<span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##   team1 team2 W_1 L_1 P_1 W_2 L_2 P_2
## 1   ATL   ARI  68  93 42%  69  93 43%
## 2   BAL   ARI  89  73 55%  69  93 43%
## 3   BOS   ARI  93  69 57%  69  93 43%
## 4   CHA   ARI  78  84 48%  69  93 43%
## 5   CHN   ARI 103  58 64%  69  93 43%
## 6   CIN   ARI  68  94 42%  69  93 43%</code></pre>
</figure>
<p>Now, we can use our <code>pdiff</code> function to calculate the probability that team 1’s win percentage is greater than team 2’s win percentage.</p>
<figure>
<pre><code data-lang="r">teams_compared <span>&lt;-</span> teams_compared <span>%>%</span>
  mutate<span>(</span>p_1_greater <span>=</span> pmap<span>(</span><span>list</span><span>(</span>k1 <span>=</span> W_1<span>,</span> k2 <span>=</span> W_2<span>,</span> n1 <span>=</span> W_1 <span>+</span> L_1<span>,</span> n2 <span>=</span> W_2 <span>+</span> L_2<span>),</span> pdiff<span>))</span> <span>%>%</span>
  unnest<span>(</span>p_1_greater<span>)</span></code></pre>
</figure>
<p>Let’s see which teams are most likely to have a higher win percentage than other teams:</p>
<figure>
<pre><code data-lang="r">teams_compared <span>%>%</span>
  filter<span>(</span>P_1 <span>></span> P_2<span>)</span> <span>%>%</span>
  arrange<span>(</span><span>-</span>p_1_greater<span>)</span> <span>%>%</span>
  <span>head</span><span>()</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##   team1 team2 W_1 L_1 P_1 W_2 L_2 P_2 p_1_greater
## 1   CHN   MIN 103  58 64%  59 103 36%   0.9999997
## 2   TEX   MIN  95  67 59%  59 103 36%   0.9999694
## 3   WAS   MIN  95  67 59%  59 103 36%   0.9999694
## 4   CHN   CIN 103  58 64%  68  94 42%   0.9999631
## 5   CHN   SDN 103  58 64%  68  94 42%   0.9999631
## 6   CHN   TBA 103  58 64%  68  94 42%   0.9999631</code></pre>
</figure>
<p>The Cubs (CHN), Rangers, and Nationals all have extremely high posterior probabilities of having a better win percentage than the Minnesota Twins.</p>
<p>Let’s see which teams are most likely to be the same:</p>
<figure>
<pre><code data-lang="r">teams_compared <span>%>%</span>
  filter<span>(</span>P_1 <span>></span> P_2<span>)</span> <span>%>%</span>
  arrange<span>(</span>p_1_greater<span>)</span> <span>%>%</span>
  <span>head</span><span>()</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##   team1 team2 W_1 L_1 P_1 W_2 L_2 P_2 p_1_greater
## 1   TEX   CLE  95  67 59%  94  67 58%   0.5186491
## 2   WAS   CLE  95  67 59%  94  67 58%   0.5186491
## 3   NYN   DET  87  75 54%  86  75 53%   0.5206036
## 4   SFN   DET  87  75 54%  86  75 53%   0.5206036
## 5   ARI   ATL  69  93 43%  68  93 42%   0.5257282
## 6   OAK   ATL  69  93 43%  68  93 42%   0.5257282</code></pre>
</figure>
<p>Texas and Cleveland had very similar records, so the posterior probability that Texas was better than Cleveland is only 52%.</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="https://silentspring.org/devblog/2018/01/11/bayesian-binomial-test.html"> Silent Spring Institute Developer Blog</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/IsSmeHPCMTU" height="1"></a></div>