.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/model_selection/plot_grid_search_stats.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_grid_search_stats.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_grid_search_stats.py:
==================================================
Statistical comparison of models using grid search
==================================================
This example illustrates how to statistically compare the performance of models
trained and evaluated using :class:`~sklearn.model_selection.GridSearchCV`.
.. GENERATED FROM PYTHON SOURCE LINES 12-16
We will start by simulating moon shaped data (where the ideal separation
between classes is non-linear), adding to it a moderate degree of noise.
Datapoints will belong to one of two possible classes to be predicted by two
features. We will simulate 50 samples for each class:
.. GENERATED FROM PYTHON SOURCE LINES 16-29
.. code-block:: default
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import make_moons
X, y = make_moons(noise=0.352, random_state=1, n_samples=100)
sns.scatterplot(
x=X[:, 0], y=X[:, 1], hue=y, marker="o", s=25, edgecolor="k", legend=False
).set_title("Data")
plt.show()
.. image-sg:: /auto_examples/model_selection/images/sphx_glr_plot_grid_search_stats_001.png
:alt: Data
:srcset: /auto_examples/model_selection/images/sphx_glr_plot_grid_search_stats_001.png
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
.. code-block:: none
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1498: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1498: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1498: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1498: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1498: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
.. GENERATED FROM PYTHON SOURCE LINES 30-38
We will compare the performance of :class:`~sklearn.svm.SVC` estimators that
vary on their `kernel` parameter, to decide which choice of this
hyper-parameter predicts our simulated data best.
We will evaluate the performance of the models using
:class:`~sklearn.model_selection.RepeatedStratifiedKFold`, repeating 10 times
a 10-fold stratified cross validation using a different randomization of the
data in each repetition. The performance will be evaluated using
:class:`~sklearn.metrics.roc_auc_score`.
.. GENERATED FROM PYTHON SOURCE LINES 38-55
.. code-block:: default
from sklearn.model_selection import GridSearchCV, RepeatedStratifiedKFold
from sklearn.svm import SVC
param_grid = [
{"kernel": ["linear"]},
{"kernel": ["poly"], "degree": [2, 3]},
{"kernel": ["rbf"]},
]
svc = SVC(random_state=0)
cv = RepeatedStratifiedKFold(n_splits=10, n_repeats=10, random_state=0)
search = GridSearchCV(estimator=svc, param_grid=param_grid, scoring="roc_auc", cv=cv)
search.fit(X, y)
.. raw:: html
<div class="output_subarea output_html rendered_html output_result">
<style>#sk-container-id-57 {color: black;}#sk-container-id-57 pre{padding: 0;}#sk-container-id-57 div.sk-toggleable {background-color: white;}#sk-container-id-57 label.sk-toggleable__label {cursor: pointer;display: block;width: 100%;margin-bottom: 0;padding: 0.3em;box-sizing: border-box;text-align: center;}#sk-container-id-57 label.sk-toggleable__label-arrow:before {content: "▸";float: left;margin-right: 0.25em;color: #696969;}#sk-container-id-57 label.sk-toggleable__label-arrow:hover:before {color: black;}#sk-container-id-57 div.sk-estimator:hover label.sk-toggleable__label-arrow:before {color: black;}#sk-container-id-57 div.sk-toggleable__content {max-height: 0;max-width: 0;overflow: hidden;text-align: left;background-color: #f0f8ff;}#sk-container-id-57 div.sk-toggleable__content pre {margin: 0.2em;color: black;border-radius: 0.25em;background-color: #f0f8ff;}#sk-container-id-57 input.sk-toggleable__control:checked~div.sk-toggleable__content {max-height: 200px;max-width: 100%;overflow: auto;}#sk-container-id-57 input.sk-toggleable__control:checked~label.sk-toggleable__label-arrow:before {content: "▾";}#sk-container-id-57 div.sk-estimator input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-57 div.sk-label input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-57 input.sk-hidden--visually {border: 0;clip: rect(1px 1px 1px 1px);clip: rect(1px, 1px, 1px, 1px);height: 1px;margin: -1px;overflow: hidden;padding: 0;position: absolute;width: 1px;}#sk-container-id-57 div.sk-estimator {font-family: monospace;background-color: #f0f8ff;border: 1px dotted black;border-radius: 0.25em;box-sizing: border-box;margin-bottom: 0.5em;}#sk-container-id-57 div.sk-estimator:hover {background-color: #d4ebff;}#sk-container-id-57 div.sk-parallel-item::after {content: "";width: 100%;border-bottom: 1px solid gray;flex-grow: 1;}#sk-container-id-57 div.sk-label:hover label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-57 div.sk-serial::before {content: "";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: 0;}#sk-container-id-57 div.sk-serial {display: flex;flex-direction: column;align-items: center;background-color: white;padding-right: 0.2em;padding-left: 0.2em;position: relative;}#sk-container-id-57 div.sk-item {position: relative;z-index: 1;}#sk-container-id-57 div.sk-parallel {display: flex;align-items: stretch;justify-content: center;background-color: white;position: relative;}#sk-container-id-57 div.sk-item::before, #sk-container-id-57 div.sk-parallel-item::before {content: "";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: -1;}#sk-container-id-57 div.sk-parallel-item {display: flex;flex-direction: column;z-index: 1;position: relative;background-color: white;}#sk-container-id-57 div.sk-parallel-item:first-child::after {align-self: flex-end;width: 50%;}#sk-container-id-57 div.sk-parallel-item:last-child::after {align-self: flex-start;width: 50%;}#sk-container-id-57 div.sk-parallel-item:only-child::after {width: 0;}#sk-container-id-57 div.sk-dashed-wrapped {border: 1px dashed gray;margin: 0 0.4em 0.5em 0.4em;box-sizing: border-box;padding-bottom: 0.4em;background-color: white;}#sk-container-id-57 div.sk-label label {font-family: monospace;font-weight: bold;display: inline-block;line-height: 1.2em;}#sk-container-id-57 div.sk-label-container {text-align: center;}#sk-container-id-57 div.sk-container {/* jupyter's `normalize.less` sets `[hidden] { display: none; }` but bootstrap.min.css set `[hidden] { display: none !important; }` so we also need the `!important` here to be able to override the default hidden behavior on the sphinx rendered scikit-learn.org. See: https://github.com/scikit-learn/scikit-learn/issues/21755 */display: inline-block !important;position: relative;}#sk-container-id-57 div.sk-text-repr-fallback {display: none;}</style><div id="sk-container-id-57" class="sk-top-container"><div class="sk-text-repr-fallback"><pre>GridSearchCV(cv=RepeatedStratifiedKFold(n_repeats=10, n_splits=10, random_state=0),
estimator=SVC(random_state=0),
param_grid=[{'kernel': ['linear']},
{'degree': [2, 3], 'kernel': ['poly']},
{'kernel': ['rbf']}],
scoring='roc_auc')</pre><b>In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. <br />On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.</b></div><div class="sk-container" hidden><div class="sk-item sk-dashed-wrapped"><div class="sk-label-container"><div class="sk-label sk-toggleable"><input class="sk-toggleable__control sk-hidden--visually" id="sk-estimator-id-241" type="checkbox" ><label for="sk-estimator-id-241" class="sk-toggleable__label sk-toggleable__label-arrow">GridSearchCV</label><div class="sk-toggleable__content"><pre>GridSearchCV(cv=RepeatedStratifiedKFold(n_repeats=10, n_splits=10, random_state=0),
estimator=SVC(random_state=0),
param_grid=[{'kernel': ['linear']},
{'degree': [2, 3], 'kernel': ['poly']},
{'kernel': ['rbf']}],
scoring='roc_auc')</pre></div></div></div><div class="sk-parallel"><div class="sk-parallel-item"><div class="sk-item"><div class="sk-label-container"><div class="sk-label sk-toggleable"><input class="sk-toggleable__control sk-hidden--visually" id="sk-estimator-id-242" type="checkbox" ><label for="sk-estimator-id-242" class="sk-toggleable__label sk-toggleable__label-arrow">estimator: SVC</label><div class="sk-toggleable__content"><pre>SVC(random_state=0)</pre></div></div></div><div class="sk-serial"><div class="sk-item"><div class="sk-estimator sk-toggleable"><input class="sk-toggleable__control sk-hidden--visually" id="sk-estimator-id-243" type="checkbox" ><label for="sk-estimator-id-243" class="sk-toggleable__label sk-toggleable__label-arrow">SVC</label><div class="sk-toggleable__content"><pre>SVC(random_state=0)</pre></div></div></div></div></div></div></div></div></div></div>
</div>
<br />
<br />
.. GENERATED FROM PYTHON SOURCE LINES 56-58
We can now inspect the results of our search, sorted by their
`mean_test_score`:
.. GENERATED FROM PYTHON SOURCE LINES 58-68
.. code-block:: default
import pandas as pd
results_df = pd.DataFrame(search.cv_results_)
results_df = results_df.sort_values(by=["rank_test_score"])
results_df = results_df.set_index(
results_df["params"].apply(lambda x: "_".join(str(val) for val in x.values()))
).rename_axis("kernel")
results_df[["params", "rank_test_score", "mean_test_score", "std_test_score"]]
.. 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>params</th>
<th>rank_test_score</th>
<th>mean_test_score</th>
<th>std_test_score</th>
</tr>
<tr>
<th>kernel</th>
<th></th>
<th></th>
<th></th>
<th></th>
</tr>
</thead>
<tbody>
<tr>
<th>rbf</th>
<td>{'kernel': 'rbf'}</td>
<td>1</td>
<td>0.9400</td>
<td>0.079297</td>
</tr>
<tr>
<th>linear</th>
<td>{'kernel': 'linear'}</td>
<td>2</td>
<td>0.9300</td>
<td>0.077846</td>
</tr>
<tr>
<th>3_poly</th>
<td>{'degree': 3, 'kernel': 'poly'}</td>
<td>3</td>
<td>0.9044</td>
<td>0.098776</td>
</tr>
<tr>
<th>2_poly</th>
<td>{'degree': 2, 'kernel': 'poly'}</td>
<td>4</td>
<td>0.6852</td>
<td>0.169106</td>
</tr>
</tbody>
</table>
</div>
</div>
<br />
<br />
.. GENERATED FROM PYTHON SOURCE LINES 69-92
We can see that the estimator using the `'rbf'` kernel performed best,
closely followed by `'linear'`. Both estimators with a `'poly'` kernel
performed worse, with the one using a two-degree polynomial achieving a much
lower performance than all other models.
Usually, the analysis just ends here, but half the story is missing. The
output of :class:`~sklearn.model_selection.GridSearchCV` does not provide
information on the certainty of the differences between the models.
We don't know if these are **statistically** significant.
To evaluate this, we need to conduct a statistical test.
Specifically, to contrast the performance of two models we should
statistically compare their AUC scores. There are 100 samples (AUC
scores) for each model as we repreated 10 times a 10-fold cross-validation.
However, the scores of the models are not independent: all models are
evaluated on the **same** 100 partitions, increasing the correlation
between the performance of the models.
Since some partitions of the data can make the distinction of the classes
particularly easy or hard to find for all models, the models scores will
co-vary.
Let's inspect this partition effect by plotting the performance of all models
in each fold, and calculating the correlation between models across folds:
.. GENERATED FROM PYTHON SOURCE LINES 92-114
.. code-block:: default
# create df of model scores ordered by performance
model_scores = results_df.filter(regex=r"split\d*_test_score")
# plot 30 examples of dependency between cv fold and AUC scores
fig, ax = plt.subplots()
sns.lineplot(
data=model_scores.transpose().iloc[:30],
dashes=False,
palette="Set1",
marker="o",
alpha=0.5,
ax=ax,
)
ax.set_xlabel("CV test fold", size=12, labelpad=10)
ax.set_ylabel("Model AUC", size=12)
ax.tick_params(bottom=True, labelbottom=False)
plt.show()
# print correlation of AUC scores across folds
print(f"Correlation of models:\n {model_scores.transpose().corr()}")
.. image-sg:: /auto_examples/model_selection/images/sphx_glr_plot_grid_search_stats_002.png
:alt: plot grid search stats
:srcset: /auto_examples/model_selection/images/sphx_glr_plot_grid_search_stats_002.png
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
.. code-block:: none
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1498: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1498: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1498: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1498: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1498: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1498: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1498: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1498: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1498: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1498: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1498: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1498: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1498: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1498: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1498: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1119: FutureWarning:
use_inf_as_na option is deprecated and will be removed in a future version. Convert inf values to NaN before operating instead.
/home/circleci/mambaforge/envs/testenv/lib/python3.9/site-packages/seaborn/_oldcore.py:1119: FutureWarning:
use_inf_as_na option is deprecated and will be removed in a future version. Convert inf values to NaN before operating instead.
Correlation of models:
kernel rbf linear 3_poly 2_poly
kernel
rbf 1.000000 0.882561 0.783392 0.351390
linear 0.882561 1.000000 0.746492 0.298688
3_poly 0.783392 0.746492 1.000000 0.355440
2_poly 0.351390 0.298688 0.355440 1.000000
.. GENERATED FROM PYTHON SOURCE LINES 115-126
We can observe that the performance of the models highly depends on the fold.
As a consequence, if we assume independence between samples we will be
underestimating the variance computed in our statistical tests, increasing
the number of false positive errors (i.e. detecting a significant difference
between models when such does not exist) [1]_.
Several variance-corrected statistical tests have been developed for these
cases. In this example we will show how to implement one of them (the so
called Nadeau and Bengio's corrected t-test) under two different statistical
frameworks: frequentist and Bayesian.
.. GENERATED FROM PYTHON SOURCE LINES 128-164
Comparing two models: frequentist approach
------------------------------------------
We can start by asking: "Is the first model significantly better than the
second model (when ranked by `mean_test_score`)?"
To answer this question using a frequentist approach we could
run a paired t-test and compute the p-value. This is also known as
Diebold-Mariano test in the forecast literature [5]_.
Many variants of such a t-test have been developed to account for the
'non-independence of samples problem'
described in the previous section. We will use the one proven to obtain the
highest replicability scores (which rate how similar the performance of a
model is when evaluating it on different random partitions of the same
dataset) while maintaining a low rate of false positives and false negatives:
the Nadeau and Bengio's corrected t-test [2]_ that uses a 10 times repeated
10-fold cross validation [3]_.
This corrected paired t-test is computed as:
.. math::
t=\frac{\frac{1}{k \cdot r}\sum_{i=1}^{k}\sum_{j=1}^{r}x_{ij}}
{\sqrt{(\frac{1}{k \cdot r}+\frac{n_{test}}{n_{train}})\hat{\sigma}^2}}
where :math:`k` is the number of folds,
:math:`r` the number of repetitions in the cross-validation,
:math:`x` is the difference in performance of the models,
:math:`n_{test}` is the number of samples used for testing,
:math:`n_{train}` is the number of samples used for training,
and :math:`\hat{\sigma}^2` represents the variance of the observed
differences.
Let's implement a corrected right-tailed paired t-test to evaluate if the
performance of the first model is significantly better than that of the
second model. Our null hypothesis is that the second model performs at least
as good as the first model.
.. GENERATED FROM PYTHON SOURCE LINES 164-222
.. code-block:: default
import numpy as np
from scipy.stats import t
def corrected_std(differences, n_train, n_test):
"""Corrects standard deviation using Nadeau and Bengio's approach.
Parameters
----------
differences : ndarray of shape (n_samples,)
Vector containing the differences in the score metrics of two models.
n_train : int
Number of samples in the training set.
n_test : int
Number of samples in the testing set.
Returns
-------
corrected_std : float
Variance-corrected standard deviation of the set of differences.
"""
# kr = k times r, r times repeated k-fold crossvalidation,
# kr equals the number of times the model was evaluated
kr = len(differences)
corrected_var = np.var(differences, ddof=1) * (1 / kr + n_test / n_train)
corrected_std = np.sqrt(corrected_var)
return corrected_std
def compute_corrected_ttest(differences, df, n_train, n_test):
"""Computes right-tailed paired t-test with corrected variance.
Parameters
----------
differences : array-like of shape (n_samples,)
Vector containing the differences in the score metrics of two models.
df : int
Degrees of freedom.
n_train : int
Number of samples in the training set.
n_test : int
Number of samples in the testing set.
Returns
-------
t_stat : float
Variance-corrected t-statistic.
p_val : float
Variance-corrected p-value.
"""
mean = np.mean(differences)
std = corrected_std(differences, n_train, n_test)
t_stat = mean / std
p_val = t.sf(np.abs(t_stat), df) # right-tailed t-test
return t_stat, p_val
.. GENERATED FROM PYTHON SOURCE LINES 223-236
.. code-block:: default
model_1_scores = model_scores.iloc[0].values # scores of the best model
model_2_scores = model_scores.iloc[1].values # scores of the second-best model
differences = model_1_scores - model_2_scores
n = differences.shape[0] # number of test sets
df = n - 1
n_train = len(list(cv.split(X, y))[0][0])
n_test = len(list(cv.split(X, y))[0][1])
t_stat, p_val = compute_corrected_ttest(differences, df, n_train, n_test)
print(f"Corrected t-value: {t_stat:.3f}\nCorrected p-value: {p_val:.3f}")
.. rst-class:: sphx-glr-script-out
.. code-block:: none
Corrected t-value: 0.750
Corrected p-value: 0.227
.. GENERATED FROM PYTHON SOURCE LINES 237-238
We can compare the corrected t- and p-values with the uncorrected ones:
.. GENERATED FROM PYTHON SOURCE LINES 238-247
.. code-block:: default
t_stat_uncorrected = np.mean(differences) / np.sqrt(np.var(differences, ddof=1) / n)
p_val_uncorrected = t.sf(np.abs(t_stat_uncorrected), df)
print(
f"Uncorrected t-value: {t_stat_uncorrected:.3f}\n"
f"Uncorrected p-value: {p_val_uncorrected:.3f}"
)
.. rst-class:: sphx-glr-script-out
.. code-block:: none
Uncorrected t-value: 2.611
Uncorrected p-value: 0.005
.. GENERATED FROM PYTHON SOURCE LINES 248-257
Using the conventional significance alpha level at `p=0.05`, we observe that
the uncorrected t-test concludes that the first model is significantly better
than the second.
With the corrected approach, in contrast, we fail to detect this difference.
In the latter case, however, the frequentist approach does not let us
conclude that the first and second model have an equivalent performance. If
we wanted to make this assertion we need to use a Bayesian approach.
.. GENERATED FROM PYTHON SOURCE LINES 259-300
Comparing two models: Bayesian approach
---------------------------------------
We can use Bayesian estimation to calculate the probability that the first
model is better than the second. Bayesian estimation will output a
distribution followed by the mean :math:`\mu` of the differences in the
performance of two models.
To obtain the posterior distribution we need to define a prior that models
our beliefs of how the mean is distributed before looking at the data,
and multiply it by a likelihood function that computes how likely our
observed differences are, given the values that the mean of differences
could take.
Bayesian estimation can be carried out in many forms to answer our question,
but in this example we will implement the approach suggested by Benavoli and
colleagues [4]_.
One way of defining our posterior using a closed-form expression is to select
a prior conjugate to the likelihood function. Benavoli and colleagues [4]_
show that when comparing the performance of two classifiers we can model the
prior as a Normal-Gamma distribution (with both mean and variance unknown)
conjugate to a normal likelihood, to thus express the posterior as a normal
distribution.
Marginalizing out the variance from this normal posterior, we can define the
posterior of the mean parameter as a Student's t-distribution. Specifically:
.. math::
St(\mu;n-1,\overline{x},(\frac{1}{n}+\frac{n_{test}}{n_{train}})
\hat{\sigma}^2)
where :math:`n` is the total number of samples,
:math:`\overline{x}` represents the mean difference in the scores,
:math:`n_{test}` is the number of samples used for testing,
:math:`n_{train}` is the number of samples used for training,
and :math:`\hat{\sigma}^2` represents the variance of the observed
differences.
Notice that we are using Nadeau and Bengio's corrected variance in our
Bayesian approach as well.
Let's compute and plot the posterior:
.. GENERATED FROM PYTHON SOURCE LINES 300-306
.. code-block:: default
# initialize random variable
t_post = t(
df, loc=np.mean(differences), scale=corrected_std(differences, n_train, n_test)
)
.. GENERATED FROM PYTHON SOURCE LINES 307-308
Let's plot the posterior distribution:
.. GENERATED FROM PYTHON SOURCE LINES 308-319
.. code-block:: default
x = np.linspace(t_post.ppf(0.001), t_post.ppf(0.999), 100)
plt.plot(x, t_post.pdf(x))
plt.xticks(np.arange(-0.04, 0.06, 0.01))
plt.fill_between(x, t_post.pdf(x), 0, facecolor="blue", alpha=0.2)
plt.ylabel("Probability density")
plt.xlabel(r"Mean difference ($\mu$)")
plt.title("Posterior distribution")
plt.show()
.. image-sg:: /auto_examples/model_selection/images/sphx_glr_plot_grid_search_stats_003.png
:alt: Posterior distribution
:srcset: /auto_examples/model_selection/images/sphx_glr_plot_grid_search_stats_003.png
:class: sphx-glr-single-img
.. GENERATED FROM PYTHON SOURCE LINES 320-325
We can calculate the probability that the first model is better than the
second by computing the area under the curve of the posterior distribution
from zero to infinity. And also the reverse: we can calculate the probability
that the second model is better than the first by computing the area under
the curve from minus infinity to zero.
.. GENERATED FROM PYTHON SOURCE LINES 325-337
.. code-block:: default
better_prob = 1 - t_post.cdf(0)
print(
f"Probability of {model_scores.index[0]} being more accurate than "
f"{model_scores.index[1]}: {better_prob:.3f}"
)
print(
f"Probability of {model_scores.index[1]} being more accurate than "
f"{model_scores.index[0]}: {1 - better_prob:.3f}"
)
.. rst-class:: sphx-glr-script-out
.. code-block:: none
Probability of rbf being more accurate than linear: 0.773
Probability of linear being more accurate than rbf: 0.227
.. GENERATED FROM PYTHON SOURCE LINES 338-344
In contrast with the frequentist approach, we can compute the probability
that one model is better than the other.
Note that we obtained similar results as those in the frequentist approach.
Given our choice of priors, we are essentially performing the same
computations, but we are allowed to make different assertions.
.. GENERATED FROM PYTHON SOURCE LINES 346-365
Region of Practical Equivalence
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Sometimes we are interested in determining the probabilities that our models
have an equivalent performance, where "equivalent" is defined in a practical
way. A naive approach [4]_ would be to define estimators as practically
equivalent when they differ by less than 1% in their accuracy. But we could
also define this practical equivalence taking into account the problem we are
trying to solve. For example, a difference of 5% in accuracy would mean an
increase of $1000 in sales, and we consider any quantity above that as
relevant for our business.
In this example we are going to define the
Region of Practical Equivalence (ROPE) to be :math:`[-0.01, 0.01]`. That is,
we will consider two models as practically equivalent if they differ by less
than 1% in their performance.
To compute the probabilities of the classifiers being practically equivalent,
we calculate the area under the curve of the posterior over the ROPE
interval:
.. GENERATED FROM PYTHON SOURCE LINES 365-374
.. code-block:: default
rope_interval = [-0.01, 0.01]
rope_prob = t_post.cdf(rope_interval[1]) - t_post.cdf(rope_interval[0])
print(
f"Probability of {model_scores.index[0]} and {model_scores.index[1]} "
f"being practically equivalent: {rope_prob:.3f}"
)
.. rst-class:: sphx-glr-script-out
.. code-block:: none
Probability of rbf and linear being practically equivalent: 0.432
.. GENERATED FROM PYTHON SOURCE LINES 375-376
We can plot how the posterior is distributed over the ROPE interval:
.. GENERATED FROM PYTHON SOURCE LINES 376-388
.. code-block:: default
x_rope = np.linspace(rope_interval[0], rope_interval[1], 100)
plt.plot(x, t_post.pdf(x))
plt.xticks(np.arange(-0.04, 0.06, 0.01))
plt.vlines([-0.01, 0.01], ymin=0, ymax=(np.max(t_post.pdf(x)) + 1))
plt.fill_between(x_rope, t_post.pdf(x_rope), 0, facecolor="blue", alpha=0.2)
plt.ylabel("Probability density")
plt.xlabel(r"Mean difference ($\mu$)")
plt.title("Posterior distribution under the ROPE")
plt.show()
.. image-sg:: /auto_examples/model_selection/images/sphx_glr_plot_grid_search_stats_004.png
:alt: Posterior distribution under the ROPE
:srcset: /auto_examples/model_selection/images/sphx_glr_plot_grid_search_stats_004.png
:class: sphx-glr-single-img
.. GENERATED FROM PYTHON SOURCE LINES 389-393
As suggested in [4]_, we can further interpret these probabilities using the
same criteria as the frequentist approach: is the probability of falling
inside the ROPE bigger than 95% (alpha value of 5%)? In that case we can
conclude that both models are practically equivalent.
.. GENERATED FROM PYTHON SOURCE LINES 395-405
The Bayesian estimation approach also allows us to compute how uncertain we
are about our estimation of the difference. This can be calculated using
credible intervals. For a given probability, they show the range of values
that the estimated quantity, in our case the mean difference in
performance, can take.
For example, a 50% credible interval [x, y] tells us that there is a 50%
probability that the true (mean) difference of performance between models is
between x and y.
Let's determine the credible intervals of our data using 50%, 75% and 95%:
.. GENERATED FROM PYTHON SOURCE LINES 405-418
.. code-block:: default
cred_intervals = []
intervals = [0.5, 0.75, 0.95]
for interval in intervals:
cred_interval = list(t_post.interval(interval))
cred_intervals.append([interval, cred_interval[0], cred_interval[1]])
cred_int_df = pd.DataFrame(
cred_intervals, columns=["interval", "lower value", "upper value"]
).set_index("interval")
cred_int_df
.. 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>lower value</th>
<th>upper value</th>
</tr>
<tr>
<th>interval</th>
<th></th>
<th></th>
</tr>
</thead>
<tbody>
<tr>
<th>0.50</th>
<td>0.000977</td>
<td>0.019023</td>
</tr>
<tr>
<th>0.75</th>
<td>-0.005422</td>
<td>0.025422</td>
</tr>
<tr>
<th>0.95</th>
<td>-0.016445</td>
<td>0.036445</td>
</tr>
</tbody>
</table>
</div>
</div>
<br />
<br />
.. GENERATED FROM PYTHON SOURCE LINES 419-423
As shown in the table, there is a 50% probability that the true mean
difference between models will be between 0.000977 and 0.019023, 70%
probability that it will be between -0.005422 and 0.025422, and 95%
probability that it will be between -0.016445 and 0.036445.
.. GENERATED FROM PYTHON SOURCE LINES 425-441
Pairwise comparison of all models: frequentist approach
-------------------------------------------------------
We could also be interested in comparing the performance of all our models
evaluated with :class:`~sklearn.model_selection.GridSearchCV`. In this case
we would be running our statistical test multiple times, which leads us to
the `multiple comparisons problem
<https://en.wikipedia.org/wiki/Multiple_comparisons_problem>`_.
There are many possible ways to tackle this problem, but a standard approach
is to apply a `Bonferroni correction
<https://en.wikipedia.org/wiki/Bonferroni_correction>`_. Bonferroni can be
computed by multiplying the p-value by the number of comparisons we are
testing.
Let's compare the performance of the models using the corrected t-test:
.. GENERATED FROM PYTHON SOURCE LINES 441-467
.. code-block:: default
from itertools import combinations
from math import factorial
n_comparisons = factorial(len(model_scores)) / (
factorial(2) * factorial(len(model_scores) - 2)
)
pairwise_t_test = []
for model_i, model_k in combinations(range(len(model_scores)), 2):
model_i_scores = model_scores.iloc[model_i].values
model_k_scores = model_scores.iloc[model_k].values
differences = model_i_scores - model_k_scores
t_stat, p_val = compute_corrected_ttest(differences, df, n_train, n_test)
p_val *= n_comparisons # implement Bonferroni correction
# Bonferroni can output p-values higher than 1
p_val = 1 if p_val > 1 else p_val
pairwise_t_test.append(
[model_scores.index[model_i], model_scores.index[model_k], t_stat, p_val]
)
pairwise_comp_df = pd.DataFrame(
pairwise_t_test, columns=["model_1", "model_2", "t_stat", "p_val"]
).round(3)
pairwise_comp_df
.. 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>model_1</th>
<th>model_2</th>
<th>t_stat</th>
<th>p_val</th>
</tr>
</thead>
<tbody>
<tr>
<th>0</th>
<td>rbf</td>
<td>linear</td>
<td>0.750</td>
<td>1.000</td>
</tr>
<tr>
<th>1</th>
<td>rbf</td>
<td>3_poly</td>
<td>1.657</td>
<td>0.302</td>
</tr>
<tr>
<th>2</th>
<td>rbf</td>
<td>2_poly</td>
<td>4.565</td>
<td>0.000</td>
</tr>
<tr>
<th>3</th>
<td>linear</td>
<td>3_poly</td>
<td>1.111</td>
<td>0.807</td>
</tr>
<tr>
<th>4</th>
<td>linear</td>
<td>2_poly</td>
<td>4.276</td>
<td>0.000</td>
</tr>
<tr>
<th>5</th>
<td>3_poly</td>
<td>2_poly</td>
<td>3.851</td>
<td>0.001</td>
</tr>
</tbody>
</table>
</div>
</div>
<br />
<br />
.. GENERATED FROM PYTHON SOURCE LINES 468-473
We observe that after correcting for multiple comparisons, the only model
that significantly differs from the others is `'2_poly'`.
`'rbf'`, the model ranked first by
:class:`~sklearn.model_selection.GridSearchCV`, does not significantly
differ from `'linear'` or `'3_poly'`.
.. GENERATED FROM PYTHON SOURCE LINES 475-483
Pairwise comparison of all models: Bayesian approach
----------------------------------------------------
When using Bayesian estimation to compare multiple models, we don't need to
correct for multiple comparisons (for reasons why see [4]_).
We can carry out our pairwise comparisons the same way as in the first
section:
.. GENERATED FROM PYTHON SOURCE LINES 483-506
.. code-block:: default
pairwise_bayesian = []
for model_i, model_k in combinations(range(len(model_scores)), 2):
model_i_scores = model_scores.iloc[model_i].values
model_k_scores = model_scores.iloc[model_k].values
differences = model_i_scores - model_k_scores
t_post = t(
df, loc=np.mean(differences), scale=corrected_std(differences, n_train, n_test)
)
worse_prob = t_post.cdf(rope_interval[0])
better_prob = 1 - t_post.cdf(rope_interval[1])
rope_prob = t_post.cdf(rope_interval[1]) - t_post.cdf(rope_interval[0])
pairwise_bayesian.append([worse_prob, better_prob, rope_prob])
pairwise_bayesian_df = pd.DataFrame(
pairwise_bayesian, columns=["worse_prob", "better_prob", "rope_prob"]
).round(3)
pairwise_comp_df = pairwise_comp_df.join(pairwise_bayesian_df)
pairwise_comp_df
.. 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>model_1</th>
<th>model_2</th>
<th>t_stat</th>
<th>p_val</th>
<th>worse_prob</th>
<th>better_prob</th>
<th>rope_prob</th>
</tr>
</thead>
<tbody>
<tr>
<th>0</th>
<td>rbf</td>
<td>linear</td>
<td>0.750</td>
<td>1.000</td>
<td>0.068</td>
<td>0.500</td>
<td>0.432</td>
</tr>
<tr>
<th>1</th>
<td>rbf</td>
<td>3_poly</td>
<td>1.657</td>
<td>0.302</td>
<td>0.018</td>
<td>0.882</td>
<td>0.100</td>
</tr>
<tr>
<th>2</th>
<td>rbf</td>
<td>2_poly</td>
<td>4.565</td>
<td>0.000</td>
<td>0.000</td>
<td>1.000</td>
<td>0.000</td>
</tr>
<tr>
<th>3</th>
<td>linear</td>
<td>3_poly</td>
<td>1.111</td>
<td>0.807</td>
<td>0.063</td>
<td>0.750</td>
<td>0.187</td>
</tr>
<tr>
<th>4</th>
<td>linear</td>
<td>2_poly</td>
<td>4.276</td>
<td>0.000</td>
<td>0.000</td>
<td>1.000</td>
<td>0.000</td>
</tr>
<tr>
<th>5</th>
<td>3_poly</td>
<td>2_poly</td>
<td>3.851</td>
<td>0.001</td>
<td>0.000</td>
<td>1.000</td>
<td>0.000</td>
</tr>
</tbody>
</table>
</div>
</div>
<br />
<br />
.. GENERATED FROM PYTHON SOURCE LINES 507-520
Using the Bayesian approach we can compute the probability that a model
performs better, worse or practically equivalent to another.
Results show that the model ranked first by
:class:`~sklearn.model_selection.GridSearchCV` `'rbf'`, has approximately a
6.8% chance of being worse than `'linear'`, and a 1.8% chance of being worse
than `'3_poly'`.
`'rbf'` and `'linear'` have a 43% probability of being practically
equivalent, while `'rbf'` and `'3_poly'` have a 10% chance of being so.
Similarly to the conclusions obtained using the frequentist approach, all
models have a 100% probability of being better than `'2_poly'`, and none have
a practically equivalent performance with the latter.
.. GENERATED FROM PYTHON SOURCE LINES 522-541
Take-home messages
------------------
- Small differences in performance measures might easily turn out to be
merely by chance, but not because one model predicts systematically better
than the other. As shown in this example, statistics can tell you how
likely that is.
- When statistically comparing the performance of two models evaluated in
GridSearchCV, it is necessary to correct the calculated variance which
could be underestimated since the scores of the models are not independent
from each other.
- A frequentist approach that uses a (variance-corrected) paired t-test can
tell us if the performance of one model is better than another with a
degree of certainty above chance.
- A Bayesian approach can provide the probabilities of one model being
better, worse or practically equivalent than another. It can also tell us
how confident we are of knowing that the true differences of our models
fall under a certain range of values.
- If multiple models are statistically compared, a multiple comparisons
correction is needed when using the frequentist approach.
.. GENERATED FROM PYTHON SOURCE LINES 543-567
.. topic:: References
.. [1] Dietterich, T. G. (1998). `Approximate statistical tests for
comparing supervised classification learning algorithms
<http://web.cs.iastate.edu/~jtian/cs573/Papers/Dietterich-98.pdf>`_.
Neural computation, 10(7).
.. [2] Nadeau, C., & Bengio, Y. (2000). `Inference for the generalization
error
<https://papers.nips.cc/paper/1661-inference-for-the-generalization-error.pdf>`_.
In Advances in neural information processing systems.
.. [3] Bouckaert, R. R., & Frank, E. (2004). `Evaluating the replicability
of significance tests for comparing learning algorithms
<https://www.cms.waikato.ac.nz/~ml/publications/2004/bouckaert-frank.pdf>`_.
In Pacific-Asia Conference on Knowledge Discovery and Data Mining.
.. [4] Benavoli, A., Corani, G., Demšar, J., & Zaffalon, M. (2017). `Time
for a change: a tutorial for comparing multiple classifiers through
Bayesian analysis
<http://www.jmlr.org/papers/volume18/16-305/16-305.pdf>`_.
The Journal of Machine Learning Research, 18(1). See the Python
library that accompanies this paper `here
<https://github.com/janezd/baycomp>`_.
.. [5] Diebold, F.X. & Mariano R.S. (1995). `Comparing predictive accuracy
<http://www.est.uc3m.es/esp/nueva_docencia/comp_col_get/lade/tecnicas_prediccion/Practicas0708/Comparing%20Predictive%20Accuracy%20(Dielbold).pdf>`_
Journal of Business & economic statistics, 20(1), 134-144.
.. rst-class:: sphx-glr-timing
**Total running time of the script:** (0 minutes 1.466 seconds)
.. _sphx_glr_download_auto_examples_model_selection_plot_grid_search_stats.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/1.3.X?urlpath=lab/tree/notebooks/auto_examples/model_selection/plot_grid_search_stats.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_grid_search_stats.ipynb
:alt: Launch JupyterLite
:width: 150 px
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_grid_search_stats.py <plot_grid_search_stats.py>`
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_grid_search_stats.ipynb <plot_grid_search_stats.ipynb>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_