Repository: scikit-learn/scikit-learn
Pull request: FEA add zero_division to matthews_corrcoef
#### Reference Issues/PRs

Partially address: https://github.com/scikit-learn/scikit-learn/issues/29048
Fixes #25258
See also #19977

#### What does this implement/fix? Explain your changes.
The Matthews correlation coefficient is ill-defined (due to zero division) when only one class is present in either the true or predicted labels.

If only one of the true or predicted labels contains a single class, the limit value is 0 (can be shown using polar coordinates). This is sensible as it suggests that the model either provided constant predictions on non-constant data, or variable predictions on single-class data. In such cases, the metric should return a value of 0.

However, if both the true and predicted labels contain only a single class, the limit does not exist, rendering the metric undefined. Also, if both true and model labels are single class, the model succeeded in trivial task and we genuinely can't tell  if the correlation is good or poor. Consequently, in this scenario, the metric should return a nan value. This behavior was chosen to avoid returning 0 for perfect predictions on single-class data.

#### Any other comments?

I noticed that for some other metrics the handling of zero division was based on `_prf_divide()`, but it didn't seem a good fit for this case as its default behavior is to return 0 on both cases and IMHO directly returning nan or 0 yields more readable code and is simpler for users. 

@marctorsoc @glemaitre WDYT?

Problem statement:
Matthews correlation coefficient gives wrong results for single-class result
### Describe the bug

for a single class input, (ie perfect correlation) MCC returns 0, when it should return 1.

since MCC is a special case of pearson correlation, it should fall back to pearson when appropriate

### Steps/Code to Reproduce

from sklearn.metrics import matthews_corrcoef
y_true = [1, 1, 1, 1]
y_pred = [1, 1, 1, 1] 
matthews_corrcoef(y_true, y_pred)
Out  [9]: 0.0

### Expected Results

output should be 1

### Actual Results

output is 0

### Versions

```shell
import sklearn; sklearn.show_versions()

System:
    python: 3.9.12 (main, Apr  5 2022, 01:53:17)  [Clang 12.0.0 ]
   machine: macOS-10.16-x86_64-i386-64bit

Python dependencies:
          pip: 22.3.1
   setuptools: 65.5.1
      sklearn: 1.0.2
        numpy: 1.22.3
        scipy: 1.9.3
       Cython: 0.29.32
       pandas: 1.5.2
   matplotlib: 3.6.2
       joblib: 1.2.0
threadpoolctl: 3.1.0

Built with OpenMP: True
/Users/phild1/opt/anaconda3/lib/python3.9/site-packages/_distutils_hack/__init__.py:33: UserWarning: Setuptools is replacing distutils.
  warnings.warn("Setuptools is replacing distutils.")
```

Hints from issue discussion:
First, I thought that we could quickly fix the problem you pointed out with something like:
```python
if all(y_true == y_pred):
    return 1
```
But this problem is more general: MCC and Pearson Correlation (as far as I understand) are not defined if we have a constant array/variable because the standard deviation would be zero. This would lead to a division by zero, as we can see from the [Wikipedia definition below](https://en.wikipedia.org/wiki/Pearson_correlation_coefficient):
![image](https://user-images.githubusercontent.com/55899543/210184465-897b3915-7510-4e25-99ef-052b67d1798d.png)

This happens in other implementations of the Pearson correlation coefficient as well:
```python
import numpy as np
from scipy.stats import pearsonr

y_true = 10*[1] + [0]
y_pred = 10*[1] + [1]

np.corrcoef(y_true, y_pred)
>>>> RuntimeWarning: invalid value encountered in divide c /= stddev[None, :]
>>>> array([[ 1., nan],
       [nan, nan]])

pearsonr(y_true, y_pred)
>>>> ConstantInputWarning: An input array is constant; the correlation coefficient is not defined.
>>>> PearsonRResult(statistic=nan, pvalue=nan)

# Note that matthews_corrcoef = 0 here for the same reason:
matthews_corrcoef(y_true, y_pred)
>>>> 0
```

But I agree with you that this can be misleading. Especially if you are doing some `sklearn.model_selection.BaseSearchCV` for hyperparameter optimization, you will get some 0 scores (that actually doesn't mean anything). Comparing them with another model with a -0.1 MCC score is not straightforward. Which one is better?
<details>

Indeed MCC is one "available by default" metric in the `sklearn.model_selection.BaseSearchCV` classes.
```python
from sklearn.metrics import get_scorer_names
'matthews_corrcoef' in get_scorer_names()
>>>> True
```

Also, the example I'm talking about is not extreme, as y_pred = 0 for all instances can happen when you have a highly imbalanced problem (and you are using default `.predict`, i.e., with 0.5 as the threshold, which is the case for `sklearn.model_selection.BaseSearchCV`).
</details>

Anyway, I don't see a straightforward change to fix this issue, as using Pearson Correlation doesn't do the trick as far as I can see. Is there another way of using it that you were thinking @dsphilli1?

(Also, you can check #19977 for historical decision/context.)
I think that we are in a corner case. We should probably add a parameter `zero_division` to control the output value and provide a warning in the default case. Regarding if we should output 0.0 or 1.0, we should be looking at the other metric to have something consistent.
@dsphilli1 I added support for this in my related PR: https://github.com/scikit-learn/scikit-learn/pull/23183, as per Guillaume's suggestion
The decision in https://github.com/scikit-learn/scikit-learn/pull/19977 was based on a Wikipedia article that stated:
> If any of the four sums in the denominator is zero, the denominator can be arbitrarily set to one; this results in a Matthews correlation coefficient of zero, which can be shown to be the correct limiting value

While this is correct (and yields favorable metric behavior) when **only one** of the sums is zero, in the case two of them is zero (as in the case in this issue), the (multivariable calculus) limit does not exist and the MCC is 0 for perfect prediction. I've corrected the Wikipedia article accordingly. This might be especially problematic in case someone calculates the metric in batches and aggregates (not the best practice here, but people might do it), then this can change the results even is the data is not really a single class.

I'd suggest returning zero in the cases the limit is zero and a `nan` when it doesn't exist.
Feel free to contribute :)
The same way we introduced the `zero_division` in some metric, we can add the option and being consistent with other metrics.
I think that we did step one of the bullet list given in: https://github.com/scikit-learn/scikit-learn/pull/23183#issuecomment-1411118385

A way forward is to propose small individual PR that makes it easier to review and be merged.
> Feel free to contribute :)

Sure, I'll open a little PR over the weekend

Commit message: Fixed a typo

Patch to review:
diff --git a/sklearn/metrics/_classification.py b/sklearn/metrics/_classification.py
index b6d59d25696ce..a5db319cea4cb 100644
--- a/sklearn/metrics/_classification.py
+++ b/sklearn/metrics/_classification.py
@@ -959,6 +959,24 @@ def matthews_corrcoef(y_true, y_pred, *, sample_weight=None):
         prediction, 0 an average random prediction and -1 and inverse
         prediction).
 
+    Notes
+    -----
+    :func:`matthews_corrcoef` is ill-defined (due to zero division)
+    when only one class is present in either the true or predicted
+    labels.
+
+    If only one of the true or predicted labels contains a single
+    class, the limit value is 0. This is sensible as it suggests
+    that the model either provided constant predictions on
+    non-constant data, or variable predictions on single-class data.
+    In such cases, the metric will return a value of 0.
+
+    However, if both the true and predicted labels contain only a
+    single class, the limit does not exist, rendering the metric
+    undefined. Consequently, in this scenario, the metric will
+    return a nan value. This behavior was chosen to avoid
+    returning 0 for perfect predictions on single-class data.
+
     References
     ----------
     .. [1] :doi:`Baldi, Brunak, Chauvin, Andersen and Nielsen, (2000). Assessing the
@@ -1003,8 +1021,12 @@ def matthews_corrcoef(y_true, y_pred, *, sample_weight=None):
     cov_ypyp = n_samples**2 - np.dot(p_sum, p_sum)
     cov_ytyt = n_samples**2 - np.dot(t_sum, t_sum)
 
-    if cov_ypyp * cov_ytyt == 0:
+    # Zero division, metric limit value is 0
+    if bool(cov_ypyp) ^ bool(cov_ytyt):
         return 0.0
+    # Zero division, metric limit does not exist
+    elif cov_ypyp * cov_ytyt == 0:
+        return np.nan
     else:
         return cov_ytyp / np.sqrt(cov_ytyt * cov_ypyp)
 
diff --git a/sklearn/metrics/tests/test_classification.py b/sklearn/metrics/tests/test_classification.py
index 8f79874e72162..bc7f40eeb3729 100644
--- a/sklearn/metrics/tests/test_classification.py
+++ b/sklearn/metrics/tests/test_classification.py
@@ -774,7 +774,7 @@ def test_cohen_kappa():
 
 
 def test_matthews_corrcoef_nan():
-    assert matthews_corrcoef([0], [1]) == 0.0
+    assert np.isnan(matthews_corrcoef([0], [1]))
     assert matthews_corrcoef([0, 0], [0, 1]) == 0.0
 
 
@@ -886,12 +886,13 @@ def test_matthews_corrcoef():
     y_true_inv2 = np.where(y_true_inv2, "a", "b")
     assert_almost_equal(matthews_corrcoef(y_true, y_true_inv2), -1)
 
-    # For the zero vector case, the corrcoef cannot be calculated and should
-    # output 0
-    assert_almost_equal(matthews_corrcoef([0, 0, 0, 0], [0, 0, 0, 0]), 0.0)
+    # For if both the true and predicted labels contain only a single class,
+    # the metric is undefined and the limit doesn't exist, should return a nan.
+    assert np.isnan(matthews_corrcoef([0, 0, 0, 0], [0, 0, 0, 0]))
 
-    # And also for any other vector with 0 variance
-    assert_almost_equal(matthews_corrcoef(y_true, ["a"] * len(y_true)), 0.0)
+    # If only one of the true or predicted labels contains a single
+    # class is undefined, the limit is 0, should return a 0.
+    assert matthews_corrcoef(y_true, ["a"] * len(y_true)) == 0.0
 
     # These two vectors have 0 correlation and hence mcc should be 0
     y_1 = [1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1]
@@ -952,14 +953,12 @@ def test_matthews_corrcoef_multiclass():
         matthews_corrcoef(y_true, y_pred, sample_weight=sample_weight), -1
     )
 
-    # For the zero vector case, the corrcoef cannot be calculated and should
-    # output 0
+    # For if both the true and predicted labels contain only a single class,
+    # the metric is undefined and the limit doesn't exist, should return a nan.
     y_true = [0, 0, 1, 2]
     y_pred = [0, 0, 1, 2]
     sample_weight = [1, 1, 0, 0]
-    assert_almost_equal(
-        matthews_corrcoef(y_true, y_pred, sample_weight=sample_weight), 0.0
-    )
+    assert np.isnan(matthews_corrcoef(y_true, y_pred, sample_weight=sample_weight))
 
 
 @pytest.mark.parametrize("n_points", [100, 10000])


Code review comments:
diff_hunk: @@ -886,12 +886,13 @@ def test_matthews_corrcoef():
     y_true_inv2 = np.where(y_true_inv2, "a", "b")
     assert_almost_equal(matthews_corrcoef(y_true, y_true_inv2), -1)
 
-    # For the zero vector case, the corrcoef cannot be calculated and should
-    # output 0
-    assert_almost_equal(matthews_corrcoef([0, 0, 0, 0], [0, 0, 0, 0]), 0.0)
+    # For if both the true and predicted labels contain only a single class,
+    # the metric is undefined and the limit doesn't exist, should return a nan.
path: sklearn/metrics/tests/test_classification.py
text: @user1:
```suggestion
    # If both the true and predicted labels contain only a single class,
    # the metric is undefined and the limit doesn't exist, thus should return a nan.
```

diff_hunk: @@ -959,6 +959,24 @@ def matthews_corrcoef(y_true, y_pred, *, sample_weight=None):
         prediction, 0 an average random prediction and -1 and inverse
         prediction).
 
+    Notes
+    -----
+    :func:`matthews_corrcoef` is ill-defined (due to zero division)
+    when only one class is present in either the true or predicted
+    labels.
+
+    If only one of the true or predicted labels contains a single
+    class, the limit value is 0. This is sensible as it suggests
+    that the model either provided constant predictions on
+    non-constant data, or variable predictions on single-class data.
+    In such cases, the metric will return a value of 0.
+
+    However, if both the true and predicted labels contain only a
+    single class, the limit does not exist, rendering the metric
+    undefined. Consequently, in this scenario, the metric will
+    return a nan value. This behavior was chosen to avoid
+    returning 0 for perfect predictions on single-class data.
path: sklearn/metrics/_classification.py
text: @user1:
> and returning a 1 for perfect predictions on single-class data  

But the definition above prays
```
The Matthews correlation coefficient (+1 represents a perfect
        prediction, 0 an average random prediction and -1 and inverse
        prediction).
```

I think we should add a `zero_division` param as proposed in the Issue, so the user decides what's returned for this case. For the case where only one of them has single class, I'm with you on what's sensible. But I could see people arguing otherwise 🤷‍♂️ 

@author:
I agree that the desired behavior for perfect correlation on single class data is arguable. 

Some will say it's perfect prediction and should return 1, some will say it's zero division so should return 0 and raise a warning with an option to choose other return values via `zero_division`, to be in line with the other metrics.

My logic was:
- The two options of returning 0 or 1 are both mathematically (actual value and limit value) not correct.
- Differentiation between the non-existing limit case (perfect correlation for trivial task - we can't tell if the model is good) and the 0 limit case (non uniform predictions for single class data or the opposite - the model is not good) seems important.
-  `scikit-learn` users come from different backgrounds, while returning an arbitrary value and raising a warning might be better for SW engineers, returning nan for an undetermined value seems more like a lingua franca between developers, scientists, statisticians, teachers, students, etc.

The implementation is simple in either case, it's more about style and about how opinionated the module should be.   

@user1:
precisely because not everyone will want the same, I'd argue a `zero_division` param gives flexibility. But let's find some other opinion, esp. from a code owner

@user2:
I also agree that we should have the `zero_division` parameter as we previously decided. @author You can actually refer to the work of @user1 in #23183 where the changes were already done.

@author:
Alright, so the conclusion is that we want the zero limit cases (correspond to bad predictions) and the non-existing limit case (perfect prediction on single class data) to have the same behavior, determined by `zero_division` as in the other metrics?

@user2:
Yes indeed.

@author:
ok 👍

diff_hunk: @@ -886,12 +886,13 @@ def test_matthews_corrcoef():
     y_true_inv2 = np.where(y_true_inv2, "a", "b")
     assert_almost_equal(matthews_corrcoef(y_true, y_true_inv2), -1)
 
-    # For the zero vector case, the corrcoef cannot be calculated and should
-    # output 0
-    assert_almost_equal(matthews_corrcoef([0, 0, 0, 0], [0, 0, 0, 0]), 0.0)
+    # For if both the true and predicted labels contain only a single class,
+    # the metric is undefined and the limit doesn't exist, should return a nan.
+    assert np.isnan(matthews_corrcoef([0, 0, 0, 0], [0, 0, 0, 0]))
 
-    # And also for any other vector with 0 variance
-    assert_almost_equal(matthews_corrcoef(y_true, ["a"] * len(y_true)), 0.0)
+    # If only one of the true or predicted labels contains a single
+    # class is undefined, the limit is 0, should return a 0.
path: sklearn/metrics/tests/test_classification.py
text: @user1:
```suggestion
    # If either the true or predicted labels contains a single
    # class (but not both), the metric is undefined but the limit is 0, thus should return a 0.
```

diff_hunk: @@ -952,14 +953,12 @@ def test_matthews_corrcoef_multiclass():
         matthews_corrcoef(y_true, y_pred, sample_weight=sample_weight), -1
     )
 
-    # For the zero vector case, the corrcoef cannot be calculated and should
-    # output 0
+    # For if both the true and predicted labels contain only a single class,
+    # the metric is undefined and the limit doesn't exist, should return a nan.
path: sklearn/metrics/tests/test_classification.py
text: @user1:
```suggestion
    # If both the true and predicted labels contain only a single class,
    # the metric is undefined and the limit doesn't exist, thus should return a nan.
```