Регрессия Tweedie на страховых претензиях

Загрузка наборов данных, базовая обработка функций и определения целевых значений

Мы создаем набор данных 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, 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.756746e+03  1.744042e+03
mean squared error  5.801770e+07  5.030677e+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

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

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

Мы сравним производительность обоих подходов. Для количественной оценки производительности обеих моделей можно вычислить среднюю девиацию обучающей и тестовой выборки, предполагая распределение «Пуассон-Гамма» для общей суммы претензий. Это эквивалентно распределению Tweedie с параметром 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.669930e+01  7.617050e+01     7.640770e+01  7.640880e+01
mean Tweedie dev p=1.7000  3.695740e+01  3.683980e+01     3.682880e+01  3.692270e+01
mean Tweedie dev p=1.8000  3.046010e+01  3.040530e+01     3.037600e+01  3.045390e+01
mean Tweedie dev p=1.9000  3.387580e+01  3.385000e+01     3.382120e+01  3.387830e+01
mean Tweedie dev p=1.9900  2.015716e+02  2.015414e+02     2.015347e+02  2.015587e+02
mean Tweedie dev p=1.9990  1.914573e+03  1.914370e+03     1.914538e+03  1.914387e+03
mean Tweedie dev p=1.9999  1.904751e+04  1.904556e+04     1.904747e+04  1.904558e+04
mean abs. error            2.730119e+02  2.722128e+02     2.739865e+02  2.731249e+02
mean squared error         3.295040e+07  3.212197e+07     3.295505e+07  3.213056e+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.916555e+07  1.313276e+07
predicted, tweedie, power=1.90       3.951751e+07  1.325198e+07

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

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

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

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

Наконец, следует отметить, что модель Пуассона-Гамма, которая напрямую подгоняется под чистый тариф, проще в разработке и обслуживании, так как она состоит из одного оценщика 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 минут 9,391 секунды)

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

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

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

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

Отображение данных в нормальное распределение

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