poisson models for proportions
The trick to this problem is that it is more a proportion than it is a count. You may have been told that a Poisson model is "used for count data". In fact, you could use Poisson data to model proportions (or binomial) outcomes. The only issue is that the Poisson model tends to over estimate the variance, (Binomial: var = p(1-p), Poisson: var = p) leading to wider confidence intervals and larger p-values, on average. It can also predict beyond the range of normal values. Is that a bad thing? Just truncate the over predictions, and wide confidence intervals just mean reduced power.
independent data models for over/under-dispersed binomial outcomes
Diving into the problem more, binomial variables come from independent and identically distributed Bernoulli events. It's important to verify those assumptions. One way is to inspect the dispersion of Bernoulli outcomes. For instance, let's simulate 3 meetings and their attendance, but 1st not so popular, and 3rd is very popular.
p <- c(0.3, 0.5, 0.7)
x <- matrix(rbinom(100*3, 1, p), 100, 3, byrow = T)
y <- rowSums(x)
summary(glm(cbind(y, 3) ~ 1, family=quasibinomial))
glm(formula = cbind(y, 3) ~ 1, family = quasibinomial)
Min 1Q Median 3Q Max
-1.5555 -0.3563 0.3188 0.3188 0.8486
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.69984 0.05662 -12.36 <2e-16 ***
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for quasibinomial family taken to be 0.3191032)
Null deviance: 41.543 on 99 degrees of freedom
Residual deviance: 41.543 on 99 degrees of freedom
Number of Fisher Scoring iterations: 4
A dispersion parameter of 0.32 which is highly under-dispersed. Another way to check in this univariate case is simply:
Note: 0.08 (actual variance) / 0.25 (expected variance under binomial probability model) = 0.32. I consider the
glm approach superior because it allows for adjusting for covariates.
One way of conceptualizing the problem of underdispersion (in this case) is that if you simulated outcomes from a binomial GLM, their distributions would not match the target distributions. The calibration of the probability model for the response is:
> # expected
> dbinom(x = 0:3, size = 3, prob = mean(y/3))
 0.1275167 0.3774832 0.3724834 0.1225166
> # observed
0 1 2 3
0.11 0.37 0.44 0.08
(see figure at bottom of the page for a calibration plot)
Here, using a quasibinomial GLM will naturally widen confidence intervals to compensate for the miscalibration noted above. That way we can get correct p-values and confidence intervals. Predicting and simulating outcomes remains a problem, however.
proportional odds models for over/under-dispersed binomial outcomes.
Another approach is to use a proportional odds model. These "cumulative link" models allow for adjacent counts to have a unique latent log-odds difference determining the likelihood of achieving at least one higher outcome. This is why this type of model is preferred for ordinal data.
A "novariate" (intercept only) model for the response gives the following log-odds parameters for endorsing an at-least-one-higher response
polr(factor(y) ~ 1)
polr(formula = factor(y) ~ 1)
0|1 1|2 2|3
-2.09049882 -0.07990206 2.44238850
Residual Deviance: 234.7927
Importantly, this leads to a probability model with a much higher agreement with the empirical frequencies (calibration)
> predict(polr(factor(y) ~ 1), newdata = data.frame('1'=1), type = 'probs')
0 1 2 3
0.11002372 0.37001138 0.43996795 0.07999695
So I find that this method is superior for handling the case when the Bernoulli outcomes contributing the binomial event are possibly non-independent and possible non-identical in their probability.