.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/model_selection/plot_likelihood_ratios.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_model_selection_plot_likelihood_ratios.py>`
        to download the full example code or to run this example in your browser via JupyterLite or Binder

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_model_selection_plot_likelihood_ratios.py:


=============================================================
Class Likelihood Ratios to measure classification performance
=============================================================

This example demonstrates the :func:`~sklearn.metrics.class_likelihood_ratios`
function, which computes the positive and negative likelihood ratios (`LR+`,
`LR-`) to assess the predictive power of a binary classifier. As we will see,
these metrics are independent of the proportion between classes in the test set,
which makes them very useful when the available data for a study has a different
class proportion than the target application.

A typical use is a case-control study in medicine, which has nearly balanced
classes while the general population has large class imbalance. In such
application, the pre-test probability of an individual having the target
condition can be chosen to be the prevalence, i.e. the proportion of a
particular population found to be affected by a medical condition. The post-test
probabilities represent then the probability that the condition is truly present
given a positive test result.

In this example we first discuss the link between pre-test and post-test odds
given by the :ref:`class_likelihood_ratios`. Then we evaluate their behavior in
some controlled scenarios. In the last section we plot them as a function of the
prevalence of the positive class.

.. GENERATED FROM PYTHON SOURCE LINES 27-30

.. code-block:: Python


    # Authors:  Arturo Amor <david-arturo.amor-quiroz@inria.fr>
    #           Olivier Grisel <olivier.grisel@ensta.org>







.. GENERATED FROM PYTHON SOURCE LINES 31-38

Pre-test vs. post-test analysis
===============================

Suppose we have a population of subjects with physiological measurements `X`
that can hopefully serve as indirect bio-markers of the disease and actual
disease indicators `y` (ground truth). Most of the people in the population do
not carry the disease but a minority (in this case around 10%) does:

.. GENERATED FROM PYTHON SOURCE LINES 38-44

.. code-block:: Python


    from sklearn.datasets import make_classification

    X, y = make_classification(n_samples=10_000, weights=[0.9, 0.1], random_state=0)
    print(f"Percentage of people carrying the disease: {100*y.mean():.2f}%")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Percentage of people carrying the disease: 10.37%




.. GENERATED FROM PYTHON SOURCE LINES 45-48

A machine learning model is built to diagnose if a person with some given
physiological measurements is likely to carry the disease of interest. To
evaluate the model, we need to assess its performance on a held-out test set:

.. GENERATED FROM PYTHON SOURCE LINES 48-53

.. code-block:: Python


    from sklearn.model_selection import train_test_split

    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)








.. GENERATED FROM PYTHON SOURCE LINES 54-57

Then we can fit our diagnosis model and compute the positive likelihood
ratio to evaluate the usefulness of this classifier as a disease diagnosis
tool:

.. GENERATED FROM PYTHON SOURCE LINES 57-66

.. code-block:: Python


    from sklearn.linear_model import LogisticRegression
    from sklearn.metrics import class_likelihood_ratios

    estimator = LogisticRegression().fit(X_train, y_train)
    y_pred = estimator.predict(X_test)
    pos_LR, neg_LR = class_likelihood_ratios(y_test, y_pred)
    print(f"LR+: {pos_LR:.3f}")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    LR+: 12.617




.. GENERATED FROM PYTHON SOURCE LINES 67-77

Since the positive class likelihood ratio is much larger than 1.0, it means
that the machine learning-based diagnosis tool is useful: the post-test odds
that the condition is truly present given a positive test result are more than
12 times larger than the pre-test odds.

Cross-validation of likelihood ratios
=====================================

We assess the variability of the measurements for the class likelihood ratios
in some particular cases.

.. GENERATED FROM PYTHON SOURCE LINES 77-97

.. code-block:: Python


    import pandas as pd


    def scoring(estimator, X, y):
        y_pred = estimator.predict(X)
        pos_lr, neg_lr = class_likelihood_ratios(y, y_pred, raise_warning=False)
        return {"positive_likelihood_ratio": pos_lr, "negative_likelihood_ratio": neg_lr}


    def extract_score(cv_results):
        lr = pd.DataFrame(
            {
                "positive": cv_results["test_positive_likelihood_ratio"],
                "negative": cv_results["test_negative_likelihood_ratio"],
            }
        )
        return lr.aggregate(["mean", "std"])









.. GENERATED FROM PYTHON SOURCE LINES 98-100

We first validate the :class:`~sklearn.linear_model.LogisticRegression` model
with default hyperparameters as used in the previous section.

.. GENERATED FROM PYTHON SOURCE LINES 100-106

.. code-block:: Python


    from sklearn.model_selection import cross_validate

    estimator = LogisticRegression()
    extract_score(cross_validate(estimator, X, y, scoring=scoring, cv=10))






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>positive</th>
          <th>negative</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>mean</th>
          <td>16.661086</td>
          <td>0.724702</td>
        </tr>
        <tr>
          <th>std</th>
          <td>4.383973</td>
          <td>0.054045</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 107-113

We confirm that the model is useful: the post-test odds are between 12 and 20
times larger than the pre-test odds.

On the contrary, let's consider a dummy model that will output random
predictions with similar odds as the average disease prevalence in the
training set:

.. GENERATED FROM PYTHON SOURCE LINES 113-119

.. code-block:: Python


    from sklearn.dummy import DummyClassifier

    estimator = DummyClassifier(strategy="stratified", random_state=1234)
    extract_score(cross_validate(estimator, X, y, scoring=scoring, cv=10))






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>positive</th>
          <th>negative</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>mean</th>
          <td>1.108843</td>
          <td>0.986989</td>
        </tr>
        <tr>
          <th>std</th>
          <td>0.268147</td>
          <td>0.034278</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 120-125

Here both class likelihood ratios are compatible with 1.0 which makes this
classifier useless as a diagnostic tool to improve disease detection.

Another option for the dummy model is to always predict the most frequent
class, which in this case is "no-disease".

.. GENERATED FROM PYTHON SOURCE LINES 125-129

.. code-block:: Python


    estimator = DummyClassifier(strategy="most_frequent")
    extract_score(cross_validate(estimator, X, y, scoring=scoring, cv=10))






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>positive</th>
          <th>negative</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>mean</th>
          <td>NaN</td>
          <td>1.0</td>
        </tr>
        <tr>
          <th>std</th>
          <td>NaN</td>
          <td>0.0</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 130-145

The absence of positive predictions means there will be no true positives nor
false positives, leading to an undefined `LR+` that by no means should be
interpreted as an infinite `LR+` (the classifier perfectly identifying
positive cases). In such situation the
:func:`~sklearn.metrics.class_likelihood_ratios` function returns `nan` and
raises a warning by default. Indeed, the value of `LR-` helps us discard this
model.

A similar scenario may arise when cross-validating highly imbalanced data with
few samples: some folds will have no samples with the disease and therefore
they will output no true positives nor false negatives when used for testing.
Mathematically this leads to an infinite `LR+`, which should also not be
interpreted as the model perfectly identifying positive cases. Such event
leads to a higher variance of the estimated likelihood ratios, but can still
be interpreted as an increment of the post-test odds of having the condition.

.. GENERATED FROM PYTHON SOURCE LINES 145-150

.. code-block:: Python


    estimator = LogisticRegression()
    X, y = make_classification(n_samples=300, weights=[0.9, 0.1], random_state=0)
    extract_score(cross_validate(estimator, X, y, scoring=scoring, cv=10))






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>positive</th>
          <th>negative</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>mean</th>
          <td>17.8000</td>
          <td>0.373333</td>
        </tr>
        <tr>
          <th>std</th>
          <td>8.5557</td>
          <td>0.235430</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 151-168

Invariance with respect to prevalence
=====================================

The likelihood ratios are independent of the disease prevalence and can be
extrapolated between populations regardless of any possible class imbalance,
**as long as the same model is applied to all of them**. Notice that in the
plots below **the decision boundary is constant** (see
:ref:`sphx_glr_auto_examples_svm_plot_separating_hyperplane_unbalanced.py` for
a study of the boundary decision for unbalanced classes).

Here we train a :class:`~sklearn.linear_model.LogisticRegression` base model
on a case-control study with a prevalence of 50%. It is then evaluated over
populations with varying prevalence. We use the
:func:`~sklearn.datasets.make_classification` function to ensure the
data-generating process is always the same as shown in the plots below. The
label `1` corresponds to the positive class "disease", whereas the label `0`
stands for "no-disease".

.. GENERATED FROM PYTHON SOURCE LINES 168-194

.. code-block:: Python


    from collections import defaultdict

    import matplotlib.pyplot as plt
    import numpy as np

    from sklearn.inspection import DecisionBoundaryDisplay

    populations = defaultdict(list)
    common_params = {
        "n_samples": 10_000,
        "n_features": 2,
        "n_informative": 2,
        "n_redundant": 0,
        "random_state": 0,
    }
    weights = np.linspace(0.1, 0.8, 6)
    weights = weights[::-1]

    # fit and evaluate base model on balanced classes
    X, y = make_classification(**common_params, weights=[0.5, 0.5])
    estimator = LogisticRegression().fit(X, y)
    lr_base = extract_score(cross_validate(estimator, X, y, scoring=scoring, cv=10))
    pos_lr_base, pos_lr_base_std = lr_base["positive"].values
    neg_lr_base, neg_lr_base_std = lr_base["negative"].values








.. GENERATED FROM PYTHON SOURCE LINES 195-198

We will now show the decision boundary for each level of prevalence. Note that
we only plot a subset of the original data to better assess the linear model
decision boundary.

.. GENERATED FROM PYTHON SOURCE LINES 198-228

.. code-block:: Python


    fig, axs = plt.subplots(nrows=3, ncols=2, figsize=(15, 12))

    for ax, (n, weight) in zip(axs.ravel(), enumerate(weights)):
        X, y = make_classification(
            **common_params,
            weights=[weight, 1 - weight],
        )
        prevalence = y.mean()
        populations["prevalence"].append(prevalence)
        populations["X"].append(X)
        populations["y"].append(y)

        # down-sample for plotting
        rng = np.random.RandomState(1)
        plot_indices = rng.choice(np.arange(X.shape[0]), size=500, replace=True)
        X_plot, y_plot = X[plot_indices], y[plot_indices]

        # plot fixed decision boundary of base model with varying prevalence
        disp = DecisionBoundaryDisplay.from_estimator(
            estimator,
            X_plot,
            response_method="predict",
            alpha=0.5,
            ax=ax,
        )
        scatter = disp.ax_.scatter(X_plot[:, 0], X_plot[:, 1], c=y_plot, edgecolor="k")
        disp.ax_.set_title(f"prevalence = {y_plot.mean():.2f}")
        disp.ax_.legend(*scatter.legend_elements())




.. image-sg:: /auto_examples/model_selection/images/sphx_glr_plot_likelihood_ratios_001.png
   :alt: prevalence = 0.22, prevalence = 0.34, prevalence = 0.45, prevalence = 0.60, prevalence = 0.76, prevalence = 0.88
   :srcset: /auto_examples/model_selection/images/sphx_glr_plot_likelihood_ratios_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 229-230

We define a function for bootstrapping.

.. GENERATED FROM PYTHON SOURCE LINES 230-245

.. code-block:: Python



    def scoring_on_bootstrap(estimator, X, y, rng, n_bootstrap=100):
        results_for_prevalence = defaultdict(list)
        for _ in range(n_bootstrap):
            bootstrap_indices = rng.choice(
                np.arange(X.shape[0]), size=X.shape[0], replace=True
            )
            for key, value in scoring(
                estimator, X[bootstrap_indices], y[bootstrap_indices]
            ).items():
                results_for_prevalence[key].append(value)
        return pd.DataFrame(results_for_prevalence)









.. GENERATED FROM PYTHON SOURCE LINES 246-247

We score the base model for each prevalence using bootstrapping.

.. GENERATED FROM PYTHON SOURCE LINES 247-267

.. code-block:: Python


    results = defaultdict(list)
    n_bootstrap = 100
    rng = np.random.default_rng(seed=0)

    for prevalence, X, y in zip(
        populations["prevalence"], populations["X"], populations["y"]
    ):
        results_for_prevalence = scoring_on_bootstrap(
            estimator, X, y, rng, n_bootstrap=n_bootstrap
        )
        results["prevalence"].append(prevalence)
        results["metrics"].append(
            results_for_prevalence.aggregate(["mean", "std"]).unstack()
        )

    results = pd.DataFrame(results["metrics"], index=results["prevalence"])
    results.index.name = "prevalence"
    results






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead tr th {
            text-align: left;
        }

        .dataframe thead tr:last-of-type th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr>
          <th></th>
          <th colspan="2" halign="left">positive_likelihood_ratio</th>
          <th colspan="2" halign="left">negative_likelihood_ratio</th>
        </tr>
        <tr>
          <th></th>
          <th>mean</th>
          <th>std</th>
          <th>mean</th>
          <th>std</th>
        </tr>
        <tr>
          <th>prevalence</th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0.2039</th>
          <td>4.507943</td>
          <td>0.113516</td>
          <td>0.207667</td>
          <td>0.009778</td>
        </tr>
        <tr>
          <th>0.3419</th>
          <td>4.443238</td>
          <td>0.125140</td>
          <td>0.198766</td>
          <td>0.008915</td>
        </tr>
        <tr>
          <th>0.4809</th>
          <td>4.421087</td>
          <td>0.123828</td>
          <td>0.192913</td>
          <td>0.006360</td>
        </tr>
        <tr>
          <th>0.6196</th>
          <td>4.409717</td>
          <td>0.164009</td>
          <td>0.193949</td>
          <td>0.005861</td>
        </tr>
        <tr>
          <th>0.7578</th>
          <td>4.334795</td>
          <td>0.175298</td>
          <td>0.189267</td>
          <td>0.005840</td>
        </tr>
        <tr>
          <th>0.8963</th>
          <td>4.197666</td>
          <td>0.238955</td>
          <td>0.185654</td>
          <td>0.005027</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 268-271

In the plots below we observe that the class likelihood ratios re-computed
with different prevalences are indeed constant within one standard deviation
of those computed with on balanced classes.

.. GENERATED FROM PYTHON SOURCE LINES 271-326

.. code-block:: Python


    fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(15, 6))
    results["positive_likelihood_ratio"]["mean"].plot(
        ax=ax1, color="r", label="extrapolation through populations"
    )
    ax1.axhline(y=pos_lr_base + pos_lr_base_std, color="r", linestyle="--")
    ax1.axhline(
        y=pos_lr_base - pos_lr_base_std,
        color="r",
        linestyle="--",
        label="base model confidence band",
    )
    ax1.fill_between(
        results.index,
        results["positive_likelihood_ratio"]["mean"]
        - results["positive_likelihood_ratio"]["std"],
        results["positive_likelihood_ratio"]["mean"]
        + results["positive_likelihood_ratio"]["std"],
        color="r",
        alpha=0.3,
    )
    ax1.set(
        title="Positive likelihood ratio",
        ylabel="LR+",
        ylim=[0, 5],
    )
    ax1.legend(loc="lower right")

    ax2 = results["negative_likelihood_ratio"]["mean"].plot(
        ax=ax2, color="b", label="extrapolation through populations"
    )
    ax2.axhline(y=neg_lr_base + neg_lr_base_std, color="b", linestyle="--")
    ax2.axhline(
        y=neg_lr_base - neg_lr_base_std,
        color="b",
        linestyle="--",
        label="base model confidence band",
    )
    ax2.fill_between(
        results.index,
        results["negative_likelihood_ratio"]["mean"]
        - results["negative_likelihood_ratio"]["std"],
        results["negative_likelihood_ratio"]["mean"]
        + results["negative_likelihood_ratio"]["std"],
        color="b",
        alpha=0.3,
    )
    ax2.set(
        title="Negative likelihood ratio",
        ylabel="LR-",
        ylim=[0, 0.5],
    )
    ax2.legend(loc="lower right")

    plt.show()



.. image-sg:: /auto_examples/model_selection/images/sphx_glr_plot_likelihood_ratios_002.png
   :alt: Positive likelihood ratio, Negative likelihood ratio
   :srcset: /auto_examples/model_selection/images/sphx_glr_plot_likelihood_ratios_002.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 3.589 seconds)


.. _sphx_glr_download_auto_examples_model_selection_plot_likelihood_ratios.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: binder-badge

      .. image:: images/binder_badge_logo.svg
        :target: https://mybinder.org/v2/gh/scikit-learn/scikit-learn/main?urlpath=lab/tree/notebooks/auto_examples/model_selection/plot_likelihood_ratios.ipynb
        :alt: Launch binder
        :width: 150 px

    .. container:: lite-badge

      .. image:: images/jupyterlite_badge_logo.svg
        :target: ../../lite/lab/?path=auto_examples/model_selection/plot_likelihood_ratios.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_likelihood_ratios.ipynb <plot_likelihood_ratios.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_likelihood_ratios.py <plot_likelihood_ratios.py>`


.. include:: plot_likelihood_ratios.recommendations


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_