{"id":1409,"date":"2024-01-21T21:11:13","date_gmt":"2024-01-21T20:11:13","guid":{"rendered":"https:\/\/lorentzen.ch\/?p=1409"},"modified":"2024-01-21T21:11:13","modified_gmt":"2024-01-21T20:11:13","slug":"ml-xai-strong-glm","status":"publish","type":"post","link":"https:\/\/lorentzen.ch\/index.php\/2024\/01\/21\/ml-xai-strong-glm\/","title":{"rendered":"ML + XAI -> Strong GLM"},"content":{"rendered":"\n

My last post<\/a> was using {hstats}, {kernelshap} and {shapviz} to explain a binary classification random forest. Here, we use the same package combo to improve a Poisson GLM with insights from a boosted trees model.<\/p>\n\n\n\n

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

This time, we work with a synthetic, but quite realistic dataset. It describes 1 Mio insurance policies and their corresponding claim counts. A reference for the data is:<\/p>\n\n\n\n

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

library(OpenML)\nlibrary(lightgbm)\nlibrary(splines)\nlibrary(ggplot2)\nlibrary(patchwork)\nlibrary(hstats)\nlibrary(kernelshap)\nlibrary(shapviz)\n\n#===================================================================\n# Load and describe data\n#===================================================================\n\ndf <- getOMLDataSet(data.id = 45106)$data\n\ndim(df)  # 1000000       7\nhead(df)\n\n# year town driver_age car_weight car_power car_age claim_nb\n# 2018    1         51       1760       173       3        0\n# 2019    1         41       1760       248       2        0\n# 2018    1         25       1240       111       2        0\n# 2019    0         40       1010        83       9        0\n# 2018    0         43       2180       169       5        0\n# 2018    1         45       1170       149       1        1\n\nsummary(df)\n\n# Response\nggplot(df, aes(claim_nb)) +\n  geom_bar(fill = "chartreuse4") +\n  ggtitle("Distribution of the response")\n\n# Features\nxvars <- c("year", "town", "driver_age", "car_weight", "car_power", "car_age")\n\ndf[xvars] |> \n  stack() |> \nggplot(aes(values)) +\n  geom_histogram(fill = "chartreuse4", bins = 19) +\n  facet_wrap(~ind, scales = "free", ncol = 2) +\n  ggtitle("Distribution of the features")\n\n# car_power and car_weight are correlated 0.68, car_age and driver_age 0.28\ndf[xvars] |> \n  cor() |> \n  round(2)\n#            year  town driver_age car_weight car_power car_age\n# year          1  0.00       0.00       0.00      0.00    0.00\n# town          0  1.00      -0.16       0.00      0.00    0.00\n# driver_age    0 -0.16       1.00       0.09      0.10    0.28\n# car_weight    0  0.00       0.09       1.00      0.68    0.00\n# car_power     0  0.00       0.10       0.68      1.00    0.09\n# car_age       0  0.00       0.28       0.00      0.09    1.00\n\n<\/pre><\/div>\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 naive additive linear GLM and a tuned Boosted Trees model.<\/li>\n\n\n\n
  2. We combine the models and specify their predict function.<\/li>\n<\/ol>\n\n\n\n
    # Train\/test split\nset.seed(8300)\nix <- sample(nrow(df), 0.9 * nrow(df))\ntrain <- df[ix, ]\nvalid <- df[-ix, ]\n\n# Naive additive linear Poisson regression model\n(fit_glm <- glm(claim_nb ~ ., data = train, family = poisson()))\n\n# Boosted trees with LightGBM. The parameters (incl. number of rounds) have been \n# by combining early-stopping with random search CV (not shown here)\n\ndtrain <- lgb.Dataset(data.matrix(train[xvars]), label = train$claim_nb)\n\nparams <- list(\n  learning_rate = 0.05, \n  objective = "poisson", \n  num_leaves = 7, \n  min_data_in_leaf = 50, \n  min_sum_hessian_in_leaf = 0.001, \n  colsample_bynode = 0.8, \n  bagging_fraction = 0.8, \n  lambda_l1 = 3, \n  lambda_l2 = 5\n)\n\nfit_lgb <- lgb.train(params = params, data = dtrain, nrounds = 300)  \n\n# {hstats} works for multi-output predictions,\n# so we can combine all models to a list, which simplifies the XAI part.\nmodels <- list(GLM = fit_glm, LGB = fit_lgb)\n\n# Custom predictions on response scale\npf <- function(m, X) {\n  cbind(\n    GLM = predict(m$GLM, X, type = "response"),\n    LGB = predict(m$LGB, data.matrix(X[xvars]))\n  )\n}\npf(models, head(valid, 2))\n#       GLM        LGB\n# 0.1082285 0.08580529\n# 0.1071895 0.09181466\n\n# And on log scale\npf_log <- function(m, X) {\n  log(pf(m = m, X = X))\n}\npf_log(models, head(valid, 2))\n#       GLM       LGB\n# -2.223510 -2.455675\n# -2.233157 -2.387983 -2.346350<\/pre><\/div>\n\n\n\n
    Traditional XAI<\/h2>\n\n\n\n
    Performance<\/h4>\n\n\n\n
    Comparing average Poisson deviance on the validation data shows that the LGB model is clearly better than the naively built GLM, so there is room for improvent!<\/p>\n\n\n\n
    perf <- average_loss(\n  models, X = valid, y = "claim_nb", loss = "poisson", pred_fun = pf\n)\nperf\n#       GLM       LGB \n# 0.4362407 0.4331857\n<\/pre><\/div>\n\n\n\n
    Feature importance<\/h4>\n\n\n\n
    Next, we calculate permutation importance on the validation data with respect to mean Poisson deviance loss. The results make sense, and we note that year<\/code> and car_weight<\/code> seem to be negligile.<\/p>\n\n\n\n
    imp <- perm_importance(\n  models, v = xvars, X = valid, y = "claim_nb", loss = "poisson", pred_fun = pf\n)\nplot(imp)<\/pre><\/div>\n\n\n\n
    \"\"<\/figure>\n\n\n\n
    Main effects<\/h4>\n\n\n\n
    Next, we visualize estimated main effects by partial dependence plots on log link scale. The differences between the models are quite small, with one big exception: Investing more parameters into driver_age<\/code> via spline will greatly improve the performance and usefulness of the GLM.<\/p>\n\n\n\n
    partial_dep(models, v = "driver_age", train, pred_fun = pf_log) |> \n  plot(show_points = FALSE)\n\npdp <- function(v) {\n  partial_dep(models, v = v, X = train, pred_fun = pf_log) |> \n    plot(show_points = FALSE)\n}\nwrap_plots(lapply(xvars, pdp), guides = "collect") &\n  ylim(-2.8, -1.7)<\/pre><\/div>\n\n\n\n
    \"\"<\/figure>\n\n\n\n
    Interaction effects<\/h4>\n\n\n\n
    Friedman’s H-squared (per feature and feature pair) and on log link scale shows that – unsurprisingly – our GLM does not contain interactions, and that the strongest relative interaction happens between town<\/code> and car_power<\/code>. The stratified PDP visualizes this interaction. Let’s add a corresponding interaction effect to our GLM later.<\/p>\n\n\n\n
    system.time(  # 5 sec\n  H <- hstats(models, v = xvars, X = train, pred_fun = pf_log)\n)\nH\nplot(H)\n\n# Visualize strongest interaction by stratified PDP\npartial_dep(models, v = "car_power", X = train, pred_fun = pf_log, BY = "town") |> \n  plot(show_points = FALSE)<\/pre><\/div>\n\n\n\n
    \"\"<\/figure>\n\n\n\n
    \"\"<\/figure>\n\n\n\n
    SHAP<\/h2>\n\n\n\n
    As an elegant alternative to studying feature importance, PDPs and Friedman’s H, we can simply run a SHAP analysis on the LGB model.  <\/p>\n\n\n\n
    set.seed(22)\nX_explain <- train[sample(nrow(train), 1000), xvars]\n \nshap_values_lgb <- shapviz(fit_lgb, data.matrix(X_explain))\nsv_importance(shap_values_lgb)\nsv_dependence(shap_values_lgb, v = xvars) &\n  ylim(-0.35, 0.8)<\/pre><\/div>\n\n\n\n
    \"\"<\/figure>\n\n\n\n
    \"\"<\/figure>\n\n\n\n
    Here, we would come to the same conclusions:<\/p>\n\n\n\n
      \n
    1. car_weight<\/code> and year<\/code> might be dropped.<\/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
      Pimp the GLM<\/h2>\n\n\n\n
      In the final section, we apply the three insights from above with very good results. <\/p>\n\n\n\n
      fit_glm2 <- glm(\n  claim_nb ~ car_power * town + ns(driver_age, df = 7) + car_age, \n  data = train, \n  family = poisson()\n  \n# Performance now as good as LGB\nperf_glm2 <- average_loss(\n  fit_glm2, X = valid, y = "claim_nb", loss = "poisson", type = "response"\n)\nperf_glm2  # 0.432962\n\n# Effects similar as LGB, and smooth\npartial_dep(fit_glm2, v = "driver_age", X = train) |> \n  plot(show_points = FALSE)\n\npartial_dep(fit_glm2, v = "car_power", X = train, BY = "town") |> \n  plot(show_points = FALSE)<\/pre><\/div>\n\n\n\n
      \"\"<\/figure>\n\n\n\n
      \"\"<\/figure>\n\n\n\n
      Or even via permutation or kernel SHAP:<\/p>\n\n\n\n
      set.seed(1)\nbg <- train[sample(nrow(train), 200), ]\nxvars2 <- setdiff(xvars, c("year", "car_weight"))\n\nsystem.time(  # 4 sec\n  ks_glm2 <- permshap(fit_glm2, X = X_explain[xvars2], bg_X = bg)\n)\nshap_values_glm2 <- shapviz(ks_glm2)\nsv_dependence(shap_values_glm2, v = xvars2) &\n  ylim(-0.3, 0.8)<\/pre><\/div>\n\n\n\n
      \"\"<\/figure>\n\n\n\n
      Final words<\/h2>\n\n\n\n