<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/simulation-of-empirical-bayesian-methods-using-baseball-statistics/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://varianceexplained.org/r/simulation-bayes-baseball/"> Variance Explained</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><strong>Previously in this series:</strong></p>
<ul>
<li><a target="_blank" rel="nofollow" href="http://varianceexplained.org/statistics/beta_distribution_and_baseball/">The beta distribution</a></li>
<li><a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/empirical_bayes_baseball/">Empirical Bayes estimation</a></li>
<li><a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/credible_intervals_baseball/">Credible intervals</a></li>
<li><a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/bayesian_fdr_baseball/">The Bayesian approach to false discovery rates</a></li>
<li><a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/bayesian_ab_baseball/">Bayesian A/B testing</a></li>
<li><a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/beta_binomial_baseball/">Beta-binomial regression</a></li>
<li><a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/hierarchical_bayes_baseball/">Understanding empirical Bayesian hierarchical modeling</a></li>
<li><a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/mixture-models-baseball/">Mixture models and expectation-maximization</a></li>
<li><a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/ebbr-package/">The ebbr package</a></li>
</ul>
<p>We’re approaching the end of this series on empirical Bayesian methods, and have touched on many statistical approaches for analyzing binomial (success / total) 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 <em>know</em> it, so we can’t see if our methods estimated it accurately. For example, we think that empirical Bayes shrinkage gets closer to the true probabilities than raw batting averages do, but we can’t actually measure the mean-squared error. This means we can’t test our methods, or examine when they work well and when they don’t.</p>
<p>In this 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 get to the true solution. 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. In particular, this post demonstrates the tidyverse approach to simulation, which takes advantage of the <code>dplyr</code>, <code>tidyr</code>, <code>purrr</code> and <code>broom</code> packages to examine many combinations of input parameters.</p>
<h2>Setup</h2>
<p>Most of our posts started by assembling some per-player batting data. We’re going to be simulating (i.e. making up) our data for this analysis, so you might think we don’t need to look at real data at all. However, data is still necessary to estimating the parameters we’ll use in the simulation, which keeps the experiment realistic and ensures that our conclusions will be useful.</p>
<p>(Note that all the code in this post can be found <a target="_blank" rel="nofollow" href="https://github.com/dgrtwo/dgrtwo.github.com/blob/master/_R/2017-01-10-simulation-bayes-baseball.Rmd">here</a>).</p>
<figure>
<pre><code data-lang="r"><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># include the &quot;bats&quot; (handedness) and &quot;year&quot; column for later
</span><span>career</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></code></pre>
</figure>
<h3>Choosing a distribution of p and AB</h3>
<p>In the beta-binomial model we’ve been using for most of these posts, there are two values for each player :</p>
<p></p>
<p></p>
<p> are “hyperparameters”: two unobserved values that describe the entire distribution.  is the true batting average for each player- we don’t observe this, but it’s the “right answer” for each batter that we’re trying to estimate.  is the number of at-bats the player had, which <em>is</em> observed. (You might recall we had a more complicated model in the <a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/beta_binomial_baseball/">beta-binomial regression post</a> that had  depend on : we’ll get back to that).</p>
<p>Our approach is going to be to pick some “true” , then simulate  for each player. Since we’re just picking any  to start with, we may as well estimate them from our data, since we know those are plausible values (though if we wanted to be more thorough, we could try a few other values and see how our accuracy changes).</p>
<p>To do this estimation, we can use our new <a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/ebbr-package/">ebbr</a> package to fit the empirical Bayes prior.</p>
<figure>
<pre><code data-lang="r"><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>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>
</figure>
<figure>
<pre><code data-lang="text">## Empirical Bayes binomial fit with method mle 
## Parameters:
## # A tibble: 1 × 2
##   alpha  beta
##   &lt;dbl> &lt;dbl>
## 1  72.1   215</code></pre>
</figure>
<p>These two hyperparameters are all we need to simulate a few thousand values of , using the <code>rbeta</code> function:</p>
<figure>
<pre><code data-lang="r"><span>alpha0</span><span> </span><span>&lt;-</span><span> </span><span>tidy</span><span>(</span><span>prior</span><span>)</span><span>$</span><span>alpha</span><span>
</span><span>beta0</span><span> </span><span>&lt;-</span><span> </span><span>tidy</span><span>(</span><span>prior</span><span>)</span><span>$</span><span>beta</span><span>
</span><span>qplot</span><span>(</span><span>rbeta</span><span>(</span><span>10000</span><span>,</span><span> </span><span>alpha0</span><span>,</span><span> </span><span>beta0</span><span>))</span></code></pre>
</figure>
<p><img alt="center" src="https://i1.wp.com/varianceexplained.org/figs/2017-01-11-simulation-bayes-baseball/alpha0-1.png?w=456" data-recalc-dims="1"></p>
<p>There’s another component to this model: , the distribution of the number of at-bats. This is a much more unusual distribution:</p>
<figure>
<pre><code data-lang="r"><span>ggplot</span><span>(</span><span>career</span><span>,</span><span> </span><span>aes</span><span>(</span><span>AB</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_histogram</span><span>()</span><span> </span><span>+</span><span>
</span><span>scale_x_log10</span><span>()</span></code></pre>
</figure>
<p><img alt="center" src="https://i2.wp.com/varianceexplained.org/figs/2017-01-11-simulation-bayes-baseball/unnamed-chunk-3-1.png?w=456" data-recalc-dims="1"></p>
<p>The good news is, we don’t <em>need</em> to simulate these  values, since we’re not trying to estimate them with empirical Bayes. We can just use the observed values we have! (In a different study, we may be interested in how the success of empirical Bayes depends on the distribution of the s).</p>
<p>Thus, to recap, we will:</p>
<ul>
<li><strong>Estimate</strong> , which works because the parameters are not observed, but there are only a few and we can predict them with confidence.</li>
<li><strong>Simulate</strong> , based on a beta distribution, so that we can test our ability to estimate them.</li>
<li><strong>Use observed</strong> , since we know the true values and we might as well.</li>
</ul>
<h2>Shrinkage on simulated data</h2>
<p>The beta-binomial model is easy to simulate, with applications of the <code>rbeta</code> and <code>rbinom</code> functions:</p>
<figure>
<pre><code data-lang="r"><span># always set a seed when simulating
</span><span>set.seed</span><span>(</span><span>2017</span><span>)</span><span>
</span><span>career_sim</span><span> </span><span>&lt;-</span><span> </span><span>career</span><span> </span><span>%>%</span><span>
</span><span>mutate</span><span>(</span><span>p</span><span> </span><span>=</span><span> </span><span>rbeta</span><span>(</span><span>n</span><span>(),</span><span> </span><span>alpha0</span><span>,</span><span> </span><span>beta0</span><span>),</span><span>
</span><span>H</span><span> </span><span>=</span><span> </span><span>rbinom</span><span>(</span><span>n</span><span>(),</span><span> </span><span>AB</span><span>,</span><span> </span><span>p</span><span>))</span><span>
</span><span>career_sim</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## # A tibble: 10,388 × 4
##     playerID     H    AB     p
##        &lt;chr> &lt;int> &lt;int> &lt;dbl>
## 1  aaronha01  3817 12364 0.299
## 2  aaronto01   240   944 0.249
## 3   abadan01     5    21 0.274
## 4  abadijo01     9    49 0.198
## 5  abbated01   755  3044 0.249
## 6  abbeych01   449  1751 0.265
## 7  abbotda01     2     7 0.191
## 8  abbotfr01   138   513 0.251
## 9  abbotje01   116   596 0.244
## 10 abbotku01   570  2044 0.261
## # ... with 10,378 more rows</code></pre>
</figure>
<p>Just like that, we’ve generated a “true”  for each player, and then a new value of  based on it. (This means there is no relationship between how good a particular player is in our simulation and how good they are in reality).</p>
<p>Our working theory has been that our raw  estimates have had a large amount of noise when  is small, and that empirical Bayes helps moderate it. Now, since we know the true value of  for each player, we can finally examine whether that’s true: and we can see what the empirical Bayes method is giving up as well.</p>
<p>Let’s visualize the true values of  versus the estimates, which we’ll call , using either raw estimation or empirical Bayes shrinkage. (Again, we couldn’t have made this plot with the real data since we don’t know the true : it’s possible only in our simulation).</p>
<figure>
<pre><code data-lang="r"><span>career_sim_eb</span><span> </span><span>&lt;-</span><span> </span><span>career_sim</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>career_sim_gathered</span><span> </span><span>&lt;-</span><span> </span><span>career_sim_eb</span><span> </span><span>%>%</span><span>
</span><span>rename</span><span>(</span><span>Shrunken</span><span> </span><span>=</span><span> </span><span>.fitted</span><span>,</span><span> </span><span>Raw</span><span> </span><span>=</span><span> </span><span>.raw</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>Shrunken</span><span>,</span><span> </span><span>Raw</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>career_sim_gathered</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>p</span><span>,</span><span> </span><span>estimate</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>geom_smooth</span><span>(</span><span>method</span><span> </span><span>=</span><span> </span><span>&quot;lm&quot;</span><span>,</span><span> </span><span>color</span><span> </span><span>=</span><span> </span><span>&quot;white&quot;</span><span>,</span><span> </span><span>lty</span><span> </span><span>=</span><span> </span><span>2</span><span>,</span><span> </span><span>se</span><span> </span><span>=</span><span> </span><span>FALSE</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>10</span><span>,</span><span> </span><span>100</span><span>,</span><span> </span><span>1000</span><span>,</span><span> </span><span>10000</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>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>&quot;True batting average (p)&quot;</span><span>,</span><span>
</span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;Raw or shrunken batting average&quot;</span><span>,</span><span>
</span><span>title</span><span> </span><span>=</span><span> </span><span>&quot;Empirical Bayes shrinkage reduces variance, but causes bias&quot;</span><span>,</span><span>
</span><span>subtitle</span><span> </span><span>=</span><span> </span><span>&quot;Red line is x = y; dashed white line is a linear fit&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="center" src="https://i1.wp.com/varianceexplained.org/figs/2017-01-11-simulation-bayes-baseball/career_sim_gathered_plot-1.png?w=456" data-recalc-dims="1"></p>
<p>Our method works: the the Raw (H / AB) estimates have a <em>lot</em> more noise than the shrunken estimates, just as we expected. (We filtered out cases where  in this graph: if we hadn’t, the difference would have been even starker). So</p>
<p>However, notice the dashed line representing the best-fit slope. One property that we’d prefer an estimate to have is that it’s equally likely to be an overestimate or an underestimate (that is, that ), and that’s true for the raw average. However, the shrunken estimate tends to be too high for low values of , and too low for high values of . The empirical Bayes method has introduced <strong>bias</strong> into our estimate, in exchange for drastically reducing the <strong>variance</strong>. This is a <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Bias%E2%80%93variance_tradeoff">classic tradeoff in statistics and machine learning</a>.</p>
<h3>Mean-squared error and bias relative to AB</h3>
<p>Typically, when statisticians are facing a tradeoff between bias and variance, we use mean squared error (or MSE) as a balance, which is computed as . We can easily compute that for both the raw and shrunken methods:</p>
<figure>
<pre><code data-lang="r"><span>career_sim_gathered</span><span> </span><span>%>%</span><span>
</span><span>group_by</span><span>(</span><span>type</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>summarize</span><span>(</span><span>mse</span><span> </span><span>=</span><span> </span><span>mean</span><span>((</span><span>estimate</span><span> </span><span>-</span><span> </span><span>p</span><span>)</span><span> </span><span>^</span><span> </span><span>2</span><span>))</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## # A tibble: 2 × 2
##       type     mse
##      &lt;chr>   &lt;dbl>
## 1      Raw 0.01516
## 2 Shrunken 0.00035</code></pre>
</figure>
<p>The MSE of the shrunken estimate was <em>much</em> lower than the raw estimate, as we probably could have guessed by eyeballing the graph. So by this standard, the method succeeded!</p>
<p>We’ve seen in the graph how the variance depends on , so we may want to compute the MSE within particular bins. We should use logarithmic bins (<code>10 ^ (round(log10(AB)))</code> is a handy shortcut).</p>
<figure>
<pre><code data-lang="r"><span>metric_by_bin</span><span> </span><span>&lt;-</span><span> </span><span>career_sim_gathered</span><span> </span><span>%>%</span><span>
</span><span>group_by</span><span>(</span><span>type</span><span>,</span><span> </span><span>AB</span><span> </span><span>=</span><span> </span><span>10</span><span> </span><span>^</span><span> </span><span>(</span><span>round</span><span>(</span><span>log10</span><span>(</span><span>AB</span><span>))))</span><span> </span><span>%>%</span><span>
</span><span>summarize</span><span>(</span><span>mse</span><span> </span><span>=</span><span> </span><span>mean</span><span>((</span><span>estimate</span><span> </span><span>-</span><span> </span><span>p</span><span>)</span><span> </span><span>^</span><span> </span><span>2</span><span>))</span><span>
</span><span>ggplot</span><span>(</span><span>metric_by_bin</span><span>,</span><span> </span><span>aes</span><span>(</span><span>AB</span><span>,</span><span> </span><span>mse</span><span>,</span><span> </span><span>color</span><span> </span><span>=</span><span> </span><span>type</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_line</span><span>()</span><span> </span><span>+</span><span>
</span><span>scale_x_log10</span><span>()</span><span> </span><span>+</span><span>
</span><span>scale_y_log10</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;Number of at-bats (AB)&quot;</span><span>,</span><span>
</span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;Mean-squared-error within this bin (note log scale)&quot;</span><span>,</span><span>
</span><span>title</span><span> </span><span>=</span><span> </span><span>&quot;Mean squared error is higher with raw estimate, especially for low AB&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="center" src="https://i1.wp.com/varianceexplained.org/figs/2017-01-11-simulation-bayes-baseball/metric_by_bin-1.png?w=456" data-recalc-dims="1"></p>
<p>We could also examine the bias within each bin, measured as the slope between the estimate and the true value of .</p>
<figure>
<pre><code data-lang="r"><span>career_sim_gathered</span><span> </span><span>%>%</span><span>
</span><span>mutate</span><span>(</span><span>AB</span><span> </span><span>=</span><span> </span><span>10</span><span> </span><span>^</span><span> </span><span>(</span><span>round</span><span>(</span><span>log10</span><span>(</span><span>AB</span><span>))))</span><span> </span><span>%>%</span><span>
</span><span>filter</span><span>(</span><span>AB</span><span> </span><span>></span><span> </span><span>1</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>nest</span><span>(</span><span>-</span><span>type</span><span>,</span><span> </span><span>-</span><span>AB</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>unnest</span><span>(</span><span>map</span><span>(</span><span>data</span><span>,</span><span> </span><span>~</span><span> </span><span>tidy</span><span>(</span><span>lm</span><span>(</span><span>estimate</span><span> </span><span>~</span><span> </span><span>p</span><span>,</span><span> </span><span>.</span><span>))))</span><span> </span><span>%>%</span><span>
</span><span>filter</span><span>(</span><span>term</span><span> </span><span>==</span><span> </span><span>&quot;p&quot;</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>color</span><span> </span><span>=</span><span> </span><span>type</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_line</span><span>()</span><span> </span><span>+</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>10</span><span>,</span><span> </span><span>100</span><span>,</span><span> </span><span>1000</span><span>,</span><span> </span><span>10000</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>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>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>&quot;Number of at-bats (AB)&quot;</span><span>,</span><span>
</span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;Slope of estimate/p within this bin&quot;</span><span>,</span><span>
</span><span>title</span><span> </span><span>=</span><span> </span><span>&quot;Shrunken estimates introduce bias for low AB&quot;</span><span>,</span><span>
</span><span>subtitle</span><span> </span><span>=</span><span> </span><span>&quot;Note that an unbiased estimate would have a slope of 0&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="center" src="https://i2.wp.com/varianceexplained.org/figs/2017-01-11-simulation-bayes-baseball/bias_plot-1.png?w=456" data-recalc-dims="1"></p>
<p>Another way to visualize how this tradeoff between bias and variance happens for varying AB is to recreate the above graph of the true batting average versus the estimate, this time binning by .</p>
<p><img alt="center" src="https://i0.wp.com/varianceexplained.org/figs/2017-01-11-simulation-bayes-baseball/career_sim_gathered_bin_plot-1.png?w=456" data-recalc-dims="1"></p>
<p>Notice how the variance around the true (red) line shrinks in the raw estimate, and the bias (the flatness of the gray dashed line) in the shrunken estimate decreases, until both look quite similar in the 1000+ bin.</p>
<h3>Credible intervals</h3>
<p>Besides the shrunken empirical Bayes estimates, the <code>add_ebb_estimate</code> function also adds <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Credible_interval">credible intervals</a> for each of our playesr. For example, we</p>
<figure>
<pre><code data-lang="r"><span>set.seed</span><span>(</span><span>2017</span><span>)</span><span>
</span><span>career_sim_eb</span><span> </span><span>%>%</span><span>
</span><span>sample_n</span><span>(</span><span>20</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>mutate</span><span>(</span><span>playerID</span><span> </span><span>=</span><span> </span><span>reorder</span><span>(</span><span>playerID</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>playerID</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_point</span><span>()</span><span> </span><span>+</span><span>
</span><span>geom_point</span><span>(</span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>p</span><span>),</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>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>theme</span><span>(</span><span>axis.text.y</span><span> </span><span>=</span><span> </span><span>element_blank</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% credible interval)&quot;</span><span>,</span><span>
</span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;Player&quot;</span><span>,</span><span>
</span><span>title</span><span> </span><span>=</span><span> </span><span>&quot;Credible intevals for 20 randomly selected players&quot;</span><span>,</span><span>
</span><span>subtitle</span><span> </span><span>=</span><span> </span><span>&quot;The true batting average of each player is shown in red&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="center" src="https://i0.wp.com/varianceexplained.org/figs/2017-01-11-simulation-bayes-baseball/career_sim_eb-1.png?w=456" data-recalc-dims="1"></p>
<p>Notice that out of 20 randomly selected players, the credible interval contained the true batting average (shown in red) in 19 case. This is a 95% coverage rate, which is just what we’d hoped for! Indeed, we can examine this across all players and see that 95% of the intervals contained the true probability:</p>
<figure>
<pre><code data-lang="r"><span>career_sim_eb</span><span> </span><span>%>%</span><span>
</span><span>summarize</span><span>(</span><span>coverage</span><span> </span><span>=</span><span> </span><span>mean</span><span>(</span><span>.low</span><span> </span><span>&lt;=</span><span> </span><span>p</span><span> </span><span>&amp;</span><span> </span><span>p</span><span> </span><span>&lt;=</span><span> </span><span>.high</span><span>))</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## # A tibble: 1 × 1
##   coverage
##      &lt;dbl>
## 1    0.949</code></pre>
</figure>
<p>We could also have set the threshold of the credible interval to 90%, or 75%. Does the probability that the parameter is contained within the interval change accordingly?</p>
<figure>
<pre><code data-lang="r"><span>library</span><span>(</span><span>purrr</span><span>)</span><span>
</span><span># fit the prior once
</span><span>sim_prior</span><span> </span><span>&lt;-</span><span> </span><span>ebb_fit_prior</span><span>(</span><span>career_sim</span><span>,</span><span> </span><span>H</span><span>,</span><span> </span><span>AB</span><span>)</span><span>
</span><span># find the coverage probability for each level
</span><span>estimate_by_cred_level</span><span> </span><span>&lt;-</span><span> </span><span>data_frame</span><span>(</span><span>level</span><span> </span><span>=</span><span> </span><span>seq</span><span>(</span><span>.5</span><span>,</span><span> </span><span>.98</span><span>,</span><span> </span><span>.02</span><span>))</span><span> </span><span>%>%</span><span>
</span><span>unnest</span><span>(</span><span>map</span><span>(</span><span>level</span><span>,</span><span> </span><span>~</span><span> </span><span>augment</span><span>(</span><span>sim_prior</span><span>,</span><span> </span><span>career_sim</span><span>,</span><span> </span><span>cred_level</span><span> </span><span>=</span><span> </span><span>.</span><span>)))</span><span>
</span><span>estimate_by_cred_level</span><span> </span><span>%>%</span><span>
</span><span>group_by</span><span>(</span><span>level</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>mutate</span><span>(</span><span>cover</span><span> </span><span>=</span><span> </span><span>.low</span><span> </span><span>&lt;=</span><span> </span><span>p</span><span> </span><span>&amp;</span><span> </span><span>p</span><span> </span><span>&lt;=</span><span> </span><span>.high</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>summarize</span><span>(</span><span>coverage</span><span> </span><span>=</span><span> </span><span>mean</span><span>(</span><span>cover</span><span>))</span><span> </span><span>%>%</span><span>
</span><span>ggplot</span><span>(</span><span>aes</span><span>(</span><span>level</span><span>,</span><span> </span><span>coverage</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>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>&quot;Level of credible interval&quot;</span><span>,</span><span>
</span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;Probability credible interval contains the true value&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="center" src="https://i0.wp.com/varianceexplained.org/figs/2017-01-11-simulation-bayes-baseball/unnamed-chunk-6-1.png?w=456" data-recalc-dims="1"></p>
<p>Notice that the probability (the points) hugs the red  line almost precisely. This shows that these per-observation credible intervals are generally <em>well-calibrated</em>: if you ask for a X% credible interval, you get a region that contains the true parameter about X% of the time.</p>
<h3>FDR control</h3>
<p>In <a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/bayesian_fdr_baseball/">another post</a> we considered the problem of Bayesian hypothesis testing and FDR control. In particular, we considered the problem of constructing a list of players whose true batting average was above .300, and controlling such that only (say) 10% of the players on the list were included incorrectly.</p>
<p>The q-value, which controls FDR, can be calculated with the <code>add_ebb_prop_test</code> function:</p>
<figure>
<pre><code data-lang="r"><span>pt</span><span> </span><span>&lt;-</span><span> </span><span>career_sim_eb</span><span> </span><span>%>%</span><span>
</span><span>add_ebb_prop_test</span><span>(</span><span>.3</span><span>,</span><span> </span><span>sort</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>)</span><span>
</span><span># Control for FDR of 10%
</span><span>hall_of_fame</span><span> </span><span>&lt;-</span><span> </span><span>pt</span><span> </span><span>%>%</span><span>
</span><span>filter</span><span>(</span><span>.qvalue</span><span> </span><span>&lt;=</span><span> </span><span>.1</span><span>)</span><span>
</span><span>nrow</span><span>(</span><span>hall_of_fame</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## [1] 74</code></pre>
</figure>
<p>If the FDR control were successful, we’d expect 10% of the true batting averages (<code>p</code>) to be false discoveries, and therefore below .300. Did it work?</p>
<figure>
<pre><code data-lang="r"><span>mean</span><span>(</span><span>hall_of_fame</span><span>$</span><span>p</span><span> </span><span>&lt;</span><span> </span><span>.3</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## [1] 0.0946</code></pre>
</figure>
<p>Yes- almost exactly 10% of the players included in this “hall of fame” were included incorrectly, indicating that the q-value succeeded in controlling FDR. We could instaed try this for all q-values, using the <code>cummean</code> function:</p>
<figure>
<pre><code data-lang="r"><span>pt</span><span> </span><span>%>%</span><span>
</span><span>mutate</span><span>(</span><span>true_fdr</span><span> </span><span>=</span><span> </span><span>cummean</span><span>(</span><span>p</span><span> </span><span>&lt;</span><span> </span><span>.3</span><span>))</span><span> </span><span>%>%</span><span>
</span><span>ggplot</span><span>(</span><span>aes</span><span>(</span><span>.qvalue</span><span>,</span><span> </span><span>true_fdr</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_line</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>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>&quot;q-value&quot;</span><span>,</span><span>
</span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;True FDR at this q-value threshold&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="center" src="https://i0.wp.com/varianceexplained.org/figs/2017-01-11-simulation-bayes-baseball/unnamed-chunk-8-1.png?w=456" data-recalc-dims="1"></p>
<p>Notice that the FDR was often a little bit higher than we aimed for with the q-value, which could be due to random noise. Later in this post, we’ll perform many replications of this simulation and confirm whether the FDR method was successful on average.</p>
<h3>Beta-binomial regression</h3>
<p>Most simulation analyses start with a simple model, than gradually add complications. In <a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/beta_binomial_baseball/">a post on beta-binomial regression</a>, we discovered that there is a relationship between  and the true batting average  that we need to incorporate into our model. Let’s add that complication to our simulation, and see if the method we used to account for it actually works.</p>
<p>The model described in that post had three hyperparameters: , \mu_{\mbox{AB}} and . Then each of the probabilities  was computed as:</p>
<p></p>
<p></p>
<p></p>
<p></p>
<p></p>
<p>Much as we estimated  and  from the data before using them in the simulation, we would estimate , , and  from the data:</p>
<figure>
<pre><code data-lang="r"><span>bb_reg</span><span> </span><span>&lt;-</span><span> </span><span>career</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>log10</span><span>(</span><span>AB</span><span>))</span><span>
</span><span>tidy</span><span>(</span><span>bb_reg</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## # A tibble: 3 × 6
##   parameter        term estimate std.error statistic p.value
##      &lt;fctr>       &lt;chr>    &lt;dbl>     &lt;dbl>     &lt;dbl>   &lt;dbl>
## 1        mu (Intercept)   -1.688   0.00907    -186.1       0
## 2        mu   log10(AB)    0.192   0.00280      68.7       0
## 3     sigma (Intercept)   -6.299   0.02316    -272.0       0</code></pre>
</figure>
<p>It turns out this step is pretty easy with the <code>augment</code> method of the beta-binomial prior.</p>
<figure>
<pre><code data-lang="r"><span>set.seed</span><span>(</span><span>2017</span><span>)</span><span>
</span><span>career_sim_ab</span><span> </span><span>&lt;-</span><span> </span><span>augment</span><span>(</span><span>bb_reg</span><span>,</span><span> </span><span>career</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>select</span><span>(</span><span>playerID</span><span>,</span><span> </span><span>AB</span><span>,</span><span> </span><span>true_alpha0</span><span> </span><span>=</span><span> </span><span>.alpha0</span><span>,</span><span> </span><span>true_beta0</span><span> </span><span>=</span><span> </span><span>.beta0</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>mutate</span><span>(</span><span>p</span><span> </span><span>=</span><span> </span><span>rbeta</span><span>(</span><span>n</span><span>(),</span><span> </span><span>true_alpha0</span><span>,</span><span> </span><span>true_beta0</span><span>),</span><span>
</span><span>H</span><span> </span><span>=</span><span> </span><span>rbinom</span><span>(</span><span>n</span><span>(),</span><span> </span><span>AB</span><span>,</span><span> </span><span>p</span><span>))</span></code></pre>
</figure>
<h4>Performance of beta-binomial regression method</h4>
<p>First question: are we able to extract the right hyperparameters through beta-binomial regression? We’ll fit the prior and then compare:</p>
<figure>
<pre><code data-lang="r"><span>career_ab_prior</span><span> </span><span>&lt;-</span><span> </span><span>career_sim_ab</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>log10</span><span>(</span><span>AB</span><span>))</span></code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>tidy</span><span>(</span><span>bb_reg</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## # A tibble: 3 × 6
##   parameter        term estimate std.error statistic p.value
##      &lt;fctr>       &lt;chr>    &lt;dbl>     &lt;dbl>     &lt;dbl>   &lt;dbl>
## 1        mu (Intercept)   -1.688   0.00907    -186.1       0
## 2        mu   log10(AB)    0.192   0.00280      68.7       0
## 3     sigma (Intercept)   -6.299   0.02316    -272.0       0</code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>tidy</span><span>(</span><span>career_ab_prior</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## # A tibble: 3 × 6
##   parameter        term estimate std.error statistic p.value
##      &lt;fctr>       &lt;chr>    &lt;dbl>     &lt;dbl>     &lt;dbl>   &lt;dbl>
## 1        mu (Intercept)   -1.688   0.00908    -185.9       0
## 2        mu   log10(AB)    0.193   0.00281      68.7       0
## 3     sigma (Intercept)   -6.279   0.02564    -244.9       0</code></pre>
</figure>
<p>That’s sure pretty close! It looks like beta-binomial regression was able to estimate the true parameters accurately, which suggests the resulting prior (which depends on ) will be accurate.</p>
<p>How did this prior, which depends on , affect our shrunken estimates? Again, since we’re working from a simulation we can compare the true values to the estimates, and do so within each model.</p>
<figure>
<pre><code data-lang="r"><span>career_flat_prior</span><span> </span><span>&lt;-</span><span> </span><span>career_sim_ab</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></code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>data_frame</span><span>(</span><span>method</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>&quot;Flat prior&quot;</span><span>,</span><span> </span><span>&quot;Prior depending on AB&quot;</span><span>),</span><span>
</span><span>model</span><span> </span><span>=</span><span> </span><span>list</span><span>(</span><span>career_flat_prior</span><span>,</span><span> </span><span>career_ab_prior</span><span>))</span><span> </span><span>%>%</span><span>
</span><span>unnest</span><span>(</span><span>map</span><span>(</span><span>model</span><span>,</span><span> </span><span>augment</span><span>,</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>career_sim_ab</span><span>))</span><span> </span><span>%>%</span><span>
</span><span>ggplot</span><span>(</span><span>aes</span><span>(</span><span>p</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>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>+</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>facet_wrap</span><span>(</span><span>~</span><span> </span><span>method</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;True batting average (p)&quot;</span><span>,</span><span>
</span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;Shrunken batting average estimate&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="center" src="https://i0.wp.com/varianceexplained.org/figs/2017-01-11-simulation-bayes-baseball/career_flat_prior_plot-1.png?w=456" data-recalc-dims="1"></p>
<p>Look at the bias when we don’t take the AB to batting average relationship into account: batters with low AB and low averages were universally overestimated. This is exactly the issue we had expected in the <a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/beta_binomial_baseball/">beta-binomial regression post</a>:</p>
<blockquote>
<p>Since low-AB batters are getting overestimated, and high-AB batters are staying where they are, we’re working with a biased estimate that is systematically <em>overestimating</em> batter ability.</p>
</blockquote>
<p>If you’re interested, you could take this more complex model and perform the same examinations of credible intervals and priors that we did for the simple model. (Indeed, you could incorporate some of the other trends that could affect your prior, such as year and handedness, that were considered in our <a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/hierarchical_bayes_baseball/">hierarchical model</a>).</p>
<h2>Replications</h2>
<p>So far, we ran a single simulation of our players’ batting averages, used it to perform estimation, and then examined whether our results were accurate. This is a valuable way to sanity-check the accuracy of the empirical Bayes method.</p>
<p>But what if we just got lucky? What if empirical Bayes shrinkage works about half the time, and on some datasets it gives terrible results? This is an important concern if we want to trust the method on our real data. As the next step, rather than simulating a single example, let’s create <strong>50 simulations</strong>, and run the method on each of them. We can then examine how the method performs on real data.</p>
<figure>
<pre><code data-lang="r"><span>set.seed</span><span>(</span><span>2017</span><span>)</span><span>
</span><span>sim_replications</span><span> </span><span>&lt;-</span><span> </span><span>career</span><span> </span><span>%>%</span><span>
</span><span>crossing</span><span>(</span><span>replication</span><span> </span><span>=</span><span> </span><span>1</span><span>:</span><span>50</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>mutate</span><span>(</span><span>p</span><span> </span><span>=</span><span> </span><span>rbeta</span><span>(</span><span>n</span><span>(),</span><span> </span><span>alpha0</span><span>,</span><span> </span><span>beta0</span><span>),</span><span>
</span><span>H</span><span> </span><span>=</span><span> </span><span>rbinom</span><span>(</span><span>n</span><span>(),</span><span> </span><span>AB</span><span>,</span><span> </span><span>p</span><span>))</span><span>
</span><span>sim_replications</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## # A tibble: 519,400 × 5
##     playerID     H    AB replication     p
##        &lt;chr> &lt;int> &lt;int>       &lt;int> &lt;dbl>
## 1  aaronha01  3624 12364           1 0.299
## 2  aaronha01  3090 12364           2 0.249
## 3  aaronha01  3407 12364           3 0.274
## 4  aaronha01  2418 12364           4 0.198
## 5  aaronha01  3112 12364           5 0.249
## 6  aaronha01  3253 12364           6 0.265
## 7  aaronha01  2421 12364           7 0.191
## 8  aaronha01  3075 12364           8 0.251
## 9  aaronha01  2974 12364           9 0.244
## 10 aaronha01  3186 12364          10 0.261
## # ... with 519,390 more rows</code></pre>
</figure>
<p>(This is the slowest step of our simulation: if you’re following along and want to speed it up, you could turn down the number of replications from 50).</p>
<figure>
<pre><code data-lang="r"><span>sim_replication_models</span><span> </span><span>&lt;-</span><span> </span><span>sim_replications</span><span> </span><span>%>%</span><span>
</span><span>nest</span><span>(</span><span>-</span><span>replication</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>mutate</span><span>(</span><span>prior</span><span> </span><span>=</span><span> </span><span>map</span><span>(</span><span>data</span><span>,</span><span> </span><span>~</span><span> </span><span>ebb_fit_prior</span><span>(</span><span>.</span><span>,</span><span> </span><span>H</span><span>,</span><span> </span><span>AB</span><span>)))</span></code></pre>
</figure>
<h3>Estimations of prior parameters</h3>
<figure>
<pre><code data-lang="r"><span>sim_replication_priors</span><span> </span><span>&lt;-</span><span> </span><span>sim_replication_models</span><span> </span><span>%>%</span><span>
</span><span>unnest</span><span>(</span><span>map</span><span>(</span><span>prior</span><span>,</span><span> </span><span>tidy</span><span>),</span><span> </span><span>.drop</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>)</span><span>
</span><span>sim_replication_priors</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## # A tibble: 50 × 4
##    replication alpha  beta  mean
##          &lt;int> &lt;dbl> &lt;dbl> &lt;dbl>
## 1            1  73.9   220 0.251
## 2            2  71.3   212 0.252
## 3            3  71.9   214 0.251
## 4            4  68.1   204 0.251
## 5            5  73.8   220 0.251
## 6            6  76.3   228 0.251
## 7            7  74.1   221 0.251
## 8            8  72.3   216 0.250
## 9            9  70.1   209 0.251
## 10          10  72.3   216 0.251
## # ... with 40 more rows</code></pre>
</figure>
<p>Much like our earlier individual simulation, it looks like most of these are pretty close to each other, and to the true  and . We could visualize to confirm:</p>
<p><img alt="center" src="https://i2.wp.com/varianceexplained.org/figs/2017-01-11-simulation-bayes-baseball/sim_replication_priors_plot-1.png?w=456" data-recalc-dims="1"></p>
<p>We thus notice that our estimates of , , and the mean  are mostly unbiased: generally they’re equally likely to be above or below the true parameter. We particularly note that the mean is almost always between .251 and .252; since this is what every player is being shrunk towards, it’s good that we’re so precise.</p>
<p>Our accuracy gives us confidence in the empirical Bayesian approach for this problem: we have enough data that we can feel good about estimating the hyperparameters from the data, and then using those hyperparameters as our prior.</p>
<h3>Estimates, intervals, and hypothesis testing across replications</h3>
<p>Now that we’ve seen that the prior parameters are generally estimated accurately, we can examine whether the empirical Bayes shrinkage and credible intervals worked on average.</p>
<p>One of the first metrics we examined in the simulation was the mean squared error between the true . This was . Did we just get lucky with that run?</p>
<figure>
<pre><code data-lang="r"><span>sim_replication_au</span><span> </span><span>&lt;-</span><span> </span><span>sim_replication_models</span><span> </span><span>%>%</span><span>
</span><span>unnest</span><span>(</span><span>map2</span><span>(</span><span>prior</span><span>,</span><span> </span><span>data</span><span>,</span><span> </span><span>augment</span><span>))</span></code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>sim_replication_mse</span><span> </span><span>&lt;-</span><span> </span><span>sim_replication_au</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>group_by</span><span>(</span><span>type</span><span>,</span><span> </span><span>replication</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>summarize</span><span>(</span><span>mse</span><span> </span><span>=</span><span> </span><span>mean</span><span>((</span><span>estimate</span><span> </span><span>-</span><span> </span><span>p</span><span>)</span><span> </span><span>^</span><span> </span><span>2</span><span>))</span><span>
</span><span>ggplot</span><span>(</span><span>sim_replication_mse</span><span>,</span><span> </span><span>aes</span><span>(</span><span>type</span><span>,</span><span> </span><span>mse</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_boxplot</span><span>()</span><span> </span><span>+</span><span>
</span><span>ylab</span><span>(</span><span>&quot;Mean squared error across 50 replications&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="center" src="https://i0.wp.com/varianceexplained.org/figs/2017-01-11-simulation-bayes-baseball/sim_replication_mse_plot-1.png?w=456" data-recalc-dims="1"></p>
<p>No, it looks like the MSE was always much lower than the raw estimates. This is a good sign: even in 50 examples, it never fails “catastrophically” and gives terrible estimates. This is not true of all methods!</p>
<p>We earlier saw that the credible intervals had good coverage probabilities: is this usually true?</p>
<figure>
<pre><code data-lang="r"><span>sim_replication_au</span><span> </span><span>%>%</span><span>
</span><span>mutate</span><span>(</span><span>cover</span><span> </span><span>=</span><span> </span><span>.low</span><span> </span><span>&lt;=</span><span> </span><span>p</span><span> </span><span>&amp;</span><span> </span><span>p</span><span> </span><span>&lt;=</span><span> </span><span>.high</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>group_by</span><span>(</span><span>replication</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>summarize</span><span>(</span><span>coverage</span><span> </span><span>=</span><span> </span><span>mean</span><span>(</span><span>cover</span><span>))</span><span> </span><span>%>%</span><span>
</span><span>ggplot</span><span>(</span><span>aes</span><span>(</span><span>coverage</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_histogram</span><span>(</span><span>binwidth</span><span> </span><span>=</span><span> </span><span>.001</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;% of time true value was in a 95% confidence interval&quot;</span><span>,</span><span>
</span><span>title</span><span> </span><span>=</span><span> </span><span>&quot;95% credible interval is well calibrated across replications&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="center" src="https://i1.wp.com/varianceexplained.org/figs/2017-01-11-simulation-bayes-baseball/sim_replication_au_plot-1.png?w=456" data-recalc-dims="1"></p>
<p>Yes, it looks like the coverage of a 95% credible interval was generally between 94.4% and 95.5%. Is it well calibrated at other levels: does an 80% credible interval contain the true value about 80% of the time?</p>
<p>We earlier created a plot that compares credibility level to the % of intervals containing the true parameter. We can now recreate that plot, but do it across all fifty replications:</p>
<figure>
<pre><code data-lang="r"><span>sim_replication_intervals</span><span> </span><span>&lt;-</span><span> </span><span>sim_replication_models</span><span> </span><span>%>%</span><span>
</span><span>crossing</span><span>(</span><span>cred_level</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>seq</span><span>(</span><span>.5</span><span>,</span><span> </span><span>.9</span><span>,</span><span> </span><span>.05</span><span>),</span><span> </span><span>.95</span><span>))</span><span> </span><span>%>%</span><span>
</span><span>unnest</span><span>(</span><span>pmap</span><span>(</span><span>list</span><span>(</span><span>prior</span><span>,</span><span> </span><span>data</span><span>,</span><span> </span><span>cred_level</span><span> </span><span>=</span><span> </span><span>cred_level</span><span>),</span><span> </span><span>augment</span><span>))</span><span> </span><span>%>%</span><span>
</span><span>select</span><span>(</span><span>replication</span><span>,</span><span> </span><span>cred_level</span><span>,</span><span> </span><span>p</span><span>,</span><span> </span><span>.low</span><span>,</span><span> </span><span>.high</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>sim_replication_intervals</span><span> </span><span>%>%</span><span>
</span><span>mutate</span><span>(</span><span>cover</span><span> </span><span>=</span><span> </span><span>.low</span><span> </span><span>&lt;=</span><span> </span><span>p</span><span> </span><span>&amp;</span><span> </span><span>p</span><span> </span><span>&lt;=</span><span> </span><span>.high</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>group_by</span><span>(</span><span>replication</span><span>,</span><span> </span><span>cred_level</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>summarize</span><span>(</span><span>coverage</span><span> </span><span>=</span><span> </span><span>mean</span><span>(</span><span>cover</span><span>))</span><span> </span><span>%>%</span><span>
</span><span>ggplot</span><span>(</span><span>aes</span><span>(</span><span>cred_level</span><span>,</span><span> </span><span>coverage</span><span>,</span><span> </span><span>group</span><span> </span><span>=</span><span> </span><span>replication</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_line</span><span>(</span><span>alpha</span><span> </span><span>=</span><span> </span><span>.3</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>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>&quot;Credibility level&quot;</span><span>,</span><span>
</span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;% of credible intervals in this replication that contain the true parameter&quot;</span><span>,</span><span>
</span><span>title</span><span> </span><span>=</span><span> </span><span>&quot;Credible intervals are well calibrated across 50 replications&quot;</span><span>,</span><span>
</span><span>subtitle</span><span> </span><span>=</span><span> </span><span>&quot;Red line is x = y&quot;</span><span>)</span></code></pre>
</figure>
<p><img alt="center" src="https://i2.wp.com/varianceexplained.org/figs/2017-01-11-simulation-bayes-baseball/unnamed-chunk-10-1.png?w=456" data-recalc-dims="1"></p>
<p>Each of these lines is one replication tracing from a “50% credible interval” to a “95% credible interval”: since all the lines are close to the red  line, we can see that an X% credible interval contains the true value about X% of the time. This is an important lesson of tidy simulations: whenever you can make a plot to check one simulation to check accuracy or calibration, you can also recreate the plot across many replications.</p>
<p>We can also examine our method for false discovery rate control, and see whether we can trust a q-value of (say) .05 to keep the FDR below 5%. The approach is similar to the one for credible interval coverage (so I won’t show the code here): group by each replication, then perform the same analysis we did on a single replication.</p>
<p><img alt="center" src="https://i0.wp.com/varianceexplained.org/figs/2017-01-11-simulation-bayes-baseball/sim_replication_prop_tests_plot-1.png?w=456" data-recalc-dims="1"></p>
<p>Each of these lines represents a replication tracing along every possible q-value threshold. We see that the proportion of false discoveries below a q-value is sometimes higher than the q-value promises, and sometimes lower. That’s OK: the promise of FDR control isn’t that the false discovery rate will always be exactly 5% (that would be impossible due to random noise), but that it would be on average.</p>
<p>These replications are just a start in terms of the random simulations we could perform to examine the method. For example, we could have varied the size of each replication (by randomly subsampling batters out of our total). Would our method still have been accurate if we had only 1,000 batters, or 100? We could have varied the algorithm used to estimate our hyperparameters: what if we’d used the (much faster) <a target="_blank" rel="nofollow" href="http://stats.stackexchange.com/questions/12232/calculating-the-parameters-of-a-beta-distribution-using-the-mean-and-variance">method of moments</a> to compute the beta prior, rather than maximum likelihod? I encourage you to explore other simulations you might be interested in.</p>
<h3>Conclusion: “I have only proved it correct, not simulated it”</h3>
<p>Computer scientist Donald Knuth has a <a target="_blank" rel="nofollow" href="https://staff.fnwi.uva.nl/p.vanemdeboas/knuthnote.pdf">famous quote</a>: <strong>“Beware of bugs in the above code; I have only proved it correct, not tried it.”</strong> I feel the same way about statistical methods.</p>
<p>When I look at mathematical papers about statistical methods, the text tends to look something like:</p>
<blockquote>
<p>Smith et al (2008) proved several asymptotic properties of empirical Bayes estimators of the exponential family under regularity assumptions i, ii, and iii. We extend this to prove the estimator of Jones et al (2001) is inadmissable, in the case that  is an unbiased estimator and  is convex…</p>
</blockquote>
<p>This kind of paper is an important part of the field, but it does almost nothing for me. I’m not particularly good at manipulating equations, and I get rustier every year out of grad school. If I’m considering applying a statistical method to a dataset, papers and proofs like this won’t help me judge whether it will work. (“Oh- I should have known that my data didn’t follow regularity assumption ii!”) What does help me is the approach we’ve used here, where I can see for myself just how accurate the method tends to be.</p>
<p>For example, last month I was working on a problem of logistic regression that I suspected had mislabeled outcomes (some zeroes turned to ones, and vice versa), and read up on <a target="_blank" rel="nofollow" href="https://www.jstor.org/stable/2345763?seq=1#page_scan_tab_contents">some robust logistic regression methods</a>, implemented in the <a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/robust/robust.pdf">robust package</a>. But I wasn’t sure they would be effective on my data, so I <a target="_blank" rel="nofollow" href="http://rpubs.com/dgrtwo/235656">did some random simulation</a> of mislabeled outcomes and applied the method. The method didn’t work as well as I needed it to, which saved me from applying it to my data and thinking I’d solved the problem.</p>
<p>For this reason, no matter how much math and proofs there are that show a method is , I really only feel comfortable with a method once I’ve worked with it simulated data. It’s also a great way to teach myself about the statistical method.</p>
<h3>What’s Next</h3>
<p>This post is the last of my ten-part series about empirical Bayesian methods applied to baseball batting data. This is the first time I’ve extended an analysis and discussion across a long series of posts, and it was a very positive experience. I hope you’ve enjoyed it and learned something useful.</p>
<p>But <strong>stay tuned</strong>, because this isn’t the last of my plans for the empirical Bayes series. I’ve got an upcoming announcement about how it will be expanded into a new format. (Be sure to subscribe in the form on the left if you’re interested).</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://varianceexplained.org/r/simulation-bayes-baseball/"> Variance Explained</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 title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/Oqm-Mp2FI9g"></a></div>