Статистическое сравнение моделей с использованием поиска по сетке

Сравнение двух моделей: частотный подход

Можно начать с вопроса: «Является ли первая модель статистически значимо лучше второй модели (при ранжировании по mean_test_score)?»

Для ответа на этот вопрос с помощью частотного подхода можно выполнить парный t-тест и вычислить p-значение. Это также известно как тест Диболда-Мариано в литературе по прогнозированию [5]. Было разработано много вариантов такого t-теста, чтобы учесть проблему «независимости выборок», описанную в предыдущем разделе. Мы будем использовать тот, который показал лучшие результаты по воспроизводимости (что оценивает насколько похожа производительность модели при оценке её на разных случайных разбиениях того же набора данных), сохраняя низкий уровень ложноположительных и ложноотрицательных результатов: скорректированный t-тест Надьо и Бенжио [2], который использует 10-кратную перекрёстную проверку, повторяемую 10 раз [3].

Этот скорректированный парный t-тест вычисляется как:

где \(k\) — число итераций, \(r\) — число повторений в перекрестной проверке, \(x\) — разность в производительности моделей, \(n_{test}\) — число выборок, используемых для тестирования, \(n_{train}\) — число выборок, используемых для обучения, и \(\hat{\sigma}^2\) представляет собой дисперсию наблюдаемых различий.

Давайте реализуем скорректированный односторонний парный t-тест, чтобы оценить, является ли производительность первой модели статистически значимо лучше, чем у второй модели. Наша нулевая гипотеза заключается в том, что вторая модель работает по крайней мере так же хорошо, как и первая.

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

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}")

Corrected t-value: 0.750
Corrected p-value: 0.227

Давайте сравним скорректированные значения t и p с нескорректированными:

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}"
)

Uncorrected t-value: 2.611
Uncorrected p-value: 0.005

Используя обычный уровень значимости альфа в p=0.05, мы видим, что нескорректированный t-тест говорит о том, что первая модель значимо лучше второй.

В отличие от этого, скорректированный подход не позволяет сделать этого вывода.

В последнем случае, однако, частотный подход не позволяет сделать вывод о том, что первая и вторая модели имеют эквивалентную производительность. Если мы хотим сделать это утверждение, нам нужно использовать байесовский подход.

Сравнение двух моделей: байесовский подход

Мы можем использовать байесовскую оценку для расчета вероятности того, что первая модель лучше второй. Байесовская оценка выведет распределение, за которым следует среднее значение \(\mu\) разницы в производительности двух моделей.

Для получения апостериорного распределения нам нужно определить априорное распределение, которое моделирует наши убеждения о том, как распределяется среднее значение до просмотра данных, и умножить его на функцию правдоподобия, которая вычисляет, насколько вероятны наши наблюдаемые различия, учитывая значения, которые может принимать среднее значение различий.

Байесовская оценка может быть выполнена в различных формах для ответа на наш вопрос, но в этом примере мы будем использовать подход, предложенный Бенаволи и коллегами [4].

Один из способов определения нашего апостериорного распределения с помощью замкнутой формы — выбрать априорное распределение, сопряженное с функцией правдоподобия. Бенаволи и коллеги [4] показывают, что при сравнении производительности двух классификаторов мы можем смоделировать априорное распределение как нормально-гамма-распределение (с неизвестным средним значением и дисперсией), сопряженное с нормальным правдоподобием, чтобы выразить апостериорное распределение как нормальное распределение. Выполняя интегрирование по дисперсии из этого нормального апостериорного распределения, мы можем определить апостериорное распределение параметра среднего значения как распределение Стьюдента. В частности:

где \(n\) — общее количество образцов, \(\overline{x}\) представляет собой среднее значение разницы в баллах, \(n_{test}\) — количество образцов, используемых для тестирования, \(n_{train}\) — количество образцов, используемых для обучения, а \(\hat{\sigma}^2\) — дисперсия наблюдаемых различий.

Обратите внимание, что мы также используем скорректированную дисперсию Надо и Бенжио в нашем байесовском подходе.

Давайте вычислим и построим апостериорное распределение:

# initialize random variable
t_post = t(
    df, loc=np.mean(differences), scale=corrected_std(differences, n_train, n_test)
)

Давайте построим график апостериорного распределения:

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()

Мы можем вычислить вероятность того, что первая модель лучше второй, вычислив площадь под кривой апостериорного распределения от нуля до бесконечности. И наоборот, мы можем вычислить вероятность того, что вторая модель лучше первой, вычислив площадь под кривой от минус бесконечности до нуля.

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}"
)

Probability of rbf being more accurate than linear: 0.773
Probability of linear being more accurate than rbf: 0.227

В отличие от частотного подхода, мы можем вычислить вероятность того, что одна модель лучше другой.

Обратите внимание, что мы получили аналогичные результаты, как и в частотном подходе. Учитывая наш выбор априорных распределений, мы фактически выполняем те же вычисления, но нам разрешено делать разные утверждения.

Область практического эквивалента

Иногда мы заинтересованы в определении вероятностей того, что наши модели имеют эквивалентную производительность, где «эквивалентность» определяется практическим способом. Примитивный подход [4] заключался бы в определении оценщиков как практически эквивалентных, если они отличаются менее чем на 1% по точности. Но мы также можем определить эту практическую эквивалентность, учитывая задачу, которую мы пытаемся решить. Например, различие в 5% по точности означало бы увеличение продаж на 1000 долларов, и мы считаем любые значения выше этого значения релевантными для нашего бизнеса.

В этом примере мы определим область практического эквивалента (ROPE) как \([-0.01, 0.01]\). То есть мы будем считать две модели практически эквивалентными, если они отличаются менее чем на 1% по производительности.

Для вычисления вероятностей практической эквивалентности классификаторов мы вычисляем площадь под кривой апостериорного распределения на интервале ROPE:

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}"
)

Probability of rbf and linear being practically equivalent: 0.432

Мы можем построить график распределения апостериорного распределения по интервалу ROPE:

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()

Как предлагается в [4], мы можем далее интерпретировать эти вероятности, используя те же критерии, что и частотный подход: является ли вероятность попадания в ROPE больше 95% (значение альфа 5%)? В этом случае мы можем заключить, что обе модели практически эквивалентны.

Байесовский подход к оценке также позволяет нам вычислить, насколько не уверены мы в нашей оценке разницы. Это можно вычислить, используя доверительные интервалы. Для заданной вероятности они показывают диапазон значений, которые может принимать оцениваемая величина, в нашем случае среднее значение разницы в производительности. Например, 50%-ный доверительный интервал [x, y] говорит нам о том, что существует 50%-ная вероятность того, что истинное (среднее) различие в производительности между моделями находится между x и y.

Давайте определим доверительные интервалы наших данных, используя 50%, 75% и 95%:

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

	lower value	upper value
interval
0.50	0.000977	0.019023
0.75	-0.005422	0.025422
0.95	-0.016445	0.036445

Как показано в таблице, существует 50%-ная вероятность того, что истинное среднее различие между моделями будет находиться между 0,000977 и 0,019023, 70%-ная вероятность того, что оно будет находиться между -0,005422 и 0,025422, и 95%-ная вероятность того, что оно будет находиться между -0,016445 и 0,036445.

Парное сравнение всех моделей: частотный подход

Мы также можем быть заинтересованы в сравнении производительности всех наших моделей, оцененных с помощью GridSearchCV. В этом случае мы будем многократно запускать наш статистический тест, что приводит нас к проблеме множественного сравнения.

Существует много возможных способов решения этой проблемы, но стандартным подходом является применение поправки Бонферрони. Поправка Бонферрони может быть вычислена путем умножения p-значения на количество проводимых сравнений.

Давайте сравним производительность моделей, используя скорректированный t-тест:

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

Мы наблюдаем, что после поправки на множественные сравнения единственной моделью, которая значительно отличается от других, является '2_poly'. 'rbf', модель, занявшая первое место по результатам GridSearchCV, не имеет значимых различий с 'linear' или '3_poly'.

Парное сравнение всех моделей: байесовский подход

При использовании байесовской оценки для сравнения нескольких моделей нам не нужно вносить поправки на множественные сравнения (причины, по которым см. [4]).

Мы можем провести наши парные сравнения так же, как и в первой части:

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

Используя байесовский подход, мы можем вычислить вероятность того, что одна модель имеет лучшую, худшую или практически эквивалентную производительность по сравнению с другой.

Результаты показывают, что модель, занявшая первое место по результатам GridSearchCV 'rbf', имеет примерно 6,8% шанс быть хуже, чем 'linear', и 1,8% шанс быть хуже, чем '3_poly'. 'rbf' и 'linear' имеют 43% вероятность быть практически эквивалентными, а 'rbf' и '3_poly' имеют 10% шанс на это.

Аналогично выводам, полученным с помощью частотного подхода, все модели имеют 100% вероятность быть лучше, чем '2_poly', и ни одна из них не имеет практически эквивалентной производительности с последней.

Основные выводы

• Небольшие различия в показателях производительности могут быть случайными, а не потому, что одна модель систематически предсказывает лучше, чем другая. Как показано в этом примере, статистика может показать, насколько это вероятно.
• При статистическом сравнении производительности двух моделей, оцененных с помощью GridSearchCV, необходимо скорректировать рассчитанную дисперсию, которая может быть занижена, так как оценки моделей не являются независимыми друг от друга.
• Частотный подход, использующий (исправленный на дисперсию) парный t-тест, может показать, лучше ли производительность одной модели, чем другой, с уверенностью, превышающей случайность.
• Баесовский подход может предоставить вероятности того, что одна модель лучше, хуже или практически эквивалентна другой. Он также может показать, насколько мы уверены в том, что истинные различия в наших моделях попадают в определенный диапазон значений.
• Если сравнивается несколько моделей, при использовании частотного подхода требуется корректировка по множественному сравниванию.

Ссылки

Общее время выполнения скрипта: (0 минут 2.032 секунды)

Download Jupyter notebook: plot_grid_search_stats.ipynb

Download Python source code: plot_grid_search_stats.py

Download zipped: plot_grid_search_stats.zip

Связанные примеры

Тестирование с перестановками значимости классификационной оценки

Сравнение калибровки классификаторов

Общие ловушки при интерпретации коэффициентов линейных моделей

Пример конвейера для извлечения и оценки признаков текста