In this article, I test the fitness of Poisson and Poisson-Gamma stochastic models to the provided data. The data consist of the counts of deaths, which occurred in each one of 60 Game of Thrones episodes. The goal is to find models’ parameters that maximise the maximum likelihood.
Poisson model
Let \(y_i \sim \text{Poisson}(\lambda)\), where \(y_i\) is the number of deaths which occurred in episode \(i\). Then:
\[\begin{align} f(y_i) &= \frac{\lambda^{y_i}}{y_i!} e^{-\lambda} \\ \text{Likelihood: }&\prod_i f(y_i) = \frac{\lambda^{\Sigma y_i}}{\Pi y_i!}e^{-n\lambda} \\ \text{Loglikelihood: }&L = \sum_i y_i \log\lambda - \sum_i \log(y_i!) -n\lambda \\ \text{Maximise Loglikelihood: }& \frac{\partial L}{\partial \lambda} = \frac{\sum_i y_i}{\lambda} - n = 0 \\ &\hat{\lambda} = \frac{\sum y_i}{n} = \bar{y} \\ \text{Second derivative: }& \frac{\partial^2 L}{\partial \lambda^2} = -\frac{\sum y_i}{\lambda^2} > 0 \text{ (i.e. maximum)}\\ \text{Fisher information: }& \Rightarrow I(\hat\lambda) = \frac{\sum y_i}{\lambda^2} = \frac{n^2}{\sum y_i} \\ \text{Standard error: }& \Rightarrow \text{SE}(\hat\lambda) = \sqrt{I(\hat\lambda)^{-1}} \end{align}\]
By plugging in the observed data, we get the maximum likelihood estimate \(\hat\lambda = 3.27\) and \(\text{SE}(\hat\lambda)=0.23\). The results are exactly the same as those achieved by the optim function:
lambda= 3.266667
SE(lambda)= 0.2333333 
Because \(E[Y]=Var[Y]=\lambda\), then this model may not be able to explain all the variance within the data.
Poisson-Gamma model
Derivation of closed-form likelihood for Poisson-Gamma model is challenging. Instead, we can ask optim function to use its variation of the Newton-Raphson method (BFGS) to get the numerical approximation to the parameters (\(\alpha\) and \(\beta\)), which maximise the likelihood.
alpha= 4.386928 , SE(alpha)= 1.920489
beta= 0.7447539 , SE(beta)= 0.3335913 
Because this model’s expected value and variance are different functions of its parameters, it can explain more variance than the first model.
Model Comparision
The Poisson model cannot model the variance in the data as good as the Poisson-Gamma model, but it has only one parameter, i.e. \(\lambda\). Poisson-Gamma model has two parameters, i.e. \(\alpha\) and \(\beta\), and the chances for it to overfit the data is greater. To compare these two models, we can use, e.g. Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC) or leave-one-out score.
| Leave.one.out | AIC | BIC | |
|---|---|---|---|
| Poisson | -133.6948 | -273.9216 | -276.0159 |
| Poisson-Gamma | -128.0138 | -264.3671 | -268.5558 |
The Poisson-Gamma model has better (i.e. greater) scores than the Poisson model, despite its increased complexity. It is considered a better model for the observed data.
Goodness of fit
To assess the goodness of fit for the data, I propose sampling model parameters from their 95% CI ranges and calculating the predicted frequencies of deaths using these parameters. Then, we can inspect if the observed data lay within the 95% CI of the predictions.
Only one observation is located outside 95% CI in the Poisson-Gamma model. The vanilla Poisson model has eight such observations. Thanks to the increased degree of freedom in modelling the variance, the former model is a better fit to the observed data.
Conclusion
According to the Poisson model, the expected value of deaths in an episode of the Game of Thrones is 3.27 with SE=1.81. The Poisson-Gamma model estimates the expected value to be nearly the same as in the Poisson model, i.e. 3.27, but thanks to the additional parameter, it can model the variance to fit better the observed data (SE=5.7). The model comparison shows that the Poisson-Gamma model gets better information criteria scores despite its increased complexity. The goodness of fit plots show that the observed data is located within 95% CI of predictions for the Poisson-Gamma model. The Poisson model, with it’s reduced variability fails to include observed data in its confidence intervals more often. Both models fail to predict a large number of deaths in an episode. According to data, the chance of observing 11 deaths in an episode is 3.33%. The Poisson and Poisson-Gamma models predict it to be 0.04% and 0.45% respectively.