Сравнение метрик двух моделей

Представим, что мы хотим выбрать лучшую регрессионную модель для неких данных. Для этого мы обучим, например, катбуст и линейную регрессию. А затем проведем кросс-валидацию на 5 фолдах и получим следующие значения mse:

На первый взгляд может показаться, что у линейной регрессии лучшие показатели. Так ли это на самом деле? Если мы вспомним наш предыдущий пост, то нам станет понятно: пока рановато говорить о том, что catboost хуже, ведь мы провалидировали нашу модель не на генеральной совокупности, а лишь на подвыборках. Следовательно, нам нужно привлечь статистику!

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

Возьмем критическое значение \( α=0.05 \), с которым будем сравнивать \( P_{value} \) для опровержения/подтверждения гипотезы. А выбор статистического теста будет зависеть от того, является ли распределение метрик нормальным. Кстати, в следующих статьях мы обязательно расскажем вам о возможностях проверки нормального распределения метрик 🙂 При нормальном распределении метрик мы можем использовать параметрические тесты. Они подразумевают, что выборка была получена из распределения с конечным числом параметров. Например, из \( N(\alpha, \sigma) \). Отсюда и их название. Большинство параметрических тестов ожидают нормально распределенную выборку, соответственно, и все дальнейшие параметрические тесты (о которых будет сказано далее) для нормальных распределений. Параметрические тесты с большей вероятностью обнаружат различия, если они есть, поэтому лучше использовать их, когда это возможно.

Т-тесты

Итак, давайте считать распределение нормальным. Тогда мы можем использовать параметрический Т-тест Стюдента.

1. Вычислим попарные разницы между моделями (строка diff в табличке выше).
2. Вычислим среднюю разницу: \( M=-0,01 \).
3. Вычислим стандартное отклонение разниц: \( sd=0,007 \).
4. Вычислим стандартную ошибку: \( SE =\frac{sd}{sqrt(n)} \). \( SE=0,003 \).
5. Посчитаем T-статистику: \( t=\frac{M}{SE}=-3.33 \).
6. Посчитаем число степеней свободы: \( df=n-1 = 4 \).

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

7. Воспользуемся калькулятором T-распределения и получим \( Pvalue = 0.029024 \).
8. Сравним Pvalue и α и поймем: мы можем опровергнуть нашу нулевую гипотезу и сказать, что модели предсказывают не одинаково, так как \( P_{value}<0.05 \). Учитывая, что мы опровергли нулевую гипотезу и то, что в среднем у линейной регрессии mse ниже, можно сделать вывод, что она действительно лучше! В этом и дальнейших примерах нам повезло, разброс между метриками на фолдах был не сильно большой, но обычно это не так. В реальности лучше делать несколько разбиений на фолды (как в этой статье) либо использовать 5x2cv.

Стоит заметить, что T-тестом Стьюдента не стоит пользоваться, если есть разница в количестве наблюдений в выборке. В таких случаях лучше подходит T-тест Уэлча.

Тест Стьюдента в питоне

Теперь давайте посмотрим, как эти тесты можно выполнить в python.

Сначала возьмем реальный датасет предсказания степени заболевания диабетом, создадим два регрессора — катбуст и случайный лес, запустим кроссвалидацию.

from sklearn.datasets import load_diabetes
from sklearn.ensemble import RandomForestRegressor
from catboost import CatBoostRegressor
from sklearn.model_selection import cross_val_score, KFold

X, y = load_diabetes(return_X_y=True)
rf = RandomForestRegressor()
cb = CatBoostRegressor(verbose=False)
folds = KFold(n_splits=5)
scores_rf = cross_val_score(rf, X, y, scoring='neg_mean_squared_error', cv=folds)
scores_cb = cross_val_score(cb, X, y, scoring='neg_mean_squared_error', cv=folds)

Далее посчитаем значение T-статистики с помощью библиотеки stats.

import scipy.stats as stats

stats.ttest_rel(scores_rf, scores_cb)
# output: Ttest_relResult(statistic=1.9353919550011929, pvalue=0.12503778196479817)

Так как \( P_{value}>0.05 \) мы не можем опровергнуть нашу нулевую гипотезу.

Тест Уэлча в питоне

Давайте предположим, что мы по проводили кроссвалидацию с разным числом фолдов:

scores_rf = cross_val_score(rf, X, y, cv=5, scoring='neg_mean_squared_error')
scores_cb = cross_val_score(cb, X, y, cv=10, scoring='neg_mean_squared_error')

В таком случае нам нужно использовать тест Уэлча, потому что мы уже не можем вычислить попарные разницы.

stats.ttest_ind(scores_rf, scores_cb, equal_var=False)
# output: Ttest_indResult(statistic=0.8564592491178395, pvalue=0.40778645901968114)

Так как \( P_{value}>0.05 \), мы не можем отклонить нулевую гипотезу.

Критерий Вилкоксона

Теперь давайте представим, что наши разницы распределены ненормально.

В таком случае нам поможет критерий Вилкоксона. Его еще называют критерием знаковых рангов. Сейчас станет понятно почему 🙂

Важное замечание: в данном тесте наша нулевая гипотеза заключается в том, что медианы метрик моделей одинаковые!

1. Вычислим попарные разницы между моделями (строка diff в табличке выше).
2. Проранжируем в порядке возрастания их модулей (строка rank abs). Можно заметить, что для повторяющихся значений мы выбрали средний ранг, а ноль пропустили. В нашем примере три повторяющихся значения -0.01, они получат значения рангов 1, 2, 3 — среднее равно 2, поэтому для всех значения -0.01 в таблице установили ранг 2.
3. Посчитаем сумму рангов, соответствующих разнице положительной \( W+=17,5 \) и отрицательной \( W-=10,5 \).
4. Сравним \( W=min(W-,W+)=10,5 \) c табличным значением = 3. Это значит, что мы мы не можем опровергнуть нулевую гипотезу.

Критерий Вилкоксона в питоне

Теперь давайте попробуем сделать тоже самое на Python.

Снова возьмем датасет про диабет, создадим два регрессора — катбуст и случайный лес, запустим кроссвалидацию.

from sklearn.datasets import load_diabetes
from sklearn.ensemble import RandomForestRegressor
from catboost import CatBoostRegressor
from sklearn.model_selection import cross_val_score, KFold

X, y = load_diabetes(return_X_y=True)
rf = RandomForestRegressor()
cb = CatBoostRegressor(verbose=False)
folds = KFold(n_splits=5)
scores_rf = cross_val_score(rf, X, y, scoring='neg_mean_squared_error', cv=folds)
scores_cb = cross_val_score(cb, X, y, scoring='neg_mean_squared_error', cv=folds)

Опять же обратимся к stats.

from scipy.stats import wilcoxon

wilcoxon(scores_rf, scores_cb)
# output: WilcoxonResult(statistic=1.0, pvalue=0.125)

И снова не сможем отклонить нашу нулевую гипотезу, так как \( P_{value}>0.05 \).