Dear esteemed colleagues. I’m trying to compare, given binomial data, the p-values obtained from running a Gaussian-family identity-link GLM vs a Binomial-family logit-link fit. I’ve generated 1000 data sets from H0, and counted how many times I reject H0 for an alpha of 0.05. It seems that the Gaussian version rejects H0 in around 5% of the data sets, and the binomial version rejects around 1% of the data sets. This confuses me profusely, shouldn’t an alpha of 0.05 reject 5% of H0 data when using the correct binomial model rather than the incorrect Gaussian model? R-code follows

`base_p = 0.2 # probability in the baseline drug_p = 0.2 # probability in the drug N = 20 # number of observations n = 50 # number of bernoulli experiments p_gauss = rep(NaN, 1000) p_binom = rep(NaN, length(p_gauss)) for (i in 1:length(p_gauss)) { df = data.frame(y = c(rbinom(N/2, n, base_p), rbinom(N/2, n, drug_p)), a = c(rep(0,N/2), rep(1,N/2)), n = n) p_gauss[i] = coef(summary(glm( y ~ a, family = gaussian, data = df)))[,'Pr(>|t|)'][2] p_binom[i] = coef(summary(glm(cbind(y,n) ~ a, family = binomial, data = df)))[,'Pr(>|z|)'][2] } cat(sprintf('gauss reject H0: %.1f%%n', 100 * mean(p_gauss < 0.05))) cat(sprintf('binom reject H0: %.1f%%n', 100 * mean(p_binom < 0.05))) `