# Appropriate application of Poisson regression?

I am looking to predict the number of shots that are going to be made by soccer players in their next game.

Currently, I use multiple linear regression, where I regress my entire data set (shotdata) for all c.300 players at once, with the following explanatory variables

av = average (mean) shots in last 5 games

po = possession during the games

ha = dummy variable whether the player is playing home or away

such that:

lm(shots ~ av + po + ha, data = shotdata)


I then predict one observation , for each player, after plugging in relevant values for (predicted) possession and whether the player is playing home/away.

I am looking to improve my model, however, given that shots are count data and seem to exert a Poisson distribution (see http://imgur.com/LDSZSDE):

I have tried using a Poisson regression model, such that:

glm(shots ~ av + po + ha, data = shotdata, family = “poisson”)


On an initial, layman, inspection of the results, however, this model does not seem to offer an improvement in prediction. The response residuals it is generating are about 15% higher than those from the multiple linear regression and simple eyeballing of the predictions suggests they are not as strong as those from the multiple linear regression.

Might this not be an appropriate application of the Poisson regression? Is there something I’m obviously doing wrong with the Poisson distribution? Are there other methods for determining the comparative strengths of the model I should look at?

I've tried to keep this question concise with relevant information only, but please let me know if I can expand on any of it.

• What metric are you using to assess model quality? What is your loss function? (AKA, what do you want to achieve with your prediction?) Are you doing out-of-sample testing, or cross-validation? – Stephan Kolassa Sep 30 '15 at 16:06
• When you generate the predictions, are you using type='response' in your call to predict? (I've left that out before myself, much to my subsequent confusion, hence my question.) – jbowman Sep 30 '15 at 17:24
• @StephanKolassa To be honest, I'm not doing out of sample testing at the moment, as I feel my data set is currently too small. But maybe that's flawed thinking? I also just looked into cross-validation and it seems that that is better at dealing with small datasets than conventional validation (70/30 partitioning). My current 'validation' is just to look at RMSE and to eyeball the predictions it is giving me for the next soccer match. The Poisson based predictions seem to generate more obvious bad predictions, particularly for players with extremely high or extremely low av values. – Will T-E Oct 1 '15 at 10:16
• @jbowman yup thanks, I am, although I initially made that mistake too! A quicker way around it I found is just to do modelname\$fitted as opposed to predict(modelname, type = "response") – Will T-E Oct 1 '15 at 10:18

A couple of thoughts that may or may not help below. It's a bit hard for us to be helpful without seeing your actual data...

• First of all, your data rather obviously does not follow a standard regression model form: observations are integer, residuals will certainly not be normally distributed, and so forth. This is a textbook case of where to use count data models.

So if you find that an OLS model performs better than a count data model, something seems to be badly broken.

• I'd strongly suggest running lots of diagnostics for both models, which should show you ways in which you could improve either model. Plot your raw shots response against your predictors. Plot your predicted response against the predictors. Plot your residuals against the predictors, and plot actuals and residuals against predictions. Look at prediction distributions, even if you are most interested in point predictions - your bad point forecasts may be due to high and unmodeled variance. If so, you may want to think about this in the context of your specific application for predictions. Unfortunately, diagnostics for count data models are not as well researched as for OLS models.

The vignette "Regression models for count data in R" for the pscl package may be helpful here. Note that you can also fit OLS models with glm(), by specifying family = gaussian, which is actually the default, so you could directly use many of the summaries explained in the vignette.

• You write in a comment that the fit is particularly bad for high or low av values. This to me suggests transforming av, for instance using splines. Look at Frank Harrell's Regression Modeling Strategies, which has just come out in a brand new second edition, with an accompanying R package called rms.

• While we are discussing the model form, do you have information on the position your players play? A striker will have a much higher number of shots on goal (which I assume you are modeling, not simply shots and passes as such) than a defender.

• I'm already happy you are using the RMSE and not the MAD, which would be minimized by the conditional median...

• I am not overly surprised that the OLS model yields a better in-sample fit than a Poisson regression - after all, it has one more degree of freedom to play around with, namely the variance. To "even the playing field", as it were, you may want to offer your count data model another parameter, as well. For instance, try looking at a negative binomial regression. Or a hurdle or a zero-inflated model (see the vignette linked above).

• In-sample fit is a notoriously misleading proxy for out-of-sample predictive accuracy, especially for models with varying numbers of predictors. Even with only 300 data points, I'd suggest you at least run cross-validation, say five-fold. This would also give you an idea of the variability of your prediction error. It's quite possible that your prediction error for both models will jump around erratically.

• @StephenKolassa Many thanks, I will work my way through your recommendations and report my findings back in this thread. And yes, I do have data on position, and I have actually be clustering my data-set based on position, and then performing regression on each cluster. This is one of the details I omitted from the question, in the interests of keeping the question concise. – Will T-E Oct 1 '15 at 16:05
• Separate analyses is one possibility. Alternatively, build one big model containing all players, with indicator variables for position - potentially in interactions with your other predictors. – Stephan Kolassa Oct 1 '15 at 16:25
• @StephenKolassa I've just been thinking about doing that! I can't find a discussion anywhere of it though, of why you should do one rather than the other (i.e. interaction effects vs clustering). Do you know anywhere I can find such a discussion? Or does it just come down to running diagnostics, again, as per your second point above? – Will T-E Oct 1 '15 at 16:29
• If you do separate analyses, this amounts to having different coefficients for each predictor within each cluster. Interaction terms in a big model will yield the same. One difference lies in whether residual variance is homogeneous between clusters - if you use one big model, this will typically assume homogeneous variance (although you can model heteroskedasticity in principle). A good guideline would be: run cross-validations for both approaches and check whether one consistently yields better out-of-bag predictions. – Stephan Kolassa Oct 1 '15 at 16:34