Year: 2021

Abstract:
Deployment of machine learning models in real high-risk settings (e.g. healthcare) often depends not only on model’s accuracy but also on its fairness, robustness and interpretability. Generalized Additive Models (GAMs) are a class of interpretable models with a long history of use in these high-risk domains, but they lack desirable features of deep learning such as differentiability and scalability. In this work, we propose a neural GAM (NODE-GAM) and neural GAM (NODE-GAM) that scale well and perform better than other GAMs on large datasets, while remaining interpretable compared to other ensemble and deep learning models. We demonstrate that our models find interesting patterns in the data. Lastly, we show that we improve model accuracy via self-supervised pre-training, an improvement that is not possible for non-differentiable GAMs.

text: Deployment of machine learning models in real high-risk settings (e.g. healthcare) often depends not only on model’s accuracy but also on its fairness, robustness and interpretability. Generalized Additive Models (GAMs) are a class of interpretable models with a long history of use in these high-risk domains, but they lack desirable features of deep learning such as differentiability and scalability. In this work, we propose a neural GAM (NODE-GAM) and neural GA2M (NODE-GA2M) that scale well and perform better than other GAMs on large datasets, while remaining interpretable compared to other ensemble and deep learning models. We demonstrate that our models find interesting patterns in the data. Lastly, we show that we improve model accuracy via self-supervised pre-training, an improvement that is not possible for non-differentiable GAMs.

heading: 1 INTRODUCTION
text: As machine learning models become increasingly adopted in everyday life, we begin to require models to not just be accurate, but also satisfy other constraints such as fairness, bias discovery, and robustness under distribution shifts for high-stakes decisions (e.g., in healthcare, finance and criminal justice). These needs call for an easier ability to inspect and understand a model’s predictions.
Generalized Additive Models (GAMs) (Hastie & Tibshirani, 1990) have a long history of being used to detect and understand tabular data patterns in a variety of fields including medicine (Hastie & Tibshirani, 1995; Izadi, 2020), business (Sapra, 2013) and ecology (Pedersen et al., 2019). Recently proposed tree-based GAMs and GA2Ms models (Lou et al., 2013) further improve on original GAMs (Spline) having higher accuracy and better ability to discover data patterns (Caruana et al., 2015). These models are increasingly used to detect dataset bias (Chang et al., 2021) or audit black-box models (Tan et al., 2018a;b). As a powerful class of models, they still lack some desirable features of deep learning that made these models popular and effective, such as differentiability and scalability.
In this work, we propose a deep learning version of GAM and GA2M that enjoy the benefits of both worlds. Our models are comparable to other deep learning approaches in performance on tabular data while remaining interpretable. Compared to other GAMs, our models can be optimized using GPUs and mini-batch training allowing for higher accuracy and more effective scaling on larger datasets. We also show that our models improve performance when labeled data is limited by self-supervised pretraining and finetuning, where other non-differentiable GAMs cannot be applied.
Several works have focused on building interpretable deep learning models that are effective for tabular data. TabNet (Arik & Pfister, 2020) achieves state-of-the-art performance on tabular data while also providing feature importance per example by its attention mechanism. Although attention seems to be correlated with input importance (Xu et al., 2015), in the worst case they might not correlate well (Wiegreffe & Pinter, 2019). Yoon et al. (2020) proposes to use self-supervised learning on tabular data and achieves state-of-the-art performance but does not address interpretability. NIT (Tsang et al., 2018) focuses on building a neural network that produces at most K-order interactions and thus include GAM and GA2M. However, NIT requires a two-stage iterative training process that requires longer computations. And their performance is slightly lower to DNNs while ours are overall on par with it. They also do not perform purification that makes GA2M graphs unique when showing them.
The most relevant approaches to our work are NODE (Popov et al., 2019) and NAM (Agarwal et al., 2020). Popov et al. (2019) developed NODE that mimics an ensemble of decision trees but permits differentiability and achieves state-of-the-art performance on tabular data. Unfortunately, NODE suffers from a lack of interpretability similarly to other ensemble and deep learning models. On the other hand, Neural Additive Model (NAM) whose deep learning architecture is a GAM, similar to our proposal, thus assuring interpretability. However, NAM can not model the pairwise interactions and thus do not allow GA2M. Also, because NAM builds a small feedforward net per feature, in high-dimensional datasets NAM may require large memory and computation. Finally, NAM requires training of 10s to 100s of models and ensemble them which incurs large computations and memory, while ours only trains once; our model is also better than NAM without the ensemble (Supp. A).
To make our deep GAM scalable and effective, we modify NODE architecture (Popov et al., 2019) to be a GAM and GA2M, since NODE achieves state-of-the-art performance on tabular data, and its tree-like nature allows GAM to learn quick, non-linear jumps that better match patterns seen in real data (Chang et al., 2021). We thus call our models NODE-GAM and NODE-GA2M respectively.
One of our key contributions is that we design several novel gating mechanisms that gradually reduce higher-order feature interactions learned in the representation. This also enables our NODE-GAM and NODE-GA2M to automatically perform feature selection via back-propagation for both marginal and pairwise features. This is a substantial improvement on tree-based GA2M that requires an additional algorithm to select which set of pairwise feature interactions to learn (Lou et al., 2013).
Overall, our contributions can be summarized as follows:
• Novel architectures for neural GAM and GA2M thus creating interpretable deep learning models. • Compared to state-of-the-art GAM methods, our NODE-GAM and NODE-GA2M achieve similar
performance on medium-sized datasets while outperforming other GAMs on larger datasets. • We demonstrate that NODE-GAM and NODE-GA2M discover interesting data patterns. • Lastly, we show that NODE-GAM benefits from self-supervised learning that improves performance
when labeled data is limited, and performs better than other GAMs.
We foresee our novel deep learning formulation of the GAMs to be very useful in high-risk domains, such as healthcare, where GAMs have already proved to be useful but stopped short from being applied to new large data collections due to scalability or accuracy issues, as well as settings where access to labeled data is limited. Our novel approach also benefits the deep learning community by adding high accuracy interpretable models to the deep learning repertoire.

heading: 2 BACKGROUND
text: GAM and GA2M: GAMs and GA2Ms are interpretable by design because of their functional forms. Given an input x ∈ RD, a label y, a link function g (e.g. g is log p1−p in binary classification), main effects fj for each feature j, and feature interactions fjj′ , GAM and GA2M are expressed as:
GAM: g(y) = f0 + D∑ j=1 fj(xj), GA2M: g(y) = f0 + D∑ j=1 fj(xj) + D∑ j=1 ∑ j′>j fjj′(xj , xj′).
Unlike full complexity models (e.g. DNNs) that have y = f(x1, ..., xj), GAMs and GA2M are interpretable because the impact of each feature fj and each feature interaction fjj′ can be visualized as a graph (i.e. for fj , x-axis shows xj and y-axis shows fj(xj)). Humans can easily simulate how they work by reading fjs and fjj′ off different features from the graph and adding them together.
GAM baselines: We compare with Explainable Boosting Machine (EBM) (Nori et al., 2019) that implements tree-based GAM and GA2M. We also compare with splines proposed in the 80s (Hastie & Tibshirani, 1990) using Cubic splines in pygam package (Servén & Brummitt, 2018).
Neural Oblivious Decision Trees (NODEs): We describe NODEs for completeness and refer the readers to Popov et al. (2019) for more details. NODE consists of L layers where each layer has m differentiable oblivious decision trees (ODT) of equal depth C. Below we describe a single ODT.
Differentiable Oblivious Decision Trees: An ODT works like a traditional decision tree except for all nodes in the same depth share the same input features and thresholds, which allows parallel computation and makes it suitable for deep learning. Specifically, an ODT of depth C compares C chosen input feature to C thresholds, and returns one of the 2C possible responses. Mathmatically, for feature functions F c which choose what features to split, splitting thresholds bc, and a response vector R ∈ R2C , the tree output h(x) is defined as:
h(x) = R · ([
I(F 1(x) ≤ b1) I(F 1(x) > b1) ] ⊗ [ I(F 2(x) ≤ b2) I(F 2(x) > b2) ] ⊗ · · · ⊗ [ I(FC(x) ≤ bC) I(FC(x) > bC) ]) (1)
Here I is the indicator function, ⊗ is the outer product, and · is the inner product. Both feature functions F c and I prevent differentiability. To make them differentiable, Popov et al. (2019) replace F c(x) as a weighted sum of features:
F c(x) = D∑ j=1 xjentmaxα(F c)j = x · entmaxα(F c). (2)
Here F c ∈ RD are the logits for which features to choose, and entmaxα (Peters et al., 2019) is the entmax transformation which works like a sparse version of softmax such that the sum of the output equals to 1. They also replace the I with entmoid which works like a sparse sigmoid that has output values between 0 and 1. Since all operations are differentiable (entmax, entmoid, outer and inner products), the ODT is differentiable.
Stacking trees into deep layers: Popov et al. (2019) follow the design similar to DenseNet where all tree outputs h(x) from previous layers (each layer consists of total I trees) become the inputs to the next layer. For input features x, the inputs xl to each layer l becomes:
x1 = x, xl = [x,h1(x1), ...,h(l−1)(x(l−1))] for l > 1. (3)
And the final output of the model ŷ(x) is the average of all tree outputs h1, ...,hL of all L layers:
ŷ(x) = 1
LI L∑ l=1 I∑ i=1 hli(x l) (4)

heading: 3 OUR MODEL DESIGN
text: GAM design: See Supp. B for a complete pseudo code. To make NODE a GAM, we make three key changes to avoid any feature interactions in the architecture (Fig. 1). First, instead of letting F c(x) be a weighted sum of features (Eq. 2), we make it only pick 1 feature. We introduce a temperature annealing parameter T that linearly decreases from 1 to 0 for the first S learning steps to make entmaxα(F c/T ) gradually become one-hot:
F c(x) = x · entmaxα(F c/T ), T S steps−−−−→ 0. (5)
Second, within each tree, we make the logits F c the same across depth C i.e. F 1 = · · · = FC = F to avoid any feature interaction within a tree. Third, we avoid the DenseNet connection between two trees that focus on different features j, j′, since they create feature interactions between features j and j′ if two trees connect. Thus we introduce a gate that only allows connections between trees that take the same features. Let Gi = entmaxα(Fi/T ) of the tree i. For tree i in layer l and another tree î in layer l̂ for l̂ < l, the gating weight gl̂ii and the feature function Fli for tree i become:
gl̂ii = Gî ·Gi, Fli(x) = x ·Gi + l−1∑ l̂=1 I∑ î=1 hl̂̂i(x)gl̂ii. (6)
Since G becomes gradually one-hot by Eq. 5, after S steps gîi would only become 1 when Gî = Gi and 0 otherwise. This enforces no feature interaction between tree connections.
Attention-based GAMs (AB-GAMs): To make the above GAM more expressive, we add an attention weight al̂ii in the feature function Fli(x) to decide which previous tree to focus on:
Fli(x) = D∑ j=1 xjGij + l−1∑ l̂=1 I∑ î=1 hl̂̂i(x)gl̂iial̂ii where l−1∑ l̂=1 I∑ î=1 gl̂iial̂ii = 1. (7)
To achieve this, we introduce attention logits Ali for each tree i that after entmax it produces al̂ii:
al̂ii = gl̂iientmaxα(log(gi) +Ali)î. (8) This forces the attention of a tree i that ∑ î al̂ii = 1 for all î that gl̂ii = 1 and al̂ii = 0 when gl̂ii = 0.
The attention logits A requires a large matrix size [I , (l−1)I] for each layer l > 1 which explodes the memory. We instead make A as the inner product of two smaller matrices such that A = BC where B is of size [I , E] and C is of size [E, (l − 1)I], where E is a hyperparameter for the embedding dimension of the attention.
Last Linear layer: Lastly, instead of averaging the outputs of all trees as the output of the model (Eq. 4), we add the last linear layer to be a weighted sum of all outputs:
ŷ(x) = L∑ l=1 I∑ i=1 hli(xl)wli. (9)
Note that in self-supervised learning, wli has multiple output heads to predict multiple tasks.
Regularization: We also include other changes that improves performance. First, we add Dropout (rate p1) on the outputs of trees hli(xl), and Dropout (rate p2) on the final weights wli. Also, to increase diversity of trees, each tree can only model on a random subset of features (η), an idea similar to Random Forest. We also add an `2 penalization (λ) on hli(xl). In binary classification task where labels y are imbalanced between class 0 and 1, we set a constant as log p(y)1−p(y) that is added to the final output of the model such that after sigmoid it becomes the p(y) if the output of the model is 0. We find it’s crucial for `2 penalization to work since `2 induces the model to output 0.
NODE-GA2Ms — extending NODE-GAMs to two-way interactions: To allow two-way interactions, for each tree we introduce two logits F 1 and F 2 instead of just one, and let F c = F (c−1) mod 2+1 for c > 2; this allows at most 2 features to interact within each tree (Fig. 7). Besides temperature annealing (Eq. 5), we make the gating weights gîi = 1 only if the combination of F 1,F 2 is the same between tree î and i (i.e. both trees î and i focus on the same 2 features). We set gîi as:
gîi = min((G 1 i ·G1î )× (G 2 i ·G2î ) + (G 1 i ·G2î )× (G 2 i ·G1î ), 1). (10)
We cap the value at 1 to avoid uneven amplifications as gîi = 2 when G 1 i = G 2 i = G 1 î = G2 î .
Data Preprocessing and Hyperparameters: We follow Popov et al. (2019) to do target encoding for categorical features, and do quantile transform for all features to Gaussian distribution (we find Gaussian works better than Uniform). We use random search to search the architecture space for NODE, NODE-GAM and NODE-GA2M. We use QHAdam (Ma & Yarats, 2018) and average the most recent 5 checkpoints (Izmailov et al., 2018). In addition, we adopt learning rate warmup (Goyal et al., 2017), and do early stopping and learning rate decay on the plateau. More details in Supp. F.
Extracting shape graphs from GAMs: We follow Chang et al. (2021) to implement a function that extracts main effects fj from any GAM model including NODE-GAM, Spline and EBM. The main idea is to take the difference between the model’s outputs of two examples (x1, x2) that have the same values except for feature j. Since the intercept and other main effects are canceled out when taking the difference, the difference f(x2)− f(x1) is equal to fj(x2j )− fj(x1j ). If we query all the unique values of xj , we get all values of fj relative to fj(x1j ). Then we center the graph of fj by setting the average of fj(xj) across the dataset as 0 and add the average to the intercept term f0.
Extracting shape graphs from GA2Ms: Designing a black box function to extract from any GA2M is non-trivial, as each changed feature xj would change not just main effect term fj but also every interactions ∀j′fjj′ that involve feature j. Instead, since we know which features each tree takes, we can aggregate the output of trees into corresponding main fj and interaction terms fjj′ .
Note that GA2M can have many representations that result in the same function. For example, for a prediction value v associated with x2, we can move v to the main effect f2(x2) = v, or the interaction effect f23(x2, ·) = v that involves x2. To solve this ambiguity, we adopt "purification" (Lengerich et al., 2020) that pushes interaction effects into main effects if possible. See Supp. C for details.

heading: 4 RESULTS
text: We first show the accuracy of models in Sec. 4.1 and interpret 2 datasets, Bikeshare and MIMIC2, in Sec. 4.2. In Sec. 4.3, we show that NODE-GAM benefits from self-supervised pre-training and often outperforms other GAMs when labels are limited.
4.1 ARE NODE-GAM AND NODE-GA2M ACCURATE?
We compare our performance on 6 popular binary classification datasets (Churn, Support2, MIMIC2, MIMIC3, Income, and Credit) and 2 regression datasets (Wine and Bikeshare). These datasets are medium-sized with 6k-300k samples and 6-57 features (Table 5). We use 5-fold cross validation to derive the mean and standard deviation for each model. We use 80-20 splits for training and val set. To compare models across datasets, we calculate 2 summary metrics: (1) Rank: we rank the performance on each dataset, and then compute the average rank across all 9 datasets (the lower the rank the better). (2) Normalized Score (NS): for each dataset, we set the worst performance for that dataset as 0 and the best as 1, and scale all other scores linearly between 0 and 1.
In Table 1, we show the performance of all GAMs, GA2Ms and full complexity models. First, we compare 4 GAMs (here NODE-GA2M-main is the purified main effect from NODE-GA2M). We find all 4 GAMs perform similarly and the best GAM in different datasets varies, with Spline as the best in Rank and NODE-GAM in NS. But the differences are often smaller than the standard deviation. Next, both NODE-GA2M and EBM-GA2M perform similarly, with NODE-GA2M better in Rank and EBM-GA2M better in NS. Lastly, within all full complexity methods, XGB performs the best
with not much difference from NODE and RF performs the worst. In summary, all GAMs perform similarly. NODE-GA2M is similar to EBM-GA2M, and slightly outperforms full-complexity models.
In Table 2, we test our methods on 6 large datasets (all have samples > 500K) used in the NODE paper, and we use the same train-test split to be comparable. Since these only provide 1 test split we report standard deviation across multiple random seeds. First, on a cluster with 32 CPU and 120GB memory, Spline goes out of memory on 5 out of 6 datasets and EBM also can not be run on dataset Epsilon with 2k features, showing their lack of ability to scale to large datasets. For 5 datasets that EBM can run, our NODE-GAMs runs slightly better than EBM. But when considering GA2M, NODE-GA2M outperforms EBM-GA2M up to 7.3% in Higgs and average relative improvement of 3.6%. NODE outperforms all GAMs and GA2Ms substantially on Higgs and Year, suggesting both datasets might have important higher-order feature interactions.
4.2 SHAPE GRAPHS OF NODE-GAM AND NODE-GA2M: BIKESHARE AND MIMIC2
In this section, we highlight our key findings, and show the rest of the plots in Supp. H.
Bikeshare dataset: Here we interpret the Bikeshare dataset. It contains the hourly count of rental bikes between the years 2011 and 2012 in Capital bikeshare system located in Washington, D.C. Note that all 4 GAMs trained on Bikeshare are equally accurate with < 0.1% error difference (Table 1).
In Fig. 2, we show the shape plots of 4 features: Hour, Temperature, Month, and Week day. First, Hour (Fig. 2a) is the strongest feature with two peaks around 9 AM and 5 PM, representing the time that people commute, and all 4 models agree. Then we show Temperature in Fig. 2b. Here temperature is normalized between 0 and 1, where 0 means -8°C and 1 means 39°C. When the weather is hot (Temp > 0.8, around 30°C), all models agree rental counts decrease which makes sense. Interestingly, when it’s getting colder (Temp < 0.4, around 11°C) there is a steady rise shown by NODE-GAM, Spline and EBM but not NODE-GA2M (blue). Since it’s quite unlikely people rent more bikes when it’s getting colder especially below 0°C, the pattern shown by GA2M seems more plausible. Similarly, in feature Month (Fig. 2c), NODE-GA2M shows a rise in summer (month 6− 8) while others indicate a strong decline of rental counts. Since we might expect more people to rent bikes during summer since it’s warmer, NODE-GA2M might be more plausible, although we might explain it due to summer vacation fewer students ride bikes to school. Lastly, for Weekday (Fig. 2d) all 4 models agree with each other that the lowest number of rentals happen at the start of the week (Sunday and Monday) and slowly increase with Saturday as the highest number.
In Fig. 3, we show the 4 feature interactions (out of 67) from our NODE-GA2M. The strongest effect happens in Hr x Working day (Fig. 3(a)): this makes sense since in working day (orange), people usually rent bikes around 9AM and 5PM to commute. Otherwise, if the working day is 0 (blue), the number peaks from 10AM to 3PM which shows people going out more often in daytime. In Hr x Weekday (Fig. 3(b)), we can see more granularly that this commute effect happens strongly on Monday to Thursday, but on Friday people commute a bit later, around 10 AM, and return earlier, around 3 or 4 PM. In Hr x Temperature (Fig. 3(c)), it shows that in the morning rental count is high when it’s cold, while in the afternoon the rental count is high when it’s hot. We also find in Hr x Humidity (Fig. 3(d)) that when humidity is high from 3-6 PM, people ride bikes less. Overall these interpretable graphs enable us to know how the model predicts and find interesting patterns.
MIMIC2: MIMIC2 is the hospital ICU mortality prediction task (Johnson et al., 2016a). We extract 17 features within the first 24 hour measurements, and we use mean imputation for missingness.
We show the shape plots (Fig. 4) of 4 features: Age, PFratio, Bilirubin and GCS. In feature Age (Fig.. 4(a)), we see the all 4 models agree the risk increases from age 20 to 90. Overall NODEGAM/GA2M are pretty similar to EBM in that they all have small jumps in a similar place at age 55, 80 and 85; spline (red) is as expected very smooth. Interestingly, we see NODE-GAM/GA2M shows risk increases a bit when age < 20. We think the risk is higher in younger people because this is generally a healthier age in the population, so their presence in ICU indicates higher risk conditions.
In Fig. 4(b), we show PFratio: a measure of how well patients oxygenate the blood. Interestingly, NODE-GAM/GA2M and EBM capture a sharp drop at 332. It turns out that PFratio is usually not measured for healthier patients, and the missing values have been imputed by the population mean 332, thus placing a group of low-risk patients right at the mean value of the feature. However, Spline (red) is unable to capture this and instead have a dip around 300-600. Another drop captured by NODE-GAM/GA2M from 400− 500 matches clinical guidelines that > 400 is healthy. Bilirubin shape plot is shown in Fig. 4(c). Bilirubin is a yellowish pigment made during the normal breakdown of red blood cells. High bilirubin indicates liver or bile duct problems. Indeed, we can see risk quickly goes up as Bilirubin is > 2, and all 4 models roughly agree with each other except for Spline which has much lower risk when Billirubin is 80, which is likely caused by Spline’s smooth inductive bias and unlikely to be true. Lastly, in Fig. 4(d) we show Glasgow Coma Scale (GCS): a bedside measurement for how conscious the patient is with 1 in a coma and 5 as conscious. Indeed, we find the risk is higher for patients with GCS=1 than 5, and all 4 models agree.
In Fig. 5, we show the 4 of 154 feature interactions learned in the NODE-GA2M. First, in Age x Bilirubin (Fig. 5(a)), when Billirubin is high (>2), we see an increase of risk (blue) in people with age 18− 70. Risk decreases (red) when age > 80. Combined with the shape plots of Age (Fig. 4(a)) and Bilirubin (Fig. 4(c)), we find this interaction works as a correction effect: if patients have Bilirubin > 2 (high risk) but are young (low risk), they should have a higher risk than what their main (univariate) effects suggest. On the other hand, if patients have age > 80 (high risk) and Bilirubin > 2 (high risk), they already get very high risk from main effects, and in fact the interaction effect is negative to correct for the already high main effects. It suggests that Billirubin=2 is an important threshold that should affect risk adjustments.
Also in GCS x Bilirubin plot (Fig. 5(b)), we find similar effects: if Bilirubin >2, the risk of GCS is correctly lower for GCS=1,2 and higher for 3-5. In Fig. 5(c) we find patients with GCS=1-3 (high risk) and age>80 (high risk), surprisingly, have even higher risk (blue) for these patients; it shows models think these patients are more in danger than their main effects suggest. Finally, in Fig. 5(d) we show interaction effect GCS x PFratio. We find PFratio also has a similar threshold effect: if PFratio > 400 (low risk), and GCS=1,2 (high risk), model assigns higher risk for these patients while decreasing risks for patients with GCS=3,4,5.

heading: 4.3 SELF-SUPERVISED PRE-TRAINING
text: By training GAMs with neural networks, it enables self-supervised learning that learns representations from unlabeled data which improves accuracy in limited labeled data scenarios. We first learn a NODE-GAM that reconstructs input features under randomly masked inputs (we use 15% masks). Then we remove and re-initialize the last linear weight and fine-tune it on the original targets under limited labeled data. For fine-tuning, we freeze the embedding and only train the last linear weight for the first 500 steps; this helps stabilize the training. We also search smaller learning rates [5e−5, 1e−4, 3e−4, 5e−4] and choose the best model by validation set. We compare our self-supervised model (SS) with 3 other baselines: (1) NODE-GAM without self-supervision (No-SS), (2) EBM and (3) Spline. We randomly search 15 attention based AB-GAM architectures for both SS and No-SS.
In Figure 6, we show the relative improvement over the AUC of No-SS under variously labeled data ratio 0.5%, 1%, 2%, 5%, and 10% (except Credit which 0.5% has too few positive samples and thus crashes). And we run 3 different train-test split folds to derive mean and standard deviation. Here the relative improvement means improvement over No-SS baselines. First, we find NODE-GAM with self-supervision (SS, blue) outperforms No-SS in 6 of 7 datasets (except Income) with MIMIC2 having the most improvement (6%). This shows our NODE-GAM benefits from self-supervised pre-training. SS also outperforms EBM in 3 out of 6 datasets (Churn, MIMIC2 and Credit) up to 10% improvement in Churn, demonstrating the superiority of SS when labeled data is limited.

heading: 5 LIMITATIONS, DISCUSSIONS AND CONCLUSIONS
text: Although we interpret and explain the shape graphs in this paper, we want to emphasize that the shown patterns should be treated as an association not causation. Any claim based on the graphs requires a proper causal study to validate it.
In this paper, we assumed that the class of GAMs is interpretable and proposed a new deep learning model in this class, so our NODE-GAM is as interpretable to other GAMs. But readers might wonder if GAMs are interpretable to human. Hegselmann et al. (2020) show GAM is interpretable for doctors in a clinical user study; Tan et al. (2019) find GAM is more interpretable than the decision tree helping users discover more patterns and understand feature importance better; Kaur et al. (2020) compare GAM to post-hoc explanations (SHAP) and find that GAM significantly makes users answer questions more accurately, have higher confidence in explanations, and reduce cognitive load.
In this paper we propose a deep-learning version of GAM and GA2M that automatically learn which main and pairwise interactions to focus without any two-stage training and model ensemble. Our GAM is also more accurate than traditional GAMs in both large datasets and limited-labeled settings. We hope this work can further inspire other interpretable design in the deep learning models.

heading: ETHICS STATEMENT
text: A potential misuse of this tool is to claim causal relationships of what GAM models learned. For example, one might wrongfully claim that having asthma is good for pneumonia patients (Caruana et al., 2015). This is obviously a wrong conclusion predicated by the bias in the dataset.
The interpretable models are a great tool to discover potential biases hiding in the data. Especially in high-dimensional datasets when it’s hard to pinpoint where the bias is, GAM provides easy visualization that allows users to confirm known biases and examine hidden biases. We believe that such models foster ethical adoption of machine learning in high-stakes settings. By providing a deep learning version of GAMs, our work enables the use of GAM models on larger datasets thus increasing GAM adoption in real-world settings.

heading: REPRODUCIBILITY STATEMENT
text: We released our code anonymously in the supplementary with instructions and hyperparameters to reproduce our final results. Our datasets can be automatically downloaded to the directory with train and test fold using included code. We report the details of our datasets in Supp. E, and hyperparameters for both our model and baselines in Supp. F. We list the hyperparameters corresponding to the best performance in Supp. G.

heading: A COMPARISON TO NAM (AGARWAL ET AL., 2020) AND THE UNIVARIATE
text: NETWORK IN NID (TSANG ET AL., 2017)
First, we compare with NAM’s (Agarwal et al., 2020) performance without and with the ensemble. Since NAM requires an extensive hyperparameter search and to be fair with NAM, we focus on 3 of their 5 datasets (MIMIC2, COMPAS, Credit) and use their reported best hyperparameters to compare. We find that NODE-GAM is better than NAM without ensemble consistently across 3 datasets.
As the reviewer points out, in the NID paper (Tsang et al., 2017) they also have a similar idea to NAM that trains a univariate network alongside an MLP to model the main effect. But there are two key differences between the univariate network in NID and NAM: (1) NAM proposes a new jumpy activation function - ExU - that can model quick, non-linear changes of the inputs as part of their hyperparameter selection, and (2) NAM uses multiple networks to ensemble. To be thorough, we compare with NAM that only uses normal units to resemble the univariate network in NID. We also considered whether removing the ensemble would disproportionately impact ExU activations since ExU activations are more prone to overfitting.
We show the performance in Table 4. We show the results in 2 of their 5 datasets (MIMIC-II, Credit) that find ExU perform better than normal units (in other 3 datasets normal units perform better). First, we find that after ensemble normal units perform quite similar to ExU in MIMIC2 but worse in Credit. And given that in 3 other datasets NAM already finds normal units to perform better, we think normal units and ExU probably have similar accuracy. Besides, ensemble helps improve performance much more for ExU units but not so much for normal units, since ExU is a more low-bias high-variance unit that benefits more from ensembles.

heading: B PSEUDO-CODE FOR NODE-GAM
text: Here we provide the pseudo codes for our model in Alg. 1-4. We highlight our key changes that make NODE as a GAM in red, and the new architectures or regularization in blue. We show a single GAM decision tree in Alg. 1, a single GA2M tree in Alg. 2, model algorithm in Alg. 3, and the model update in Alg. 4.
Algorithm 1 A GAM differentiable oblivious decision tree (ODT) Input: Input X ∈ RD, Temperature T (T −→ 0), Previous layers’ outputs Xp ∈ RP , Previous layers’ feature selection Gp ∈ RP×D, Attention Matrix A ∈ RP Hyperparameters: Tree Depth C, Column subsample ratio η Trainable Parameters: Feature Selection Logits F ∈ RD , Split Thresholds b ∈ RC , Split slope S ∈ RC , Node weights W ∈ R2C
if η < 1 then n = max(Dη, 1) . Number of subset features i = shuffle(range(D))[int(n):] . Randomly choose features to exclude F [i] = − inf . Exclude features end if G = EntMax(F /T ) . Go through EntMax to generate soft one-hot vector (Eq. 5) K = X ·G . Pick 1 feature softly if XP is not None then
g = GPG ∈ RP . Calculate gating weights g (Eq. 6) g′ = g/ ∑ g if A is None else g·entmaxα(log(g) +A) . Attention-based GAM (Eq.8)
K = K +Xp · g′ . Add previous outputs with g normalized to 1 end if H = EntMoid((K − b)/S) ∈ RC . Generate soft binary value
e =
([ H1
(1−H1)
] ⊗ · · · ⊗ [ (HC)
(1−HC)
]) ∈ R2C . Go through the decision tree
h = e ·W . Select one weight value softly as the output
Return: h, G . Return tree response h and feature selection G
Algorithm 2 A GA2M differentiable oblivious decision tree (ODT) Input: Input X ∈ RD, Temperature T (T −→ 0), Previous layers’ outputs Xp ∈ RP , Previous layers’ feature selection G1p,G 2 p ∈ RP×D, Attention Matrix A ∈ RP
Hyperparameters: Tree Depth C, Column subsample ratio η Trainable Parameters: Feature Selection Logits F 1,F 2 ∈ RD , Split Thresholds b ∈ RC , Split slope S ∈ RC , Node weights W ∈ R2C
if η < 1 and first time running then n = max(Dη, 1) . Number of subset features i = shuffle(range(D))[n:] . Randomly exclude features F 1[i] = − inf , F 2[i] = − inf . Exclude features end if G1 = EntMax(F 1/T ), G2 = EntMax(F 2/T ) . Get soft one-hot vector (Eq. 5) K1 = X ·G1, K2 = X ·G2 . Pick 1 feature softly if XP is not None then
g = min((G1 ·G1P )× (G2 ·G2P ) + (G1 ·G2P )× (G2 ·G1P ), 1) . Gating weights (Eq. 10) g′ = g/ ∑ g if A is None else g·entmaxα(log(g) +A) . Attention-based GAM (Eq.8) K1 = K1 +Xp · g′ . Add previous outputs with g normalized to 1 K2 = K2 +Xp · g′ K = [K1,K2,K1, ...,K2] ∈ RC . Alternating between K1 and K2 end if H = EntMoid((K − b)/S) ∈ RC . Generate soft binary value
e =
([ H1
(1−H1)
] ⊗ · · · ⊗ [ (HC)
(1−HC)
]) ∈ R2C . Go through the decision tree
h = e ·W . Select one weight value softly as the output
Return: h, [G1,G2] . Return tree response h and feature selection [G1,G2]
Algorithm 3 The NODE-GAM / NODE-GA2M algorithm 1: Input: Input X ∈ RD 2: Hyperparameters: number of layers L, number of trees I per layer, tree depth C, current
optimization step s, temperature annealing step S, Attention Embedding E 3: Trainable Parameters: the decision trees Ml in each layer l (either GAM trees (Alg. 1) or
GA2M trees (Alg. 2)), the final output weights WL ∈ R(LI) and bias w0 4: if E > 0 then . Use attention-based GAM 5: Initialize Bl ∈ R(l−1)I×E and Cl ∈ RE×I for l = 2...L 6: end if 7: T = 10−2(s/S) if s ≤ S else 0 . Slowly decrease temperature to 0 8: Xp = None, Gp = None . Initialize previous trees’ outputs Xp and feature selections Gp 9: for l = 1 to L do 10: Al = None if E = 0 or l = 1 else BlCl . Calculate attention matrix 11: hl,Gl = Ml(X, T ,Xp,Gp,Al) . Run total I trees in Ml by Alg. 1 or 2 12: hl = Dropout(hl) . Dropout rate p1 13: Xp = hl if Xp is None else [Xp,hl] . Concatenate outputs hl 14: Gp = Gl if Gp is None else [Gp,Gl] . Concatenate feature selection Gl 15: end for 16: WL = Dropout(WL) . Dropout rate p2 17: R = Xp ·WL + w0 . Go through last linear layer 18: Return: R, XP . Return model response R and all trees’ outputs XP
Algorithm 4 Our model’s update 1: Input: An input X ∈ RD, target y, Node-GAM modelM 2: R,XP =M(X) . Run Node-GAM (Alg. 3) 3: L = BCELoss(y,R) if binary classification else MSELoss(y,R) 4: L = L+ λ ∑ (XP )
2 . Add the `2 on the output of trees 5: Optimize L by Adam optimizer
C PURIFICATION OF GA2M
Note that GA2M can have many representations that result in the same function. For example, for a prediction value v associated with x2, we can move v to the main effect f2(x2) = v, or the interaction effect f23(x2, ·) = v that involves x2. To solve this ambiguity, we adopt "purification" (Lengerich et al., 2020) that pushes interaction effects into main effects if possible.
To purify an interaction fjj′ , we first bin continuous feature xj into at most K quantile bins with K unique values x1j , ...x K j and for xj′ as well. Then for every x k j , we move the average a k j of interactions fjj′ to main effects fj :
∀Kk=1xkj , akj = 1
NK′ K′∑ k′=1 fjj′(x k j , x k′ j′ ), fjj′(x k j , x k′ j′ ) = fjj′(x k j , x k′ j′ )−akj , fj(xkj ) = fj(xkj )+akj
This is one step to purify fjj′ to fj . Then we purify fjj′ to fj′ , and so on until all akj and a k′
j′ are close to 0.
D NODE-GA2M FIGURES
Here we show the architecutre of NODE-GA2M in Figure 7.

heading: E DATASET DESCRIPTIONS
text: Here we describe all 8 datasets we use and we summarize them in Table 5.
• Churn: this is to predict which user is a potential subscription churner for telecom company. https://www.kaggle.com/blastchar/telco-customer-churn
• Support2: this is to predict mortality in the hospital by several lab values. http:// biostat.mc.vanderbilt.edu/DataSets
• MIMIC-II and MIMIC-III dataset (Johnson et al., 2016b): this is an ICU patient datasets to predict mortality of patients in a tertiary academic medical center in Boston, MA, USA.
• Income: UCI Dua & Graff (2017). This is a dataset from census collected in 1994, and the goal is to predict who has income >50K/year. https://archive.ics.uci.edu/ ml/datasets/adult
• Credit: this is to predict which transaction is a fraud. The features provided are the coefficient of PCA components to protect privacy. https://www.kaggle.com/mlg-ulb/ creditcardfraud
• Bikeshare (Dua & Graff, 2017): this is the hourly bikeshare rental counts in Washington D.C., USA. https://archive.ics.uci.edu/ml/datasets/bike+ sharing+dataset
• Wine (Dua & Graff, 2017): this is to predict the wine quality based on a variety of lab values. https://archive.ics.uci.edu/ml/datasets/wine+quality
For 6 datasets used in NODE, we use the scripts from NODE paper (https://github.com/ Qwicen/node) which directly downloads the dataset. Here we still cite and list their sources:
• Click: https://www.kaggle.com/c/kddcup2012-track2
• Higgs: UCI (Dua & Graff, 2017) https://archive.ics.uci.edu/ml/ datasets/HIGGS
• Epsilon: https://www.k4all.org/project/large-scale-learning-challenge/
• Microsoft: https://www.microsoft.com/en-us/research/project/ mslr/

heading: E.1 PREPROCESSING
text: For NODE and NODE-GAM/GA2M, we follow Popov et al. (2019) to do target encoding for categorical features, and do quantile transform1 with 2000 bins for all features to Gaussian distribution (we find Gaussian performs better than Uniform). We find adding small gaussian noise (e.g. 1e-5) when fitting quantile transformation (but no noise in transformation stage) is crucial to have mean 0 and standard deviation close to 1 after transformation.

heading: F HYPERPARAMETERS SELECTION
text: In order to tune the hyperparameters, we performed a random stratified split of full training data into train set (80%) and validation set (20%) for all datasets. For datasets we compile of medium-sized (Income, Churn, Credit, Mimic2, Mimic3, Support2, Bikeshare), we do a 5-fold cross validation for 5 different test splits. For datasets in NODE paper (Click, Epsilon, Higgs, Microsoft, Yahoo, Year), we use train/val/test split provided by the NODE paper author. Since they only provide 1 test split, we report standard deviation by different random seeds on these datasets. For medium-sized datasets, we only tune hyperparameters on the first train-val-test fold split, and fix the hyperparameters to run the rest of 4 folds. This means that we do not search hyperparameters per fold to avoid computational overheads. All NODE, NODE-GAM/GA2M are run with 1 TITAN Xp GPU, 4 CPU and 8GB memory. For EBM and Spline, they are run with a machine with 32 CPUs and 120GB memory.
Below we describe the hyperparameters we use for each method:

heading: F.1 EBM
text: For EBM, we set inner_bags=100 and outer_bags=100 and set the maximum rounds as 20k to make sure it converges; we find EBM performs very stable out of this choice probably because we set total bagging as 10k that makes it stable; other parameters have little effect on final performance.
For EBM GA2M, we search the number of interactions for 16, 32, 64, 128 and choose the best one on validation set. On large datasets we set the number of iterations as 64 as we find it performs quite well on medium-sized datasets.
1sklearn.preprocessing.quantile_transform

heading: F.2 SPLINE
text: We use the cubic spline in PyGAM package (Servén & Brummitt, 2018) that we follow Chang et al. (2021) to set the number of knots per feature to a large number 50 (we find setting it larger would crash the model), and search the best lambda penalty between 1e-3 to 1e3 for 15 times and return the best model.
F.3 NODE, NODE-GA2M AND NODE
We follow NODE to use QHAdam (Ma & Yarats, 2018) and average the most recent 5 checkpoints. In addition, we adopt learning rate warmup at first 500 steps. And we early stop our training for no improvement for 11k steps and decay learning rate to 1/5 if no improvement happens in 5k steps.
Here we list the hyperparameters we find works quite well and we do not do random search on these hyperparameters:
• optimizer: QHAdam (Ma & Yarats, 2018) (same as NODE paper)
• lr_warmup_steps: 500
• num_checkpoints_avged: 5
• temperature_annealing_steps (S): 4k
• min_temperature: 0.01 (0.01 is small enough for making one-hot vector. And after S steps we set the function to produce one-hot vector exactly.)
• batch_size: 2048, or the max batch size that fits in GPU memory with minimum 128 Maximum training time: 20 hours. This is just to avoid model training for too long.
We use random search to find the best hyperparameters which we list the range in below. We list the random search range for NODE:
• num_layers: {2, 3, 4, 5}
• total tree counts (= num_trees × num_layers): {500, 1000, 2000, 4000} • depth: {2, 4, 6}
• tree_dim: {1, 2, 3}
• output_dropout (p1): {0, 0.1, 0.2}
• colsample_bytree: {1, 0.5, 0.1, 1e-5}
• lr: {0.01, 0.005}
• l2_lambda: {0., 1e-7, 1e-6, 1e-5}
• add_last_linear (to add last linear weight or not): {0, 1}
• last_dropout (p2, only if add_last_linear=1): {0, 0.1, 0.2, 0.3}
• seed: uniform distribution [1, 100]
For NODE-GAM and NODE-GA2M, we have additional parameters:
• arch: {GAM, AB-GAM}
• dim_att (dimension of attention embedding E): {8, 16, 32}
We show the best hyperparameters for each dataset in Section G.
F.4 XGBOOST
For large datasets in NODE, we directly report the performance from the original NODE paper. For medium-sized data, we set the depth of xgboost as 3, and learning rate as 0.1 with n_estimators=50k and set early stopping for 50 rounds to make sure it converges.

heading: F.5 RANDOM FOREST (RF)
text: We use the default hyperparameters from sklearn and set the number of trees to a large number 1000.

heading: G BEST HYPERPARAMETERS FOUND IN EACH DATASET
text: Here we report the best hyperparameters we find for 9 medium-sized datasets in Table 6 (NODEGAM), Table 7 (NODE-GA2M), and Table 8 (NODE). We report the best hyperparameters for large datasets in Table 9 (NODE-GAM) and Table 10 (NODE-GA2M).

heading: H COMPLETE SHAPE GRAPHS IN BIKESHARE AND MIMIC2
text: We list all main effects of Bikeshare in Fig. 8 and top 16 interactions effects in Fig. 9. We list all main effects of MIMIC2 in Fig. 10 and top 16 interactions effects in Fig. 11.

Review summaries:
Generalized Additive Models (GAMs) are a class of interpretable models with a long history of use in these high-risk domains, but they lack desirable features of deep learning such as differentiability and scalability.   The authors propose a neural GAM (NODE-GAM) and neural GA2M (NODE-GA2M) that scale well and perform better than other GAMs on large datasets, while remaining interpretable.  Popov et al. (2019) developed NODE that mimics an ensemble of decision trees, but suffers from lack of interpretability similarly to other ensemble and deep learning models. Agarwal et al. (2020) developed a Neural Additive Model (NAM) whose deep learning architecture is a GAM, but does not scale well.  The authors combine the above mentioned 2 approaches via a 3 step process: 1. Instead of letting the ODT feature function be a weighted sum of features, they make it pick only one feature. 2. Within each tree, they make a logit layer the same across all depths to remove interactions between features. 3. Finally, they avoid the DenseNet connection between two trees that focus on different features.  The authors then output performance on 6 popular binary classification datasets (Churn, Support2, MIMIC2, MIMIC3, Income, and Credit) and 2 regression datasets (Wine and Bikeshare). Looking at the experiments, their NODE GA^2M version seems to perform similar to the compared state of the art.  The authors additionally perform an ablation experiment by testing their method's semi supervised learning capabilities, show that it is able to learn on mask data and be fine tuned on a limited labeled dataset.  Finally, the conclusion has arguments to support the interpretability of GAMs and to defend their model's interpretability

 - The authors develop the NODE-GAM and NODE-GA$^2$M architectures by modifying the previously-developed NODE architecture with constraints and a gating mechanism to ensure that, both within each tree and across layers, the model is only allowed to learn feature interactions of order 1 (for NODE-GAM) or 2 for (NODE-GA$^2$M). - They show that their proposed NODE-GAM and NODE-GA$^2$M architectures perform comparably to the baselines of Explainable Boosting Machines, tree-based GAM/GA$^2$M and the traditional spline-based GAMs on medium-sized datasets while outperforming them on larger datasets. They also show that NODE-GAM can benefit from self-supervised learning by pretraining a model to reconstruct the original input from a masked input. - They apply NODE-GAM and NODE-GA$^2$M to real-world datasets and showcase the patterns uncovered.

The authors proposed novel architectures for neural GAM and GA2M, which preserves the interpretability of GAM and leveraging the deep learning architectures for performance gains. The new method was assessed on 14 different dataset which covers a wide variety of prediction tasks. The performance of methods (NODE-GAM and NODE-GA2M) were comparable with the other GAM methods. Compared to other similar methods, the proposed method has a better scale ability, and better performance on larger datasets.