Rating college basketball teams (with probabilistic programming!)

So lately I have been getting very interested in probabilistic programming, specifically probabilistic graphical models, and wondering how I can apply them to college basketball ratings. I really like the idea of being to able to specify data problems graphically through probability models. So I started looking around the internet for similar examples of probabilistic graphical models (PGMs) applied to sports ratings, specifically implementations using Python interfaces like PyMC.

With some success, I found a really nice write-up on building a hierarchical Bayesian model to rate 20 soccer teams of the Premier league. The post is actually (mostly) a reconstruction of work performed in a similar paper found here. Anyway, they both use relatively simple, intuitive models to describe how points (goals in this case) are produced in a soccer game. I'll adapt it here to build a similar process for college basketball games.

A scoring model for college basketball¶

There are many, many different approaches we could take to rating sports team's and, indeed, many have been applied in practice. Being able to better predict the outcomes of games that people care about and that people bet money on, it turns out, is quite appealing. One challenge that college basketball has, moreso than other sports, is that it is more difficult to measure a team's strength based off of traditional statistics that can be directly observed. We cannot simply look at wins/losses or points scored/allowed directly because college basketball team's play vastly different strengths of schedule. A team's observed output depends on the true skill of the team's that they have played. And the true skill of each team that they have played depends, in turn, on the true skill of the teams those teams have played. So for college basketball ratings we can't just take some stats ("advanced" or otherwise) from the box score, shove them into a feature matrix, and run some scikit-learn algorithm (well, we could, but it wouldn't be very good). We need to think a bit deeper.

To solve our issues here, let's not think about how what we observe determines team ratings, but about how team ratings determine what we observe. This is a latent variable problem. When two team's play their ratings, which are hidden to us, combine to determine (approximately) what the outcome will be.

To formulate the process, we say that each team possesses an "offensive skill" rating and a "defensive skill" rating. When team A playes team B, the expected points that team A will score is determined by team A's offensive skill and team B's defensive skill. But, as most college basketball fans (or sports fans in general) can tell you, if two team's play multiple games under exactly the same conditions, the actual score for team A will vary game-to-game. We in the data game like to call that randomness. So, a reasonable model for this randomness might be to assume that the true score is drawn from a normal distribution centered at the expected points scored by team A. This model will be slightly more interpretable if we add an intercept, such that the offensive and defensive "skills" can be interpreted as points above or below average.

For example, say the Bayes Ballers matchup versus Newton Netters where the two team's skills are:

We can interpret this as the Bayes Ballers scoring 5 points more than the league average and allowing 2 points more to opposing offenses than the league average. For the Newton Netters, they score 1 point less than the league average but allow 3 points less on defense. Assuming the league point average is 70 points, the expected points for each team are:

• Bayes Ballers expected points: 5 + -3 + 70 = 72
• Newton Netters expected points: -1 + 2 + 70 = 71

So, on average, we expect the Bayes Ballers to defeat the Newton Netters by a score of 72-71, a margin of one point. Still, that does not mean that the Newton Netters won't win from time to time. After all, the two team's seem pretty evenly matched and we know there is some randomness in the scores. Above, we said that the expected points are fully determined by the team's ratings, but the true points, i.e. the points that we actually observe, are noisy or have some randomness. We chose to model that by saying that the true points are drawn from a constant mean Gaussian distribution centered at each team's expected point totals.

We now have a model that allows us to explain the scores we see in college basketball games, and we know they are partly due to skill, and partly due to randomness. We can show this as a probabilistic graphical model via a DAG structure.

Note: We'll also add a home court advantage parameter to our model, which will come in handy using real data.

TODO

Cool, we have now described a Bayesian network structure and all we should have left to do is to dump it into PyMC(3) and hit the inference button!

Not so fast.

An old friend¶

Let's restate what we have here: we have a model where a linear combination of explanatory variables (offensive/defensive skills, intercept, home advantage) influences the expected value of the response variable, which is a random variable drawn from a normal distribution. Sound familiar? That's right, this is just the description of generalized linear model with a Gaussian family and identity link function! We got all excited to use our fancy probabilistic programming tools, equipped with our pretty graphical model, and all the while the dusty old linear model was sitting in the corner just waiting for us to notice.

In fact, it turns out that this type of ratings model was made prominent in the college basketball landscape by Kenneth Massey, who details such a model in this mathematical paper from 1997. Despite the wonky formatting, it's definitely worth a read for those looking for a gentle, mathematical introduction to sports ratings.

How shall we proceed then? We got all excited to use our PyMC inference button, only to realize that we had created a glorified linear regression that's been sitting in some pre-LaTeX era word document since 1997?! Despite our initial disappointment, it will be informative and even interesting to carry out this exercise using both the boring, old GLM and our PyMC inference button. To start, we'll work through an example with some made-up data using the Statsmodels Python package for the GLM, and PyMC3 for the Bayesian network. We'll use some helper functions to synthesize data for a number of games between 5 teams.

Validating models using simulated data¶

For each team we made up an offensive and defensive skill rating, and then randomly created pairings of the teams as games. In each game, we compute an expected home score E[Home] and an expected away score E[Away] and then use these values as the means of normal distributions from which the true scores (Home Score and Away Score) are drawn. We chose the standard deviation for this normal distribution to be a constant (i.e. homoscedastic) value of 3. To verify this, we can plot the errors between actual and expected scores, which should be distributed according to $\mathcal{N}(0, 3)$

The data is pretty close to the normal distribution we'd expect, give or take some randomness. Now, to actually use this in a ratings model, we need to do some data wrangling. There are actually two sample points per row in the current dataframe, one for the home team and one for the away team. Let's "stack" the data which will double the number of rows, and drop some columns we don't need for the model.

In our new feature dataframe, the first two rows correspond to game 1, which paired team 0 (home) vs team 1 (away), while the second game was team 4 (home) vs team 3 (away), and so on. Still, for use in a GLM model, we cannot simply use the "Off. Team" and "Def. Team" columns as features because they are categorical. Though these columns have numeric values, they are just team identifiers and carry no true numeric meaning. We can use pandas to generate dummy variables, or binary encodings, for our team columns

Taking the subset of these columns that we care about, we end up with the following linear system:

$$ \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 1\\ 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & -1 & 1\\ &&&&&&\vdots\\ \end{bmatrix} \begin{bmatrix} \beta^o_0\\ \beta^o_1\\ \vdots\\ \beta^d_0\\ \beta^d_1\\ \vdots\\ home\\ intercept\\ \end{bmatrix} = \begin{bmatrix} 76.223\\ 61.633\\ \vdots\\ \end{bmatrix} $$

A system of this form can be solved easily via the classical method of least squares using whatever one of the many solver packages you like. However, the astute reader will have detected a small (but absolutely critical!) problem with the system we have constructed. When we encoded the dummy variables, we did not drop any of the columns and so we are left with linearly dependent columns in our feature matrix. To be more explicit, one can determine, exactly, the fifth offensive rating column given any of the other four (likewise with the defense).

We can also view this slightly differently--since we have an intercept, we could add an arbitrary offset to all the offensive ratings and subtract the same from the intercept and still have a valid solution. So, we need to impose some identifiability constraints. We can do this simply by mandating that the offensive ratings all sum to zero (i.e. the average team has a zero offensive rating) and do the same for the defense. This can be interpreted as a constraint on the last columns of the offense and defense:

$$ \sum_{i=1}^{N} \beta_i^o = 0, \sum_{i=1}^{N} \beta_i^d = 0\\ -\sum_{i=1}^{N-1} \beta^o_i = \beta^o_N, \: -\sum_{i=1}^{N-1} \beta^d_i = \beta^d_N $$

Where a one appeared in the last column of our feature matrix before, we can replace it with -1 in all the other columns and we will have effectively enforced the sum to zero constraints. Our new dataframe is shown below:

We can now solve the linear system to find the maximum likelihood estimates of the team ratings. We'll use the statsmodels GLM API to solve the system and get some summary statistics about the model fit.

Hmmmmm, well they certainly aren't dead on are they? But we shouldn't really expect them to be since there is some randomness and we only sampled 200 games (see also: central limit theorem). Since we're inferring the coefficients from random sample, we would expect them to vary from sample to sample. In fact, it is well established in theory that the regression coefficients have a sampling distribution, which can be computed analytically, and what we see here is just a single sample from that distribution (please do read more about this here). Statsmodels is so nice as to provide us with 95% confidence interval bounds as a part of the returned model summary. Instead of comparing directly to the sampled values, let's compare the true values to the confidence intervals.

The true coefficient is contained within the confidence interval bounds that we computed from our linear model fit for each offensive coefficient. We should not, in general, expect every confidence interval to contain the true value of its respective parameter. These are 95% confidence intervals, so a non-negligible 5% of the time the value will be outside of the interval.

We now have a data model that we have validated with synthesized data and provides us with a sampling distribution for each of our parameters--I am quite certain we will find interesting ways to use this on real data. Still, I gotta admit, that wasn't much of a thrill. A linear model? C'mon, I'm a data adrenaline junkie and it's Bayesian networks that get my pulse pounding!

A Bayesian network ratings model in PyMC3¶

Well, luckily for us Pythonistas, PyMC3 happens to be a wonderful library for implementing directed graphical models in Python. Even though we might not need to implement this same model in PyMC3, it will be useful in future endeavors where we may expand our data model. It will also serve as a great way to get started with PyMC3 and probabilistic programming. Without further ado, we'll declare our model using PyMC3 below.

What just happened here? First, we used PyMC's find_MAP() method to find the maximum a posteriori point estmate of the regression coefficients. That is, given the data we observed and our priors on the coefficients, this MAP estimate is the most likely point estimate for the regression parameters. We then used that MAP estimate as the starting point in our parameter space to begin taking samples from the posterior. The "how" of MCMC sampling is a deep subject and I couldn't even begin to do it justice in this small paragraph. For an excellent introduction to the subject of MCMC sampling, please take a look at Chapter 3 of "Bayesian Methods for Hackers" by Cam Davidson-Pilon; to take what you learn there a bit further, Thomas Wiecki has a really nice write-up on his blog, summarizing the inner workings of the Metropolis sampling method.

All we need to know here is that PyMC is generating 3000 points from the posterior - 3000 points in our parameter space where the true parameters are most likely to be. This will provide us with a histogram for each parameter, which we can interpret as an approximation to the true posterior probability distribution of each parameter.

PyMC provides the raw samples, usually called the "trace", from the model fit, but it also provides us with a nice summary table, seen here:

We see that for each parameter, we get the mean and standard deviation of the posterior distribution, as well as the 95% highest posterior density region, which for unimodal distributions corresponds to the 95% credible region. Let's check that the true values are within the credible regions for the offensive parameters.

The results here and the interval bounds matchup nicely with those from the GLM model which is a good sign that our models agree. However, we still want to check the posterior distributions from PyMC and compare them with the sampling distributions from the GLM.

Validating the posterior distribution¶

As Bayesians, we seek out the posterior distribution of our model's parameters given the data observed. Combining our prior beliefs with actual data provides with a posterior belief about where the true values of the parameters might lie. In fact, we get an entire k-dimensional probability distribution for the parameter values which tells us not only the most likely posterior values of the parameter (the MAP estimate) but also provides us with measures of uncertainty. This is better than simply having a point estimate for a variety of reasons. In a college basketball context, the Bayesian posterior distribution can tell us the associated uncertainty in a prediction. A gambling man could use this certainty to select less risky bets to place on games.

But you might say, the GLM model above also provided us with a distribution on the parameters. Wrong! Err... sort of. We should be quite careful not to confuse the sampling distribution of the coefficients in a linear regression with the posterior distribution of the coefficients in a Bayesian linear regression. This is a matter of frequentist vs Bayesian viewpoints and the topic can become quite obscure rather quickly, but the difference mostly comes down to a fundamental difference in the interpretation of what probability actually means. I really must implore you to, at this point, run on over to Quora and read this nice explanation by Michael Hochster, "What is a confidence interval in layman's terms?". There's rather a lot to say on this subject, but I will say none of it here and leave it as a most gruesome exercise for the interested reader. (I also rather enjoyed reading this paper here on the matter)

Despite all this, the Bayesian posterior and the sampling distribution of the GLM coefficients do happen to match in this case. Let's now compare the sampling distributions from the GLMs with the posterior distributions from the Bayesian model. First, let's take a look at the offensive and defensive ratings.

Ah, warm fuzzies. The Bayesian model matches the GLM, and we cry tears of unadultered joy at the purity of our universe. Anyway, let's plug on and check the intercept and home court advantage as well.

We nailed it! At this point, we're feeling pretty good. We came up with a nice model for rating sports teams, implemented it in two separate ways that ought to provide equivalent results, and we validated those results against each other and the true, known values that we made up. I highly suggest starting any analysis, especially with probabilistic models, by generating synthetic data and validating that the results are as expected against known parameters. One point of caution, however, is that you should not expect, in general, to recover 100% of the true, known parameters inside of Bayesian credible intervals. Remember, these are 95% credible intervals such that if you have 100 parameters, you should expect that ~5 of those credible regions will not contain the true value.

Well, we've done our due diligence and we've thoroughly checked our models for completeness. I think we're finally ready to move on to using some real college basketball data in our analysis. First we'll rate. Then, we'll profit. I think.

Rating real college basketball teams in 2016¶

Ideally, I'd like to provide a way to get the data I've collected for anyone to use, but I don't have an easy solution in mind at the moment. For now, I'll load data from a local PostgreSQL database which contains box score information for recent seasons of college basketball. I'll also use a helper function from a utility class that will return the data in the format that we need. Additionally, I'll load in Ken Pomeroy's publicly available ratings for the 2016 season so that we have something to compare to. In general, we should not assume that these ratings will be better than the model we've come up with, but given that Kenpom's ratings are the result of years of work, we should probably suspect they are better than this simple model we'll use here.

Also note that instead of using actual point scores as our rating metric, we'll use points per possession (ppp) instead. You can read more on this topic here, but basically ppp is a metric that normalizes a team's score by the number of opportunities it actually had to score, which is a more apt representation of their skill.

Let's inspect the games dataframe.

Some important columns:

• i_hteam, i_ateam - a unique index for home and away teams, ranging from 0 to num_teams - 1
• hppp, appp - home and away points per possessions for each game. Computed according to the description here.
• neutral - boolean indicating whether or not game was played on a neutral court (this will be useful in assigning home court advantage)
• hteam_id, ateam_id - home and away team ids (these correspond to ids used on the website from which the data is collected)

Below, we perform similar data wrangling as before so that we can get the data in the form we need for both the linear model and the Bayesian model.

With our data prepared in nice formats for both models (GLM and Bayesian), we can now train each model on the data.

With both models trained, let's check the GLM sampling distributions versus the Bayesian posterior densities as we did with synthesized data.

Once again, the Bayesian posteriors are in agreement with the sampling distributions from the GLM, this time with a vastly increased number of parameters. Still, just because these two models agree with one another does not mean that they are correct. Before, we had true values to compare to because we generated the data ourselves, but this time we don't know the actual ratings (that's the entire point). Instead, we can compare to Ken Pomeroy's ratings - while these are not "true" values, they are a pretty good time-tested set of ratings to compare with.

Kenpom rankings validation¶

As a final inspection in this exercise, we will check that our ranking of each of the team's aligns with the Kenpom rankings. However, since we should not expect these to match up perfectly, it will help to compare against a baseline model. For now, we can just check that our rankings are better than simply ranking teams on their raw average output in each game. First, let's compute each team's raw average.

It turns out that Kenpom's ratings do not contain one of our team's, Texas Pan American, and since we are comparing rankings we'll just drop it from our dataframe.

First, take a look at the summary statistics of the ranking distances dataframe. The 25% and 75% quantiles are in the 20-30 ranking spots range for the average ratings to kenpom ratings, while our GLM and Bayesian ratings have improved these to the 10-15 range. Additionally, the standard deviations of the distances are greatly reduced. Finally, we can inspect the scatter plots above to see that our ratings are much more tightly coupled to the Kenpom ratings than are the raw averages. We could go further still, but this should be enough to convince ourselves that our ratings are significantly better than our average baseiline ratings.

On a related note, there are many reasons that our ratings will not exactly match Ken Pomeroy's. To be more specific, there are many reasons why his ratings are better than this first version. For one, he applies a temporal weighting on games that essentially phases out really old games that aren't likely to reflect a team's current skill as well as recent games. He also reduces the weight of enormous blowouts such that a team that beats an inferior opponent by 60 does not get (much) more credit than a team that beats an inferior opponent by 20. There are other good reasons as well (hey, the guy's been doing this a while), but with some time and experimentation, our rating system here may get better. This brings me to my next point: what does "better" mean?

Most likely, we are rating sports teams because we'd like to be able to intelligently predict the outcome of a matchup before it happens. We might care about predicting the actual score, or just predicting the winner, but either way we will evaluate our ratings on how well they predict future data. In this exercise, we rated each of the team's in Division 1 college basketball after all the games had been played... this doesn't sound very useful! In a follow up post, I'll take a look at generating rankings throughout the season and applying them for prediction of past "future" games.

Conclusions¶

Saying that you used probabilistic programming to create a rating system for college basketball teams using a Bayesian network definitely sounds cool. But as we saw, a linear regression gave us all we needed to create a solid, effective rating system. So why did we bother with the Bayesian network? For one, it was a nice way to learn PyMC and probabilistic programming. Second, the Bayesian network will give us nice flexibility to extend our model in the future, and to also address some problems we haven't discussed. How should we rate teams in the beginning of the season when there is really very little data to work with? The Bayesian model will allow us to specify priors on the team ratings, which we could base off of a variety of things, the simplest being the previous year's ratings. Bayes nets also allow us to reason about the uncertainty in our ratings more clearly. We should be clear not to confuse the sampling distributions of the coefficients provided by the GLM with the posterior distribution on the parameters we get from the Bayesian model. Among other things, the Bayes net gives us the ability to say things like "team A's rating is between 0.9 and 1.1 with 95% probability," something we cannot say using frequentist methods.