Hierarchical Bayesian Batting Ability, with Multiple Comparisons

by

Just in time for the last game(s) of this year’s World Series, the final installment of my example of Bayesian stats using baseball batting ability. I took us off-topic (text analytics) with a general intro to Bayesian stats (What is “Bayesian” Statistical Inference?). As a first example, I compared Bayesian calculations of binomial posteriors with point estimates (Batting Averages: Bayesian vs. MLE Estimates). In exploring the question of “where do priors come from?”, I started with simple non-Bayesian point estimates (Moment Matching for Empirical Bayes Beta Priors for Batting Averages). Finally, I realized I’d meant “ability” where I’d said “average” (former is a latent parameter, latter is a statistic calculated from at bats), when I considered Bayesian point estimates (Bayesian Estimators for the Beta-Binomial Model of Batting Ability).

Hierarchical Model Recap

For batter j, the number of hits n_j in N_j at bats is modeled as a binomial,

n_j \sim \mbox{\sf Binom}(\theta_j,N_j) for j \in 1:J,

where the ability, or chance of getting a hit, for batter j is \theta_j. Ability is modeled as having a beta distribution

\theta_j \sim \mbox{\sf Beta}(\alpha,\beta) for j \in 1:J

with prior number of hits \alpha-1 and prior number of outs \beta-1. These parameters, which act as priors for the for the binomial parameter, are themselves given priors. The mean of the beta, \alpha/(\alpha+\beta), is given a uniform prior,

\alpha/(\alpha+\beta) \sim \mbox{\sf Beta}(1,1).

The prior for the scale \alpha+\beta is modeled by a Pareto distribution, which has probability proportional to 1/(\alpha+\beta)^{2.5},

\alpha + \beta \sim \mbox{Pareto}(1.5).

Fitting

I used the 2006 AL position player data (given in this previous blog post). That fixes the number of players J and the hits n_j and at-bats N_j for each player j \in 1:J.

I then used BUGS encoding of the model above (as shown in this previous post). BUGS automaticaly implements Gibbs sampling, a form of Markov Chain Monte Carlo (MCMC) inference. I ran 3 chains of 1000 samples each, retaining the last 500 samples, for 1500 posterior samples total. All parameters had \hat{R} values very near 1, indicating convergence of the Markov chains. As usual, we read off statistics from the samples and used the sampled values for inference, where they allow the integrals involved in Bayesian inference (as descirbed in this blog post) to be computed.

Beta Posterior

We can inspect the posterior for the beta mean \alpha/(\alpha+\beta) and scale \alpha+\beta parameters. We plot them as a 2D scatterplot with their means picked out as lines.

beta parameters posterior

As usual, the scale is much harder to estimate than the mean.

Ability Posteriors

We can also plot the ability for particular players. Here’s histograms of sampled posteriors for the players with the best average posterior ability estimates, with their hits and at-bats provided as labels above the plot:

Notice how the variance of the posterior is determined by the number of at bats. The player with 60 at bats has a much wider posterior than the player with 695 at bats.

Multiple Comparisons: Who’s the Best Batter?

We’re going to do what Andrew, Jennifer and Masanao suggest, which is to use our hierarchical model to make a posterior comparison about player’s ability. In particular, we’ll estimate the posterior probabability that a particular player is the best player. We simply look at all 1500 posterior samples, which include ability samples as plotted above, and count how many times a player has the highest ability in a sample. Then divide by 1500, and we’re done. It’s a snap in R, and here’s the results, for the same batters as the plot above:

Average At-Bats Pr(best ability)
.367 60 0.02
.342 482 0.08
.347 521 0.12
.322 695 0.02
.330 648 0.04
.343 623 0.11
.330 607 0.04

The .347 batter with 521 at bats not only has the highest estimated chance in our model of having the highest ability, he also won the batting crown (Joe Mauer of Minnesota). The beta-binomial hierarchical model estimate is only 12% that this batter has the highest ability. The estimate is very close for the .343 batter with 623 at bats (Derek Jeter of the Yankees). [It turns out the race to the batting crown came down to the wire.]

The larger number of at bats provides more evidence that the batter has a higher ability than average, thus pulling the posterior ability estimate further away from the prior mean. Finally, note that we’re assigning some chance that the .367 batter with only 60 at bats is best in the league. That’s because when the samples are on the high side of the distribution, this batter’s best.

4 Responses to “Hierarchical Bayesian Batting Ability, with Multiple Comparisons”

  1. Gregor Gorjanc Says:

    You can also obtain ranks directly with BUGS – take a look at functions rank*.

  2. Probability Measures and Random Variables « LingPipe Blog Says:

    [...] multivariate density functions explicitly and reason through integration. For instance, in the hierarchical batting ability model I discussed as an example of Bayesian inference, the model is fully described by the [...]

  3. Falk Says:

    Do you think it is possible to derive analytical update equations for alpha and beta?

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s


Follow

Get every new post delivered to your Inbox.

Join 797 other followers