{"id":1440,"date":"2024-02-02T19:15:49","date_gmt":"2024-02-02T18:15:49","guid":{"rendered":"https:\/\/lorentzen.ch\/?p=1440"},"modified":"2024-04-12T08:43:50","modified_gmt":"2024-04-12T06:43:50","slug":"ml-xai-strong-glm-in-python","status":"publish","type":"post","link":"https:\/\/lorentzen.ch\/index.php\/2024\/02\/02\/ml-xai-strong-glm-in-python\/","title":{"rendered":"Building Strong GLMs in Python via ML + XAI"},"content":{"rendered":"\n

In our latest post<\/a>, we explained how to use ML + XAI to build strong generalized linear models with R. Let’s do the same with Python.<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n

Insurance pricing data<\/h2>\n\n\n\n

We will use again a synthetic dataset with 1 Mio insurance policies, with reference:<\/p>\n\n\n\n

Mayer, M., Meier, D. and Wuthrich, M.V. (2023),
SHAP for Actuaries: Explain any Model.
https:\/\/<\/a>doi.org\/10.2139\/ssrn.4389797<\/a><\/p>\n\n\n\n

Let’s start by loading and describing the data:<\/p>\n\n\n\n

import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport shap\nfrom sklearn.datasets import fetch_openml\nfrom sklearn.inspection import PartialDependenceDisplay\nfrom sklearn.metrics import mean_poisson_deviance\nfrom sklearn.dummy import DummyRegressor\nfrom lightgbm import LGBMRegressor\n# We need preview version of glum that adds formulaic API\n# !pip install git+https:\/\/github.com\/Quantco\/glum@glum-v3#egg=glum\nfrom glum import GeneralizedLinearRegressor\n\n# Load data\n\ndf = fetch_openml(data_id=45106, parser="pandas").frame\ndf.head()\n\n# Continuous features\ndf.hist(["driver_age", "car_weight", "car_power", "car_age"])\n_ = plt.suptitle("Histograms of continuous features", fontsize=15)\n\n# Response and discrete features\nfig, axes = plt.subplots(figsize=(8, 3), ncols=3)\n\nfor v, ax in zip(["claim_nb", "year", "town"], axes):\n    df[v].value_counts(sort=False).sort_index().plot(kind="bar", ax=ax, rot=0, title=v)\nplt.suptitle("Barplots of response and discrete features", fontsize=15)\nplt.tight_layout()\nplt.show()\n\n# Rank correlations\ncorr = df.corr("spearman")\nmask = np.triu(np.ones_like(corr, dtype=bool))\nplt.suptitle("Rank-correlogram", fontsize=15)\n_ = sns.heatmap(\n    corr, mask=mask, vmin=-0.7, vmax=0.7, center=0, cmap="vlag", square=True\n)\n\n<\/pre><\/div>\n\n\n\n
\"\"<\/figure>\n\n\n\n
\"\"<\/figure>\n\n\n\n
\"\"<\/figure>\n\n\n\n
Modeling<\/h2>\n\n\n\n
    \n
  1. We fit a tuned Boosted Trees model to model log(E(claim count)) via Poisson deviance loss.<\/li>\n\n\n\n
  2. And perform a SHAP analysis to derive insights.<\/li>\n<\/ol>\n\n\n\n
    from sklearn.model_selection import train_test_split\n\nX_train, X_test, y_train, y_test = train_test_split(\n    df.drop("claim_nb", axis=1), df["claim_nb"], test_size=0.1, random_state=30\n)\n\n# Tuning step not shown. Number of boosting rounds found via early stopping on CV performance\nparams = dict(\n    learning_rate=0.05,\n    objective="poisson",\n    num_leaves=7,\n    min_child_samples=50,\n    min_child_weight=0.001,\n    colsample_bynode=0.8,\n    subsample=0.8,\n    reg_alpha=3,\n    reg_lambda=5,\n    verbose=-1,\n)\n\nmodel_lgb = LGBMRegressor(n_estimators=360, **params)\nmodel_lgb.fit(X_train, y_train)\n\n# SHAP analysis\nX_explain = X_train.sample(n=2000, random_state=937)\nexplainer = shap.Explainer(model_lgb)\nshap_val = explainer(X_explain)\n\nplt.suptitle("SHAP importance", fontsize=15)\nshap.plots.bar(shap_val)\n\nfor s in [shap_val[:, 0:3], shap_val[:, 3:]]:\n    shap.plots.scatter(s, color=shap_val, ymin=-0.5, ymax=1)<\/pre><\/div>\n\n\n\n
    \"\"<\/figure>\n\n\n\n
    \"\"<\/figure>\n\n\n\n
    Here, we would come to the conclusions:<\/p>\n\n\n\n
      \n
    1. car_weight<\/code> and year<\/code> might be dropped, depending on the specify aim of the model.<\/li>\n\n\n\n
    2. Add a regression spline for driver_age<\/code>.<\/li>\n\n\n\n
    3. Add an interaction between car_power<\/code> and town<\/code>.<\/li>\n<\/ol>\n\n\n\n
      Build strong GLM<\/h2>\n\n\n\n
      Let’s build a GLM with these insights. Two important things:<\/p>\n\n\n\n
        \n
      1. Glum<\/a> is an extremely powerful GLM implementation that was inspired by a pull request of our Christian Lorentzen.<\/li>\n\n\n\n
      2. In the upcoming version 3.0, it adds a formula API based of formulaic<\/a>, a very performant formula parser. This gives a very easy way to add interaction effects, regression splines, dummy encodings etc.<\/li>\n<\/ol>\n\n\n\n
        model_glm = GeneralizedLinearRegressor(\n    family="poisson",\n    l1_ratio=1.0,\n    alpha=1e-10,\n    formula="car_power * C(town) + bs(driver_age, 7) + car_age",\n)\nmodel_glm.fit(X_train, y=y_train)  # 1 second on old laptop\n\n# PDPs of both models\nfig, ax = plt.subplots(2, 2, figsize=(7, 5))\ncols = ("tab:blue", "tab:orange")\nfor color, name, model in zip(cols, ("GLM", "LGB"), (model_glm, model_lgb)):\n    disp = PartialDependenceDisplay.from_estimator(\n        model,\n        features=["driver_age", "car_age", "car_power", "town"],\n        X=X_explain,\n        ax=ax if name == "GLM" else disp.axes_,\n        line_kw={"label": name, "color": color},\n    )\nfig.suptitle("PDPs of both models", fontsize=15)\nfig.tight_layout()\n\n# Stratified PDP of car_power\nfor color, town in zip(("tab:blue", "tab:orange"), (0, 1)):\n    mask = X_explain.town == town\n    disp = PartialDependenceDisplay.from_estimator(\n        model_glm,\n        features=["car_power"],\n        X=X_explain[mask],\n        ax=None if town == 0 else disp.axes_,\n        line_kw={"label": town, "color": color},\n    )\nplt.suptitle("PDP of car_power stratified by town (0 vs 1)", fontsize=15)\n_ = plt.ylim(0, 0.2)<\/pre><\/div>\n\n\n\n
        \"\"<\/figure>\n\n\n\n
        \"\"<\/figure>\n\n\n\n
        In this relatively simple situation, the mean Poisson deviance of our models are very simlar now:<\/p>\n\n\n\n
        model_dummy = DummyRegressor().fit(X_train, y=y_train)\ndeviance_null = mean_poisson_deviance(y_test, model_dummy.predict(X_test)) \n\ndev_imp = []\nfor name, model in zip(("GLM", "LGB", "Null"), (model_glm, model_lgb, model_dummy)):\n    dev_imp.append((name, mean_poisson_deviance(y_test, model.predict(X_test))))\npd.DataFrame(dev_imp, columns=["Model", "Mean_Poisson_Deviance"])<\/pre><\/div>\n\n\n\n
        \"\"<\/figure>\n\n\n\n
        Final words<\/h2>\n\n\n\n