I’m confused about the operation of binom.test.
Say I want to test a sample of 4/10 success against p=0.5. The P value should be:
P(X <= 4) + P(X >=6) or P(X <= 4) + 1-P(X <= 5)
and indeed:
>pbinom(4,10,p=0.5) + 1-pbinom(5,10,0.5)
[1] 0.7539063
or:
>binom.test(4,10,p=0.5)
Exact binomial test
data: 4 and 10
number of successes = 4, number of trials = 10, p-value = 0.7539
But now I want to test a sample of 95/150 against p=0.66 Here, the expected value is 99, so the P value should be
P(X <= 95) + P(X >= 103) or P(X <= 95) + 1-P(X <= 102)
which is
>pbinom(95,150,.66) + 1-pbinom(102,150,.66)
[1] 0.5464849
but
>binom.test(95,150,.66)
Exact binomial test
data: 95 and 150
number of successes = 95, number of trials = 150, p-value = 0.4914
In fact, the difference in the two P-values is exactly dbinom(103,150,.66). So it seems R has failed to include X=103.
The only explanation I can guess for this is that there is roundoff error due to the inexact representation of .66 causing R to just miss X=103. Is this all it is, or is there something else going on?
Here is the code for computing the p-value in binom.test(x = 95, n = 150, p= 0.66)
relErr <- 1 + 1e-07
d <- dbinom(x, n, p)
m <- n * p
i <- seq.int(from = ceiling(m), to = n)
y <- sum(dbinom(i, n, p) <= d * relErr)
pbinom(x, n, p) + pbinom(n - y, n, p, lower.tail = FALSE)
So, binom.test doesn't look symmetrically about the expected value. It looks for the first integer C such that C is bigger than or equal to the expected value and the probability of getting exactly C successes is less than or equal to the probability of getting exactly x successes, up to the fudge factor in relErr. So, instead of saying that p is the probability of getting "at least that far away from the expected value", they say that p is the probability that the probability is at least as small as the value that you obtained.
In this case,
dbinom(95,n,p)
is 0.05334916. So, binom.test looks for the values of x such that dbinom(x,n,p) is less than 0.05334916. It turns out that those are 0:95 and 104:150. So, binom.test returns the value of
sum(dbinom(0:95,n,p)) + sum(dbinom(104:150,n,p))
which is 0.4914044.