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

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

Для начала можно задаться вопросом: «Является ли первая модель статистически значимо лучше второй (при ранжировании по 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(next(iter(cv.split(X, y)))[0])
n_test = len(next(iter(cv.split(X, y)))[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

При общепринятом уровне значимости alpha, равном 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] говорит о том, что вероятность того, что истинная (средняя) разность в качестве моделей находится между x и y, равна 50%.

Определим доверительные интервалы для наших данных с уровнями 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

	нижняя граница	верхняя граница
интервал
0.50	0.000977	0.019023
0.75	-0.005422	0.025422
0.95	-0.016445	0.036445

Как показано в таблице, вероятность того, что истинная средняя разность между моделями находится между 0.000977 и 0.019023, равна 50%; вероятность того, что она находится между -0.005422 и 0.025422, равна 70%, а вероятность того, что она находится между -0.016445 и 0.036445, равна 95%.

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

Нас также может интересовать сравнение качества всех наших моделей, оценённых с помощью 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', вероятность оказаться хуже 'linear' составляет примерно 6.8%, а вероятность оказаться хуже '3_poly' — 1.8%. Вероятность практической эквивалентности 'rbf' и 'linear' составляет 43%, тогда как для 'rbf' и '3_poly' она равна 10%.

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

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

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

Литература

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

Download Jupyter notebook: plot_grid_search_stats.ipynb

Download Python source code: plot_grid_search_stats.py

Download zipped: plot_grid_search_stats.zip

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

Проверка значимости оценки классификации с помощью перестановок

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

Распространённые ошибки при интерпретации коэффициентов линейных моделей

Параметры SVM с ядром RBF