Using `glm()` in R for correlated count data (data provided) I have data from one group ($n = 25$) of subjects pre- and post- tested on an $15$-item dichotomously scored test (see below).
Given that the data are from the same subjects (i.e., correlated), how can I correctly use count regression to compare the performance of the group at the two testing points (i.e., pre and post)?
I appreciate an R demonstration.
Here is what I have done (ignoring the correlated data) in R:
pre = c(9,10,8,9,6,7,9,7,8,6,8,4,8,11,6,7,9,9,6,6,9,9,10,9,8)

post = c(13,10,10,12,9,11,12,9,7,9,7,3,7,9,6,10,10,12,8,9,8,9,10,10,7)


y <- c(pre, post) 
group <- rep(0:1, each = 25)
summary(model <- glm(cbind(y, 15 - y) ~ group, binomial) )

 A: You should use the Wilcoxon signed-rank test, which is essentially a non-parametric version of the paired t-test. In R, you would run
wilcox.test(pre, post, paired=TRUE)

You could also run a paired t-test (ignoring the count nature of the data) using t.test(), which will yield the same conclusion (although in this case, the non-normality of the difference scores and the small sample size might steer you away from this method).
A: This isn't what's typically meant by "count data" in statistics.  When statisticians talk about count data, they are referring to non-negative whole numbers that have no upper limit ($[0, \infty)$).  Your data are bounded counts (counts out of a finite total).  They are binomial data.  These would typically be modeled with logistic regression (or probit, etc.), not "count regression", which typically implies Poisson or negative binomial regression.
The next issue is that you have repeated measures, which means the data are not independent.  The proposed model treats the data as independent—that will lead to incorrect standard errors, p-values, etc.
There are several ways to model data like these.  The most common would probably be to use a GLMM to model scores as a function of time:
library(lme4)                # we'll need this package
y     = c(pre, post) 
time  = rep(0:1, each = 25)
id    = rep(1:25, times=2)   # we need an indicator for the participants
model = glmer(cbind(y, 15-y)~time+(1|id), family=binomial)  # this fits the GLMM
summary(model, correlation=FALSE)
# Generalized linear mixed model fit by maximum likelihood (Laplace
#   Approximation) [glmerMod]
#  Family: binomial  ( logit )
# Formula: cbind(y, 15 - y) ~ time + (1 | id)
# 
#      AIC      BIC   logLik deviance df.resid 
#    209.4    215.2   -101.7    203.4       47 
# 
# Scaled residuals: 
#     Min      1Q  Median      3Q     Max 
# -2.2694 -0.5616  0.1626  0.3804  1.6842 
# 
# Random effects:
#  Groups Name        Variance Std.Dev.
#  id     (Intercept) 0.07019  0.2649  
# Number of obs: 50, groups:  id, 25
# 
# Fixed effects:
#             Estimate Std. Error z value Pr(>|z|)  
# (Intercept)   0.1138     0.1170   0.973   0.3307  
# time          0.3211     0.1492   2.152   0.0314 *
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

