GAMformer: In-Context Learning for Generalized Additive Models

## Abstract

Abstract Generalized Additive Models (GAMs) are widely recognized for their ability to create fully interpretable machine learning models for tabular data. Traditionally, training GAMs involves iterative learning algorithms, such as splines, boosted trees, or neural networks, which refine the additive components through repeated error reduction. In this paper, we introduce GAMformer , the first method to leverage in-context learning to estimate shape functions of a GAM in a single forward pass, representing a significant departure from the conventional iterative approaches to GAM fitting. Building on previous research applying in-context learning to tabular data, we exclusively use complex, synthetic data to train GAMformer, yet find it extrapolates well to real-world data. Our experiments show that GAMformer performs on par with other leading GAMs across various classification benchmarks while generating highly interpretable shape functions.

## 1 Introduction

The growing importance of interpretability in machine learning is evident, especially in areas where transparency, fairness, and accountability are critical (Barocas and Selbst, ; Rudin et al., ). Interpretable models are essential for building trust between humans and AI systems by allowing users to understand the reasoning behind the model’s predictions and decisions (Ribeiro et al., ). This is crucial in safety-critical fields like healthcare, where incorrect or biased decisions can have severe consequences (Caruana et al., ). Additionally, interpretability is vital for regulatory compliance in sectors like finance and hiring, where explaining and justifying model outcomes is necessary (Arun et al., ; Dattner et al., ). Interpretable models also help detect and mitigate bias by revealing the factors influencing predictions, ensuring fair and unbiased decisions across different population groups (Mehrabi et al., ).

Generalized Additive Models (GAMs) have proven a popular choice for interpretable modeling due to their high accuracy and interpretability. In GAMs, the target variable is expressed as a sum of non-linearly transformed features. This approach strikes a balance between the interpretability of linear models and the flexibility of capturing non-linear relationships between features and the target variable (Hastie and Tibshirani, ). A wide variety of GAMs exist, differing in the non-linear functions used to transform features and the methods employed to fit these functions to training data. Traditionally, GAMs have used splines in conjunction with the backfitting algorithm (Hastie and Tibshirani, ), while Explainable Boosting Machines (EBMs) utilize decision trees and cyclic gradient boosting (Lou et al., ; ; Caruana et al., ). More recently, Neural Additive Models (NAMs) have employed multilayer perceptrons (MLPs) optimized via gradient descent (Agarwal et al., ). All existing GAM variants share the need for an iterative optimization algorithm to fit the shape functions, which introduces additional hyperparameters for optimization and regularization that require tuning (Siems et al., ; Kovács, ).

Recently, in-context learning (ICL) has emerged as a powerful paradigm for eliminating explicit optimization in models. This breakthrough was first observed in large language models (Brown et al., ), where a model trained in an unsupervised manner on vast amounts of unlabeled data can learn to execute a new task when presented with examples, without any further optimization or updates to its parameters. Since then, ICL has been applied to various domains, including multi-modal foundation models (Li et al., ) and time-series forecasting (Dooley et al., ). Of particular relevance to our work is TabPFN (Hollmann et al., ; Müller et al., ), a transformer model pretrained on complex, synthetic tabular data. This pretraining enables TabPFN to generalize to real-world data when presented with a dataset in the form of in-context examples, demonstrating the potential of ICL.

We introduce GAMformer (see  Figure 1), the first GAM method to estimate shape functions using ICL in a single forward pass. GAMformer distinguishes itself from existing GAM methods by employing a non-parametric, binned representation of shape functions, thus eliminating the need to impose a specific model class. Similar to TabPFN, our model is trained exclusively on large-scale synthetic datasets, yet demonstrates robust performance on real-world data. During training, GAMformer estimates shape functions for each feature based on the training data’s features and labels. These estimated functions are then utilized to generate predictions for test data points by summing the shape function values across features. The model is trained end-to-end based on the GAM’s predictions, ensuring that it learns to accurately construct shape functions for reliable predictions.

Figure: Figure 1: GAMformer’s forward pass on a new dataset with three features ($x_{1}$, $x_{2}$, $x_{3}$) and label $y$ and two data points: (1) For each data point, we bin all features, one-hot encode them, embed the resulting vectors and add the label of the data point. (2) We alternate between applying attention across the features and the data points, allowing us to handle varying numbers of each. (3) We decode per-feature shape functions using a shared MLP decoder. (4) We infer the prediction for test data points by looking up and adding each feature’s shape function value (red bins) forming the GAM prediction. (5) Finally, we compute the loss based on the prediction allowing the end-to-end training of the shape function estimation based on (in our case, *synthetic*) training datasets.
Refer to caption: x1.png

Our main contributions can be summarized as follows:

- •
We introduce GAMformer, the first method to utilize in-context learning with sequence-to-sequence models to form shape functions in a single forward pass, eliminating the need for iterative learning and hyperparameter tuning.
- •
Our experimental results demonstrate GAMformer’s capacity to match the accuracy of leading GAMs on various classification benchmarks.
- •
Our case study on MIMIC-II demonstrates how GAMformer can be applied to real-world data to generate interpretable models and insights of that data.

## 2 Background and Related Work

In this section, we provide some background and related work on generalized additive models and in-context learning.

### 2.1 Generalized Additive Models.

Generalized Additive Models (GAMs) (Hastie and Tibshirani, ) emerged as a generalization of Generalized Linear Models (Nelder and Wedderburn, ) which include non-linear transformations of the input features. The structure of a GAM is given by:

$$ $g(\mathbb{E}[y|x])=\beta+\sum\nolimits_{i=1}^{p}f_{i}(x_{i}),$ (1) $$

where $x=(x_{1},\dots x_{p})\in\mathcal{X}\subseteq\mathbb{R}^{p}$ is the input with $p$ features, $y\in\mathcal{Y}\subseteq\mathbb{R}^{m}$ is the response variable, and $f_{i}:\mathbb{R}\to\mathbb{R}$ are univariate functions termed shape functions that capture the individual contributions of each feature. The intercept $\beta\in\mathbb{R}$ is a learnable bias term, and $g:\mathbb{R}\to\mathbb{R}$ is the link function that connects the expected outcome to the linear predictor, examples of which include the logit or softmax function for binary or multiclass classification or the identity function for linear regression. The shape functions $f_{i}$ in GAMs, also sometimes called partial dependence plots, allow for an interpretable representation of each feature’s effect, akin to the role of coefficients in linear regression, thus enabling practitioners to inspect the learned potentially non-linear relationships.

Traditional GAMs often use splines and backfitting (Hastie and Tibshirani, ), enhanced by penalized regression splines (Wood, ) and fast fitting algorithms (Wood, ). Spline-based GAMs use the backfitting algorithm, iteratively updating each shape function to fit the residuals of others until convergence. More recent advances include Explainable Boosting Machines (EBMs) (Lou et al., ; ; Caruana et al., ), which use decision trees to model shape functions via cyclic gradient boosting. This approach learns each feature’s contribution iteratively in a round-robin manner, mitigating collinearity effects and accurately modeling steps in the data, which is crucial for capturing discontinuities like treatment effects in medical data.
On the other hand, Neural Additive Models (NAMs) (Agarwal et al., ) and follow up works (Chang et al., ; Dubey et al., ; Radenovic et al., ; Xu et al., ; Enouen and Liu, ; Bouchiat et al., ) use multilayer perceptrons (MLPs) as non-linear transformations to model the shape functions $f_{i}$. As a result, NAMs can be optimized using variants of gradient descent by leveraging automatic differentiation frameworks.
Finally, GAMs have also found applications in time-series forecasting, with models such as Prophet (Taylor and Letham, ) and NeuralProphet (Triebe et al., ). For a more comprehensive related work refer to Appendix A.

### 2.2 In-Context Learning & Prior-Data Fitted Networks

In-Context Learning (ICL) was first demonstrated alongside the introduction of GPT-3 (Brown et al., ), where the authors showed that Transformer models (Vaswani et al., ) could learn to perform tasks solely from input examples, without explicit training or fine-tuning, after self-supervised pre-training. This capability marks a significant paradigm shift from the traditional machine learning paradigm of in-weights learning, where the parameters of a model are adjusted in order to learn a new task. The discovery of ICL has led to numerous investigations into the mechanisms used by trained transformers that enable ICL. Olsson et al. () found that a two-layer attention-only network can develop “induction heads”, a mechanism that outputs the token succeeding a previous instance of the current token, precisely when its ICL performance increases. Chan et al. () investigated the properties of the data that contribute to the emergence of ICL abilities, while Reddy () identified factors responsible for the abrupt emergence of induction heads.

Of particular relevance to this paper are Prior-Data-Fitted Networks (PFNs) (Müller et al., ; Hollmann et al., ), which showed that a transformer trained on complex synthetic data generated using random causal graphs can be used for tabular classification. From a Bayesian perspective, such causal graphs $\phi$ sampled from a hypothesis space $\Phi$ (the prior), define a mechanism that describes the relationship between the input and output variables. In TabPFNs (Hollmann et al., ), a synthetic dataset $D\sim p(D)=\mathbb{E}_{\phi\sim p(\phi)}[p(D|\phi)]$ is repeatedly constructed by propagating samples $x\sim p(\mathcal{X})$ from the input space through a randomly sampled structural causal model (SCM), $\phi\sim p(\phi)$, to obtain the corresponding $y$ values. We denote the dataset containing $N$ such examples as the set $D:=\{(x^{(n)},y^{(n)})\}_{n=1}^{N}$. To simulate practical inference scenarios, the dataset $D$ is split into $D_{\text{train}}$ and the context dataset $D_{\text{test}}=D\setminus D_{\text{train}}$. The transformer model parses the pairs $(x_{\text{train}},y_{\text{train}})\in D_{\text{train}}$, as well as $x_{\text{test}}$, as single input tokens and its parameters $\theta$ are updated to minimize the negative log likelihood on the test held-out examples:

$$ $\mathbb{E}_{(D_{\text{train}}\cup(x_{\text{test}},y_{\text{test}}))\sim p(D)}[ -\log q_{\theta}(y_{\text{test}}|x_{\text{test}},D_{\text{train}})].$ (2) $$

Müller et al. () showed that by minimizing this loss, TabPFN approximates the true posterior predictive distribution

$$ $p(y_{\text{test}}|x_{\text{test}},D_{\text{train}})=\int_{\Phi}p(y_{\text{test }}|x_{\text{test}},\phi)p(\phi|D_{\text{train}})d\phi\propto\int_{\Phi}p(y_{ \text{test}}|x_{\text{test}},\phi)p(D_{\text{train}}|\phi)p(\phi)d\phi$ (3) $$

on a new input point from the test set $x_{\text{test}}$ up to an additive constant. This paradigm has since been extended to time-series forecasting (Dooley et al., ), hyperparameter optimization (Müller et al., ; Adriaensen et al., ; Rakotoarison et al., ) and the prediction of neural network weights (Müller et al., ). Similarly, Conditional Neural Processes (Garnelo et al., ) also perform a form of ICL, using a neural architecture with weights meta-learned on real data. (Nguyen and Grover, ) extended Neural Processes to a transformer architecture, leading to an architecture similar to PFNs. *GAMformer* builds on top of TabPFN by training a transformer on synthetically generated datasets to estimate the shape function per feature and computing predictions by adding the individual shape function values.

## 3 GAMformer

We first provide a high-level overview of how GAMformer works before delving into the details of each of its components. GAMformer follows a two-step approach that first fits a GAM on training data $D_{\text{train}}$ and then predicts on test data $x_{\text{test}}$, as illustrated in Figure 1. Initially, a transformer estimates shape functions using ICL on the training dataset $D_{\text{train}}$. Next, predictions are computed by aggregating the shape function values for each test data point $x_{\text{test}}$. This methodology replaces the traditional data fitting process of GAM variants with a single forward pass of a pre-trained transformer model, eliminating the need for optimization and regularization hyperparameters. We now describe each model component in more detail.

### 3.1 Shape Estimation and Predictions

We obtain the shape functions with ICL by applying a transformer on the training input points and labels:

$$ $\tilde{f}=\mathcal{T}_{\theta}(x_{\text{train}},y_{\text{train}})\in\mathbb{R} ^{p\times n_{\text{bins}}\times m},$ (4) $$

where $p,m$ and $n_{\text{bins}}$ are respectively the numbers of features, classes and bins. To get predictions on a new test point $x_{\text{test}}$, we first bin each feature value and then apply the estimated shape function:

$$ $g\left(\tilde{y}_{\text{test}}\right)=\sum_{i=1}^{p}\tilde{f}_{ij_{x_{i}}}\in \mathbb{R}^{m},$ (5) $$

where $j_{x_{i}}\in[n_{\text{bins}}]$ denotes the bin index corresponding to the $i-$th feature of $x_{\text{test}}$. We now give more details on the binning and the architecture used for $\mathcal{T}_{\theta}$ in Eq. 4 before discussing the pre-training.

### 3.2 Model Architecture

#### Feature Preprocessing.

Prior to being passed through the transformer, all features of each data point are binned, one-hot encoded, and finally embedded using an MLP. We use $n_{\text{bins}}=64$ bins for each feature, allocating bins based on the quantiles of the feature in the training dataset. Similarly to TabPFN, we embed the label of each datapoint and add it to the embedding of each feature. Categorical features are equally distributed across the 64 bins according to their ratios.

#### Representation of the shape functions.

To accurately represent the shape functions, we chose to predict a discrete representation for each feature by discretizing it into 64 bins. An alternative approach would have been to predict the weights of a Neural Additive Model (NAM), similar to the method employed by Mothernet (Müller et al., ). However, we decided against this approach to more naturally represent sudden discontinuities in the shape functions(^1^11We refer to our case study on MIMIC-II for an illustration of this effect.).

#### Transformer.

The preprocessed training datapoints are processed by a transformer architecture consisting of 12 layers, each with a dual-module design that sequentially applies self-attention—first over the features and then over the data points. This design, inspired by (Lorch et al., ), ensures that our model is agnostic to the number of features and data points, and is equivariant with respect to the order of both. As a result, unlike TabPFN (Hollmann et al., ), our approach does not require padding to a fixed maximum number of features.

After the transformer layers, we compute the average embeddings for each class based on training labels enabling multi-class classification (limited to 10 classes in our experiments). This averaging yields one embedding per class per feature which we denote $h\in\mathbb{R}^{p\times d\times m}$ where $d$ denotes the embedding dimension of the transformer(^2^22This embedding is equivariant with respect to input features but not invariant to class ordering due to distinct class encodings in the input layer.). Each embedding is then passed through a shared decoder MLP to produce the binned shape functions $\tilde{f}\in\mathbb{R}^{p\times n_{\text{bins}}\times m}$. This architecture is parameter-efficient as it allows sharing of parameters across features and classes. The model comprises 40k parameters in the encoder layer, 50.5M parameters in the transformer layers, and 0.3M parameters in the decoder, resulting in a total of 50.8M parameters. Note that while the shape function estimation scales quadratically in the number of features and datapoints, the inference only scales linearly in both.

### 3.3 Training Procedure

We train with SGD on synthetic data priors, a method introduced in Prior-Data Fitted Networks (PFNs) (Müller et al., ; Hollmann et al., ). These priors are designed to be diverse, facilitating the generation of realistic tabular datasets and enabling extrapolation to real-world data. We utilize two types of priors for training: (1) Structural Causal Models, which involve sampling random causal graphs and generating data from them, and (2) Gaussian Processes, where random Gaussian Processes are sampled and used to generate data. For more details on the synthetic data generation process, we refer to Appendix D. During training, the synthetic data is randomly split into train and test datasets. To obtain the parameters $\theta$ of Eq. 4 we minimize a cross-entropy loss between the estimated GAM prediction and ground truth labels on the test dataset $D_{\text{test}}$:

$$ $\theta^{*}\in\text{argmin}_{\theta}\mathbb{E}_{(D_{\text{train}}\cup(x_{\text{ test}},y_{\text{test}}))\sim p(D)}\left[\mathcal{L}(\tilde{y}_{\text{test}},y_ {\text{test}})\right]$ (6) $$

Additional details on the training are given in Appendix E.

GAMformer’s core contribution is the substitution of the data fitting process of traditional GAM variants with a single forward pass of a pre-trained transformer model, which is presented with data through in-context examples. Consequently, GAMformer replaces the manually crafted fitting procedures used in methods like EBMs (Caruana et al., ), where the boosting procedure is restricted to one feature at a time in a round-robin manner, or the joint optimization of all shape functions in NAMs (Agarwal et al., ) using SGD. Note that in both traditional GAM fitting and GAMformer, the output of the processes remains the same; a main effects GAM fitted to a given dataset represented by its shape functions.

### 3.4 Higher-order effects

We now describe how GAMformer can be extended to handle higher-orders effects. We extend GAMformer to model higher-order effects, specifically pairwise interactions, by incorporating feature products, resulting in up to $\mathcal{O}(p^{2})$ potential features. GAMformer can accommodate this by performing ICL on concatenated original data and higher-order effects, represented as feature vectors in $\mathbb{R}^{p+P}$, where $P$ denotes the number of pair interactions. However, increasing feature dimensions beyond the $10$ used in pretraining is problematic and adds complexity to shape function estimation. To mitigate this, we rank the most informative pairs via the FAST method (Lou et al., ) and the optimal number of pairs is determined as a hyperparameter through cross-validation during inference.

## 4 Experiments

After pretraining GAMformer on the synthetic datasets, we evaluate it on both illustrative and real-world tasks in 4.1 and 4.2, respectively. Moreover, in 4.3, we highlight its potential to assist in decision-making in a clinical setting by predicting the mortality rate of patients in the intensive care unit (ICU). We compare to Explainable Boosting Machines (EBMs) (Lou et al., ; ; Caruana et al., ) in terms of estimated shape function quality, as well as to other state-of-the-art tabular classification models such as XGBoost (Chen and Guestrin, ) and TabPFN (Hollmann et al., ) in terms of predictive performance. On the downstream datasets, differently from EBM and the other baselines, GAMformer requires *only a single forward pass* of the transformer model to estimate the shape functions and construct prediction on the entire test set, without any parameter updates.

### 4.1 Illustrative Examples

Figure: Figure 2: Shape functions derived from GAMformer and EBMs applied to the linear, binary classification problem $f(x_{1},x_{2},x_{3})=\mathbb{I}((-1)x_{1}+0x_{2}+x_{3}>0)$. We use a twin y axis with GAMformer and EBM on left and right, respectively. All models shown result from a 30-fold cross-validation over 1500 data points.
Refer to caption: x2.png

Before demonstrating GAMformer on real-world tabular data, we first investigate its behavior on synthetic data where the data-generation process is known. This allows us to validate the effectiveness of GAMformer in capturing the underlying relationships between features and the target variable. All considered examples are binary classification and hence we only show one shape function per class and per feature. In the context of GAMs with a logit link function (used for binary classification), log-odds is the unit of the predictors. Therefore, the shape functions’ output values are on the log-odds scale, which are then transformed to overall prediction probabilities after summing via the logistic function. For all metrics reported in the paper, we use ROC-AUC (Receiver Operating Characteristic - Area Under the Curve).

Figure: Figure 3: Robustness analysis (linear, binary classification): GAMformer consistently outperforms or matches EBM across various sample sizes and feature counts, showcasing its efficiency
Refer to caption: x5.png

Linear, binary classification. We begin by evaluating GAMformer and, for comparison, EBMs on data generated by the linear, binary classification problem $f(x_{1},x_{2},x_{3})=\mathbb{I}((-1)x_{1}+0x_{2}+x_{3}>0)$, where $\mathbb{I}$ is the indicator function. We sample 2000 data points uniformly and independently from the interval [-2, 2] and split the data into 1500 training points and 500 test points. The results, shown in Figure 2, demonstrate that both GAMformer and EBMs accurately estimate the slopes for each feature and achieve an ROC AUC of 1.0 on the test dataset. However, the shape functions learned by GAMformer are noticeably smoother, suggesting that it may have captured some bias towards smoother models during pretraining. Additionally, we compared the effect of varying the number of datapoints or features in this example on EBMs and GAMformer in Figure 3. Our findings indicate that GAMformer consistently outperforms EBMs across various sample sizes and feature counts.

Figure: Figure 4: (a) Shape functions derived from GAMformer and EBMs applied to the polynomial, binary classification problem $f(x_{1},x_{2})=\mathbb{I}(x_{1}+x_{2}^{2}>0)$. All models result from a 30-fold cross-validation over 1500 data points are shown.
Refer to caption: x6.png

Polynomial, binary classification. To further validate the robustness of GAMformer, we evaluate it on data generated by a more complex function $f(x_{1},x_{2})=\mathbb{I}(x_{1}+x_{2}^{2}>0)$. The experimental setup remains the same as for the logistic regression case.
The results, presented in Figure 4, show that both GAMformer and EBMs successfully capture the quadratic relationship in $x_{2}$ and the linear contribution of $x_{1}$ up to $x_{1}\leq 0$. For $x_{1}>0$, $f$ always predicts true, resulting in a constant contribution. Consistent with the previous experiment, GAMformer produces smoother shape functions. Again both models achieve an ROC AUC of 1.0 on the test dataset.

Classification Boundaries.
We visualize the classification boundaries of GAMformer compared to TabPFN and EBM on the scikit-learn (Pedregosa et al., ) test datasets in Figure 5. We find that GAMformer performs similarly to TabPFN and EBMs on most of the example datasets. LA-NAM (Bouchiat et al., ) (main effects only), a Bayesian version of NAMs (Agarwal et al., ), provides good uncertainty estimates despite exhibiting slightly worse predictive performance. It is worth noting that GAMformer, EBM and LA-NAM struggle with accurately modeling the ‘XOR’ dataset (bottom row) due to the absence of higher-order feature interaction terms in these models. This is resolved by incorporating second-order effects (EBM^∗ and GAMformer^∗; see Section 3.4 for details), allowing them to effectively learn the non-linear decision boundary of the ‘XOR’ function.

Figure: Figure 5: Visualization of classification boundaries for various baseline classifiers and GAMformer on scikit-learn dataset examples (Pedregosa et al., ), in the lower right corner we show the ROC-AUC on a validation split. Due to the absence of higher-order feature interaction terms in both GAMformer and EBM (main effects), the ’XOR’ dataset (bottom row) is not accurately modeled by them. Incorporating second-order effects solves the problem (EBM^∗ and GAMformer^∗).
Refer to caption: x8.png

### 4.2 Multi-class Classification on OpenML Tabular Datasets

To assess the transferability of pretraining on synthetic data to real-world tabular data, we evaluate GAMformer’s performance on the test datasets from TabPFN (Hollmann et al., ), which include up to 2000 datapoints (see Section B.1 for dataset details).
Figure 6 reports Critical Diagrams (CD) from Demšar () showing the average rank across datasets for each method, with statistically tied methods grouped by horizontal bars. Our method outperforms EBM when using only main effects. With pair effects, both GAMformer* and EBM* show slight improvements, matching XGBoost’s performance. We also compare against GAMs from the mgcv R (^3^33[https://www.rdocumentation.org/packages/mgcv/versions/1.9-1/topics/gam](https://www.rdocumentation.org/packages/mgcv/versions/1.9-1/topics/gam)) library. mgcv GAM models the relationships between features and output variables by combining parametric and non-parametric terms. The non-parametric components are represented by splines, thus capturing nonlinear relationships. In mgcv GAM the degree of smoothness in every spline is automatically selected using Restricted Maximum Likelihood (REML) (Wood, ).

We note that the small difference in performance between XGBoost and GAMformer suggests that the trade-offs in model capacity when choosing a main effects only GAM are often less significant than expected. As a result, the substantial interpretability benefits offered by the GAM model class become even more appealing, making it a viable choice for many applications. We present additional results on five binary classification datasets used by Chang et al. () in Section B.2. Despite these datasets falling outside the recommended range of 2000 datapoints, GAMformer still demonstrates comparable performance to more complex models.

### 4.3 Case Study: Intensive Care Unit Mortality Risk

In this case study, we examine shape functions derived from GAMformer and EBMs (main effects only) using the MIMIC-II dataset (Lee et al., ), a publicly available critical care dataset for predicting mortality risk based on various demographic and biophysical indicators. Our analysis focuses on four key clinical variables: Age, Heart Rate (HR), PFratio (PaO2/FiO2 ratio), and Glasgow Coma Scale (GCS), as shown in Figure 7 (remaining variables in Section G.1). Further results on the MIMIC-III dataset are available in Section G.2.

Figure: Figure 6: Critical Difference diagram demonstrating GAMformer’s competitive performance against state-of-the-art baselines across diverse datasets. Lower ranks indicate superior performance; connected algorithms are not statistically significantly different ($p=0.05$).
Refer to caption: x9.png

Figure: Figure 7: Shape functions derived from GAMformer and EBMs applied to the MIMIC-II dataset for critical clinical variables. The data density plot is shown above each figure. The results are based on 30 models for both GAMformer and EBMs, each fitted on 10,000 randomly selected data points.
Refer to caption: x10.png

For Age, the GAMformer shape function shows a steady increase in the log-odds of adverse outcomes with advancing age, stabilizing at older ages. The data density plot reveals a higher concentration of data points in middle age, with fewer at the extremes. The shape function exhibits less variance where data is denser, indicating the model’s reliability in these regions. Overall, the shape function highlights increased risk in elderly patients due to declining physiological reserves and multiple chronic conditions. Heart Rate (HR) exhibits a complex relationship with adverse outcomes. Both GAMformer and EBMs capture a U-shaped risk profile, indicating increased risk at very high and very low heart rates, underscoring the importance of maintaining HR within a normal range. PFratio, a lung function and oxygenation efficiency measure, shows a steep risk increase as values decrease. Lower PFratio values, critical in diagnosing and managing conditions like Acute Respiratory Distress Syndrome (ARDS), indicate worse lung function. Notably, both models display a sharp drop in risk at a PFratio of approximately 325, likely an artifact from data preprocessing where missing values were imputed at the mean, previously pointed out by Chen et al. () for MIMIC-2. In healthcare, missing values often suggest healthier patients, as data collection was deemed unnecessary by professionals. Here, patients with missing PFratio values, representing the majority, have lower risk than those with collected values. GAMformer more precisely isolates these missing value patients, demonstrating its potential to detect data processing artifacts better than prior GAM algorithms. For the Glasgow Coma Scale (GCS), which measures the level consciousness, there is a strong negative correlation with adverse outcomes. Lower GCS scores, indicating reduced consciousness, are associated with significantly higher mortality risk. Our findings show that GAMformer effectively handles categorical data, identifying patterns similar to those detected by EBMs.

## 5 Limitations & Broader Impact

Limitations. While GAMformer introduces a novel approach to estimating Generalized Additive Models (GAMs), it is important to acknowledge its current limitations. This work primarily focuses on main and second-order effect GAMs and does not account for higher-order interactions, which are addressed in other GAM implementations, such as EBMs (Lou et al., ; Nori et al., ; Chang et al., ). Future research could explore incorporating these interactions to enhance the model’s expressiveness and predictive capabilities.
Another limitation of the current GAMformer model is its difficulty in improving predictions when presented with datasets that exceed twice the size of the data it saw during training (c.f. Figure 8). This issue is related to the well-known challenge of length extrapolation in sequence-to-sequence models, including transformers (Grazzi et al., ; Zhou et al., ). Addressing this limitation may require exposing the model to a larger variety of number of examples during training. However, due to computational constraints, the experiments in this work were limited to a maximum of 500 datapoints during training. Future studies with increased computational resources could investigate the model’s performance on larger datasets and develop strategies to mitigate the length extrapolation problem.
The GAMformer model’s transformer-based architecture scales quadratically with both the number of training data points and features, posing a similar challenge to handling large datasets as faced by TabPFN (Hollmann et al., ). Novel, scalable transformer alternatives, such as the recently proposed Mamba (Gu and Dao, ) or Gated Linear Attention (Yang et al., ), may prove useful in overcoming this issue.

Broader Impact. As a versatile machine learning model for tabular data, GAMformer offers both positive and negative societal impacts. Positively, it can generate novel insights in fields like medicine, enhancing disease diagnosis and treatment. However, it can also be misused to not mitigate but exploit biases, such as adjusting insurance premiums based on ethnicity, leading to discrimination.

## 6 Conclusion

In this paper, we introduce GAMformer, a novel approach to creating GAMs using in-context learning with transformer models. By leveraging a single forward pass to form shape functions, GAMformer overcomes the limitations of traditional GAM algorithms that require iterative learning processes and hence hyperparameter tuning. Our approach uses non-parametric, binned representations of shape functions, resulting in significant improvements in efficiency and accuracy across various classification benchmarks. Extensive experiments demonstrate that GAMformer approaches the accuracy of leading GAM variants while exhibiting robustness to label noise and class imbalance. The model’s ability to generalize beyond the number of examples seen during training highlights its adaptability and potential for practical applications.

GAMformer is fundamentally different from the iterative optimization methods traditionally used, and offers a new research direction for interpretable models on tabular data. Further, our case study on the MIMIC-II dataset showcases that interpreting GAMformer’s shape functions can yield qualitative insights and uncover flaws in datasets similar to state of the art GAM methods. This work contributes to the development of more transparent and explainable AI systems, with potential applications in various domains where interpretability is crucial. Future research can expand on this initial new paradigm, and explore scalable alternatives to transformers to handle larger datasets.

## Acknowledgments

This research was partially supported by the following sources: TAILOR, a project funded by EU Horizon 2020 research and innovation programme under GA No 952215; the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant number 417962828; the European Research Council (ERC) Consolidator Grant “Deep Learning 2.0” (grant no. 101045765). Frank Hutter acknowledges financial support by the Hector Foundation. The authors acknowledge support from ELLIS and ELIZA. Funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the ERC. Neither the European Union nor the ERC can be held responsible for them.

.

## Appendix A Generalized Additive Models: Extended Related Work

As with many families of machine learning algorithms, the differences among GAM algorithms lie in (a) the functional form of the shape functions $f_{i}$, (b) the learning algorithm used for their estimation and (c) regularity assumptions and regularization. Two important properties that all GAMs share are (1) the ability to learn non-linear transformations for each feature and (2) additively combining these shape functions (prior to applying the link function) to create modularity that aids interpretability by allowing users to examine shape functions one-at-a-time.

Typically, GAMs have relied on splines and backfitting algorithms for estimation [Hastie and Tibshirani, ], with subsequent works focusing on improving efficiency and stability through penalized regression splines [Wood, ] and fast, stable fitting algorithms [Wood, ]. Spline-based GAMs are typically fitted using the backfitting algorithm, an iterative procedure that starts with initial estimates of the smooth functions for each predictor variable. The algorithm then repeatedly updates each function by fitting a weighted additive model to the residuals of the other functions until convergence is achieved. The weights are determined by the current estimates of the other functions and the link function in the case of generalized additive models.

Modern approaches leverage machine learning advances. Explainable Boosting Machines (EBMs) [Lou et al., , , Caruana et al., ] model the shape functions using decision trees, which are fitted using a variant of gradient boosting called cyclic gradient boosting. The model iteratively learns the contribution of each feature and interaction term in a round-robin fashion, using a low learning rate to ensure that the order of features does not affect the final model. This cyclic training procedure helps mitigate the effects of colinearity among predictors by providing opportunity for data-driven credit attribution among the features while preventing multiple counting of evidence. EBMs are also popular because they can accurately capture steps in the shape functions, which is important for modeling discontinuities in data, such as treatment effects in medical data.

More recently, Neural Additive Models (NAMs) [Agarwal et al., ] and follow up works [Chang et al., , Dubey et al., , Radenovic et al., , Xu et al., , Enouen and Liu, , Bouchiat et al., ] use multilayer perceptrons (MLPs), as non-linear transformations, to model the shape functions $f_{i}$. As a result, NAMs can be optimized using variants of gradient descent by leveraging automatic differentiation frameworks.

Finally, GAMs have also found applications in time-series forecasting, with models such as Prophet [Taylor and Letham, ] and NeuralProphet [Triebe et al., ]. Interestingly, the 1-layer versions of the recently proposed Kolmogorov-Arnold Networks (KANs) [Liu et al., ] may be viewed as GAMs with spline based shape functions.

## Appendix B Dataset Details

In this section, we provide details on the datasets used in our empirical evaluations of GAMformer and other baselines in Section 4 of the main paper.

### B.1 TabPFN Test Datasets

As test dataset, we used the 30 datasets used in Hollmann et al. [] which were obtained from OpenML [Vanschoren et al., ]. These were chosen because they contain up to 2000 samples, 100 features and 10 classes, show in Table 1.

**Table 1: Test dataset names and properties, taken from Hollmann et al. []. Here *did* is the OpenML Dataset ID, *d* the number of features, *n* the number of instances, and *k* the number of classes in each dataset.**
| did | name | d | n | k |
| --- | --- | --- | --- | --- |
| 11 | balance-scale | 5 | 625 | 3 |
| 14 | mfeat-fourier | 77 | 2000 | 10 |
| 15 | breast-w | 10 | 699 | 2 |
| 16 | mfeat-karhunen | 65 | 2000 | 10 |
| 18 | mfeat-morphological | 7 | 2000 | 10 |
| 22 | mfeat-zernike | 48 | 2000 | 10 |
| 23 | cmc | 10 | 1473 | 3 |
| 29 | credit-approval | 16 | 690 | 2 |
| 31 | credit-g | 21 | 1000 | 2 |
| 37 | diabetes | 9 | 768 | 2 |
| 50 | tic-tac-toe | 10 | 958 | 2 |
| 54 | vehicle | 19 | 846 | 4 |
| 188 | eucalyptus | 20 | 736 | 5 |
| 458 | analcatdata_authorship | 71 | 841 | 4 |
| 469 | analcatdata_dmft | 5 | 797 | 6 |

### B.2 Binary Classification

**Table 2: Comparison of GAMformer with other GAM variants and full complexity models on various datasets. We report ROC-AUC (%) (higher is better) and the standard error over 10 fold cross-validation. We also report results by pyGAM [Servén and Brummitt, ].**
|  | GAMs | Full Complexity |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
|  | GAMformer (ours) | EBM (Main effects) | Logistic Regression | pyGAM (Main effects) | EBM | XGBoost | Random Forest |
| Churn | 81.69 $\pm$ 0.1 | 83.59 $\pm$ 0.1 | 81.66 $\pm$ 0.1 | 82.03 $\pm$ 0.0 | 83.68 $\pm$ 0.1 | 83.53 $\pm$ 0.0 | 82.07 $\pm$ 0.0 |
| Support2 | 80.84 $\pm$ 0.1 | 82.36 $\pm$ 0.0 | 81.1 $\pm$ 0.0 | 81.74 $\pm$ 0.2 | 83.51 $\pm$ 0.0 | 84.03 $\pm$ 0.0 | 83.93 $\pm$ 0.0 |
| Adult | 90.05 $\pm$ 0.0 | 93.05 $\pm$ 0.0 | 90.73 $\pm$ 0.0 | 91.55 $\pm$ 0.0 | 93.07 $\pm$ 0.0 | 93.16 $\pm$ 0.0 | 91.8 $\pm$ 0.0 |
| MIMIC-2 | 82.22 $\pm$ 0.0 | 85.15 $\pm$ 0.0 | 81.62 $\pm$ 0.0 | 83.89 $\pm$ 0.1 | 86.36 $\pm$ 0.1 | 87.29 $\pm$ 0.0 | 87.31 $\pm$ 0.0 |
| MIMIC-3 | 74.41 $\pm$ 0.1 | 81.14 $\pm$ 0.0 | 78.05 $\pm$ 0.0 | 79.95 $\pm$ 0.1 | 82.52 $\pm$ 0.1 | 83.32 $\pm$ 0.0 | 81.28 $\pm$ 0.1 |

#### Churn dataset.

The [Telco Customer Churn Dataset](https://www.kaggle.com/datasets/blastchar/telco-customer-churn) is a binary classification dataset for predicting potential subscription churners in a telecom company, containing customer information and churn-related features.

#### Adult dataset.

The Adult dataset Dua and Graff [], also known as the “Census Income” dataset, is a widely-used benchmark for binary classification, predicting whether an individual’s annual income exceeds $50,000 based on 14 attributes from the 1994 United States Census Bureau data.

#### MIMIC-II dataset.

The MIMIC-II dataset Lee et al. [] is a publicly-available database of clinical data from diverse ICU patients, integrating demographics, vital signs, lab results, medications, procedures, notes, and imaging reports, along with mortality outcomes.

#### MIMIC-III dataset.

The MIMIC-III dataset Johnson et al. [] expands on MIMIC-II, with a larger patient cohort, more recent records, enhanced data granularity, and the inclusion of free-text imaging report interpretations.

#### SUPPORT2 dataset.

The SUPPORT2 dataset Connors Jr et al. [] contains medical information from critically ill hospitalized adults, compiled to study the relationships between medical decision-making, patient preferences, and treatment outcomes, with variables spanning demographics, physiology, diagnostics, treatments, and survival/quality of life outcomes.

## Appendix C Properties of GAMformer

### C.1 Data Scaling

To assess GAMformer’s ability to generalize to datasets containing more datapoints than it saw during training, i.e. larger context sizes, we conducted an experiment that varied the number of training data points and evaluated the impact on ROC-AUC performance using a consistent validation split. To ensure the robustness of our findings, we sampled training datasets three times with replacement for each training size. The results in Figure 8 demonstrate that GAMformer’s ROC-AUC improves across datasets when the number of training examples is up to twice the number of training examples seen during training. For comparison, we also evaluated the performance of EBMs under the same conditions. While EBMs also exhibited improvements in ROC-AUC with increased training data, they achieved higher accuracy when provided with a larger number of examples. This observation highlights a limitation of GAMformer in its ability to fully leverage additional training samples.

Figure: (a) MIMIC-II
Refer to caption: x14.png

### C.2 Class Imbalance

Figure: (a) Imbalanced Data
Refer to caption: x18.png

To compare GAMformer’s sensitivity to class imbalance with that of EBMs, we conduct the following analysis. First, we sample 300 data points from two centroids in a 20-dimensional feature space, creating a binary classification problem. We then vary the ratio of the two classes to introduce increasing levels of imbalance in the sampled data. Next, we split the data into train and test sets using a 75% to 25% split and evaluate the performance using the AUC-ROC metric. We repeat the experiment 10 times for each data ratio. Our results are shown in Figure 9(a), the shaded area are the 5%, 95% confidence intervals estimated using 1000 bootstrap samples. We see that GAMformer performs on average better than EBMs in this setting and shows no inherent sensitivity to class imbalance.

### C.3 Noise Robustness

To gain a deeper understanding of GAMformers’ sensitivity to noisy or incorrect labels, we conducted an experiment similar to the one described in Section C.2. We generated 300 data points and randomly perturbed the labels in the train split with increasing probability (75%, 25% train/test split), repeating each experiment 10 times. Figure 9(b) illustrates our findings. Once again, we observed that GAMformer exhibits a sensitivity to noisy labels comparable to that of EBMs.

## Appendix D Synthetic Data Priors

We use the same synthetic data generation process proposed in Prior-Data-Fitted Networks (PFNs) [Hollmann et al., , Müller et al., ] and provide a brief summary of the process.

TabPFN is trained on two synthetic data priors, which are mixed during training.TabPFN introduced a synthetic data prior based on Structural Causal Models (SCMs). SCMs are particularly suitable for modeling tabular data as they capture causal relationships between columns, a strong prior in human reasoning. An SCM comprises a set of structural assignments (mechanisms) where each mechanism is defined by a deterministic function and a noise variable, structured within a Directed Acyclic Graph (DAG). The causal relationships are represented by directed edges from causes to effects, facilitating the modeling of complex dependencies within the data.
To instantiate a PFN prior based on SCMs, one defines a sampling procedure to create supervised learning tasks. Each dataset is generated from a randomly sampled SCM, including its DAG structure and deterministic functions. Nodes in the causal graph are selected to represent features and targets, and samples are generated by propagating noise variables through the graph. This process results in features and targets that are conditionally dependent through the DAG structure, capturing both forward and backward causation [Hollmann et al., ]. This allows for the generation of diverse datasets.

The second prior samples of synthetic data using Gaussian Processes (GPs) [Rasmussen and Williams, ] with a constant mean function and a radial basis function (RBF) kernel to define the covariance structure. Hyperparameters such as noise level, output scale, and length scale are sampled from predefined distributions to introduce variability. Depending on the configuration, input data points can be sampled uniformly, normally, or as equidistant points and the target column is generated by passing the input data through the GP. This prior gives the model the ability to learn smoother functions.

For multi-class prediction, scalar labels are transformed into discrete class labels by partitioning the scalar values into intervals corresponding to different classes, ensuring the synthetic data is suitable for imbalanced multi-class classification tasks.

Finally, both priors are combined by sampling batches of data from each prior with different probabilities during training. In all of our experiments we sampled from the SCM and GP prior with probability $0.96$ and $0.04$, respectively.

## Appendix E Training Details

In GAMformer, we used a transformer model with 12 hidden layers, 512 embedding size and 4 heads per attention. To bin the shape functions and all features we used 64 bins.
For training, we use the AdamW [Loshchilov and Hutter, ] optimizer ($\beta_{1}=0.9$) and cosine learning rate schedule with initial learning rate of 3e-5, 20 warm up epochs and minimum learning rate of 1e-8 for 25 days on a A100 GPU with 80Gb of memory. We used mixed precision training. Each epoch (arbitrarily) consists of 65536 synthetic datasets; the model trained for 1800 epochs, meaning it saw over 100M synthetic datasets. We used a batch size of 8, that we doubled at epoch 20, 50, 200 and 1000. Each synthetic dataset consisted of 500 samples that were split into training and test portions using using a uniform sampling of the training fraction, and used a number of features drawn uniformly between 1 and 10.

## Appendix F Higher-order effects

To handle higher-order effects, we compute the best pairs with the FAST algorithm [Lou et al., ] and evaluate GAMformer on the top pairs using the following ratios of features:

$$ $\mathcal{P}=[0.01p,0.05p,0.1p,0.2p,0.4p,0.8p,0.9p]$ $$

where we recall that $p$ denotes the number of features. We round off each ratio to determine the number of target pair features, evaluate performance on hold-out validation data from the training set, and select the number of pairs with the best validation performance. The model is then fitted on the entire training dataset. This involves doing $|\mathcal{P}|+1$ forward passes, which is unproblematic as doing one forward pass is very fast, even on a CPU. One could also vectorialize all computations which we do not do given the low fitting time.

## Appendix G Shape Functions

In this section, we show complementary results on the shape functions estimates from GAMformer and EBM (main effects only) on the MIMIC-II [Lee et al., ] (complementary to the plots in Figure 7) and on the MIMIC-III datasets.

### G.1 MIMIC-II dataset

Figure: Figure 10: The remaining shape functions derived from GAMformer and EBMs on the MIMIC-II dataset for critical clinical variables. The plot above each figure shows the data density. There are interesting differences between the EBM and GAMformer shape plots for several of the categorical variables. Although different GAM algorithms do not usually learn identical functions, we are investigating to better understand these differences.
Refer to caption: x20.png

### G.2 MIMIC-III dataset

Figure: Figure 11: The shape functions derived from GAMformer and EBMs on the MIMIC-III dataset for critical clinical variables. The plot above each figure shows the data density. The results are based on 30 models for both GAMformer and EBMs, each fitted on 10,000 randomly selected data points. There are interesting differences between the EBM and GAMformer shape plots for several of the categorical variables. Although different GAM algorithms do not usually learn identical functions, we are investigating to better understand these differences.
Refer to caption: x33.png

Figure: Figure 12: The remaining shape functions derived from GAMformer and EBMs applied to the MIMIC-III dataset for critical clinical variables. The plot above each figure shows the data density in the training set. The results are based on 30 models for both GAMformer and EBMs, each fitted on 10,000 randomly selected data points.
Refer to caption: x61.png

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: Please see the weaknesses section.
soundness: 3
strengths: GAMformer is a contribution to GAMs, leveraging ICL and transformer models to eliminate iterative optimization, thereby simplifying the modeling process and reducing the computational overhead associated with traditional GAMs.

The model maintains high interpretability—crucial for critical fields like healthcare—while matching the performance of established methods like Explainable Boosting Machines (EBMs).

GAMformer’s training on synthetic data enables it to generalize to real-world data effectively, a challenging task for many models, especially in interpretability-driven applications.

The use of a non-parametric, binned representation for shape functions allows for flexibility, particularly for capturing discontinuities or sudden shifts in feature impacts.

The model was rigorously tested across various benchmark datasets, and a case study on ICU mortality in the MIMIC-II dataset demonstrated its clinical interpretability potential, which is well-aligned with the paper’s goals.
summary: The paper presents GAMformer, a novel model for fitting Generalized Additive Models (GAMs) using in-context learning (ICL) within a transformer-based framework. Unlike traditional GAMs that rely on iterative methods such as splines or gradient boosting, GAMformer uses a single forward pass to estimate shape functions for each feature, eliminating the need for hyperparameter tuning and iterative optimization. This approach is trained exclusively on synthetic data but performs competitively with existing GAMs on real-world tasks. GAMformer’s non-parametric, binned approach to shape function estimation enables high interpretability of feature impacts. Experimental results show that GAMformer matches or surpasses other interpretable machine learning methods on both synthetic and real-world tabular datasets, including clinical applications on the MIMIC dataset for ICU mortality prediction. Additionally, the model’s adaptability to real-world data demonstrates its potential for scalable, interpretable applications without extensive tuning.
weaknesses: GAMformer currently only supports main effects and second-order feature interactions, limiting its applicability for datasets where higher-order interactions are significant. 

The Transformer architecture in GAMformer scales quadratically with the number of data points, leading to potential performance bottlenecks for very large datasets. Exploring scalable attention mechanisms, as the authors suggest, would strengthen the model’s practical use.

While the clinical case study is insightful, further empirical evaluations across diverse fields (e.g., finance, manufacturing) would provide a clearer picture of GAMformer’s interpretability and performance across different domains.

There is a lack of quantitative results tables comparing the model with recent baselines.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 4
questions: Do you only consider discrete target variables in the experiments? Given that the features are binned and discretized, could the method be applied for regression with continuous variables?
soundness: 3
strengths: The method appears to be a novel approach for learning generalized additive models.

The paper is well-written, ideas and goals are clearly stated, background work is acknowledged and limitations are addressed.

Experiments are done on synthetic and real examples with an extensive public health case study interpreting the learned shape functions and their implications. 

The paper discusses the limitations of the model which are 1) lack of accounting of higher-order interactions 2) lack of improvement over datasets larger than seen during training and 3) quadratic complexity of the transformer.

Also propose an extension to model higher-order effects by concatenating data and high-order effects.
summary: The paper proposes an in-context learning approach for learning generalized additive models for tabular data building on prior work (PFN and TabPFN) for tabular classification.
The training procedure executes on synthetic data by sampling a random causal graph and generating data from an initial random sample. The data is split into training and test datasets to simulate inference. 
A transformer model applies attention across the data points and features and handles tabular data of varying sizes. A single forward pass of the transformer estimates the shape functions for the given in-context training data which are then applied to the test example. The shape functions themselves are represented as discrete functions which apply to discretized and binned features. The method is demonstrated experimentally on synthetic and real data including a mortality risk case study where the shape functions are used to interpret model predictions.
weaknesses: The approach appears to be limited to discrete target values. Shape functions are learned as discretized functions over discretized features which could be limiting.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: Please, refer to the main weaknesses.
soundness: 2
strengths: 1. The paper is clear and well-written **Clarity**
2. The paper addresses an important and relevant problem, specifically how to leverage deep learning to learn a prior for predictive tasks on tabular data. **Relevance**
3. The code is available, and a Jupyter Notebook is provided to demonstrate how the proposed model and explainable boosted machines generate the predictive functions **Code Availability**. However, no checks have been performed to verify the reproducibility of the experiments.
summary: The paper addresses the problem of supervised learning for tabular data and proposes a solution based on generalized additive models. A key feature is the use of an attention-based neural network (Transformer) to process the training data and provide a prior over the parameters of the non-linear predictive functions. The learning process involves splitting the training data into a training set and a holdout set. A predictive likelihood over the holdout set is used to learn the prior based on the training set. Experiments on synthetic data and OpenML datasets are conducted to compare the proposed solution with explainable boosted machines, demonstrating its ability to achieve comparable predictive performance
weaknesses: 1. The novelty of the paper is limited and incremental. **Novelty**. The main ideas have already appeared in two previous works [1,2], and the primary difference seems to be the use of a different classifier/regressor. In other words, instead of considering Bayesian neural networks or structural causal models like in [1,2], the authors focus on generalised additive models. In essence, the work can be seen as an application of existing ideas within the context of generalised additive models.
2. There are several vague and overstated claims that are not properly supported. For instance, the abstract mentions that the proposed solution generates highly interpretable predictive functions. **Soundness** However, this is also true for the competitors, and it is unclear what the real advantage of the proposed solution is over existing generalised additive models and other interpretable models (such as XGBoost). In the experiments (e.g., lines 304-305), it is stated that the proposed solution outperforms explainable boosted machines (EBMs), but these claims seem exaggerated. Firstly, in the low-data regime (32 samples) with a larger number of features (64), the proposed solution clearly underperforms compared to EBMs by 14 points, suggesting a possible blind spot and indicating that sufficient data is required for the proposed solution to perform on par. Secondly, it is unclear whether the differences in the results are statistically significant, as no standard deviation is provided. Similarly, for Figures 2 and 3, it is claimed that the proposed solution clearly learns smoother predictive functions. However, this is subjective and not consistently true (only the 1st and 3rd plots in Figure 2 support the authors' claim).
3. The experimental analysis lacks a consistent comparison across datasets and tasks with other interpretable models. Additionally, the analysis focuses on the case where the ground truth classifier lies within the hypothesis space. What about the agnostic case? **Quality**
4. The experiments are conducted on small datasets, reflecting the poor scalability of the approach. While the idea of synthesizing data may be reasonable for small datasets, it may not be tractable or feasible for higher-dimensional data, given the potential for combinatorial explosion. Scalability and feasibility are currently overlooked, which is a significant limitation of the proposed solution. As a result, it is unclear why one should prefer this approach over existing interpretable models that are more scalable. **Quality/Significance**

**References** \
[1] Müller, Hollmann, Arango, Grabocka, Hutter. “Transformers can do Bayesian Inference”. ICLR 2022 \
[2] Hollmann, Müller, Eggensperger, Hutter. “TabPFN: A Transformer that Solves Small Tabular Classification Problems in A Second”. ICLR 2023

---

confidence: 4
contribution: 1
initial_rating: 3
presentation: 2
questions: - **Q1**: I would be very happy if the authors could address the weaknesses I have mentioned above
- **Q2**: Are there any insights of GAMformers that I might have missed?
- **Q3**: Have the authors thought about analyzing the smoothness of GAMformers and whether this could provide more interpretable functions compared to the jagged functions of NAMs / EBMs?
- **Q4**: Have the authors thought about extending the class of GAMs? I would assume that this model could also learn a combination of GAMs, trees, NODEs, etc., and still remain interpretable.
soundness: 2
strengths: - **Originality**: The idea of learning shape functions in-context is novel and interesting. 
- **Numerical Experiments** The comparisons cover a variety of models from different classes, thus giving a bigger picture of GAMformer's capabilities.
summary: The authors use prior fitted networks to train a large transform-type architecture that can then learn the shape functions of additive models in a single forward pass. The resulting model is competitive to other approaches with similar interpretability and capacity.
weaknesses: ## Major

### Contributions

- **[C1]** The claim that "experimental results demonstrate GAMformer's capacity to match the accuracy of leading GAMs" might be accurate regarding performance, but the interpretability has not been sufficiently scrutinized. See comments on Experiments below.
- **[C2]** In light of GAMs requiring no tuning (at least in the `mgcv` package), the claim "... to form shape functions ... eliminating the need for ... iterative learning and hyperparameter tuning" does not seem particularly significant.
- **[C3]** The contribution claiming the model was applied to the MIMIC-II dataset lacks significance. This dataset has been analyzed previously. The current study does not add any new insights. The dataset itself is also not particularly challenging, yet the modeling approach seems to have missed a key property of the dataset (see **E3** below).

### Technical soundness/correctness

- **[T1]** The introduction to GAMs is missing a distributional assumption (a GAM consists of both structural and distributional assumptions; see Wood, 2017).
- **[T2]** The simulated functions are not GAMs but deterministic functions. As correctly noted by the authors, a GAM is defined by a link function, yet they use a simple indicator function without induced noise or distributional assumptions for the simulation, which does not correspond to the data-generating process of a GAM. The code uses `sigmoid`, but there is no distribution involved (same for the regression task).
- **[T3]** "Allocating bins based on the quantiles of the feature in the training dataset" → This approach is likely inferior to equidistant binning, as quantile-based binning alters the data distribution of the feature (see Li and Wood, 2017).
- **[T4]** The comparisons with `mgcv::gam` appear incomplete (see below).

### Significance

- **[S1]** The computational costs of:  
    + fitting a GAM are $O(N_{train} \cdot (K \cdot p)^2)$ (see Wood, 2020), where $p$ is the number of features and $K$ the number of basis functions (in `mgcv`, often set to 10). For the data used by the authors, this would amount to 50-800 parameters;  
    + predicting with a fitted GAM is $O(N_{test} \cdot K \cdot p)$.

  In contrast, ICL requires millions of parameters, if I understand Sec. 3.2 correctly, and even with fast inference, it is slow compared to GAMs, where typically $N > p$ and hence the quadratic scaling of $N$ in the transformer is still the bottleneck. Moreover, the authors report that the model required 25 days on a high-performance GPU, whereas all the analyzed datasets could be fit within seconds using GAMs. GAMs can also be applied to datasets of size $10^8$ using `mgcv::bam` (see Wood et al., 2017).
  
- **[S2]** The method does not seem to outperform other models in prediction accuracy and appears to be inferior to TabPFN. TabPFN itself could also be analyzed using SHAP after computing the predictions, raising the question if a specific architecture is even necessary.
- **[S3]** I could not identify any other significant insights of theoretical nature or similar that could be derived from a GAMformer. In particular, I would assume that the SCMs in TabPFN likely already cover GAM-type models (related to **S2**). This raises the question of what additional benefit is gained by making them explicit, as done here.

### Experiments

- **[E1]** The experiments do not show the shape functions of the GAM method, which would be particularly useful for illustrative examples. 
- **[E2]** Simulations should be designed to correspond to an actual GAM to see whether the GAMformer can actually recover those (see **T2**).
- **[E2/T4]** The `mgcv::gam` should not be inferior to logistic regression if used correctly (see Figure 6).
- **[E3]** Isn't there censoring in the MIMIC datasets? A time-to-event model might be more appropriate in this case then.
- **[E4]** In the Appendix experiments, the authors switch to `pyGAM`, which is known to be inferior to `mgcv::gam`, and do not report the latter's performance.

### Reproducibility

- **[R1]** The code does not provide competitor models.

### Writing

- **[W1]** There is some redundancy between Sections 1 and 2, which disrupts the flow of reading.
- **[W2]** The notation $j_{x_i}$ is somewhat confusing, as $j$ is an index in the bins and $x_i$ represents the $i$th feature in $x$.

## Minor / Technical soundness

- **[M1]** The $f$ functions are typically referred to as *smooth terms* or *smooth functions* in the GAM literature, not *shape functions* (a term seemingly invented by the NAM community). They are also not *partial dependence plots* (as these are plots, not functions; in GAM literature, they are referred to as *partial effects*).
- **[M2]** The $g$ function typically does not map to $\mathbb{R}$ but to a subspace ($\mathcal{Y}$, or more specifically, e.g., (0,1) for the logistic function).
- **[M3]** What is $q_\theta$ in equation (2)?
- **[M4]** "Spline-based GAMs use the backfitting algorithm" $\rightarrow$ Backfitting was proposed by Hastie and Tibshirani. More recent approaches, like those from Wood, use PIRLS or alternatives like INLA (see Wood, 2017).
- **[M5]** The citation in footnote 3 (and for the `mgcv` package in general) seems incorrect.
- **[M6]** No "shape functions" for pairwise smooth interactions are shown.

## References

- Wood, 2017: https://www.taylorfrancis.com/books/mono/10.1201/9781315370279/generalized-additive-models-simon-wood  
- Li and Wood, 2017: https://link.springer.com/article/10.1007/s11222-019-09864-2  
- Wood et al., 2017: https://www.tandfonline.com/doi/full/10.1080/01621459.2016.1195744  
- Wood, 2020: https://www.maths.ed.ac.uk/~swood34/test-gam.pdf

## Suggestions for Improving the Paper

Here are some suggestions for improving the paper:

1. **Writing/Novelty/Significance**: 1) Clearly articulate how GAMformer advances beyond current GAM implementations and/or what additional insights they provide for PFNs or in-context learning. If you consider not changing your listed contributions, provide a more rigorous theoretical comparison with GAMs (computational complexity, etc.). I would, however, suggest changing your argumentation and thinking about what other aspects a PFN-type model can provide that a GAM cannot.

2. **Technical Correctness and Clarity**: Revise the GAM background to more formally introduce GAMs, and merge the redundant parts of Sections 1 and 2.

3. **Numerical Experiments**: Consider 1) simulating datasets that follow a GAM data-generating process, 2) comparing against ``mgcv::gam`` in performance, 3) showing shape functions also for GAMs.

4. **Application**: If the authors' approach allows the inclusion of censoring, consider modifying your application. Alternatively, consider a different and more challenging dataset.

5. **Reproducibility**: Include the code for competitor models.

Meta-review: The paper introduces GAMformer, a method that uses in-context learning to fit Generalized Additive Models in a single step, rather than traditional iterative methods. The models presented are trained on synthetic data but demonstrates good performance on real-world datasets. The authors claim competitive performance with leading GAM implementations while maintaining interpretability.

The method is limited to first and second-order feature interactions. As well, the neural network based approach loses many of the theoretical guarantees compared to traditional GAMs. Reviewers were concerned by the cost of training the model compared to ad-hoc GAM training. Reviewers also had concerns about the experimental validation and the novelty compared to recent works.

Based on the reviews and discussion, this paper appears to be marginally below the acceptance threshold. While the approach is interesting, I believe this is outweighed by reviewers' concerns so I will recommend rejection.

Decision: Reject