This week Ellis Hughes and I continue to discuss Bayesian analysis and this time we transition to talking about count data. We build an example of using the gamma-poisson conjugate to take observations of a basketball players points per game and update our prior knowledge about the points per game for our given population.
To watch the screen cast, CLICK HERE.
To access the code, CLICK HERE.
Some friends were discussing Philadelphia 76er’s point guard, Tyrese Maxey’s, three point% today. They were discussing how well he has performed over 72 games with a success rate of 43% behind the arc (at the time this data was scraped, 4/6/2022). While his percentage from 3pt range is very impressive I did notice that he has 294 attempts, which is less than 3 out of the 4 player’s that are ahead of him (Kyrie only has 214 attempts and he is ranked 3rd at the time of this writing) and Steph Curry is just behind Maxey in the ranking (42.4% success) with nearly 70 more attempts.
The question becomes, how can we be of Maxey’s three point percentage relative to those with more attempts? We will take a Bayesian approach, using a beta conjugate, to consider the success rate of these players relative to what we believe the average three point success rate is for an NBA shooter (our prior), which we will determine from observing 3 point shooting over previous 3 seasons.
NOTE: On basketball-reference.com, they have a nice check box that automatically will filter out players that are non-qualifiers for rate stats. After playing around with this, it appears that 200 attempts is their cut off. So, I will keep that and filter the data down to only those with 200 or more 3pt attempts.
All of the code, web scrapping, and csv files of the data (if you are looking to run it prior to when I scrapped it) are available on my GITHUB PAGE.
Exploratory Data Analysis
First, let’s view the top 10 three point shooters this season (size of the dot represents the number of three point attempts taken).
Visualize the distribution of three point attempts and three point% for the 2022 season, so far.
Establishing Our Prior
Since we are dealing with a binary outcome of successes (made the shot) and failures (missed the shot) we will use the beta distribution, which is the conjugate prior for the binomial distribution.
The beta distribution has two parameters, alpha and beta. To determine what these parameters should be, we will use the method of moments with the data from the previous three seasons.
To do this, we need to first find the mean and variance for the previous three seasons.
Next, we create a function that calculates alpha and beta based on the mean and variance from our observed data.
# function for calculating alpha and beta
beta_parameters <- function(dist_avg, dist_var){
alpha <- dist_avg * (dist_avg * (1 - dist_avg)/dist_var - 1)
beta <- alpha * (1 - dist_avg)/dist_avg
list(alpha = alpha,
beta = beta)
}
The function works to produce the two parameters we need. The data is returned in list format, so we will call each element of the list and store the respective values in their own variable.
The function works to produce the two parameters we need. The data is returned in list format, so we will call each element of the list and store the respective values in their own variable.
The alpha and beta parameters derived from our method of moments function appear to produce the mean and standard deviation correctly. We can plot this distribution to see what it looks like.
Update the 3pt% of the players in the 2022 season with our meta prior
We calculate our Bayes adjusted three point percentage for the players by adding their successes to `alpha` and their failures to `beta` and then calculating the new posterior percentage as
alpha / (alpha + beta)
and the posterior standard deviation as
sqrt((alpha * beta) / ((alpha + beta)^2 * (alpha + beta + 1)))
tbl2022 <- tbl2022 %>%
mutate(three_pt_missed = three_pt_att - three_pt_made,
posterior_alpha = three_pt_made + alpha,
posterior_beta = three_pt_missed + beta,
posterior_three_pt_pct = posterior_alpha / (posterior_alpha + posterior_beta),
posterior_three_pt_sd = sqrt((posterior_alpha * posterior_beta) / ((posterior_alpha + posterior_beta)^2 * (posterior_alpha + posterior_beta + 1))))
Have any of the players in the top 10 changes following in the adjustment?
We can plot the results of these top 10 players showing the posterior Bayes three point% relative to their observed three point%.
Show the uncertainty in Tyrese Maxies Performance versus Luke Kennard, who has 409 attempts
kennard <- tbl2022 %>%
filter(player == "Luke Kennard")
maxey <- tbl2022 %>%
filter(player == "Tyrese Maxey")
plot(density(rbeta(n = 1e6, shape1 = maxey$posterior_alpha, shape2 = maxey$posterior_beta)),
col = "blue",
lwd = 4,
ylim = c(0, 20),
xlab = "3pt %",
main = "Bayes Adjusted 3pt%\nBlue = Tyrese Maxey | Red = Luke Kennard")
lines(density(rbeta(n = 1e6, shape1 = kennard$posterior_alpha, shape2 = kennard$posterior_beta)),
col = "red",
lwd = 4)
If we sample from the posterior for both players, how much better is Kennard?
maxey_sim <- rbeta(n = 1e6, shape1 = maxey$posterior_alpha, shape2 = maxey$posterior_beta)
kennard_sim <- rbeta(n = 1e6, shape1 = kennard$posterior_alpha, shape2 = kennard$posterior_beta)
plot(density(kennard_sim - maxey_sim),
lwd = 4,
col = "black",
main = "Kennard Posterior Sim - Maxie Posterior Sim",
xlab = "Difference between Kennard & Maxie")
abline(v = 0,
lwd = 4,
lty = 2,
col = "red")
On average, Kennard was better in ~74% of the 1,000,000 simulations.
Long story short, Tyrese Maxey has been a solid 3pt shooter, he just happens to play on a team where James Harden takes many of the shots (maybe he should distribute the ball more?).
One last thing…..Shrinkage
So, what happened? Basically, the Bayesian adjustment created “shrinkage” whereby the players that are above average are pulled down slightly towards the population average and the players below average are pulled up slightly towards the population average. The amount of shrinkage depends on the number of attempts the player has had (the size of their sample). More attempts leads to less shrinkage (more certainty about their performance) and smaller attempts leads to more shrinkage (more certainty about their). Basically, if we haven’t seen you shoot very much then our best guess is that you are probably closer to average until we are provided more evidence to believe otherwise.
Since we were originally dealing with only players that have had 200 or more three point attempts, let’s scrape all players from the 2022 season and apply the same approach to see what shrinkage looks like.
url2022 <- read_html("https://www.basketball-reference.com/leagues/NBA_2022_totals.html")
tbl2022a <- html_nodes(url2022, 'table') %>%
html_table(fill = TRUE) %>%
pluck(1) %>%
janitor::clean_names() %>%
select("player", three_pt_att = "x3pa", three_pt_made = "x3p", three_pt_pct = "x3p_percent") %>%
filter(player != "Player") %>%
mutate(across(.cols = three_pt_att:three_pt_pct,
~as.numeric(.x))) %>%
filter(!is.na(three_pt_pct)) %>%
arrange(desc(three_pt_pct)) %>%
mutate(three_pt_missed = three_pt_att - three_pt_made,
posterior_alpha = three_pt_made + alpha,
posterior_beta = three_pt_missed + beta,
posterior_three_pt_pct = posterior_alpha / (posterior_alpha + posterior_beta),
posterior_three_pt_sd = sqrt((posterior_alpha * posterior_beta) / ((posterior_alpha + posterior_beta)^2 * (posterior_alpha + posterior_beta + 1))))
tbl2022a %>%
mutate(pop_avg = alpha / (alpha + beta)) %>%
ggplot(aes(x = three_pt_pct, y = posterior_three_pt_pct, size = three_pt_att)) +
geom_point(color = "black",
alpha = 0.8) +
geom_hline(aes(yintercept = pop_avg),
color = "green",
size = 1.2,
linetype = "dashed") +
geom_abline(intercept = 0,
slope = 1,
size = 1.2,
color = "green") +
labs(x = "Observed 3pt%",
y = "Bayesian Adjusted 3pt%",
size = "Attempts",
title = "Shirnkage of 3pt% using Beta-Conjugate",
caption = "Data Source: https://www.basketball-reference.com/leagues/NBA_2022_totals.html")
What does this tell us?
This week, Ellis Hughes and I shift our initial intro to Bayes towards solving an actual problem using the beta distribution as a conjugate prior for a binomial likelihood. We discuss what a conjugate prior is, we cover updating your prior knowledge in line, as new data becomes available, and show how the posterior distribution, produced by our Bayesian approach, works to change our beliefs about a specific outcome.
Aside from this episode, we also discussed this approach in TidyX Episode 11, where we used the a beta prior to update our knowledge about the winning percentage of professional beach volleyball teams across a season. That episode also discussed how to use the method of moments (sometimes referred to as moment matching) to select the alpha and beta parameters for the beta distribution, if you don’t have a good sense of what they should be in order to establish your prior.
To watch our screen cast, CLICK HERE.
To access our code, CLICK HERE.
To access the screen cast and code for Episode 11, CLICK HERE.
I enjoyed this short article by Andrew Gelman on the blog he keeps with several of his colleagues. The main point, which I agree 100% with, is that the examples in our stats and data analysis textbooks never seem to match what we see in the real world. The examples always seem to work! The data is clean and looks perfectly manicured for the analysis. I get it! The idea is to convey the concept of how different analytical approaches work. The rub is that once you get to the real world and look at your data you end up being like, “Uh. Wait….what is this? What do I do now?!”
The blog got me thinking about something else, though. Something that really frustrates me. If the examples we see in textbooks don’t reflect the data problems we face in the real world, what about the examples we read about in applied sport science research? How much do those examples reflect what we see in the real world?
At the risk of upsetting some colleagues, I’ll go ahead and say it:
I’m not convinced that the research we read in publications completely represents the real world either!
How can that be? This is applied science! Isn’t the real word THE research?
Well, yes and no. Yes, the data was collected in the real world, with real athletes, and in a sport setting. But often, reading a paper, looking at the aim and the conclusions, and then parsing through the methods section to see how they handled the data leaves me scratching my head.
My entire day revolves around looking at data and I can tell you; the real world is very messy.
There are potential hiccups all along the way, no matter how stringent you are in trying to keep sound data collection practices!
In research, the problem is that the data is often touched up in some manner to create an analysis sufficient for publication. Missing values need to be handled a certain way (sometimes those rows get dropped, sometimes values get imputed), class imbalance can be an issue, errant values and outliers from technology flub-ups are a very real thing, data entry issues arise, etc. These things are all problematic and, if not identified prior to analysis, can be a major issue with the findings. I do believe that most people recognize these problems and would agree with me that they are very real issues. However, it is less about knowing that there are problems but rather, figuring out what to do about them. Often the methods sections gloss over these details (I get it, word counts for journals can be a pain in the butt) and simply produce a result that, on paper at least, seems overly optimistic. As I read through the results section, without details about data processing, I frequently say to myself, “No way. There is no way this effect is as real as they report. I can’t reproduce this without knowing how they cleaned up their data to observe this effect.”
Maybe someone should post a paper about how crappy their data is in the applied setting? Maybe we should just be more transparent about the data cleaning processes we go through so that we aren’t overly bullish on our findings and more realistic about the things we can say with confidence in the applied setting?
Does anyone else feel this way? Maybe I’m wrong and I’m being pessimistic, and this isn’t as big of an issue as I believe it to be? Is the data we see in publication truly representative of the real world?
Ellis Hughes and I continue our series on simulation and sampling by moving into a basic introduction to Bayesian thinking. We code two of the common Bayesian problems that get discussed in intro courses/textbooks: (1) The probability of disease given a positive test; and, (2) the probability that a coin is biased towards heads or tails.
To watch our screen cast, CLICK HERE.
To access our code, CLICK HERE.
If you’d like to see the code in Python (I wrote a script for that too), CLICK HERE.
