Repository: scipy/scipy
Pull request: Handle constant input for scipy.stats.f_oneway
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
https://docs.scipy.org/doc/numpy/dev/development_workflow.html#writing-the-commit-message

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
<!--Example: Closes gh-WXYZ.-->
Ref. #11669 
This PR does not closes #11669 fully, because there remains question of input which is "almost constant" e.g. up to some small epsilon

#### What does this implement/fix?
<!--Please explain your changes.-->
This PR deals with degenerated cases of input provided to scipy.stats.f_oneway. When lists on input are all of constant values, calculating F and p-value makes a little sense, and some unexpected values are being returned now (before PR). 
#### Additional information
<!--Any additional information you think is important.-->

Problem statement:
p-value varies with the order of the data
<!-- 
Thank you for taking the time to file a bug report. 
Please fill in the fields below, deleting the sections that 
don't apply to your issue. You can view the final output
by clicking the preview button above.
Note: This is a comment, and won't appear in the output.
-->

My issue is about the results of f_oneway function are not consistent. The size of the p-value varies with the order of the data. Is this a bug?

#### Reproducing code example:
<!-- 
If you place your code between the triple backticks below, 
it will be rendered as a code block. 
-->

```
import scipy.stats as stats

a = [2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
b = [1000, 1000, 1000, 1000, 1000]
c = [100, 100, 100, 100, 100, 100]

print(stats.f_oneway(a, c, b))   
print(stats.f_oneway(a, b, c))   
print(stats.f_oneway(b, c, a))   
```

#### Result:
<!-- If any, paste the *full* error message inside a code block
as above (starting from line Traceback)
-->

```
F_onewayResult(statistic=-6.872847162465134e+16, pvalue=nan)
F_onewayResult(statistic=6.872847162465134e+16, pvalue=1.1322762786915646e-143)
F_onewayResult(statistic=-3.436423581232567e+16, pvalue=nan)
```

Hints from issue discussion:
Your data is degenerate, in that the values in each group are constant.  If you work through the analysis of variance calculation, you'll find that this data leads to division by zero when computing the F statistic, so in theory, F should be infinite, and the p-value should be 0.   In the actual numerical calculation, however, the denominator that should be 0 ends up being a small value that is basically numerical noise, and therefore the computed F statistic is, to use the technical term, garbage.  Note that the value of F reported in two of those three results is *negative*.  F should never be negative.

So I guess we can call this a bug, but you should be aware that this is an unusual edge case.
Thank you for your interpretation, we’ll consider this situation!

Thank you at all,
Gangli


开 2020年3月15日 在 下午3:05:52, Warren Weckesser (notifications@github.com) 写到:

Your data is degenerate, in that the values in each group are constant. If you work through the analysis of variance calculation, you'll find that this data leads to division by zero when computing the F statistic, so in theory, F should be infinite, and the p-value should be 0. In the actual numerical calculation, however, the denominator that should be 0 ends up being a small value that is basically numerical noise, and therefore the computed F statistics is, to use the technical term, garbage. Note that the value of F reported in two of those three results is negative. F should never be negative.

So I guess we can call this a bug, but you should be aware that this is an unusual edge case.
@WarrenWeckesser , is there something which can be done to fix this?
Maybe some `UserWarning` about degenerated data? I would be interested in tackling this
Thank you, because the f value cannot be negative, we changed a line in the source code of the f_oneway function, that is, f = msb / msw to f = abs (msb / msw).

Through this we could get the results near of our intention, p-value is near to zero not equal to ‘nan'.

Thank you very much for your enthusiasm and inspiration.

Best regards,
Gangli

开 2020年3月15日 在 下午8:35:33, Vladyslav Rachek (notifications@github.com) 写到:

@WarrenWeckesser , is there something which can be done to fix this?
Maybe some UserWarning about degenerated data? I would be interested in tackling this

—
You are receiving this because you modified the open/close state.
Reply to this email directly, view it on GitHub, or unsubscribe.
> @WarrenWeckesser , is there something which can be done to fix this?
> Maybe some `UserWarning` about degenerated data? I would be interested in tackling this

@erheron 
Not sure what Warren's opinion is here but we have similar behavior for Spearman correlation:
https://github.com/scipy/scipy/blob/c1372d8aa90a73d8a52f135529293ff4edb98fc8/scipy/stats/stats.py#L3223
which is thrown/raised here:
https://github.com/scipy/scipy/blob/c1372d8aa90a73d8a52f135529293ff4edb98fc8/scipy/stats/stats.py#L3373-L3376

so something like that may be in order. I'd be willing to review the PR, unless Warren has a different approach he would prefer in this case (there could be some reasons I'm not thinking why we'd want to do something different with anova than we do with correlations). 

Edit (added): There is also what seems to be a numeric issue here so some investigation might be in order before proceeding directly with raising an error or a warning. 
This bug still exists in latest SciPy :thinking: 
Maybe I don't understand something, but I can't see a linked PR or anything. Shouldn't this issue be re-opened? @Lee-nju  @rlucas7 

We could check for the values within each group being constant.  If they are, and if the common value in at least one group is different from that of another group, the function could return `inf` for the statistic and 0 for the p-value.  If all the values in all the groups are the same, then F is 0/0 and the function should probably return `nan` for F and the p-value.

In `pearsonr`, we generate a warning if an input is "almost constant", because that can result in significant loss of precision.  I haven't thought about whether we need that for `f_oneway`.  I suspect it would be possible to find some nearly constant inputs that generate inaccurate results.



>> On Mar 24, 2020, at 7:03 PM, Warren Weckesser <notifications@github.com> wrote:
> ﻿
> We could check for the values within each group being constant. If they are, and if the common value in at least one group is different from that of another group, the function could return inf for the statistic and 0 for the p-value. If all the values in all the groups are the same, then F is 0/0 and the function should probably return nan for F and the p-value.
> 
This does seem like the most sensible thing to do for the constant case. In Anova terms check that all “within” variances are constant and based on “between” variances either (0, inf) for between variances not all same and 0 and (nan, nan) for between variances all the same and all 0.

> In pearsonr, we generate a warning if an input is "almost constant", because that can result in significant loss of precision. I haven't thought about whether we need that for f_oneway. I suspect it would be possible to find some nearly constant inputs that generate inaccurate results.
> 
Yeah so that is why I was saying that it might require some investigation. Seems likely there are some behaviors like this but not clear to me what would generate them. Handling the constant input case alone would be good for a single PR I’d think. 
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub, or unsubscribe.

Through debugging, I found that the reason for the inconsistency of the floating point number calculation is about e -10. 

> Thank you, because the f value cannot be negative, we changed a line in the source code of the f_oneway function, that is, f = msb / msw to f = abs (msb / msw).

Since we finally got the results we needed, we didn't pay much attention to this issue afterwards.
I'll be now working on the PR, based on the comment of @WarrenWeckesser . Hope to produce it soon

Commit message: Add warning for constant input case

And bunch of minor changes

Patch to review:
diff --git a/scipy/stats/stats.py b/scipy/stats/stats.py
index 80615f5e6b97..ee364697f959 100644
--- a/scipy/stats/stats.py
+++ b/scipy/stats/stats.py
@@ -192,7 +192,8 @@
            'scoreatpercentile', 'percentileofscore',
            'cumfreq', 'relfreq', 'obrientransform',
            'sem', 'zmap', 'zscore', 'iqr', 'gstd', 'median_absolute_deviation',
-           'sigmaclip', 'trimboth', 'trim1', 'trim_mean', 'f_oneway',
+           'sigmaclip', 'trimboth', 'trim1', 'trim_mean',
+           'f_oneway', 'F_onewayConstantInputWarning',
            'PearsonRConstantInputWarning', 'PearsonRNearConstantInputWarning',
            'pearsonr', 'fisher_exact', 'SpearmanRConstantInputWarning',
            'spearmanr', 'pointbiserialr',
@@ -3116,6 +3117,17 @@ def trim_mean(a, proportiontocut, axis=0):
 F_onewayResult = namedtuple('F_onewayResult', ('statistic', 'pvalue'))
 
 
+class F_onewayConstantInputWarning(RuntimeWarning):
+    """Warning generated by `f_oneway` when an input is constant, e.g.
+       each of the samples provided is a constant array"""
+
+    def __init__(self, msg=None):
+        if msg is None:
+            msg = ("Each of the input arrays is constant;"
+                    "the F-value and p-value are not defined or infinite")
+        self.args = (msg,)
+
+
 def f_oneway(*args):
     """
     Perform one-way ANOVA.
@@ -3136,6 +3148,13 @@ def f_oneway(*args):
     pvalue : float
         The associated p-value from the F-distribution.
 
+    Warns 
+    -------
+    F_onewayConstantInputWarning
+        Raised if each of the input arrays is constant array.
+        In this case F-value should be either infinite or isn't defined, so
+        ``np.inf`` or ``np.nan`` is returned for F-value
+
     Notes
     -----
     The ANOVA test has important assumptions that must be satisfied in order
@@ -3150,8 +3169,14 @@ def f_oneway(*args):
     possible to use the Kruskal-Wallis H-test (`scipy.stats.kruskal`) although
     with some loss of power.
 
+    If each group is made of constant values, and
+    
+    - There exists 2 groups with different values, the function returns (``np.inf``, 0)
+    - All values in all groups are the same, function returns (``np.nan``, ``np.nan``)
+
     The algorithm is from Heiman[2], pp.394-7.
 
+
     References
     ----------
     .. [1] R. Lowry, "Concepts and Applications of Inferential Statistics",
@@ -3192,6 +3217,22 @@ def f_oneway(*args):
     alldata = np.concatenate(args)
     bign = len(alldata)
 
+    # Check if the values within each group are constant
+    # And the common value in at least one group is different
+    # from that in another group - special cases
+    # Based on https://github.com/scipy/scipy/issues/11669
+    const_groups = True
+    for group in args:
+        if not all(x == group[0] for x in group):
+            const_groups = False
+            break
+    if const_groups:
+        warnings.warn(F_onewayConstantInputWarning())
+        if len(set(group[0] for group in args)) > 1:
+            return F_onewayResult(np.inf, 0)
+        else:
+            return F_onewayResult(np.nan, np.nan)
+
     # Determine the mean of the data, and subtract that from all inputs to a
     # variance (via sum_of_sq / sq_of_sum) calculation.  Variance is invariance
     # to a shift in location, and centering all data around zero vastly
diff --git a/scipy/stats/tests/test_stats.py b/scipy/stats/tests/test_stats.py
index c739d9de02ea..98d2253c5a4a 100644
--- a/scipy/stats/tests/test_stats.py
+++ b/scipy/stats/tests/test_stats.py
@@ -4450,6 +4450,16 @@ def test_nist(self):
             assert_allclose(res[0], f, rtol=rtol,
                             err_msg='Failing testcase: %s' % test_case)
 
+    def test_constant_input(self):
+        # For more details, look on https://github.com/scipy/scipy/issues/11669
+        a = np.array([42, 42, 42])
+        b = np.array([7, 7, 7, 7])
+        with assert_warns(stats.F_onewayConstantInputWarning):
+            res = stats.f_oneway(a, b)
+            assert_equal(res, (np.inf, 0))
+            res = stats.f_oneway(a,a)
+            assert_equal(res, (np.nan, np.nan))
+
 
 class TestKruskal(object):
     def test_simple(self):


Code review comments:
diff_hunk: @@ -3116,6 +3117,17 @@ def trim_mean(a, proportiontocut, axis=0):
 F_onewayResult = namedtuple('F_onewayResult', ('statistic', 'pvalue'))
 
 
+class F_onewayConstantInputWarning(RuntimeWarning):
+    """Warning generated by `f_oneway` when an input is constant, e.g.
+       each of the samples provided is a constant array"""
+
+    def __init__(self, msg=None):
+        if msg is None:
+            msg = ("Each of the input arrays is constant;"
+                    "the F-value and p-value are not defined or infinite")
path: scipy/stats/stats.py
text: @user1:
suggest to change F-value to F-statistic

@author:
Shouldn't we also change to ``F-statistic`` here? "F-value" is in documentation now, and I think the change must be kind of consistent
https://github.com/scipy/scipy/blob/b818978f3580e05a3885006cae337f28bd2aaf7d/scipy/stats/stats.py#L3137-L3142

UPD: after quick Google search, that page suggests that "F-value" is a good choice https://www.statisticshowto.datasciencecentral.com/probability-and-statistics/f-statistic-value-test/
They say that "F-statistic" denotes something different

@user1:
Ok, let's stick with what has been in the doc rather than change.

diff_hunk: @@ -3136,6 +3148,13 @@ def f_oneway(*args):
     pvalue : float
         The associated p-value from the F-distribution.
 
+    Warns 
+    -------
+    F_onewayConstantInputWarning
+        Raised if each of the input arrays is constant array.
+        In this case F-value should be either infinite or isn't defined, so
path: scipy/stats/stats.py
text: @user1:
here maybe state: 
    In this case the F-statistic is either infinite or not defined, so ...

@author:
Did, but left "F-value"

diff_hunk: @@ -3152,8 +3171,8 @@ def f_oneway(*args):
 
     If each group is made of constant values, and
     
-    - Exists 2 groups with different values, function returns (np.inf, 0)
-    - All values in all groups are the same, function returns (np.nan, np.nan)
+    - There exists 2 groups with different values, the function returns (``np.inf``, 0)
path: scipy/stats/stats.py
text: @user1:
sorry to nitpick this but change: 
exists 2 groups => exist at least 2 groups
 

@author:
Yeeees... but in "mathematical" sense these two sentences are equal, and the former was shorter that's why I went for it. Will do :)