how to analyze overdispersed binary data In my substantive research, I often use dichotomous scoring (1 correct, 0 wrong) for my tests (tests with $15~yes/no$ items). My goal is often to compare the the proportion of correct answers to all items across two groups (e.g., $n1 = n2 = 25$).
Each test-taker's sum of $15$ item scores, will be binomially distributed. But we want to compare $25$ such test-takers from group1 to $25$ other independent test-takers in group2.
What is a reasonable plan to do such a comparison (preferably via Bayesian estimation)?
Note: I believe such a data may be 'overdispersed' due to the fact that $p$ (probability of success) may not be common within each group.
Here is my data comprised of sum of $15$ $yes/no$ item scores for each person in the R code:
group1 = c(7, 11, 10, 9, 7, 11, 7, 7, 9, 9,12,11,11, 9,10, 8,10, 9, 5,10, 8, 7,11, 9, 12)

group2 = c(6, 7, 3, 7, 10, 7, 7, 8, 7, 6, 9, 8, 9, 5, 9, 7, 9, 8, 9, 8, 8,  9, 4, 7, 9)

 A: If you analyse the data using an ordinary generalised linear model this is what you get:
> y <- c(group1, group2) 
> group <- factor(rep(1:2, each=25))
> summary(glm(cbind(y, 15 - y) ~ group, binomial))

Call:
glm(formula = cbind(y, 15 - y) ~ group, family = binomial)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.37472  -0.51306  -0.08459   0.80810   1.57631  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.4501     0.1059   4.250 2.14e-05 ***
group2       -0.4661     0.1479  -3.151  0.00163 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 51.401  on 49  degrees of freedom
Residual deviance: 41.403  on 48  degrees of freedom
AIC: 200.93

Number of Fisher Scoring iterations: 3

As you can see, the deviance of 41.4 is smaller than it's expected value of 48 under $H_0$, so if anything, this indicates a slight tendency for these data to be underdispersed.
Such underdispersion would be expected if there is little heterogeneity between test takers within each group (say, no heterogeneity for the sake of the following argument) but heterogeneity instead between questions.  The variance of the score would then be
$$
\operatorname{Var}(\sum_{j=1}^{15} I_j)=\sum_{j=1}^{15}p_j(1-p_j) \le 15\bar p(1-\bar p),
$$
where $\bar p=\frac1{15}\sum_{j=1}^{15}p_j$ (the average probability of correct answer across all questions), that is, smaller than expected based on the binomial model.  This assumes that all test takes are given the same set of questions.
Also note the significant group difference.
