mSHAP

Step 1: Prepare the Data

Let’s take a look at the data we will be using.

## R
dataOhlsson %>% head()
#>   agarald kon zon mcklass fordald bonuskl duration antskad skadkost
#> 1       0   M   1       4      12       1 0.175342       0        0
#> 2       4   M   3       6       9       1 0.000000       0        0
#> 3       5   K   3       3      18       1 0.454795       0        0
#> 4       5   K   4       1      25       1 0.172603       0        0
#> 5       6   K   2       1      26       1 0.180822       0        0
#> 6       9   K   3       3       8       1 0.542466       0        0

We will first rename the columns, and look at summaries of each of the variables. Scikit-learn does not accept non-numeric covariates, so we will also convert the gender variable to an is_male indicator.

## R
cleaned <- dataOhlsson %>%
  mutate(severity = skadkost / antskad) %>%
  select(
    severity,
    claims = antskad,
    exposure = duration,
    age = agarald,
    gender = kon,
    geographic_zone = zon,
    vehicle_age = fordald
  ) %>%
  mutate(is_male = as.numeric(gender == "M")) %>%
  select(-gender)

summary(cleaned$severity)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>      16    3008    8724   23793   26788  211254   63878
summary(cleaned$exposure)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.0000  0.4630  0.8274  1.0107  1.0000 31.3397
summary(cleaned$age)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>    0.00   31.00   44.00   42.42   52.00   92.00
summary(cleaned$vehicle_age)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>    0.00    5.00   12.00   12.54   16.00   99.00
cleaned %>% count(is_male)
#>   is_male     n
#> 1       0  9853
#> 2       1 54695
cleaned %>% count(geographic_zone)
#>   geographic_zone     n
#> 1               1  8582
#> 2               2 11794
#> 3               3 12722
#> 4               4 24816
#> 5               5  2377
#> 6               6  3884
#> 7               7   373

Our next step will be to create a train/test split on the overall data. This will allow us to train models with the train set while having a holdout set for prediction and explanation. We will use 90% of our data to train the model and only 10% to test it, for faster run times on the SHAP explainers.

## R
idx <- createDataPartition(cleaned$claims, p = 0.9, list = FALSE)
train <- cleaned[idx,]
test <- cleaned[-idx,]

Now the data for the two models must be created. The first is for the severity model. We filter for only policies with a claim to create a severity specific training set. Also, we remove the exposure variable, as it should have no bearing on severity

## R
train_sev <- train %>%
  filter(claims > 0) %>%
  select(-exposure, -claims)

Next is frequency. The data will be used to create a multinomial classification model, and in this case, the possible values are 0, 1, and 2. Technically, it is possible to have more than 2 claims, but we will consider the probability so small as to be negligible. In order to create this model, we will use the same variables as before, but weight by the exposure column. Furthermore, we will downsample the rows with 0 claims while upsampling the rows with 1 and 2 claims, so we have a balanced data set.

## R
freq_0 <- train %>%
  filter(claims == 0) %>%
  sample_n(8000, replace = TRUE)
freq_1 <- train %>%
  filter(claims == 1) %>%
  sample_n(8000, replace = TRUE)
freq_2 <- train %>%
  filter(claims == 2) %>%
  sample_n(8000, replace = TRUE)

train_freq <- freq_0 %>% 
  bind_rows(freq_1) %>%
  bind_rows(freq_2) %>%
  mutate(claims = as.factor(claims)) %>%
  select(-severity)

To conclude the data preparation step, we will split our data into the predictors and the response, for ease of model fitting in python.

## R
X_sev <- train_sev %>%
  select(-severity)
y_sev <- train_sev %>%
  pull(severity)

X_freq <- train_freq %>%
  select(-claims)
y_freq <- train_freq %>%
  pull(claims)

Step 2: Train the Models

Our first model will predict the severity, or the cost per claim. It will be trained in python.

## Python
mod_dat_sev = r.X_sev
sev = r.y_sev

sev_mod = sk.RandomForestRegressor(n_estimators = 100, max_depth = 10, max_features = 2)
sev_mod.fit(mod_dat_sev, sev)
#> RandomForestRegressor(max_depth=10, max_features=2)

The next model will predict the frequency and will also be trained in python.

## Python
mod_dat_freq = r.X_freq
claims = r.y_freq

freq_mod = sk.RandomForestClassifier(n_estimators = 100, max_depth = 10, max_features = 2)
freq_mod.fit(mod_dat_freq, claims)
#> RandomForestClassifier(max_depth=10, max_features=2)

Step 3: Get Model Predictions

We will now take our test set predict on it with the two created models. The subsequent model outputs will be multiplied together and then averaged for each row to obtain an expected cost per row.

## R
test_sev <- test %>%
  select(-exposure, -severity, -claims)
test_freq <- test %>%
  select(-severity, -claims)

## Python
test_sev = r.test_sev
test_freq = r.test_freq

preds_sev = sev_mod.predict(test_sev)
preds_freq = freq_mod.predict_proba(test_freq)

## R
preds_sev <- py$preds_sev
preds_freq <- py$preds_freq %>%
  as.data.frame()
expected_values <- map2_dfc(
  .x = preds_freq,
  .y = 0:2,
  .f = ~{
    .x * .y * preds_sev
  }
) %>%
  rowSums()

Step 4: Explain the Predictions

With the goal of explaining this “average value” prediction, we will eventually use mSHAP. However, it is necessary to first calculate the SHAP values for the two models separately.

## Python
sev_ex = shap.Explainer(sev_mod)
sev_expected_val = sev_ex.expected_value
sev_preds_explained = sev_ex.shap_values(test_sev)

freq_ex = shap.Explainer(freq_mod)
freq_expected_val = freq_ex.expected_value
freq_preds_explained = freq_ex.shap_values(test_freq)

## R
freq_shap <- py$freq_preds_explained
sev_shap <- py$sev_preds_explained

Note that we can take these raw SHAP values from python and plug them straight into the function with no additional manipulation, but that requires that we specify the arguments shap_1_names and shap_2_names. Recall that our models do not use exactly the same predictors, so specifying the names will alert the algorithm of this and create a column of zeros for all non-matching column names.

Note also that passing in a list as one of the shap* arguments will cause a nested list to be returned, instead of the normal list.

## R
mshap_res <- mshap(
  shap_1 = freq_shap,
  shap_2 = sev_shap,
  ex_1 = py$freq_expected_val,
  ex_2 = py$sev_expected_val,
  shap_1_names = colnames(test_freq),
  shap_2_names = colnames(test_sev)
)

Since we want the expected value, we would like to combine the SHAP values in the same way, multiplying the respective values by 0, 1, and 2. This can be done in the following code.

## R
ev_explained <- mshap_res[[1]]$shap_vals * 0 + 
  mshap_res[[2]]$shap_vals * 1 + 
  mshap_res[[3]]$shap_vals * 2

shap_expected_values <- mshap_res[[1]]$expected_value * 0 + 
  mshap_res[[2]]$expected_value * 1 + 
  mshap_res[[3]]$expected_value * 2

Step 5: Visualize the Results

## R
summary_plot(variable_values = test_freq, shap_values = ev_explained)

## R
observation_plot(variable_values = test_freq[1,], shap_values = ev_explained[1,], expected_value = shap_expected_values[1])
#> Warning in min(x): no non-missing arguments to min; returning Inf
#> Warning in max(x): no non-missing arguments to max; returning -Inf