Регрессия Пуассона и функция потерь, отличная от нормальной

Набор данных по страховым требованиям по автострахованию гражданской ответственности во Франции

Загрузим набор данных по автостраховым требованиям из OpenML: https://www.openml.org/d/41214

from sklearn.datasets import fetch_openml

df = fetch_openml(data_id=41214, as_frame=True).frame
df

Количество страховых требований (ClaimNb) — положительное целое число, которое можно моделировать распределением Пуассона. Предполагается, что это число дискретных событий, происходящих с постоянной интенсивностью за заданный интервал времени (Exposure, в годах).

Здесь мы хотим моделировать частоту y = ClaimNb / Exposure условно на X с помощью (масштабированного) распределения Пуассона и использовать Exposure как sample_weight.

df["Frequency"] = df["ClaimNb"] / df["Exposure"]

print(
    "Average Frequency = {}".format(np.average(df["Frequency"], weights=df["Exposure"]))
)

print(
    "Fraction of exposure with zero claims = {0:.1%}".format(
        df.loc[df["ClaimNb"] == 0, "Exposure"].sum() / df["Exposure"].sum()
    )
)

fig, (ax0, ax1, ax2) = plt.subplots(ncols=3, figsize=(16, 4))
ax0.set_title("Number of claims")
_ = df["ClaimNb"].hist(bins=30, log=True, ax=ax0)
ax1.set_title("Exposure in years")
_ = df["Exposure"].hist(bins=30, log=True, ax=ax1)
ax2.set_title("Frequency (number of claims per year)")
_ = df["Frequency"].hist(bins=30, log=True, ax=ax2)

Average Frequency = 0.10070308464041305
Fraction of exposure with zero claims = 93.9%

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

Поэтому для обучения линейных моделей с такими предикторами необходимо выполнить стандартные преобразования признаков:

from sklearn.compose import ColumnTransformer
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import (
    FunctionTransformer,
    KBinsDiscretizer,
    OneHotEncoder,
    StandardScaler,
)

log_scale_transformer = make_pipeline(
    FunctionTransformer(np.log, validate=False), StandardScaler()
)

linear_model_preprocessor = ColumnTransformer(
    [
        ("passthrough_numeric", "passthrough", ["BonusMalus"]),
        (
            "binned_numeric",
            KBinsDiscretizer(
                n_bins=10, quantile_method="averaged_inverted_cdf", random_state=0
            ),
            ["VehAge", "DrivAge"],
        ),
        ("log_scaled_numeric", log_scale_transformer, ["Density"]),
        (
            "onehot_categorical",
            OneHotEncoder(),
            ["VehBrand", "VehPower", "VehGas", "Region", "Area"],
        ),
    ],
    remainder="drop",
)

Базовый уровень с постоянным прогнозом

Стоит отметить, что у более чем 93% страхователей нет страховых требований. Если преобразовать эту задачу в задачу бинарной классификации, классы будут сильно несбалансированы, и даже простая модель, которая всегда предсказывает среднее значение, сможет достичь точности 93%.

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

from sklearn.dummy import DummyRegressor
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline

df_train, df_test = train_test_split(df, test_size=0.33, random_state=0)

dummy = Pipeline(
    [
        ("preprocessor", linear_model_preprocessor),
        ("regressor", DummyRegressor(strategy="mean")),
    ]
).fit(df_train, df_train["Frequency"], regressor__sample_weight=df_train["Exposure"])

Вычислим качество этого базового прогноза с постоянным значением, используя 3 различные метрики регрессии:

from sklearn.metrics import (
    mean_absolute_error,
    mean_poisson_deviance,
    mean_squared_error,
)


def score_estimator(estimator, df_test):
    """Score an estimator on the test set."""
    y_pred = estimator.predict(df_test)

    print(
        "MSE: %.3f"
        % mean_squared_error(
            df_test["Frequency"], y_pred, sample_weight=df_test["Exposure"]
        )
    )
    print(
        "MAE: %.3f"
        % mean_absolute_error(
            df_test["Frequency"], y_pred, sample_weight=df_test["Exposure"]
        )
    )

    # Ignore non-positive predictions, as they are invalid for
    # the Poisson deviance.
    mask = y_pred > 0
    if (~mask).any():
        n_masked, n_samples = (~mask).sum(), mask.shape[0]
        print(
            "WARNING: Estimator yields invalid, non-positive predictions "
            f" for {n_masked} samples out of {n_samples}. These predictions "
            "are ignored when computing the Poisson deviance."
        )

    print(
        "mean Poisson deviance: %.3f"
        % mean_poisson_deviance(
            df_test["Frequency"][mask],
            y_pred[mask],
            sample_weight=df_test["Exposure"][mask],
        )
    )


print("Constant mean frequency evaluation:")
score_estimator(dummy, df_test)

Constant mean frequency evaluation:
MSE: 0.564
MAE: 0.189
mean Poisson deviance: 0.625

(Обобщённые) линейные модели

Сначала смоделируем целевую переменную с помощью линейной регрессии наименьших квадратов (со штрафом l2), более известной как регрессия Ridge. Мы используем небольшой коэффициент регуляризации alpha, поскольку ожидаем, что такая линейная модель будет недообучаться на столь большом наборе данных.

from sklearn.linear_model import Ridge

ridge_glm = Pipeline(
    [
        ("preprocessor", linear_model_preprocessor),
        ("regressor", Ridge(alpha=1e-6)),
    ]
).fit(df_train, df_train["Frequency"], regressor__sample_weight=df_train["Exposure"])

Девианс Пуассона нельзя вычислить для неположительных предсказаний модели. Для моделей, которые действительно выдают несколько неположительных предсказаний (например, Ridge), мы игнорируем соответствующие образцы, то есть полученный девианс Пуассона является приближённым. Альтернативный подход — использовать метаоцениватель TransformedTargetRegressor, чтобы преобразовать y_pred в строго положительную область.

print("Ridge evaluation:")
score_estimator(ridge_glm, df_test)

Ridge evaluation:
MSE: 0.560
MAE: 0.186
WARNING: Estimator yields invalid, non-positive predictions  for 595 samples out of 223745. These predictions are ignored when computing the Poisson deviance.
mean Poisson deviance: 0.597

Затем обучим регрессор Пуассона на целевой переменной. Мы задаём силу регуляризации alpha примерно равной 1e-6, делённому на число образцов (то есть 1e-12), чтобы имитировать регрессор Ridge, член штрафа L2 которого масштабируется иначе в зависимости от числа образцов.

Поскольку регрессор Пуассона внутри модели использует логарифм ожидаемого значения целевой переменной, а не само ожидаемое значение (логарифмическую, а не тождественную функцию связи), зависимость между X и y уже не является строго линейной. Поэтому регрессор Пуассона называют обобщённой линейной моделью (GLM), а не обычной линейной моделью, как в случае регрессии Ridge.

from sklearn.linear_model import PoissonRegressor

n_samples = df_train.shape[0]

poisson_glm = Pipeline(
    [
        ("preprocessor", linear_model_preprocessor),
        ("regressor", PoissonRegressor(alpha=1e-12, solver="newton-cholesky")),
    ]
)
poisson_glm.fit(
    df_train, df_train["Frequency"], regressor__sample_weight=df_train["Exposure"]
)

print("PoissonRegressor evaluation:")
score_estimator(poisson_glm, df_test)

PoissonRegressor evaluation:
MSE: 0.560
MAE: 0.186
mean Poisson deviance: 0.594

Градиентный бустинг деревьев регрессии для регрессии Пуассона

Наконец, рассмотрим нелинейную модель — градиентный бустинг деревьев регрессии. Для моделей на основе деревьев не требуется выполнять one-hot-кодирование категориальных данных: вместо этого можно закодировать каждую метку категории произвольным целым числом с помощью OrdinalEncoder. При таком кодировании деревья будут считать категориальные признаки упорядоченными, что не всегда желательно. Однако для достаточно глубоких деревьев этот эффект ограничен: они способны восстановить категориальную природу признаков. Главное преимущество OrdinalEncoder перед OneHotEncoder заключается в ускорении обучения.

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

from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.preprocessing import OrdinalEncoder

tree_preprocessor = ColumnTransformer(
    [
        (
            "categorical",
            OrdinalEncoder(),
            ["VehBrand", "VehPower", "VehGas", "Region", "Area"],
        ),
        ("numeric", "passthrough", ["VehAge", "DrivAge", "BonusMalus", "Density"]),
    ],
    remainder="drop",
)
poisson_gbrt = Pipeline(
    [
        ("preprocessor", tree_preprocessor),
        (
            "regressor",
            HistGradientBoostingRegressor(loss="poisson", max_leaf_nodes=128),
        ),
    ]
)
poisson_gbrt.fit(
    df_train, df_train["Frequency"], regressor__sample_weight=df_train["Exposure"]
)

print("Poisson Gradient Boosted Trees evaluation:")
score_estimator(poisson_gbrt, df_test)

Poisson Gradient Boosted Trees evaluation:
MSE: 0.561
MAE: 0.184
mean Poisson deviance: 0.575

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

Оценка моделей на единственном разделении на обучающую и тестовую выборки подвержена случайным колебаниям. Если позволяют вычислительные ресурсы, следует проверить, приводят ли метрики качества при кросс-валидации к аналогичным выводам.

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

fig, axes = plt.subplots(nrows=2, ncols=4, figsize=(16, 6), sharey=True)
fig.subplots_adjust(bottom=0.2)
n_bins = 20
for row_idx, label, df in zip(range(2), ["train", "test"], [df_train, df_test]):
    df["Frequency"].hist(bins=np.linspace(-1, 30, n_bins), ax=axes[row_idx, 0])

    axes[row_idx, 0].set_title("Data")
    axes[row_idx, 0].set_yscale("log")
    axes[row_idx, 0].set_xlabel("y (observed Frequency)")
    axes[row_idx, 0].set_ylim([1e1, 5e5])
    axes[row_idx, 0].set_ylabel(label + " samples")

    for idx, model in enumerate([ridge_glm, poisson_glm, poisson_gbrt]):
        y_pred = model.predict(df)

        pd.Series(y_pred).hist(
            bins=np.linspace(-1, 4, n_bins), ax=axes[row_idx, idx + 1]
        )
        axes[row_idx, idx + 1].set(
            title=model[-1].__class__.__name__,
            yscale="log",
            xlabel="y_pred (predicted expected Frequency)",
        )
plt.tight_layout()

Экспериментальные данные для y имеют распределение с длинным хвостом. Во всех моделях мы предсказываем ожидаемую частоту случайной величины, поэтому экстремальных значений неизбежно будет меньше, чем среди наблюдаемых реализаций этой случайной величины. Это объясняет, почему мода гистограмм предсказаний модели не обязательно соответствует наименьшему значению. Кроме того, нормальное распределение, используемое в Ridge, имеет постоянную дисперсию, тогда как для распределения Пуассона, используемого в PoissonRegressor и HistGradientBoostingRegressor, дисперсия пропорциональна предсказанному ожидаемому значению.

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

Оцениватель HistGradientBoostingRegressor обладает наибольшей гибкостью и способен предсказывать более высокие ожидаемые значения.

Обратите внимание, что для модели HistGradientBoostingRegressor можно было бы использовать функцию потерь наименьших квадратов. Это ошибочно предполагало бы нормальное распределение целевой переменной, как и модель Ridge, и, возможно, приводило бы к слегка отрицательным прогнозам. Однако градиентный бустинг деревьев всё равно показал бы относительно хорошие результаты и, в частности, был бы лучше модели PoissonRegressor благодаря сочетанию гибкости деревьев с большим количеством обучающих образцов.

Оценка калибровки прогнозов

Чтобы убедиться, что оцениватели выдают разумные прогнозы для разных типов страхователей, можно разбить тестовые образцы на интервалы в соответствии с y_pred, возвращаемым каждой моделью. Затем для каждого интервала сравним среднее предсказанное значение y_pred со средним наблюдаемым значением целевой переменной:

from sklearn.utils import gen_even_slices


def _mean_frequency_by_risk_group(y_true, y_pred, sample_weight=None, n_bins=100):
    """Compare predictions and observations for bins ordered by y_pred.

    We order the samples by ``y_pred`` and split it in bins.
    In each bin the observed mean is compared with the predicted mean.

    Parameters
    ----------
    y_true: array-like of shape (n_samples,)
        Ground truth (correct) target values.
    y_pred: array-like of shape (n_samples,)
        Estimated target values.
    sample_weight : array-like of shape (n_samples,)
        Sample weights.
    n_bins: int
        Number of bins to use.

    Returns
    -------
    bin_centers: ndarray of shape (n_bins,)
        bin centers
    y_true_bin: ndarray of shape (n_bins,)
        average y_pred for each bin
    y_pred_bin: ndarray of shape (n_bins,)
        average y_pred for each bin
    """
    idx_sort = np.argsort(y_pred)
    bin_centers = np.arange(0, 1, 1 / n_bins) + 0.5 / n_bins
    y_pred_bin = np.zeros(n_bins)
    y_true_bin = np.zeros(n_bins)

    for n, sl in enumerate(gen_even_slices(len(y_true), n_bins)):
        weights = sample_weight[idx_sort][sl]
        y_pred_bin[n] = np.average(y_pred[idx_sort][sl], weights=weights)
        y_true_bin[n] = np.average(y_true[idx_sort][sl], weights=weights)
    return bin_centers, y_true_bin, y_pred_bin


print(f"Actual number of claims: {df_test['ClaimNb'].sum()}")
fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(12, 8))
plt.subplots_adjust(wspace=0.3)

for axi, model in zip(ax.ravel(), [ridge_glm, poisson_glm, poisson_gbrt, dummy]):
    y_pred = model.predict(df_test)
    y_true = df_test["Frequency"].values
    exposure = df_test["Exposure"].values
    q, y_true_seg, y_pred_seg = _mean_frequency_by_risk_group(
        y_true, y_pred, sample_weight=exposure, n_bins=10
    )

    # Name of the model after the estimator used in the last step of the
    # pipeline.
    print(f"Predicted number of claims by {model[-1]}: {np.sum(y_pred * exposure):.1f}")

    axi.plot(q, y_pred_seg, marker="x", linestyle="--", label="predictions")
    axi.plot(q, y_true_seg, marker="o", linestyle="--", label="observations")
    axi.set_xlim(0, 1.0)
    axi.set_ylim(0, 0.5)
    axi.set(
        title=model[-1],
        xlabel="Fraction of samples sorted by y_pred",
        ylabel="Mean Frequency (y_pred)",
    )
    axi.legend()
plt.tight_layout()

Actual number of claims: 11935
Predicted number of claims by Ridge(alpha=1e-06): 11933.4
Predicted number of claims by PoissonRegressor(alpha=1e-12, solver='newton-cholesky'): 11932.0
Predicted number of claims by HistGradientBoostingRegressor(loss='poisson', max_leaf_nodes=128): 12158.8
Predicted number of claims by DummyRegressor(): 11931.2

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

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

PoissonRegressor и HistGradientBoostingRegressor демонстрируют лучшее соответствие между предсказанными и наблюдаемыми целевыми значениями, особенно при низких предсказанных значениях.

Сумма всех прогнозов также подтверждает проблему калибровки модели Ridge: она недооценивает общее число страховых требований в тестовой выборке более чем на 3%, тогда как остальные три модели могут приблизительно восстановить общее число требований в тестовом портфеле.

Оценка способности ранжирования

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

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

Этот график называется кривой Лоренца, и его можно суммировать с помощью индекса Джини:

from sklearn.metrics import auc


def lorenz_curve(y_true, y_pred, exposure):
    y_true, y_pred = np.asarray(y_true), np.asarray(y_pred)
    exposure = np.asarray(exposure)

    # order samples by increasing predicted risk:
    ranking = np.argsort(y_pred)
    ranked_frequencies = y_true[ranking]
    ranked_exposure = exposure[ranking]
    cumulated_claims = np.cumsum(ranked_frequencies * ranked_exposure)
    cumulated_claims /= cumulated_claims[-1]
    cumulated_exposure = np.cumsum(ranked_exposure)
    cumulated_exposure /= cumulated_exposure[-1]
    return cumulated_exposure, cumulated_claims


fig, ax = plt.subplots(figsize=(8, 8))

for model in [dummy, ridge_glm, poisson_glm, poisson_gbrt]:
    y_pred = model.predict(df_test)
    cum_exposure, cum_claims = lorenz_curve(
        df_test["Frequency"], y_pred, df_test["Exposure"]
    )
    gini = 1 - 2 * auc(cum_exposure, cum_claims)
    label = "{} (Gini: {:.2f})".format(model[-1], gini)
    ax.plot(cum_exposure, cum_claims, linestyle="-", label=label)

# Oracle model: y_pred == y_test
cum_exposure, cum_claims = lorenz_curve(
    df_test["Frequency"], df_test["Frequency"], df_test["Exposure"]
)
gini = 1 - 2 * auc(cum_exposure, cum_claims)
label = "Oracle (Gini: {:.2f})".format(gini)
ax.plot(cum_exposure, cum_claims, linestyle="-.", color="gray", label=label)

# Random Baseline
ax.plot([0, 1], [0, 1], linestyle="--", color="black", label="Random baseline")
ax.set(
    title="Lorenz curves by model",
    xlabel="Cumulative proportion of exposure (from safest to riskiest)",
    ylabel="Cumulative proportion of claims",
)
ax.legend(loc="upper left")

<matplotlib.legend.Legend object at 0x784128f60440>

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

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

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

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

Линейные модели предполагают отсутствие взаимодействий между входными переменными, что, вероятно, приводит к недообучению. Добавление извлекателя полиномиальных признаков (PolynomialFeatures) действительно повышает их способность к дискриминации на 2 пункта индекса Джини. В частности, это улучшает способность моделей выявлять 5% профилей с наибольшим риском.

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

• Качество моделей можно оценивать по их способности выдавать хорошо откалиброванные прогнозы и обеспечивать хорошее ранжирование.
• Калибровку модели можно оценить, построив график среднего наблюдаемого значения относительно среднего предсказанного значения для групп тестовых образцов, разбитых на интервалы по предсказанному риску.
• Функция потерь наименьших квадратов (вместе с неявным использованием тождественной функции связи) в модели регрессии Ridge, по-видимому, приводит к плохой калибровке. В частности, модель склонна недооценивать риск и даже может предсказывать недопустимые отрицательные частоты.
• Использование функции потерь Пуассона с логарифмической функцией связи позволяет исправить эти проблемы и получить хорошо откалиброванную линейную модель.
• Индекс Джини отражает способность модели ранжировать прогнозы независимо от их абсолютных значений, а потому оценивает только способность к ранжированию.
• Несмотря на улучшение калибровки, способность обеих линейных моделей к ранжированию сопоставима и значительно уступает способности градиентного бустинга деревьев регрессии.
• Девианс Пуассона, используемый в качестве метрики оценки, отражает и калибровку, и способность модели к ранжированию. Он также предполагает линейную зависимость между ожидаемым значением и дисперсией целевой переменной. Для краткости мы не проверяли, выполняется ли это предположение.
• Традиционные метрики регрессии, такие как среднеквадратичная и средняя абсолютная ошибки, трудно содержательно интерпретировать для счётных данных с большим количеством нулей.

plt.show()

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

Download Jupyter notebook: plot_poisson_regression_non_normal_loss.ipynb

Download Python source code: plot_poisson_regression_non_normal_loss.py

Download zipped: plot_poisson_regression_non_normal_loss.zip

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

Регрессия Твиди для страховых требований

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

Основные нововведения scikit-learn 0.23

Объединение предикторов с помощью стекинга