I have the following data structure (sorry, I can't share actual data due to privacy concerns).
A dataframe with:
counts category exposure
There are many samples, each having a category, a number for exposure (e.g. number of attempts), and a number of counts (e.g. successes)
Additionally, the data is quite zero-inflated.
I fit the following two models, for example in R code:
A Poisson model:
glm('count ~ category + offset(log(exposure))', family = poisson())
and a Negative binomial model:
glm.nb('count ~ category + offset(log(exposure))')
By all measures I looked at, AIC, looking at deviance, etc the negative binomial model is a much better fit.
However, I'm having trouble with some of the coefficients its generating.
For a given level of the categorical variable, 'B', and where 'A' is the baseline level, if I understand the coefficients correctly, I expect the fitted coefficient, C, for level 'B' to satisfy the following equation:
exp(C) = (sum(counts[category==B])*sum(exposure[category==A]))/(sum(counts[category==A])*sum(exposure[category==B])
So basically the ratio of B's effect size to the baseline level A.
This comes out very accurately with the poisson model, which is nevertheless the worse fit. On the other hand, the negative binomial model, which is a much better fit, does not produce coefficients that are correct, and 2 of the 15 or so coefficients are not in the right order of size, relative to my expectation when I do
sum(counts[category==X])/sum(exposure[category==X]) for a given category.
To me this suggests that basically the Negaive Binomial model is fitting the individual samples better, thereby producing less error in goodness-of-fit tests, yet simultaneously is not predicting the 'mean' effect as well. Why might this be the case? What should I check? And which would be the 'better' model?
I made a plot using this code, and where
dat is a dataframe holding my data:
plot(dat$exposure[dat$category == 'B'], dat$count[dat2$category == 'B']) lines(dat$exposure[dat$category=='B'], predict(mPoisson,dat[dat$category == 'B',], type = 'response')) lines(dat$exposure[dat$category=='B'], predict(mNB,dat[dat$category == 'B',], type = 'response'), col = 'red')
So for a single category level, I plot exposure vs counts, and then for the two models I plot exposure vs predicted counts for that same category.
It turns out that this relationship is not at all linear but more like counts ~= exposure^(1.5) (very roughly, but you get the idea)
The NB model has less error in goodness-of-fit tests and in the individuals because it more closely fits the points when exposure is low, and I have more of those points. The poisson is more affected by these fewer but 'large exposure' points. In fact, the high exposure points are more influential in determining the
sum(counts)/sum(exposure) quantity I'm interested in, and that is why the poisson is doing such a good job. It also explains why the zero-inflated poisson, although an improvement on poisson, does not beat the negative binomial, since it only accounts for low exposure points where the count is 0, but NOT low exposure points where the count is not 0.
I guess the best model should be something where the offset(log(exposure)) is not a linear term. Is there a standard way of modeling this such that the coefficient you get back will still be the mean? I'm pretty sure the answer is 'no'. I think this basically points to the model that best fits the individuals not being a good fit for the 'mean' due to the shape of the data, so I should probably just pick a model that best fits the mean and ignore the AIC and error on the individuals.
Any feedback on this and whether this is the right track is appreciated.
I think what's going on is also related to this:
In their data they report that negative binomial fits the lower mean observations better and my observations from various plots is the same for me - that the negative bionomial fits the lower rates better, ignoring the larger rates. Poisson models (including quasi-poisson) seem to fit the low and high rates in such a way to produce exactly the coefficients I expect, but the lower rates seem to pull the negative binomial away from those estimates.
I guess in that sense poisson-based models are better for this data and question in spite of AIC measures.
I'm still unsure of whether this holds up when I start throwing in variables I need to control for such that the regression will no longer only consist of orthogonal variables. Perhaps I can use some sort of stratification to evaluate model fits vs expected effect sizes.