[R] The difference between chisq.test binom.test and pbinom

[Peter Ehlers]

On 2012-08-20 12:24, Tal Galili wrote:
> Hello all,
> I am trying to understand the different results I am getting from the
> following 3 commands:
>
> chisq.test(c(62,50), p = c(0.512,1-0.512), correct = F) # p-value = 0.3788
> binom.test(x=62,n=112, p= 0.512) # p-value = 0.3961
> 2*(1-pbinom(62,112, .512)) # p-value = 0.329
>
> Well, the binom.test was supposed to be "exact" and give the same results
> as the pbinom, while the chisq.test relies on the normal asymptotics.  So I
> would imagine the binom.test should be equal to one of the other two lines,
> but it is not.
>
> The same happens for larger numbers as well:
>
> chisq.test(c(1395,1278), p = c(0.512,1-0.512), correct = F)
> # chisq.test(c(1395,1278), p = c(0.512,1-0.512), simulate.p.value= T)
> binom.test(x=1395,n=2673, p= 0.512)
> 2*(1-pbinom(1395,2673, .512))

Let's first dispense with the chisq.test value; that uses the Normal
approximation and the calculation is trivial; the other two don't.

Let's do a one-sided test:

  binom.test(x=62, n=112, p=.512, alternative="greater")$p.value
  # 0.2161936

Compare with:

  pbinom(61, 112, .512, lower.tail=FALSE)
  # 0.2161936

Why '61' instead of '62'? Because we want to include 62 in the
upper tail probability.

Now it's just a matter of how to calculate the lower tail probability
in binom.test. Hint: it's not just set equal to the upper tail
probability (hence no '2*xxx').
The expected value (under H0) is 112*0.51 = 57.12.
Thus 62 is 4.88 ~ 5 units higher than the EV.
Now calculate the probability of an outcome that is 5 units *less*
(or lower) than the EV:

  pbinom(52, 112, .512)
  # 0.1799121

Add the two probabilities to get:

  0.2161936 + 0.1799121
  # 0.3961057

and that's what binom.test() reports.
For details, have a look at the binom.test code.

Peter Ehlers

>
>
> I'd be happy to know what I might be missing.
>
> Thanks,
> Tal
>
>
> ----------------Contact
> Details:-------------------------------------------------------
> Contact me: Tal.Galili at gmail.com |  972-52-7275845
> Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) |
> www.r-statistics.com (English)
> ----------------------------------------------------------------------------------------------
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>