5
$\begingroup$

In a previous question, I asked how to test if a coin is fair. Now I want to empirically test if this test works.

One answer is that programs like R and python have built-in "binomial p-tests" that can be called to do this.

Here is an example of some python code to do such a p-test for a single case of flipping 1000 fair coins:

import numpy as np
from numpy import random
import scipy
from scipy import stats


def flipAndDoPTest(numberOfFlips, weight):
    flippedCoinSimulation = np.random.binomial(1, weight, numberOfFlips) #first input changes sum of coins
    numberOfHeads = np.sum(flippedCoinSimulation==1)
    numberOfTails = np.sum(flippedCoinSimulation==0)
    pvalue = stats.binom_test(numberOfHeads, numberOfFlips, weight)
    return pvalue

numberOfFlips = 1000
weight = .50
ptestvalue = flipAndDoPTest(numberOfFlips, weight)
if ptestvalue>.05:
    print("the ptest has a value of:", ptestvalue)
    print("The null hypothesis cannot be rejected at the 5% level of significance because the returned p-value is greater than the critical value of 5%.")
if ptestvalue<.05:
    print("the ptest has a value of:", ptestvalue)
    print("The null hypothesis can be rejected at the 5% level of significance because the returned p-value is less than the critical value of 5%.")

Now I want to empirically test what this "5% level of significance" means. It seems like there is a lot of disagreement over the interpretation of p-values, so I want to just simulate what occurs in my case.

First I would like to test if 5% of the time a p-value less than .05 will occur for a fair coin. To do this, I repeat this p-test 1000 times (and each p-test is for the event of flipping a fair coin 10000 times). Now I collect all of the times the p-value is less than .05. The code for this is here:

numberOfFlips = 10000
weight = .50
numberOfTests = 1000
StatisticalSignificanceProbability = .05
pTestList = np.zeros(numberOfTests) #initialization
for i in range(numberOfTests):
    #for each i in the loop, do a p-test of 10,000 fair coin flips and add it to a list
    ithPTest = flipAndDoPTest(numberOfFlips, weight)
    pTestList[i] = ithPTest
#take this list and count all of the times there are cases below .05
numberOfSheerCoincidences = sum(pTestList<StatisticalSignificanceProbability)
expectedNumberOfSheerCoincidences = numberOfTests*StatisticalSignificanceProbability

print("numberOfSheerCoincidences: ", numberOfSheerCoincidences)
print("expectedNumberOfSheerCoincidences: ", expectedNumberOfSheerCoincidences)

Now I expect that 5% of my 1000 p-tests are going to be less than .05 (so .05*1000 = 50). But every time I run it I get a number that is significantly smaller than 50. Now this result has a random distribution, so I then wrote code to repeat this procedure to get a histogram distribution of the resultant data:

numberOfFlips = 100
weight = .50   
numberOfDataPoints = 1000
pTestResultsDataPoints = np.zeros(numberOfDataPoints) #initialization
for j in range(numberOfDataPoints):
    #repeating this collection of p-test to get a range of different values
    numberOfTests = 1000
    StatisticalSignificanceProbability = .05
    pTestList = np.zeros(numberOfTests) #initialization
    for i in range(numberOfTests):
        ithPTest = flipAndDoPTest(numberOfFlips, weight)
        pTestList[i] = ithPTest
    numberOfSheerCoincidences = sum(pTestList<StatisticalSignificanceProbability)
    expectedNumberOfSheerCoincidences = numberOfTests*StatisticalSignificanceProbability
    pTestResultsDataPoints[j] = numberOfSheerCoincidences

n, bins, patches = plt.hist(pTestResultsDataPoints, 50)
plt.show()

And with that result I get a distribution that is centered around 35 and not 50. enter image description here

Is this result expected? I was expecting a normal distribution around 50.

$\endgroup$

2 Answers 2

3
$\begingroup$

In general, there is no such binomial test with significance level exactly $\alpha = 0.05,$ on account of the discreteness of binomial distributions.

For an exact test at $\alpha = 0.05$ based on a continuous test statistic, the distribution of P-value when $H_0$ is true would be standard uniform and the probability that the P-value is below $0.05$ would be exactly $0.05.$

If $n = 100,$ testing $H_0: p = .5$ against $H_a: p \ne 0.5,$ the closest one can get to a test at the 5% level (without going over 5%) is $0.0352 = 3.52\%.$

2*(1 - pbinom(60, 100, .5))
[1] 0.0352002
2*(1 - pbinom(59, 100, .5))
[1] 0.05688793

[It does not help to use a normal approximation with nominal 5% level because z values close to $\pm 1.96$ cannot be achieved. For $\mathsf{Binom}(100,.5)$ the normal approximation is very accurate so it hardly matters whether one does an 'exact' binomial test or an approximate normal one.]

Below I simulate 100,000 tests in R for $n=100$ observations, and summarize and plot histograms of P-values. Probabilities should be accurate to about two places. The binomial test 'binom.test' and approximate normal test 'prop.test' both have the anticipated P-values.

set.seed(2021);  n = 100
pv.b = replicate(10^5, 
        binom.test(rbinom(1,n,.5),n,.5)$p.val)
mean(pv.b < 0.05)
[1] 0.03605        # aprx 0.0352
2*sd(pv.b < 0.05)/sqrt(10^5)
[1] 0.001178995    # aprx 95% margin of sim error

set.seed(2021);  n = 100
pv.n = replicate(10^5, 
        prop.test(rbinom(1,n,.5),n,.5)$p.val)
mean(pv.n < 0.05)
[1] 0.03605

The figures below show the simulated distributions of the P-values under $H_0.$

enter image description here

R code for figure:

par(mfrow=c(1,2))
hist(pv.b, prob=T, xlim=c(-.01,1.01), col="skyblue2")
 abline(v = .05, col="red")
 curve(dunif(x), add=T, n=10001, lwd=2)
hist(pv.n, prob=T, xlim=c(-.01,1.01), col="skyblue2")
 abline(v = .05, col="red")
 curve(dunif(x), add=T, n=10001, lwd=2)
par(mfrow=c(1,1))

Notes: (1) At the resolution of the plot, the two histograms look the same, but there are miniscule differences in the distributions; at the borderline the binomial and approximate normal tests had slightly different P-values in 8070 tests out of 100,000. But they never disagreed about rejection at the 5% level.

sum(pv.b == pv.n)
[1] 8070
mean(abs(pv.b - pv.n))
[1] 0.0001792384
sum((pv.n <= .05) == (pv.n <=.05))
[1] 100000

(2) Without using a randomized test, the closest $\alpha$ to 5% without going over for $n=1000$ is $0.046 = 4.6\%,$

2*(1-pbinom(531,1000,.5))
[1] 0.0462912
2*(1-pbinom(530,1000,.5))
[1] 0.05367785
$\endgroup$
4
$\begingroup$

It is (nearly) impossible for a test statistic that has a discrete distribution to attain exactly a 5% (or any other) significance level. So instead, many so-called "exact" tests use the largest attainable significance level that is less than $\alpha$ (or 5% in your case). That is what you are seeing.

Here is a reference that elaborates:

https://www.jstor.org/stable/2684683?origin=crossref

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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