Michael's and Christian's Blog

Category: Machine Learning

  • Model Diagnostics: Statistics vs Machine Learning

    In this post, we show how different use cases require different model diagnostics. In short, we compare (statistical) inference and prediction.

    As an example, we use a simple linear model for the Munich rent index dataset, which was kindly provided by the authors of Regression – Models, Methods and Applications 2nd ed. (2021). This dataset contains monthy rents in EUR (rent) for about 3000 apartments in Munich, Germany, from 1999. The apartments have several features such as living area (area) in squared meters, year of construction (yearc), quality of location (location, 0: average, 1: good, 2: top), quality of bath rooms (bath, 0:standard, 1: premium), quality of kitchen (kitchen, 0: standard, 1: premium), indicator for central heating (cheating).

    The target variable is Y=rentY=\text{rent} and the goal of our model is to predict the mean rent, E[Y]E[Y] (we omit the conditioning on X for brevity).

    Disclaimer: Before presenting the use cases, let me clearly state that I am not in the apartment rent business and everything here is merely for the purpose of demonstrating statistical good practice.

    Inference

    The first use case is about inference of the effect of the features. Imagine the point of view of an investor who wants to know whether the installation of a central heating is worth it (financially). To lay the ground on which to base a decision, a statistician must have answers to:

    • What is the effect of the variable cheating on the rent.
    • Is this effect statistically significant?

    Prediction

    The second use case is about prediction. This time, we take the point of view of someone looking out for a new apartment to rent. In order to know whether the proposed rent by the landlord is about right or improper (too high), a reference value would be very convenient. One can either ask the neighbors or ask a model to predict the rent of the apartment in question.

    Model Fit

    Before answering the above questions and doing some key diagnostics, we must load the data and fit a model. We choose a simple linear model and directly model rent.

    Notes:

    • For rent indices as well as house prices, one often log-transforms the target variable before modelling or one uses a log-link and an appropriate loss function (e.g. Gamma deviance).
    • Our Python version uses GeneralizedLinearRegressor from the package glum. We could as well have chosen other implementations like statsmodels.regression.linear_model.OLS. This way, we have to implement the residual diagnostics ourselves which makes it clear what is plotted.

    For brevity, we skip imports and data loading. Our model is then fit by:

    Python
    R

    Diagnostics for Inference

    The coefficient table will already tell us the effect of the cheating variable. For more involved models like gradient boosted trees or neural nets, one can use partial dependence and shap values to assess the effect of features.

    Python
    R
    Variablecoefsep_valueci_lowerci_upper
    intercept-3682.5327.00.0-4323-3041
    bs(area, ..)[1]88.531.34.6e-0327150
    bs(area,..)[2]316.824.50.0269365
    bs(area, ..)[3]547.762.80.0425671
    bs(area, ..)[4]733.791.71.3e-15554913
    yearc1.90.20.01.62.3
    C(location)[2]48.25.94.4e-163760
    C(location)[3]137.927.76.6e-0784192
    C(bath)[1]50.016.52.4e-031882
    C(kitchen)[1]98.218.51.1e-0762134
    C(cheating)[1]107.810.60.087.0128.6

    We see that ceteris paribus, meaning all else equal, a central heating increases the monthly rent by about 108 EUR. Not the size of the effect of 108 EUR, but the fact that there is an effect of central heating on the rent seems statistically significant:
    This is indicated by the very low probability, i.e. p-value, for the null-hypothesis of cheating having a coefficient of zero.
    We also see that the confidence interval with the default confidence level of 95%: [ci_lower, ci_upper] = [87, 129].
    This shows the uncertainty of the estimated effect.

    For a building with 10 apartments and with an investment horizon of about 10 years, the estimated effect gives roughly a budget of 13000 EUR (range is roughly 10500 to 15500 with 95% confidence).

    A good statistician should ask several further questions:

    • Is the dataset at hand a good representation of the population?
    • Are there confounders or interaction effects, in particular between cheating and other features?
    • Are the assumptions for the low p-value and the confidence interval of cheating valid?

    Here, we will only address the last question, and even that one only partially. Which assumptions were made? The error term, ϵ=Y−E[Y]\epsilon = Y - E[Y], should be homoscedastic and Normal distributed. As the error is not observable (because the true model for E[Y]E[Y] is unknown), one replaces E[Y]E[Y] by the model prediction E^[Y]\hat{E}[Y], this gives the residuals, ϵ^=Y−E^[Y]=y−fitted values\hat{\epsilon} = Y - \hat{E}[Y] = y - \text{fitted values}, instead. For homoscedasticity, the residuals should look like white (random) noise. Normality, on the other hand, becomes less of a concern with larger data thanks to the central limit theorem. With about 3000 data points, we are far away from small data, but it might still be a good idea to check for normality.

    The diagnostic tools to check that are residual and quantile-quatile (QQ) plots.

    Python
    R
    Residual plots on the training data.

    The more data points one has the less informative is a scatter plot. Therefore, we put a contour plot on the right.

    Visual insights:

    • There seems to be a larger variability for larger fitted values. This is a hint that the homoscedasticity might be violated.
    • The residuals seem to be centered around 0. This is a hint that the model is well calibrated (adequate).
    Python
    R

    The QQ-plot shows the quantiles of the theoretical assumed distribution of the residuals on the x-axis and the ordered values of the residuals on the y-axis. In the Python version, we decided to use the studentized residuals because normality of the error implies a student (t) distribution for these residuals.

    Concluding remarks:

    • We might do similar plots on the test sample, but we don’t necessarily need a test sample to answer the inference questions.
    • It is good practice to plot the residuals vs each of the features as well.

    Diagnostics for Prediction

    If we are only interested in predictions of the mean rent, E^[Y]\hat{E}[Y], we don’t care much about the probability distribution of YY. We just want to know if the predictions are close enough to the real mean of the rent E[Y]E[Y]. In a similar argument as for the error term and residuals, we have to accept that E[Y]E[Y] is not observable (it is the quantity that we want to predict). So we have to fall back to the observations of YY in order to judge if our model is well calibrated, i.e., close the the ideal E[Y]E[Y].

    Very importantly, here we make use of the test sample in all of our diagnostics because we fear the in-sample bias.

    We start simple by a look at the unconditional calibration, that is the average (negative) residual 1n∑(E^[Yi]−Yi)\frac{1}{n}\sum(\hat{E}[Y_i]-Y_i).

    Python
    R
    setmean biascountstderrp-value
    train-3.2e-1224652.81.0
    test2.16175.80.72

    It is no surprise that bias_mean in the train set is almost zero.
    This is the balance property of (generalized) linear models (with intercept term). On the test set, however, we detect a small bias of about 2 EUR per apartment on average.

    Next, we have a look a reliability diagrams which contain much more information about calibration and bias of a model than the unconditional calibration above. In fact, it assesses auto-calibration, i.e. how well the model uses its own information.
    An ideal model would lie on the dotted diagonal line.

    Python
    R

    Visual insights:

    • Graphs on train and test set look very similar.
      The larger uncertainty intervals on the test set stem from the fact that is has a smaller sample size.
    • The model seems to lie around the diagonal indicating good auto-calibration for the largest part of the range.
    • Very high predicted values seem to be systematically too low, i.e. the graph is above the diagonal.

    Finally, we assess conditional calibration, i.e. the calibration with respect to the features. Therefore, we plot one of our favorite graphs for each feature. It consists of:

    • average observed value of YY for each (binned) value of the feature
    • average predicted value
    • partial dependence
    • histogram of the feature (grey, right y-axis)
    Python
    R

    Visual insights:

    • On the train set, the categorical features seem to have perfect calibration as average observed equals average predicted. This is again a result of the balance property. On the test set, we see a deviation, especially for the categorical level with smaller sample size. This is a good demonstration why plotting on both train and test set is a good idea.
    • The numerical features area and year of construction seem fine, but a closer look can’t hurt.

    We next perform a bias plot, which is plotting the average difference of predicted minus observed per feature value. The values should be around zero, so we can zoom in on the y-axis.
    This is very similar to the residual plot, but the information is better condensed for its purpose.

    Python
    R

    Visual insights:

    • For large values of area and yearc in the 1940s and 1950s, there are only few observations available. Still, the model might be improved for those regions.
    • The bias of yearc shows a parabolic curve. The simple linear effect in our model seems too simplistic. A refined model could use splines instead, as for area.

    Concluding remarks:

    • The predictions for area larger than around 120 square meters and for year of construction around the 2nd world war are less reliable.
    • For all the rest, the bias is smaller than 50 EUR on average.
      This is therefore a rough estimation of the prediction uncertainty.
      It should be enough to prevent improperly high (or low) rents (on average).

    The full Python and R code is available under:

  • Fast Grouped Counts and Means in R

    Edited on 2025-05-01: Multiple improvements by Christian, especially on making Polars neater, DuckDB faster, and the plot easier to read.

    From time to time, the following questions pop up:

    1. How to calculate grouped counts and (weighted) means?
    2. What are fast ways to do it in R?

    This blog post presents a couple of approaches and then compares their speed with a naive (non-scientific!) benchmark.

    Base R

    There are many ways to calculate grouped counts and means in base R, e.g., aggregate(), tapply(), by(), split() + lapply(). In my experience, the fastest way is a combination of tabulate() and rowsum().

    R

    But: tabulate() ignores missing values. To avoid problems, create an explicit missing level via factor(x, exclude = NULL).

    Let’s turn to some other approaches.

    dplyr

    Not optimized for speed or memory, but the de-facto standard in data processing with R. I love its syntax.

    R

    data.table

    Does not need an introduction. Since 2006 the package for fast data manipulation written in C.

    R

    DuckDB

    Extremely powerful query engine / database system written in C++, with initial release in 2019, and R bindings since 2020. Allows larger-than-RAM calculations.

    R

    collapse

    C/C++-based package for data transformation and statistical computing. {collapse} was initially released on CRAN in 2020. It can do much more than grouped calculations, check it out!

    R

    Polars

    R bindings of the fantastic Polars project that started in 2020. First R release in 2022. Currently under heavy revision.

    The current package is not up-to-date with the main project, thus we expect the revised version (available in this branch) to be faster.

    R

    Naive Benchmark

    Let’s compare the speed of these approaches for sample sizes up to 10^8 using a Windows system with an Intel i7-13700H CPU.

    This is not at all meant as a scientific benchmark!

    R

    Memory

    What about memory? {dplyr}, {data.table}, and rowsum() require a lot of it, as does collapse::fcount(). For the other approaches, almost no memory is required, or profmem can’ t measure it.

    Final words

    • {duckdb} is increadibly fast for large data.
    • {collapse} is increadibly fast for all sample sizes. In other benchmarks, it is slower because there, the grouping has to be a string rather than a factor.
    • {polars} looks really cool.
    • rowsum() and tabulate() provide fast solutions with base R.

    R script

  • Dictionary for Data Scientists and Statisticians

    During my journey through machine learning (ML) and statistics, I was faced some many times with surprisingly different usage of terms. To improve the understanding of data scientists and statisticians, I present a dictionary and hope the humour does not get unnoticed.

    data scientiststatisticiancomment
    sampleobservation
    (training) setsample
    featurecovariate, predictormany more terms
    labelcategorical response
    inferenceprediction, forecast
    statisticsinference
    trainingfitting
    training errorin-sample error
    test/validation sethold-out sample
    regressionregression
    classificationregression (on categorical response) + decision makingthus the name logistic / multinomial regression!
    supervised machine learningregression
    AIAI for funding, else regressionsee EU AI Act article 3
    confidence scorepredicted probabilityconfidence scores might not represent probabilities
    (binary/multiclass) cross-entropy(binomial/multinomial) log likelihooda.k.a. log loss
    unbalanced data problem🤷‍♂️what problem?if any, a small data problem
    SMOTEdevil’s work

    Statistics is about the honest interpretation of data, which is much less appealing than less honest interpretation.

    by Prof. Simon Wood, a.k.a. Mr GAM/mgcv

  • Converting arbitrarily large CSVs to Parquet with R

    In this recent post, we have used Polars and DuckDB to convert a large CSV file to Parquet in steaming mode – and Python.

    Different people have contacted me and asked: “and in R?”

    Simple answer: We have DuckDB, and we have different Polars bindings. Here, we are using {polars} which is currently being overhauled into {neopandas}.

    So let’s not wait any longer!


    Run times are on a Windows system with an Intel i7-13700H CPU.

    Generate 2.2 GB csv file

    We use {data.table} to dump a randomly generated dataset with 100 Mio rows into a csv file.

    R

    DuckDB

    Then, we use DuckDB to fire a query to the file and stream the result into Parquet.

    Threads and RAM can be set on the fly, which is very convenient. Setting a low memory limit (e.g., 500 MB) will work – try it out!

    R

    3.5 seconds – wow! The resulting file looks good. It is 125 MB large.

    Polars

    Let’s do the same with Polars.

    R

    With nine seconds, it is slower than DuckDB. But the output looks as expected and has the same size as with DuckDB.

    Final words

    • With DuckDB or Polars, conversion of CSVs to Parquet is easy and fast, even in larger-than-RAM situations.
    • We can apply filters, selects, sorts etc. on the fly.
    • Let’s keep an eye on Polars in R. It looks really interesting.

    R script

  • Converting arbitrarily large CSVs to Parquet with Python

    Conversion from CSV to Parquet in streaming mode? No problem for the two power houses Polars and DuckDB. We can even throw in some data preprocessing steps in-between, like column selection, data filters, or sorts.

    Edit: Streaming writing (or “lazy sinking”) of data with Polars was introduced with release 1.25.2 in March 2025, thanks Christian for pointing this out.

    pip install polars

    pip install duckdb


    Run times are on a normal laptop, dedicating 8 threads to the crunching.

    Let’s generate a 2 GB csv file first

    Python

    Polars

    Let’s use Polars in Lazy mode to connect to the CSV, apply some data operations, and stream the result into a Parquet file.

    Python

    In case you prefer to write SQL code, you can alternatively use the SQL API of Polars. Curiously, run time is substantially longer:

    Python

    In both cases, the result looks as expected, and the resulting Parquet file is about 170 MB large.

    Python

    DuckDB

    As an alternative, we use DuckDB. Thread pool size and RAM limit can be set on the fly. Setting a low memory limit (e.g., 500 MB) will lead to longer run times, but it works.

    Python

    Again, the output looks as expected. The Parquet file is again 170 MB large, thanks to using the same compression (“zstd”) as with Polars..

    Python

    Final words

    • With Polars or DuckDB, conversion of CSVs to Parquet is easy and fast, even in larger-than-RAM situations.
    • We can apply filters, selects, sorts etc. on the fly.

    Python notebook

  • Effect Plots in Python and R

    Christian and me did some code magic: Highly effective plots that help to build and inspect any model:

    The functionality is best described by its output:

    Python
    R

    The plots show different types of feature effects relevant in modeling:

    • Average observed: Descriptive effect (also interesting without model).
    • Average predicted: Combined effect of all features. Also called “M Plot” (Apley 2020).
    • Partial dependence: Effect of one feature, keeping other feature values constant (Friedman 2001).
    • Number of observations or sum of case weights: Feature value distribution.
    • R only: Accumulated local effects, an alternative to partial dependence (Apley 2020).

    Both implementations…

    • are highly efficient thanks to {Polars} in Python and {collapse} in R, and work on datasets with millions of observations,
    • support case weights with all their statistics, ideal in insurance applications,
    • calculate average residuals (not shown in the plots above),
    • provide standard deviations/errors of average observed and bias,
    • allow to switch to Plotly for interactive plots, and
    • are highly customizable (the R package, e.g., allows to collapse rare levels after calculating statistics via the update() method or to sort the features according to main effect importance).

    In the spirit of our “Lost In Translation” series, we provide both high-quality Python and R code. We will use the same data and models as in one of our latest posts on how to build strong GLMs via ML + XAI.

    Example

    Let’s build a Poisson LightGBM model to explain the claim frequency given six traditional features in a pricing dataset on motor liability claims. 80% of the 1 Mio rows are used for training, the other 20% for evaluation. Hyper-parameters have been slightly tuned (not shown).

    R
    Python

    Let’s inspect the (main effects) of the model on the test data.

    R
    Python

    The output can be seen at the beginning of this blog post.

    Here some model insights:

    • Average predictions closely match observed frequencies. No clear bias is visible.
    • Partial dependence shows that the year and the car weight almost have no impact (regarding their main effects), while the driver_age and car_power effects seem strongest. The shared y axes help to assess these.
    • Except for car_weight, the partial dependence curve closely follows the average predictions. This means that the model effect seems to really come from the feature on the x axis, and not of some correlated other feature (as, e.g., with car_weight which is actually strongly correlated with car_power).

    Final words

    • Inspecting models has become much relaxed with above functions.
    • The packages used offer much more functionality. Try them out! Or we will show them in later posts ;).

    References

    1. Apley, Daniel W., and Jingyu Zhu. 2020. Visualizing the Effects of Predictor Variables in Black Box Supervised Learning Models. Journal of the Royal Statistical Society Series B: Statistical Methodology, 82 (4): 1059–1086. doi:10.1111/rssb.12377.
    2. Friedman, Jerome H. 2001. Greedy Function Approximation: A Gradient Boosting Machine. Annals of Statistics 29 (5): 1189–1232. doi:10.1214/aos/1013203451.

    R script , Python notebook

  • Explaining a Causal Forest

    We use a causal forest [1] to model the treatment effect in a randomized controlled clinical trial. Then, we explain this black-box model with usual explainability tools. These will reveal segments where the treatment works better or worse, just like a forest plot, but multivariately.

    Data

    For illustration, we use patient-level data of a 2-arm trial of rectal indomethacin against placebo to prevent post-ERCP pancreatitis (602 patients) [2]. The dataset is available in the package {medicaldata}.

    The data is in fantastic shape, so we don’t need to spend a lot of time with data preparation.

    1. We integer encode factors.
    2. We select meaningful features, basically those shown in the forest plot of [2] (Figure 4) without low-information features and without hospital.

    The marginal estimate of the treatment effect is -0.078, i.e., indomethacin reduces the probability of post-ERCP pancreatitis by 7.8 percentage points. Our aim is to develop and interpret a model to see if this value is associated with certain covariates.

    R

    The model

    We use the {grf} package to fit a causal forest [1], a tree-ensemble trying to estimate conditional average treatment effects (CATE) E[Y(1) – Y(0) | X = x]. As such, it can be used to study treatment effect inhomogeneity.

    In contrast to a typical random forest:

    • Honest trees are grown: Within trees, part of the data is used for splitting, and the other part for calculating the node values. This anti-overfitting is implemented for all random forests in {grf}.
    • Splits are selected to produce child nodes with maximally different treatment effects (under some additional constraints).

    Note: With about 13%, the complication rate is relatively low. Thus, the treatment effect (measured on absolute scale) can become small for certain segments simply because the complication rate is close to 0. Ideally, we could model relative treatment effects or odds ratios, but I have not found this option in {grf} so far.

    R

    Explain the model with “classic” techniques

    After looking at tree split importance, we study the effects via partial dependence plots and Friedman’s H. These only require a predict() function and a reference dataset.

    R

    Variable importance

    Variable importance of the causal forest can be measured by the relative counts each feature had been used to split on (in the first 4 levels). The most important variable is age.

    Main effects

    To study the main effects on the CATE, we consider partial dependence plots (PDP). Such plot shows how the average prediction depends on the values of a feature, keeping all other feature values constant (can be unnatural.)

    We can see that the treatment effect is strongest for persons up to age 35, then reduces until 45. For older patients, the effect increases again.

    Remember: Negative values mean a stronger (positive) treatment effect.

    Interaction strength

    Between what covariates are there strong interactions?

    A model agnostic way to assess pairwise interaction strength is Friedman’s H statistic [3]. It measures the error when approximating the two-dimensional partial dependence function of the two features by their univariate partial dependence functions. A value of zero means there is no interaction. A value of α means that about 100α% of the joint effect (variability) comes from the interaction.

    This measure is shown on the right hand side of the plot. More than 15% of the joint effect variability of age and biliary sphincterotomy (bsphinc) comes from their interaction.

    Typically, pairwise H-statistics are calculated only for the most important variables or those with high overall interaction strength. Overall interaction strength (left hand side of the plot) can be measured by a version of Friedman’s H. It shows how much of the prediction variability comes from interactions with that feature.

    Visualize strong interaction

    Interactions can be visualized, e.g., by a stratified PDP. We can see that the treatment effect is associated with age mainly for persons with biliary sphincterotomy.

    SHAP Analysis

    A “modern” way to explain the model is based on SHAP [4]. It decomposes the (centered) predictions into additive contributions of the covariates.

    Because there is no TreeSHAP shipped with {grf}, we use the much slower Kernel SHAP algorithm implemented in {kernelshap} that works for any model.

    First, we explain the prediction of a single data row, then we decompose many predictions. These decompositions can be analysed by simple descriptive plots to gain insights about the model as a whole.

    R

    Explain one CATE

    Explaining the CATE corresponding to the feature values of the first patient via waterfall plot.

    SHAP importance plot

    The bars show average absolute SHAP values. For instance, we can say that biliary sphincterotomy impacts the treatment effect on average by more than +- 0.01 (but we don’t see how).

    SHAP summary plot

    One-dimensional plot of SHAP values with scaled feature values on the color scale, sorted in the same order as the SHAP importance plot. Compared to the SHAP importance barplot, for instance, we can additionally see that biliary sphincterotomy weakens the treatment effect (positive SHAP value).

    SHAP dependence plots

    Scatterplots of SHAP values against corresponding feature values. Vertical scatter (at given x value) indicates presence of interactions. A candidate of an interacting feature is selected on the color scale. For instance, we see a similar pattern in the age effect on the treatment effect as in the partial dependence plot. Thanks to the color scale, we also see that the age effect depends on biliary sphincterotomy.

    Remember that SHAP values are on centered prediction scale. Still, a positive value means a weaker treatment effect.

    Wrap-up

    • {grf} is a fantastic package. You can expect more on it here.
    • Causal forests are an interesting way to directly model treatment effects.
    • Standard explainability methods can be used to explain the black-box.

    References

    1. Athey, Susan, Julie Tibshirani, and Stefan Wager. “Generalized Random Forests”. Annals of Statistics, 47(2), 2019.
    2. Elmunzer BJ et al. A randomized trial of rectal indomethacin to prevent post-ERCP pancreatitis. N Engl J Med. 2012 Apr 12;366(15):1414-22. doi: 10.1056/NEJMoa1111103.
    3. Friedman, Jerome H., and Bogdan E. Popescu. Predictive Learning via Rule Ensembles. The Annals of Applied Statistics 2, no. 3 (2008): 916-54.
    4. Scott M. Lundberg and Su-In Lee. A Unified Approach to Interpreting Model Predictions. Advances in Neural Information Processing Systems 30 (2017).

    The full R notebook

  • Out-of-sample Imputation with {missRanger}

    {missRanger} is a multivariate imputation algorithm based on random forests, and a fast version of the original missForest algorithm of Stekhoven and Buehlmann (2012). Surprise, surprise: it uses {ranger} to fit random forests. Especially combined with predictive mean matching (PMM), the imputations are often quite realistic.

    Out-of-sample application

    The newest CRAN release 2.6.0 offers out-of-sample application. This is useful for removing any leakage between train/test data or during cross-validation. Furthermore, it allows to fill missing values in user provided data. By default, it uses the same number of PMM donors as during training, but you can change this by setting pmm.k = nice value.

    We distinguish two types of observations to be imputed:

    1. Easy case: Only a single value is missing. Here, we simply apply the corresponding random forest to fill the one missing value.
    2. Hard case: Multiple values are missing. Here, we first fill the values univariately, and then repeatedly apply the corresponding random forests, with the hope that the effect of univariate imputation vanishes. If values of two highly correlated features are missing, then the imputations can be non-sensical. There is no way to mend this.

    Example

    To illustrate the technique with a simple example, we use the iris data.

    1. First, we randomly add 10% missing values.
    2. Then, we make a train/test split.
    3. Next, we “fit” missRanger() to the training data.
    4. Finally, we use its new predict() method to fill the test data.

    R

    The results look reasonable, in this case even for the “hard case” row 6 with missing values in two variables. Here, it is probably the strong association with Species that helped to create good values.

    The new predict() also works with single row input.

    Learn more about {missRanger}

    The full R script

  • SHAP Values of Additive Models

    Within only a few years, SHAP (Shapley additive explanations) has emerged as the number 1 way to investigate black-box models. The basic idea is to decompose model predictions into additive contributions of the features in a fair way. Studying decompositions of many predictions allows to derive global properties of the model.

    What happens if we apply SHAP algorithms to additive models? Why would this ever make sense?

    In the spirit of our “Lost In Translation” series, we provide both high-quality Python and R code.

    The models

    Let’s build the models using a dataset with three highly correlated covariates and a (deterministic) response.

    R
    Python

    SHAP

    For both models, we use exact permutation SHAP and exact Kernel SHAP. Furthermore, the linear model is analyzed with “additive SHAP”, and the tree-based model with TreeSHAP.

    Do the algorithms provide the same?

    R
    Python
    SHAP dependence plot of the additive linear model and the additive explainer (Python).

    Yes – the three algorithms within model provide the same SHAP values. Furthermore, the SHAP values reconstruct the additive components of the features.

    Didactically, this is very helpful when introducing SHAP as a method: Pick a white-box and a black-box model and compare their SHAP dependence plots. For the white-box model, you simply see the additive components, while the dependence plots of the black-box model show scatter due to interactions.

    Remark: The exact equivalence between algorithms is lost, when

    • there are too many features for exact procedures (~10+ features), and/or when
    • the background data of Kernel/Permutation SHAP does not agree with the training data. This leads to slightly different estimates of the baseline value, which itself influences the calculation of SHAP values.

    Final words

    • SHAP algorithms applied to additive models typically give identical results. Slight differences might occur because sampling versions of the algos are used, or a different baseline value is estimated.
    • The resulting SHAP values describe the additive components.
    • Didactically, it helps to see SHAP analyses of white-box and black-box models side by side.

    R script , Python notebook

  • A Tweedie Trilogy — Part III: From Wrights Generalized Bessel Function to Tweedie’s Compound Poisson Distribution

    TLDR: The scipy 1.7.0 release introduced Wright’s generalized Bessel function in the Python ecosystem. It is an important ingredient for the density and log-likelihood of Tweedie probabilty distributions. In this last part of the trilogy I’d like to point out why it was important to have this function and share the endeavor of implementing this inconspicuous but highly intractable special function. The fun part is exploiting a free parameter in an integral representation, which can be optimized by curve fitting to the minimal arc length.

    This trilogy celebrates the 40th birthday of Tweedie distributions in 2024 and highlights some of their very special properties.

    See part i and part ii.

    Tweedie Distributions

    As pointed out in part I and part II, the family of Tweedie distributions is a very special one with outstanding properties. They are central for estimating expectations with GLMs. The probability distributions have mainly positive (non-negative) support and are skewed, e.g. Poisson, Gamma, Inverse Gaussian and compound Poisson-Gamma.

    As members of the exponential dispersion family, a slight extension of the exponential family, the probability density can be written as

    f(y;θ,ϕ)=c(y,ϕ)exp⁡(yθ−κ(θ)ϕ)κ(θ)=κp(θ)=12−p((1−p)θ)2−p1−p\begin{align*} f(y; \theta, \phi) &= c(y, \phi) \exp\left(\frac{y\theta - \kappa(\theta)}{\phi}\right) \\ \kappa(\theta) &= \kappa_p(\theta) = \frac{1}{2-p}((1-p)\theta)^{\frac{2-p}{1-p}} \end{align*}

    It is often more instructive to parametrise the distribution with pp, μ\mu and ϕ\phi, using

    θ={μ1−p1−p ,p≠1log⁡(μ) ,p=1κ(θ)={μ2−p2−p ,p≠2log⁡(μ) ,p=2\begin{align*} \theta &= \begin{cases} \frac{\mu^{1-p}}{1-p}\,,\quad p\neq 1\\ \log(\mu)\,,\quad p=1 \end{cases} \\ \kappa(\theta) &= \begin{cases} \frac{\mu^{2-p}}{2-p}\,,\quad p\neq 2\\ \log(\mu)\,,\quad p=2 \end{cases} \end{align*}

    and write

    Y∼Twp(μ,ϕ)\begin{align*} Y &\sim \mathrm{Tw}_p(\mu, \phi) \end{align*}
    Probability density of several Tweedie distributions.

    Compound Poisson Gamma

    A very special domain for the power parameter is between Poisson and Gamma: 1<p<21<p<2. This range results in the Compound Poisson distribution which is suitable if you have a random count process and if each count itself has a random amount. A well know example is insurance claims. Typically, there is a random number of insurance claims, and each and every claim has a random amount of claim costs.

    N∼Poisson(λ)Xi∼Gamma(a,b)Y=∑i=0NXi∼CompPois(λ,a,b)\begin{align*} N &\sim \mathrm{Poisson}(\lambda)\\ X_i &\sim \mathrm{Gamma}(a, b)\\ Y &= \sum_{i=0}^N X_i \sim \mathrm{CompPois}(\lambda, a, b) \end{align*}

    For Poisson count we have E⁡[N]=λ\operatorname{E}[N]=\lambda and Var⁡[N]=λ=E⁡[N]\operatorname{Var}[N]=\lambda=\operatorname{E}[N], for the Gamma amount E⁡[X]=ab\operatorname{E}[X]=\frac{a}{b} and Var⁡[X]=ab2=1aE⁡[X]2\operatorname{Var}[X]=\frac{a}{b^2}=\frac{1}{a}\operatorname{E}[X]^2. For the compound Poisson-Gamma variable, we obtain

    E⁡[Y]=E⁡[N]E⁡[X]=λab=μVar⁡[Y]=Var⁡[N]E⁡[X]2+E⁡[N]Var⁡[X]=ϕμpp=a+2a+1∈(1,2)ϕ=(λa)1−p(p−1)b2−p\begin{align*} \operatorname{E}[Y] &= \operatorname{E}[N] \operatorname{E}[X] = \lambda\frac{a}{b}=\mu\\ \operatorname{Var}[Y] &= \operatorname{Var}[N] \operatorname{E}[X]^2 + \operatorname{E}[N] \operatorname{Var}[X] = \phi \mu^p\\ p &= \frac{a + 2}{a+1} \in (1, 2)\\ \phi &= \frac{(\lambda a)^{1-p}}{(p-1)b^{2-p}} \end{align*}

    What’s so special here is that there is a point mass at zero, i.e., P(Y=0)=exp⁡(−μ2−pϕ(2−p))>0P(Y=0)=\exp(-\frac{\mu^{2-p}}{\phi(2-p)}) > 0. Hence, it is a suitable distribution for non-negative quantities with some exact zeros.

    Probability density for compound Poisson Gamma, point masses at zero are marked as points.

    Code

    The rest of this post is about how to compute the density for this parameter range. The easy part is exp⁡(yθ−κ(θ)ϕ)\exp\left(\frac{y\theta - \kappa(\theta)}{\phi}\right) which can be directly implemented. The real obstacle is the term c(y,ϕ)c(y, \phi) which is given by

    c(y,ϕ)=Φ(−α,0,t)yα=2−p1−pt=((p−1)ϕy)α(2−p)ϕ\begin{align*} c(y, \phi) &= \frac{\Phi(-\alpha, 0, t)}{y} \\ \alpha &= \frac{2 - p}{1 - p} \\ t &= \frac{\left(\frac{(p - 1)\phi}{y}\right)^{\alpha}}{(2-p)\phi} \end{align*}

    This depends on Wright’s (generalized Bessel) function Φ(a,b,z)\Phi(a, b, z) as introduced in a 1933 paper by E. Wright.

    Wright’s Generalized Bessel Function

    According to DLMF 10.46, the function is defined as

    Φ(a,b,z)=∑k=0∞zkk!Γ(ak+b),a>−1,b∈R,z∈C\begin{equation*} \Phi(a, b, z) = \sum_{k=0}^{\infty} \frac{z^k}{k!\Gamma(ak+b)}, \quad a > -1, b \in R, z \in C \end{equation*}

    which converges everywhere because it is an entire function. We will focus on the positive real axis z=x≥0z=x\geq 0 and the range a≥0,b≥0a\geq 0, b\geq 0 (note that a=−α∈(0,∞)a=-\alpha \in (0,\infty) for 1<p<21<p<2). For the compound Poisson-Gamma, we even have b=0b=0.

    Implementation of such a function as done in scipy.stats.wright_bessel, even for the restricted parameter range, poses tremendous challenges. The first one is that it has three parameters which is quite a lot. Then the series representation above, for instance, can always be used, but depending on the parameters, it will require a huge amount of terms, particularly for large xx. As each term involves the Gamma function, this becomes expensive very fast. One ends up using different representations and strategies for different parameter regions:

    • Small xx: Taylor series according to definition
    • Small aa: Taylor series in a=0a=0
    • Large xx: Asymptotic series due to Wright (1935)
    • Large aa: Taylor series according to definition for a few terms around the approximate maximum term kmaxk_{max} due to Dunn & Smyth (2005)
    • General: Integral represantation due to Luchko (2008)

    Dunn & Smyth investigated several evaluation strategies for the simpler Tweedie density which amounts to Wright’s functions with b=0b=0, see Dunn & Smyth (2005). Luchko (2008) lists most of the above strategies for the full Wright’s function.

    Note that Dunn & Smyth (2008) provide another strategy to evaluate the Tweedie distribution function by means of the inverse Fourier transform. This does not involve Wright’s function, but also encounters complicated numerical integration of oscillatory functions.

    The Integral Representation

    This brings us deep into complex analysis: We start with Hankel’s contour integral representation of the reciprocal Gamma function.

    1Γ(z)=12πi∫Ha−ζ−zeζ  dζ\begin{equation*} \frac{1}{\Gamma(z)} = \frac{1}{2\pi i} \int_{Ha^-} \zeta^{-z} e^\zeta \; d\zeta \end{equation*}

    with the Hankel path Ha−Ha^- from negative infinity (A) just below the real axis, counter-clockwise with radius ϵ>0\epsilon>0 around the origin and just above the real axis back to minus infinity (D).

    Hankel contour Ha– in the complex plane.

    In principle, one is free to choose any such path with the same start (A) and end point (D) as long as one does not cross the negative real axis. One usually lets the AB and CD be infinitesimal close to the negative real line. Very importantly, the radius ϵ>0\epsilon>0 is a free parameter! That is real magic🪄

    By interchanging sum and integral and using the series of the exponential, Wright’s function becomes

    Φ(a,b,z)=∑k=0∞zkk!12πi∫Ha−ζ−(ak+b)eζ  dζ=12πi∫Ha−ζ−beζ+zζ−a  dζ\begin{align*} \Phi(a, b, z) &= \sum_{k=0}^{\infty} \frac{z^k}{k!} \frac{1}{2\pi i} \int_{Ha^-} \zeta^{-(ak+b)} e^\zeta \; d\zeta \\ &= \frac{1}{2\pi i} \int_{Ha^-} \zeta^{-b} e^{\zeta + z\zeta^{-a}} \; d\zeta \end{align*}

    Now, one needs to do the tedious work and split the integral into the 3 path sections AB, BC, CD. Putting AB and CD together gives an integral over K, the circle BC gives an integral over P:

    Φ(a,b,x)=1π∫ϵ∞K(a,b,x,r)  dr+ϵ1−bπ∫0πP(ϵ,a,b,x,φ)  dφK(a,b,x,r)=r−bexp⁡(−r+xr−acos⁡(πa))sin⁡(x⋅r−asin⁡(πa)+πb)P(ϵ,a,b,x,φ)=exp⁡(ϵcos⁡(φ)+xϵ−acos⁡(aφ))cos⁡(ϵsin⁡(φ)−x⋅ϵ−asin⁡(aφ)+(1−b)φ)\begin{align*} \Phi(a, b, x) &= \frac{1}{\pi} \int_{\epsilon}^\infty K(a, b, x, r) \; dr \\ &+ \frac{\epsilon^{1-b}}{\pi} \int_0^\pi P(\epsilon, a, b, x, \varphi) \; d\varphi \\ K(a, b, x, r) &= r^{-b}\exp(-r + x r^{-a} \cos(\pi a)) \\ &\quad \sin(x \cdot r^{-a} \sin(\pi a) + \pi b) \\ P(\epsilon, a, b, x, \varphi) &= \exp(\epsilon \cos(\varphi) + x \epsilon^{-a}\cos(a \varphi)) \\ &\quad \cos(\epsilon \sin(\varphi) - x \cdot \epsilon^{-a} \sin(a \varphi) + (1-b) \varphi) \end{align*}

    What remains is to carry out the numerical integration, also known as quadrature. While this is an interesting topic in its own, let’s move to the magic part.

    Arc Length Minimization

    If you have come so far and say, wow, puh, uff, crazy, 🤯😱 Just keep on a little bit because here comes the real fun part🤞

    It turns out that most of the time, the integral over P is the most difficult. The worst behaviour an integrand can have is widely oscillatory. Here is one of my favorite examples:

    Integrands for a=5, b=1, x=100 and two choices of epsilon.

    With the naive choice of ϵ=1\epsilon=1, both integrands (blue) are—well—crazy. There is basically no chance the most sophisticated quadrature rule will work. And then look at the other choice of ϵ≈4\epsilon\approx 4. Both curves seem well behaved (for P, we would need a closer look).

    So the idea is to find a good choice of ϵ\epsilon to make P well behaved. Well behaved here means most boring, if possible a straight line. What makes a straight line unique? In flat space, it is the shortest path between two points. Therefore, well behaved integrands have minimal arc length. That is what we want to minimize.

    The arc length SS from x=ax=a to x=bx=b of a 1-dimensional function ff is given by

    S=∫ab1+f′(x)2  dx\begin{equation*} S = \int_a^b \sqrt{1 + f^\prime(x)^2} \; dx \end{equation*}

    Instead of f=Pf=P, we only take the oscillatory part of P and approximate the arc length as f(φ)=f(φ)=ϵsin⁡(φ)−xϵ−ρsin⁡(ρφ)+(1−β)φf(\varphi)=f(\varphi) = \epsilon \sin(\varphi) - x \epsilon^{-\rho} \sin(\rho \varphi) + (1-\beta) \varphi. For a single parameter point a,b,za, b, z this looks like

    Arc length and integrand P for different epsilon, given a=0.1, b=5, x=100.

    Note the logarithmic y-scale for the right plot of P. The optimal ϵ=10\epsilon=10 is plotted in red and behaves clearly better than smaller values of ϵ\epsilon.

    What remains to be done for an actual implementation is

    • Calculate minimal ϵ\epsilon for a large grid of values a,b,xa, b, x.
    • Choose a function with some parameters.
    • Curve fitting (so again optimisation): Fit this function to the minimal ϵ\epsilon of the grid via minimising least squares.
    • Implement some quadrature rules and use this choice of ϵ\epsilon in the hope that it intra- and extrapolates well.

    This strategy turns out to work well in practice and is implemented in scipy. As the parameter space of 3 variables is huge, the integral representation breaks down in certain areas, e.g. huge values of ϵ\epsilon where the integrands just overflow numerically (in 64-bit floating point precision). But we have our other evaluation strategies for that.

    Conclusion

    An extensive notebook for Wright’s function, with all implementation strategies can be found here.

    After an adventurous journey, we arrived at one implementation strategy of Wright’s generalised Bessel function, namely the integral representation. The path went deep into complex analysis and contour integration, then further to the arc length of a function and finally curve fitting via optimisation. I am really astonished how connected all those different areas of mathematics can be.

    Wright’s function is the missing piece to compute full likelihoods and probability functions of the Tweedie distribution family and is now available in the Python ecosystem via scipy.

    We are at the very end of this Tweedie trilogy. I hope it has been entertaining and it has become clear why Tweedie deserves to be celebrated.

    Further references:

  • A Tweedie Trilogy — Part II: Offsets

    TLDR: This second part of the trilogy will have a deeper look at offsets and sample weights of a GLM. Their non-equivalence stems from the mean-variance relationship. This time, we not only have a Poisson frequency but also a Gamma severity model.

    This trilogy celebrates the 40th birthday of Tweedie distributions in 2024 and highlights some of their very special properties.

    See part I.

    From Mean-Variance Relation to Score Equations

    In the part I, we have already introduced the mean-variance relation of a Tweedie random variable Y∼Twp(μ,ϕ)Y\sim Tw_p(\mu, \phi) with Tweedie power pp, mean μ\mu and dispersion parameter ϕ\phi:

    E⁡[Y]=μVar⁡[Y]=ϕμp=ϕv(μ)\begin{align*} \operatorname{E}[Y] &= \mu \\ \operatorname{Var}[Y] &= \phi \mu^p = \phi v(\mu) \end{align*}

    with variance function v(μ)v(\mu).

    This variance function directly impacts the estimation of GLMs. Assume the task is to estimate the expectation of a random variableYi∼Twp(μi,ϕ/wi)Y_i\sim Tw_p(\mu_i, \phi/w_i), given observations of the target yiy_i and of explanatories variables, aka features, xi∈Rkx_i\in R^k. A GLM then assumes a link functiong(μi)=∑j=1kxijβjg(\mu_i) = \sum_{j=1}^k x_{ij}\beta_jwith coefficients β\beta to be estimated via an optimization procedure, of which the first order condition, also called score equation, reads

    ∑iwiyi−μiv(μi)g′(μi)xij=0∀j=1,…,k\begin{equation*} \sum_i w_i \frac{y_i - \mu_i}{v(\mu_i)g'(\mu_i)} x_{ij}= 0 \quad \forall j =1, \ldots, k \end{equation*}

    This shows that the higher the Tweedie power pp, entering via v(μ)v(\mu) only, the less weight is given to deviations of large values. In other words, higher Tweedie powers result in GLMs that are less and less sensitive to what happens at large (expected) values.

    This is also reflected in the deviance loss function. They can be derived from the negative log-likelihood and are given by

    dp(y,μ)=2⋅{max⁡(0,y2−p)(1−p)(2−p)−yμ1−p1−p+μ2−p2−pp∈R∖(0,1]∪{2}ylog⁡yμ−y+μp=0yμ−log⁡yμ−1p=2\begin{equation*} d_p(y, \mu) = 2 \cdot \begin{cases} \frac{\max(0, y^{2-p})}{(1-p)(2-p)}-\frac{y\mu^{1-p}}{1-p}+\frac{\mu^{2-p}}{2-p} & p \in \mathrm{R}\setminus (0,1] \cup \{2\} \\ y\log\frac{y}{\mu} - y + \mu & p=0 \\ \frac{y}{\mu} - \log\frac{y}{\mu} - 1 & p=2 \end{cases} \end{equation*}

    These are the only strictly consistent scoring functions for the expectation (up to one multiplicative and one additive constant) that are homogeneous functions (of degree 2−p2-p), see, e.g., Fissler et al (2022). The Poisson deviance (p=1p=1), for example, has a degree of homogeneity of 1 and the same unit as the target variable. The Gamma deviance (p=2p=2), on the other side, is zero-homogeneous and is completely agnostic to the scale of its arguments. This is another way of stating the above: the higher the Tweedie power the less it cares about large values.

    It is also connected to the fact that Tweedie distributions are the only distributions from the exponential dispersion family that are closed under scale transformations:

    Y∼Twp(μ,ϕ)cY∼Twp(cμ,c2−pϕ)∀c>0\begin{align*} Y &\sim Tw_p(\mu, \phi) \\ cY &\sim Tw_p(c\mu, c^{2-p}\phi) \quad \forall c>0 \end{align*}

    Offsets and Sample Weights

    Poisson GLM

    When estimating counts with a Poisson GLM, there is often an exposure measure like time under consideration or underlying number of things (insurance policies, trees in a forest, radioactive atoms). One then often finds two different, but equivalent formulations of a Poisson GLM with log-link.

    • Sample weights: Model frequency y=Nwy=\frac{N}{w} and fit with sample weights ww to estimate E⁡[y]=μy=exp⁡(xβ)\operatorname{E}[y] = \mu_y = \exp(x \beta).
    • Offsets: Model counts NN, but account for the exposure ww via an offset as E⁡[N]=μN=exp⁡(xβ+log⁡(w))=wμy\operatorname{E}[N]=\mu_N = \exp(x \beta + \log(w)) = w \mu_y.

    Note that each way models a different target, so we had to use subscripts to distinguish the mean parameters μ\mu.

    In this special case of a Poisson GLM with (canonical) log-link, both models are equivalent and will result in the exact same parameters β\beta. You can plug it into the score equation to convince yourself.

    Tweedie GLM

    Very importantly, this simple equivalence of GLM fomulations with offsets and with sample weights does only hold for the Poisson GLM with log-link. It does not hold for any other Tweedie parameter or even other distributions from the exponential dispersion family.

    One can show that a Tweedie GLM with log-link and offset (additive in link space) log⁡(u)\log(u) on target yy with weights ww is equivalent to the same Tweedie GLM but with target yu\frac{y}{u} and weights wu2−pw u^{2-p}.

    So one can construct an equivalence between unweighted with offsets and weighted without offsets by setting u=w2−pu = \sqrt[2-p]{w}. But note that this does not work for a Gamma GLM which as p=2p=2.

    Example

    We continue with the same dataset and model as in part I and show this (non-) equivalence with the offsets.

    Python

    As next, we model the severity Y=lossNY = \frac{loss}{N} with claim counts NN as weights. As is standard, we use a Gamma GLM with log-link (which is not canonical this time).

    Python

    The deviations might seem small, but they are there and add up:

    Python

    Here, it becomes evident that the two models are quite different.

    Outlook

    The full notebook can be found here.

    The final part III of the Tweedie trilogy will follow in one week and go into details of the probability density function.

  • A Tweedie Trilogy — Part I: Frequency and Aggregration Invariance

    TLDR: In this first part of the Tweedie Trilogy, we will take a look at what happens to a GLM if we aggregate the data by a group-by operation. A frequency model for insurance pricing will serve as an example.

    This trilogy celebrates the 40th birthday of Tweedie distributions in 2024 and highlights some of their very special properties.

    Intro

    Tweedie distributions and Generalised Linear Models (GLM) have an intertwined relationship. While GLMs are, in my view, one of the best reference models for estimating expectations, Tweedie distributions lie at the heart of expectation estimation. In fact, basically all applied GLMs in practice use Tweedie distributions with three notable exceptions: the binomial, the multinomial and the negative binomial distribution.

    Mean-Variance Relation

    “An index which distinguishes between some important exponential families” is the original publication title of Maurice Charles Kenneth Tweedie in 1984—but note that Shaul K. Bar-Lev and Peter Enis published around the same time; as their 1986 paper was received November 1983, the distribution could also be named Bar-Lev & Enis distribution.1 This index is meanwhile called the Tweedie power parameter pp. Recall that distributions of the exponential dispersion family always fulfil a mean-variance relationship. Its even a way to define them. For the Tweedie distribution, denoted Twp(μ,ϕ)Tw_p(\mu, \phi), the relation reads

    E⁡[Y]=μVar⁡[Y]=ϕμp\begin{align*} \operatorname{E}[Y] &= \mu \\ \operatorname{Var}[Y] &= \phi \mu^p \end{align*}

    with dispersion parameter ϕ\phi. Some very common members are given in the following table.

    ppdistributiondomain YYdomain μ\mu
    0Normal / GaussianR\mathrm{R}R\mathrm{R}
    1Poisson0,1,2,…0, 1, 2, \ldotsR+\mathrm{R}_+
    (1,2)(1,2)Compound Poisson-GammaR+∪{0}\mathrm{R}_+ \cup \{0\}R+\mathrm{R}_+
    2GammaR+\mathrm{R}_+R+\mathrm{R}_+
    3inverse GaussianR+\mathrm{R}_+R+\mathrm{R}_+

    Insurance Pricing Models

    In non-life insurance pricing, most claims happen somewhat randomly, typically the occurrence as well as the size. Take the theft of your bike or a water damage of your basement due to flooding as an example. Pricing actuaries usually want to predict the expected loss E[Y∣X]E[Y|X] given some features XX of a policy. The set of features could contain the purchasing price of your bike or the proximity of your house to a river.

    Instead of directly modelling the expected loss per exposure ww, e.g. the time duration of the insurance contract, the most used approach is the famous frequency-severity split:

    E⁡[Yw]=E⁡[Nw]⏟frequency⋅E⁡[Yn∣N=n]⏟severity\begin{align*} \operatorname{E}\left[\frac{Y}{w}\right] = \underbrace{\operatorname{E}\left[\frac{N}{w}\right]}_{frequency} \cdot \underbrace{\operatorname{E}\left[\left. \frac{Y}{n}\right| N=n\right]}_{severity} \end{align*}

    For simplicity, the conditioning on XX is suppressed, it would occur in every expectation. The first part E⁡[Nw]\operatorname{E}\left[\frac{N}{w}\right]is the (expected) frequency, i.e. the number of claims per exposure (time). The second term E⁡[YN∣N]\operatorname{E}\left[\left.\frac{Y}{N}\right| N\right] is the (expected) severity, i.e. the average claim size (per claim) given a fixed number of claims. Here, we focus on the frequency part.

    Convolution and Aggregation Invariance

    This property might first seem very theoretical, but it may be one of the most important properties for the estimation of expectations E[Y∣X]E[Y|X] with GLMs. It is in fact a property valid for the whole exponential dispersion family: The weighted mean of i.i.d. random variables has (almost) the same distribution!

    If

    Yi∼i.i.dTwp(μ,ϕ/wi) ,w+=∑iwiwith wi>0 ,\begin{align*} Y_i &\overset{i.i.d}{\sim} \mathrm{Tw}_p(\mu, \phi/w_i) \,, \\ w_+ &= \sum_i w_i \quad\text{with } w_i >0 \,, \end{align*}

    then

    Y=∑inwiYiw+∼Twp(μ,ϕ/w+) .\begin{align*} Y &=\sum_i^n \frac{w_i Y_i}{w_+} \sim \mathrm{Tw}_p(\mu, \phi/w_+) \,. \end{align*}

    It is obvious that the mean of YYis again μ\mu. But is is remarkable that it has the same distribution with the same power parameter, only the 2nd argument with the dispersion parameter differs. But the dispersion parameter cancels out in GLM estimations. In fact, we will show that two GLMs, one on aggregated data, give identical results. Another way of saying the same in statistical terms is that (weighted) averages are the sufficient statistic for the expectation within the exponential dispersion family.

    This is quite an essential property for data aggregation. It means that one can aggregate rows with identical features and still do an analysis (of the conditional expectation) without loss of information.

    The weighted average above can be written a bit more intuitive. For instance, a frequency Yi=NiwiY_i=\frac{N_i}{w_i} has weighted average Y=∑iNi∑iwiY=\frac{\sum_i N_i}{\sum_i w_i}.

    Poisson Distribution

    When modelling counts, the Poisson distribution is by far the easiest distribution one can think of. It only has a single parameter, is a member of the Tweedie family, and fulfils the mean-variance relation

    E⁡[N]=μ=Var⁡[N] .\begin{equation*} \operatorname{E}[N] = \mu = \operatorname{Var}[N] \,.\end{equation*}

    In particular, p=1p=1. While the distribution is strictly speaking only for counts, i.e. NN takes on non-negative integer values, Poisson regression also works for any non-negative response variable like N/w∈RN/w \in \mathrm{R}.

    Frequency Example

    For demonstration, we fit a Poisson GLM on the french motor third-party liability claims dataset, cf. the corresponding scikit-learn example and the case study 1 of the Swiss Association of Actuaries on the same dataset.

    Python

    In fact, both models have the same intercept term and same coefficients, they are really identical models (up to numerical precision):

    Python

    Outlook

    The full notebook can be found here.

    In the next week, part II of this trilogy will follow. There, we will meet some more of its quite remarkable properties.

    Further references:

    • Tweedie M.C.K. 1984. “An index which distinguishes between some important exponential families”. Statistics: Applications and New Directions. Proceedings of the Indian Statistical Institute Golden Jubilee International Conference, Indian Statistical Institute, Cal- cutta, pp. 579–604.
    • Bar-Lev, S.K., Enis, P. Reproducibility in the one-parameter exponential family. Metrika 32, 391–394 (1985). https://doi.org/10.1007/BF01897827
    • Shaul K. Bar-Lev. Peter Enis. “Reproducibility and Natural Exponential Families with Power Variance Functions.” Ann. Statist. 14 (4) 1507 – 1522, December, 1986. https://doi.org/10.1214/aos/1176350173
    1. A great thanks to Prof. Mario Wüthrich for pointing out the references of Bar-Lev and Enis. ↩︎

  • Building Strong GLMs in Python via ML + XAI

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

    Insurance pricing data

    We will use again a synthetic dataset with 1 Mio insurance policies, with reference:

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

    Let’s start by loading and describing the data:

    Python

    Modeling

    1. We fit a tuned Boosted Trees model to model log(E(claim count)) via Poisson deviance loss.
    2. And perform a SHAP analysis to derive insights.
    Python

    Here, we would come to the conclusions:

    1. car_weight and year might be dropped, depending on the specify aim of the model.
    2. Add a regression spline for driver_age.
    3. Add an interaction between car_power and town.

    Build strong GLM

    Let’s build a GLM with these insights. Two important things:

    1. Glum is an extremely powerful GLM implementation that was inspired by a pull request of our Christian Lorentzen.
    2. In the upcoming version 3.0, it adds a formula API based of formulaic, a very performant formula parser. This gives a very easy way to add interaction effects, regression splines, dummy encodings etc.
    Python

    In this relatively simple situation, the mean Poisson deviance of our models are very simlar now:

    Python

    Final words

    • Glum is an extremely powerful GLM implementation – we have only scratched its surface. You can expect more blogposts on Glum…
    • Having a formula interface is especially useful for adding interactions. Fingers crossed that the upcoming version 3.0 will soon be released.
    • Building GLMs via ML + XAI is so smooth, especially when you work with large data. For small data, you need to be careful to not add hidden overfitting to the model.

    Click here for the full Python notebook

  • ML + XAI -> Strong GLM

    My last post 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.

    Insurance pricing data

    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:

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

    R

    Modeling

    1. We fit a naive additive linear GLM and a tuned Boosted Trees model.
    2. We combine the models and specify their predict function.
    R

    Traditional XAI

    Performance

    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!

    R

    Feature importance

    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 and car_weight seem to be negligile.

    R

    Main effects

    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 via spline will greatly improve the performance and usefulness of the GLM.

    R

    Interaction effects

    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 and car_power. The stratified PDP visualizes this interaction. Let’s add a corresponding interaction effect to our GLM later.

    R

    SHAP

    As an elegant alternative to studying feature importance, PDPs and Friedman’s H, we can simply run a SHAP analysis on the LGB model.

    R

    Here, we would come to the same conclusions:

    1. car_weight and year might be dropped.
    2. Add a regression spline for driver_age
    3. Add an interaction between car_power and town.

    Pimp the GLM

    In the final section, we apply the three insights from above with very good results.

    R

    Or even via permutation or kernel SHAP:

    R

    Final words

    • Improving naive GLMs with insights from ML + XAI is fun.
    • In practice, the gap between GLM and a boosted trees model can’t be closed that easily. (The true model behind our synthetic dataset contains a single interaction, unlike real data/models that typically have much more interactions.)
    • {hstats} can work with multiple regression models in parallel. This helps to keep the workflow smooth. Similar for {kernelshap}.
    • A SHAP analysis often brings the same qualitative insights as multiple other XAI tools together.

    The full R script

  • Explain that tidymodels blackbox!

    Let’s explain a {tidymodels} random forest by classic explainability methods (permutation importance, partial dependence plots (PDP), Friedman’s H statistics), and also fancy SHAP.

    Disclaimer: {hstats}, {kernelshap} and {shapviz} are three of my own packages.

    Diabetes data

    We will use the diabetes prediction dataset of Kaggle to model diabetes (yes/no) as a function of six demographic features (age, gender, BMI, hypertension, heart disease, and smoking history). It has 100k rows.

    Note: The data additionally contains the typical diabetes indicators HbA1c level and blood glucose level, but we wont use them to avoid potential causality issues, and to gain insights also for people that do not know these values.

    R
    “yes” proportion of binary variables (including the response)
    Distribution of numeric variables
    Distribution of smoking_history

    Modeling

    Let’s fit a random forest via tidymodels with {ranger} backend.

    We add a predict function pf() that outputs only the probability of the “Yes” class.

    R

    Classic explanation methods

    R

    Feature importance

    Permutation importance measures by how much the average test loss (in our case log loss) increases when a feature is shuffled before calculating the losses. We repeat the process four times and also show standard errors.

    Permutation importance: Age and BMI are the two main risk factors.

    Main effects

    Main effects are estimated by PDP. They show how the average prediction changes with a feature, keeping every other feature fixed. Using a fixed vertical axis helps to grasp the strenght of the effect.

    PDPs: The diabetes risk tends to increase with age, high (and very low) BMI, presence of heart disease/hypertension, and it is a bit lower for females and non-smoker.

    Interaction strength

    Interaction strength can be measured by Friedman’s H statistics, see the earlier blog post. A specific interaction can then be visualized by a stratified PDP.

    Friedman’s H statistics: Left: BMI and age are the two features with clearly strongest interactions. Right: Their pairwise interaction explains about 10% of their joint effect variability.
    Stratified PDP: The strong interaction between age and BMI is clearly visible. A high BMI makes the age effect on diabetes stronger.

    SHAP

    What insights does a SHAP analysis bring?

    We will crunch slow exact permutation SHAP values via kernelshap::permshap(). If we had more features, we could switch to

    • kernelshap::kernelshap()
    • Brandon Greenwell’s {fastshap}, or to the
    • {treeshap} package of my colleages from TU Warsaw.
    R

    SHAP importance

    SHAP importance: On average, the age increases or decreases the diabetes probability by 4.7% etc. In this case, the top three features are the same as in permutation importance.

    SHAP “summary” plot

    SHAP “summary” plot: Additionally to the bar plot, we see that higher age, higher BMI, hypertension, smoking, males, and having a heart disease are associated with higher diabetes risk.

    SHAP dependence plots

    SHAP dependence plots: We see similar shapes as in the PDPs. Thanks to the vertical scatter, we can, e.g., spot that the BMI effect strongly depends on the age. As in the PDPs, we have selected a common vertical scale to also see the effect strength.

    Final words

    • {hstats}, {kernelshap} and {shapviz} can explain any model with XAI methods like permutation importance, PDPs, Friedman’s H, and SHAP. This, obviously, also includes models developed with {tidymodels}.
    • They would actually even work for multi-output models, e.g., classification with more than two categories.
    • Studying a blackbox with XAI methods is always worth the effort, even if the methods have their issues. I.e., an imperfect explanation is still better than no explanation.
    • Model-agnostic SHAP takes a little bit of time, but it is usually worth the effort.

    The full R script