0
$\begingroup$

Problem

I am reading the Wikipedia article on the binomial test and I have been unable to duplicate the purported R and SAS code showing they are equivalent.

The p-value I get for R code is 0.04375, but the p-value I get for SAS code is < .0001. Clearly, these results are not equivalent.

Code

Suppose you suspect your dice is not fair and favors 6. Therefore, you run a trial and flip the dice 235 times. You observe 51 instances of getting a 6. In order to assess whether the dice is fair or not, you run the following tests in SAS and R.

Here is the R code:

binom.test(51, 235, (1/6), alternative = "two.sided")

The SAS is as follows:

data diceRoll;
input success n;
datalines;
0 235
1 51
;

The code for the diceroll is as follows:

proc freq data=diceRoll;
    tables success / binomial (p=0.1666667) ALPHA = 0.05;
    exact binomial; 
    weight n; *number of obs for [0,1];
    run;

The above has been taken directly from the Wikipedia article, but they did not specify the definition of diceroll. Therefore, I have done my best to reconstruct it, but clearly failed.

I'm primarily an R user, but I want to bring up my SAS chops, but clearly missing something obvious here. Please assist.

Reference

$\endgroup$
2
$\begingroup$

To the extent that your question is about the proper method of computing the P-value of a binomial test, I will answer, illustrating results with R code. To the extent that this is a question about debugging SAS code, I believe it's 'off topic' for this site.

Your null hypothesis is that the die is fair so that the probability of getting a 6 is $1/6.$ Then with $n = 235$ rolls of the die, the null distribution of then number of 6's observed is $X \sim \mathsf{Binom}(n=236, p=1/6).$

According to the null hypothesis, you would 'expect' to see $E(X) = np = 236/6 = 39.1667$ 6s. However, you have seen $51$ 6s. If your alternative were $H_a: p > 1/6$ then the P-value would be $$P(X \ge 51|n=236,p=1/6)\\ = 1 - P(X \le 50|n=236,p=1/6) =0.02654.$$

1 - pbinom(50, 235, 1/6)
[1] 0.02654425

If we use binom.test in R to do the same thing, we get exactly the same P-value:

binom.test(51,235,1/6, alt="g")

        Exact binomial test

data:  51 and 235
number of successes = 51, number of trials = 235, 
  p-value = 0.02654
alternative hypothesis: 
  true probability of success is greater than 0.1666667
95 percent confidence interval:
 0.1735253 1.0000000
sample estimates:
probability of success 
             0.2170213 

In the figure below, the P-value is the sum of the heights of the black bars to the right of the vertical red line.

x = 20:70;  pdf = dbinom(x, 235, 1/6)
plot(x, pdf, type="h", lwd=2)
 abline(h = 0, col="green2");  abline(v = 51.5, col="red")

enter image description here

For testing $H_0: p = 1/6$ against the two-sided alternative $H_a: p \ne 1/6,$ the P-value must include area that is as extreme or more extreme than what we observed, in the lower tail. Using the 'center' to be integer 39 (mode), this additional probability is

$$P(X \le 27|n = 235, p = 1/6) = 0.01720.$$

pbinom(27, 235, 1/6)
[1] 0.01720372

Thus the total P-value for the two-sided test is 0.0265 + 0.0172 = 0.0437.

From binom.test we obtain the the same result, consistent with rounding error,

binom.test(51,235,1/6, alt="t")

        Exact binomial test

data:  51 and 235
number of successes = 51, number of trials = 235, 
  p-value = 0.04375
alternative hypothesis: 
  true probability of success is not equal to 0.1666667
95 percent confidence interval:
 0.1660633 0.2752684
sample estimates:
probability of success 
             0.2170213 

enter image description here

So the implementation of the binomial test in R's binom.test seems to perform as it should.

Note: There can be quibbles how to compute the P-value for a two-sided test when the null distribution is not precisely symmetrical. One possibility is to double the P-value from the appropriate one-sided test, another is to use some scheme to interpolate when finding the probability in the opposite tail. Also, some programs approximate the P-value by using a (symmetrical) normal distribution.

None of these alternatives seem to account for the P-value you got in SAS. Although debugging SAS (or other) programs is off-topic here, maybe someone will explain the difference or point out a mistake in your use of the SAS procedure.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.