distribution of p-values under the null for a binomial test via GLM

binom_p.R

## adapted from Richard McElreath's code at
## https://gist.github.com/rmcelreath/4f7010e8d5688c69bbeb7008f0aabe65

binom_p <- function(N=10, M=5, p=0.25, b=0,
                    test = c("Wald", "LRT")) {
    test <- match.arg(test)
    x <- rep(0:1, each = N/2)
    p <- plogis(qlogis(p)+b*x)
    y <- rbinom(N, size=M, prob=p)
    z <- glm( cbind(y, M-y) ~ x , family=binomial )
    pval <- switch(test,
                   Wald = coef(summary(z))["x", "Pr(>|z|)"],
                   LRT = anova(z)[["Pr(>Chi)"]][2]
                   )
    return(pval)
}

binom_p()
binom_p(test="LRT")

p1 <- replicate(10000, binom_p())
p2 <- replicate(10000, binom_p(test = "LRT"))

plot(density(p1, adjust = 0.5))
par(las = 1)
pfun <- function(x, ...) {
    plot(prop.table(table(x)), axes = FALSE,
         ann = FALSE, ylim = c(0, 0.07), xlim = c(0, 1), lwd = 5,
         lend = 1, yaxs = "i", ...)
}
png("binom_p.png")
pfun(p1)
par(new=TRUE)
tred <- adjustcolor("red", alpha.f = 0.2)
pfun(p2, col = tred)
axis(side = 1)
axis(side = 2)
mtext(side=2, at = 0.035, line = 3, "proportion", las = 3)
mtext(side=1, at = 0.5, line = 3, "p-value")
legend("topleft", lty =1 ,lwd = 5, col = c("black", tred),
       legend = c("Wald", "LRT"))
dev.off()