Регрессия Пуассона и ненормальное потерь

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

Загрузим набор данных о требованиях по страховым случаям из 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.10070308464041304
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, 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-регуляризацией, более известной как регрессия Риджа. Мы используем слабую регуляризацию 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) для имитации регрессора Риджа, чьи штрафные члены L2 масштабируются по-разному в зависимости от количества выборок.

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

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.566
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): 12196.1
Predicted number of claims by DummyRegressor(): 11931.2

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

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

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

Сумма всех предсказаний также подтверждает проблему калибровки модели Ridge: она недооценивает общее количество претензий в тестовой выборке более чем на 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 0x7656aa17e940>

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

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

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

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

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

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

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

plt.show()

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

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

Пример обычного метода наименьших квадратов