Repository: scikit-learn/scikit-learn
Pull request: DOC note on calibration impact on ranking
<!--
Thanks for contributing a pull request! Please ensure you have taken a look at
the contribution guidelines: https://github.com/scikit-learn/scikit-learn/blob/main/CONTRIBUTING.md
-->

#### Reference Issues/PRs
<!--
Example: Fixes #1234. See also #3456.
Please use keywords (e.g., Fixes) to create link to the issues or pull requests
you resolved, so that they will automatically be closed when your pull request
is merged. See https://github.com/blog/1506-closing-issues-via-pull-requests
-->
Fixes #16321.


#### What does this implement/fix? Explain your changes.


#### Any other comments?


<!--
Please be aware that we are a loose team of volunteers so patience is
necessary; assistance handling other issues is very welcome. We value
all user contributions, no matter how minor they are. If we are slow to
review, either the pull request needs some benchmarking, tinkering,
convincing, etc. or more likely the reviewers are simply busy. In either
case, we ask for your understanding during the review process.
For more information, see our FAQ on this topic:
http://scikit-learn.org/dev/faq.html#why-is-my-pull-request-not-getting-any-attention.

Thanks for contributing!
-->


Problem statement:
Isotonic calibration changes rank-based test metrics values
[As discussed with @ogrisel]

#### Describe the bug
Using isotonic calibration changes metrics values. This is because it is a non-strictly monotonic calibration. Sigmoid calibration being strictly monotonic doesn't suffer from this.

#### Steps/Code to Reproduce
Here is quick example where we split in three train/calibration/test sets and compare ROC AUC on the test set before and after calibrating for both isotonic and sigmoid..

```python
import numpy as np
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.calibration import CalibratedClassifierCV
from sklearn.metrics import roc_auc_score

X, y = datasets.make_classification(n_samples=100000, n_features=20,
                                    n_informative=18, n_redundant=0,
                                    random_state=42)

X_train, X_, y_train, y_ = train_test_split(X, y, test_size=0.5,
                                                    random_state=42)
X_test, X_calib, y_test, y_calib = train_test_split(X_, y_, test_size = 0.5,
                                                    random_state = 42)

clf = LogisticRegression(C=1.)
clf.fit(X_train, y_train)

y_pred = clf.predict_proba(X_test)
print(roc_auc_score(y_test, y_pred[:,1]))
```
The ROC AUC is then:    0.88368

```python
isotonic = CalibratedClassifierCV(clf, method='isotonic', cv='prefit')
isotonic.fit(X_calib, y_calib)

y_pred_calib = isotonic.predict_proba(X_test)
print(roc_auc_score(y_test, y_pred_calib[:,1]))
```
After isotonic calibration, the ROC AUC becomes:    0.88338

```python
isotonic = CalibratedClassifierCV(clf, method='sigmoid', cv='prefit')
isotonic.fit(X_calib, y_calib)

y_pred_calib = isotonic.predict_proba(X_test)
print(roc_auc_score(y_test, y_pred_calib[:,1]))
```
As expected for sigmoid calibration, the ROC AUC is constant.    0.88368

#### Versions
System:
    python: 3.8.1 | packaged by conda-forge | (default, Jan 29 2020, 15:06:10)  [Clang 9.0.1 ]
executable: /Users/leodreyfusschmidt/opt/miniconda2/envs/isotonic/bin/python
   machine: macOS-10.13.4-x86_64-i386-64bit

Python dependencies:
       pip: 20.0.2
setuptools: 45.1.0.post20200119
   sklearn: 0.22.1
     numpy: 1.16.5
     scipy: 1.4.1
    Cython: None
    pandas: None
matplotlib: 3.1.2
    joblib: 0.14.1

Built with OpenMP: True

Hints from issue discussion:
Thank you very much @dsleo for the detailed report. I started from your code and pushed the analysis further.

If one introspect the calibrator, on can observe the following:

```python

y_pred_calib = isotonic.predict_proba(X_test)
print(roc_auc_score(y_test, y_pred_calib[:, 1]))

calibrator = isotonic.calibrated_classifiers_[0].calibrators_[0]

import matplotlib.pyplot as plt
plt.figure(figsize=(16, 10))
plt.plot(calibrator._necessary_X_, calibrator._necessary_y_)
```
![fig1](https://user-images.githubusercontent.com/89061/73534572-b9cc3980-4421-11ea-9070-e71fadaf9357.png)

Note that this kind of plot is very interesting and there is another related PR (#16289) to make it possible to expose those thresholds for inspection reasons here.

This is a bit weird to have a constant step wise function to map `y_raw_pred` to `y_cal_pred`. 

Furthermore a close look at the (x, y) pairs of thresholds look as follows (for the first 10 pairs):
```python
for x, y in zip(calibrator._necessary_X_[:10],
                calibrator._necessary_y_[:10]):
    print(x, y)
```
```
-8.31123478939124 0.0
-5.600958281618134 0.0
-5.590794325886634 0.00980392156862745
-5.044899575456449 0.00980392156862745
-5.042974497498811 0.017804154302670624
-4.303066186168206 0.017804154302670624
-4.301859876807849 0.024193548387096774
-4.151532225028019 0.024193548387096774
-4.149327735774874 0.027548209366391185
-3.763637782629291 0.027548209366391185
```

The y values are monotonic (but not strictly) <del>but what is weird is that the x values are not even monotonic and go a bit back in time at each step</del>. This is a bit fishy but fair enough, once we build the `scipy.interpolation.interp1d` function on the threshold in `IsotonicRegression` these artifacts can be ignored.


I don't know why we want to use such a piece-wise constant mapping. One could instead use a piece-wise linear mapping that would be trivially strictly monotonic. Here is a quick hack to show that this would trivially fix the issue when using isotonic calibration for classifier calibration:

```python
X_thresholds_fixed = np.concatenate((calibrator._necessary_X_[::2],
                                     calibrator._necessary_X_[-1:]))
y_thresholds_fixed = np.concatenate((calibrator._necessary_y_[::2],
                                     calibrator._necessary_y_[-1:]))


plt.figure(figsize=(16, 10))
plt.plot(calibrator._necessary_X_, calibrator._necessary_y_)
plt.plot(X_thresholds_fixed, y_thresholds_fixed)
```
![fig2](https://user-images.githubusercontent.com/89061/73535600-39f39e80-4424-11ea-88a9-2ef052298209.png)

Linearly interpolating using those thresholds makes the mapping strictly monotonic and the ROC-AUC of the original model is recovered:

```python

calibrator._build_f(X_thresholds_fixed, y_thresholds_fixed)
y_pred_calib = isotonic.predict_proba(X_test)

print(roc_auc_score(y_test, y_pred[:, 1]))
print(roc_auc_score(y_test, y_pred_calib[:, 1]))
```
```
0.8836876399954915
0.8836876399954915
```

Instead of the max values of each steps, one could have used mid-points. This could be explored in a fix in `IsotonicRegression` itself (maybe as a new option to keep the piece-wise constant mapping as the default behavior for backward compat reasons).
The above analysis is wrong: I made wrong expectations on the structure of the steps. They are not always exact steps. There are already piece-wise linear components in the default prediction function.

Still it would be nice to have an option to make it possible to enforce a strict monotonicity option, maybe by adding a small eps on one of the edges whenever y is constant on a segment.
Metrics that should not be impacted by Isotonic Calibration but are:

- ROC AUC
- Average Precision
- NDCG...
> Linearly interpolating using those thresholds makes the mapping strictly monotonic and the ROC-AUC of the original model is recovered:

This would solve the problem but we would have to make clear that we are not strictly doing isotonic regression.

I had a look at what R does. While the core stats package and a number of external packages implement isotonic regression but none seem to deal with this issue.

There is also few mentions of this in literature (that I could fine). ['Predicting accurate probabilities with a ranking loss'](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4180410/) mentions:

> Observe that isotonic regression preserves the ordering of the input model’s scores, although potentially introducing ties i.e. f(ŝ(·)) is not injective. To break ties on training examples, we may simply refer to the corresponding original model’s scores.

This would be a simple way to ensure rank stays consistent and we would still strictly be using isotonic regression but the API would be difficult to work out.

The only other paper I found that touches on this ['Smooth Isotonic Regression: A New Method to Calibrate Predictive Models'](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3248752/). Their method would solve the rank problem but I think it is too complex for our purposes. Summary of their method:

> 1. Obtain f* = argminf∑i (ci− f(pi))2, subject to f(pi) ≤ f(pi+1), ∀i (Isotonic Regression).
> 2. Sample s points from f(γ), γ ∈ (0, 1), one point per flat region of f(·). Denote samples as P′, and their corresponding class labels as C′.
> 3. Construct a Piecewise Cubic Hermite Interpolating Polynomial function t(f(P′), C′)) to obtain a monotone smoothing spline as the final transformation function for calibration.
Should we just add a note in the UG, something like

> It is generally expected that calibration does not affect ranking metrics such as ROC-AUC or Average Prevision. However, these metrics might differ after calibration when using method='isotonic' since isotonic regression introduces ties in the predicted probabilities. Using method='logistic' will not suffer from this because the sigmoid is strictly monotonic.


Slightly related, it would be interesting to discuss in the UG when one might prefer sigmoid vs isotonic. I personally have no idea.
> Slightly related, it would be interesting to discuss in the UG when one might prefer sigmoid vs isotonic. I personally have no idea.

I am happy to add this to the doc. This was explained in some of the papers I looked at.


please ping me on the PR :)
> This would be a simple way to ensure rank stays consistent and we would still strictly be using isotonic regression but the API would be difficult to work out.

I am curious if this can be worked into our `IsotonicRegression`. We can have it as an option if it adds too much computational cost.
> it would be interesting to discuss in the UG when one might prefer sigmoid vs isotonic. I personally have no idea.

sigmoid is parametric and strongly biased (assume the calibration curve of the un-calibrated model has a sigmoid shape) while isotonic calibration is non-parametric and does not make this assumption.

In general, the non-parametric calibration would tend to overfit on dataset with a small number of samples while the parametric sigmoid calibration can underfit if the assumption is wrong.

TL;DR: sigmpoid on small datasets, isotonic on large datasets.
> I had a look at what R does. While the core stats package and a number of external packages implement isotonic regression but none seem to deal with this issue.

Not all application of isotonic regression are for probabilistic calibration of classifiers. For the calibration use case it might make sense to further impose strict monotonicity.

This could be a non-default option for the `IsotonicRegression` class it-self but that would be enabled by default in `CalibratedClassifierCV` when using the isotonic method.

Commit message: DOC fix doi

Patch to review:
diff --git a/doc/developers/contributing.rst b/doc/developers/contributing.rst
index 9a3efd3d2c8db..a0680ed87887e 100644
--- a/doc/developers/contributing.rst
+++ b/doc/developers/contributing.rst
@@ -714,7 +714,7 @@ Building the documentation requires installing some additional packages:
 
     pip install sphinx sphinx-gallery numpydoc matplotlib Pillow pandas \
                 scikit-image packaging seaborn sphinx-prompt \
-                sphinxext-opengraph plotly
+                sphinxext-opengraph plotly pooch
 
 To build the documentation, you need to be in the ``doc`` folder:
 
@@ -736,7 +736,7 @@ To also generate the example gallery you can use:
 
 .. prompt:: bash $
 
-    make html
+    make doc
 
 This will run all the examples, which takes a while. If you only want to
 generate a few examples, you can use:
diff --git a/doc/modules/calibration.rst b/doc/modules/calibration.rst
index e35a4a12cbb8a..9828c5bf77708 100644
--- a/doc/modules/calibration.rst
+++ b/doc/modules/calibration.rst
@@ -20,10 +20,24 @@ prediction.
 Well calibrated classifiers are probabilistic classifiers for which the output
 of the :term:`predict_proba` method can be directly interpreted as a confidence
 level.
-For instance, a well calibrated (binary) classifier should classify the samples
-such that among the samples to which it gave a :term:`predict_proba` value
-close to 0.8,
-approximately 80% actually belong to the positive class.
+For instance, a well calibrated (binary) classifier should classify the samples such
+that among the samples to which it gave a :term:`predict_proba` value close to, say,
+0.8, approximately 80% actually belong to the positive class.
+
+Before we show how to re-calibrate a classifier, we first need a way to detect how
+good a classifier is calibrated.
+
+.. note::
+    Strictly proper scoring rules for probabilistic predictions like
+    :func:`sklearn.metrics.brier_score_loss` and
+    :func:`sklearn.metrics.log_loss` assess calibration (reliability) and
+    discriminative power (resolution) of a model, as well as the randomness of the data
+    (uncertainty) at the same time. This follows from the well-known Brier score
+    decomposition of Murphy [1]_. As it is not clear which term dominates, the score is
+    of limited use for assessing calibration alone (unless one computes each term of
+    the decomposition). A lower Brier loss, for instance, does not necessarily
+    mean a better calibrated model, it could also mean a worse calibrated model with much
+    more discriminatory power, e.g. using many more features.
 
 .. _calibration_curve:
 
@@ -33,7 +47,7 @@ Calibration curves
 Calibration curves, also referred to as *reliability diagrams* (Wilks 1995 [2]_),
 compare how well the probabilistic predictions of a binary classifier are calibrated.
 It plots the frequency of the positive label (to be more precise, an estimation of the
-*conditional event probability* :math:`P(Y=1|\text{predict\_proba})`) on the y-axis
+*conditional event probability* :math:`P(Y=1|\text{predict_proba})`) on the y-axis
 against the predicted probability :term:`predict_proba` of a model on the x-axis.
 The tricky part is to get values for the y-axis.
 In scikit-learn, this is accomplished by binning the predictions such that the x-axis
@@ -62,7 +76,7 @@ by showing the number of samples in each predicted probability bin.
 
 :class:`LogisticRegression` returns well calibrated predictions by default as it has a
 canonical link function for its loss, i.e. the logit-link for the :ref:`log_loss`.
-This leads to the so-called **balance property**, see [7]_ and
+This leads to the so-called **balance property**, see [8]_ and
 :ref:`Logistic_regression`.
 In contrast to that, the other shown models return biased probabilities; with
 different biases per model.
@@ -79,7 +93,7 @@ case in this dataset which contains 2 redundant features.
 :class:`RandomForestClassifier` shows the opposite behavior: the histograms
 show peaks at probabilities approximately 0.2 and 0.9, while probabilities
 close to 0 or 1 are very rare. An explanation for this is given by
-Niculescu-Mizil and Caruana [1]_: "Methods such as bagging and random
+Niculescu-Mizil and Caruana [3]_: "Methods such as bagging and random
 forests that average predictions from a base set of models can have
 difficulty making predictions near 0 and 1 because variance in the
 underlying base models will bias predictions that should be near zero or one
@@ -99,7 +113,7 @@ to 0 or 1 typically.
 .. currentmodule:: sklearn.svm
 
 :class:`LinearSVC` (SVC) shows an even more sigmoid curve than the random forest, which
-is typical for maximum-margin methods (compare Niculescu-Mizil and Caruana [1]_), which
+is typical for maximum-margin methods (compare Niculescu-Mizil and Caruana [3]_), which
 focus on difficult to classify samples that are close to the decision boundary (the
 support vectors).
 
@@ -167,29 +181,18 @@ fit the regressor. It is up to the user to
 make sure that the data used for fitting the classifier is disjoint from the
 data used for fitting the regressor.
 
-:func:`sklearn.metrics.brier_score_loss` may be used to assess how
-well a classifier is calibrated. However, this metric should be used with care
-because a lower Brier score does not always mean a better calibrated model.
-This is because the Brier score metric is a combination of calibration loss
-and refinement loss. Calibration loss is defined as the mean squared deviation
-from empirical probabilities derived from the slope of ROC segments.
-Refinement loss can be defined as the expected optimal loss as measured by the
-area under the optimal cost curve. As refinement loss can change
-independently from calibration loss, a lower Brier score does not necessarily
-mean a better calibrated model.
-
-:class:`CalibratedClassifierCV` supports the use of two 'calibration'
-regressors: 'sigmoid' and 'isotonic'.
+:class:`CalibratedClassifierCV` supports the use of two regression techniques
+for calibration via the `method` parameter: `"sigmoid"` and `"isotonic"`.
 
 .. _sigmoid_regressor:
 
 Sigmoid
 ^^^^^^^
 
-The sigmoid regressor is based on Platt's logistic model [3]_:
+The sigmoid regressor, `method="sigmoid"`` is based on Platt's logistic model [4]_:
 
 .. math::
-       p(y_i = 1 | f_i) = \frac{1}{1 + \exp(A f_i + B)}
+       p(y_i = 1 | f_i) = \frac{1}{1 + \exp(A f_i + B)} \,,
 
 where :math:`y_i` is the true label of sample :math:`i` and :math:`f_i`
 is the output of the un-calibrated classifier for sample :math:`i`. :math:`A`
@@ -200,10 +203,10 @@ The sigmoid method assumes the :ref:`calibration curve <calibration_curve>`
 can be corrected by applying a sigmoid function to the raw predictions. This
 assumption has been empirically justified in the case of :ref:`svm` with
 common kernel functions on various benchmark datasets in section 2.1 of Platt
-1999 [3]_ but does not necessarily hold in general. Additionally, the
+1999 [4]_ but does not necessarily hold in general. Additionally, the
 logistic model works best if the calibration error is symmetrical, meaning
 the classifier output for each binary class is normally distributed with
-the same variance [6]_. This can be a problem for highly imbalanced
+the same variance [7]_. This can be a problem for highly imbalanced
 classification problems, where outputs do not have equal variance.
 
 In general this method is most effective for small sample sizes or when the
@@ -213,7 +216,7 @@ high and low outputs.
 Isotonic
 ^^^^^^^^
 
-The 'isotonic' method fits a non-parametric isotonic regressor, which outputs
+The `method="isotonic"` fits a non-parametric isotonic regressor, which outputs
 a step-wise non-decreasing function, see :mod:`sklearn.isotonic`. It minimizes:
 
 .. math::
@@ -226,10 +229,20 @@ calibrated classifier for sample :math:`i` (i.e., the calibrated probability).
 This method is more general when compared to 'sigmoid' as the only restriction
 is that the mapping function is monotonically increasing. It is thus more
 powerful as it can correct any monotonic distortion of the un-calibrated model.
-However, it is more prone to overfitting, especially on small datasets [5]_.
+However, it is more prone to overfitting, especially on small datasets [6]_.
 
 Overall, 'isotonic' will perform as well as or better than 'sigmoid' when
-there is enough data (greater than ~ 1000 samples) to avoid overfitting [1]_.
+there is enough data (greater than ~ 1000 samples) to avoid overfitting [3]_.
+
+.. note:: Impact on ranking metrics like AUC
+
+    It is generally expected that calibration does not affect ranking metrics such as
+    ROC-AUC. However, these metrics might differ after calibration when using
+    `method="isotonic"` since isotonic regression introduces ties in the predicted
+    probabilities. This can be seen as within the uncertainty of the model predictions.
+    In case, you strictly want to keep the ranking and thus AUC scores, use
+    `method="logistic"` which is a strictly monotonic transformation and thus keeps
+    the ranking.
 
 Multiclass support
 ^^^^^^^^^^^^^^^^^^
@@ -239,7 +252,7 @@ support 1-dimensional data (e.g., binary classification output) but are
 extended for multiclass classification if the `base_estimator` supports
 multiclass predictions. For multiclass predictions,
 :class:`CalibratedClassifierCV` calibrates for
-each class separately in a :ref:`ovr_classification` fashion [4]_. When
+each class separately in a :ref:`ovr_classification` fashion [5]_. When
 predicting
 probabilities, the calibrated probabilities for each class
 are predicted separately. As those probabilities do not necessarily sum to
@@ -254,36 +267,42 @@ one, a postprocessing is performed to normalize them.
 
 .. topic:: References:
 
-    .. [1] `Predicting Good Probabilities with Supervised Learning
-           <https://www.cs.cornell.edu/~alexn/papers/calibration.icml05.crc.rev3.pdf>`_,
-           A. Niculescu-Mizil & R. Caruana, ICML 2005
+    .. [1] Allan H. Murphy (1973).
+           :doi:`"A New Vector Partition of the Probability Score"
+           <10.1175/1520-0450(1973)012%3C0595:ANVPOT%3E2.0.CO;2>`
+           Journal of Applied Meteorology and Climatology
 
     .. [2] `On the combination of forecast probabilities for
            consecutive precipitation periods.
            <https://journals.ametsoc.org/waf/article/5/4/640/40179>`_
            Wea. Forecasting, 5, 640–650., Wilks, D. S., 1990a
 
-    .. [3] `Probabilistic Outputs for Support Vector Machines and Comparisons
+    .. [3] `Predicting Good Probabilities with Supervised Learning
+           <https://www.cs.cornell.edu/~alexn/papers/calibration.icml05.crc.rev3.pdf>`_,
+           A. Niculescu-Mizil & R. Caruana, ICML 2005
+
+
+    .. [4] `Probabilistic Outputs for Support Vector Machines and Comparisons
            to Regularized Likelihood Methods.
            <https://www.cs.colorado.edu/~mozer/Teaching/syllabi/6622/papers/Platt1999.pdf>`_
            J. Platt, (1999)
 
-    .. [4] `Transforming Classifier Scores into Accurate Multiclass
+    .. [5] `Transforming Classifier Scores into Accurate Multiclass
            Probability Estimates.
            <https://dl.acm.org/doi/pdf/10.1145/775047.775151>`_
            B. Zadrozny & C. Elkan, (KDD 2002)
 
-    .. [5] `Predicting accurate probabilities with a ranking loss.
+    .. [6] `Predicting accurate probabilities with a ranking loss.
            <https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4180410/>`_
            Menon AK, Jiang XJ, Vembu S, Elkan C, Ohno-Machado L.
            Proc Int Conf Mach Learn. 2012;2012:703-710
 
-    .. [6] `Beyond sigmoids: How to obtain well-calibrated probabilities from
+    .. [7] `Beyond sigmoids: How to obtain well-calibrated probabilities from
            binary classifiers with beta calibration
            <https://projecteuclid.org/euclid.ejs/1513306867>`_
            Kull, M., Silva Filho, T. M., & Flach, P. (2017).
 
-    .. [7] Mario V. Wüthrich, Michael Merz (2023).
+    .. [8] Mario V. Wüthrich, Michael Merz (2023).
            :doi:`"Statistical Foundations of Actuarial Learning and its Applications"
            <10.1007/978-3-031-12409-9>`
            Springer Actuarial


Code review comments:
diff_hunk: @@ -736,7 +736,7 @@ To also generate the example gallery you can use:
 
 .. prompt:: bash $
 
-    make html
+    make doc
path: doc/developers/contributing.rst
text: @user1:
I don't think this is correct for the `Makefile` of the `/doc` folder. This is only if you build the doc from the root directory. Above we state that you should be in the `/doc` folder.

@author:
I read too fast. But I never go into that path, I always use make doc:smirk: