# How to test whether two binomially distributed variables come from the same distribution?

I have data where I have a set of pairs of integer, each binomially distributed, and I want to be able to test whether the observations came from different distributions, though I don't care what actual distribution each comes from. For example, if n1 = 50, k1 = 20, n2 = 40, k2 = 10, what is the probability of observing the difference between k1 and k2 by chance (that is, if p1 = p2)? I've seen it suggested that one use Fisher's exact test, with the contingency table being:

contingency.tab <- matrix(c(k1[i], size - k1[i], k2[i], size - k2[i]),
nrow=2, ncol=2)


This sort of works for large size parameters, but other wise seems to break down. For example, in the following null case, where p1 = p2 for all observations:

p <- .3
N <- 1000
size1 <- round(runif(N, 50, 150))
size2 <- round(runif(N, 50, 150))
k1 <- rbinom(N, prob = p, size = size1)
k2 <- rbinom(N, prob = p, size = size2)
fisher.pvals_h0 <- sapply(1:N,function(i){fisher.test(matrix(c(k1[i],
size1[i] - k1[i], k2[i], size2[i] - k2[i]), nrow=2, ncol=2), alternative = 'two.sided')$p.value}) hist(fisher.pvals_h0,breaks=40)  There's a serious rightward skew. Here's a case where there is a difference in between p1 and p2: p1 <- .3 p2 <- .36 k3 <- rbinom(N, prob = p1, size = size1) k4 <- rbinom(N, prob = p2, size = size2) fisher.pvals_h1 <- sapply(1:N,function(i){fisher.test(matrix(c(k3[i], size1 - k3[i], k4[i], size2 - k4[i]), nrow=2, ncol=2), alternative = 'two.sided')$p.value})
hist(fisher.pvals_h1,breaks=40)


There are clearly 'significant' differences in some observations, but also a peak near 1.

How would I test for this kind of difference? Thanks.

• If $n_1=50$ and $n_2=40,$ then the distributions definitely differ: there's no question about that and no probability associated with this assertion. Are you trying to ask whether $p_1=p_2$? If so, you're in the textbook setting of a two-sample Binomial test for a difference of means. I wonder, though, because you write of "pairs of integers," which doesn't sound like something that would be described by a Binomial distribution at all. Could you clarify your setting and what you are trying to do? – whuber Oct 15 '19 at 21:48
• Because the test statistic of the Fisher exact test is discrete, you cannot expect the P-value (under the null hypothesis) to be uniform. Your first histogram is roughly uniform. The important part for a test at roughly the 5% level is the probability that the P-value is below 5%. In your second histogram there is more probability below 0.05 because that represents the power of the test (probability of rejection when $H_0$ is false. So it is a good thing that the area in the 2nd histogram has more probability < .05. // BTW, the label of the 2nd histogram should be something like pvals_ha – BruceET Oct 16 '19 at 2:48
• @whuber: the 'pair of integers' refers to k1 and k2. I am indeed asking whether p1 = p2. I tried using prop.test() in R which I think does what your suggesting, binomial two sample test, but it yields a big peak at 1 in the p-value histogram. – user78453 Oct 16 '19 at 21:13
• In your plots the reason for the "peak" at p=1 is the absence of any values between 0.9 and 1, likely because of the discreteness of the sampling statistic. In effect, one bar serves the role of four bars elsewhere in the plot and thereby is four times taller. This makes your plots visually deceptive. To get a more accurate sense of the distribution, consider using either (a) a variable-width histogram or (b) a constrained KDE. – whuber Oct 16 '19 at 21:18

The Fisher Exact Test is based on a hypergeometric distribution, which is discrete. This means that it is not usually possible to do a test at exactly the 5% level. It also means that the distribution of the P-value under the null hypothesis $$H_0: p_1 = p_2$$ does not have a (continuous) uniform distribution on $$(0,1).$$