I'm cross validating a model that's trying to predict a count. If this was a binary classification problem, I'd calculate out-of-fold AUC, and if this was a regression problem I'd calculate out-of-fold RMSE or MAE.

For a Poisson model, what error metrics can I use to evaluate the "accuracy" of the out-of-sample predictions? Is there a Poisson extension of AUC that looks at how well the predictions order the actual values?

It seems that a lot of Kaggle competitions for counts (e.g. number of useful votes a yelp review will get, or number of days a patient will spend in the hospital) use root mean log squared error, or RMLSE.

/Edit: One thing I've been doing is calculating deciles of the predicted values, and then looking at the actual counts, binned by decile. If decile 1 is low, decile 10 is high, and the deciles in between are strictly increasing, I've been calling the model "good," but I've been having trouble quantifying this process, and I've convinced there's a better approach.

/Edit 2: I'm looking for a formula that takes predicted and actual values and returns some "error" or "accuracy" metric. My plan is to calculate this function on the out-of-fold data during cross-validation, and then use it to compare a wide variety of models (e.g. a poisson regression, a random forest and a GBM).

For example, one such function is RMSE = sqrt(mean((predicted-actual)^2)). Another such function would be AUC. Neither function seems to be right for poisson data.

  • 1
    $\begingroup$ For Poisson models you could use the deviance; which is akin to a MSE but better suited to a Poisson. If you samples sizes aren't small, a weighted MSE would be quite similar. $\endgroup$ – Glen_b Oct 2 '13 at 22:22
  • $\begingroup$ @Glen_b what's the formula for deviance? $\endgroup$ – Zach Oct 2 '13 at 23:54
  • 1
    $\begingroup$ Deviance. How are you fitting your Poisson model? $\endgroup$ – Glen_b Oct 3 '13 at 0:02
  • $\begingroup$ A couple of different ways, ranging from a penalized poisson regression to a gbm. I'm looking for a good error metric to compare different models. Thanks for the advice. $\endgroup$ – Zach Oct 3 '13 at 2:59
  • $\begingroup$ Poisson regression at least should give you a deviance automatically $\endgroup$ – Glen_b Oct 3 '13 at 3:02

There are a couple of proper and strictly proper scoring rules for count data you can use. Scoring rules are penalties $s(y,P)$ introduced with $P$ being the predictive distribution and $y$ the observed value. They have a number of desirable properties, first and foremost that a forecast that is closer to the true probability will always receive less penalty and there is a (unique) best forecast and that one is when the predicted probability coincides with the true probability. Thus minimizing the expectation of $s(y,P)$ means reporting the true probabilities. See also Wikipedia.

Often one takes an average of those over all predicted values as

$S=\frac{1}{n}\sum_{i=1}^n s(y^{(i)},P^{(i)})$

Which rule to take depends on your objective but I'll give a rough characterization when each is good to be used.

In what follows I use $f(y)$ for the predictive probability mass function $\Pr(Y=y)$ and $F(y)$ the predictive cumulative distribution function. A $\sum_k$ runs over the whole support of the count distribution (i.e, $0,1,\dots, \infty$). $I$ denotes an indicator function. $\mu$ and $\sigma$ are the mean and standard deviation of the predictive distribution (which are usually directly estimated quantities in count data models).

Strictly proper scoring rules

  • Brier Score: $s(y,P)=-2 f(y) + \sum_k f^2(k)$ (stable for size imbalance in categorical predictors)
  • Dawid-Sebastiani score: $s(y,P)=(\frac{y-\mu}{\sigma})^2+2\log\sigma$ (good for general predictive model choice; stable for size imbalance in categorical predictors)
  • Deviance score: $s(y,P)=-2\log f(y) + g_y$ ($g_y$ is a normalization term that only depends on $y$, in Poisson models it is usually taken as the saturated deviance; good for use with estimates from an ML framework)
  • Logarithmic score: $s(y,P)=-\log f(y)$ (very easily calculated; stable for size imbalance in categorical predictors)
  • Ranked probability score: $s(y,P)=\sum_k \{F(k)-I(y\leq k)\}^2$ (good for contrasting different predictions of very high counts; susceptible to size imbalance in categorical predictors)
  • Spherical score: $s(y,P)=\frac{f(y)}{\sqrt{\sum_k f^2(k)}}$ (stable for size imbalance in categorical predictors)

Other scoring rules (not so proper but often used)

  • Absolute error score: $s(y,P)=|y-\mu|$ (not proper)
  • Squared error score: $s(y,P)=(y-\mu)^2$ (not strictly proper; susceptible to outliers; susceptible to size imbalance in categorical predictors)
  • Pearson normalized squared error score: $s(y,P)=(\frac{y-\mu}{\sigma})^2$ (not strictly proper; susceptible to outliers; can be used for checking if model checking if the averaged score is very different from 1; stable for size imbalance in categorical predictors)

Example R code for the strictly proper rules:

m1 <- glm(Freq ~ mental, family=poisson, data=Mental) 

# scores for the first observation
mu <- predict(m1, type="response")[1]
x  <- Mental$Freq[1]

# logarithmic (equivalent to deviance score up to a constant) 
-log(dpois(x, lambda=mu))

# quadratic (brier)
-2*dpois(x,lambda=mu) + sapply(mu, function(x){ sum(dpois(1:1000,lambda=x)^2) })

# spherical
- dpois(x,mu) / sqrt(sapply(mu, function(x){ sum(dpois(1:1000,lambda=x)^2) }))

# ranked probability score
sum(ppois((-1):(x-1), mu)^2) + sum((ppois(x:10000,mu)-1)^2)

# Dawid Sebastiani
(x-mu)^2/mu + log(mu)
| cite | improve this answer | |
  • $\begingroup$ @Momo, it's a old thread but very good and useful. Question however about the logarithmic score. You used function -log(f(y)). Is the - sign really should be there? In your scoring rule wikipedia link (en.wikipedia.org/wiki/Scoring_rule#Logarithmic_scoring_rule), the logarithmic score as no negative sign: L(r,i)=ln(ri) is it normal? Finally, in that case is a higher score better or worst? $\endgroup$ – Bastien Nov 13 '18 at 14:24
  • $\begingroup$ Is it better (or at least more conservative and more realistic) to calculate these measures on a validation data-set that wasn't a part of the data used for estimating the models? $\endgroup$ – Fred May 17 '19 at 19:20
  • $\begingroup$ Given that GLMs are fit using iteratively reweighted least squares, as in bwlewis.github.io/GLM, what would be the objection actually of calculating a weighted R2 on the GLM link scale, using 1/variance weights as weights (which glm gives back in the slot weights in a glm fit)? This would also work for a Poisson glm right? $\endgroup$ – Tom Wenseleers Jun 11 '19 at 13:17
  • $\begingroup$ See stats.stackexchange.com/questions/412580/… for a reproducible example... $\endgroup$ – Tom Wenseleers Jun 11 '19 at 22:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.