Paper deep dive
Interpretability-by-Design with Accurate Locally Additive Models and Conditional Feature Effects
Vasilis Gkolemis, Loukas Kavouras, Dimitrios Kyriakopoulos, Konstantinos Tsopelas, Dimitrios Rontogiannis, Giuseppe Casalicchio, Theodore Dalamagas, Christos Diou
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 91%
Last extracted: 7/21/2026, 1:34:41 AM
Summary
The paper introduces Conditionally Additive Local Models (CALMs), a novel interpretable-by-design model class that balances the interpretability of Generalized Additive Models (GAMs) with the accuracy of Generalized Additive Models with pairwise interactions (GA2Ms). CALMs achieve this by allowing multiple univariate shape functions per feature, each active in distinct regions of the input space defined by simple logical conditions (thresholds) on interacting features. A distillation-based training pipeline is proposed to identify these homogeneous regions and fit shape functions via region-aware backfitting, demonstrating superior accuracy compared to GAMs and comparable accuracy to GA2Ms while maintaining strict interpretability.
Entities (6)
Relation Signals (6)
CalM → outperforms → GAM
confidence 95% · Experiments on diverse classification and regression tasks show that CALMs consistently outperform GAMs
CalM → uses → Shape Function
confidence 95% · CALMs allow multiple univariate shape functions per feature, each active in different regions of the input space.
CalM → achievescomparableaccuracyto → GA2M
confidence 90% · CALMs consistently outperform GAMs and achieve accuracy comparable with GA2Ms.
CalM → balances → Interpretability
confidence 90% · CALMs offer a compelling trade-off between predictive accuracy and interpretability.
CalM → trainedby → Distillation Pipeline
confidence 90% · We further propose a principled distillation-based training pipeline that identifies homogeneous regions... and fits interpretable shape functions
Distillation Pipeline → usesmetric → Heterogeneity
confidence 85% · To identify such regions, we use the interaction-related heterogeneity measure... making it an effective criterion for our purpose.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Generalized additive models (GAMs) offer interpretability through independent univariate feature effects but underfit when interactions are present in data. GA$^2$Ms add selected pairwise interactions which improves accuracy, but sacrifices interpretability and limits model auditing. We propose \emph{Conditionally Additive Local Models} (CALMs), a new model class, that balances the interpretability of GAMs with the accuracy of GA$^2$Ms. CALMs allow multiple univariate shape functions per feature, each active in different regions of the input space. These regions are defined independently for each feature as simple logical conditions (thresholds) on the features it interacts with. As a result, effects remain locally additive while varying across subregions to capture interactions. We further propose a principled distillation-based training pipeline that identifies homogeneous regions with limited interactions and fits interpretable shape functions via region-aware backfitting. Experiments on diverse classification and regression tasks show that CALMs consistently outperform GAMs and achieve accuracy comparable with GA$^2$Ms. Overall, CALMs offer a compelling trade-off between predictive accuracy and interpretability.
Tags
Links
- Source: https://arxiv.org/abs/2602.16503v1
- Canonical: https://arxiv.org/abs/2602.16503v1
Trouble viewing inline? Open PDF directly →
Full Text
146,789 characters extracted from source content.
Expand or collapse full text
Interpretability-by-Design with Accurate Locally Additive Models and Conditional Feature Effects Vasilis Gkolemis Loukas Kavouras Dimitrios Kyriakopoulos Konstantinos Tsopelas Dimitrios Rontogiannis Giuseppe Casalicchio Theodore Dalamagas Christos Diou Abstract Generalized additive models (GAMs) offer interpretability through independent univariate feature effects but underfit when interactions are present in data. GA2Ms add selected pairwise interactions which improves accuracy, but sacrifices interpretability and limits model auditing. We propose Conditionally Additive Local Models (CALMs), a new model class, that balances the interpretability of GAMs with the accuracy of GA2Ms. CALMs allow multiple univariate shape functions per feature, each active in different regions of the input space. These regions are defined independently for each feature as simple logical conditions (thresholds) on the features it interacts with. As a result, effects remain locally additive while varying across subregions to capture interactions. We further propose a principled distillation-based training pipeline that identifies homogeneous regions with limited interactions and fits interpretable shape functions via region-aware backfitting. Experiments on diverse classification and regression tasks show that CALMs consistently outperform GAMs and achieve accuracy comparable with GA2Ms. Overall, CALMs offer a compelling trade-off between predictive accuracy and interpretability. Machine Learning, ICML 1 Introduction In high-stakes decision making, machine learning models must provide not only reliable predictions but also human-understandable explanations (Murdoch et al., 2019; Doshi-Velez and Kim, 2017). This requirement has motivated the development of interpretable-by-design models, whose structure permits direct inspection of their behavior (Molnar, 2020). However, interpretability is a continuum rather than a binary property: interpretable-by-design models span a spectrum where gains in predictive accuracy often entail greater structural complexity (Rudin, 2019; Lipton, 2018). This trade-off is particularly evident in GAMs (Lou et al., 2012) and, their extension, GA2Ms (Caruana et al., 2015). GAMs achieve high interpretability by modeling predictions as a sum of independent univariate effects, but this strict additivity prevents them from capturing feature interactions, limiting their predictive accuracy. GA2Ms add selected pairwise interactions, improving performance at the cost of crucial interpretability properties (Radenovic et al., 2022). Specifically, interaction terms (i) obscure the unique attribution of the prediction to individual features because the pairwise interaction terms are generally not additively separable (i) complicate global auditing by requiring the simultaneous inspection of multiple (univariate and bivariate) effects. Section 3.3 analyzes these limitations in detail. ⋯·sxdx_dx1|x3x_1|x_3x2|x1x_2|x_1∑Σβ0 _0g−1g^-1y y Figure 1: CALMs use conditional feature effects: every feature effect is expressed by a collection of 1D functions, each associated with a different region of the input space. In the example, the effect of x1x_1 conditions on x3x_3, the effect of x2x_2 on x1x_1, while xdx_d does not interact with any other feature and thus has a single plot. CALMs are interpretable-by-design—summarized in d figures of 1D plots—and accurate, as they can model feature interactions. Code available at: https://github.com/givasile/CALM. We introduce Conditionally Additive Local Models (CALMs), a new model class that strikes a balance between the predictive accuracy of GA2Ms and the interpretability of GAMs. The core idea is to learn conditional feature effects: multiple univariate effects per feature, each active within a distinct region of the input space where the feature exhibits nearly additive behavior—that is, where it minimally interacts with others. These regions are defined through threshold-based conditions on interacting features. For example, if xix_i interacts with xjx_j, a CALM may learn separate effects for xix_i; one when xj<τx_j<τ and another when xj≥τx_j≥τ for some threshold τ. Figure 1 illustrates this representation. Conditional effects thus enable interaction-aware modeling while maintaining (i) transparent feature contributions and (i) straightforward global model auditing. To fit CALMs to data, we propose a three-step distillation-based pipeline: (1) Train a black-box reference model to capture complex interactions present in the data; (2) Partition the input space independently for each feature using CART-based splitters (Herbinger et al., 2022) optimized to minimize heterogeneity—a criterion that, when minimized, provably reduces feature interactions within the resulting regions (Herbinger et al., 2024); (3) Fit region-specific shape functions using a region-aware extension of the standard backfitting algorithm. This procedure efficiently identifies interacting features (e.g., xix_i interacts with xjx_j), optimal split points (e.g., τ such that separate effects are learned for xj<τx_j<τ and xj≥τx_j≥τ), and requires minimal tuning—primarily the maximum allowed number of interactions per feature. Extensive evaluation on diverse regression and classification benchmarks demonstrates that CALMs consistently outperform GAMs in predictive accuracy, and often match or exceed GA2Ms. Furthermore, we formally show that CALMs satisfy key interpretability properties unsupported by GA2Ms. Code for reproducing all experiments is provided in the supplementary material. Contributions. (i) We introduce Conditionally Additive Local Models (CALMs), a novel interpretable-by-design model class that balances GAM’s interpretability and GA2Ms accuracy via conditional feature effects. (i) We develop a robust training algorithm based on distillation, heterogeneity-minimization and region-aware backfitting. (i) We provide a formal analysis of the interpretability properties inherent in CALMs. (iv) We present extensive empirical evidence demonstrating improved accuracy–interpretability trade-offs over GAMs and GA2Ms. 2 Background and Related Work Interpretable-by-design models enforce transparency in their decision-making process. Examples include decision trees (Breiman et al., 1984), rule lists (Letham et al., 2015; Angelino et al., 2018), and prototype-based classifiers (Chen et al., 2019; Bien and Tibshirani, 2011). Among them, GAMs (Hastie and Tibshirani, 1990; Wood, 2017) stand as a particularly popular candidate, mainly due to their simple global interpretability. GAMs model the output as g([y])=β0+∑ifi(xi)g(E[y])= _0+ _if_i(x_i), where each fif_i is a univariate shape function, β0 _0 is the global intercept, and g is a link function. The key benefit of GAMs is that feature effects can be independently visualized through a simple 1D plot. Methods for learning fi∀if_i\>∀ i include spline-based approaches (Lou et al., 2012; Wood, 2017), gradient boosting (Friedman, 2001), and neural-based variants (Agarwal et al., 2021; Kraus et al., 2024; Radenovic et al., 2022). However, standard GAMs assume that features contribute independently to the output, which limits their accuracy when feature interactions are present in data. To capture these interactions, GA2Ms add pairwise terms: g([y])=β0+∑ifi(xi)+∑i<jfij(xi,xj)g(E[y])= _0+ _if_i(x_i)+ _i<jf_ij(x_i,x_j) (Lou et al., 2013). GA2Ms instances vary in (i) how they select the top-K interactions to maintain sparsity and (i) how they learn fif_i and fijf_ij. For example, (Lou et al., 2013) greedily selects interactions based on loss reduction, while, neural-based approaches (Yang et al., 2020; Ibrahim et al., 2023; Chang et al., 2022) integrate pairwise interactions into standard neural-based architectures. While GA2Ms improve accuracy, they obscure individual feature contributions and complicate model auditing (see Section 3.3). A key open question remains: Can we achieve GA2M-level accuracy while maintaining GAM-level interpretability? To this end, we draw inspiration from regional effect methods, which handle interactions by partitioning the feature space into subregions where interactions are weak. To understand regional effects, consider first how standard global effect work. Global effect plots, such as PDP, ALE or SHAP-DP, explain a black box model by decomposing its complex d-dimensional function f()→yf(x)→ y into d one-dimensional plots xi→yx_i→ y. However, when strong interactions exist between xix_i and other features, these plots can yield misleading explanations (Gkolemis et al., 2023a, b). Regional feature effect plots address that by partitioning the feature space into subregions where interactions are minimal (Herbinger et al., 2022, 2024). Within each subregion, simple univariate plots accurately represent feature effects without being confounded by interactions with other features. We adopt this strategy to identify low-interaction subregions and then we fit shape functions within each subregion. Our approach relates to model distillation (Hinton et al., 2015; Tan et al., 2018) and surrogate modeling (Guidotti et al., 2018). These approaches train an interpretable “student” model to mimic a complex “teacher”, either locally (Ribeiro et al., 2016) or globally (Guidotti et al., 2018). However, instead of explaining the teacher, we use it to identify subregions with minimal feature interactions, on which we then fit shape functions. In a rough analogy, CALM serves as the interpretable-by-design “student” that replaces the black-box “teacher” as the final predictor. Other works that share the above idea are: SLIM (Hu et al., 2020) and related analyses (Herbinger et al., 2023) propose tree-based models with simple predictors in each region, mainly in the context of model distillation and explanation. In contrast, our goal is not to approximate a black-box model, but to construct an interpretable-by-design predictor with GAM-like structure. Related ideas also appear in (Gkolemis et al., 2023c), which leverages regional feature effects to build interpretable models; however, that work does not explicitly characterize the interpretability guarantees of the learned model or evaluate them systematically. 3 CALM: Conditionally Additive Local Model CALM captures feature interactions via conditional feature effects, a set of univariate shape functions per feature, each active in a different region of the input space. 3.1 Model Formulation Let =(X1,…,Xd)∈⊆ℝdX=(X_1,…,X_d) ^d be a random vector with joint distribution P_X over the input data, =(x1,…,xd)x=(x_1,…,x_d) where =1×⋯×dX=X_1×·s×X_d. We use the subscript −i-i to denote quantities excluding the i-th feature (e.g., −iX_-i defined on space −iX_-i). Let also Y be the output random variable, with Y∈ℝY for regression and Y∈0,1Y∈\0,1\ for binary classification. CALM is then defined as g([Y∣=])=β0+∑i=1dfi(ri(−i))(xi)g (E[Y =x] )= _0+ _i=1^df_i^(r_i(x_-i))(x_i) (1) where β0∈ℝ _0 is an intercept and g is a link function (identity for regression and logit for classification). For each feature i, CALM learns a set of univariate shape functions fi(r)r=1Ri\f_i^(r)\_r=1^R_i. At prediction time, a region selection function ri(−i):ℝd−1→1,…,Rir_i(x_-i):R^d-1→\1,…,R_i\ chooses the specific shape function for xix_i based on the context of other features −ix_-i. This allows the effect of xix_i to vary across regions, effectively modeling feature interactions, while maintaining univariate interpretability within each region. The final prediction is given by: fCALM()=g−1(β0+∑i=1dfi(ri(−i))(xi))f_ CALM(x)=g^-1 ( _0+ _i=1^df_i^(r_i(x_-i))(x_i) ) (2) For the rest of the paper, we may use y^() y(x) to denote fCALMf_ CALM for brevity. The selection function ri(−i)r_i(x_-i) partitions the input space (excluding xix_i) into i:=ℛi(r)r=1RiP_i:=\R_i^(r)\_r=1^R_i, independently for each feature i, where each ℛi(r)⊆−iR_i^(r) _-i. The goal of these partitions is to minimize interactions between xix_i and other features within each region, ensuring the effect of xix_i is accurately captured by a single univariate shape function. To maintain interpretability, each partition is represented by a binary decision tree TiT_i of maximum depth dmaxd_max, built over −ix_-i. Internal nodes apply axis-aligned splits, i.e., inequality thresholds xj<τx_j<τ or xj≥τx_j≥τ for continuous features, and equality tests xj=τx_j=τ or xj≠1x_j≠ 1 for categorical ones. Each leaf of the tree defines a region ℛi(r)R_i^(r). Formally, each region is defined by a conjnuction of at most mi(r)≤dmaxm_i^(r)≤ d_max rules: ℛi(r)=−i|∧k=1mi(r)(xjkopkτk).R_i^(r)= \x_-i\; |\; _k=1^m_i^(r)(x_j_k\;op_k\; _k) \. (3) where jk∈1,…,d∖ij_k∈\1,…,d\ \i\ indexes an interacting feature, opk∈<,≥,=,≠op_k∈\<,≥,=\,,≠\ is a comparison operator, and τk∈ℝ _k is a threshold. CALM is able to capture feature interactions involving up to (dmax+1)(d_max+1) features, as each region conditions on up to dmaxd_max features. While increasing dmaxd_max can improve accuracy, it comes at the cost of interpretability. To balance this accuracy-interpretability tradeoff, CALM exposes two hyperparameters: dmaxd_max (tree depth) bounds the number of shape functions per feature to Ri≤2dmaxR_i≤ 2^d_max, while K limits the total number of shape functions across all features via ∑iRi≤K _iR_i≤ K. In our experiments, we set dmax=2d_max=2 (yielding Ri=4R_i=4 regions per feature) and leave K unconstrained. Notably, CALM learns, at most, d⋅2dmaxd· 2^d_max univariate shape functions (up to 2dmax2^d_max per feature), yet these combine to express up to 2d⋅dmax2^d· d_max distinct additive models–an exponential increase in expressive power. 3.2 Training algorithm for fitting CALM Given a dataset =((i),y(i))i=1ND=\(x^(i),y^(i))\_i=1^N, we fit a CALM predictor by minimizing the empirical risk (,Y)∼P^[ℒ(Y,y^())]E_(X,Y) P_D [L(Y, y(X)) ], where y^() y(X) denotes a CALM model, ℒL a loss function and P P_D the empirical distribution over D. Fitting a CALM requires estimating: (a) a set of d partitions iP_i, i∈1,…,di∈\1,…,d\, represented by binary trees TiT_i (one per feature); (b) region-specific shape functions fi(r)r=1Ri\f_i^(r)\_r=1^R_i and a global intercept β0 _0. We propose a three-step distillation-based pipeline summarized in Algorithm 1. Algorithm 1 Training a CALM model 0: Training data =((i),y(i))i=1ND=\(x^(i),y^(i))\_i=1^N 0: CALM predictor fCALM()f_ CALM(x) 1: Step 1: Train reference model freff_ ref on D. 2: Step 2: Learn feature-specific trees TiT_i, i=1,…,di=1,…,d. 3: Step 3: Estimate shape functions fi(r)r=1Rii=1d\\f_i^(r)\_r=1^R_i\_i=1^d. Step 1: Train a Reference Model. We train a high-capacity black-box predictor fref:→ℝf_ ref:X on D. This model serves as functional proxy for [Y|]E[Y|X] and is used exclusively to detect interactions. The reference model is discarded after training. The choice of freff_ ref is independent of later stages; any accurate predictor can be used, such as gradient-boosted trees (our default), neural networks, random forests or foundation models, like TabPFN (Hollmann et al., 2022). Step 2: Learn feature-specific partitioning trees. This step learns d independent, feature-specific partitions i:=ℛi(r)r=1RiP_i:=\R_i^(r)\_r=1^R_i, each represented by a binary tree TiT_i defined over −ix_-i. The objective is to identify near-additive regions ℛi(r)R_i^(r), in which the effect of xix_i on the output is well approximated by a univariate function fi(r)f_i^(r), due to locally weak higher-order interactions. To identify such regions, we use the interaction-related heterogeneity measure from (Herbinger et al., 2024). Feature effect methods (e.g., PDP or ALE) explain a black-box model (like fref()f_ ref(x)) by decomposing it into univariate effects fi(xi)∀if_i(x_i)\>∀ i. To do so, they first define the local effects h(xi,−i(j))h(x_i,x_-i^(j)) that quantify the contribution of feature xix_i on a specific instance (j)x^(j). Then they compute fi(xi)f_i(x_i) by averaging the local effects. Heterogeneity is the variability of these local effects around their average, directly measuring how much xix_i’s effect depends on other features. High heterogeneity signals strong interactions, while low heterogeneity identifies near-additive, interaction-free regions—making it an effective criterion for our purpose. Following the above, we define the pointwise heterogeneity of xix_i in a region ℛR (e.g., ℛi(r)R_i^(r)) as: Hiℛ(xi) H_i^R(x_i) =−i∣−i∈ℛ[(h(xi,−i)−μiℛ(xi))2] =E_X_-i _-i [ (h(x_i,X_-i)- _i^R(x_i) )^2 ] (4) where μiℛ(xi)=−i∣−i∈ℛ[h(xi,−i)] _i^R(x_i)=E_X_-i _-i [h(x_i,X_-i) ] The corresponding feature-level heterogeneity is: Hiℛ=Xi|−i∈ℛ[Hiℛ(Xi)]H^R_i=E_X_i|X_-i [H_i^R(X_i) ] (5) In our experiments, heterogeneity is computed using PDP-based formulas. As shown in (Herbinger et al., 2024), ALE- and SHAP-DP variants also provide valid alternatives; we refer to (Herbinger et al., 2022, 2024) for details. As additional theoretical justification for using heterogeneity as the splitting criterion, we show that, under specific assumptions, the approximation error of CALM is bounded by the sum of the expected feature heterogeneities. Consequently, reducing heterogeneity of each feature through partitioning yields a CALM model fCALMf_ CALM with lower mean squared error, assuming perfect fit in Step 3 of Algorithm 1. The proof is provided in Appendix A (Proposition A.1). Therefore, we use HiℛH_i^R as the splitting criterion and learn a binary decision tree TiT_i for each feature xix_i using a greedy CART-style algorithm. At each node, we consider splits over all conditioning features xj≠xix_j≠ x_i and select the split that maximizes the reduction in heterogeneity. For numerical features, splits take the form xj≤,>τx_j\≤,>\τ, while for categorical features we use xj=,≠τx_j\=,≠\τ. A split is accepted only if the relative reduction in HiℛH_i^R exceeds a threshold ϵε; otherwise, the node becomes a leaf. This procedure yields a partition i=ℛi(r)r=1RiP_i=\R_i^(r)\_r=1^R_i where heterogeneity is significantly reduced, ensuring the effect of xix_i is locally stable and less distorted by interactions than in the global space. Additional implementation details are provided in Appendix A. In all experiments, we use PDP-based heterogeneity, set the maximum tree depth to dmax=2d_ =2, and fix ϵ=0.2ε=0.2. Step 3: Estimate Shape Functions. Given the partitions ii=1d\P_i\_i=1^d, we estimate the region-specific shape functions fi(r)\f_i^(r)\ by minimizing the empirical loss of the CALM predictor. We use a modified gradient boosting procedure in which, at each iteration, a single shape function fi(r)f_i^(r) is updated using only the observations whose −ix_-i fall into the corresponding region ℛi(r)R_i^(r) (see Appendix A for details). Although each update is restricted to a single region, the resulting optimization problem is inherently coupled. The regions are defined separately for each feature, so the subsets of samples used to update different shape functions generally overlap. Consequently, updating one fi(r)f_i^(r) changes the residuals seen by all other shape functions. Therefore it can be understood as a coordinated optimization of a single additive predictor with region-gated components. The following proposition formalizes this intuition by characterizing the target of Step 3 at the population level and its convergence behavior under idealized updates. The proof is in Appendix A. In practice, Step 3 is implemented using gradient boosting to approximate the exact regional updates. Proposition 3.1 (Optimality and convergence). Let m()=[Y∣=]m(x)=E[Y =x] denote the true regression function. Assume regression with squared loss, fixed partition trees Tii=1d\T_i\_i=1^d (as learned by Step 2) and [Y2]<∞E[Y^2]<∞. Let ℋ(Ti)H(\T_i\) denote the fixed-tree CALM class (Appendix A.3.1), and assume ℋ(Ti)⊂L2(PX)H(\T_i\)⊂ L_2(P_X) is nonempty, closed, and convex. Then: 1. Optimality: Any s⋆∈argmins∈ℋ(Ti)[(Y−s())2]s ∈ _s (\T_i\)E[(Y-s(X))^2] is the L2(P)L_2(P_X)-best approximation of m within ℋ(Ti)H(\T_i\). 2. Convergence: The idealized exact cyclic regional backfitting converges to an empirical risk minimizer over ℋ(Ti)H(\T_i\). 3.3 Interpretability of a CALM CALMs offer interpretability similar to GAMs, as both require inspecting d univariate plots; one per feature. CALMs are slightly more complex: unlike GAMs, each plot can contain up to 2dmax2^d_max curves corresponding to a different region of the input space and with interaction-induced discontinuities marked by vertical lines. However, because regions are interpretable (specified by simple threshold conditions), these region-specific curves provide explanations which are suitable for model auditing and decision support. For example: The effect of age on _ mortality\_rate follows this curve when the patient is male and has a BMI above 30. Figure 2: CALM plot for x1x_1. Each curve shows the contribution of x1x_1 to y (P1) in a different region of the input space (x3≶0x_3 0). Vertical lines indicate interaction-induced discontinuities, e.g., at x1≈−0.4x_1≈-0.4 there exist a positive hidden jump of [0,0.37][0,0.37] due to interaction of x1x_1 with x2x_2, which must be considered when assessing regional sensitivity (P2) or global properties (P3). Interpreting a CALM plot. In Figure 2 each curve gives the contribution of x1x_1 to y in a specific region; the blue curve when x3>0x_3>0 and the orange curve when x3≤0x_3≤ 0. For example, at x1=−0.5x_1=-0.5, the contribution is approximately −0.2-0.2 (blue) or −0.75-0.75 (orange), depending on x3x_3. The plots also illustrate how altering x1x_1 to x1→x1+Δx_1→ x_1+ x impacts the prediction. Vertical dotted lines mark points of a hidden discontinuity which is due to x1x_1 participating as an interaction term for feature x2x_2. As shown in Figure 1, the effect of x2x_2 is conditioned by x1≤−0.4x_1≤-0.4, −0.4≤x1≤0.4-0.4≤ x_1≤ 0.4 and x1>0.4x_1>0.4, therefore in Figure 2 we observe vertical lines in x1±0.4x_1± 0.4. If a change in x1x_1 does not cross a vertical line, the change in the output (Δy)( y) equals the curve difference (Δfi)( f_i). Crossing a line signifies a hidden jump, in the range [α,β][α,β], so Δfi+α≤Δy≤Δfi+β f_i+α≤ y≤ f_i+β. Arrows provide a fast understanding of the jump: ↑ means Δy>Δfi y> f_i, ↓ means Δy<Δfi y< f_i, ↕ means it depends. Below, we outline three crucial interpretability properties, along with discussion about their satisfiability by GAM , GA2M, and CALM. Complete proofs are provided in Appendix B. P1. Local Feature Contribution: What is the contribution of each feature to the prediction? Formally: Given an input x, how much does each xix_i contribute to y^() y(x)? The explanation is local—it concerns a specific input x. In GAM , the contribution is fi(xi)f_i(x_i). In GA2M, the contribution is neither explicit nor unique; it must be inferred post-hoc (e.g., via SHAP or LIME), with each method relying on different assumptions and yielding different results. In CALM, the contribution is fi(r(−i))(xi)f_i^(r(x_-i))(x_i). P2. Regional Feature Sensitivity: How does changing xix_i change the prediction? Formally: Given xix_i and Δx>0 x>0, what is Δy^=y^(+iΔx)−y^() y= y(x+e_i x)- y(x), assuming only xix_i is perturbed and −ix_-i is fixed? The explanation is regional—it characterizes the effect of xix_i on a specific region ([xi,xi+Δx][x_i,x_i+ x]) independently of the values of other features. In GAM , the change is Δy^=fi(xi+Δx)−fi(xi) y=f_i(x_i+ x)-f_i(x_i). In GA2M, the change cannot be determined without knowing the values of all features that interact with xix_i. In CALM, an exact answer is possible only when the perturbation does not cross a vertical line. In this case, the resulting change is a set of values Δy:=Δfi(r)r=1Ri y:=\ f_i^(r)\_r=1^R_i, one for each curve in the plot: Δf(r)=fi(r)(xi+Δx)−fi(r)(xi) f^(r)=f_i^(r)(x_i+ x)-f_i^(r)(x_i). If a crossing occurs, an exact answer is not attainable, however, for less precise questions, such as whether the Δx x change in xix_i will have a positive impact on y, an answer is still feasible (see P3. below). P3. Global Feature Property: Is the model globally monotonic increasing with respect to xix_i? Formally: For all x and Δx>0 x>0, is Δy^=y^(+iΔx)−y^()>0 y= y(x+e_i x)- y(x)>0? The explanation is global—it assesses the monotonicity of xix_i across the entire input space. In GAM , monotonicity is easily determined by the shape of fi(xi)f_i(x_i). In GA2M, verifying monotonicity requires a concurrent examination of the 1D shape of fi(xi)f_i(x_i) along with all pairwise interactions fijf_ij along all j, an inspection which is infeasible for a human. In CALM, if no vertical lines are present, simply check whether all curves fi(r)(xi)f_i^(r)(x_i) for r=1,…,Rir=1,…,R_i are monotonically increasing. If vertical lines exist, simply verify that all arrows are positive. 3.4 Efficiency The efficiency of fitting a CALM to a dataset is mainly affected by the number of instances N and the feature dimensionality d. Step 1 involves fitting a black box model, which, in general, is computationally efficient. Our default setup uses an XGBoost model that fits in a few seconds. Step 2 is the computational bottleneck as it requires fitting d binary trees TiT_i, resulting in a worst-case cost of (dmaxd2NlogN)O(d_maxd^2N N). Since the trees are shallow (dmaxd_max is small), this reduces to (d2NlogN)O(d^2N N), which is acceptable in practice given that d is typically on the order of tens for tabular datasets. Crucially, the local effects needed for evaluating candidate splits are computed once for the entire dataset and stored. Then, at each candidate split, the algorithm checks the precomputed effects that relate to the samples in the node (indexed lookups), avoiding repeated evaluations of the black-box model. While this does not change the asymptotic complexity, it makes each split evaluation computationally inexpensive and is critical for practical efficiency; in our experiments, this entire step completes in a few seconds. Step 3 simply requires fitting the regional shape functions, which depends on the underlying fitting approach. Our default setup uses gradient boosting which is computationally efficient. Overall our method is efficient, typically fitting most tabular datasets within a few seconds, as shown by the runtimes reported in the Appendix C. 4 Empirical Evaluation (a) GAM (b) CALM (c) GA2M Figure 3: Explanatory plots for the three synthetic regression tasks: top row (Task 1) shows two-region interaction; middle row (Task 2) shows four-region interaction; bottom row (Task 3) shows general interactions. The empirical evaluation of CALM includes three synthetic regression datasets and 25 public real-world tabular datasets—10 for classification and 15 for regression. Across all experiments, CALM is applied with its default configuration: Step 1 uses XGBoost as the black-box model. Step 2 uses PDP-based heterogeneity with dmax=2d_max=2, ϵ=0.2ε=0.2, and K remains unconstrained. Step 3 uses standard gradient boosting. All predictive results are reported as mean ± standard deviation over standard 5-fold cross-validation. Full experimental details are in Appendix C. Implementation details. We used well-known open-source implementations for all baseline models: scikit-learn (random forests), XGBoost, tensorflow (neural networks and NAM), interpretml (EBM, EB2M), pygam (SPLINE), and official repositories for NODE-GA2M(Chang et al., 2022) and GAMI‐Net(Yang et al., 2020). 4.1 Synthetic example We compare CALM against EBM (as the GAM baseline) and EB2M (as the GA2M baseline) on three synthetic datasets. Both EBM and EB2M are used with their default parameters. Evaluation is based on R2R^2 performance. Each dataset contains 1000 samples with features drawn from xi∼(−1,1)x_i (-1,1). Table 1 summarizes the accuracy of each approach. Table 1: Synthetic examples: R2R^2 comparison (mean ± std. dev. over 5-fold CV). Dataset GAM GA2M CALM Case 1 0.737±0.0280.737 $±$0.028 0.974±0.0110.974 $±$0.011 0.995±0.0020.995 $±$0.002 Case 2 0.479±0.0520.479 $±$0.052 0.712±0.0360.712 $±$0.036 0.949±0.0240.949 $±$0.024 Case 3 0.527±0.0230.527 $±$0.023 0.961±0.0090.961 $±$0.009 0.975±0.0060.975 $±$0.006 Table 2: Classification accuracy vs. baselines (mean ± std. dev. over 5-fold CV). ✓/✗ denotes higher/lower accuracy than the average GAM (or GA2M) baseline. Dataset B GAM CALM GA2M CALM vs. XGB NAM EBM EB2M NODE GAMI GAM GA2M Adult 0.870±0.0020.870 $±$0.002 0.851±0.0020.851 $±$0.002 0.870±0.0020.870 $±$0.002 0.870±0.0020.870 $±$0.002 0.871±0.0030.871 $±$0.003 0.850±0.0020.850 $±$0.002 0.857±0.0030.857 $±$0.003 ✓ ✓ COMPAS 0.661±0.0070.661 $±$0.007 0.682±0.0160.682 $±$0.016 0.681±0.0120.681 $±$0.012 0.684±0.0130.684 $±$0.013 0.684±0.0110.684 $±$0.011 0.667±0.0130.667 $±$0.013 0.686±0.0130.686 $±$0.013 ✓ ✓ HELOC 0.717±0.0130.717 $±$0.013 0.723±0.0110.723 $±$0.011 0.728±0.0140.728 $±$0.014 0.728±0.0120.728 $±$0.012 0.730±0.0120.730 $±$0.012 0.715±0.0080.715 $±$0.008 0.725±0.0120.725 $±$0.012 ✓ ✓ MIMIC2 0.890±0.0010.890 $±$0.001 0.886±0.0010.886 $±$0.001 0.886±0.0030.886 $±$0.003 0.886±0.0030.886 $±$0.003 0.886±0.0030.886 $±$0.003 0.888±0.0020.888 $±$0.002 0.884±0.0030.884 $±$0.003 ✓ ✓ Appendicitis 0.868±0.0690.868 $±$0.069 0.848±0.0560.848 $±$0.056 0.877±0.0770.877 $±$0.077 0.878±0.0630.878 $±$0.063 0.869±0.0790.869 $±$0.079 0.849±0.0160.849 $±$0.016 0.764±0.0740.764 $±$0.074 ✓ ✓ Phoneme 0.898±0.0060.898 $±$0.006 0.808±0.0020.808 $±$0.002 0.821±0.0070.821 $±$0.007 0.861±0.0110.861 $±$0.011 0.863±0.0070.863 $±$0.007 0.869±0.0030.869 $±$0.003 0.876±0.0090.876 $±$0.009 ✓ ✗ SPECTF 0.862±0.0290.862 $±$0.029 0.839±0.0300.839 $±$0.030 0.894±0.0150.894 $±$0.015 0.894±0.0150.894 $±$0.015 0.882±0.0170.882 $±$0.017 0.820±0.0670.820 $±$0.067 0.874±0.0160.874 $±$0.016 ✓ ✓ Magic 0.885±0.0040.885 $±$0.004 0.850±0.0060.850 $±$0.006 0.857±0.0050.857 $±$0.005 0.864±0.0040.864 $±$0.004 0.872±0.0030.872 $±$0.003 0.876±0.0040.876 $±$0.004 0.874±0.0030.874 $±$0.003 ✓ ✗ Bank 0.908±0.0030.908 $±$0.003 0.901±0.0030.901 $±$0.003 0.902±0.0020.902 $±$0.002 0.905±0.0020.905 $±$0.002 0.909±0.0020.909 $±$0.002 0.903±0.0010.903 $±$0.001 0.908±0.0030.908 $±$0.003 ✓ ✗ Churn 0.958±0.0040.958 $±$0.004 0.885±0.0090.885 $±$0.009 0.886±0.0050.886 $±$0.005 0.946±0.0030.946 $±$0.003 0.957±0.0090.957 $±$0.009 0.954±0.0090.954 $±$0.009 0.952±0.0070.952 $±$0.007 ✓ ✗ Table 3: Regression RMSE vs. baselines (mean ± std. dev. over 5-fold CV). ✓/✗ denotes lower/higher RMSE than the average GAM (or GA2Ms) baseline. Dataset B GAM CALM GA2M CALM vs. XGB NAM EBM EB2M NODE GAMI GAM GA2M Bike Sharing 39.35±1.3839.35 $±$1.38 101.97±1.31101.97 $±$1.31 100.21±1.01100.21 $±$1.01 55.67±1.2255.67 $±$1.22 54.80±0.9654.80 $±$0.96 54.47±1.5854.47 $±$1.58 53.44±1.9053.44 $±$1.90 ✓ ✗ California Housing 0.45±0.010.45 $±$0.01 0.61±0.010.61 $±$0.01 0.55±0.010.55 $±$0.01 0.51±0.010.51 $±$0.01 0.49±0.010.49 $±$0.01 0.50±0.010.50 $±$0.01 0.51±0.040.51 $±$0.04 ✓ ✗ Parkinsons Motor 1.44±0.091.44 $±$0.09 6.11±0.166.11 $±$0.16 4.20±0.094.20 $±$0.09 2.24±0.132.24 $±$0.13 2.35±0.062.35 $±$0.06 3.49±0.313.49 $±$0.31 2.76±0.312.76 $±$0.31 ✓ ✓ Parkinsons Total 1.86±0.081.86 $±$0.08 7.90±0.107.90 $±$0.10 4.85±0.114.85 $±$0.11 2.97±0.092.97 $±$0.09 2.77±0.062.77 $±$0.06 4.60±0.534.60 $±$0.53 3.81±0.663.81 $±$0.66 ✓ ✓ Seoul Bike 209.6±3.47209.6 $±$3.47 320.2±4.38320.2 $±$4.38 303.7±3.86303.7 $±$3.86 238.9±1.64238.9 $±$1.64 235.2±1.58235.2 $±$1.58 231.0±4.23231.0 $±$4.23 245.82±8.45245.82 $±$8.45 ✓ ✗ Wine 0.62±0.010.62 $±$0.01 0.72±0.010.72 $±$0.01 0.70±0.010.70 $±$0.01 0.69±0.020.69 $±$0.02 0.68±0.010.68 $±$0.01 0.69±0.010.69 $±$0.01 0.71±0.000.71 $±$0.00 ✓ ✗ Energy 67.97±2.5867.97 $±$2.58 89.80±2.6789.80 $±$2.67 85.23±2.4985.23 $±$2.49 83.09±1.9783.09 $±$1.97 77.33±2.4677.33 $±$2.46 82.18±2.5382.18 $±$2.53 180.1±84.09180.1 $±$84.09 ✓ ✓ CCPP 3.09±0.093.09 $±$0.09 4.21±0.054.21 $±$0.05 3.44±0.083.44 $±$0.08 3.42±0.073.42 $±$0.07 3.28±0.073.28 $±$0.07 3.97±0.073.97 $±$0.07 3.92±0.073.92 $±$0.07 ✓ ✓ Electrical 0.04±0.020.04 $±$0.02 0.04±0.010.04 $±$0.01 0.02±0.010.02 $±$0.01 0.02±0.010.02 $±$0.01 0.02±0.010.02 $±$0.01 0.01±0.010.01 $±$0.01 0.02±0.010.02 $±$0.01 ✓ ✗ Elevators 0.00±0.000.00 $±$0.00 0.00±0.000.00 $±$0.00 0.00±0.000.00 $±$0.00 0.00±0.000.00 $±$0.00 0.00±0.000.00 $±$0.00 0.00±0.000.00 $±$0.00 0.00±0.000.00 $±$0.00 ✓ ✓ No2 0.47±0.030.47 $±$0.03 0.50±0.020.50 $±$0.02 0.49±0.030.49 $±$0.03 0.49±0.040.49 $±$0.04 0.47±0.030.47 $±$0.03 0.52±0.020.52 $±$0.02 0.50±0.030.50 $±$0.03 ✓ ✓ Sensory 0.51±0.030.51 $±$0.03 0.48±0.010.48 $±$0.01 0.48±0.010.48 $±$0.01 0.45±0.020.45 $±$0.02 0.45±0.020.45 $±$0.02 0.49±0.010.49 $±$0.01 0.46±0.030.46 $±$0.03 ✓ ✓ Airfoil 1.54±0.101.54 $±$0.10 4.80±0.204.80 $±$0.20 4.57±0.154.57 $±$0.15 2.50±0.152.50 $±$0.15 2.17±0.092.17 $±$0.09 2.13±0.112.13 $±$0.11 4.79±0.214.79 $±$0.21 ✓ ✓ Skill Craft 0.94±0.020.94 $±$0.02 0.93±0.030.93 $±$0.03 0.90±0.020.90 $±$0.02 0.90±0.020.90 $±$0.02 0.90±0.030.90 $±$0.03 0.91±0.020.91 $±$0.02 1.38±0.801.38 $±$0.80 ✓ ✓ Ailerons 0.00±0.000.00 $±$0.00 0.00±0.000.00 $±$0.00 0.00±0.000.00 $±$0.00 0.00±0.000.00 $±$0.00 0.00±0.000.00 $±$0.00 0.00±0.000.00 $±$0.00 0.00±0.000.00 $±$0.00 ✓ ✓ Case 1. We set y=x12+log(|x2|)+2sin(π2x3)x2≥0+2cos(π2x3)x2<0y=x_1^2+ (|x_2|)+2 ( π2x_3 )1_x_2≥ 0+2 ( π2x_3 )1_x_2<0. Since GAM cannot model the x2−x3x_2-x_3 interaction, it simplifies the x3→yx_3→ y effect by averaging the sine and cosine modes into one curve (Fig. 3(a), top row), resulting in R2=0.737R^2=0.737. Both GA2M and CALM capture the interaction correctly, achieving near-perfect accuracy (0.9740.974 and 0.9950.995). However, their interpretability differs. GA2M requires three views to interpret x3x_3’s effect: (i) the main effect averaging the sine and cosine cases, (i) a near-zero x1x_1–x3x_3 heatmap, and (i) an x2x_2–x3x_3 heatmap capturing the sign-dependent interaction (Fig. 3(c), top row). In contrast, CALM simplifies the interpretation with two 1D plots and a clear separation of the regimes by conditioning on x2x_2’s sign (Fig. 3(b), top row). Case 2. We set: y=x02+log(|x1|)+2[sin(π2x2)x0≥0,x1≥0+cos(π2x2)x0≥0,x1<0+sin(2πx2)x0<0,x1≥0+cos(2πx2)x0<0,x1<0]y=x_0^2+ (|x_1|)+2 [ ( π2x_2 )1_x_0≥ 0,\,x_1≥ 0+ ( π2x_2 )1_x_0≥ 0,\,x_1<0+ (2π x_2)1_x_0<0,\,x_1≥ 0+ (2π x_2)1_x_0<0,\,x_1<0 ]. GAM (Fig. 3(a), middle row) merges all four modes into a single curve, limiting its accuracy to R2=0.479R^2=0.479. GA2M (Fig. 3(c), middle row) captures partial structure through two 2D interactions (x1x_1–x3x_3 and x2x_2–x3x_3) but misses the full 3-way interaction, reaching R2=0.712R^2=0.712. Furthermore, the model’s behavior is hard to grasp, as it requires integrating information from three distinct plots: the main effect and two interaction heatmaps. CALM (Fig. 3(b), middle row) separates the four regimes into distinct 1D curves, achieving both high accuracy (R2=0.949R^2=0.949) and clear interpretability. Case 3. We set y=x12+log(|x2|)sin(π2x3)y=x_1^2+ (|x_2|) ( π2x_3), where the logarithmic effect of x2x_2, log|x2| |x_2|, is modulated by the sinusoidal term sin(π2x3) ( π2x_3), resulting in a general interaction structure. GAM (Fig. 3(a), bottom row) fails to capture the interaction and instead fits a blurred sinusoidal curve, leading to low accuracy (R2=0.527R^2=0.527). GA2M (Fig. 3(c), bottom row) achieves high accuracy (R2=0.961R^2=0.961) by capturing the interaction in the x2x_2–x3x_3 heatmap. Still, interpretation remains difficult, as understanding the x2→yx_2→ y effect requires integrating three views: the main effect plot and two interaction heatmaps. CALM (Fig. 3(b), bottom row) cannot fully express the continuous x2x_2–x3x_3 interaction, but approximates it by partitioning x2x_2 and assigning each segment an average logarithmic response. While simplified, this yields a high accuracy (R2=0.975R^2=0.975) and a concise, transparent explanation that remains close to the underlying behavior. 4.2 Evaluation on Real Datasets On each real dataset, we evaluate CALM against three black-box baselines (neural networks, random forests, and XGBoost), three GAM models (NAM, EBM, and SPLINE), and three GA2M models: EB2M (EBM with pairwise interactions enabled), NODE-GA2M, and GAMI‐Net. Evaluation is based on accuracy for classification tasks and RMSE for regression datasets. Classification Results. On classification datasets (Table 2), CALM: (i) outperforms the average performance of GAM s in 9/10 datasets and matches it on the remaining one; (i) outperforms the average performance of GA2Ms on 5/10 datasets and matches it on one case; (i) outperforms the black-box model on 4/10 datasets of the datasets and matches it on one. Regression Results. On regression datasets (Table 3), CALM: (i) outperforms the average performance of GAM s in 14/15 datasets and matches it on the remaining one; (i) outperforms the average performance of GA2Ms on 8/15 datasets and matches it on two; (i) outperforms the black-box model on 3/10 datasets of the datasets and matches it on two. Model Complexity. In addition to its strong predictive performance, CALM achieves its results using remarkably few pairwise feature interactions, 6.2 on average for classification and 15.1 for regression. On classification tasks, this is fewer than EB2M (15.2), GAMI-Net (17.6), and NODE-GA2M (97.7). On regression tasks, CALM remains sparser than GAMI-Net (16.3) and NODE-GA2M (87.3), and close to EB2M (14.1), as shown in Tables 6 and 7 of Appendix C. While most GA2M methods allow explicit control over the number of interactions, we apply them using their default configurations. In contrast, CALM does not impose an interaction cap, but still discovers a sparse structure through region-based conditioning. Runtime. On classification tasks, CALM is slightly slower than EBM and EB2M, but noticeably faster than the other two GA2M models (GAMI-Net and NODE-GA2M), and even faster than NAM, despite the fact that NAM does not model interactions. In regression tasks, CALM exhibits a runtime similar to NAM and EB2M, and remains significantly faster than the remaining GA2M baselines. Average runtimes across datasets are reported in Tables 8 and 9 (Appendix C). Summary. Across both classification and regression tasks, CALM consistently outperforms the average GAM models and achieves comparable performance to GAM models in most cases. Remarkably, it surpasses the best individual GA2M on selected datasets. These gains are achieved with fewer interactions (highlighting the model’s inherent sparsity) and with practical runtimes, significantly lower than those of complex GA2M architectures. Together, these results demonstrate that CALM strikes an effective balance between accuracy, simplicity, and efficiency. 5 Conclusions We introduced CALM, a novel class of interpretable-by-design models that bridge the gap between the interpretability of GAMs and the accuracy of GA2Ms. CALMs model feature interactions using conditional feature effects, a set of univariate shape functions per feature, conditioned on its interacting features. Extensive evaluation shows that conditional effects are sufficient to preserve, or even exceed, the accuracy of GA2Ms (which only model bivariate interactions), without resorting to 3D plots or heatmaps. CALM main limitations is that interpretability may harden as the number of interactions increases. Although each feature is linked to a maximum of 2dmax=42^d_max=4 shape functions, extensive cross-feature dependencies can clutter plots with vertical lines, making them harder to read. Furthermore, the method’s accuracy is strictly tied to the quality of the underlying black-box and additive models; if either underperforms, the entire predictive output suffers. References R. Agarwal, L. Melnick, N. Frosst, X. Zhang, B. Lengerich, R. Caruana, and G. E. Hinton (2021) Neural additive models: interpretable machine learning with neural nets. In Advances in Neural Information Processing Systems, Vol. 34, p. 4699–4711. Cited by: §2. E. Angelino, N. Larus-Stone, D. Alabi, M. Seltzer, and C. Rudin (2018) Learning certifiably optimal rule lists for categorical data. External Links: 1704.01701 Cited by: §2. J. Bien and R. Tibshirani (2011) Prototype selection for interpretable classification. The Annals of Applied Statistics 5 (4). External Links: ISSN 1932-6157, Document Cited by: §2. L. Breiman, J. Friedman, R.A. Olshen, and C.J. Stone (1984) Classification and regression trees. Chapman and Hall/CRC. External Links: Document Cited by: §2. L. M. Candanedo, V. Feldheim, and D. Deramaix (2017) Data driven prediction models of energy use of appliances in a low-energy house. Energy and Buildings 140, p. 81–97. External Links: Document Cited by: §C.1. L. M. Candanedo (2017) Appliances energy prediction dataset. Cited by: §C.1. R. Caruana, Y. Lou, J. Gehrke, P. Koch, M. Sturm, and N. Elhadad (2015) Intelligible models for healthcare: predicting pneumonia risk and hospital 30-day readmission. In Proceedings of the 21th ACM SIGKDD international conference on knowledge discovery and data mining, p. 1721–1730. Cited by: §1. C. Chang, R. Caruana, and A. Goldenberg (2022) NODE-gam: neural generalized additive model for interpretable deep learning. In International Conference on Learning Representations, Cited by: §2, §4. C. Chen, O. Li, C. Tao, A. J. Barnett, J. Su, and C. Rudin (2019) This looks like that: deep learning for interpretable image recognition. In NeurIPS, External Links: 1806.10574 Cited by: §2. F. Doshi-Velez and B. Kim (2017) Towards a rigorous science of interpretable machine learning. arXiv preprint arXiv:1702.08608. Cited by: §1. J. H. Friedman (2001) Greedy function approximation: a gradient boosting machine. Annals of Statistics 29 (5), p. 1189–1232. Cited by: §A.3, §2. V. Gkolemis, T. Dalamagas, and C. Diou (2023a) DALE: differential accumulated local effects for efficient and accurate global explanations. In Asian Conference on Machine Learning, p. 375–390. Cited by: §2. V. Gkolemis, T. Dalamagas, E. Ntoutsi, and C. Diou (2023b) RHALE: robust and heterogeneity-aware accumulated local effects. In ECAI 2023, p. 859–866. Cited by: §2. V. Gkolemis, A. Tzerefos, T. Dalamagas, E. Ntoutsi, and C. Diou (2023c) Regionally additive models: explainable-by-design models minimizing feature interactions. In Joint European Conference on Machine Learning and Knowledge Discovery in Databases, p. 433–447. Cited by: §2. R. Guidotti, A. Monreale, S. Ruggieri, F. Turini, F. Giannotti, and D. Pedreschi (2018) A survey of methods for explaining black box models. ACM Computing Surveys 51 (5), p. 1–42. Cited by: §2. T. Hastie and R. Tibshirani (1990) Generalized additive models. Chapman and Hall/CRC. Cited by: §2. J. Herbinger, B. Bischl, and G. Casalicchio (2022) REPID: regional effect plots with implicit interaction detection. In Proceedings of The 25th International Conference on ArtificialIntelligence and Statistics, G. Camps-Valls, F. J. R. Ruiz, and I. Valera (Eds.), Proceedings of Machine Learning Research, Vol. 151, p. 10209–10233. Cited by: §A.2.1, §1, §2, §3.2. J. Herbinger, S. Dandl, F. K. Ewald, S. Loibl, and G. Casalicchio (2023) Leveraging model-based trees as interpretable surrogate models for model distillation. In European Conference on Artificial Intelligence, p. 232–249. Cited by: §2. J. Herbinger, M. N. Wright, T. Nagler, B. Bischl, and G. Casalicchio (2024) Decomposing global feature effects based on feature interactions. Journal of Machine Learning Research 25 (381), p. 1–65. Cited by: §A.2.1, §1, §2, §3.2, §3.2. G. Hinton, O. Vinyals, and J. Dean (2015) Distilling the knowledge in a neural network. In NeurIPS Deep Learning and Representation Learning Workshop, Cited by: §2. N. Hollmann, S. Müller, K. Eggensperger, and F. Hutter (2022) Tabpfn: a transformer that solves small tabular classification problems in a second. arXiv preprint arXiv:2207.01848. Cited by: §A.1, §3.2. L. Hu, J. Chen, V. N. Nair, and A. Sudjianto (2020) Surrogate locally-interpretable models with supervised machine learning algorithms. arXiv preprint arXiv:2007.14528. Cited by: §2. S. Ibrahim, G. Afriat, K. Behdin, and R. Mazumder (2023) GRAND-SLAMIN’ Interpretable Additive Modeling with Structural Constraints. Advances in Neural Information Processing Systems 36, p. 61158–61186 (en). Cited by: §2. M. Kelly, R. Longjohn, and K. Nottingham (2025) The uci machine learning repository. Note: Accessed: [Insert Date Here] Cited by: §C.1. M. Kraus, D. Tschernutter, S. Weinzierl, and P. Zschech (2024) Interpretable generalized additive neural networks. European Journal of Operational Research 317 (2), p. 303–316. External Links: ISSN 0377-2217, Document Cited by: §2. B. Letham, C. Rudin, T. H. McCormick, and D. Madigan (2015) Interpretable classifiers using rules and bayesian analysis: building a better stroke prediction model. The Annals of Applied Statistics 9 (3). External Links: ISSN 1932-6157, Document Cited by: §2. Z. C. Lipton (2018) The mythos of model interpretability: in machine learning, the concept of interpretability is both important and slippery.. Queue 16 (3), p. 31–57. Cited by: §1. Y. Lou, R. Caruana, J. Gehrke, and G. Hooker (2013) Accurate intelligible models with pairwise interactions. In Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining, p. 623–631. Cited by: §2. Y. Lou, R. Caruana, and J. Gehrke (2012) Intelligible models for classification and regression. In Proceedings of the 18th ACM SIGKDD international conference on Knowledge discovery and data mining, p. 150–158. Cited by: §A.3, §1, §2. C. Molnar (2020) Interpretable machine learning. Cited by: §1. W. J. Murdoch, C. Singh, K. Kumbier, R. Abbasi-Asl, and B. Yu (2019) Interpretable machine learning: definitions, methods, and applications. arXiv preprint arXiv:1901.04592. Cited by: §1. A. Oliabev (2018) Home equity line of credit (HELOC). Note: Accessed: [Insert Date Here] Cited by: §C.1. R. K. Pace and R. Barry (1997) Sparse spatial autoregressions. Statistics and Probability Letters 33, p. 291–297. Cited by: §C.1. ProPublica (2016) COMPAS data and analysis for ’machine bias’. Cited by: §C.1. F. Radenovic, A. Dubey, and D. Mahajan (2022) Neural basis models for interpretability. Advances in Neural Information Processing Systems 35, p. 8414–8426. Cited by: §1, §2. M. T. Ribeiro, S. Singh, and C. Guestrin (2016) “Why should i trust you?”: explaining the predictions of any classifier. In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, p. 1135–1144. Cited by: §2. J. D. Romano, T. T. Le, W. La Cava, J. T. Gregg, D. J. Goldberg, P. Chakraborty, N. L. Ray, D. Himmelstein, W. Fu, and J. H. Moore (2021) PMLB v1.0: an open source dataset collection for benchmarking machine learning methods. arXiv preprint arXiv:2012.00058v2. Cited by: §C.1. C. Rudin (2019) Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead. Nature machine intelligence 1 (5), p. 206–215. Cited by: §1. M. Saeed, M. Villarroel, A. T. Reisner, G. Clifford, L. Lehman, G. Moody, T. Heldt, T. H. Kyaw, B. Moody, and R. G. Mark (2011) Multiparameter intelligent monitoring in intensive care i (mimic-i): a public-access intensive care unit database. Critical Care Medicine 39 (5), p. 952–960. External Links: Document Cited by: §C.1. S. Tan, R. Caruana, G. Hooker, and Y. Lou (2018) Distill-and-compare: auditing black-box models using transparent model distillation. In Proceedings of the 2018 AAAI/ACM Conference on AI, Ethics, and Society, p. 303–310. Cited by: §2. J. Vanschoren, J. N. van Rijn, B. Bischl, and L. Torgo (2013) OpenML: networked science in machine learning. SIGKDD Explorations 15 (2), p. 49–60. External Links: Document Cited by: §C.1. S. N. Wood (2017) Generalized additive models: an introduction with r. CRC press. Cited by: §2, §2. Z. Yang, A. Zhang, and A. Sudjianto (2020) GAMI-net: an explainable neural network based on generalized additive models with structured interactions. arXiv preprint arXiv:2003.07132. Cited by: §2, §4. Appendix A Conditionally Additive Local Models (CALMs) In the main paper (Algorithm 1), we summarize the procedure of fitting a Conditionally Additive Local Model (CALM): fCALM()=g−1(β0+∑i=1dfiri(−i)(xi)),f_ CALM(x)=g^-1 ( _0+ _i=1^df_i^r_i(x_-i)(x_i) ), as defined in Eq.(1), to a dataset =(i),y(i)i=1ND=\x^(i),y^(i)\_i=1^N in three main steps: (i) fit a high‑capacity reference model f_ ref on D, (i) fit a partition tree TiT_i using an interaction-related heterogeneity measure HiH_i for each feature i=1,…,Di=1,…,D, and (i) estimate region‑specific effects fi(r)r=1Ri\f_i^(r)\_r=1^R_i for each feature i=1,…,Di=1,…,D. We provide additional details for each step below. A.1 Step 1: Fit a high‑capacity reference model f_ ref on D The initial stage involves fitting a high-capacity predictive model, denoted as f:ℝD→ℝf_ ref:R^D to dataset =((i),y(i))i=1ND=\(x^(i),y^(i))\_i=1^N. The reference model serves as an accurate functional approximation of the underlying relationship between the input features x and the response variable y. The choice of f_ ref is flexible; any sufficiently expressive and accurate model can be employed. Commonly used options include gradient-boosted decision trees (e.g., XGBoost, which we adopt by default), deep neural networks, random forests, or more recent architectures such as TabPFN (Hollmann et al., 2022). A.2 Step 2: Fit a partition tree TiT_i using heterogeneity HiH_i for each feature We first define a suitable heterogeneity measure HiH_i (see Section A.2.1), and then explore the relationship between heterogeneity and CALM’s error in approximating freff_ ref. Motivated by these results, we fit a partition tree TiT_i for each feature i (see Section A.2.4). Finally, if the user specifies a constraint K on the maximum number of interactions, we apply a pruning step to retain the K most significant interactions (see Section A.2.5). A.2.1 Heterogeneity Measures HiH_i Interaction-related heterogeneity HiH_i quantifies the extent to which a feature xix_i interacts with all other features xjx_j, where j∈1,…,d∖ij∈\1,…,d\ \i\. The computation of HiH_i is based on the variance of local effects: h(xi,−i(j))h(x_i,x_-i^(j)). Local Effect. We define h(xi,−i(j))h(x_i,x_-i^(j)) as the local effect of feature xix_i when applied to the j-th background observation −i(j)x_-i^(j). That is, if we take the j-th observation and substitute the i-th feature with the value xix_i, the resulting change (local effect) in the model output is captured by h(xi,−i(j))h(x_i,x_-i^(j)). The computation of h(⋅)h(·) relies on two components: (i) a reference teacher model f_ ref and (i) a background dataset (i)i=1N\x^(i)\_i=1^N containing input features only. The specific form of h(⋅)h(·) depends on the chosen feature effect method. Several methods can be used, including Partial Dependence Plot (PDP), SHAP-Dependence Plot (SHAP-DP) and Robust and Heterogeneity-aware ALE (RHALE). Below, we provide the definition of h in the case of PDP, which we adopt as the default method in our experiments. Let f_ ref be the reference model and (j)j=1N\x^(j)\_j=1^N the background dataset. The PDP-based heterogeneity is defined as: h(xi,−i(j))=f(xi,−i(j))−c(−i)h (x_i,\,x_-i^(j) )=f_ ref (x_i,\,x_-i^(j) )-c(x_-i) (6) where c(−i(j))=Xi[f(Xi,−i(j))]c(x^(j)_-i)=E_X_i [f_ ref (X_i,\,x_-i^(j) ) ] is a centering constant. For details on alternative approaches, we refer the reader to (Herbinger et al., 2024, 2022). Point-wise heterogeneity Hi(xi)H_i(x_i) and Heterogeneity HiH_i. The pointwise heterogeneity is then defined as Hi(xi)=−i[(h(xi,−i)−μi(xi))2]H_i(x_i)=E_X_-i [ (h(x_i,X_-i)- _i(x_i) )^2 ] (7) where μi(xi)=−i[h(xi,−i)] _i(x_i)=E_X_-i[h(x_i,X_-i)] is the mean local effect. Inside a region, this heterogeneity takes the form of Eq. (4). Feature-level heterogeneity is Hi=Xi[Hi(Xi)]=Xi[Var−i[h(Xi,−i)]]H_i=E_X_i [H_i(X_i) ]=E_X_i [Var_X_-i[h(X_i,X_-i)] ] (8) Inside a region this definition takes the form of Eq. (5). A more formal motivation for the use of heterogeneity in CALM comes from the proposition of the following section. A.2.2 Heterogeneity and CALM approximation error Proposition A.1. Let ℰ=[(fref()−fCALM())2]E=E_X [ (f_ ref(X)-f_ CALM(X) )^2 ] be the mean squared error of the CALM approximation of freff_ ref. Assume that fref:→ℝf_ ref:X admits a multivariate regionally additive approximation fref()≈β0+∑i=1dhiri(−i)()f_ ref(x)≈ _0+ _i=1^dh_i^r_i(x_-i)(x) such that each hiri(−i)h_i^r_i(x_-i) captures the contribution of the i-th feature to the function’s output (including interactions). Also, assume that there is no error in fitting of the univariate CALM shape functions. Then, minfCALMℰ≤d∑i=1d(∑r=1RiPℛi(r)Hiℛi(r)) _f_ CALME≤ d _i=1^d ( _r=1^R_iP_R_i^(r)H^R_i^(r)_i ) (9) where region ℛi(r)R_i^(r) is the r-th region of the i-th feature, Pℛi(r)P_R_i^(r) is the probability mass of the region in the data and Hiℛi(r)H^R_i^(r)_i is the heterogeneity of hiℛi(r)()h_i^R_i^(r)(x). Proof. For the approximation error we have ℰ =[(fref()−fCALM())2]=[((β0+∑i=1dhiri(−i)())−(β0+∑i=1dfiri(−i)(xi)))2] =E_X [ (f_ ref(X)-f_ CALM(X) )^2 ]=E_X [ ( ( _0+ _i=1^dh_i^r_i(x_-i)(x) )- ( _0+ _i=1^df_i^r_i(x_-i)(x_i) ) )^2 ] =[(∑i=1d(hiri(−i)()−firi(−i)(xi)))2]≤d∑i=1d[(hiri(−i)()−firi(−i)(xi))2] =E_X [ ( _i=1^d (h_i^r_i(x_-i)(x)-f_i^r_i(x_-i)(x_i) ) )^2 ]≤ d _i=1^dE_X [ (h_i^r_i(x_-i)(x)-f_i^r_i(x_-i)(x_i) )^2 ] where the last step results from the Jensen’s inequality. We therefore have ℰ≤d∑i=1dℰiri(−i)E≤ d _i=1^dE^r_i(x_-i)_i (10) Next, we show that the minimum error achieved for each region is equal to the heterogeneity of the local effects. For any measurable function gi(xi)g_i(x_i) and region ℛR, decompose, using μiℛ(xi)=−i∣−i∈ℛ[hiℛ(xi,−i)]μ^R_i(x_i)=E_X_-i _-i [h^R_i(x_i,X_-i) ]: hi(Xi,−i)−gi(Xi)=(hi(Xi,−i)−μiℛ(Xi))+(μiℛ(Xi)−gi(Xi))h_i(X_i,X_-i)-g_i(X_i)= (h_i(X_i,X_-i)-μ^R_i(X_i) )+ (μ^R_i(X_i)-g_i(X_i) ) Taking expectations inside ℛR and squaring: Xi,−i∣−i∈ℛ[(hiℛ(Xi,−i)−gi(Xi))2] _X_i,X_-i _-i [ (h_i^R(X_i,X_-i)-g_i(X_i) )^2 ] =Xi,−i∣−i∈ℛ[(hiℛ(Xi,−i)−μiℛ(Xi))2] =E_X_i,X_-i _-i [ (h_i^R(X_i,X_-i)-μ^R_i(X_i) )^2 ] +Xi[(μiℛ(Xi)−gi(Xi))2] +E_X_i [ (μ^R_i(X_i)-g_i(X_i) )^2 ] +2Xi[−i∣−i∈ℛ[hiℛ(Xi,−i)−μiℛ(Xi)](μiℛ(Xi)−gi(Xi))] +2E_X_i [E_X_-i _-i [h_i^R(X_i,X_-i)-μ^R_i(X_i) ] (μ^R_i(X_i)-g_i(X_i) ) ] The cross-term equals zero because: −i∣−i∈ℛ[hiℛ(xi,−i)−μiℛ(xi)]=μiℛ(xi)−μiℛ(xi)=0E_X_-i _-i [h_i^R(x_i,X_-i)-μ^R_i(x_i) ]=μ^R_i(x_i)-μ^R_i(x_i)=0 by definition of μiℛμ^R_i. Therefore: Xi,−i∣−i∈ℛ[(hiℛ(Xi,−i)−gi(Xi))2]=ℰiℛ+Xi[(μiℛ(Xi)−gi(Xi))2]≥ℰiℛE_X_i,X_-i _-i [ (h_i^R(X_i,X_-i)-g_i(X_i) )^2 ]=E^R_i+E_X_i [ (μ^R_i(X_i)-g_i(X_i) )^2 ] ^R_i with equality if and only if gi=μiℛg_i=μ^R_i almost surely. This indicates that within each region, μiℛ _i^R minimizes the approximation error. Moreover, from the definition of pointwise heterogeneity: Hiℛ(xi):=−i∣−i∈ℛ[(hℛ(xi,−i)−μiℛ(xi))2]H^R_i(x_i):=E_X_-i _-i [ (h^R(x_i,X_-i)-μ^R_i(x_i) )^2 ] Taking the expectation over XiX_i: Hiℛ=Xi[Hiℛ(Xi)]=Xi[−i∣−i∈ℛ[(hℛ(Xi,−i)−μiℛ(Xi))2]]H^R_i=E_X_i [H^R_i(X_i) ]=E_X_i [E_X_-i _-i [ (h^R(X_i,X_-i)-μ^R_i(X_i) )^2 ] ] which gives Hiℛ=Xi,−i|−i∈ℛ[(hℛ(Xi,−i)−μiℛ(Xi))2]=minℰiℛH^R_i=E_X_i,X_-i|X_-i [ (h^R(X_i,X_-i)-μ^R_i(X_i) )^2 ]= _i^R i.e., heterogeneity is the minimum approximation error. Since heterogeneity minimizes the error for each region we can estimate the feature-level error across regions as minℰi=∑r=1RiPℛi(r)Hiℛi(r) E_i= _r=1^R_iP_R_i^(r)H^R_i^(r)_i where Pℛi(r)P_R_i^(r) is the probability mass of region ℛi(r)R_i^(r). Using this result with Eq. (10), we have minfCALMℰ≤d∑i=1d(∑r=1RiPℛi(r)Hiℛi(r)) _f_ CALME≤ d _i=1^d ( _r=1^R_iP_R_i^(r)H^R_i^(r)_i ) ∎ This result links heterogeneity with the estimation error of the local effect inside a region and motivates the CALM approach for region splitting. The initial assumption about the additive approximation of freff_ ref is a common approach followed by feature effect methods (where the effect of each feature is computed independently). Moreover, the selected local effect function and heterogeneity are solely used for identifying the splitting regions and not as estimators of freff_ ref. Using this result, and especially Eq. (9), the following two sections present the CALM approach for using heterogeneity to define a partition for each feature. A.2.3 Heterogeneity estimation For heterogeneity estimation, denote with ℐ⊆1,…,NI \1,…,N\ the index set of active background instances considered in the heterogeneity calculation i.e., those that reside in region ℛi(ri(−i))R^(r_i(x_-i))_i of the i-th instance. Then, the point-wise heterogeneity at feature value xix_i over ℐI can be estimated by H^iℐ(xi)=1|ℐ|∑j∈ℐ(h(xi,−i(j))−μ(xi))2, where μ(xi)=1|ℐ|∑j∈ℐh(xi,−i(j)) H_i^I(x_i)= 1|I| _j (h(x_i,x_-i^(j))-μ(x_i) )^2, where μ(x_i)= 1|I| _j h(x_i,x_-i^(j)) (11) Eq.(11) quantifies the strength of the interactions of the i-th feature, at position xix_i, as the variance of the local effects at xix_i. To obtain a global measure of these interactions, we average he pointwise heterogeneity over a grid of M values x~i(m)m=1M\ x_i^(m)\_m=1^M sampled from the domain of xix_i: H^iℐ=1M∑m=1MHiℐ(x~i(m)). H_i^I= 1M _m=1^MH_i^I( x_i^(m)). (12) A.2.4 Computing the partition tree TiT_i for each feature. Algorithm 2 Fitting a Partition Tree TiT_i for Feature xix_i 1: Input: freff_ref, =((k),y(k))k=1ND=\(x^(k),y^(k))\_k=1^N; Parameters: max depth dmaxd_ , threshold ϵε, grid size M=20M=20 2: Output: Binary tree TiT_i 3: 4: function BuildTree(node ν, depth ℓ , active indices ℐI) 5: if ℓ=dmax =d_ then 6: return ν Stop splitting 7: end if 8: Compute HiℐH_i^I using Eq.(12) 9: Initialize: ΔHmax←0 H_ ← 0 10: for each feature j≠ij≠ i do 11: Determine candidate thresholds jT_j Default: Unique values if xjx_j categorical, else M-grid (default: M=20M=20) 12: for each τ∈jτ _j do 13: if xjx_j is numerical then 14: ℐL←k∈ℐ:xj(k)≤τI_L←\k :x_j^(k)≤τ\, ℐR←ℐ∖ℐLI_R _L 15: else xjx_j is categorical 16: ℐL←k∈ℐ:xj(k)=τI_L←\k :x_j^(k)=τ\, ℐR←ℐ∖ℐLI_R _L 17: end if 18: Compute ΔHi H_i using Eq. (13) 19: if ΔHi>ΔHmax H_i> H_ then 20: ΔHmax←ΔHi H_ ← H_i, Store (j,τ)(j,τ) as optimal split 21: end if 22: end for 23: end for 24: if ΔHmax>ϵ H_ >ε then 25: Split node ν into νL,νR _L, _R using optimal (j,τ)(j,τ) 26: νL← _L← BuildTree(νL,ℓ+1,ℐL _L, +1,I_L) 27: νR← _R← BuildTree(νR,ℓ+1,ℐR _R, +1,I_R) 28: end if 29: return ν 30: end function For each feature xix_i, we construct a binary tree TiT_i of maximum depth dmaxd_ . At each internal node of the tree, we evaluate candidate binary splits based on all features xjx_j for j∈1,…,d∖ij∈\1,…,d\ \i\. For each candidate splitting feature xjx_j, we consider a fixed number T of candidate threshold values τ. These candidate thresholds are selected as T equally spaced values over the range of xjx_j in the training data, i.e., [minkxj(k),maxkxj(k)][ _kx_j^(k), _kx_j^(k)], where xj(k)x_j^(k) denotes the value of feature xjx_j for the k-th training instance. Alternatively, users may choose to evaluate all unique observed values of xjx_j, i.e., xj(k)k=1N\x_j^(k)\_k=1^N. To determine the best split at each node, we compute the heterogeneity drop, which quantifies the decrease in interaction-related heterogeneity of the target feature xix_i after performing the split. Specifically, for a candidate split, we partition the current set of instances ℐI into left and right subsets ℐLI_L and ℐRI_R, based on whether the splitting feature xjx_j falls below or above the threshold τ. The heterogeneity drop is defined as: ΔHi=Hiℐ−(|ℐL||ℐ|HiℐL+|ℐR||ℐ|HiℐR)Hiℐ H_i^= H_i^I- ( |I_L||I|H_i^I_L+ |I_R||I|H_i^I_R )H_i^I (13) where HiℐH_i^I denotes the interaction-related heterogeneity of feature xix_i over the current set of instances ℐI; ℐLI_L and ℐRI_R are the subsets of ℐI resulting from the candidate split; HiℐLH_i^I_L and HiℐRH_i^I_R are the heterogeneities of xix_i computed on the left and right subsets, respectively; and |ℐL|/|ℐ||I_L|/|I| and |ℐR|/|ℐ||I_R|/|I| are the proportions of instances in each subset, which serve as weights in the weighted average. This normalized metric enhances interpretability by allowing users to define a meaningful threshold, ϵε, for split acceptance. For instance, a user may require a candidate split to achieve a minimum heterogeneity reduction of ϵ=0.2ε=0.2 (representing a 20%20\% drop) to be considered significant. In our experimental setup, we utilize a default threshold of ϵ=0.2ε=0.2. The procedure for fitting the tree TiT_i is given in Algorithm 2. A.2.5 Prune trees to keep top K interactions If the user sets an upper threshold K on the number of interactions, we want to select the K most important splits. Given the trees TiT_i for i=1,…,di=1,…,d, we denote by ΔHiν H_i^ν the heterogeneity reduction associated with node ν in tree TiT_i. We collect all such splits ΔHiν H_i^ν across all features i and nodes ν, and then sort them in descending order. From this ordered list, we select the top K nodes that correspond to the largest heterogeneity decreases, while ensuring that no child node is retained without its parent node also being included. Based on this selection, we prune each tree TiT_i by removing nodes outside the retained set, thereby preserving only the most significant interactions. A.3 Step 3: Estimating the region‑specific effects fi(r)r=1Ri\f_i(r)\_r=1^R_i for each feature. The original gradient boosting procedure for fitting a standard, global GAM, fGAM()=g−1(β0+∑ifi(xi))f_ GAM(x)=g^-1 ( _0+ _if_i(x_i) ), follows the round-robin approach of (Friedman, 2001). At each boosting iteration, one univariate shape function fif_i is updated by fitting to the current residuals over the entire dataset, thereby greedily reducing the overall loss. The procedure as describe by (Lou et al., 2012) is summarized in Algorithm 3. Algorithm 3 Gradient Boosting for GAM 1: Initialize fi←0f_i← 0 for all i=1,…,di=1,…,d 2: for m=1m=1 to M do 3: for i=1i=1 to d do 4: ℰ←(xj(i),y(j)−fGAM((j))j=1NE←\(x^(i)_j,y^(j)-f_GAM(x^(j))\_j=1^N Compute partial residuals 5: Learn shaping function S:xi→ℝS:x_i using ℰE Fit shape function 6: fi←fi+ηSf_i← f_i+η S Update with learning rate η 7: end for 8: end for To fit a CALM: fCALM()=g−1(β0+∑i=1dfiri(−i)(xi))f_ CALM(x)=g^-1 ( _0+ _i=1^df_i^r_i(x_-i)(x_i) ) we customize the standard gradient boosting to accommodate for region-specific effects, i.e., we adapt this scheme so that each shape function fi(r)f_i^(r) is trained only on the subset of instances that lie within its assigned subregion ℛi(r)R_i^(r). This preserves the additive, boosting framework while enforcing that each fi(r)f_i^(r) captures only the behavior of xix_i within its corresponding region. The procedure is summarized in Algorithm 4. Algorithm 4 Gradient Boosting for CALM 1: Initialize fi(r)←0f_i^(r)← 0 for all i=1,…,di=1,…,d and r=1,…,Rir=1,…,R_i 2: for m=1m=1 to M do 3: for i=1i=1 to d do 4: for r=1r=1 to RiR_i do 5: ℰ←(xi(j),y(j)−fCALM((j))):(j)∈Regionr(Ti)E←\(x^(j)_i,y^(j)-f_CALM(x^(j))):x^(j) _r(T_i)\ Filter by region ℛi(r)R_i^(r) 6: Learn shaping function S:xi→ℝS:x_i using ℰE 7: fi(r)←fi(r)+ηSf_i^(r)← f_i^(r)+η S Update regional shape function 8: end for 9: end for 10: end for A.3.1 Theoretical analysis of Step 3 This section formalizes the claims made in Step 3 of the main text. We show that, under squared loss and fixed region-selection trees, (i) the population target of Step 3 is an L2(P)L_2(P_X)-best approximation of the true regression function onto the CALM function class, and (i) an idealized exact version of the algorithm converges to the corresponding empirical risk minimizer. Consider the CALM score function s()=β0+∑i=1dfi(ri(−i))(xi),s(x)= _0+ _i=1^df_i^(r_i(X_-i))(x_i), (14) where the region-selection functions ri\r_i\ (or equivalently the trees Ti\T_i\) are treated as fixed. Let m()=[Y∣=]m(x)=E[Y =x] denote the true regression function. We assume regression with identity link and squared loss, and [Y2]<∞E[Y^2]<∞. We restrict attention to scores s∈L2(P)s∈ L_2(P_X) (hence m∈L2(P)m∈ L_2(P_X)). To ensure identifiability, we impose the centering condition [fi(r)(Xi)|ri(−i)=r]=0for all i,r with ℙ(ri(−i)=r)>0.E\! [f_i^(r)(X_i)\, |\,r_i(X_-i)=r ]=0 all i,r with P(r_i(X_-i)=r)>0. (15) This convention assigns all region-wise constants to the intercept β0 _0 and ensures uniqueness of the decomposition. Let ℋ(Ti)H(\T_i\) denote the class of all CALM score functions of the form (14) satisfying (15), and assume ℋ(Ti)⊂L2(P)H(\T_i\)⊂ L_2(P_X) is a nonempty closed convex set. A.3.2 Population target Proposition A.2 (Population optimality). There exists at least one minimizer s⋆∈argmins∈ℋ(Ti)[(Y−s())2],s ∈ _s (\T_i\)E[(Y-s(X))^2], and it satisfies s⋆∈argmins∈ℋ(Ti)[(m()−s())2].s ∈ _s (\T_i\)E[(m(X)-s(X))^2]. In particular, [(Y−s⋆())2]−[(Y−m())2]=[(m()−s⋆())2].E[(Y-s (X))^2]-E[(Y-m(X))^2]=E[(m(X)-s (X))^2]. Proof. For any measurable function s, Y−s()=(Y−m())+(m()−s()).Y-s(X)=(Y-m(X))+(m(X)-s(X)). Taking squares and expectations yields [(Y−s())2]=[(Y−m())2]+[(m()−s())2],E[(Y-s(X))^2]=E[(Y-m(X))^2]+E[(m(X)-s(X))^2], since [Y−m()∣]=0E[Y-m(X) ]=0 eliminates the cross-term. The first term does not depend on s, so minimizing prediction error over ℋ(Ti)H(\T_i\) is equivalent to minimizing [(m()−s())2]E[(m(X)-s(X))^2]. Existence of s⋆s follows since ℋ(Ti)⊂L2(P)H(\T_i\)⊂ L_2(P_X) is nonempty, closed, and convex. Condition (15) ensures uniqueness of the representation (β0,fi(r))( _0,\f_i^(r)\) for any given s∈ℋ(Ti)s (\T_i\). ∎ Interpretation. Proposition A.2 shows that, with fixed regions, Step 3 computes an L2(P)L_2(P_X)-best approximation of the true regression function m within the CALM function class. Any remaining error is therefore purely due to the structural limitations of the chosen regions and univariate components, not to the optimization procedure. A.3.3 Empirical convergence of exact Step 3 Given i.i.d. data (x(k),y(k))k=1N\(x^(k),y^(k))\_k=1^N, define the empirical risk ℛ^N(s)=1N∑k=1N(y(k)−s(x(k)))2. R_N(s)= 1N _k=1^N(y^(k)-s(x^(k)))^2. Proposition A.3 (Convergence of exact cyclic regional backfitting). Assume fixed trees and squared loss. Consider an idealized version of Step 3 that cyclically updates each fi(r)f_i^(r) by exactly minimizing ℛ^N R_N over that function using only samples with ri(−i(k))=r_i(x^(k)_-i)=r, followed by empirical centering implemented by shifting the region-wise sample mean from fi(r)f_i^(r) into the intercept β0 _0, so that fitted values (and hence ℛ^N R_N) are unchanged. Then ℛ^N R_N is non-increasing along the iterates, and the procedure converges (in empirical L2L_2) to a limit s^∈argmins∈ℋ(Ti)ℛ^N(s) s∈ _s (\T_i\) R_N(s). Moreover, the fitted-value vector (s^((1)),…,s^((N)))( s(x^(1)),…, s(x^(N))) is unique. Proof. With fixed trees, the CALM score is linear in the collection of shape functions fi(r)\f_i^(r)\ evaluated on the data, and ℛ^N R_N is a convex quadratic function of these parameters. On the finite sample, ℛ^N R_N depends on each fi(r)f_i^(r) only through the finite vector fi(r)(xi(k)):ri(−i(k))=r\f_i^(r)(x_i^(k)):r_i(x_-i^(k))=r\, so the problem can be viewed as a finite-dimensional least-squares problem. Each regional update solves the exact least-squares problem for one block of parameters while holding the others fixed, and therefore cannot increase ℛ^N R_N. Although regions overlap across features, all updates jointly optimize a single global objective. Standard results for cyclic block coordinate descent on convex quadratics imply convergence to a global minimizer, and the minimizing fitted values are unique by strict convexity of the quadratic loss in the fitted-value vector. The centering condition (15) ensures a unique decomposition into components. ∎ Interpretation. This result shows that Step 3 is not a collection of independent local fits. Instead, it performs a coupled optimization of a single additive predictor with region-gated components, analogous to classical backfitting for GAMs. Under exact updates, this procedure is guaranteed to converge to the best CALM fit for the given data and fixed regions. Appendix B Formal Proofs for Interpretability Properties This appendix provides formal characterizations and proofs for the interpretability properties discussed in Section 3.3. For each property, we define it precisely and analyze whether it can be answered exactly under GAM , GA2M, and CALM. Proposition B.1 (Local Feature Contribution). Let y^:ℝd→ℝ y:R^d be a prediction function, and fix an input =(x1,…,xd)∈ℝdx=(x_1,…,x_d) ^d. For each index i∈1,…,di∈\1,…,d\, we define the local contribution ϕi()∈ℝ _i(x) as a function intended to represent the contribution of feature xix_i to the output y^() y(x). We analyze whether such a decomposition y^()=∑i=1dϕi() y(x)= _i=1^d _i(x) is uniquely determined by the structure of the model. (A) In GAM . Assume the model has the additive form y^()=∑j=1dfj(xj), y(x)= _j=1^df_j(x_j), with each fj:ℝ→ℝf_j:R . Then for all i∈1,…,di∈\1,…,d\, ϕi():=fi(xi) _i(x):=f_i(x_i) is well-defined, and y^()=∑i=1dϕi(). y(x)= _i=1^d _i(x). (B) In GA2M. Assume the model includes pairwise interactions: y^()=∑j=1dfj(xj)+∑j<kfjk(xj,xk), y(x)= _j=1^df_j(x_j)+ _j<kf_jk(x_j,x_k), where each fj:ℝ→ℝf_j:R and fjk:ℝ2→ℝf_jk:R^2 . Then in general, there does not exist a unique additive decomposition y^()=∑i=1dϕi() y(x)= _i=1^d _i(x) with each ϕi() _i(x) depending only on xix_i. (C) In CALM. Assume the model has the form y^()=∑j=1dfj(rj(−j))(xj), y(x)= _j=1^df_j^(r_j(x_-j))(x_j), where for each j, rj:ℝd−1→1,…,Rjr_j:R^d-1→\1,…,R_j\ is a region selector and fj(r):ℝ→ℝf_j^(r):R is a univariate shape function for region r. Here −j:=∖xj∈ℝd−1x_-j:=x x_j ^d-1. Then for all i∈1,…,di∈\1,…,d\, ϕi():=fi(ri(−i))(xi) _i(x):=f_i^(r_i(x_-i))(x_i) is well-defined, and y^()=∑i=1dϕi(). y(x)= _i=1^d _i(x). Proof. We examine each case separately. (A) In GAM . By direct substitution, the prediction is an additive sum of univariate functions, each depending only on a single feature xjx_j. Thus the contribution of feature xix_i is uniquely defined as fi(xi)f_i(x_i), and the decomposition is exact. (B) In GA2M. Suppose, for contradiction, that there exists a decomposition y^()=∑i=1dϕi(), y(x)= _i=1^d _i(x), where each ϕi:ℝd→ℝ _i:R^d represents the contribution of feature xix_i and depends only on xix_i, i.e., ϕi()=ϕi(xi) _i(x)= _i(x_i). Then each interaction term fjk(xj,xk)f_jk(x_j,x_k) must be split between ϕj _j and ϕk _k in such a way that the total sum remains correct and additive over individual features. This is only possible if fjkf_jk is additively separable, i.e., if there exist functions gj,gk:ℝ→ℝg_j,g_k:R such that fjk(xj,xk)=gj(xj)+gk(xk).f_jk(x_j,x_k)=g_j(x_j)+g_k(x_k). However, the model does not constrain fjkf_jk to be separable. In the general case, fjkf_jk is non-separable and depends jointly on xjx_j and xkx_k. Therefore, no decomposition ∑iϕi(xi) _i _i(x_i) can reproduce the prediction y^() y(x) using only univariate terms. The contributions ϕi() _i(x) are not uniquely determined by the model, and must rely on external assumptions or post-hoc attribution methods. (C) In CALM. For a fixed input x, each region function rj(−j)r_j(x_-j) returns a unique region index in 1,…,Rj\1,…,R_j\. The corresponding function fj(rj(−j))f_j^(r_j(x_-j)) is applied to xjx_j, yielding a scalar. Thus, the contribution of feature xix_i is precisely defined as: ϕi():=fi(ri(−i))(xi), _i(x):=f_i^(r_i(x_-i))(x_i), and the model prediction is recovered by summing over all features: y^()=∑i=1dfi(ri(−i))(xi)=∑i=1dϕi(). y(x)= _i=1^df_i^(r_i(x_-i))(x_i)= _i=1^d _i(x). ∎ The second property asks: How does the prediction change if we perturb a single feature xix_i while keeping all others fixed? The answer differs across GAM , GA2M, and CALM, as shown below. Proposition B.2 (Regional Feature Sensitivity)). Let y^:ℝd→ℝ y:R^d be a prediction function. Fix an input ∈ℝdx ^d, an index i∈1,…,di∈\1,…,d\, and a perturbation ε>0 >0. Define the prediction change under perturbation of xix_i as Δy^:=y^(+εi)−y^(). y:= y(x+ _i)- y(x). We analyze how Δy y can be computed in GAM , GA2M, and CALM. (A) In GAM . Assume the model has the additive form y^()=∑j=1dfj(xj), y(x)= _j=1^df_j(x_j), with each fj:ℝ→ℝf_j:R univariate. Then Δy^=fi(xi+ε)−fi(xi). y=f_i(x_i+ )-f_i(x_i). (B) In GA2M. We show that in GA2M, the change Δy y caused by perturbing xix_i depends on the values of all interacting features xjx_j for j≠ij≠ i, and thus cannot be computed from xix_i alone. (C) In CALM: We show that the change in prediction Δy^:=y^(+εi)−y^() y:= y(x+ _i)- y(x) is computable from xix_i alone if and only if no region transitions are triggered in other features. Otherwise, the change depends on −ix_-i, and regional sensitivity is not determined solely by xix_i. Proof. (A) In GAM . Since all other features remain unchanged, and their contributions are independent of xix_i, we have y^(+εi)=∑j≠ifj(xj)+fi(xi+ε), y(x+ _i)= _j≠ if_j(x_j)+f_i(x_i+ ), y^()=∑j≠ifj(xj)+fi(xi), y(x)= _j≠ if_j(x_j)+f_i(x_i), and therefore Δy^=fi(xi+ε)−fi(xi). y=f_i(x_i+ )-f_i(x_i). (B) In GA2M. Let :=j≠i∣fij(xi,xj) is present in the modelJ:=\j≠ i f_ij(x_i,x_j) is present in the model\, i.e., the set of features that interact with xix_i. Since only xix_i is perturbed, all other features remain unchanged. The only affected terms are the univariate function fi(xi)f_i(x_i), and the interaction terms fij(xi,xj)f_ij(x_i,x_j) for each j∈j . All other terms in the model remain constant. The prediction after perturbation is: y^(+εi)=fi(xi+ε)+∑j∈fj(xj)+∑j∈fij(xi+ε,xj)+C, y(x+ _i)=f_i(x_i+ )+ _j f_j(x_j)+ _j f_ij(x_i+ ,x_j)+C, and before perturbation: y^()=fi(xi)+∑j∈fj(xj)+∑j∈fij(xi,xj)+C, y(x)=f_i(x_i)+ _j f_j(x_j)+ _j f_ij(x_i,x_j)+C, where C:=∑j<k,j,k≠ifjk(xj,xk)C:= _j<k,\,j,k≠ if_jk(x_j,x_k) denotes the interaction terms not involving xix_i, which are unaffected. Subtracting, we obtain: Δy^=fi(xi+ε)−fi(xi)+∑j∈[fij(xi+ε,xj)−fij(xi,xj)]. y=f_i(x_i+ )-f_i(x_i)+ _j [f_ij(x_i+ ,x_j)-f_ij(x_i,x_j) ]. Since this expression depends on interacting xjx_j for j≠ij≠ i, it cannot be computed from xix_i alone. (C) In CALM. Assume the model has the form y^()=∑j=1dfj(rj(−j))(xj), y(x)= _j=1^df_j^(r_j(x_-j))(x_j), where rj:ℝd−1→1,…,Rjr_j:R^d-1→\1,…,R_j\ assigns a region index to each feature j based on the context −j:=∖xjx_-j:=x x_j, and fj(r):ℝ→ℝf_j^(r):R is the shape function for region r. Let −i=∖xix_-i=x x_i, and define the perturbation xi↦xi+εx_i x_i+ . Then Δy^=fi(ri(−i))(xi+ε)−fi(ri(−i))(xi)+∑j≠i[fj(rj(−j+ε))(xj)−fj(rj(−j))(xj)], y=f_i^(r_i(x_-i))(x_i+ )-f_i^(r_i(x_-i))(x_i)+ _j≠ i [f_j^(r_j(x_-j^+ ))(x_j)-f_j^(r_j(x_-j))(x_j) ], where −j+ε:=−jx_-j^+ :=x_-j with xi↦xi+εx_i x_i+ (since xi∈−jx_i _-j for all j≠ij≠ i). Case 1: No region transitions in other features. Assume that for all j≠ij≠ i, we have rj(−j+ε)=rj(−j)r_j(x_-j^+ )=r_j(x_-j). Then: Δy^=fi(ri(−i))(xi+ε)−fi(ri(−i))(xi). y=f_i^(r_i(x_-i))(x_i+ )-f_i^(r_i(x_-i))(x_i). This case is fully interpretable from the function fi(r)f_i^(r) alone. Case 2: Region transitions occur in other features. If for some j≠ij≠ i, the region function changes: rj(−j+ε)≠rj(−j)r_j(x_-j^+ )≠ r_j(x_-j), then the additional terms contribute: Δfj:=fj(rj(−j+ε))(xj)−fj(rj(−j))(xj). f_j:=f_j^(r_j(x_-j^+ ))(x_j)-f_j^(r_j(x_-j))(x_j). Hence, the total change includes not only the direct shift in fif_i, but also discrete region-based changes in other features and the value of Δy y cannot be recovered from xix_i alone. ∎ The third property asks whether increasing a feature always increases the model’s prediction, regardless of the values of other features. In other words, does the model treat the feature as globally monotonic? Proposition B.3 (Global Feature Monotonicity). Let y^:ℝd→ℝ y:R^d be a prediction function, and fix a feature index i∈1,…,di∈\1,…,d\. Define the global monotonicity condition for feature xix_i as follows: the model is said to be globally increasing in xix_i if for all ∈ℝdx ^d and all δ>0δ>0, y^(+δi)≥y^(), y(x+ _i)≥ y(x), and strictly increasing if the inequality is strict. We analyze whether this property can be verified from the structure of the model. (A) In GAM . Assume the model has the form: y^()=∑j=1dfj(xj), y(x)= _j=1^df_j(x_j), with each fj:ℝ→ℝf_j:R . Then the model is globally increasing in xix_i if and only if fif_i is monotonically increasing. (B) In GA2M. Assume the model has the form: y^()=∑j=1dfj(xj)+∑j<kfjk(xj,xk). y(x)= _j=1^df_j(x_j)+ _j<kf_jk(x_j,x_k). Then global monotonicity in xix_i cannot be determined from fif_i alone. (C) In CALM. Assume the model has the form: y^()=∑j=1dfj(rj(−j))(xj), y(x)= _j=1^df_j^(r_j(x_-j))(x_j), where rj:ℝd−1→1,…,Rjr_j:R^d-1→\1,…,R_j\ is a region selector, and −j:=∖xjx_-j:=x x_j. Then global monotonicity in xix_i holds if the following two conditions are satisfied: 1. For all regions r∈1,…,Rir∈\1,…,R_i\, the function fi(r):ℝ→ℝf_i^(r):R is increasing. 2. For all j≠ij≠ i, and all xj∈ℝx_j , the function xi↦fj(rj(−j))(xj)x_i f_j^(r_j(x_-j))(x_j) is non-decreasing; i.e., region transitions caused by changing xix_i do not decrease fjf_j’s contribution. Proof. (A) In GAM . Since the model is additive and each term depends only on a single variable, we have: y^(+δi)−y^()=fi(xi+δ)−fi(xi). y(x+ _i)- y(x)=f_i(x_i+δ)-f_i(x_i). Thus, y^(+δi)≥y^() y(x+ _i)≥ y(x) if and only if fif_i is increasing. (B) In GA2M. As in case (A), fi(xi+δ)−fi(xi)f_i(x_i+δ)-f_i(x_i) gives the direct contribution. However, for each j≠ij≠ i, the term fij(xi,xj)f_ij(x_i,x_j) also contributes. The total change is: y^(+δi)−y^()=fi(xi+δ)−fi(xi)+∑j≠i[fij(xi+δ,xj)−fij(xi,xj)]. y(x+ _i)- y(x)=f_i(x_i+δ)-f_i(x_i)+ _j≠ i [f_ij(x_i+δ,x_j)-f_ij(x_i,x_j) ]. The sign of this expression depends on the values of xjx_j, and therefore cannot be determined from fif_i alone. Consequently, verifying global monotonicity in xix_i requires knowing the full interaction structure and input values. (C) In CALM Let ∈ℝdx ^d and δ>0δ>0. Define: Δy^:=y^(+δi)−y^(). y:= y(x+ _i)- y(x). Then: Δy^=fi(ri(−i))(xi+δ)−fi(ri(−i))(xi)+∑j≠i[fj(rj(−j+δ))(xj)−fj(rj(−j))(xj)], y=f_i^(r_i(x_-i))(x_i+δ)-f_i^(r_i(x_-i))(x_i)+ _j≠ i [f_j^(r_j(x_-j^+δ))(x_j)-f_j^(r_j(x_-j))(x_j) ], where −j+δx_-j^+δ is obtained by replacing xi↦xi+δx_i x_i+δ in −jx_-j. The first term is non-negative if fi(r)f_i^(r) is increasing for all regions r, and the region ri(−i)r_i(x_-i) is fixed. The second term is a sum over changes in other features’ contributions due to possible changes in their regions rjr_j. For the model to be globally increasing in xix_i, all such contributions must be non-negative. This requires that increasing xix_i does not cause a decrease in the contribution of any other feature xjx_j, i.e., region transitions must preserve or increase each fjf_j. Therefore, global monotonicity in xix_i requires both conditions above. Conversely, if both conditions are satisfied, then each term in the sum is non-negative, implying Δy^≥0 y≥ 0. ∎ Appendix C Detailed Experiments C.1 Experimental Setup Details Reporting convention. As in Section 4, we report mean ± standard deviation over 5-fold cross-validation. Table 4: Classification Datasets Dataset Size Attributes Adult 45222 13 COMPAS 6167 9 HELOC 10459 23 MIMIC2 24508 17 Appendicitis 106 7 Phoneme 5404 5 SPECTF 349 44 Magic 19020 10 Bank 45211 16 Churn 5000 19 Table 5: Regression Datasets Dataset Size Attributes Bike Sharing 17379 11 California Housing 20640 8 Parkinsons Motor 5875 19 Parkinsons Total 5875 19 Seoul Bike 8465 14 Wine 6497 11 Energy 19735 28 CCPP 9568 4 Electrical 10000 13 Elevators 16599 18 No2 500 7 Sensory 576 11 Airfoil 1503 5 Skill Craft 3338 19 Ailerons 13750 39 Datasets details. Tables 5 and 5 summarize the 10 classification and 15 regression datasets used in our experiments. These datasets are sourced from various publicly available repositories. The UCI Machine Learning Repository (Kelly et al., 2025) provides datasets including Adult, Magic, Bank, Bike Sharing, Parkinson’s Motor, Parkinson’s Total, Seoul Bike, Wine, CCPP and Skill Craft. OpenML (Vanschoren et al., 2013) offers datasets such as Electrical, Elevators, No2, Sensory, Airfoil, and Ailerons. The Penn Machine Learning Benchmarks (PMLB) (Romano et al., 2021) includes datasets like Appendicitis, Phoneme, SPECTF, and Churn. Additional datasets used in our experiments include California Housing (Pace and Barry, 1997), MIMIC-I (Saeed et al., 2011), COMPAS (ProPublica, 2016), HELOC (Oliabev, 2018), and Energy (Candanedo, 2017; Candanedo et al., 2017). Data Preprocessing. All input features are standardized using z-score normalization. For regression tasks, the target variable y is also standard scaled. The only exception is the GAMI-Net model, which requires input features to be scaled using min-max normalization. For RMSE, predictions and targets are inverse-transformed to the original scale prior to evaluation. Model Configurations. Our black box models comprise a fully connected neural network, a random forest and an XGBoost ensemble. The deep network contains two hidden layers with 50 units each; ReLU activations are used for regression whereas a final sigmoid unit closes the classification variant. Training is carried out for 200 epochs with the Adam optimiser (learning rate=0.001), a batch size of 200 and either mean-squared error or binary cross-entropy loss, depending on the task. The random-forest baseline consists of 500 trees grown to a maximum depth of 25 with a minimum of three samples required at each leaf. For boosted trees we employ XGBoost with 300 boosting rounds, a learning rate of 0.1 and the log-loss evaluation metric for classification For GAM models we consider NAM and EBM without interactions and Spline. The Neural Additive Model (NAM) builds one independent multilayer perceptron per feature; each sub-network has three hidden layers of 100,100 and 10 ReLU units followed by a linear output, and the additive sum is passed through a sigmoid for binary classification. Training uses Adam (learning rate 0.001), for 10 epochs and a mini-batch size of 32. Spline-GAMs are implemented with PyGAM’s LinearGAM or LogisticGAM, allocating a single spline term to every feature while relying on PyGAM’s automatic smoothing. Finally, the Explainable Boosting Machine (EBM) is used without interaction terms by explicitly setting n_interactions=0n\_interactions=0. To capture pairwise interactions we experiment with three GA2M variants. EB2M extends the EBM by enabling interactions with a default strength of 0.9; NODE‑GA2M activates its interaction mode via the ga2m=1 flag and restricts training to a five‑minute time limit to ensure parity with the other models; GAMI‑Net is run with the default hyper‑parameters. For CALM , we configure the heterogeneity drop threshold to 0.2, determining whether a split is considered statistically significant. The region detector, which measures heterogeneity , relies on Partial Dependence Plots (PDP) for RF and XGBoost black-box models, while for DNNs we consider both PDP and RHALE to measure heterogeneity. The detector evaluates 20 candidate split points per feature. The GAM used within regions is an EBM without interactions and the CALM is applied on top of all three black-box models (DNN, RF, XGBoost). Compatibility. The GAMI-Net baseline depends on a native binary (lib_ebmcore_mac_x64.dylib) compiled for x86_64 architecture. As a result, it is not compatible with Apple Silicon Macs, and running it on such systems will lead to an architecture mismatch error. This issue is limited to this external model and does not affect any of the proposed methods or other baselines. We provide instructions in the README file to guide such users on running all other methods except GAMI-Net. Computer Resources. All experiments were conducted on an in-house server with cloud infrastructure equipped with an Intel(R) Core(TM) i9-10900X CPU @ 3.70GHz, 128 GB of RAM. No GPU acceleration was utilized during these experiments. C.2 Number of Interactions Number of Interactions (Classification). Table 6 reports the number of feature interactions selected or used by each model across classification datasets. CALM consistently uses significantly fewer interactions than GA2M-style baselines, often by a wide margin. In most datasets, CALM activates less than a third of the interactions compared to NodeGA2M, and even fewer than GAMINet’s default of 20. The numbers shown for CALM reflect actual detected regions with interaction-specific behavior, confirming that its compact regionalization mechanism results in sparse and interpretable models. Number of Interactions (Regression). As shown in Table 7, CALM activates fewer feature interactions than other GA2M-based models in regression tasks. While GAMI-Net and EB2M use a fixed or near-fixed interaction set (e.g., 20 by default), and NodeGA2M often selects over 100 interactions, CALM remains sparse across datasets. Table 6: Number of Interactions per Model in Classification Datasets Dataset CALM EB2M NodeGA2M GAMINet Adult 4.2±1.44.2 $±$1.4 12.0±0.012.0 $±$0.0 72.6±3.272.6 $±$3.2 20.0±0.020.0 $±$0.0 COMPAS 1.4±1.71.4 $±$1.7 9.0±0.09.0 $±$0.0 36.0±0.036.0 $±$0.0 20.0±0.020.0 $±$0.0 HELOC 1.1±1.61.1 $±$1.6 21.0±0.021.0 $±$0.0 163.8±4.4163.8 $±$4.4 20.0±0.020.0 $±$0.0 MIMIC2 11.5±2.111.5 $±$2.1 16.0±0.016.0 $±$0.0 115.0±3.3115.0 $±$3.3 20.0±0.020.0 $±$0.0 Appendicitis 4.0±1.74.0 $±$1.7 7.0±0.07.0 $±$0.0 21.0±0.021.0 $±$0.0 6.2±8.16.2 $±$8.1 Phoneme 8.6±1.68.6 $±$1.6 5.0±0.05.0 $±$0.0 10.0±0.010.0 $±$0.0 10.0±0.010.0 $±$0.0 SPECTF 0.0±0.00.0 $±$0.0 40.0±0.040.0 $±$0.0 273.0±14.4273.0 $±$14.4 20.0±0.020.0 $±$0.0 Magic 8.5±2.48.5 $±$2.4 9.0±0.09.0 $±$0.0 44.0±0.644.0 $±$0.6 20.0±0.020.0 $±$0.0 Bank 9.2±2.09.2 $±$2.0 15.0±0.015.0 $±$0.0 104.2±3.2104.2 $±$3.2 20.0±0.020.0 $±$0.0 Churn 14.0±0.014.0 $±$0.0 18.0±0.018.0 $±$0.0 137.4±2.9137.4 $±$2.9 20.0±0.020.0 $±$0.0 Avg. 6.2 15.2 97.7 17.6 Table 7: Number of Interactions per Model in Regression Datasets Dataset CALM EB2M NodeGA2M GAMINet Bike Sharing 19.3±2.719.3 $±$2.7 10.0±0.010.0 $±$0.0 53.6±1.453.6 $±$1.4 20.0±0.020.0 $±$0.0 California Housing 10.7±2.010.7 $±$2.0 8.0±0.08.0 $±$0.0 28.0±0.028.0 $±$0.0 20.0±0.020.0 $±$0.0 Parkinsons Motor 13.8±4.413.8 $±$4.4 18.0±0.018.0 $±$0.0 126.2±1.7126.2 $±$1.7 20.0±0.020.0 $±$0.0 Parkinsons Total 14.4±5.814.4 $±$5.8 18.0±0.018.0 $±$0.0 129.2±3.9129.2 $±$3.9 20.0±0.020.0 $±$0.0 Seoul Bike 20.4±1.020.4 $±$1.0 13.0±0.013.0 $±$0.0 83.2±0.783.2 $±$0.7 20.0±0.020.0 $±$0.0 Wine 13.9±3.313.9 $±$3.3 10.0±0.010.0 $±$0.0 54.6±0.554.6 $±$0.5 20.0±0.020.0 $±$0.0 Energy 49.2±4.449.2 $±$4.4 26.0±0.026.0 $±$0.0 192.0±4.5192.0 $±$4.5 4.0±8.04.0 $±$8.0 CCPP 1.9±2.21.9 $±$2.2 4.0±0.04.0 $±$0.0 6.0±0.06.0 $±$0.0 5.8±0.45.8 $±$0.4 Electrical 0.0±0.00.0 $±$0.0 12.0±0.012.0 $±$0.0 68.6±3.068.6 $±$3.0 12.0±0.012.0 $±$0.0 Elevators 11.7±3.011.7 $±$3.0 17.0±0.017.0 $±$0.0 115.4±2.1115.4 $±$2.1 20.0±0.020.0 $±$0.0 No2 12.8±1.912.8 $±$1.9 7.0±0.07.0 $±$0.0 21.0±0.021.0 $±$0.0 20.0±0.020.0 $±$0.0 Sensory 6.9±4.26.9 $±$4.2 10.0±0.010.0 $±$0.0 54.0±0.954.0 $±$0.9 16.0±4.916.0 $±$4.9 Airfoil 9.3±0.59.3 $±$0.5 5.0±0.05.0 $±$0.0 10.0±0.010.0 $±$0.0 10.0±0.010.0 $±$0.0 Skill Craft 0.0±0.00.0 $±$0.0 18.0±0.018.0 $±$0.0 135.6±2.4135.6 $±$2.4 16.0±8.016.0 $±$8.0 Ailerons 41.7±2.641.7 $±$2.6 36.0±0.036.0 $±$0.0 232.4±5.1232.4 $±$5.1 20.0±0.020.0 $±$0.0 Avg. 15.1 14.1 87.3 16.3 C.3 Practical Runtime Table 8: Runtime (Seconds) for Classification Datasets Dataset BlackBox GAM CALM GA2M XGB NAM EBM EB2M NodeGA2M GAMINet Adult 0.2±0.010.2 $±$0.01 38±238 $±$2 10±110 $±$1 36±136 $±$1 13±213 $±$2 97±997 $±$9 718±122718 $±$122 COMPAS 0.1±0.0040.1 $±$0.004 10±0.110 $±$0.1 2±22 $±$2 2±0.12 $±$0.1 1±0.041 $±$0.04 51±151 $±$1 129±5129 $±$5 HELOC 0.2±0.010.2 $±$0.01 31±631 $±$6 2±22 $±$2 17±117 $±$1 2±0.22 $±$0.2 57±0.157 $±$0.1 249±60249 $±$60 MIMIC2 0.2±0.010.2 $±$0.01 31±0.331 $±$0.3 4±34 $±$3 26±326 $±$3 5±15 $±$1 69±969 $±$9 359±32359 $±$32 Appendicitis 0.1±0.0010.1 $±$0.001 5±0.15 $±$0.1 2±42 $±$4 2±0.12 $±$0.1 1±0.51 $±$0.5 24±0.124 $±$0.1 6±56 $±$5 Phoneme 0.1±0.010.1 $±$0.01 6±16 $±$1 3±43 $±$4 3±33 $±$3 4±14 $±$1 61±461 $±$4 166±38166 $±$38 SPECTF 0.1±0.010.1 $±$0.01 32±0.332 $±$0.3 3±43 $±$4 3±33 $±$3 4±54 $±$5 25±0.225 $±$0.2 180±13180 $±$13 Magic 0.2±0.0030.2 $±$0.003 28±1828 $±$18 5±55 $±$5 17±317 $±$3 7±47 $±$4 98±798 $±$7 489±63489 $±$63 Bank 0.2±0.040.2 $±$0.04 44±144 $±$1 5±15 $±$1 42±342 $±$3 13±113 $±$1 98±1798 $±$17 994±265994 $±$265 Churn 0.2±0.010.2 $±$0.01 19±319 $±$3 2±12 $±$1 10±410 $±$4 3±23 $±$2 56±256 $±$2 257±32257 $±$32 Avg. 0.2 24 4 16 5 64 355 Table 9: Runtime (Seconds) for Regression Datasets Dataset BlackBox GAM CALM GA2M XGB NAM EBM EB2M NodeGA2M GAMINet Bike Sharing 0.1±0.010.1 $±$0.01 18±118 $±$1 2±12 $±$1 15±115 $±$1 19±319 $±$3 187±31187 $±$31 299±33299 $±$33 California Housing 0.2±0.010.2 $±$0.01 15±0.115 $±$0.1 4±24 $±$2 13±113 $±$1 16±0.516 $±$0.5 200±46200 $±$46 677±147677 $±$147 Parkinsons Motor 0.3±0.0030.3 $±$0.003 21±621 $±$6 5±25 $±$2 17±217 $±$2 44±444 $±$4 193±34193 $±$34 536±88536 $±$88 Parkinsons Total 0.3±0.010.3 $±$0.01 18±0.218 $±$0.2 5±35 $±$3 16±216 $±$2 40±540 $±$5 227±76227 $±$76 543±79543 $±$79 Seoul Bike 0.2±0.010.2 $±$0.01 15±0.215 $±$0.2 3±33 $±$3 13±313 $±$3 12±412 $±$4 127±18127 $±$18 268±53268 $±$53 Wine 0.2±0.010.2 $±$0.01 12±112 $±$1 2±32 $±$3 9±39 $±$3 3±0.43 $±$0.4 62±1062 $±$10 187±20187 $±$20 Energy 0.3±0.010.3 $±$0.01 43±0.443 $±$0.4 8±38 $±$3 82±582 $±$5 97±1597 $±$15 206±57206 $±$57 253±214253 $±$214 CCPP 0.2±0.0040.2 $±$0.004 6±0.16 $±$0.1 5±45 $±$4 8±58 $±$5 14±214 $±$2 133±18133 $±$18 88±888 $±$8 Electrical 0.1±0.0020.1 $±$0.002 25±2025 $±$20 9±49 $±$4 11±511 $±$5 8±38 $±$3 112±28112 $±$28 279±43279 $±$43 Elevators 0.3±0.10.3 $±$0.1 25±0.425 $±$0.4 7±47 $±$4 35±535 $±$5 43±843 $±$8 181±38181 $±$38 773±269773 $±$269 No2 0.1±0.0030.1 $±$0.003 5±0.045 $±$0.04 2±0.12 $±$0.1 2±0.12 $±$0.1 2±0.12 $±$0.1 26±0.126 $±$0.1 78±178 $±$1 Sensory 0.1±0.0010.1 $±$0.001 10±110 $±$1 4±64 $±$6 7±57 $±$5 2±0.42 $±$0.4 27±0.527 $±$0.5 74±1674 $±$16 Airfoil 0.1±0.0010.1 $±$0.001 4±0.034 $±$0.03 3±63 $±$6 5±65 $±$6 7±17 $±$1 51±1251 $±$12 76±3476 $±$34 Skill Craft 0.3±0.0030.3 $±$0.003 17±0.217 $±$0.2 4±64 $±$6 9±69 $±$6 1±0.11 $±$0.1 44±0.144 $±$0.1 141±52141 $±$52 Ailerons 0.3±0.020.3 $±$0.02 53±153 $±$1 5±75 $±$7 70±870 $±$8 3±0.13 $±$0.1 126±36126 $±$36 619±134619 $±$134 Avg. 0.2 19 5 21 21 127 326 We report the total runtime of each method across all datasets in Tables 8 and 9. For CALM, the reported time includes not only the regions detection fitting steps, but also the training time of the underlying black-box model and the GAM component used within regions. While this inclusion gives a complete view of end-to-end cost, it somewhat disadvantages CALM in comparison to other models, whose runtimes reflect only their own training processes. In practice, the black-box or GAM components can be selected to be lightweight and the blackbox model can be reused or pretrained independently, making CALM’s region discovery step lightweight and modular. Despite this conservative accounting, CALM’s runtime remains competitive. For example, it often trains faster than or comparably to GA2M methods such as NODE-GA2M and GAMI-Net, which tend to have high computational overhead. On small and medium-sized datasets like COMPAS, HELOC, California Housing, and Wine, CALM runs in under a minute, showing that its interpretability gains come with reasonable computational cost. On larger datasets such as Bank or Ailerons, CALM’s runtime scales moderately but remains practical. C.4 Performance Metrics for all Experiments Performance Summary. We evaluated CALM on 25 datasets using three model classes: DNN, XGB, and RF. For each model and dataset, we compare CALM to both its corresponding GAM baseline and the original black-box model. In comparisons with GAMs, we match the model structure: for example, CALM-NAM is compared directly to NAM, so the performance difference reflects only the benefit of introducing regional modeling. In comparisons with the black-box, we report the accuracy of the best CALM-GAM. For XGB and RF, this is selected across multiple GAM types (NAM, EBM, Spline) using a fixed heterogeneity threshold of 0.2. For DNNs, we first identify the best heterogeneity modeling method (PDP or RHALE) for each GAM, and then select the best-performing GAM among all types. While the threshold is fixed in our experiments, we note that alternative thresholds could yield further improvements. The results cover six tables: Tables 10–12 report performance on regression datasets using DNN, XGB, and RF respectively, while Tables 13–15 report on the same three model classes for classification datasets. Across these settings, CALM consistently outperforms the corresponding GAM baseline. It matches or exceeds the performance of its GAM counterpart in 97.8% and 100% of cases for DNNs (Tables 10, 13), 97.8% and 96.7% for XGB (Tables 11, 14), and 95.6% and 96.7% for RF (Tables 12, 15). These results reflect direct, structure-matched comparisons between CALM and GAMs of the same architecture. When compared to the black-box models, CALM also achieves strong results. The best CALM-GAM variant outperforms the original black-box in 73.3% of DNN regression cases (Table 10), 33.3% for XGB regression (Table 11), and 26.7% for RF regression (Table 12). In classification, it beats the black-box in 90.0% of DNN cases (Table 13), 50.0% for XGB (Table 14), and 60.0% for RF (Table 15). Overall, CALM outperforms or matches the corresponding GAM in 219 out of 225 cases (97.3%) and surpasses the black-box in 40 out of 75 cases (53.3%). The rare cases where CALM does not improve over the GAM occur in both regression and classification tasks and almost always correspond to scenarios where the GAM itself performs as well as or better than the black-box. These cases suggest that the underlying function is already well captured by global additive effects, leaving little room for further refinement through regional modeling. Table 10: RMSE Score for Regression Datasets, DNN Dataset BlackBox GAM CALM - NAM CALM - EBM CALM - Spline DNN NAM EBM Spline PDP RHALE PDP RHALE PDP RHALE Bike Sharing 42.952±1.46042.952 $±$1.460 101.968±1.309101.968 $±$1.309 100.211±1.010100.211 $±$1.010 100.876±0.868100.876 $±$0.868 63.034±0.88363.034 $±$0.883 64.941±2.56664.941 $±$2.566 58.444±1.63758.444 $±$1.637 60.966±2.98360.966 $±$2.983 59.868±1.56459.868 $±$1.564 62.108±2.98962.108 $±$2.989 California Housing 0.519±0.0170.519 $±$0.017 0.611±0.0110.611 $±$0.011 0.554±0.0080.554 $±$0.008 0.632±0.0090.632 $±$0.009 0.580±0.0120.580 $±$0.012 0.591±0.0160.591 $±$0.016 0.522±0.0100.522 $±$0.010 0.526±0.0090.526 $±$0.009 0.593±0.0180.593 $±$0.018 0.597±0.0060.597 $±$0.006 Parkinsons Motor 3.654±0.1623.654 $±$0.162 6.112±0.1586.112 $±$0.158 4.204±0.0874.204 $±$0.087 6.027±0.0706.027 $±$0.070 4.220±0.1924.220 $±$0.192 5.303±0.3125.303 $±$0.312 2.681±0.1052.681 $±$0.105 3.658±0.3243.658 $±$0.324 4.149±0.2114.149 $±$0.211 5.318±0.4305.318 $±$0.430 Parkinsons Total 5.139±0.1025.139 $±$0.102 7.897±0.1027.897 $±$0.102 4.847±0.1074.847 $±$0.107 7.519±0.0727.519 $±$0.072 5.456±0.1825.456 $±$0.182 6.776±0.3816.776 $±$0.381 3.587±0.0303.587 $±$0.030 4.319±0.2224.319 $±$0.222 5.408±0.3655.408 $±$0.365 6.827±0.2096.827 $±$0.209 Seoul Bike 244.967±5.081244.967 $±$5.081 320.225±4.377320.225 $±$4.377 303.685±3.855303.685 $±$3.855 302.020±4.701302.020 $±$4.701 269.277±6.414269.277 $±$6.414 273.241±8.379273.241 $±$8.379 249.457±2.407249.457 $±$2.407 245.880±5.153245.880 $±$5.153 246.855±3.936246.855 $±$3.936 243.926±4.433243.926 $±$4.433 Wine 0.700±0.0150.700 $±$0.015 0.719±0.0090.719 $±$0.009 0.703±0.0120.703 $±$0.012 0.736±0.0390.736 $±$0.039 0.712±0.0050.712 $±$0.005 0.711±0.0120.711 $±$0.012 0.695±0.0110.695 $±$0.011 0.695±0.0120.695 $±$0.012 0.729±0.0270.729 $±$0.027 0.734±0.0350.734 $±$0.035 Energy 77.139±2.91677.139 $±$2.916 89.799±2.66589.799 $±$2.665 85.225±2.48985.225 $±$2.489 86.498±2.38486.498 $±$2.384 89.304±2.48089.304 $±$2.480 88.600±2.17988.600 $±$2.179 84.347±2.35884.347 $±$2.358 83.972±2.55783.972 $±$2.557 84.929±1.86884.929 $±$1.868 83.526±2.65383.526 $±$2.653 CCPP 3.916±0.0793.916 $±$0.079 4.213±0.0504.213 $±$0.050 3.435±0.0763.435 $±$0.076 4.144±0.0544.144 $±$0.054 4.205±0.0674.205 $±$0.067 4.188±0.0794.188 $±$0.079 3.427±0.0583.427 $±$0.058 3.379±0.0833.379 $±$0.083 3.995±0.0613.995 $±$0.061 4.006±0.0584.006 $±$0.058 Electrical 0.051±0.0090.051 $±$0.009 0.040±0.0050.040 $±$0.005 0.018±0.0110.018 $±$0.011 0.095±0.0050.095 $±$0.005 0.040±0.0050.040 $±$0.005 0.040±0.0050.040 $±$0.005 0.018±0.0110.018 $±$0.011 0.012±0.0070.012 $±$0.007 0.094±0.0050.094 $±$0.005 0.088±0.0030.088 $±$0.003 Elevators 0.002±0.0000.002 $±$0.000 0.003±0.0000.003 $±$0.000 0.002±0.0000.002 $±$0.000 0.002±0.0000.002 $±$0.000 0.002±0.0000.002 $±$0.000 0.002±0.0000.002 $±$0.000 0.002±0.0000.002 $±$0.000 0.002±0.0000.002 $±$0.000 0.002±0.0000.002 $±$0.000 0.002±0.0000.002 $±$0.000 No2 0.522±0.0120.522 $±$0.012 0.503±0.0210.503 $±$0.021 0.492±0.0340.492 $±$0.034 0.480±0.0320.480 $±$0.032 0.507±0.0160.507 $±$0.016 0.501±0.0250.501 $±$0.025 0.492±0.0300.492 $±$0.030 0.498±0.0340.498 $±$0.034 0.485±0.0280.485 $±$0.028 0.492±0.0330.492 $±$0.033 Sensory 0.522±0.0220.522 $±$0.022 0.483±0.0140.483 $±$0.014 0.477±0.0090.477 $±$0.009 0.482±0.0120.482 $±$0.012 0.461±0.0250.461 $±$0.025 0.449±0.0190.449 $±$0.019 0.447±0.0210.447 $±$0.021 0.444±0.0120.444 $±$0.012 0.453±0.0250.453 $±$0.025 0.446±0.0200.446 $±$0.020 Airfoil 2.131±0.1312.131 $±$0.131 4.801±0.1984.801 $±$0.198 4.565±0.1484.565 $±$0.148 4.542±0.1214.542 $±$0.121 3.275±0.1963.275 $±$0.196 3.143±0.1503.143 $±$0.150 2.647±0.1942.647 $±$0.194 2.449±0.2272.449 $±$0.227 2.697±0.1612.697 $±$0.161 2.588±0.2482.588 $±$0.248 Skill Craft 1.231±0.1951.231 $±$0.195 0.925±0.0250.925 $±$0.025 0.901±0.0210.901 $±$0.021 2.515±3.0922.515 $±$3.092 0.923±0.0280.923 $±$0.028 0.918±0.0230.918 $±$0.023 0.901±0.0210.901 $±$0.021 0.905±0.0180.905 $±$0.018 1.652±1.3231.652 $±$1.323 2.345±2.7922.345 $±$2.792 Ailerons 0.0002±0.0000.0002 $±$0.000 0.0002±0.0000.0002 $±$0.000 0.0002±0.0000.0002 $±$0.000 0.0002±0.0000.0002 $±$0.000 0.0002±0.0000.0002 $±$0.000 0.0002±0.0000.0002 $±$0.000 0.0002±0.0000.0002 $±$0.000 0.0002±0.0000.0002 $±$0.000 0.0002±0.0000.0002 $±$0.000 0.001±0.0010.001 $±$0.001 Table 11: RMSE Score for Regression Datasets, XGB Dataset BlackBox GAM CALM XGB NAM EBM Spline NAM EBM Spline Bike Sharing 39.346±1.37539.346 $±$1.375 101.968±1.309101.968 $±$1.309 100.211±1.010100.211 $±$1.010 100.876±0.868100.876 $±$0.868 59.732±1.41459.732 $±$1.414 55.667±1.22055.667 $±$1.220 56.656±1.32056.656 $±$1.320 California Housing 0.454±0.0100.454 $±$0.010 0.611±0.0110.611 $±$0.011 0.554±0.0080.554 $±$0.008 0.632±0.0090.632 $±$0.009 0.566±0.0100.566 $±$0.010 0.511±0.0100.511 $±$0.010 0.594±0.0140.594 $±$0.014 Parkinsons Motor 1.443±0.0941.443 $±$0.094 6.112±0.1586.112 $±$0.158 4.204±0.0874.204 $±$0.087 6.027±0.0706.027 $±$0.070 4.574±0.1374.574 $±$0.137 2.243±0.1262.243 $±$0.126 4.432±0.1854.432 $±$0.185 Parkinsons Total 1.863±0.0791.863 $±$0.079 7.897±0.1027.897 $±$0.102 4.847±0.1074.847 $±$0.107 7.519±0.0727.519 $±$0.072 6.109±0.3616.109 $±$0.361 2.968±0.0892.968 $±$0.089 5.638±0.2885.638 $±$0.288 Seoul Bike 209.587±3.474209.587 $±$3.474 320.225±4.377320.225 $±$4.377 303.685±3.855303.685 $±$3.855 302.020±4.701302.020 $±$4.701 270.365±7.049270.365 $±$7.049 238.912±1.643238.912 $±$1.643 237.705±2.658237.705 $±$2.658 Wine 0.622±0.0120.622 $±$0.012 0.719±0.0090.719 $±$0.009 0.703±0.0120.703 $±$0.012 0.736±0.0390.736 $±$0.039 0.703±0.0110.703 $±$0.011 0.693±0.0160.693 $±$0.016 0.737±0.0430.737 $±$0.043 Energy 67.967±2.57967.967 $±$2.579 89.799±2.66589.799 $±$2.665 85.225±2.48985.225 $±$2.489 86.498±2.38486.498 $±$2.384 87.416±2.89487.416 $±$2.894 83.088±1.97483.088 $±$1.974 82.747±2.42182.747 $±$2.421 CCPP 3.086±0.0903.086 $±$0.090 4.213±0.0504.213 $±$0.050 3.435±0.0763.435 $±$0.076 4.144±0.0544.144 $±$0.054 4.171±0.0424.171 $±$0.042 3.419±0.0663.419 $±$0.066 4.039±0.0364.039 $±$0.036 Electrical 0.037±0.0150.037 $±$0.015 0.040±0.0050.040 $±$0.005 0.018±0.0110.018 $±$0.011 0.095±0.0050.095 $±$0.005 0.040±0.0050.040 $±$0.005 0.018±0.0110.018 $±$0.011 0.095±0.0050.095 $±$0.005 Elevators 0.002±0.0000.002 $±$0.000 0.003±0.0000.003 $±$0.000 0.002±0.0000.002 $±$0.000 0.002±0.0000.002 $±$0.000 0.002±0.0000.002 $±$0.000 0.002±0.0000.002 $±$0.000 0.002±0.0000.002 $±$0.000 No2 0.473±0.0260.473 $±$0.026 0.503±0.0210.503 $±$0.021 0.492±0.0340.492 $±$0.034 0.480±0.0320.480 $±$0.032 0.498±0.0220.498 $±$0.022 0.492±0.0360.492 $±$0.036 0.479±0.0310.479 $±$0.031 Sensory 0.508±0.0280.508 $±$0.028 0.483±0.0140.483 $±$0.014 0.477±0.0090.477 $±$0.009 0.482±0.0120.482 $±$0.012 0.462±0.0160.462 $±$0.016 0.450±0.0200.450 $±$0.020 0.451±0.0290.451 $±$0.029 Airfoil 1.544±0.1001.544 $±$0.100 4.801±0.1984.801 $±$0.198 4.565±0.1484.565 $±$0.148 4.542±0.1214.542 $±$0.121 3.146±0.1373.146 $±$0.137 2.503±0.1472.503 $±$0.147 2.582±0.1302.582 $±$0.130 Skill Craft 0.938±0.0190.938 $±$0.019 0.925±0.0250.925 $±$0.025 0.901±0.0210.901 $±$0.021 2.515±3.0922.515 $±$3.092 0.925±0.0250.925 $±$0.025 0.901±0.0210.901 $±$0.021 2.167±2.2612.167 $±$2.261 Ailerons 0.0002±0.0000.0002 $±$0.000 0.0002±0.0000.0002 $±$0.000 0.0002±0.0000.0002 $±$0.000 0.0002±0.0000.0002 $±$0.000 0.0002±0.0000.0002 $±$0.000 0.0002±0.0000.0002 $±$0.000 0.0002±0.0000.0002 $±$0.000 Table 12: RMSE Score for Regression Datasets, RF Dataset BlackBox GAM CALM RF NAM EBM Spline NAM EBM Spline Bike Sharing 43.247±1.03743.247 $±$1.037 101.968±1.309101.968 $±$1.309 100.211±1.010100.211 $±$1.010 100.876±0.868100.876 $±$0.868 61.174±1.54461.174 $±$1.544 57.158±1.28857.158 $±$1.288 57.823±1.39357.823 $±$1.393 California Housing 0.502±0.0120.502 $±$0.012 0.611±0.0110.611 $±$0.011 0.554±0.0080.554 $±$0.008 0.632±0.0090.632 $±$0.009 0.570±0.0150.570 $±$0.015 0.515±0.0120.515 $±$0.012 0.611±0.0250.611 $±$0.025 Parkinsons Motor 1.477±0.1551.477 $±$0.155 6.112±0.1586.112 $±$0.158 4.204±0.0874.204 $±$0.087 6.027±0.0706.027 $±$0.070 4.515±0.0944.515 $±$0.094 2.230±0.0812.230 $±$0.081 4.349±0.0614.349 $±$0.061 Parkinsons Total 1.869±0.2431.869 $±$0.243 7.897±0.1027.897 $±$0.102 4.847±0.1074.847 $±$0.107 7.519±0.0727.519 $±$0.072 5.761±0.2595.761 $±$0.259 3.044±0.0683.044 $±$0.068 5.393±0.1625.393 $±$0.162 Seoul Bike 225.572±3.086225.572 $±$3.086 320.225±4.377320.225 $±$4.377 303.685±3.855303.685 $±$3.855 302.020±4.701302.020 $±$4.701 270.514±8.146270.514 $±$8.146 243.910±2.580243.910 $±$2.580 240.471±4.337240.471 $±$4.337 Wine 0.618±0.0120.618 $±$0.012 0.719±0.0090.719 $±$0.009 0.703±0.0120.703 $±$0.012 0.736±0.0390.736 $±$0.039 0.708±0.0060.708 $±$0.006 0.691±0.0110.691 $±$0.011 0.726±0.0360.726 $±$0.036 Energy 69.406±2.78169.406 $±$2.781 89.799±2.66589.799 $±$2.665 85.225±2.48985.225 $±$2.489 86.498±2.38486.498 $±$2.384 88.448±2.60788.448 $±$2.607 84.741±2.52084.741 $±$2.520 84.690±2.54684.690 $±$2.546 CCPP 3.378±0.0783.378 $±$0.078 4.213±0.0504.213 $±$0.050 3.435±0.0763.435 $±$0.076 4.144±0.0544.144 $±$0.054 4.129±0.0574.129 $±$0.057 3.422±0.1123.422 $±$0.112 4.029±0.0524.029 $±$0.052 Electrical 0.009±0.0090.009 $±$0.009 0.040±0.0050.040 $±$0.005 0.018±0.0110.018 $±$0.011 0.095±0.0050.095 $±$0.005 0.040±0.0050.040 $±$0.005 0.018±0.0110.018 $±$0.011 0.095±0.0050.095 $±$0.005 Elevators 0.003±0.0000.003 $±$0.000 0.003±0.0000.003 $±$0.000 0.002±0.0000.002 $±$0.000 0.002±0.0000.002 $±$0.000 0.002±0.0000.002 $±$0.000 0.002±0.0000.002 $±$0.000 0.002±0.0000.002 $±$0.000 No2 0.466±0.0330.466 $±$0.033 0.503±0.0210.503 $±$0.021 0.492±0.0340.492 $±$0.034 0.480±0.0320.480 $±$0.032 0.481±0.0230.481 $±$0.023 0.482±0.0310.482 $±$0.031 0.484±0.0290.484 $±$0.029 Sensory 0.446±0.0070.446 $±$0.007 0.483±0.0140.483 $±$0.014 0.477±0.0090.477 $±$0.009 0.482±0.0120.482 $±$0.012 0.451±0.0220.451 $±$0.022 0.439±0.0200.439 $±$0.020 0.442±0.0250.442 $±$0.025 Airfoil 2.196±0.1702.196 $±$0.170 4.801±0.1984.801 $±$0.198 4.565±0.1484.565 $±$0.148 4.542±0.1214.542 $±$0.121 3.073±0.1433.073 $±$0.143 2.533±0.1372.533 $±$0.137 2.627±0.1082.627 $±$0.108 Skill Craft 0.916±0.0190.916 $±$0.019 0.925±0.0250.925 $±$0.025 0.901±0.0210.901 $±$0.021 2.515±3.0922.515 $±$3.092 0.926±0.0370.926 $±$0.037 0.901±0.0220.901 $±$0.022 2.320±2.6962.320 $±$2.696 Ailerons 0.0002±0.0000.0002 $±$0.000 0.0002±0.0000.0002 $±$0.000 0.0002±0.0000.0002 $±$0.000 0.0002±0.0000.0002 $±$0.000 0.0002±0.0000.0002 $±$0.000 0.0002±0.0000.0002 $±$0.000 0.0002±0.0000.0002 $±$0.000 Table 13: Accuracy Score for Classification Datasets, DNN Dataset BlackBox GAM CALM - NAM CALM - EBM CALM - Spline DNN NAM EBM Spline PDP RHALE PDP RHALE PDP RHALE Adult 0.849±0.0020.849 $±$0.002 0.851±0.0020.851 $±$0.002 0.870±0.0020.870 $±$0.002 0.856±0.0010.856 $±$0.001 0.853±0.0010.853 $±$0.001 0.854±0.0030.854 $±$0.003 0.870±0.0030.870 $±$0.003 0.871±0.0030.871 $±$0.003 0.857±0.0030.857 $±$0.003 0.859±0.0030.859 $±$0.003 COMPAS 0.682±0.0110.682 $±$0.011 0.682±0.0160.682 $±$0.016 0.681±0.0120.681 $±$0.012 0.685±0.0120.685 $±$0.012 0.682±0.0070.682 $±$0.007 0.681±0.0140.681 $±$0.014 0.683±0.0110.683 $±$0.011 0.679±0.0120.679 $±$0.012 0.683±0.0150.683 $±$0.015 0.683±0.0110.683 $±$0.011 HELOC 0.721±0.0090.721 $±$0.009 0.723±0.0110.723 $±$0.011 0.728±0.0140.728 $±$0.014 0.726±0.0100.726 $±$0.010 0.723±0.0110.723 $±$0.011 0.723±0.0130.723 $±$0.013 0.727±0.0150.727 $±$0.015 0.728±0.0140.728 $±$0.014 0.721±0.0120.721 $±$0.012 0.726±0.0110.726 $±$0.011 MIMIC2 0.885±0.0030.885 $±$0.003 0.886±0.0010.886 $±$0.001 0.886±0.0030.886 $±$0.003 0.886±0.0030.886 $±$0.003 0.887±0.0010.887 $±$0.001 0.886±0.0010.886 $±$0.001 0.886±0.0030.886 $±$0.003 0.886±0.0020.886 $±$0.002 0.886±0.0030.886 $±$0.003 0.886±0.0030.886 $±$0.003 Appendicitis 0.877±0.0720.877 $±$0.072 0.848±0.0560.848 $±$0.056 0.877±0.0770.877 $±$0.077 0.887±0.0640.887 $±$0.064 0.839±0.0250.839 $±$0.025 0.848±0.0560.848 $±$0.056 0.897±0.0550.897 $±$0.055 0.877±0.0770.877 $±$0.077 0.887±0.0640.887 $±$0.064 0.887±0.0640.887 $±$0.064 Phoneme 0.815±0.0090.815 $±$0.009 0.808±0.0020.808 $±$0.002 0.821±0.0070.821 $±$0.007 0.832±0.0060.832 $±$0.006 0.839±0.0120.839 $±$0.012 0.834±0.0080.834 $±$0.008 0.858±0.0100.858 $±$0.010 0.856±0.0080.856 $±$0.008 0.855±0.0090.855 $±$0.009 0.852±0.0060.852 $±$0.006 SPECTF 0.814±0.0160.814 $±$0.016 0.839±0.0300.839 $±$0.030 0.894±0.0150.894 $±$0.015 0.845±0.0400.845 $±$0.040 0.848±0.0360.848 $±$0.036 0.882±0.0250.882 $±$0.025 0.885±0.0210.885 $±$0.021 0.903±0.0140.903 $±$0.014 0.848±0.0320.848 $±$0.032 0.782±0.0360.782 $±$0.036 Magic 0.870±0.0040.870 $±$0.004 0.850±0.0060.850 $±$0.006 0.857±0.0050.857 $±$0.005 0.855±0.0040.855 $±$0.004 0.860±0.0070.860 $±$0.007 0.858±0.0040.858 $±$0.004 0.862±0.0060.862 $±$0.006 0.860±0.0060.860 $±$0.006 0.865±0.0050.865 $±$0.005 0.860±0.0050.860 $±$0.005 Bank 0.903±0.0030.903 $±$0.003 0.901±0.0030.901 $±$0.003 0.902±0.0020.902 $±$0.002 0.903±0.0020.903 $±$0.002 0.902±0.0010.902 $±$0.001 0.903±0.0050.903 $±$0.005 0.905±0.0030.905 $±$0.003 0.904±0.0030.904 $±$0.003 0.906±0.0030.906 $±$0.003 0.905±0.0050.905 $±$0.005 Churn 0.938±0.0050.938 $±$0.005 0.885±0.0090.885 $±$0.009 0.886±0.0050.886 $±$0.005 0.889±0.0050.889 $±$0.005 0.957±0.0100.957 $±$0.010 0.949±0.0040.949 $±$0.004 0.953±0.0060.953 $±$0.006 0.946±0.0060.946 $±$0.006 0.959±0.0060.959 $±$0.006 0.941±0.0070.941 $±$0.007 Table 14: Accuracy Score for Classification Datasets, XGB Dataset BlackBox GAM CALM XGB NAM EBM Spline NAM EBM Spline Adult 0.870±0.0020.870 $±$0.002 0.851±0.0020.851 $±$0.002 0.870±0.0020.870 $±$0.002 0.856±0.0010.856 $±$0.001 0.853±0.0020.853 $±$0.002 0.870±0.0020.870 $±$0.002 0.858±0.0010.858 $±$0.001 COMPAS 0.661±0.0070.661 $±$0.007 0.682±0.0160.682 $±$0.016 0.681±0.0120.681 $±$0.012 0.685±0.0120.685 $±$0.012 0.686±0.0140.686 $±$0.014 0.684±0.0130.684 $±$0.013 0.686±0.0160.686 $±$0.016 HELOC 0.717±0.0130.717 $±$0.013 0.723±0.0110.723 $±$0.011 0.728±0.0140.728 $±$0.014 0.726±0.0100.726 $±$0.010 0.724±0.0110.724 $±$0.011 0.728±0.0120.728 $±$0.012 0.726±0.0110.726 $±$0.011 MIMIC2 0.890±0.0010.890 $±$0.001 0.886±0.0010.886 $±$0.001 0.886±0.0030.886 $±$0.003 0.886±0.0030.886 $±$0.003 0.886±0.0010.886 $±$0.001 0.886±0.0030.886 $±$0.003 0.886±0.0030.886 $±$0.003 Appendicitis 0.868±0.0690.868 $±$0.069 0.848±0.0560.848 $±$0.056 0.877±0.0770.877 $±$0.077 0.887±0.0640.887 $±$0.064 0.868±0.0690.868 $±$0.069 0.878±0.0630.878 $±$0.063 0.878±0.0630.878 $±$0.063 Phoneme 0.898±0.0060.898 $±$0.006 0.808±0.0020.808 $±$0.002 0.821±0.0070.821 $±$0.007 0.832±0.0060.832 $±$0.006 0.843±0.0060.843 $±$0.006 0.861±0.0110.861 $±$0.011 0.860±0.0100.860 $±$0.010 SPECTF 0.862±0.0290.862 $±$0.029 0.839±0.0300.839 $±$0.030 0.894±0.0150.894 $±$0.015 0.845±0.0400.845 $±$0.040 0.857±0.0430.857 $±$0.043 0.894±0.0150.894 $±$0.015 0.845±0.0400.845 $±$0.040 Magic 0.885±0.0040.885 $±$0.004 0.850±0.0060.850 $±$0.006 0.857±0.0050.857 $±$0.005 0.855±0.0040.855 $±$0.004 0.863±0.0070.863 $±$0.007 0.864±0.0040.864 $±$0.004 0.868±0.0060.868 $±$0.006 Bank 0.908±0.0030.908 $±$0.003 0.901±0.0030.901 $±$0.003 0.902±0.0020.902 $±$0.002 0.903±0.0020.903 $±$0.002 0.903±0.0020.903 $±$0.002 0.905±0.0020.905 $±$0.002 0.907±0.0020.907 $±$0.002 Churn 0.958±0.0040.958 $±$0.004 0.885±0.0090.885 $±$0.009 0.886±0.0050.886 $±$0.005 0.889±0.0050.889 $±$0.005 0.954±0.0070.954 $±$0.007 0.946±0.0030.946 $±$0.003 0.953±0.0080.953 $±$0.008 Table 15: Accuracy Score for Classification Datasets, RF Dataset BlackBox GAM CALM RF NAM EBM Spline NAM EBM Spline Adult 0.843±0.0010.843 $±$0.001 0.851±0.0020.851 $±$0.002 0.870±0.0020.870 $±$0.002 0.856±0.0010.856 $±$0.001 0.853±0.0020.853 $±$0.002 0.871±0.0030.871 $±$0.003 0.859±0.0010.859 $±$0.001 COMPAS 0.662±0.0160.662 $±$0.016 0.682±0.0160.682 $±$0.016 0.681±0.0120.681 $±$0.012 0.685±0.0120.685 $±$0.012 0.682±0.0160.682 $±$0.016 0.683±0.0120.683 $±$0.012 0.685±0.0120.685 $±$0.012 HELOC 0.724±0.0140.724 $±$0.014 0.723±0.0110.723 $±$0.011 0.728±0.0140.728 $±$0.014 0.726±0.0100.726 $±$0.010 0.723±0.0140.723 $±$0.014 0.728±0.0140.728 $±$0.014 0.726±0.0100.726 $±$0.010 MIMIC2 0.888±0.0020.888 $±$0.002 0.886±0.0010.886 $±$0.001 0.886±0.0030.886 $±$0.003 0.886±0.0030.886 $±$0.003 0.886±0.0020.886 $±$0.002 0.887±0.0020.887 $±$0.002 0.886±0.0030.886 $±$0.003 Appendicitis 0.859±0.0640.859 $±$0.064 0.848±0.0560.848 $±$0.056 0.877±0.0770.877 $±$0.077 0.887±0.0640.887 $±$0.064 0.848±0.0560.848 $±$0.056 0.877±0.0770.877 $±$0.077 0.887±0.0640.887 $±$0.064 Phoneme 0.898±0.0050.898 $±$0.005 0.808±0.0020.808 $±$0.002 0.821±0.0070.821 $±$0.007 0.832±0.0060.832 $±$0.006 0.838±0.0050.838 $±$0.005 0.859±0.0050.859 $±$0.005 0.859±0.0090.859 $±$0.009 SPECTF 0.877±0.0270.877 $±$0.027 0.839±0.0300.839 $±$0.030 0.894±0.0150.894 $±$0.015 0.845±0.0400.845 $±$0.040 0.848±0.0260.848 $±$0.026 0.894±0.0180.894 $±$0.018 0.842±0.0480.842 $±$0.048 Magic 0.878±0.0040.878 $±$0.004 0.850±0.0060.850 $±$0.006 0.857±0.0050.857 $±$0.005 0.855±0.0040.855 $±$0.004 0.861±0.0040.861 $±$0.004 0.864±0.0040.864 $±$0.004 0.866±0.0050.866 $±$0.005 Bank 0.901±0.0040.901 $±$0.004 0.901±0.0030.901 $±$0.003 0.902±0.0020.902 $±$0.002 0.903±0.0020.903 $±$0.002 0.905±0.0030.905 $±$0.003 0.907±0.0020.907 $±$0.002 0.907±0.0020.907 $±$0.002 Churn 0.957±0.0050.957 $±$0.005 0.885±0.0090.885 $±$0.009 0.886±0.0050.886 $±$0.005 0.889±0.0050.889 $±$0.005 0.952±0.0050.952 $±$0.005 0.951±0.0070.951 $±$0.007 0.955±0.0060.955 $±$0.006