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

Загрузка наборов данных, базовое выделение признаков и определение целевых переменных

Мы создаём набор данных freMTPL2, объединяя таблицу freMTPL2freq, содержащую число страховых требований (ClaimNb), с таблицей freMTPL2sev, содержащей суммы страховых требований (ClaimAmount), для одинаковых идентификаторов полисов (IDpol).

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

df = load_mtpl2()


# Correct for unreasonable observations (that might be data error)
# and a few exceptionally large claim amounts
df["ClaimNb"] = df["ClaimNb"].clip(upper=4)
df["Exposure"] = df["Exposure"].clip(upper=1)
df["ClaimAmount"] = df["ClaimAmount"].clip(upper=200000)
# If the claim amount is 0, then we do not count it as a claim. The loss function
# used by the severity model needs strictly positive claim amounts. This way
# frequency and severity are more consistent with each other.
df.loc[(df["ClaimAmount"] == 0) & (df["ClaimNb"] >= 1), "ClaimNb"] = 0

log_scale_transformer = make_pipeline(
    FunctionTransformer(func=np.log), StandardScaler()
)

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

# Insurances companies are interested in modeling the Pure Premium, that is
# the expected total claim amount per unit of exposure for each policyholder
# in their portfolio:
df["PurePremium"] = df["ClaimAmount"] / df["Exposure"]

# This can be indirectly approximated by a 2-step modeling: the product of the
# Frequency times the average claim amount per claim:
df["Frequency"] = df["ClaimNb"] / df["Exposure"]
df["AvgClaimAmount"] = df["ClaimAmount"] / np.fmax(df["ClaimNb"], 1)

with pd.option_context("display.max_columns", 15):
    print(df[df.ClaimAmount > 0].head())

       ClaimNb  Exposure Area  VehPower  VehAge  DrivAge  BonusMalus VehBrand  \
IDpol
139          1      0.75    F         7       1       61          50      B12
190          1      0.14    B        12       5       50          60      B12
414          1      0.14    E         4       0       36          85      B12
424          2      0.62    F        10       0       51         100      B12
463          1      0.31    A         5       0       45          50      B12

          VehGas  Density Region  ClaimAmount   PurePremium  Frequency  \
IDpol
139    'Regular'    27000    R11       303.00    404.000000   1.333333
190     'Diesel'       56    R25      1981.84  14156.000000   7.142857
414    'Regular'     4792    R11      1456.55  10403.928571   7.142857
424    'Regular'    27000    R11     10834.00  17474.193548   3.225806
463    'Regular'       12    R73      3986.67  12860.225806   3.225806

       AvgClaimAmount
IDpol
139            303.00
190           1981.84
414           1456.55
424           5417.00
463           3986.67

Модель частоты — распределение Пуассона

Число страховых требований (ClaimNb) — это положительное целое число (включая 0). Поэтому эту целевую переменную можно моделировать распределением Пуассона. Предполагается, что оно описывает число дискретных событий, происходящих с постоянной интенсивностью в заданном временном интервале (Exposure, в годах). Здесь мы моделируем частоту y = ClaimNb / Exposure, которая по-прежнему имеет (масштабированное) распределение Пуассона, и используем Exposure в качестве sample_weight.

from sklearn.linear_model import PoissonRegressor
from sklearn.model_selection import train_test_split

df_train, df_test, X_train, X_test = train_test_split(df, X, random_state=0)

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

len(df_test)

169504

len(df_test[df_test["ClaimAmount"] > 0])

6237

В результате мы ожидаем значительную изменчивость результатов оценки при случайном повторном разбиении на обучающую и тестовую выборки.

Параметры модели оцениваются минимизацией девианса Пуассона на обучающей выборке с помощью метода Ньютона. Некоторые признаки коллинеарны (например, поскольку мы не удалили ни одного категориального уровня в OneHotEncoder), поэтому для устранения численных проблем используется слабая L2-регуляризация.

glm_freq = PoissonRegressor(alpha=1e-4, solver="newton-cholesky")
glm_freq.fit(X_train, df_train["Frequency"], sample_weight=df_train["Exposure"])

scores = score_estimator(
    glm_freq,
    X_train,
    X_test,
    df_train,
    df_test,
    target="Frequency",
    weights="Exposure",
)
print("Evaluation of PoissonRegressor on target Frequency")
print(scores)

Evaluation of PoissonRegressor on target Frequency
subset               train    test
metric
D² explained        0.0448  0.0427
mean abs. error     0.1379  0.1378
mean squared error  0.2441  0.2246

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

Можно визуально сравнить наблюдаемые и предсказанные значения, сгруппированные по возрасту водителей (DrivAge), возрасту автомобилей (VehAge) и страховому бонусу/малусу (BonusMalus).

fig, ax = plt.subplots(ncols=2, nrows=2, figsize=(16, 8))
fig.subplots_adjust(hspace=0.3, wspace=0.2)

plot_obs_pred(
    df=df_train,
    feature="DrivAge",
    weight="Exposure",
    observed="Frequency",
    predicted=glm_freq.predict(X_train),
    y_label="Claim Frequency",
    title="train data",
    ax=ax[0, 0],
)

plot_obs_pred(
    df=df_test,
    feature="DrivAge",
    weight="Exposure",
    observed="Frequency",
    predicted=glm_freq.predict(X_test),
    y_label="Claim Frequency",
    title="test data",
    ax=ax[0, 1],
    fill_legend=True,
)

plot_obs_pred(
    df=df_test,
    feature="VehAge",
    weight="Exposure",
    observed="Frequency",
    predicted=glm_freq.predict(X_test),
    y_label="Claim Frequency",
    title="test data",
    ax=ax[1, 0],
    fill_legend=True,
)

plot_obs_pred(
    df=df_test,
    feature="BonusMalus",
    weight="Exposure",
    observed="Frequency",
    predicted=glm_freq.predict(X_test),
    y_label="Claim Frequency",
    title="test data",
    ax=ax[1, 1],
    fill_legend=True,
)

Согласно наблюдаемым данным, частота аварий выше у водителей младше 30 лет и положительно коррелирует с переменной BonusMalus. Наша модель в целом корректно описывает это поведение.

Модель тяжести убытка — распределение Гамма

Эмпирически можно показать, что средняя сумма страхового требования, или тяжесть убытка (AvgClaimAmount), приблизительно подчиняется распределению Гамма. Мы обучаем GLM-модель тяжести убытка с теми же признаками, что и модель частоты.

Примечание:

• Мы отфильтровываем ClaimAmount == 0, поскольку область значений распределения Гамма — \((0, \infty)\), а не \([0, \infty)\).
• Мы используем ClaimNb в качестве sample_weight, чтобы учесть полисы, по которым было более одного страхового требования.

from sklearn.linear_model import GammaRegressor

mask_train = df_train["ClaimAmount"] > 0
mask_test = df_test["ClaimAmount"] > 0

glm_sev = GammaRegressor(alpha=10.0, solver="newton-cholesky")

glm_sev.fit(
    X_train[mask_train.values],
    df_train.loc[mask_train, "AvgClaimAmount"],
    sample_weight=df_train.loc[mask_train, "ClaimNb"],
)

scores = score_estimator(
    glm_sev,
    X_train[mask_train.values],
    X_test[mask_test.values],
    df_train[mask_train],
    df_test[mask_test],
    target="AvgClaimAmount",
    weights="ClaimNb",
)
print("Evaluation of GammaRegressor on target AvgClaimAmount")
print(scores)

Evaluation of GammaRegressor on target AvgClaimAmount
subset                     train          test
metric
D² explained        3.900000e-03  4.400000e-03
mean abs. error     1.756752e+03  1.744055e+03
mean squared error  5.801775e+07  5.030676e+07

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

from sklearn.dummy import DummyRegressor

dummy_sev = DummyRegressor(strategy="mean")
dummy_sev.fit(
    X_train[mask_train.values],
    df_train.loc[mask_train, "AvgClaimAmount"],
    sample_weight=df_train.loc[mask_train, "ClaimNb"],
)

scores = score_estimator(
    dummy_sev,
    X_train[mask_train.values],
    X_test[mask_test.values],
    df_train[mask_train],
    df_test[mask_test],
    target="AvgClaimAmount",
    weights="ClaimNb",
)
print("Evaluation of a mean predictor on target AvgClaimAmount")
print(scores)

Evaluation of a mean predictor on target AvgClaimAmount
subset                     train          test
metric
D² explained        0.000000e+00 -0.000000e+00
mean abs. error     1.756687e+03  1.744497e+03
mean squared error  5.803882e+07  5.033764e+07

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

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

print(
    "Mean AvgClaim Amount per policy:              %.2f "
    % df_train["AvgClaimAmount"].mean()
)
print(
    "Mean AvgClaim Amount | NbClaim > 0:           %.2f"
    % df_train["AvgClaimAmount"][df_train["AvgClaimAmount"] > 0].mean()
)
print(
    "Predicted Mean AvgClaim Amount | NbClaim > 0: %.2f"
    % glm_sev.predict(X_train).mean()
)
print(
    "Predicted Mean AvgClaim Amount (dummy) | NbClaim > 0: %.2f"
    % dummy_sev.predict(X_train).mean()
)

Mean AvgClaim Amount per policy:              71.78
Mean AvgClaim Amount | NbClaim > 0:           1951.21
Predicted Mean AvgClaim Amount | NbClaim > 0: 1940.95
Predicted Mean AvgClaim Amount (dummy) | NbClaim > 0: 1978.59

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

fig, ax = plt.subplots(ncols=1, nrows=2, figsize=(16, 6))

plot_obs_pred(
    df=df_train.loc[mask_train],
    feature="DrivAge",
    weight="Exposure",
    observed="AvgClaimAmount",
    predicted=glm_sev.predict(X_train[mask_train.values]),
    y_label="Average Claim Severity",
    title="train data",
    ax=ax[0],
)

plot_obs_pred(
    df=df_test.loc[mask_test],
    feature="DrivAge",
    weight="Exposure",
    observed="AvgClaimAmount",
    predicted=glm_sev.predict(X_test[mask_test.values]),
    y_label="Average Claim Severity",
    title="test data",
    ax=ax[1],
    fill_legend=True,
)
plt.tight_layout()

В целом возраст водителей (DrivAge) слабо влияет на тяжесть убытка как в наблюдаемых, так и в предсказанных данных.

Моделирование чистой страховой премии: произведение моделей или отдельный TweedieRegressor

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

В качестве альтернативы можно напрямую моделировать общие убытки с помощью обобщённой линейной модели Compound Poisson Gamma (с логарифмической функцией связи). Эта модель является частным случаем GLM-модели Твиди с параметром «степени» \(p \in (1, 2)\). Здесь мы заранее фиксируем параметр power модели Твиди на произвольном значении (1.9) из допустимого диапазона. В идеале это значение следует выбирать с помощью поиска по сетке, минимизируя отрицательное логарифмическое правдоподобие модели Твиди, но, к сожалению, текущая реализация пока этого не поддерживает.

Мы сравним эффективность обоих подходов. Чтобы оценить эффективность обеих моделей, можно вычислить средний девианс для обучающих и тестовых данных, предполагая распределение Compound Poisson-Gamma для общей суммы страховых требований. Это эквивалентно распределению Твиди с параметром power между 1 и 2.

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

from sklearn.linear_model import TweedieRegressor

glm_pure_premium = TweedieRegressor(power=1.9, alpha=0.1, solver="newton-cholesky")
glm_pure_premium.fit(
    X_train, df_train["PurePremium"], sample_weight=df_train["Exposure"]
)

tweedie_powers = [1.5, 1.7, 1.8, 1.9, 1.99, 1.999, 1.9999]

scores_product_model = score_estimator(
    (glm_freq, glm_sev),
    X_train,
    X_test,
    df_train,
    df_test,
    target="PurePremium",
    weights="Exposure",
    tweedie_powers=tweedie_powers,
)

scores_glm_pure_premium = score_estimator(
    glm_pure_premium,
    X_train,
    X_test,
    df_train,
    df_test,
    target="PurePremium",
    weights="Exposure",
    tweedie_powers=tweedie_powers,
)

scores = pd.concat(
    [scores_product_model, scores_glm_pure_premium],
    axis=1,
    sort=True,
    keys=("Product Model", "TweedieRegressor"),
)
print("Evaluation of the Product Model and the Tweedie Regressor on target PurePremium")
with pd.option_context("display.expand_frame_repr", False):
    print(scores)

Evaluation of the Product Model and the Tweedie Regressor on target PurePremium
                          Product Model               TweedieRegressor
subset                            train          test            train          test
metric
D² explained                        NaN           NaN     1.640000e-02  1.370000e-02
mean Tweedie dev p=1.5000  7.670090e+01  7.616900e+01     7.640460e+01  7.641980e+01
mean Tweedie dev p=1.7000  3.695810e+01  3.683920e+01     3.682720e+01  3.692730e+01
mean Tweedie dev p=1.8000  3.046050e+01  3.040490e+01     3.037490e+01  3.045690e+01
mean Tweedie dev p=1.9000  3.387610e+01  3.384970e+01     3.382040e+01  3.388020e+01
mean Tweedie dev p=1.9900  2.015718e+02  2.015412e+02     2.015342e+02  2.015600e+02
mean Tweedie dev p=1.9990  1.914574e+03  1.914370e+03     1.914537e+03  1.914388e+03
mean Tweedie dev p=1.9999  1.904751e+04  1.904556e+04     1.904747e+04  1.904558e+04
mean abs. error            2.730129e+02  2.722124e+02     2.740176e+02  2.731633e+02
mean squared error         3.295040e+07  3.212197e+07     3.295518e+07  3.213087e+07

В этом примере оба подхода к моделированию дают сопоставимые метрики эффективности. По причинам, связанным с реализацией, доля объяснённой дисперсии \(D^2\) для модели-произведения недоступна.

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

res = []
for subset_label, X, df in [
    ("train", X_train, df_train),
    ("test", X_test, df_test),
]:
    exposure = df["Exposure"].values
    res.append(
        {
            "subset": subset_label,
            "observed": df["ClaimAmount"].values.sum(),
            "predicted, frequency*severity model": np.sum(
                exposure * glm_freq.predict(X) * glm_sev.predict(X)
            ),
            "predicted, tweedie, power=%.2f" % glm_pure_premium.power: np.sum(
                exposure * glm_pure_premium.predict(X)
            ),
        }
    )

print(pd.DataFrame(res).set_index("subset").T)

subset                                      train          test
observed                             3.917618e+07  1.299546e+07
predicted, frequency*severity model  3.916579e+07  1.313280e+07
predicted, tweedie, power=1.90       3.952914e+07  1.325666e+07

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

Коэффициент Джини (основанный на площади между кривой и диагональю) можно использовать как метрику выбора модели для оценки её способности ранжировать страхователей. Обратите внимание, что эта метрика не отражает способность модели точно предсказывать абсолютную величину общей суммы страховых требований, а оценивает только относительные величины в качестве метрики ранжирования. Коэффициент Джини ограничен сверху значением 1.0, однако даже идеальная модель, ранжирующая страхователей по наблюдаемым суммам требований, не может набрать 1.0.

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

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

Наконец, важно отметить, что модель Compound Poisson Gamma, напрямую обученная для чистой страховой премии, проще в разработке и сопровождении, поскольку состоит из одного оценщика scikit-learn, а не пары моделей, каждая из которых имеет собственный набор гиперпараметров.

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_exposure = exposure[ranking]
    ranked_pure_premium = y_true[ranking]
    cumulative_claim_amount = np.cumsum(ranked_pure_premium * ranked_exposure)
    cumulative_claim_amount /= cumulative_claim_amount[-1]
    cumulative_exposure = np.cumsum(ranked_exposure)
    cumulative_exposure /= cumulative_exposure[-1]
    return cumulative_exposure, cumulative_claim_amount


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

y_pred_product = glm_freq.predict(X_test) * glm_sev.predict(X_test)
y_pred_total = glm_pure_premium.predict(X_test)

for label, y_pred in [
    ("Frequency * Severity model", y_pred_product),
    ("Compound Poisson Gamma", y_pred_total),
]:
    cum_exposure, cum_claims = lorenz_curve(
        df_test["PurePremium"], y_pred, df_test["Exposure"]
    )
    gini = 1 - 2 * auc(cum_exposure, cum_claims)
    label += " (Gini index: {:.3f})".format(gini)
    ax.plot(cum_exposure, cum_claims, linestyle="-", label=label)

# Oracle model: y_pred == y_test
cum_exposure, cum_claims = lorenz_curve(
    df_test["PurePremium"], df_test["PurePremium"], df_test["Exposure"]
)
gini = 1 - 2 * auc(cum_exposure, cum_claims)
label = "Oracle (Gini index: {:.3f})".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",
    xlabel=(
        "Cumulative proportion of exposure\n(ordered by model from safest to riskiest)"
    ),
    ylabel="Cumulative proportion of claim amounts",
)
ax.legend(loc="upper left")
plt.plot()

[]

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

Download Jupyter notebook: plot_tweedie_regression_insurance_claims.ipynb

Download Python source code: plot_tweedie_regression_insurance_claims.py

Download zipped: plot_tweedie_regression_insurance_claims.zip

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

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

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

Преобразование данных к нормальному распределению

Визуализация прогнозов с кросс-валидацией