I need help choosing between the binomial, negative binomial, and beta-binomial distributions and understanding the rationale of why I would choose one over the other, as well as understanding how it changes my interpretation of coefficients in the model. The most helpful links I've found have been: "Negative binomial vs binomial distribution for proportion data" (asked by Joris) and "Diagnostic plots for count regression" (asked by half-pass). Both were helpful, but I'm still quite uncertain.
Data: Cross-sectional, self-report. We gave 222 veterans using opioids a previously validated scale of 21 questions asking about days of opioid overdose risk behaviors in the last 30 days (e.g. "How many days in the last 30 did you use heroin?"). Disregarding the problems with the survey items, we assume they all represent equally risky opioid-related overdose risk behaviors and take a mean score across the items for a final Risk Behaviors score from 0-30. Huge assumptions, but that's what the authors want to do.
The item-level distributions mostly look like this:
How many days did you take opioid pain medicine that you got from some source other than your own doctor’s prescription (for example, that you bought from someone or that were given to you by someone)? (note that I've modeled the variables with 0s on left and without 0s on the right)
Which empirical distribution best reflects the observed distribution? Why?
Should I select the model by comparing the different regressions, using AIC (or something else) as my criterion?
There is an upper limit of 30, which might suggest the need for a binomial or beta binomial distribution, but the observed data only have a max of 23, and there's only 5 people with scores above 15 (not a meaningful cutoff, just what I see above). They may be outliers.
There is extreme overdispersion (Mean=3.1, Var=15.1), so no Poisson, and probably not binomial.
I have compared negative binomial and binomial regression below:
a. negative binomial regression
glm(SumofDV/#ofitems/#ofdaysinamonth ~ ., df, family='binomial', weights=#ofdaysinamonth)
b. binomial regression
and I draw approximately the same conclusions (for individual t-tests of the variable-outcome relationship).
(I haven't figured out how to implement the beta-binomial model in r because, if I take the mean of event frequency, I have non-negative integer data and it doesn't have an easy way to handle that. In the other models (bin,nnegbin), I can use an offset or weights to get the average and make the analysis work.)