<p><strong>Previously in this series:</strong></p>
<ul><li><a href="http://varianceexplained.org/statistics/beta_distribution_and_baseball/">The beta distribution</a></li>
<li><a href="http://varianceexplained.org/r/empirical_bayes_baseball/">Empirical Bayes estimation</a></li>
<li><a href="http://varianceexplained.org/r/credible_intervals_baseball/">Credible intervals</a></li>
<li><a href="http://varianceexplained.org/r/bayesian_fdr_baseball/">The Bayesian approach to false discovery rates</a></li>
<li><a href="http://varianceexplained.org/r/bayesian_ab_baseball/">Bayesian A/B testing</a></li>
<li><a href="http://varianceexplained.org/r/beta_binomial_baseball/">Beta-binomial regression</a></li>
<li><a href="http://varianceexplained.org/r/hierarchical_bayes_baseball/">Understanding empirical Bayesian hierarchical modeling</a></li>
<li><a href="http://varianceexplained.org/r/mixture-models-baseball/">Mixture models and expectation-maximization</a></li>
</ul><p>We’ve introduced a number of statistical techniques in this series: estimating a beta prior, beta-binomial regression, hypothesis testing, mixture models, and many other components of the empirical Bayes approach. These approaches are useful whenever you have many observations of success/total data.</p>
<p>Since I’ve provided the R code within each post, you can certainly apply these methods to your own data (as <a href="http://juliasilge.com/blog/Bayesian-Blues/">some already have!</a>). However, you’d probably find yourself copying and pasting a rather large amount of code, which can take you out of the flow of your own data analysis and introduces opportunities for mistakes.</p>
<p>Here, I’ll introduce the new <a href="https://github.com/dgrtwo/ebbr">ebbr package</a> for performing empirical Bayes on binomial data. The package offers convenient tools for performing almost all the analyses we’ve done during this series, along with documentation and examples. This post also serves as review of our entire empirical Bayes series so far: we’ll touch on each post briefly and recreate some of the key results.</p>
<h3>Setup</h3>
<p>We start, as always, by assembling some per-player batting data (note that all the code in this post can be found <a href="https://github.com/dgrtwo/dgrtwo.github.com/blob/master/_R/2017-01-05-ebbr-package.Rmd">here</a>). It’s worth remembering that while we’re using batting averages as an example, this package can be applied to many types of data.</p>
<pre><code><span>library</span><span>(</span><span>Lahman</span><span>)</span><span>
</span><span>library</span><span>(</span><span>dplyr</span><span>)</span><span>
</span><span>library</span><span>(</span><span>tidyr</span><span>)</span><span>
</span><span># Grab career batting average of non-pitchers
# (allow players that have pitched &lt;= 3 games, like Ty Cobb)
</span><span>pitchers</span><span> </span><span>&lt;-</span><span> </span><span>Pitching</span><span> </span><span>%>%</span><span>
</span><span>group_by</span><span>(</span><span>playerID</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>summarize</span><span>(</span><span>gamesPitched</span><span> </span><span>=</span><span> </span><span>sum</span><span>(</span><span>G</span><span>))</span><span> </span><span>%>%</span><span>
</span><span>filter</span><span>(</span><span>gamesPitched</span><span> </span><span>></span><span> </span><span>3</span><span>)</span><span>
</span><span># Add player names
</span><span>player_names</span><span> </span><span>&lt;-</span><span> </span><span>Master</span><span> </span><span>%>%</span><span>
</span><span>tbl_df</span><span>()</span><span> </span><span>%>%</span><span>
</span><span>dplyr</span><span>::</span><span>select</span><span>(</span><span>playerID</span><span>,</span><span> </span><span>nameFirst</span><span>,</span><span> </span><span>nameLast</span><span>,</span><span> </span><span>bats</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>unite</span><span>(</span><span>name</span><span>,</span><span> </span><span>nameFirst</span><span>,</span><span> </span><span>nameLast</span><span>,</span><span> </span><span>sep</span><span> </span><span>=</span><span> </span><span>&quot; &quot;</span><span>)</span><span>
</span><span># include the &quot;bats&quot; (handedness) and &quot;year&quot; column for later
</span><span>career_full</span><span> </span><span>&lt;-</span><span> </span><span>Batting</span><span> </span><span>%>%</span><span>
</span><span>filter</span><span>(</span><span>AB</span><span> </span><span>></span><span> </span><span>0</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>anti_join</span><span>(</span><span>pitchers</span><span>,</span><span> </span><span>by</span><span> </span><span>=</span><span> </span><span>&quot;playerID&quot;</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>group_by</span><span>(</span><span>playerID</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>summarize</span><span>(</span><span>H</span><span> </span><span>=</span><span> </span><span>sum</span><span>(</span><span>H</span><span>),</span><span> </span><span>AB</span><span> </span><span>=</span><span> </span><span>sum</span><span>(</span><span>AB</span><span>),</span><span> </span><span>year</span><span> </span><span>=</span><span> </span><span>mean</span><span>(</span><span>yearID</span><span>))</span><span> </span><span>%>%</span><span>
</span><span>inner_join</span><span>(</span><span>player_names</span><span>,</span><span> </span><span>by</span><span> </span><span>=</span><span> </span><span>&quot;playerID&quot;</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>filter</span><span>(</span><span>!</span><span>is.na</span><span>(</span><span>bats</span><span>))</span><span>
</span><span># We don't need all this data for every step
</span><span>career</span><span> </span><span>&lt;-</span><span> </span><span>career_full</span><span> </span><span>%>%</span><span>
</span><span>select</span><span>(</span><span>-</span><span>bats</span><span>,</span><span> </span><span>-</span><span>year</span><span>)</span></code></pre><h2>Empirical Bayes estimation</h2>
<ul><li><a href="http://varianceexplained.org/r/empirical_bayes_baseball/">Empirical Bayes estimation</a></li>
<li><a href="http://varianceexplained.org/r/credible_intervals_baseball/">Credible intervals</a></li>
</ul><p>In one of the <a href="http://varianceexplained.org/r/empirical_bayes_baseball/">first posts in the series</a>, we noticed that the distribution of player batting averages looked roughly like a beta distribution:</p>
<pre><code><span>career</span><span> </span><span>%>%</span><span>
</span><span>filter</span><span>(</span><span>AB</span><span> </span><span>>=</span><span> </span><span>100</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>ggplot</span><span>(</span><span>aes</span><span>(</span><span>H</span><span> </span><span>/</span><span> </span><span>AB</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_histogram</span><span>()</span></code></pre><p></p>
<p>We thus wanted to estimate the beta prior for the overall dataset, which is the first step of empirical Bayes analysis. The <code>ebb_fit_prior</code> function encapsulates this, taking the data along with the success/total columns and fitting the beta through maximum likelihood.</p>
<pre><code><span>library</span><span>(</span><span>ebbr</span><span>)</span><span>
</span><span>prior</span><span> </span><span>&lt;-</span><span> </span><span>career</span><span> </span><span>%>%</span><span>
</span><span>filter</span><span>(</span><span>AB</span><span> </span><span>>=</span><span> </span><span>500</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>ebb_fit_prior</span><span>(</span><span>H</span><span>,</span><span> </span><span>AB</span><span>)</span><span>
</span><span>prior</span></code></pre><pre><code>## Empirical Bayes binomial fit with method mle 
## Parameters:
## # A tibble: 1 × 2
##   alpha  beta
##   &lt;dbl> &lt;dbl>
## 1  96.9   275</code></pre><p><code>prior</code> is an <code>ebb_prior</code> object, which is a statistical model object containing the details of the beta distribution. Here we can see the overall alpha and beta parameters. (We’ll see later that it can store more complicated models).</p>
<p>The second step of empirical Bayes analysis is updating each observation based on the overall statistical model. Based on the philosophy of the <a href="https://github.com/tidyverse/broom">broom</a> package, this is achieved with the augment function:</p>
<pre><code><span>augment</span><span>(</span><span>prior</span><span>,</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>career</span><span>)</span></code></pre><pre><code>## # A tibble: 9,548 × 10
##     playerID     H    AB           name .alpha1 .beta1 .fitted   .raw
##        &lt;chr> &lt;int> &lt;int>          &lt;chr>   &lt;dbl>  &lt;dbl>   &lt;dbl>  &lt;dbl>
## 1  aaronha01  3771 12364     Hank Aaron  3867.9   8868   0.304 0.3050
## 2  aaronto01   216   944   Tommie Aaron   312.9   1003   0.238 0.2288
## 3   abadan01     2    21      Andy Abad    98.9    294   0.252 0.0952
## 4  abadijo01    11    49    John Abadie   107.9    313   0.257 0.2245
## 5  abbated01   772  3044 Ed Abbaticchio   868.9   2547   0.254 0.2536
## 6  abbeych01   492  1751  Charlie Abbey   588.9   1534   0.277 0.2810
## 7  abbotda01     1     7     Dan Abbott    97.9    281   0.259 0.1429
## 8  abbotfr01   107   513    Fred Abbott   203.9    681   0.231 0.2086
## 9  abbotje01   157   596    Jeff Abbott   253.9    714   0.262 0.2634
## 10 abbotku01   523  2044    Kurt Abbott   619.9   1796   0.257 0.2559
## # ... with 9,538 more rows, and 2 more variables: .low &lt;dbl>, .high &lt;dbl></code></pre><p>Notice we’ve now added several columns to the original data, each beginning with <code>.</code> (which is a convention of the <code>augment</code> verb to avoid rewriting e. We have the <code>.alpha1</code> and <code>.beta1</code> columns as the parameters for each player’s posterior distribution, as well as <code>.fitted</code> representing the new posterior mean (the “shrunken average”).</p>
<p>We often want to run these two steps in sequence: estimating a model, then using it as a prior for each observation. The <code>ebbr</code> package provides a shortcut, combining them into one step with <code>add_ebb_estimate</code>:</p>
<pre><code><span>eb_career</span><span> </span><span>&lt;-</span><span> </span><span>career</span><span> </span><span>%>%</span><span>
</span><span>add_ebb_estimate</span><span>(</span><span>H</span><span>,</span><span> </span><span>AB</span><span>,</span><span> </span><span>prior_subset</span><span> </span><span>=</span><span> </span><span>AB</span><span> </span><span>>=</span><span> </span><span>500</span><span>)</span><span>
</span><span>eb_career</span></code></pre><pre><code>## # A tibble: 9,548 × 10
##     playerID     H    AB           name .alpha1 .beta1 .fitted   .raw
##        &lt;chr> &lt;int> &lt;int>          &lt;chr>   &lt;dbl>  &lt;dbl>   &lt;dbl>  &lt;dbl>
## 1  aaronha01  3771 12364     Hank Aaron  3867.9   8868   0.304 0.3050
## 2  aaronto01   216   944   Tommie Aaron   312.9   1003   0.238 0.2288
## 3   abadan01     2    21      Andy Abad    98.9    294   0.252 0.0952
## 4  abadijo01    11    49    John Abadie   107.9    313   0.257 0.2245
## 5  abbated01   772  3044 Ed Abbaticchio   868.9   2547   0.254 0.2536
## 6  abbeych01   492  1751  Charlie Abbey   588.9   1534   0.277 0.2810
## 7  abbotda01     1     7     Dan Abbott    97.9    281   0.259 0.1429
## 8  abbotfr01   107   513    Fred Abbott   203.9    681   0.231 0.2086
## 9  abbotje01   157   596    Jeff Abbott   253.9    714   0.262 0.2634
## 10 abbotku01   523  2044    Kurt Abbott   619.9   1796   0.257 0.2559
## # ... with 9,538 more rows, and 2 more variables: .low &lt;dbl>, .high &lt;dbl></code></pre><p>(The <code>add_</code> prefix is inspired by <a href="https://github.com/hadley/modelr">modelr</a>’s <code>add_residuals</code> and <code>add_predictions</code>, which also fit a statistical model then append columns). Note the <code>prior_subset</code> argument, which noted that while we wanted to keep all the shrunken values in our output, we wanted to fit the prior only on individuals with at least 500 at-bats.</p>
<h3>Estimates and credible intervals</h3>
<p>Having the posterior estimates for each player lets us explore the model results using our normal tidy tools like dplyr and ggplot2. For example, we could visualize how batting averages were shrunken towards the mean of the prior:</p>
<pre><code><span>eb_career</span><span> </span><span>%>%</span><span>
</span><span>ggplot</span><span>(</span><span>aes</span><span>(</span><span>.raw</span><span>,</span><span> </span><span>.fitted</span><span>,</span><span> </span><span>color</span><span> </span><span>=</span><span> </span><span>AB</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_point</span><span>()</span><span> </span><span>+</span><span>
</span><span>geom_abline</span><span>(</span><span>color</span><span> </span><span>=</span><span> </span><span>&quot;red&quot;</span><span>)</span><span> </span><span>+</span><span>
</span><span>scale_color_continuous</span><span>(</span><span>trans</span><span> </span><span>=</span><span> </span><span>&quot;log&quot;</span><span>,</span><span> </span><span>breaks</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>1</span><span>,</span><span> </span><span>10</span><span>,</span><span> </span><span>100</span><span>,</span><span> </span><span>1000</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>tidy</span><span>(</span><span>prior</span><span>)</span><span>$</span><span>mean</span><span>,</span><span> </span><span>color</span><span> </span><span>=</span><span> </span><span>&quot;red&quot;</span><span>,</span><span> </span><span>lty</span><span> </span><span>=</span><span> </span><span>2</span><span>)</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>&quot;Raw batting average&quot;</span><span>,</span><span>
</span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;Shrunken batting average&quot;</span><span>)</span></code></pre><p></p>
<p>This was one of our first visualizations in the <a href="http://varianceexplained.org/r/empirical_bayes_baseball/">empirical Bayes estimation post</a>. I like how it captures what empirical Bayes estimation is doing: moving all batting averages towards the prior mean (the dashed red line), but moving them less if there is a lot of information about that player (high AB).</p>
<p>In the following post, we used <a href="http://varianceexplained.org/r/credible_intervals_baseball/">credible intervals</a>, to visualize our uncertainty about each player’s true batting average. The output of <code>add_ebb_estimate</code> comes with those credible intervals in the form of <code>.low</code> and <code>.high</code>. This makes it easy to visualize these intervals for particular players, such as the 1998 Yankees:</p>
<pre><code><span>yankee_1998</span><span> </span><span>&lt;-</span><span> </span><span>c</span><span>(</span><span>&quot;brosisc01&quot;</span><span>,</span><span> </span><span>&quot;jeterde01&quot;</span><span>,</span><span> </span><span>&quot;knoblch01&quot;</span><span>,</span><span>
</span><span>&quot;martiti02&quot;</span><span>,</span><span> </span><span>&quot;posadjo01&quot;</span><span>,</span><span> </span><span>&quot;strawda01&quot;</span><span>,</span><span> </span><span>&quot;willibe02&quot;</span><span>)</span><span>
</span><span>eb_career</span><span> </span><span>%>%</span><span>
</span><span>filter</span><span>(</span><span>playerID</span><span> </span><span>%in%</span><span> </span><span>yankee_1998</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>mutate</span><span>(</span><span>name</span><span> </span><span>=</span><span> </span><span>reorder</span><span>(</span><span>name</span><span>,</span><span> </span><span>.fitted</span><span>))</span><span> </span><span>%>%</span><span>
</span><span>ggplot</span><span>(</span><span>aes</span><span>(</span><span>.fitted</span><span>,</span><span> </span><span>name</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_point</span><span>()</span><span> </span><span>+</span><span>
</span><span>geom_errorbarh</span><span>(</span><span>aes</span><span>(</span><span>xmin</span><span> </span><span>=</span><span> </span><span>.low</span><span>,</span><span> </span><span>xmax</span><span> </span><span>=</span><span> </span><span>.high</span><span>))</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>&quot;Estimated batting average (w/ 95% confidence interval)&quot;</span><span>,</span><span>
</span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;Player&quot;</span><span>)</span></code></pre><p></p>
<p>Notice that once we have the output from <code>add_ebb_estimate</code>, we’re no longer relying on the <code>ebbr</code> package, only on dplyr and ggplot2.</p>
<h2>Hierarchical modeling</h2>
<ul><li><a href="http://varianceexplained.org/r/beta_binomial_baseball/">Beta-binomial regression</a></li>
<li><a href="http://varianceexplained.org/r/hierarchical_bayes_baseball/">Understanding empirical Bayesian hierarchical modeling</a></li>
</ul><p>In two subsequent posts, we examined how this beta-binomial model may not be appropriate, because of the relationship between a player’s at-bats and their batting average. Good batters tend to have long careers, while poor batters may retire quickly.</p>
<pre><code><span>career</span><span> </span><span>%>%</span><span>
</span><span>filter</span><span>(</span><span>AB</span><span> </span><span>>=</span><span> </span><span>10</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>ggplot</span><span>(</span><span>aes</span><span>(</span><span>AB</span><span>,</span><span> </span><span>H</span><span> </span><span>/</span><span> </span><span>AB</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_point</span><span>()</span><span> </span><span>+</span><span>
</span><span>geom_smooth</span><span>(</span><span>method</span><span> </span><span>=</span><span> </span><span>&quot;lm&quot;</span><span>)</span><span> </span><span>+</span><span>
</span><span>scale_x_log10</span><span>()</span></code></pre><p></p>
<p>We solved this by fitting a prior that depended on AB, through <a href="http://varianceexplained.org/r/beta_binomial_baseball/">beta-binomial regression</a>. The <code>add_ebb_estimate</code> function offers this option, by setting <code>method = &quot;gamlss&quot;</code> and providing a formula to <code>mu_predictors</code>.</p>
<pre><code><span>eb_career_ab</span><span> </span><span>&lt;-</span><span> </span><span>career</span><span> </span><span>%>%</span><span>
</span><span>add_ebb_estimate</span><span>(</span><span>H</span><span>,</span><span> </span><span>AB</span><span>,</span><span> </span><span>method</span><span> </span><span>=</span><span> </span><span>&quot;gamlss&quot;</span><span>,</span><span>
</span><span>mu_predictors</span><span> </span><span>=</span><span> </span><span>~</span><span> </span><span>log10</span><span>(</span><span>AB</span><span>))</span><span>
</span><span>eb_career_ab</span></code></pre><pre><code>## # A tibble: 9,548 × 14
##     playerID     H    AB           name   .mu  .sigma .alpha0 .beta0
##        &lt;chr> &lt;int> &lt;int>          &lt;chr> &lt;dbl>   &lt;dbl>   &lt;dbl>  &lt;dbl>
## 1  aaronha01  3771 12364     Hank Aaron 0.289 0.00181   159.6    393
## 2  aaronto01   216   944   Tommie Aaron 0.247 0.00181   136.5    416
## 3   abadan01     2    21      Andy Abad 0.193 0.00181   106.7    446
## 4  abadijo01    11    49    John Abadie 0.204 0.00181   112.9    440
## 5  abbated01   772  3044 Ed Abbaticchio 0.265 0.00181   146.7    406
## 6  abbeych01   492  1751  Charlie Abbey 0.257 0.00181   141.8    411
## 7  abbotda01     1     7     Dan Abbott 0.179 0.00181    99.1    454
## 8  abbotfr01   107   513    Fred Abbott 0.238 0.00181   131.4    421
## 9  abbotje01   157   596    Jeff Abbott 0.240 0.00181   132.6    420
## 10 abbotku01   523  2044    Kurt Abbott 0.259 0.00181   143.2    410
## # ... with 9,538 more rows, and 6 more variables: .alpha1 &lt;dbl>,
## #   .beta1 &lt;dbl>, .fitted &lt;dbl>, .raw &lt;dbl>, .low &lt;dbl>, .high &lt;dbl></code></pre><p>(You can also provide <code>sigma_predictors</code> to have the variance of the beta depend on regressors, though it’s less common that you’d want to do so).</p>
<p>The augmented output is now a bit more complicated: besides the posterior parameters such as <code>.alpha1</code>, <code>.beta1</code>, and <code>.fitted</code>, each observation also has its own prior parameters <code>.alpha0</code> and <code>.beta0</code>. These are predicted based on the regression on <code>AB</code>.</p>
<p>The other parameters, such as <code>.fitted</code> and the credible interval, are now shrinking towards a trend rather than towards a constant. We can see this by plotting <code>AB</code> against the original and the shrunken estimates:</p>
<pre><code><span>eb_career_ab</span><span> </span><span>%>%</span><span>
</span><span>filter</span><span>(</span><span>AB</span><span> </span><span>></span><span> </span><span>10</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>rename</span><span>(</span><span>Raw</span><span> </span><span>=</span><span> </span><span>.raw</span><span>,</span><span> </span><span>Shrunken</span><span> </span><span>=</span><span> </span><span>.fitted</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>gather</span><span>(</span><span>type</span><span>,</span><span> </span><span>estimate</span><span>,</span><span> </span><span>Raw</span><span>,</span><span> </span><span>Shrunken</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>ggplot</span><span>(</span><span>aes</span><span>(</span><span>AB</span><span>,</span><span> </span><span>estimate</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_point</span><span>()</span><span> </span><span>+</span><span>
</span><span>facet_wrap</span><span>(</span><span>~</span><span> </span><span>type</span><span>)</span><span> </span><span>+</span><span>
</span><span>scale_x_log10</span><span>()</span></code></pre><p></p>
<p>As we saw in the post on <a href="http://varianceexplained.org/r/hierarchical_bayes_baseball/">hierarchical modeling</a>, the model can incorporate still more useful information. For example, it could include what year they play in, and whether they are left- or right- handed, both of which tend to affect batting average.</p>
<pre><code><span>library</span><span>(</span><span>splines</span><span>)</span><span>
</span><span>eb_career_prior</span><span> </span><span>&lt;-</span><span> </span><span>career_full</span><span> </span><span>%>%</span><span>
</span><span>ebb_fit_prior</span><span>(</span><span>H</span><span>,</span><span> </span><span>AB</span><span>,</span><span> </span><span>method</span><span> </span><span>=</span><span> </span><span>&quot;gamlss&quot;</span><span>,</span><span>
</span><span>mu_predictors</span><span> </span><span>=</span><span> </span><span>~</span><span> </span><span>0</span><span> </span><span>+</span><span> </span><span>ns</span><span>(</span><span>year</span><span>,</span><span> </span><span>df</span><span> </span><span>=</span><span> </span><span>5</span><span>)</span><span> </span><span>*</span><span> </span><span>bats</span><span> </span><span>+</span><span> </span><span>log</span><span>(</span><span>AB</span><span>))</span></code></pre><p>In this case I’m fitting the prior with <code>ebb_fit_prior</code> rather than adding the estimates with <code>add_ebb_estimate</code>. This lets us feed it new data that we generate ourselves, and examine how the posterior distribution would change. This takes a bit more work, but lets us re-generate one of the more interesting plots from that post about how time and handedness relate:</p>
<pre><code><span># fake data ranging from 1885 to 2013
</span><span>fake_data</span><span> </span><span>&lt;-</span><span> </span><span>crossing</span><span>(</span><span>H</span><span> </span><span>=</span><span> </span><span>300</span><span>,</span><span>
</span><span>AB</span><span> </span><span>=</span><span> </span><span>1000</span><span>,</span><span>
</span><span>year</span><span> </span><span>=</span><span> </span><span>seq</span><span>(</span><span>1885</span><span>,</span><span> </span><span>2013</span><span>),</span><span>
</span><span>bats</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>&quot;L&quot;</span><span>,</span><span> </span><span>&quot;R&quot;</span><span>))</span><span>
</span><span># find the mean of the prior, as well as the 95% quantiles,
# for each of these combinations. This does require a bit of
# manual manipulation of alpha0 and beta0:
</span><span>augment</span><span>(</span><span>eb_career_prior</span><span>,</span><span> </span><span>newdata</span><span> </span><span>=</span><span> </span><span>fake_data</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>mutate</span><span>(</span><span>prior</span><span> </span><span>=</span><span> </span><span>.alpha0</span><span> </span><span>/</span><span> </span><span>(</span><span>.alpha0</span><span> </span><span>+</span><span> </span><span>.beta0</span><span>),</span><span>
</span><span>prior.low</span><span> </span><span>=</span><span> </span><span>qbeta</span><span>(</span><span>.025</span><span>,</span><span> </span><span>.alpha0</span><span>,</span><span> </span><span>.beta0</span><span>),</span><span>
</span><span>prior.high</span><span> </span><span>=</span><span> </span><span>qbeta</span><span>(</span><span>.975</span><span>,</span><span> </span><span>.alpha0</span><span>,</span><span> </span><span>.beta0</span><span>))</span><span> </span><span>%>%</span><span>
</span><span>ggplot</span><span>(</span><span>aes</span><span>(</span><span>year</span><span>,</span><span> </span><span>prior</span><span>,</span><span> </span><span>color</span><span> </span><span>=</span><span> </span><span>bats</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_line</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>prior.low</span><span>,</span><span> </span><span>ymax</span><span> </span><span>=</span><span> </span><span>prior.high</span><span>),</span><span> </span><span>alpha</span><span> </span><span>=</span><span> </span><span>.1</span><span>,</span><span> </span><span>lty</span><span> </span><span>=</span><span> </span><span>2</span><span>)</span><span> </span><span>+</span><span>
</span><span>ylab</span><span>(</span><span>&quot;Prior distribution (mean + 95% quantiles)&quot;</span><span>)</span></code></pre><p></p>
<p>Since the <code>ebbr</code> package makes these models convenient to fit, you can try a few variations on the model and compare them.</p>
<h2>Hypothesis testing</h2>
<ul><li><a href="http://varianceexplained.org/r/bayesian_ab_baseball/">Bayesian A/B testing</a></li>
<li><a href="http://varianceexplained.org/r/beta_binomial_baseball/">Beta-binomial regression</a></li>
</ul><p>Another pair of posts examined the problem of hypothesis testing. For example, we wanted to get a posterior probability for the statement “this player’s true batting average is greater than .300”, so that we could construct a “Hall of Fame” of such players.</p>
<p>This method is implemented in the <code>add_ebb_prop_test</code> (notice that like <code>add_ebb_estimate</code>, it adds columns to existing data). <code>add_ebb_prop_test</code> takes the output of an earlier <code>add_ebb_estimate</code> operation, which contains posterior parameters for each observation, and appends columns to it:</p>
<pre><code><span>test_300</span><span> </span><span>&lt;-</span><span> </span><span>career</span><span> </span><span>%>%</span><span>
</span><span>add_ebb_estimate</span><span>(</span><span>H</span><span>,</span><span> </span><span>AB</span><span>,</span><span> </span><span>method</span><span> </span><span>=</span><span> </span><span>&quot;gamlss&quot;</span><span>,</span><span> </span><span>mu_predictors</span><span> </span><span>=</span><span> </span><span>~</span><span> </span><span>log10</span><span>(</span><span>AB</span><span>))</span><span> </span><span>%>%</span><span>
</span><span>add_ebb_prop_test</span><span>(</span><span>.300</span><span>,</span><span> </span><span>sort</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>)</span><span>
</span><span>test_300</span></code></pre><pre><code>## # A tibble: 9,548 × 16
##     playerID     H    AB                 name   .mu  .sigma .alpha0 .beta0
##        &lt;chr> &lt;int> &lt;int>                &lt;chr> &lt;dbl>   &lt;dbl>   &lt;dbl>  &lt;dbl>
## 1   cobbty01  4189 11434              Ty Cobb 0.287 0.00181     159    394
## 2  hornsro01  2930  8173       Rogers Hornsby 0.282 0.00181     156    397
## 3  speaktr01  3514 10195         Tris Speaker 0.285 0.00181     158    395
## 4  delahed01  2596  7505         Ed Delahanty 0.280 0.00181     155    398
## 5  willite01  2654  7706         Ted Williams 0.281 0.00181     155    398
## 6  keelewi01  2932  8591        Willie Keeler 0.283 0.00181     156    397
## 7  lajoina01  3242  9589           Nap Lajoie 0.284 0.00181     157    396
## 8  jacksjo01  1772  4981 Shoeless Joe Jackson 0.274 0.00181     151    402
## 9  gwynnto01  3141  9288           Tony Gwynn 0.284 0.00181     157    396
## 10 heilmha01  2660  7787       Harry Heilmann 0.281 0.00181     155    397
## # ... with 9,538 more rows, and 8 more variables: .alpha1 &lt;dbl>,
## #   .beta1 &lt;dbl>, .fitted &lt;dbl>, .raw &lt;dbl>, .low &lt;dbl>, .high &lt;dbl>,
## #   .pep &lt;dbl>, .qvalue &lt;dbl></code></pre><p>(Note the <code>sort = TRUE</code> argument, which sorts in order of our confidence in each player). There are now too many columns to read easily, so we’ll select only a few of the most interesting ones:</p>
<pre><code><span>test_300</span><span> </span><span>%>%</span><span>
</span><span>select</span><span>(</span><span>name</span><span>,</span><span> </span><span>H</span><span>,</span><span> </span><span>AB</span><span>,</span><span> </span><span>.fitted</span><span>,</span><span> </span><span>.low</span><span>,</span><span> </span><span>.high</span><span>,</span><span> </span><span>.pep</span><span>,</span><span> </span><span>.qvalue</span><span>)</span></code></pre><pre><code>## # A tibble: 9,548 × 8
##                    name     H    AB .fitted  .low .high     .pep  .qvalue
##                   &lt;chr> &lt;int> &lt;int>   &lt;dbl> &lt;dbl> &lt;dbl>    &lt;dbl>    &lt;dbl>
## 1               Ty Cobb  4189 11434   0.363 0.354 0.371 2.41e-49 2.41e-49
## 2        Rogers Hornsby  2930  8173   0.354 0.344 0.364 2.57e-27 1.28e-27
## 3          Tris Speaker  3514 10195   0.342 0.333 0.351 6.91e-21 2.30e-21
## 4          Ed Delahanty  2596  7505   0.341 0.331 0.352 5.77e-16 1.44e-16
## 5          Ted Williams  2654  7706   0.340 0.330 0.350 1.84e-15 4.84e-16
## 6         Willie Keeler  2932  8591   0.338 0.328 0.347 3.50e-15 9.86e-16
## 7            Nap Lajoie  3242  9589   0.335 0.326 0.344 1.04e-14 2.32e-15
## 8  Shoeless Joe Jackson  1772  4981   0.348 0.335 0.360 1.40e-14 3.78e-15
## 9            Tony Gwynn  3141  9288   0.335 0.326 0.344 2.68e-14 6.34e-15
## 10       Harry Heilmann  2660  7787   0.338 0.327 0.348 6.70e-14 1.24e-14
## # ... with 9,538 more rows</code></pre><p>Notice that two columns have been added, with per-player metrics <a href="http://varianceexplained.org/r/bayesian_fdr_baseball/">described in this post</a>.</p>
<ul><li><code>.pep</code>: the posterior error probability- the probability that this player’s true batting average is less than .3.</li>
<li><code>.qvalue</code>: the q-value, which corrects for multiple testing by controlling for false discovery rate (FDR). Allowing players with a q-value below .05 would mean only 5% of the ones included would be false discoveries.</li>
</ul><p>For example, we could find how many players would be added to our “Hall of Fame” with an FDR of 5%, or 1%:</p>
<pre><code><span>sum</span><span>(</span><span>test_300</span><span>$</span><span>.qvalue</span><span> </span><span>&lt;</span><span> </span><span>.05</span><span>)</span></code></pre><pre><code>## [1] 115</code></pre><pre><code><span>sum</span><span>(</span><span>test_300</span><span>$</span><span>.qvalue</span><span> </span><span>&lt;</span><span> </span><span>.01</span><span>)</span></code></pre><pre><code>## [1] 79</code></pre><h3>Player-player A/B test</h3>
<p><a href="http://varianceexplained.org/r/bayesian_ab_baseball/">Another post</a> discussed the case where instead of comparing each observation to a single threshold (like .300) we want to compare to another player’s posterior distribution. I noted that this is similar to the problem of “A/B testing”, where we might be comparing two clickthrough rates, each represented by successes / total.</p>
<p>The post compared each player in history to Mike Piazza, and found players we were confident were better batters. We’d first find Piazza’s posterior parameters:</p>
<pre><code><span>piazza</span><span> </span><span>&lt;-</span><span> </span><span>eb_career_ab</span><span> </span><span>%>%</span><span>
</span><span>filter</span><span>(</span><span>name</span><span> </span><span>==</span><span> </span><span>&quot;Mike Piazza&quot;</span><span>)</span><span>
</span><span>piazza_params</span><span> </span><span>&lt;-</span><span> </span><span>c</span><span>(</span><span>piazza</span><span>$</span><span>.alpha1</span><span>,</span><span> </span><span>piazza</span><span>$</span><span>.beta1</span><span>)</span><span>
</span><span>piazza_params</span></code></pre><pre><code>## [1] 2281 5183</code></pre><p>This vector of two parameters, an alpha and a beta, can be passed into <code>add_ebb_prop_test</code> just like we passed in a threshold:<a href="http://varianceexplained.org/#fn:prop">1</a></p>
<pre><code><span>compare_piazza</span><span> </span><span>&lt;-</span><span> </span><span>eb_career_ab</span><span> </span><span>%>%</span><span>
</span><span>add_ebb_prop_test</span><span>(</span><span>piazza_params</span><span>,</span><span> </span><span>approx</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>,</span><span> </span><span>sort</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>)</span></code></pre><p>Again we select only a few interesting columns:</p>
<pre><code><span>compare_piazza</span><span> </span><span>%>%</span><span>
</span><span>select</span><span>(</span><span>name</span><span>,</span><span> </span><span>H</span><span>,</span><span> </span><span>AB</span><span>,</span><span> </span><span>.fitted</span><span>,</span><span> </span><span>.low</span><span>,</span><span> </span><span>.high</span><span>,</span><span> </span><span>.pep</span><span>,</span><span> </span><span>.qvalue</span><span>)</span></code></pre><pre><code>## # A tibble: 9,548 × 8
##                    name     H    AB .fitted  .low .high     .pep  .qvalue
##                   &lt;chr> &lt;int> &lt;int>   &lt;dbl> &lt;dbl> &lt;dbl>    &lt;dbl>    &lt;dbl>
## 1               Ty Cobb  4189 11434   0.363 0.354 0.371 7.00e-17 7.00e-17
## 2        Rogers Hornsby  2930  8173   0.354 0.344 0.364 4.16e-11 2.08e-11
## 3          Tris Speaker  3514 10195   0.342 0.333 0.351 1.49e-07 4.98e-08
## 4  Shoeless Joe Jackson  1772  4981   0.348 0.335 0.360 2.45e-07 9.87e-08
## 5          Ed Delahanty  2596  7505   0.341 0.331 0.352 9.41e-07 2.67e-07
## 6          Ted Williams  2654  7706   0.340 0.330 0.350 1.84e-06 5.30e-07
## 7         Willie Keeler  2932  8591   0.338 0.328 0.347 5.06e-06 1.18e-06
## 8        Harry Heilmann  2660  7787   0.338 0.327 0.348 8.64e-06 2.11e-06
## 9        Billy Hamilton  2158  6268   0.339 0.328 0.350 1.09e-05 3.09e-06
## 10           Nap Lajoie  3242  9589   0.335 0.326 0.344 1.59e-05 4.37e-06
## # ... with 9,538 more rows</code></pre><p>Just like the one-sample test, we’ve added <code>.pep</code> and <code>.qvalue</code> columns. From this we can see a few players who we’re extremely confident are better than Piazza.</p>
<h3>Mixture models</h3>
<ul><li><a href="http://varianceexplained.org/r/mixture-models-baseball/">Mixture models and expectation-maximization</a></li>
</ul><p>This brings us to my recent post on mixture models, where we noticed that when pitchers are included, the data looks a lot less like a beta distribution and more like a combination of two betas.</p>
<pre><code><span>career_w_pitchers</span><span> </span><span>&lt;-</span><span> </span><span>Batting</span><span> </span><span>%>%</span><span>
</span><span>filter</span><span>(</span><span>AB</span><span> </span><span>>=</span><span> </span><span>25</span><span>,</span><span> </span><span>lgID</span><span> </span><span>==</span><span> </span><span>&quot;NL&quot;</span><span>,</span><span> </span><span>yearID</span><span> </span><span>>=</span><span> </span><span>1980</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>group_by</span><span>(</span><span>playerID</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>summarize</span><span>(</span><span>H</span><span> </span><span>=</span><span> </span><span>sum</span><span>(</span><span>H</span><span>),</span><span> </span><span>AB</span><span> </span><span>=</span><span> </span><span>sum</span><span>(</span><span>AB</span><span>),</span><span> </span><span>year</span><span> </span><span>=</span><span> </span><span>mean</span><span>(</span><span>yearID</span><span>))</span><span> </span><span>%>%</span><span>
</span><span>mutate</span><span>(</span><span>isPitcher</span><span> </span><span>=</span><span> </span><span>playerID</span><span> </span><span>%in%</span><span> </span><span>pitchers</span><span>$</span><span>playerID</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>inner_join</span><span>(</span><span>player_names</span><span>,</span><span> </span><span>by</span><span> </span><span>=</span><span> </span><span>&quot;playerID&quot;</span><span>)</span></code></pre><pre><code><span>ggplot</span><span>(</span><span>career_w_pitchers</span><span>,</span><span> </span><span>aes</span><span>(</span><span>H</span><span> </span><span>/</span><span> </span><span>AB</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_histogram</span><span>()</span></code></pre><p></p>
<p>Fitting a mixture model, to separate out the two beta distributions so they could be shrunken separately, took a solid amount of code.</p>
<p>The <code>ebbr</code> package thus provides tools for fitting a mixture model using an iterative expectation-maximization algorithm, with the <code>ebb_fit_mixture</code> function. Like the other estimation functions, it takes a table as the first argument, followed by two arguments for the “successes” column and the “total” column:</p>
<pre><code><span>set.seed</span><span>(</span><span>2017</span><span>)</span><span>
</span><span>mm</span><span> </span><span>&lt;-</span><span> </span><span>ebb_fit_mixture</span><span>(</span><span>career_w_pitchers</span><span>,</span><span> </span><span>H</span><span>,</span><span> </span><span>AB</span><span>,</span><span> </span><span>clusters</span><span> </span><span>=</span><span> </span><span>2</span><span>)</span></code></pre><p>It returns the parameters of two (or more) beta distributions:</p>
<pre><code><span>tidy</span><span>(</span><span>mm</span><span>)</span></code></pre><pre><code>## # A tibble: 2 × 6
##   cluster alpha  beta  mean number probability
##     &lt;chr> &lt;dbl> &lt;dbl> &lt;dbl>  &lt;int>       &lt;dbl>
## 1       2 159.2   455 0.259   1981         0.5
## 2       1  26.7   158 0.144    910         0.5</code></pre><p>It also assigns each observation to the most likely cluster. Here, we can see that cluster 1 is made up of pitchers, while cluster 2 is the non-pitchers:</p>
<pre><code><span>ggplot</span><span>(</span><span>mm</span><span>$</span><span>assignments</span><span>,</span><span> </span><span>aes</span><span>(</span><span>H</span><span> </span><span>/</span><span> </span><span>AB</span><span>,</span><span> </span><span>fill</span><span> </span><span>=</span><span> </span><span>.cluster</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_histogram</span><span>(</span><span>alpha</span><span> </span><span>=</span><span> </span><span>0.8</span><span>,</span><span> </span><span>position</span><span> </span><span>=</span><span> </span><span>&quot;identity&quot;</span><span>)</span></code></pre><p></p>
<h3>Mixture model across iterations</h3>
<p>You may be interested in how the mixture model converged to its parameters. The <code>iterations</code> component of the <code>ebb_mixture</code> object contains details on the iterations, which can be visualized.</p>
<pre><code><span>fits</span><span> </span><span>&lt;-</span><span> </span><span>mm</span><span>$</span><span>iterations</span><span>$</span><span>fits</span><span>
</span><span>fits</span></code></pre><pre><code>## # A tibble: 6 × 7
##   iteration cluster alpha  beta  mean number probability
##       &lt;int>   &lt;chr> &lt;dbl> &lt;dbl> &lt;dbl>  &lt;int>       &lt;dbl>
## 1         1       2  16.2  55.1 0.228   1456         0.5
## 2         1       1  15.7  53.6 0.226   1435         0.5
## 3         2       2 156.3 446.5 0.259   1956         0.5
## 4         2       1  28.7 169.2 0.145    935         0.5
## 5         3       2 159.2 454.6 0.259   1981         0.5
## 6         3       1  26.7 158.0 0.144    910         0.5</code></pre><pre><code><span>fits</span><span> </span><span>%>%</span><span>
</span><span>gather</span><span>(</span><span>parameter</span><span>,</span><span> </span><span>value</span><span>,</span><span> </span><span>alpha</span><span>,</span><span> </span><span>beta</span><span>,</span><span> </span><span>mean</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>ggplot</span><span>(</span><span>aes</span><span>(</span><span>iteration</span><span>,</span><span> </span><span>value</span><span>,</span><span> </span><span>color</span><span> </span><span>=</span><span> </span><span>parameter</span><span>,</span><span> </span><span>lty</span><span> </span><span>=</span><span> </span><span>cluster</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_line</span><span>()</span><span> </span><span>+</span><span>
</span><span>facet_wrap</span><span>(</span><span>~</span><span> </span><span>parameter</span><span>,</span><span> </span><span>scales</span><span> </span><span>=</span><span> </span><span>&quot;free_y&quot;</span><span>)</span></code></pre><p></p>
<p>Note that it took only about one full iteration for the parameters to get pretty close to their eventual values. We can also examine the change in cluster assignments for each observation:</p>
<pre><code><span>assignments</span><span> </span><span>&lt;-</span><span> </span><span>mm</span><span>$</span><span>iterations</span><span>$</span><span>assignments</span><span>
</span><span>assignments</span></code></pre><pre><code>## # A tibble: 11,564 × 10
##    iteration  playerID     H    AB  year isPitcher               name
##        &lt;int>     &lt;chr> &lt;int> &lt;int> &lt;dbl>     &lt;lgl>              &lt;chr>
## 1          1 abbotje01    11    42  2001     FALSE        Jeff Abbott
## 2          1 abbotku01   473  1851  1997     FALSE        Kurt Abbott
## 3          1 abbotky01     2    29  1992      TRUE        Kyle Abbott
## 4          1 abercre01    86   386  2007     FALSE Reggie Abercrombie
## 5          1 abnersh01   110   531  1989     FALSE        Shawn Abner
## 6          1 abreubo01  1602  5373  2003     FALSE        Bobby Abreu
## 7          1 abreuto01   127   497  2010     FALSE         Tony Abreu
## 8          1 acevejo01     6    77  2002      TRUE       Jose Acevedo
## 9          1 ackerji01     3    28  1986      TRUE          Jim Acker
## 10         1 adamsma01   257   909  2013     FALSE         Matt Adams
## # ... with 11,554 more rows, and 3 more variables: bats &lt;fctr>,
## #   .cluster &lt;chr>, .likelihood &lt;dbl></code></pre><pre><code><span>assignments</span><span> </span><span>%>%</span><span>
</span><span>ggplot</span><span>(</span><span>aes</span><span>(</span><span>H</span><span> </span><span>/</span><span> </span><span>AB</span><span>,</span><span> </span><span>fill</span><span> </span><span>=</span><span> </span><span>.cluster</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_histogram</span><span>(</span><span>alpha</span><span> </span><span>=</span><span> </span><span>0.8</span><span>,</span><span> </span><span>position</span><span> </span><span>=</span><span> </span><span>&quot;identity&quot;</span><span>)</span><span> </span><span>+</span><span>
</span><span>facet_wrap</span><span>(</span><span>~</span><span> </span><span>iteration</span><span>)</span></code></pre><p></p>
<p>The package’s functions for mixture modeling are a bit more primitive than the rest: I’m still working out the right output format, and some of the details are likely to change (this is one of the reasons I haven’t yet submitted <code>ebbr</code> to CRAN). Still, the package makes it easier to experiment with expectation-maximization algorithms.</p>
<h2>What’s Next: Simulation</h2>
<p>We’re approaching the end of this series on empirical Bayesian methods, and have touched on many statistical approaches for analyzing binomial data, all with the goal of estimating the “true” batting average of each player. There’s one question we haven’t answered, though: <strong>do these methods actually work?</strong></p>
<p>Even if we assume each player has a “true” batting average as our model suggests, we don’t know it, so we can’t see if our methods estimated it accurately. We believe that empirical Bayes shrinkage gets closer to the true probabilities than raw batting averages do, but we can’t measure its bias or mean-squared error. We can’t see whether credible intervals actually contain the true parameters, or whether our FDR control was successful in finding a set of players above a cutoff. This means we can’t test our methods, or examine when they work well and when they don’t.</p>
<p>In the next post, we’ll <strong>simulate</strong> some fake batting average data, which will let us know the true probabilities for each player, then examine how close our statistical methods got. Simulation is a universally useful way to test a statistical method, to build intuition about its mathematical properies, and to gain confidence that we can trust its results. This will be especially easy now that we have the <code>ebbr</code> package to encapsulate the methods themselves, and will be a good chance to demonstrate the tidyverse approach to simulation.</p>
<div>
<ol><li>
<p>This is rather similar to how the built-in <a href="https://stat.ethz.ch/R-manual/R-devel/library/stats/html/prop.test.html">prop.test</a> function handles the difference between one-sample and two-sample tests. I’ve often found this kind of behavior annoying (determining what test to use based on the length of the input), but I must admit it is convenient. <a href="http://varianceexplained.org/#fnref:prop">↩</a></p>
</li>
</ol></div>