Visually compare binomial and beta-binomial appropriateness in R

Question

I'm trying to choose the best distribution family for generalized linear regression. My outcome is cross-sectional, over-dispersed proportion data (# of behaviors/20-22 possible behaviors). I used the functions fitdistr and plotdist from the package fitdistrplus to visually rule out the binomial and poisson distributions, so now I'm choosing between the beta-binomial, quasibinomial, and negative binomial distributions.

Note: I don't really understand the beta-binomial distribution, but a lot of searching on this site suggested that it's often most appropriate for overdispersed proportion data.

Key info: I'm inclined to select the negative binomial distribution because

(1) I understand glm.nb better.

(2) glm.nb has the lowest AIC, provides similar results to both the beta-binomial and quasibinomial analyses, and my outcome closely aligns visually with the theoretical negbin distribution (see graphs below).

Questions:

1. The fitdistrplus package does not support the betabinomial distribution, so I don't know how to visually compare my data with its theoretical distribution. Is it possible to visually compare the binomial and beta-binomial distributions? How would I do that in R?
2. What are the pros and cons of selecting the negative binomial in this case? I read in Wagner (2015) that the negative binomial is supposed to approximate the beta-binomial distribution. What does that mean for my analysis?

References

Wagner, B., Riggs, B., & Mikulich-Gilbertson, S. 2015. The importance of distribution-choice in modeling substance use data: a comparison of negative binomial, beta binomial, and zero-inflated distributions. American Journal of Drug and Alcohol Abuse, 41(6), 489–497.

Answer 1

It seems from the comments that (on the ground of distribution support) the alternatives are narrowed down to binomial or beta-binomial, so I will consider that.

First, the beta-binomial distribution is a compound distribution. While in a binomial experiment, the probability parameter $p$ is the same for every draw, in the beta-binomial, it is not. Before every draw, the $p$ is itself drawn from a beta distribution, independently each time. That gives rise to over-dispersion, as you have noted.

One way of doing comparisons, is to plot the data, together with simulations from a binomial model, a few times. An example is in the following post: Independence of sexes of children born to a couple

In that example, there are no covariables, so it s simple, but the same idea can be used with more complex models.