1. Intro

Poisson Regression is a tool used for the modelling of counts such as injuries, goals or touchdowns.

We use data from the 2026 NBA Season to answer the question do older players make more threes, and do they increase on teams that win more?

In this guide we’ll explore concepts such as: - Data Exploration - Fitting a Poisson Regression - Account For Exposure - Check For Overdispersion - Quasi-Poisson and Negative Binomial - Comparing/Plotting Models

2. Why Choose Poisson?

Count data often breaks the rules of linear regression. It is rightly skewed therefore does not follow the normal distribution needs of linear regression. Poisson uses a log scale which prevents any negative numbers which is ideal for count data. Linear regression assumes the variance is similar but with counts the spread grows with the average.

Poisson regression is built for this. It assumes the count \(Y\)$Y$ follows a Poisson distribution with average rate \(\lambda\), and links that rate to your predictors through a log function:

\[\log(\lambda) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots\]

3. Setup

Install these packages if you haven’t already, then load them.

install.packages("tidyverse")
install.packages("hoopR")
install.packages("performance")
installed.packages("ggeffects")
installed.packages("sjPlot")
# The MASS package is already installed with R
library(MASS)       
library(hoopR)      
library(tidyverse)
library(performance)
library(ggeffects)
library(sjPlot)

4. NBA Data

Using the hoopR package the nba dataframe was created, which included games, minutes, 3 pointers made, 3 pointers attempted, age and wins.

nba_data <- nba_leaguedashplayerstats(
  season   = "2025-26",
  per_mode = "Totals"       # the key argument: totals, not per game
)

nba <- nba_data$LeagueDashPlayerStats |>
  mutate(across(c(GP, MIN, FG3M, FG3A, AGE, W), as.numeric)) |>
  select(GP, MIN, FG3M, FG3A, AGE, W)

5. Data Exploration

Make sure you always look at the data you’re trying to model. By visualising the data using techniques such as plotting, it helps the reader and yourself to understand what is happening in the dataset.

ggplot(nba, aes(FG3M)) +
  geom_histogram(binwidth = 2, fill = "steelblue", colour = "white") +
  labs(x = "3pt Made", y = "Number of players",
       title = "Distribution of 3 Pointers") +
  theme_minimal()

The distribution of 3 pointers made is right-skewed with most players making less than 20. This is typical in count data.

This next plot shows the effect of exposure. Players that play more will naturally make more 3 pointers.

ggplot(nba, aes(MIN, FG3M)) +
  geom_jitter(width = 0.2, height = 0.2, alpha = 0.4) +
  geom_smooth(method = "loess", se = FALSE, colour = "firebrick") +
  labs(x = "Minutes in the NBA", y = "3 Pointers Made",
       title = "3 Pointers vs Minutes") +
  theme_minimal()

If exposure was ignored then an average shooter who plays significant minutes could be seen as an elite shooter.

Finally, here is a look at age vs 3 pointers.

ggplot(nba, aes(AGE, FG3M)) +
  geom_jitter(width = 0.2, height = 0.2, alpha = 0.4) +
  geom_smooth(method = "loess", se = FALSE, colour = "firebrick") +
  labs(x = "Age of Player", y = "3 Pointers Made",
       title = "Age vs 3 Pointers Made") +
  theme_minimal()

## 6. Fitting Poisson Model With An Offset

We use glm() (generalised linear model). To make this model poisson regression this needs to be added family = poisson.

To model the 3 pointers rate per minute rather than the raw count, we add offset = log(Min) to the formula. This tells the model that the expected count is proportional to the minutes played.

p_model <- glm(FG3M ~ AGE + W,
                  offset = log(MIN),
                  family = poisson,
                  data = nba)

summary(p_model)
## 
## Call:
## glm(formula = FG3M ~ AGE + W, family = poisson, data = nba, offset = log(MIN))
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.368710   0.036000 -93.576  < 2e-16 ***
## AGE          0.014625   0.001273  11.484  < 2e-16 ***
## W            0.002482   0.000392   6.332 2.43e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 11682  on 581  degrees of freedom
## Residual deviance: 11500  on 579  degrees of freedom
## AIC: 14255
## 
## Number of Fisher Scoring iterations: 5

Reading the Estimate output, it is on the log scale. The p values (Pr(>|z|)) are significantly low however we must check the Poisson assumptions and dispersion before trusting these values.

7. Check for overdispersion

Poisson assumes that mean is equal to variance. If the dataset has more variability (dispersion ratio > 1.5) then the errors are too small and the p-values cannot be trusted. This is called overdispersion.

performance::check_overdispersion(p_model)
## # Overdispersion test
## 
##        dispersion ratio =   15.939
##   Pearson's Chi-Squared = 9228.474
##                 p-value =  < 0.001

These results show that simple Poisson regression does not fit this dataset and we should used a model that will account for the extra variability.

8. Quasi-Poisson and Negative Binomial

Quasi-Poisson Same coefficients as the Poisson model, but the standard errors are scaled up by the dispersion value.

q_model <- glm(FG3M ~ AGE + W,
                  offset = log(MIN),
                  family = quasipoisson,
                  data = nba)

summary(q_model)
## 
## Call:
## glm(formula = FG3M ~ AGE + W, family = quasipoisson, data = nba, 
##     offset = log(MIN))
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.368710   0.143723 -23.439  < 2e-16 ***
## AGE          0.014625   0.005084   2.877  0.00417 ** 
## W            0.002482   0.001565   1.586  0.11330    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 15.93864)
## 
##     Null deviance: 11682  on 581  degrees of freedom
## Residual deviance: 11500  on 579  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5

Negative Binomial. This is usually the preferred model to handle overdispersion. Negative Binomial adds an extra parameter that accounts for variance.

nb_model <- glm.nb(FG3M ~ AGE + W + offset(log(MIN)), data = nba)

summary(nb_model)
## 
## Call:
## glm.nb(formula = FG3M ~ AGE + W + offset(log(MIN)), data = nba, 
##     init.theta = 1.56348415, link = log)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.348086   0.224487 -14.914   <2e-16 ***
## AGE          0.009933   0.008250   1.204   0.2286    
## W            0.004568   0.002231   2.048   0.0406 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.5635) family taken to be 1)
## 
##     Null deviance: 695.54  on 581  degrees of freedom
## Residual deviance: 689.43  on 579  degrees of freedom
## AIC: 5123
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.563 
##           Std. Err.:  0.106 
## 
##  2 x log-likelihood:  -5115.003

Dispersion parameter for Negative Binomial(1.5635) which is significantly better performing than simple Poisson regression.

9. Model Comparison

We’ll look at the AIC comparison to determine which is the better model.

tab_model(p_model, nb_model)
  FG 3 M FG 3 M
Predictors Incidence Rate Ratios CI p Incidence Rate Ratios CI p
(Intercept) 0.03 0.03 – 0.04 <0.001 0.04 0.02 – 0.05 <0.001
AGE 1.01 1.01 – 1.02 <0.001 1.01 0.99 – 1.03 0.229
W 1.00 1.00 – 1.00 <0.001 1.00 1.00 – 1.01 0.041
Observations 582 582
R2 Nagelkerke 0.269 0.015

Negative Binomial is by far the better model for this question and dataset. AIC is over 9000 less than the Poisson regression (lower is better). Additionally the results show that older players make slightly more 3s per minute however since the p values is large we can’t say age matter. Since win the p value around 0.05 wins relate to more 3s per minute however it is weak evidence of a small effect.

10. Predictions

Here we use the model to predict 3 point makes at every age with players being categorised in amount of wins.

# Predict 
predict_3 <- ggpredict(nb_model, 
                     terms=c("AGE", "W"))

# Plot - Predictions
plot(predict_3)

There is a lot of variation in the model and uncertainty in the predictions. This shows that age and wins are not ideal for 3 point make predictions and other variables should be used.