Repository: scipy/scipy
Pull request: MAINT: stats.fisher_exact: improve docs and fix bugs
#### Reference issue
Closes gh-4130
Closes gh-9231

#### What does this implement/fix?
gh-4130 was originally about inaccurate p-values due to a large fudge factor of 1e-4 used in the implementation. IIUC, the fudge factor is a relative tolerance used to compare numerically distinct but theoeretically equal probability masses. The `hypergeom` distribution is from Boost now, so it is more accurate that what we had before. Perhaps 1e-14 is too tight, but in this case, I think we're looking at small inaccuracies in the reported p-value.

Most of the discussion in that issue is actually about confusion surrounding the definition of two-sided p-value, so I added some additional information to the documentation, including a more intutitive (but also accurate) explanation of the null hypothesis.

gh-9231 exposed a bug that made `fisher_exact` very slow in some cases. This PR fixes it by replacing the binary search in `fisher_exact` with the one written for `binomtest`.


Problem statement:
BUG: stats: fisher_exact returns incorrect p-value
For some tables, `fisher_exact` incorrectly returns a p-value of 1.  For example,

```
In [1]: from scipy.stats import fisher_exact

In [2]: m = np.array([[345, 455], [260, 345]])

In [3]: fisher_exact(m)
Out[3]: (1.0061284868977176, 1.0)
```

The correct p-value for the array `m` is 0.9567.   Here, for example, is the result in R:

```
> m <- array(c(345, 455, 260, 345), c(2,2))
> m
     [,1] [,2]
[1,]  345  260
[2,]  455  345
> fisher.test(m)

    Fisher's Exact Test for Count Data

data:  m
p-value = 0.9567
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.8083048 1.2526776
sample estimates:
odds ratio 
  1.006125 
```

The culprit is this code, from the `fisher_exact` function in `scipy/stats/stats.py`:

```
    elif alternative == 'two-sided':
        mode = int(float((n + 1) * (n1 + 1)) / (n1 + n2 + 2))
        pexact = hypergeom.pmf(c[0,0], n1 + n2, n1, n)
        pmode = hypergeom.pmf(mode, n1 + n2, n1, n)

        epsilon = 1 - 1e-4
        if np.abs(pexact - pmode) / np.maximum(pexact, pmode) <= 1 - epsilon:
            return oddsratio, 1.
```

In effect, it says the p-value is 1 if `np.abs(pexact - pmode) / np.maximum(pexact, pmode) <= 1e-4`.  It isn't clear why the value 1e-4 was used here.


infinite loop in stats.fisher_exact
I'm using scipy version 1.0.0.

```
import scipy as sp
import scipy.stats
x = [[5829225, 5692693], [5760959, 5760959]]
sp.stats.fisher_exact(x)

```
For the values above scipy does not return anything but waits. What can be the reason for that? How can I fix it?

However in R it returns a p-value almost immediately.

```
a = matrix(c(5829225,5692693,5760959,5760959), nrow=2)
fisher.test(a)
```

Thanks in advance, 
Burcak

Hints from issue discussion:
If there is no unit test explaining that special case then could the extra code just be removed?

@argriffing:  Possibly.  I'll do some experimenting.

I poked around a bit more, and found that some sort of "fudge factor" is necessary in the current implementation.  For example, suppose `table = np.array([[5, 0], [1, 4]])` (an example from the unit tests).  The set of all tables with the same marginal sums are shown here, along with their exact probabilities and the probability computed using `stats.hypergeom.pmf`:

```
          Exact      hypergeom.pmf(k, 10, 5, 6)
Table  Probability      (k = table[0,0])
1 4        1/42       0.023809523809523819
5 0

2 3       10/42       0.23809523809523803
4 1

3 2       20/42       0.47619047619047605
3 2

4 1       10/42       0.23809523809523803
2 3

5 0        1/42       0.023809523809523777
1 4
```

For `table = np.array([[5, 0], [1, 4]])`, the exact value of `fisher_exact(table, alternative="two-sided")` should be 2/42.  This is the sum of the probabilities of all the tables whose probability is less than or equal to the probability of `table`.  The scipy code currently uses `stats.hypergeom.pmf` to compute the probabilities.  If a "fudge factor" is not used when searching for the tables to include in the sum, the table `[[1,4],[5,0]]` won't be included in the sum, because its probability (computed with `hypergeom.pmf` and shown above) is greater than the probability of `table`.

I'm not going to spend time tweaking the existing fudge factor.  Even if it is tweaked, other edge cases may exist where the fudge factor results in an incorrect answer.

A way to avoid the need for a fudge factor is to compute the probabilities exactly (using big integers or `fractions.Fraction`).  I have code that does this, but it is _much_ slower than using `stats.hypergeom`.  Perhaps, similar to the functions `comb` and `factorial`, an `exact` argument could be added to `fisher_exact`, and when it is True, the value computed by `fisher_exact` is truly exact.  (Admittedly, it seems a bit funny that a function called `fisher_exact` would have an `exact` argument.)

@WarrenWeckesser  I checked the implementation of this test in R https://github.com/wch/r-source/blob/trunk/src/library/stats/R/fisher.test.R, and seems that they employed a Monte Carlo method. I believe in general the current implementation in scipy won't work. Maybe a rewrite is needed.

Hello. 
I want to report again that this function is inaccurate (https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.fisher_exact.html); and maybe this is related to addition of "alternative"  keyword that doesn't work well (https://github.com/scipy/scipy/issues/4882).

I made some tests 

```
from math import factorial
from scipy.stats import fisher_exact
def hypergeometric(series):
    a,b,c,d=series
p=(factorial(a+b)*factorial(a+c)*factorial(b+d)**factorial(c+d))/(factorial(a)*factorial(b)*factorial(c)*factorial(d)*factorial(a+b+c+d))
    return p
s=np.array([6,108,37,200]).reshape((2,2))
b=hypergeometric(s)

>>b =0.00509270
(this takes so long I was feeling older, don't try this at home)

fisher_exact(s, alternative='two-sided')
>> 0.0050926977481275922)

fisher_exact(s, alternative='less')
>>0.99914550735448127
fisher_exact(s, alternative='two-sided')
>>0.0031603713904877631

```

The following calculator gives
http://www.danielsoper.com/statcalc/calculator.aspx?id=29
Exact hypergeometric probability: 	0.00230588
Exact two-tailed probability: 	0.00509270`

For reference, other calculators give this statistic:
http://www.socscistatistics.com/tests/fisher/Default2.aspx

The value for two-sided test should be double that for the one-sided "greater" one, and this greater one should be equal to 1-"less" p-value, so it seems there is  mistakes here, not only in the doubling byt also in delivered p-values. (For reference, I checked here to be sure: https://stats.stackexchange.com/questions/267192/doubling-or-halving-p-values-for-one-vs-two-tailed-tests)
For the array [[6,108][37, 6828]], values are:
2-sided : 6.2013.... e-05
less: 0.999994886
greater: equal to 2-sided.
So it seems that for some values, it is only a doubling problem. For lower values there is some deeper mess.
The problem is in these lines:
https://github.com/scipy/scipy/blob/v0.14.0/scipy/stats/stats.py#L2450

 ```
if alternative == 'less':
        pvalue = hypergeom.cdf(c[0,0], n1 + n2, n1, n)
    elif alternative == 'greater':
        # Same formula as the 'less' case, but with the second column.
        pvalue = hypergeom.cdf(c[0,1], n1 + n2, n1, c[0,1] + c[1,1])
    elif alternative == 'two-sided':
        mode = int(float((n + 1) * (n1 + 1)) / (n1 + n2 + 2))
        pexact = hypergeom.pmf(c[0,0], n1 + n2, n1, n)
        pmode = hypergeom.pmf(mode, n1 + n2, n1, n)
```
pvalue should be replaced by pexact  for all, and actually the calculation in hypergeom.pmf should not change (still, as per hypergeom docs(https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.hypergeom.html#scipy.stats.hypergeom), n value in this method parameters is N here and n1 is n in pmf. This would be less confusing to choose another name. n is actually confusing anyway, it lets thinking it should be the sum of n1+n2 and this is not the case (it is well visible above, as it is actually the total number of value of the first category, matching or not the condition (categories in columns, conditions in rows).

I could modify the code, but I am not really familiar with github process; while I read a few tutorials I can't wrap my brain around this mess, so I fear I would make a mistake with pull request...
I won't have time to get back to this for a while, but I want to note that this statement

> The value for two-sided test should be double that for the one-sided "greater" one

is not correct in general.  The probability distribution of the tables with fixed marginal sums used in Fisher's exact test is not symmetric, in general.
also, AFAICS, the `pexact = hypergeom.pmf...`  value is the probability for the specific table.
The p-value for the test needs the "cdf", the probability of getting a table at least as "extreme".
more about doubling p-value

there are different ways of defining p-values in the asymmetric case
e.g. for binom_test
https://github.com/statsmodels/statsmodels/issues/2735

fisher exact uses the first definition `minlike`.

(There are some packages in R that have all three version for some "exact" statistics, but I don't remember details right now. The choice of test statistic that defines "at least as extreme" becomes more relevant in larger than 2x2 tables.)

Thanks for the precisions, and sorry for the bad input. Anyway it seemed to me that two-sided value being equal to greater would be illogical/strange. 

regarding

> also, AFAICS, the pexact = hypergeom.pmf... value is the probability for the specific table.
> The p-value for the test needs the "cdf", the probability of getting a table at least as "extreme".

You're right, I went too fast. 
https://github.com/scipy/scipy/blob/v0.14.0/scipy/stats/stats.py#L2450 states that the method returned
p_value : is  the probability of obtaining a distribution **at least as extreme** as the one that was actually observed, assuming that the null hypothesis is true.  
(emphasis' mine)
Hence I guess that the function returns the right value for having at least such imbalanced distribution, using pexact calculated value to perform various checks before using cdf.

Looks like a bug. It gets stuck at the `while` loop here

https://github.com/scipy/scipy/blob/10233db/scipy/stats/stats.py#L3170-L3171

The values compared are 
```
0.0 7.182425487416777e-180
0.0 7.182425487416777e-180
0.0 7.182425487416777e-180
0.0 7.182425487416777e-180
0.0 7.182425487416777e-180
0.0 7.182425487416777e-180
0.0 7.182425487416777e-180
```
Yes, I was definitely suspecting from an infinite loop resulting from comparison.
In this case, how about setting the max_iteration value for the simplicity?
```python
while hypergeom.pmf(guess, n1 + n2, n1, n) < pexact * epsilon:
    # in this point
    if guess > 1000:
        break
    guess += 1
```

For example, if the loop is stopped at the `guess` value of 1000, the results is the same as that of R.

scipy results
```
(1.0239837279122552, 3.0631061987484285e-178)
```

R results
```
        Fisher's Exact Test for Count Data

data:  a
p-value < 2.2e-16
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 1.022285 1.025680
sample estimates:
odds ratio 
  1.023992 
```

Another idea is using `np.isclose()`.
```python
while hypergeom.pmf(guess, n1 + n2, n1, n) < pexact * epsilon:
    if np.isclose(pexact, 0.0):
        break
    guess += 1
```
the problem is `pexact` very close to zero. 
If the value of `pexact` is 0.0, it exits this loop, but it gets too small to get out of the loop.
In this case, inserting the part that checks whether the `pexact` value is close to 0 is solved.
AFAIU:
`guess` is limited by `n`, i.e. the last argument is the number of draws and for we cannot draw more of type I than this.
So, the loop should `break` at  `guess >= n` or add to `while guess < n`
The subtraction path includes the boundary condition `guess > 0`

I'm not sure what should happen at guess == n. AFAICS the final guess needs to be also `n`
@josef-pkt yes, `guess` is limited and controlled by `n`. the problem is `n` is so huge.
In this issue, `n` is `11590184`. this number is too big to iterate a loop.

The loop stopped at `guess=5760958`. 
the running time was roughly 30 minute in my desktop.
the result is `1.02389837279122552`.

last value 
```
(n, guess, hypergeom.pmf(), pexact) = (11590184, 5760958, 7.0148979096678186e-180, 7.183143801796957e-180)
```

In this case, the iteration number seem to be wasted. 
I think we need a separate condition to exit the loop.
If a bisection algorithm is meant by binary search, there is a much faster way of doing this. 
I present a new method. @ilayn, @josef-pkt please review this idea.
`n` value is proportional to the input number, 
in that case `x = [[5829225, 5692693], [5760959, 5760959]]`.
So I applied an adaptive iter_unit.
```python
while hypergeom.pmf(guess, n1 + n2, n1, n) < pexact * epsilon:
    max_order_digit = np.int(np.log10(np.max([n1, n1+n2, n])))
    if max_order_digit <= 3:
        iter_unit = 1
    else:
        iter_unit = np.power(10, max_order_digit-3)
    guess += iter_unit
```
the result is `(1.0239837279122552, 6.126212397496855e-178)`
and running time is 4 sec(current version 30 minutes).

Of course, a more sophisticated process is needed to determine the value of `iter_unit`.
However, applying a variable `iter_unit` can have a tremendous benefit in computation time.
Those tricks usually play well if there is a logarithmic relationship on the size estimation. I don't think it is the case here. 

I am not familiar with `hypergeom.pmf` but if it is a monotone function and  `guess` is bounded by `n` then the obvious choice is implementing smaller binary search inside those branches. Then we can combine `lower` and `upper` and simplify things.
Okay. @ilayn 
I also am not familiar with `hypergeom.pmf`. As far as i know, this function is monotonous.
I will think and design using your suggestion of combining `lower` and `upper` .
How about this process @ilayn ?
previously, regardless of the range of the interval, `guess` is alway increased by 1.
In the new process, i changed the `guess` increase amount not to 1 but to the 
intermediate value of `upper` and `lower`.

I confirmed that applying the new method converges to the results with 21 iteration(in that case).

the previous 30-minute calculation was reduced to about 3 sec.
this is 600x speed up.

```python
while hypergeom.pmf(guess, n1 + n2, n1, n) < pexact * epsilon:
    iter_unit = (np.min([n1, n, n1+n2])//2 + guess) // 2
    guess = iter_unit + 1
```

```bash
iter_index n n1 n2 iter_unit guess hypergeom.pmf pexact
0 11590184 11521918 11521918 2880480 1 0.0 7.183143801796957e-180
1 11590184 11521918 11521918 4320720 2880481 0.0 7.183143801796957e-180
2 11590184 11521918 11521918 5040840 4320721 0.0 7.183143801796957e-180
3 11590184 11521918 11521918 5400900 5040841 0.0 7.183143801796957e-180
4 11590184 11521918 11521918 5580930 5400901 0.0 7.183143801796957e-180
5 11590184 11521918 11521918 5670945 5580931 0.0 7.183143801796957e-180
6 11590184 11521918 11521918 5715952 5670946 0.0 7.183143801796957e-180
7 11590184 11521918 11521918 5738456 5715953 0.0 7.183143801796957e-180
8 11590184 11521918 11521918 5749708 5738457 0.0 7.183143801796957e-180
9 11590184 11521918 11521918 5755334 5749709 9.460732495e-315 7.183143801796957e-180
10 11590184 11521918 11521918 5758147 5755335 1.581289453283023e-242 7.183143801796957e-180
11 11590184 11521918 11521918 5759553 5758148 5.383907224495671e-210 7.183143801796957e-180
12 11590184 11521918 11521918 5760256 5759554 1.2500431856035085e-194 7.183143801796957e-180
13 11590184 11521918 11521918 5760608 5760257 3.5998486997551737e-187 7.183143801796957e-180
14 11590184 11521918 11521918 5760784 5760609 1.7189952961060262e-183 7.183143801796957e-180
15 11590184 11521918 11521918 5760872 5760785 1.1501593175993004e-181 7.183143801796957e-180
16 11590184 11521918 11521918 5760916 5760873 9.332472063069392e-181 7.183143801796957e-180
17 11590184 11521918 11521918 5760938 5760917 2.6530186231704e-180 7.183143801796957e-180
18 11590184 11521918 11521918 5760949 5760939 4.470881558334984e-180 7.183143801796957e-180
19 11590184 11521918 11521918 5760954 5760950 5.803164239342002e-180 7.183143801796957e-180
20 11590184 11521918 11521918 5760957 5760955 6.533410382778032e-180 7.183143801796957e-180
21 11590184 11521918 11521918 5760958 5760958 7.0148979096678186e-180 7.183143801796957e-180
(1.0239837279122552, 6.126212397496855e-178)

real    0m3.077s
user    0m2.296s
sys     0m0.796s

```

Commit message: Update scipy/stats/tests/test_stats.py

Co-authored-by: Pamphile Roy <roy.pamphile@gmail.com>

Patch to review:
diff --git a/scipy/stats/_stats_py.py b/scipy/stats/_stats_py.py
index c8efffaa961d..e14d31a368c8 100644
--- a/scipy/stats/_stats_py.py
+++ b/scipy/stats/_stats_py.py
@@ -4236,6 +4236,15 @@ def pearsonr(x, y, *, alternative='two-sided'):
 def fisher_exact(table, alternative='two-sided'):
     """Perform a Fisher exact test on a 2x2 contingency table.
 
+    The null hypothesis is that the true odds ratio of the populations
+    underlying the observations is one, and the observations were sampled at
+    random from these populations. The statistic returned is the unconditional
+    maximum likelihood estimate of the odds ratio, and the p-value is the
+    probability under the null hypothesis of obtaining a table at least as
+    extreme as the one that was actually observed. There are other possible
+    choices of statistic and two-sided p-value definition associated with
+    Fisher's exact test; please see the Notes for more information.
+
     Parameters
     ----------
     table : array_like of ints
@@ -4244,9 +4253,10 @@ def fisher_exact(table, alternative='two-sided'):
         Defines the alternative hypothesis.
         The following options are available (default is 'two-sided'):
 
-        * 'two-sided'
-        * 'less': one-sided
-        * 'greater': one-sided
+        * 'two-sided': the odds ratio of the underlying population is not one
+        * 'less': the odds ratio of the underlying population is less than one
+        * 'greater': the odds ratio of the underlying population is greater
+          than one
 
         See the Notes for more details.
 
@@ -4255,9 +4265,8 @@ def fisher_exact(table, alternative='two-sided'):
     oddsratio : float
         This is prior odds ratio and not a posterior estimate.
     p_value : float
-        P-value, the probability of obtaining a distribution at least as
-        extreme as the one that was actually observed, assuming that the
-        null hypothesis is true.
+        P-value, the probability under the null hypothesis of obtaining a
+        table at least as extreme as the one that was actually observed.
 
     See Also
     --------
@@ -4273,7 +4282,10 @@ def fisher_exact(table, alternative='two-sided'):
     -----
     *Null hypothesis and p-values*
 
-    The null hypothesis is that the input table is from the hypergeometric
+    The null hypothesis is that the true odds ratio of the populations
+    underlying the observations is one, and the observations were sampled at
+    random from these populations. Equivalently,
+    the null hypothesis is that the input table is from the hypergeometric
     distribution with parameters (as used in `hypergeom`)
     ``M = a + b + c + d``, ``n = a + b`` and ``N = a + c``, where the
     input table is ``[[a, b], [c, d]]``.  This distribution has support
@@ -4461,7 +4473,7 @@ def binary_search(n, n1, n2, side):
         pexact = hypergeom.pmf(c[0, 0], n1 + n2, n1, n)
         pmode = hypergeom.pmf(mode, n1 + n2, n1, n)
 
-        epsilon = 1 - 1e-4
+        epsilon = 1 - 1e-14
         if np.abs(pexact - pmode) / np.maximum(pexact, pmode) <= 1 - epsilon:
             return oddsratio, 1.
 
diff --git a/scipy/stats/tests/test_stats.py b/scipy/stats/tests/test_stats.py
index f51db828c151..1676b63502eb 100644
--- a/scipy/stats/tests/test_stats.py
+++ b/scipy/stats/tests/test_stats.py
@@ -537,6 +537,17 @@ def test_precise(self):
             np.testing.assert_almost_equal(res[1], res_r[1], decimal=11,
                                            verbose=True)
 
+    def test_gh4130(self):
+        # Previously, a fudge factor used to distinguish between theoeretically
+        # and numerically different probability masses was 1e-4; it has been
+        # tightened to fix gh4130. Accuracy checked against R fisher.test.
+        # options(digits=16)
+        # table <- matrix(c(6, 108, 37, 200), nrow = 2)
+        # fisher.test(table, alternative = "t")
+        x = [[6, 37], [108, 200]]
+        res = stats.fisher_exact(x)
+        assert_allclose(res[1], 0.005092697748126)
+
     @pytest.mark.slow
     def test_large_numbers(self):
         # Test with some large numbers. Regression test for #1401


Code review comments:
diff_hunk: @@ -4236,6 +4236,15 @@ def pearsonr(x, y, *, alternative='two-sided'):
 def fisher_exact(table, alternative='two-sided'):
     """Perform a Fisher exact test on a 2x2 contingency table.
 
+    The null hypothesis is that the true odds ratio of the populations
+    underlying the observations is one, and the observations were sampled at
+    random from these populations. The statistic returned is the unconditional
path: scipy/stats/_stats_py.py
text: @author:
Hold off on merge, actually. This is true but incomplete. I need to figure out how to work in something about the marginals of the table being held fixed.