The Poisson distribution: From basic probability theory to regression models

[This article was first published on Achim Zeileis, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Brief introduction to the Poisson distribution for modeling count data using the distributions3 package. The distribution is illustrated using the number of goals scored at the 2018 FIFA World Cup, suitable for self-study or as a classroom exercise.

The Poisson distribution

The classic basic probability distribution employed for modeling count data is the Poisson distribution. Its probability mass function \(f(y; \lambda)\) yields the probability for a random variable \(Y\) to take a count \(y \in \{0, 1, 2, \dots\}\) based on the distribution parameter \(\lambda > 0\):

[\text{Pr}(Y = y) = f(y; \lambda) = \frac{\exp\left(-\lambda\right) \cdot \lambda^y}{y!}.]

The Poisson distribution has many distinctive features, e.g., both its expectation and variance are equal and given by the parameter \(\lambda\). Thus, \(\text{E}(Y) = \lambda\) and \(\text{Var}(Y) = \lambda\). Moreover, the Poisson distribution is related to other basic probability distributions. Namely, it can be obtained as the limit of the binomial distribution when the number of attempts is high and the success probability low. Or the Poisson distribution can be approximated by a normal distribution when \(\lambda\) is large. See Wikipedia (2002) for further properties and references.

Here, we leverage the distributions3 package (Hayes et al. 2022) to work with the Poisson distribution in R. In distributions3, Poisson distribution objects can be generated with the Poisson() function. Subsequently, methods for generic functions can be used print the objects; extract mean and variance; evaluate density, cumulative distribution, or quantile function; or simulate random samples.

library("distributions3")
Y <- Poisson(lambda = 1.5)
print(Y)
## [1] "Poisson distribution (lambda = 1.5)"
mean(Y)
## [1] 1.5
variance(Y)
## [1] 1.5
pdf(Y, 0:5)
## [1] 0.22313 0.33470 0.25102 0.12551 0.04707 0.01412
cdf(Y, 0:5)
## [1] 0.2231 0.5578 0.8088 0.9344 0.9814 0.9955
quantile(Y, c(0.1, 0.5, 0.9))
## [1] 0 1 3
set.seed(0)
random(Y, 5)
## [1] 3 1 1 2 3

Using the plot() method the distribution can also be visualized which we use here to show how the probabilities for the counts \(0, 1, \dots, 15\) change when the parameter is \(\lambda = 0.5, 2, 5, 10\).

plot(Poisson(0.5), main = expression(lambda == 0.5), xlim = c(0, 15))
plot(Poisson(2),   main = expression(lambda == 2),   xlim = c(0, 15))
plot(Poisson(5),   main = expression(lambda == 5),   xlim = c(0, 15))
plot(Poisson(10),  main = expression(lambda == 10),  xlim = c(0, 15))

Probability density for Poisson distributions with means 0.5, 2, 5, and 10

In the following we will illustrate how this infrastructure can be leveraged to obtain predicted probabilities for the number of goals in soccer matches from the 2018 FIFA World Cup.

Goals in the 2018 FIFA World Cup

To investigate the number of goals scored per match in the 2018 FIFA World Cup, the FIFA2018 data set provides two rows, one for each team, for each of the 64 matches during the tournament. In the following, we treat the goals scored by the two teams in the same match as independent which is a realistic assumption for this particular data set. We just remark briefly that there are also bivariate generalizations of the Poisson distribution that would allow for correlated observations but which are not considered here.

In addition to the goals, the data set provides some basic meta-information for the matches (an ID, team name abbreviations, type of match, group vs. knockout stage) as well as some further covariates that we will revisit later in this document. The data looks like this:

data("FIFA2018", package = "distributions3")
head(FIFA2018)
##   goals team match type stage logability difference
## 1     5  RUS     1    A group     0.1531     0.8638
## 2     0  KSA     1    A group    -0.7108    -0.8638
## 3     0  EGY     2    A group    -0.2066    -0.4438
## 4     1  URU     2    A group     0.2372     0.4438
## 5     3  RUS     3    A group     0.1531     0.3597
## 6     1  EGY     3    A group    -0.2066    -0.3597

For now, we will focus on the goals variable only. A brief summary yields

summary(FIFA2018$goals)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0     0.0     1.0     1.3     2.0     6.0

showing that the teams scored between \(0\) and \(6\) goals per match with an average of \(\bar y = 1.3\) from the observations \(y_i\) (\(i = 1, \dots, 128\)). The corresponding table of observed relative frequencies is:

observed <- proportions(table(FIFA2018$goals))
observed
## 
##        0        1        2        3        4        5        6 
## 0.257812 0.375000 0.250000 0.078125 0.015625 0.015625 0.007812

This confirms that goals are relatively rare events in a soccer game with each team scoring zero to two goals per match in almost 90 percent of the matches. Below we show that this observed frequency distribution can be approximated very well by a Poisson distribution which can subsequently be used to obtain predicted probabilities for the goals scored in a match.

Basic fitted distribution

In a first step, we simply assume that goals are scored with a constant mean over all teams and matches and hence just fit a single Poisson distribution for the number of goals. To do so, we obtain a point estimate of the Poisson parameter by using the empirical mean \(\hat \lambda = \bar y = 1.3\) and set up the corresponding distribution object:

p_const <- Poisson(lambda = mean(FIFA2018$goals))
p_const
## [1] "Poisson distribution (lambda = 1.3)"

In the technical details below we show that this actually corresponds to the maximum likelihood estimation for this distribution. It could also be fitted via fit_mle(Poisson(1), FIFA2018$goals) in distributions3.

As already illustrated above, the expected probabilities of observing counts of \(0, 1, \dots, 6\) goals for this Poisson distribution can be extracted using the pdf() method. A comparison with the observed empirical frequencies yields

expected <- pdf(p_const, 0:6)
cbind(observed, expected)
##   observed expected
## 0 0.257812 0.273385
## 1 0.375000 0.354546
## 2 0.250000 0.229901
## 3 0.078125 0.099384
## 4 0.015625 0.032222
## 5 0.015625 0.008358
## 6 0.007812 0.001806

By and large, all observed and expected frequencies are rather close. However, it is not reasonable that all teams score goals with the same probabilities, which would imply that winning or losing could just be attributed to “luck” or “random variation” alone. Therefore, while a certain level of randomness will certainly remain, we should also consider that there are stronger and weaker teams in the tournament.

Poisson regression and probabilistic forecasting

To account for different expected performances from the teams in the 2018 FIFA World Cup, the FIFA2018 data provides an estimated logability for each team. These have been estimated by Zeileis et al. (2018) prior to the start of the tournament (2018-05-20) based on quoted odds from 26 online bookmakers using the bookmaker consensus model of Leitner et al. (2010). The difference in logability between a team and its opponent is a useful predictor for the number of goals scored.

Consequently, we fit a generalized linear model (GLM) to the data that links the expected number of goals per team/match \(\lambda_i\) to the linear predictor \(x_i^\top \beta\) with regressor vector \(x_i^\top = (1, \mathtt{difference}_i)\) and corresponding coefficient vector \(\beta\) using a log-link: \(\log(\lambda_i) = x_i^\top \beta\). The maximum likelihood estimator \(\hat \beta\) with corresponding inference, predictions, residuals, etc. can be obtained using the glm() function from base R with family = poisson:

m <- glm(goals ~ difference, data = FIFA2018, family = poisson)
summary(m)
## 
## Call:
## glm(formula = goals ~ difference, family = poisson, data = FIFA2018)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -2.144  -1.155  -0.175   0.528   2.327  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.2127     0.0813    2.62   0.0088 ** 
## difference    0.4134     0.1058    3.91  9.3e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 144.20  on 127  degrees of freedom
## Residual deviance: 128.69  on 126  degrees of freedom
## AIC: 359.4
## 
## Number of Fisher Scoring iterations: 5

Both parameters can be interpreted. First, the intercept corresponds to the expected log-goals per team in a match of two equally strong teams, i.e., with zero difference in log-abilities. The corresponding prediction for the number of goals can either be obtained manually from the extracted coef() by applying exp() (as the inverse of the log-link).

lambda_zero <- exp(coef(m)[1])
lambda_zero
## (Intercept) 
##       1.237

Or equivalently the predict() function can be used with type = "response" in order to get the expected \(\hat \lambda_i\) (rather than just the linear predictor \(x_i^\top \hat \beta\) that is predicted by default).

predict(m, newdata = data.frame(difference = 0), type = "response")
##     1 
## 1.237

As above, we can also set up a Poisson() distribution object and obtain the associated expected probability distribution for zero to six goals in a mathc of two equally strong teams:

p_zero <- Poisson(lambda = lambda_zero)
pdf(p_zero, 0:6)
## [1] 0.290242 0.359041 0.222074 0.091571 0.028319 0.007006 0.001445

Note that distributions3 also provides a convenience function prodist() that allows to obtain p_zero in a single step via prodist(m, newdata = data.frame(difference = 0)).

Second, the slope of \(0.413\) can be interpreted as an ability elasticity of the number of goals scored. This is because the difference of the log-abilities can also be understood as the log of the ability ratio. Thus, when the ability ratio increases by \(1\) percent, the expected number of goals increases approximately by \(0.413\) percent.

This yields a different predicted Poisson distribution for each team/match in the tournament. We can set up the vector of all \(128\) Poisson() distribution objects by extracting the vector of all fitted point estimates \((\hat \lambda_1, \dots, \hat \lambda_{128})^\top\):

p_reg <- Poisson(lambda = fitted(m))
length(p_reg)
## [1] 128
head(p_reg)
##                                       1                                       2 
## "Poisson distribution (lambda = 1.768)" "Poisson distribution (lambda = 0.866)" 
##                                       3                                       4 
## "Poisson distribution (lambda = 1.030)" "Poisson distribution (lambda = 1.486)" 
##                                       5                                       6 
## "Poisson distribution (lambda = 1.435)" "Poisson distribution (lambda = 1.066)"

Again, the convenience function prodist(m) could also be used to directly extract p_reg.

Note that specific elements from the vector p_reg of Poisson distributions can be extracted as usual, e.g., with an index like p_reg[i] or using the head() and tail() functions etc.

As an illustration, the following goal distributions could be expected for the FIFA World Cup final (in the last two rows of the data) that France won 4-2 against Croatia:

tail(FIFA2018, 2)
##     goals team match  type    stage logability difference
## 127     4  FRA    64 Final knockout     0.8866      0.629
## 128     2  CRO    64 Final knockout     0.2576     -0.629
p_final <- tail(p_reg, 2)
p_final
##                                     127                                     128 
## "Poisson distribution (lambda = 1.604)" "Poisson distribution (lambda = 0.954)"
pdf(p_final, 0:6)
##        d_0    d_1    d_2     d_3     d_4      d_5       d_6
## 127 0.2010 0.3225 0.2587 0.13836 0.05550 0.017808 0.0047618
## 128 0.3853 0.3675 0.1752 0.05572 0.01329 0.002534 0.0004029

This shows that France was expected to score more goals than Croatia but both teams scored more goals than expected, albeit not unlikely many.

Further details and extensions

Assuming independence of the number of goals scored, we can obtain the table of possible match results (after normal time) by multiplying the marginal probabilities (again only up to six goals). In R this be done using the outer() function which by default performs a multiplication of its arguments.

res <- outer(pdf(p_final[1], 0:6), pdf(p_final[2], 0:6))
round(100 * res, digits = 2)
##       [,1]  [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]  7.74  7.39 3.52 1.12 0.27 0.05 0.01
## [2,] 12.43 11.85 5.65 1.80 0.43 0.08 0.01
## [3,]  9.97  9.51 4.53 1.44 0.34 0.07 0.01
## [4,]  5.33  5.08 2.42 0.77 0.18 0.04 0.01
## [5,]  2.14  2.04 0.97 0.31 0.07 0.01 0.00
## [6,]  0.69  0.65 0.31 0.10 0.02 0.00 0.00
## [7,]  0.18  0.17 0.08 0.03 0.01 0.00 0.00

For example, we can see from this table that the expected probability for France winning against Croatia 1-0 is \(12.43\) percent while the probability that France loses 0-1 is only \(7.39\) percent.

The advantage of France can also be brought out more clearly by aggregating the probabilities for winning (lower triangular matrix), a draw (diagonal), or losing (upper triangular matrix). In R these can be computed as:

sum(res[lower.tri(res)]) ## France wins
## [1] 0.5245
sum(diag(res))           ## draw
## [1] 0.2498
sum(res[upper.tri(res)]) ## France loses
## [1] 0.2243

Note that these probabilities do not sum up to \(1\) because we only considered up to six goals per team but more goals can actually occur with a small probability.

Next, we update the expected frequencies table by averaging across the expectations per team/match from the regression model.

expected <- pdf(p_reg, 0:6)
head(expected)
##      d_0    d_1    d_2     d_3     d_4      d_5       d_6
## 1 0.1707 0.3017 0.2667 0.15721 0.06949 0.024571 0.0072403
## 2 0.4208 0.3642 0.1576 0.04548 0.00984 0.001703 0.0002457
## 3 0.3571 0.3677 0.1893 0.06498 0.01673 0.003444 0.0005911
## 4 0.2262 0.3362 0.2498 0.12377 0.04599 0.013669 0.0033857
## 5 0.2380 0.3417 0.2452 0.11732 0.04210 0.012086 0.0028914
## 6 0.3444 0.3671 0.1957 0.06954 0.01853 0.003952 0.0007022
expected <- colMeans(expected)
cbind(observed, expected)
##   observed expected
## 0 0.257812 0.294374
## 1 0.375000 0.340171
## 2 0.250000 0.214456
## 3 0.078125 0.098236
## 4 0.015625 0.036595
## 5 0.015625 0.011727
## 6 0.007812 0.003333

As before, observed and expected frequencies are reasonably close, emphasizing that the model has a good marginal fit for this data. To bring out the discrepancies graphically we show the frequencies on a square root scale using a so-called hanging rootogram (Kleiber & Zeileis 2016). The gray bars represent the square-root of the observed frequencies “hanging” from the square-root of the expected frequencies in the red line. The offset around the x-axis thus shows the difference between the two frequencies which is reasonably close to zero.

bp <- barplot(sqrt(observed), offset = sqrt(expected) - sqrt(observed),
  xlab = "Goals", ylab = "sqrt(Frequency)")
lines(bp, sqrt(expected), type = "o", pch = 19, lwd = 2, col = 2)
abline(h = 0, lty = 2)

Rootogram for the number of goals in the 2018 FIFA World Cup modeled by a Poisson regression model

Finally, we want to point out that while the log-abilities (and thus their differences) had been obtained based on bookmakers odds prior to the tournament, the calibration of the intercept and slope coefficients was done “in-sample”. This means that we have used the data from the tournament itself for estimating the GLM and the evaluation above can only be made ex post. Alternatively, one could have used previous FIFA World Cups for calibrating the coefficients so that probabilistic forecasts for the outcome of all matches (and thus the entire tournament) could have been obtained ex ante. This is the approach used by Groll et al. (2019) and Groll et al. (2021) who additionally added further explanatory variables and used flexible machine learning regression techniques rather than a simple Poisson GLM.

Technical details: Maximum likelihood estimation of \(\lambda\)

Fitting a single Poisson distribution with constant \(\lambda\) to \(n\) independent observations \(y_1, \dots, y_n\) using maximum likelihood estimation can be done analytically using basic algebra. First, we set up the log-likelihood function \(\ell\) as the sum of the log-densities per observation:

[\begin{align} \ell(\lambda; y_1, \dots, y_n) & = \sum_{i = 1}^n \log f(y_i; \lambda)
\end{align
}]

For solving the first-order condition analytically below we need the score function, i.e., the derivative of the log-likelihood with respect to the parameter \(\lambda\). The derivative of the sum is simply the sum of the derivatives:

[\begin{align} \ell^\prime(\lambda; y_1, \dots, y_n) & = \sum_{i = 1}^n \left{ \log f(y_i; \lambda) \right}^\prime
& = \sum_{i = 1}^n \left{ -\lambda + y_i \cdot \log(\lambda) - \log(y_i!) \right}^\prime
& = \sum_{i = 1}^n \left{ -1 + y_i \cdot \frac{1}{\lambda} \right}
& = -n + \frac{1}{\lambda} \sum_{i = 1}^n y_i \end{align
}]

The first-order condition for maximizing the log-likelihood sets its derivative to zero. This can be solved as follows:

[\begin{align} \ell^\prime(\lambda; y_1, \dots, y_n) & = 0
-n + \frac{1}{\lambda} \sum_{i = 1}^n y_i & = 0
n \cdot \lambda & = \sum_{i = 1}^n y_i
\lambda & = \frac{1}{n} \sum_{i = 1}^n y_i = \bar y \end{align
}]

Thus, the maximum likelihood estimator is simply the empirical mean \(\hat \lambda = \bar y.\)

Unfortunately, when the parameter \(\lambda\) is not constant but depends on a linear predictor through a log link \(\log(\lambda_i) = x_i^\top \beta\), the corresponding log-likelihood of the regression coefficients \(\beta\) can not be maximized as easily. There is no closed-form solution for the maximum likelihood estimator \(\hat \beta\) which is why the glm() function employs an iterative numerical algorithm (so-called iteratively weighted least squares) for fitting the model.

References

  • Groll A, Hvattum LM, Ley C, Popp F, Schauberger G, Van Eetvelde H, Zeileis A (2021). “Hybrid Machine Learning Forecasts for the UEFA EURO 2020.” arXiv 2106.05799. arXiv.org E-Print Archive. https://arxiv.org/abs/2106.05799
  • Groll A, Ley C, Schauberger G, Van Eetvelde H (2019). “A Hybrid Random Forest to Predict Soccer Matches in International Tournaments.” Journal of Quantitative Analysis in Sports 15(4), 271-87. https://doi.org/10.1515/jqas-2018-0060
  • Hayes A, Moller-Trane R, Jordan D, Northrop P, Lang M, Zeileis A (2022). “distributions3: Probability Distributions as S3 Objects.” R package version 0.2.0, https://CRAN.R-project.org/package=distributions3
  • Kleiber C, Zeileis A (2016). “Visualizing Count Data Regressions Using Rootograms.” The American Statistician 70(3), 296-303. https://doi.org/10.1080/00031305.2016.1173590
  • Leitner C, Zeileis A, Hornik K (2010). “Forecasting Sports Tournaments by Ratings of (Prob)abilities: A Comparison for the EURO 2008.” International Journal of Forecasting 26(3), 471-81. https://doi.org/10.1016/j.ijforecast.2009.10.001
  • Wikipedia (2022). “Poisson Distribution - Wikipedia, the Free Encyclopedia.” https://en.wikipedia.org/wiki/Poisson_distribution, accessed 2022-02-21.
  • Zeileis A, Leitner C, Hornik, K (2018). “Probabilistic Forecasts for the 2018 FIFA World Cup Based on the Bookmaker Consensus Model.” Working Paper 2018-09. Working Papers in Economics & Statistics, Research Platform Empirical & Experimental Economics, Universität Innsbruck. https://EconPapers.RePEc.org/RePEc:inn:wpaper:2018-09
To leave a comment for the author, please follow the link and comment on their blog: Achim Zeileis.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)