How Interpretable and Trustworthy are GAMs?

Introduction

As the impact of machine learning on our daily lives continues to grow, we have begun to require that ML systems used for high-stakes decisions (e.g., in healthcare, finance and criminal justice) not only be accurate but also satisfy other properties such as fairness or interpretability (Doshi-Velez and Kim, 2017; Lipton, 2018). Generalized additive models (GAMs) have emerged as a leading model class that is designed to be accurate, and yet simple enough for humans to understand and mentally simulate how a GAM model works (Hegselmann et al., 2020), and is widely used in scientific data exploration (Hastie and Tibshirani, 1995; Pedersen et al., 2019; Izadi, 2020) and model bias discovery (Tan et al., 2018a, b).

GAMs were originally trained using smoothing splines (Hastie and Tibshirani, 1990; Wahba, 1990) that enforced smoothness in the learned functions. Later, several trend-filtering based methods including fused lasso additive models were proposed to make learned functions more sparse and jumpy (Tibshirani et al., 2005; Sadhanala and Tibshirani, 2019). Lou et al. (2012) also proposed using boosted-tree-based methods to fit GAMs. Subsequent work showed the value of tree-based GAMs on two healthcare datasets (Caruana et al., 2015), and also to help audit black-box models to ensure fairness (Tan et al., 2018a).

Do GAMs trained with different algorithms agree with each other? In Fig. 1, we show that three GAMs with similar accuracy provide very different interpretations of the COMPAS recidivism dataset, a dataset in which bias is an important concern. For instance, in Fig. 2(a) EBM-BF suggests that there is no racial bias in the data, while Spline indicates that there is strong racial bias, yet is slightly less accurate. Should we believe EBM-BF because of its slightly higher accuracy and believe there is no racial bias? Probably not. Then how should we determine which GAM to believe?

Second, we examine how much we can trust each GAM to reflect true patterns in the data, a property we call data fidelity. Shmueli (2010) contrasted predictive models that seek to minimize the combination of variance and bias (defined in a statistical sense, not in the sense of unfairness) to explanatory models that aim to capture true patterns in data by minimizing bias alone. For the former, bias can be sacrified for improved variance, and Shmueli (2010) provided examples of how the “wrong” model can sometimes predict better than the right one. In this paper, we study this phenomenon across different GAM algorithms. For real data where we do not know the underlying data patterns, we use the bias term from bias-variance analysis as a proxy for data fidelity. We also experiment with simulated datasets that have different data generators, each of which may favor GAM algorithms with certain inductive biases, and measure the worst-case data fidelity of each GAM algorithm across multiple datasets. This allows us to quantify if some GAM algorithms have high accuracy but low data fidelity, which may mislead users to trust the wrong explanations.

We compare different GAM algorithms on ten classification datasets and find that the most accurate GAMs yield similar accuracy, yet learn qualitatively different explanations.

We measure which GAM algorithms lead to models that are more or less sparse, a property we call feature sparsity. We show that sparse-feature GAMs can discriminate on rare subpopulations leading to unfairness.

We examine several case studies of data anomaly discovery to see which GAMs can or cannot be trusted to discover these true patterns in the data, a property we call data fidelity.

We show that some GAMs have high accuracy but low data fidelity which can mislead users who select models by accuracy alone.

We find that inductive bias plays a crucial role in model explanations, and recommend tree-based GAMs over other GAMs for their low feature sparsity and superior data fidelity.

Related Work

While we study GAMs in this paper, GAMs are not the only interpretable model class to come under scrutiny recently. The instability of decision trees (another model class commonly considered interpretable) has been pointed out (Dwyer and Holte, 2007), and the vulnerability of post-hoc explanation methods such as LIME (Ribeiro et al., 2016) and Shapley values (Lundberg and Lee, 2017) to input perturbation has been exploited to generate adversarial attacks on model explanations (Slack et al., 2020). Hooker and Mentch (2019) also found partial dependence and feature importance metrics based on permuting inputs to be particularly misleading when inputs are highly dependent, and earlier work found feature importance metrics to be biased for certain types of models with different inductive biases, e.g., random forest feature importance is biased towards variables with many potential splits such as categorical variables with many levels (Strobl et al., 2007; Zhou and Hooker, 2021).

Our paper is not the first to compare different GAM algorithms, but to the best of our knowledge it is the first to focus on interpretability and its relationships to fairness on different GAM algorithms. Binder and Tutz (2008) compared three different spline training algorithms, including backfitting, joint optimization, and boosting, finding that boosting performed particularly well in high-dimensional settings. Lou et al. (2012) also found that boosted shallow bagged trees yielded higher accuracy than other GAM algorithms. However both papers focused on accuracy, not interpretability.

Methods

In this section we describe the different GAM algorithms used in this paper. To make it easier for readers, we defer the description of the new metrics we define in this paper – feature sparsity and data fidelity – to just before their use in Sec. 5.

GAMs are interpretable because the impact of each feature, $f_{j}$, on the prediction can be visualized as a graph (see Fig. 2 for an example), and humans can easily simulate how a GAM works by reading $f_{j}$s off different features from the graph and adding them together. We select the following six GAM algorithms to compare in this paper based on their popularity, state-of-the-art performance and availability of open source implementations.

Explainable Boosting Machine (EBM) A tree-based GAM designed for intelligibility and high accuracy (Lou et al., 2012; Caruana et al., 2015; Nori et al., 2019) where shape functions $f_{j}$ are gradient-boosted ensembles of bagged trees. Each tree operates on a single variable, preventing interactions effects from being learned. Trees are grown by repeatedly cycling through features, which forces the model to sequentially consider each feature as an explanation of the current residual rather than greedily selecting the best feature. This deliberate construction makes this model have less feature sparsity. For comparison, we create a sparse version of EBM similar to regular gradient boosted trees, ”EBM-BF” (EBM-BestFirst), that greedily grows the next tree on the best, most informative feature to reduce error as much as possible at each step. Like most gradient boosted trees, EBM-BF is likely to put most weight on a few very important features, modest weight on a larger number of moderately useful features, and little or no weight on features whose signal could be learned by other stronger, correlated features.

Spline A classic way to train GAMs is with spline basis functions (Hastie and Tibshirani, 1990). We tried a variety of spline methods in 2 popular packages, the Python pygam (Servén and Brummitt, 2018) and R mgcv package (Wood, 2011), and chose cubic splines in pygam because it has a good combination of accuracy, robustness and speed (Fig. 2(c)).

Logistic regression (LR) and other strawmen approaches We compare these other approaches to Logistic Regression (LR), a widely used linear model that cannot learn non-linear shape plots. We also compare to two other strawmen approaches: marginalized LR (mLR) and indicator LR (iLR). We first bin each feature $x_{j}$ into at most $255$ bins. In contrast to LR that assumes $f_{j}(x_{j})=w_{j}x_{j}$, mLR sets $f_{j}(x_{j})=w_{j}g(x_{j})$ where $g(x_{j})$ is the average (marginalized) value of target $y$ within the same bin as $x_{j}$ in the dataset. This is a GAM model built by applying logistic regression on top of marginalization, thus preventing shape plots from being learned in concert with each other. iLR treats each bin as a new feature (similar to one-hot encoding) and learns an LR on the transformed features. It thus ignores proximity relationships between different feature values (Fig. 2(e)).

2. Training and Hyperparameters

To fairly compare different GAM algorithms, we choose hyperparameters that perform best for each algorithm individually. Below, we briefly mention how we tune each GAM algorithm, and point the reader to the Appendix for more details.

We split each dataset into 70-15-15% train-val-test splits and repeat our training procedure run $5$ times. This allowed us to derive uncertainty estimates in the form of standard deviation across multiple runs.

Case studies: COMPAS, Adult, MIMIC-II

We start with some case studies to highlight the implications of different GAM algorithms on common interpretability tasks such as surfacing unfairness or discovering anomalies in data. In this section, we highlight our key findings with plots specifically picked to be representative of our main results. A complete set of plots can be found in Appendix B.

One key property we study in this paper is which GAM algorithm uses fewer features to make predictions i.e. feature sparsity. Although sparsity is sometimes preferred because it appears to generate simpler explanations, it can hide data bias and discriminate against minority groups. Here we examine the sparsity properties of different GAM algorithms on two datasets that have been studied in the fairness community for racial and gender bias (Zemel et al., 2013; Chouldechova, 2017; Mehrabi et al., 2019). The COMPAS dataset contains demographic, crime, and recidivism information for defendants in Broward County, Florida, in $2013$ and $2014$. Research has suggested that the COMPAS recidivism risk score may be racially biased (Angwin et al., 2019). The Adult dataset extracted demographic information, including age, race, occupation, sex, etc. from the 1994 census data to predict if an individual’s income exceeds 50k/yr. In the dataset, males have on average higher annual incomes than females (Mehrabi et al., 2019).

To motivate our analysis, we compare two GAM algorithms that are very different from each other in terms of feature sparsity: sparse EBM-BF and regular, “dense” EBM, yet achieve similar accuracy (see Table 4). Figure 3 displays the shape plots on two sensitive attributes, race and gender, on the COMPAS dataset. Since these features have modest influence compared to other features, the sparse-feature EBM-BF shows no or only a tiny effect on these sensitive attributes, while EBM shows much larger effects. Although there is no easy way to judge which GAM is more “causally” correct, the sparse EBM-BF makes users unaware of bias that may exist in the data and has been learned by other stronger, correlated features. In contrast, the dense EBM shows effects on all features. Because of this, we suggest that the dense model is better suited for surfacing potential bias in data then can then be investigated further by humans.

Next, we investigate how feature sparsity affects minority groups. Table 1 presents the predictive performance (cross entropy loss) of EBM and EBM-BF on each minority group. Although EBM and EBM-BF have negligible difference (less than $0.5\%$) in terms of overall loss, compared to EBM, EBM-BF exhibits greater loss on minority groups Other ($1.45\%$) and Asian ($6\%$) compared to majority group White ($-0.02\%$); EBM exhibits lower loss on the Native American group ($-2.26\%$). To further investigate this phenomenon, we perform an ablation study by removing the race feature from EBM thus forcing EBM to be more sparse. While this increased overall loss by $0.1\%$ compared to EBM with the race feature, the loss for minority groups was again substantially increased, with the loss increasing by $6\%$ for Asian and $1\%$ for Native American. Similarly, when we remove the sex feature from EBM, the loss for the minority group Female increased by $0.99\%$, almost four times larger than the overall loss increase ($0.23\%$). Unexpectedly, removing the sex feature improves the loss for minority group Native Americans ($-5.32\%$); this is a possible explanation for why the loss for Native Americans is smaller for EBM-BF than EBM, as EBM-BF placed little importance on sex.

We repeat the same analysis on the Adult dataset. Table 2 presents the loss of EBM and EBM-BF on each minority group in the Adult dataset. Compared to EBM, EBM-BF exhibit greater loss on minority groups Indian ($7.13\%$) and Other ($19.05\%$), much more than the overall loss ($5.27\%$) or loss on majority group White ($5.17\%$). We also find that removing race from the EBM model increased the loss more for minority groups Indian ($5.61\%$) and Other ($1.04\%$), and removing sex from the EBM model increases the loss for Female ($5.78\%$) much more than for Male ($1.54\%$).

Implications GAM algorithms with a tendency to use fewer features to make predictions (e.g. EBM-BF) showed only small effects on sensitive attributes and exhibited greater prediction loss on minority groups causing unfairness, compared to GAM algorithms that tend to use more features to make predictions (e.g. EBM).

2. Data Anomaly Discovery

Another key property we study in this paper is which GAM algorithm is better able to capture anomalies in data. To illustrate, we train different GAM algorithms on a medical dataset: ICU mortality prediction dataset MIMIC-II (Johnson et al., 2016). On this dataset, XGB and EBM have similar shape plots thus we only present the EBM plots here for simplicity.

Fig. 4(a) displays one feature, PFratio (a measure of how well patients convert oxygen in air to oxygen in blood), for the three most accurate GAM algorithms on this dataset: EBM, Spline and FLAM. Interestingly, both EBM and FLAM capture a sharp drop in mortality risk at PFratio=$332$. It turns out that PFratio is usually not measured for healthier patients, and the missing values for these patients have been imputed by its population mean $332$ (a common preprocessing fix for missing data), thus giving a group of low-risk patients the mean value of this feature. However, Spline is unable to represent the sharp drop, becoming distorted in the region 300-600, thereby underestimating the risk for patients in this region.

Fig. 4(b) for Systolic Blood Pressure (BP) shows another data anomaly that is only captured by tree-based GAM algorithms EBM and XGB. EBM captures three jumps, exhibiting dips in risk predictions near $175$, $200$ and $225$. These are likely to be human intervention artifacts, since $175$, $200$, and $225$ are treatment thresholds used by physicians. As a patient’s Systolic BP increases the mortality risk naturally increases, but when they reach the next treatment threshold, risk actually drops because most patients just above the threshold are receiving more aggressive care that is effective at reducing their risk. Both Spline and FLAM are too smooth or flat and fail to capture these anomalies.

Implications Localized data anomalies such as mean imputation and human intervention artifacts (e.g. medical treatment thresholds), often require models to learn quick, non-linear changes in risk. Tree-based methods (e.g. EBM and XGB) can detect these much better compared to GAM algorithms that are too smooth or sparse (e.g. Spline and FLAM).

Quantitative Analysis of GAMs

In the previous section, we saw examples of how different GAM algorithms revealed different insights. In this section, we study the performance differences between GAM algorithms quantitatively. We first benchmark the test accuracy of different GAMs on ten different datasets (Sec. 5.1). Then we measure feature sparsity of different GAM algorithms (Sec. 5.2). Finally, we measure data fidelity using both real (Sec. 5.3) and simulated data (Sec. 5.4, 5.5).

How do we choose which GAM to use? Accuracy is perhaps the first obvious consideration. Table 4 provides test set AUC of different GAM algorithms on ten datasets. These datasets of varying size (500 - 250k samples) and number of features ($6$ - $57$ features) span different domains such as healthcare, criminal justice, finance, and retail (see Table 3). In addition to the nine GAM algorithms described in Sec. 3, we also include two full-complexity methods: Random Forest (RF) and XGB with depth $3$ (XGB-d3). For each method, we compute three metrics, each of which is averaged over ten datasets: (1) Test AUC; (2) Rank of test AUC compared to other methods (lower rank is better); (3) Test AUC normalized compared to other methods (lowest test AUC for a dataset has value 0, highest test AUC for a dataset has value 1, with all other test AUCs scaled linearly between them). On average across ten datasets, EBM, EBM-BF, and XGB-d3 performed the best. In general, GAMs perform better than or comparably to full complexity models. Four of the GAMs (EBM, XGB, Spline and FLAM) achieve similar top performance with average AUC differences less than $0.2\%$.

Implications There exist GAM algorithms that perform comparably to full complexity models. Several GAM algorithms are similarly accurate, hence accuracy should not be the sole consideration when selecting between different GAM algorithms.

2. GAM feature sparsity

In this section we propose a new metric to quantify feature sparsity, the notion that some GAM algorithms use fewer features than others to make predictions, which we have seen in Sec. 4.1 to impact bias discovery.

Feature density metric The idea is to quantify how fast the test error of a trained model decays (i.e., how fast the model becomes more accurate) as we allow the model to have access to more features; a sparse model only requires a few important features to quickly reduce its test error, while a dense model needs more features to recover because it will have spread learned effects across more of the features. Using the GAM formulation as in Equation 1, we proceed as follows to compute this metric: first we keep only $f_{0}$ and measure the GAM’s test set error as the initial error. Then for each step out of $D$ steps, we greedily search over which feature $f_{j}(x_{j})$, when added back to the model, reduces its validation error the most. We add that feature back and measure how the model’s test error decreases. We save the test error as each subsequent feature is added, until $D$ features are added after $D$ steps, and plot test error against features. Finally we compute the feature density metric as the normalized area under this curve, treating the initial test error as $100$ and final error (with $D$ features) as . We expect an extremely sparse model to have value close to , and a dense model to have value close to $50$.

Implications The proposed feature density metric captures expected behavior. We see lower feature density for methods that greedily select the next best feature (e.g. EBM-BF) or have penalties that regularize for sparsity (e.g. FLAM). Methods that repeatedly cycle over all features (e.g. EBM) have higher feature density.

3. GAM data fidelity

In this section we propose a new metric to quantify how well a GAM is able to capture underlying data patterns, which we have seen in Sec. 4.2 to impact data anomaly discovery.

At first glance, one may think that test accuracy is a suitable metric for this purpose, since it captures how well a model generalizes to unseen data. However, we saw in Sec. 4 when comparing GAM algorithms of similar test accuracy how some were less able to represent certain data patterns. For example, smooth basis functions in Spline, while reducing variance and hopefully improving test set generalization, limited the model’s ability to capture sharp jumps in the data. As noted by Shmueli (2010), some highly accurate predictive models may actually be “wrong” in terms of capturing underlying data patterns. This notion is exactly statistical bias, which arises from model misspecification of the underlying data patterns (Hastie et al., 2009).

Data fidelity metric We use an approximation to the bias term in a bias-variance analysis to measure data fidelity. In bias-variance analysis (Bauer and Kohavi, 1999), the loss of model is composed of noise $N(x)$, bias $B(x)$ and variance $V(x)$ terms:

where $D$ is the training distribution, $t$ is the true label, $y_{*}$ is the optimal predictions, $y_{m}$ is the mean prediction of models across possible training datasets, and $y$ is the model. Since we do not know the $y_{*}$, we instead measure the empirical bias combining both noise and bias $N(x)+B(x)=E_{t}[L\left(t,y_{m}\right)]$ following Munson and Caruana (2009). We use the following sampling procedure: in each round, we split our dataset into 85-15% train-test splits. We then randomly subsample the training data to $50\%$ and train models $5$ times, and we set the average of $5$ models as $y_{m}$ to calculate empirical bias and variance once. Finally, the bias and variance estimates are averaged over eight rounds, and ranked compared to other GAM algorithms on each dataset. We take the average ranks across the ten datasets (lower rank is better).

Fig. 5 plots average variance rank vs. average bias rank for different GAM algorithms. Considering GAM algorithms closest to the bottom left corner (i.e. (0, 0) point), which are also the most accurate GAMs (see Table 4), XGB has the highest data fidelity (lowest bias rank) but has rather high variance. FLAM has the next highest data fidelity, but has even higher variance, hence it is dominated by XGB that has both higher data fidelity and lower variance. After FLAM, XGB-L2, Spline, and EBM have the next highest data fidelity, and promisingly, with significantly lower variance than FLAM or XGB.

Implications We use statistical bias as a proxy to measure data fidelity with real data. By decomposing error into bias and variance components, we see that equally accurate GAM algorithms achieve the same accuracy in different ways. Certain GAM algorithms (e.g. XGB) have lower bias which indicates better fidelity, while other GAMs (e.g. XGB-L2) have lower variance at the expense of higher bias.

4. GAM data fidelity and generator bias

We have thus far studied the data fidelity properties of different GAM algorithms on several real datasets. However, it may be that the inductive bias of a certain GAM algorithm happened to agree with the (unknown) data pattern in a particular real dataset. In this section, we experiment with semi-synthetic datasets created using known data generators. To preserve the character of real datasets as much as possible, we keep the features $X$ but change the label $y$ by training multiple ground truth GAM models (EBM, XGB, Spline, FLAM and LR) on features $X$ and then re-generating the label $y$ as each model’s predictions. Since these GAM models (except LR) are among the most accurate models on most datasets (Table 4), the generated labels capture the real-world distribution as close as possible. As these GAM algorithms are very different from each other, this should provide a diversity of ground truth data patterns.

Fig. 6(a)-(d) shows different GAMs alongside ground truth patterns from two very different generators, Spline and FLAM, on MIMIC-II for one continuous feature (Systolic BP) and one boolean feature (AIDS). Purple represents ground truth, i.e. Spline generator for Fig. 6(a) and (c), and FLAM generator for 6(b) and (d). We see an obvious generator bias: a GAM algorithm fits the ground truth better when ground truth is generated using the same algorithm. For example, on Systolic BP, the Spline GAM fits well the data generated by its own generator (Fig. 6(a)), while doing poorly for data generated by the FLAM generator (Fig. 6(b)), and vice versa for FLAM. However, tree-based methods (EBM, XGB) on Systolic BP with the Spline generator (Fig. 6a) still learn abrupt jumps at $225$ even when the underlying ground truth is smooth; similarly there is also a drop at $175$. This illustrates that it is possible for model inductive bias to dominate irrespective of the true data generator.

To mitigate the aforementioned generator bias, we perform a worst-case analysis: what is the worst performance each GAM algorithm would get across all of the different data generators? Since we do not know the underlying generators on real datasets – they could be jumpy, smooth, or even linear – this analysis is more realistic and robust to all these cases.

Worst-case data fidelity metric To measure how well a GAM can recover the ground truth generators, we calculate the mean absolute difference of each shape plot between the ground truth GAM and the GAM model. Specifically, using the GAM formulation as in Equation 1 where $f_{j}$ is the shape function for feature $j$, and taking $g_{j}$ to be the shape function for the ground truth GAM, we calculate the absolute difference $\sum_{j=1}^{D}|f_{j}(x_{j})-g_{j}(x_{j})|$ across the whole dataset. To compare between datasets, we linearly scale the absolute difference between 0 and 100 for a particular semi-synthetic dataset, with the worst GAM algorithm having value 0 and best GAM algorithm having value 100. We then take the worst score over the five different data generators that yielded five semi-synthetic datasets from each real dataset.

Table 6 provides the worst-case data fidelity for eight GAM algorithms on six real datasets, where each dataset (row) encapsulates five semi-synthetic datasets from different data generators. FLAM and XGB performed the best, then EBM and Spline.

Implications FLAM and XGB exhibit the best worst-case data fidelity. Spline and EBM are similar, and EBM-BF is the worst. Taking into account different data generators, our results are not substantively different from the results derived from the bias-variance analysis on real data in Sec. 5.3.

5. GAM accuracy vs. data fidelity

A GAM model that has high accuracy but low data fidelity may mislead users who tend to judge models solely based on accuracy. We quantify which GAM algorithm is more likely to mislead users this way, by comparing the difference between test AUC rank and data fidelity rank. For each dataset, we compute these two ranks as in Sec. 5.1 and Sec. 5.3, with lower rank being better. Then we take the rank of fidelity minus the rank of test AUC. If the result is negative, we clip it at . We call this the “positive difference” between the two ranks. Finally, we average this over all thirty semi-synthetic datasets. We expect a misleading model to have a lower test AUC rank and higher data fidelity rank.

From Table 7, Spline has the largest difference in rank over multiple datasets with different data generators. This rank difference is largest when the data generators are jumpy, which creates challenges for Spline which uses smooth basis functions.

Implications For Spline, using high test accuracy alone to select a model may be misleading, especially when the underlying data pattern may be jumpy. Other methods are more stable.

Discussion

GAMs are widely used to discover patterns in data in a variety of fields including business (Sapra, 2013), healthcare (Hunter and Prüss-Ustün, 2016), ecology (Pedersen et al., 2019), horticulture (Saw et al., 2017), air pollution (Ravindra et al., 2019), nutrition (Rostami et al., 2020) and COVID-19 (Izadi, 2020). But most of these research only experimented with a specific GAM algorithm (typically Spline) without any comparison to other GAM algorithms. In this work, we have shown that the patterns learned by GAMs are highly impacted by their own inductive biases. If the papers that used GAMs to discover patterns had used different GAM algorithms, would they have drawn different conclusions? How many of the findings are due to true patterns in the data and not due to the inductive bias of the particular GAM algorithm chosen?

While we aimed to provide a useful and fair experimental study, there are limitations to the conclusions that can be drawn from our work due to design choices we made. In terms of data sets, we considered common Kaggle datasets across several domains that are relatively large but still have a manageable amount of features. We do not explore small datasets used in the Spline literature, where a smoothing prior might help compensate for a lack of sample size. In terms of models, we only focused on a few of the most representative GAM algrithms and make additional modifications to these methods to study different characteristics of GAMs (e.g. feature sparsity and data fidelity). We leave more theoretical comparisons to future work.

Conclusion

The key findings are summarized in Table 8, where we have synthesized our findings across six different properties studied in this paper and ranked each GAM algorithm for each property (ties count for half a rank). Although a number of GAM algorithms yield similar accuracy, tree-based methods like EBM and XGB are superior when considering issues such as bias and data anomaly discovery, sparsity, fidelity, and accuracy. Tree-based methods such as XGB and EBM have higher feature density than FLAM or Spline. They also have less bias on real data, and recover data patterns with better fidelity on semi-synthetic data. We also find Spline could have high accuracy yet at the same time low data fidelity, which might mislead users who perform model selection based on test accuracy alone. Qualitatively, Spline and FLAM are not good at detecting local anomalies such as mean imputation or treatment effects, both of which are easily detected by the tree-based methods. Spline also extrapolates over-confidently in low-sample regions (Fig. 1(b), and see other examples in Appendix B).

Future development of better GAM algorithms should focus on the following: (1) GAMs that can better capture rapid non-linear change, (2) GAMs with high feature density to improve fairness and prevent bias masking, (3) GAMs having higher data fidelity on both real and simulated data. We believe our work is an important step towards making GAMs more trustworthy, and our evaluation framework will promote the development of better GAMs in the future.

References

Appendix A Reproducibility: Training details, hyperparameters, and datasets

Code can be found at https://github.com/zzzace2000/GAMs.

In this section, we further describe training details and hyperparameters to supplement the discussion in Sec. 3.

EBM, EBM-BF: we use the open-source package from https://github.com/interpretml/interpret. We set the parameters inner bagging as $100$ and outer bagging as $100$. We find that increasing the number of bags does not further improve performance. We use the default learning rate of $0.01$, default early stopping patience set to $50$, and the maximum $30000$ episodes to make sure it converges.

XGB, XGB-d3, XGB-L2: we use the open source package https://xgboost.readthedocs.io/en/latest/index.html. We also use the default learning rate with the same early stopping patience set as $50$ and number of trees as maximum $30,000$. We use bagging of $100$ times and depth $1$ for our XGB GAM model. For XGB-d3 (XGB with tree depth $3$), we find that bagging of XGB-d3 hurts the performance a bit, and thus do not apply any bagging for XGB-d3. For XGB-L2, we set the parameter ”colsample_bytree” as a small value 1e-5 to make sure each tree only sees one feature.

FLAM: we use the package from R https://cran.r-project.org/web/packages/flam/flam.pdf. We use a $15\%$ validation set to select the best $\lambda$ penalty parameter in the fused LASSO, and then refit the whole data with the best penalty parameter. We set the parameter number of lambda as $100$ and the minimum ratio as 1e-4 to increase the performance of the model.

Spline: we use the pygam package (Servén and Brummitt, 2018). We set the number of basis functions to be $50$ and the maximum iteration as $500$. We find increasing number of basis functions more than $50$ would result in instability when fitting in large datasets.

iLR, mLR: we use the EBM package’s preprocessor to quantily bin the features into $255$ bins. Then we use LR on top of it to train a linear model.

We also tried the following GAM algorithms but do not include them in the main results, for reasons detailed below:

SKGBT: we try the gradient boosting tree in scikit-learn also with tree depth set as $1$. The result is similar to EBM so we do not compare them in the paper.

Cubic spline and plate spline in R mgcv package: to our surprise, mgcv is really unstable on two datasets, Breast Cancer and Churn. After some investigation, we found a possible reason to be that mgcv does not handle numerical instability when the prediction is too close to or $1$.

A.2. Encoding categorical features

For datasets with categorical variables, the choice of encoding can affect both the shape plots and the accuracy. For gradient boosting trees, one may think that using label encoding (LE) is better than one-hot encoding, as one-hot encoding has been shown to have inferior performance in ensemble trees (Wright and König, 2019). We investigate the effects of two types of encoding on EBM and XGB. In $6$ of the datasets with categorical features, EBM with label encoding (LE) indeed shows superior performance to one-hot encoding. However, for XGB, one-hot encoding performs slightly better on average. Thus we use LE for EBM and one-hot encoding for XGB. For the rest of the methods, we use LE for mLR and one-hot encoding for FLAM, Spline, LR and iLR as these methods cannot handle inadequate numerical ordering.

A.3. Dataset sources

The datasets used in this paper can be found at:

Credit: https://www.kaggle.com/mlg-ulb/creditcardfraud

Churn: https://www.kaggle.com/blastchar/telco-customer-churn

COMPAS: https://www.kaggle.com/danofer/compass

MIMIC-II and MIMIC-III dataset (Johnson et al., 2016)

Pneumonia: we thank the authors of Caruana et al. (2015) for running our code on their dataset.

Support2: http://biostat.mc.vanderbilt.edu/DataSets

Appendix B Additional shape plots

The complete set of shape plots can be found at https://drive.google.com/file/d/1PoMRgfuHYax6xuCVU0Dbut3yFJ2ohuLX/view?usp=sharing.