2.1 Feature embedding

We start from tabular input data $\boldsymbol{x}=(x_1,\ldots , x_q)^\top$ . In a first step, we tokenize this tabular input data to an input tensor $\boldsymbol{\phi }=\boldsymbol{\phi }(\boldsymbol{x})= [\phi _1,\ldots , \phi _q ] \in {\mathbb R}^{d \times q}$ ; this step is done in complete analogy to Gorishniy et al. (Reference Gorishniy, Rubachev, Khrulkov and Babenko2021), Brauer (Reference Brauer2024), and Richman et al. (Reference Richman, Scognamiglio and Wüthrich2025). For this first step, we select a fixed embedding dimension $d \in {\mathbb N}$ ; this is a hyperparameter selected by the modeler.

Categorical features.

For categorical feature components $x_j\in \mathcal{X}_j=\{1,\ldots , n_j\} \subset {\mathbb N}$ , having $n_j$ levels, consider a $d$ -dimensional entity embedding $\phi _j\,:\, \mathcal{X}_j \to {\mathbb R}^d$

(2.1) \begin{equation} x_j \,\mapsto \,\phi _j(x_j) = \sum _{k=1}^{n_j} W^{(0)}_{j,k} \,\unicode{x1D7D9}_{\{x_j=k\}}, \end{equation}

where $W^{(0)}_j=[W^{(0)}_{j,1}, \ldots , W^{(0)}_{j,n_j}]^\top \in {\mathbb R}^{n_j \times d}$ is an embedding matrix, which is a learnable parameter.

Continuous features.

For continuous feature components $x_j \in {\mathbb R}$ , we consider a $d$ -dimensional FNN embedding $\phi _j\,:\, {\mathbb R} \to {\mathbb R}^d$

(2.2) \begin{equation} x_j \,\mapsto \,\phi _j(x_j) = W^{(2)}_j \,\tanh \left (W^{(1)}_j x_j + b^{(1)}_j\right ) + b^{(2)}_j \end{equation}

with biases $b^{(1)}_j \in {\mathbb R}^{d'}$ and $b^{(2)}_j \in {\mathbb R}^d$ , weight matrices $W^{(1)}_j \in {\mathbb R}^{d' \times 1}$ and $W^{(2)}_j \in {\mathbb R}^{d \times d'}$ , for a fixed number of units $d'$ , and where the hyperbolic tangent function is applied element-wise.

Collecting all embedded components $(\phi _j)_{j=1}^q=(\phi _j(x_j))_{j=1}^q$ of the input $\boldsymbol{x}$ gives us the tensor

(2.3) \begin{equation} \boldsymbol{\phi }=\boldsymbol{\phi }(\boldsymbol{x})= [\phi _1,\ldots , \phi _q ]= [\phi _1(x_1),\ldots , \phi _q(x_q) ] \in {\mathbb R}^{d \times q}. \end{equation}

This tensor $\boldsymbol{\phi }$ involves the embedding weights $W^{(0)}_j \in {\mathbb R}^{n_j \times d}$ from the categorical features (of dimension $n_jd$ ), and the FNN biases and weights $(b^{(1)}_j, b^{(2)}_j , W^{(1)}_j, W^{(2)}_j)$ from the continuous feature embeddings (of dimension $2 d' + (d'+1)d$ ). All these biases and weights are learned from the available training data.

The categorical entity embeddings (2.1) and the continuous feature embeddings (2.2) do not necessarily live on the same scale. If a common scale is desired, this can be achieved either through ridge regularization or by applying separate normalization layers to the two types of embeddings. In our examples, the embeddings are not directly extracted, but they are further processed by interaction tokens introduced in the next section. These interaction tokens alleviate the potential scale discrepancies. Therefore, we do not regularize or normalize the embeddings.

2.2 Pairwise interaction layer

The tensor $\boldsymbol{\phi }=\boldsymbol{\phi }(\boldsymbol{x})$ , given in (2.3) and containing the tokens $(\phi _j(x_j))_{j=1}^q$ , is expressing the individual feature components of $\boldsymbol{x} \in \mathcal{X}$ before they have been interacting with each other. In the next step, we are going to let these tokens interact in a pairwise manner. This can be interpreted in a similar way as the query–key interaction concept in the attention layer of Vaswani et al. (Reference Vaswani, Shazeer, Parmar, Uszkoreit, Jones, Gomez, Kaiser and Polosukhin2017).

We recall the (centered version of the) hard sigmoid activation

(2.4) \begin{equation} \sigma _{\textrm {hard}} \,:\,{\mathbb R} \to [0,1], \qquad x\,\mapsto \,\sigma _{\textrm {hard}}(x) = \max \left (0, \,\min \left (1, \, \frac {1+x}{2}\right )\right ). \end{equation}

This is the centered version of the hard sigmoid activation; there is also a noncentered version with $(1+x)/2$ being replaced by $x$ in (2.4). The centered version is directly comparable to the sigmoid activation function and the step function activation, see Figure 1.

Figure 1

Activation functions: (centered) hard sigmoid, sigmoid, and step function.

We now let $\phi _j$ and $\phi _k$ interact, $j \le k$ . To support these interactions, we introduce interaction tokens $\boldsymbol{e}_{j,k} \in {\mathbb R}^{d_0}$ , for a given hyper-parameter $d_0 \in {\mathbb N}$ . The tokens consist of a set of learnable parameters that are optimized jointly with the rest of the model during network training. Before interacting, these tokens do not encode any information and are randomly initialized, like the other network parameters. This is analogous to the CLS token used in Devlin et al. (Reference Devlin, Chang, Lee and Toutanova2018) and Richman et al. (Reference Richman, Scognamiglio and Wüthrich2025). After training, the tokens are calibrated to support pairwise interactions between the $j$ -th and $k$ -th covariates. We concatenate the three vectors to a row-vector $(\phi _j, \phi _k, \boldsymbol{e}_{j,k})^\top \in {\mathbb R}^{2d+d_0}$ , $j \le k$ . The goal of the interaction tokens is to learn the interaction structure within the tensor $\boldsymbol{\phi }$ for predicting the response $Y$ , in particular, different pairs $(\phi _j, \phi _k)$ may have a different interaction behavior which is reflected in different learned interaction tokens $\boldsymbol{e}_{j,k}$ , $j \le k$ . Unlike the CLS token of Devlin et al. (Reference Devlin, Chang, Lee and Toutanova2018) and Richman et al. (Reference Richman, Scognamiglio and Wüthrich2025), we do not use these interaction tokens to strip out information, but we add them to allow for the necessary modeling flexibility to be able to account for differences in the different interaction pairs.

Define the shared interaction network $f_\theta \,:\,{\mathbb R}^{2d+d_0} \to {\mathbb R}$ as a composition of three FNN layers

(2.5) \begin{equation} f_\theta (\phi _j, \phi _k, \boldsymbol{e}_{j,k}) = W^{(3)} \,\textrm {ReLU}\big (W^{(2)} \,\textrm {ReLU}\big (W^{(1)} (\phi _j, \phi _k, \boldsymbol{e}_{j,k})^\top + b^{(1)}\big ) + b^{(2)}\big ) + b^{(3)}, \end{equation}

with:

• • the 1st FNN layer has $d_1$ units, rectified linear unit (ReLU) activation, network weights $W^{(1)}\in {\mathbb R}^{d_1 \times (2d+d_0)}$ , and bias $b^{(1)}\in {\mathbb R}^{d_1}$ ;

• • the 2nd FNN layer has $d_2$ units, ReLU activation, network weights $W^{(2)}\in {\mathbb R}^{d_2 \times d_1}$ , and bias $b^{(2)}\in {\mathbb R}^{d_2}$ ;

• • the output layer has $1$ unit, linear activation, network weights $W^{(3)}\in {\mathbb R}^{1 \times d_2}$ , and bias $b^{(3)}\in {\mathbb R}$ ;

and $\theta =(W^{(1)},b^{(1)}, W^{(2)},b^{(2)}, W^{(3)},b^{(3)})$ collects all network parameters. This interaction network $f_\theta$ is going to be shared by all pairwise interactions between the components of $\boldsymbol{\phi }$ , i.e., there is one single network parameter $\theta$ that is shared by all ordered pairs $(\phi _j,\phi _k)$ , $1\le j \le k \le q$ . This can also be seen as a so-called time-distributed network as being used, e.g., in recurrent neural networks. To differentiate different forms of interactions of the different ordered pairs, $j \le k$ , we let the interaction tokens $\boldsymbol{e}_{j,k}$ learn the pairwise interaction structures.

Definition 2.1 (pairwise interaction layer). For feature tokens $\phi _j({\cdot} )$ and $\phi _k({\cdot} )$ , $1\le j \le k \le q$ , and interaction tokens $\boldsymbol{e}_{j,k}$ , we define the pairwise interaction units by

(2.6) \begin{equation} h_{j,k}(\boldsymbol{x}) = \sigma _{\mathrm {hard}} (f_\theta (\phi _j(x_j), \phi _k(x_k), \boldsymbol{e}_{j,k} ) )\,\in \,[0,1]. \end{equation}

The pairwise interaction layer for $\boldsymbol{x} \in \mathcal{X}$ is defined by the upper-right triangular matrix

(2.7) \begin{equation} \left (h_{j,k}(\boldsymbol{x})\right )_{1\le j \le k \le q} = \begin{pmatrix} h_{1,1}(\boldsymbol{x}) & h_{1,2}(\boldsymbol{x}) & \cdots & h_{1,q-1}(\boldsymbol{x})& h_{1,q}(\boldsymbol{x})\\ \mathrm {n/a} & h_{2,2}(\boldsymbol{x}) & \cdots & h_{2,q-1}(\boldsymbol{x})& h_{2,q}(\boldsymbol{x})\\ \vdots &\vdots &\ddots &\vdots &\vdots &\\ \mathrm {n/a} & \mathrm {n/a} & \cdots & h_{q-1,q-1}(\boldsymbol{x})& h_{q-1,q}(\boldsymbol{x})\\ \mathrm {n/a} & \mathrm {n/a} & \cdots & \mathrm {n/a} & h_{q,q}(\boldsymbol{x}) \end{pmatrix} , \end{equation}

where the lower-left part is undefined indicated by n/a.

We remark that the diagonal entries represent the “self-interaction” of each covariate component and can be interpreted as the corresponding individual (main) effects. Including these diagonal terms allows the model to learn nonlinear functions that capture the marginal contribution of each variable, while preserving the same unified representation used for modeling pairwise interactions.

Example 2.2. We give a simple example of the pairwise interaction layer. We assume only continuous features, that is, $\boldsymbol{x} \in {\mathbb R}^q$ . We select one-dimensional linear feature embeddings (tokenization) $x_j \mapsto \phi _j(x_j)=x_j$ for all $1\le j \le q$ . This gives us the tensor $\boldsymbol{\phi } = \boldsymbol{\phi }(\boldsymbol{x})=\boldsymbol{x} \in {\mathbb R}^{1 \times q}$ . For the shared interaction network $f_\theta$ , we select a simple additive function and the interaction tokens are set off in this example. This gives the simplified version

(2.8) \begin{equation} h^{\mathrm {add}}_{j,k}(\boldsymbol{x}) = \sigma _{\mathrm {hard}}(x_j+x_k)\qquad \text{ for $j\le k$.} \end{equation}

This is illustrated on the left-hand side of Figure 2 for $j\lt k$ .

The right-hand side of Figure 2 shows a linear additive interaction network with the first component $x_j$ being scaled with $1/3$ . We note that these linear additive interaction networks are rather similar to binary decision tree partitions of the feature space. We allow for additional flexibility by letting the interaction network $f_\theta$ be nonlinear and by letting the interaction token $\boldsymbol{e}_{j,k}$ learn different pairwise interaction structures. It is precisely this ability to learn to modify the inputs that gives the pairwise interaction layers their power.

Figure 2

Interaction layers: (lhs) $\sigma _{\mathrm {hard}}(x_j+x_k)$ and (rhs) $\sigma _{\mathrm {hard}}(x_j/3+x_k)$ for $j\lt k$ .

2.3 Pairwise interaction network

We are now ready to define the tree-like PIN.

Definition 2.3. Choose an output activation function $g$ . The PIN) is defined by $f_{\mathrm {PIN}}\,:\,\mathcal{X} \to {\mathbb R}$

(2.9) \begin{equation} \boldsymbol{x} \,\mapsto \, f_{\mathrm {PIN}}(\boldsymbol{x}) = g\left (\sum _{1\le j \le k \le q} w_{j,k} \, h_{j,k}(\boldsymbol{x})+b\right )\!, \end{equation}

with learnable (real-valued) output weights $(w_{j,k})_{1\le j \le k \le q}$ and bias $b \in {\mathbb R}$ .

There is structural similarity between a PIN and a classical gradient boosting machine (GBM)

\begin{equation*} f_{\mathrm {GBM}}(\boldsymbol{x}) = \sum _{m=1}^M \sum _{l=1}^L \gamma _{m,l}\, \unicode {x1D7D9}_{R_{m,l}}(\boldsymbol{x}), \end{equation*}

with $M\in {\mathbb N}$ GBM iterations considering binary split decision trees with $L$ leaves giving the partition $(R_{m,l})_{l=1}^L$ of the feature space $\mathcal{X}$ and assigning the leaf values $(\gamma _{m,l})_{l=1}^L$ to this partition. We start by discussing the structural correspondence between the two regression functions.

We list some key differences:

1. 1. Smooth activation: The hard binary decisions $\unicode{x1D7D9}_{R_{m,l}}(\boldsymbol{x})$ are replaced by continuous functions involving the hard sigmoid function.

2. 2. Learned representations: There is an automatic feature transformations through embeddings resulting in the tensor $\boldsymbol{\phi }$ , whereas GBMs split on the original feature values.

3. 3. Parameter sharing: Common network $f_\theta$ across all interactions, and the interaction tokens learn the differences between the pairwise interactions $(\phi _j,\phi _k)$ .

4. 4. Joint optimization: All components trained simultaneously, where as the GBM learns the parameters recursively in stage-wise adaptive way.

The other close connection of the PIN is to a GA $^{2}$ M. A GA $^{2}$ M is obtained from a GAM by adding interaction terms to the regression function

(2.10) \begin{equation} g^{-1}\left (\mu (\boldsymbol{X})\right ) = b+\sum _{j=1}^q f_j(X_j) + \sum _{1\le j\lt k \le q} f_{j,k}(X_j,X_k), \end{equation}

where $f_j$ and $f_{j,k}$ are selected from a given class of splines; see Lou et al. (Reference Lou, Caruana, Gehrke and Hooker2013). A classical GAM only considers the terms $(f_j)_{j=1}^q$ , and the GA $^{2}$ M adds the pairwise interaction terms $(f_{j,k})_{1\le j\lt k \le q}$ . We note that this GA $^{2}$ M (2.10) has the same structure as our PIN (2.9), the difference is that we replace the splines $(f_j)_{j=1}^q$ and $(f_{j,k})_{1\le j\lt k \le q}$ by pairwise interaction layers $(h_{j,k})_{1\le j\le k \le q}$ being based on FNNs.

2.4 Connection to graph neural networks

The pairwise interaction structure of the PIN can be naturally interpreted through the lens of GNNs. We formalize this connection and demonstrate how our model can be viewed as a special case of message passing on a complete graph with learned edge features.

Our PIN can be reformulated as a message passing neural network with a single update step.

Definition 2.5 (message function). The message function $\mathcal{M}\,:\,{\mathbb R}^d \times {\mathbb R}^d \times {\mathbb R}^{d_0} \to {\mathbb R}$ is defined as

\begin{equation*} \mathcal{M}(\phi _j, \phi _k, \boldsymbol{e}_{j,k}) = \sigma _{\mathrm {hard}}(f_\theta (\phi _j, \phi _k, \boldsymbol{e}_{j,k})), \end{equation*}

where $f_\theta$ is the shared interaction network (2.5) .

It is obvious that the PIN (2.9) is equivalent to the GNN defined by

\begin{equation*} f_{\mathrm {GNN}}(\boldsymbol{x}) = g\left (\sum _{(j,k) \in E} w_{j,k}\, \mathcal{M}(\phi _j, \phi _k, \boldsymbol{e}_{j,k})\right )\!. \end{equation*}

Our architecture differs from traditional GNNs in the following ways:

1. 1. Edge-centric: While most GNNs focus on learning node representations, our model primarily focuses on edge interactions and the interaction token $\boldsymbol{e}_{j,k}$ to distinguish different interaction patterns.

2. 2. Learned edge features: The interaction tokens $\boldsymbol{e}_{j,k}$ are learned parameters rather than derived features.

3. 3. Single-pass architecture: Unlike typical GNNs that use multiple layers of message passing, our model uses a single pass with a more expressive message function.

4. 4. Bounded messages: The hard sigmoid activation ensures bounded message values, similar to the binary nature of GBM binary decision tree splits.

Having this connection to GNNs provides the following additional theoretical insights.

This proposition immediately follows by analyzing permutations of the components/vertices $V=\{1,\ldots , q\}$ .

Using GNN operations, the GNN perspective suggests several natural extensions to our base model.

Definition 2.8 (multi-layer PIN). We can define a $K$ -layer model through iterative message passing, that is,

(2.11) \begin{equation} \phi _v^{(k+1)} = \operatorname {UPDATE}\left (\phi _v^{(k)}, \sum _{u \in \mathcal{N}(v)} \mathcal{M}^{(k)}\big(\phi _u^{(k)}, \phi _v^{(k)}, \boldsymbol{e}_{u,v}\big)\right ) \end{equation}

where $\phi _v^{(0)} = \phi _v(x_v)$ , $v \in V$ , $k = 0,\ldots ,K-1$ , and $\mathcal{N}(v)$ are the vertices adjacent to $v$ . The $\operatorname {UPDATE}$ can be various functions such as the maximum or the mean of the two components.

Remark 2.9 (attention mechanism). The model can be enhanced with an attention mechanism by replacing the fixed output weights $w_{j,k}$ with learned functions of the node features

\begin{equation*} w_{j,k} = \alpha _\vartheta (\phi _j, \phi _k), \end{equation*}

where $\alpha _\vartheta$ is a (self-)attention layer with parameters $\vartheta$ ; see Vaswani et al. (Reference Vaswani, Shazeer, Parmar, Uszkoreit, Jones, Gomez, Kaiser and Polosukhin2017). If we apply the attention weights only to pairs $j\lt k$ , and if we apply a thresholding, we obtain

(2.12) \begin{equation} \boldsymbol{x} \,\mapsto \, g\left (b+\sum _{j=1}^k w_{j} \, h_{j,j}(\boldsymbol{x})+\sum _{1\le j \lt k \le q} \mathrm {ReLU} (\alpha _\vartheta (\phi _j, \phi _k)-\tau )\, h_{j,k}(\boldsymbol{x})\right ), \end{equation}

for an attention threshold $\tau \in {\mathbb R}$ . The thresholding in (2.12) gives a natural way of pairwise interaction selection.

The connection to GNNs not only provides theoretical insights but also suggests practical extensions and implementation strategies for the PIN architecture. Moreover, for efficient implementations, existing GNN libraries are available that allow for batched processing (graph batching for parallel processing), optimized GNN kernels on GPUs, and optional sparsification of the interaction graph, similar to (2.12).

4.1 PIN interaction importance

A notoriously difficult problem in machine learning models is explainability. Scholars have developed various post hoc measures of variable importance, e.g., permutation variable importance by Breiman (Reference Breiman2001) and Fisher et al. (Reference Fisher, Rudin and Dominici2018), or SHAP importance of Lundberg and Lee (Reference Lundberg and Lee2017) and Mayer and Stando (Reference Mayer and Stando2025) are two popular examples. When it comes to interactions, things become even more difficult. There is Friedman and Popescu (Reference Friedman and Popescu2008)’s $H^2$ -statistics which is a post hoc explainability tool that measures the impact of bivariate components on the regression function (by partial dependence functions). Apart from that, there is little guidance on how to measure the interaction strengths at a reasonable computational effort, i.e., testing pairwise interactions can be very costly.

Coming back to the PIN (2.9), we modify this function in a way so that we do not aim at getting the maximal predictive power, but we try to evaluate which shared interaction unit $h_{j,k}(\boldsymbol{x})$ contributes most to a loss decrease. For this, we consider the multi-output mapping

(4.1) \begin{equation} \boldsymbol{x} \,\mapsto \, ( f_{\mathrm {PIN},j,k} )_{1\le j \lt k \le q} \in \, {\mathbb R}^{q(q-1)/2}, \end{equation}

where for $1\le j\lt k \le q$ we set

\begin{equation*} f_{\mathrm {PIN},j,k} = g\left (b+\sum _{l=1}^q w_{l,l} \, h_{l,l}(\boldsymbol{x})+ b_{j,k}+ w_{j,k} \, h_{j,k}(\boldsymbol{x})\right )\!, \end{equation*}

for biases $b_{j,k}\in {\mathbb R}$ and the remaining terms have the same structure as in (2.9). The first part of this multi-output regression function captures the main effects of the covariate components, this corresponds to the diagonal of the tensor (2.7). That is, the main effects are given by

(4.2) \begin{equation} b+\sum _{l=1}^q f_l(X_l) =b+\sum _{l=1}^q w_{l,l} \, h_{l,l}(\boldsymbol{x}), \end{equation}

where we make the link to the GAM (2.10). The remaining terms $h_{j,k}(\boldsymbol{x})$ in (4.1), $j\lt k$ , correspond to the off-diagonal terms of the pairwise interaction layer (2.7). Each component in (4.1) adds another interaction term to the diagonal (non-interacting) terms (4.2), and we can determine the term $h_{j,k}(\boldsymbol{x})$ , $j\lt k$ , that leads to the biggest decrease in out-of-sample loss when added to (4.2). The great thing is that we can fit all $q(q-1)/2$ interaction terms in the same gradient descent run, i.e., we apply the objective function for model fitting (simultaneously) to all components of (4.1), and jointly fit all components by just aggregating these component-wise losses

\begin{eqnarray*} L_{\mathrm {AGGR}}(\mathcal{L}, (f_{\mathrm {PIN},j,k})_{j\lt k}) &=&\frac {1}{q(q-1)/2}\sum _{1\le j \lt k \le q} L_{\mathrm {PD}}(\mathcal{L}, f_{\mathrm {PIN},j,k}), \end{eqnarray*}

the latter is the average PD loss (3.1). This is possible here because we have isolated all interaction pairs in the pairwise interaction layer (2.7).

Figure 7

(lhs) Friedman and Popescu (Reference Friedman and Popescu2008)’s $H^2$ -statistics, and (rhs) PIN interaction importance.

In our example, we consider one further simplification, namely, we replace the diagonal (non-interacting) terms (4.2) by the fitted GAM terms that are considered to be optimal under the exclusion of interactions, see Table 2. This fitted GAM is kept fixed (frozen) and only the off-diagonal interaction terms in (4.1) are trained with SGD. The training of this multi-output PIN takes approximately 15 minutes on an ordinary laptop.Footnote 2 We then compute the decrease in out-of-sample loss compared to the GAM for each individual interaction term, $j\lt k$ . Figure 7 (rhs) shows the results (ranked by size). The biggest decrease in loss compared to the GAM can be achieved by adding an interaction term for (BonusMalus:VehBrand), and the second biggest decrease is obtained by adding (BonusMalus:Region).

Figure 7 (lhs) illustrates the Friedman and Popescu (Reference Friedman and Popescu2008) $H^2$ -statistics computed for the full PIN model. In our experiments, this procedure requires around 90 minutes (typically one chooses a small background dataset to speed up computations). Not surprisingly, the same feature interactions seem to be the most impactful ones, but the order of the most important interactions is different. Friedman–Popescu’s $H^2$ -statistics studies the partial dependence function in a fully fitted model (which can be seen as a top-down view – backward selection), whereas our PIN approach considers which interaction should be added to have the biggest impact on the prediction accuracy (which is more of a bottom-up consideration – forward selection). Thus, we can interpret this as a forward selection of interaction terms, similar to the forward selection model building in GLMs. Moreover, Friedman–Popescu’s $H^2$ -statistics is based on partial dependence plots (PDPs), which have the weakness that they are not necessarily reliable in cases where the feature components are highly correlated, because PDPs extrapolate the regression function to rare or impossible covariate combinations (e.g., car drivers aged 20 with 10 years of driving experience), where the regression function is typically poorly estimated or undefined, respectively. Not surprisingly, this weakness also translates to the $H^2$ -statistics. Our method does not suffer from this weakness, since the contribution of each pairwise interaction is evaluated on the empirically observed portfolio – on the empirical population distribution – and assessing each interaction independently of all others. As a consequence, the importance measure of each pairwise interaction remains meaningful even in the presence of feature correlations. Naturally, as is also common in the GLM forward selection process, if there are two highly correlated feature components that interact with a third component, after including the first strong interaction, the second one may no longer be needed. We discuss this next.

Figure 8

(lhs) PIN interaction importance – 1st round (this plot is identical to Figure 7, rhs), and (rhs) PIN interaction importance – 2nd round.

Similar to the method of forward selection in GLMs, we add the most significant interaction term (BonusMalus:VehBrand) to the diagonal regression model (4.2). We also freeze this interaction term, and we repeat the same fitting strategy to find the second most important interaction term, after (BonusMalus:VehBrand) has been added to the model. The result of this second round of forward selection is shown in Figure 8 (rhs), and it is compared to the first round results on the left-hand side. The second interaction that should be added is (DrivAge:BonusMalus). This gives a systematic way of model building that may also be useful in improving GLMs and GAMs; this is similar to the forward selection algorithm of GLM model building. Note that this forward selected model should be refitted simultaneously on all included terms, because in this forward inclusion process, we always freeze the previously selected terms, and an additional re-fitting (fine-tuning) step may find a more optimal jointly learned model on the selected pairwise interaction components.

4.2 SHapley’s Additive exPlanation

There is another very remarkable property of the tree-like PIN architecture, namely, it allows for a very efficient computation of SHapley’s Additive exPlanation (SHAP) as introduced by Lundberg and Lee (Reference Lundberg and Lee2017). Shapley values originate from cooperative game theory. They are used to share a total payoff of a cooperative game in a fair way among $q$ players; see Shapley (Reference Shapley1953). Denote by $\mathcal{C} \subseteq \mathcal{Q}=\{1,\ldots , q\}$ any coalition of $q$ players. Based on Shapley’s fairness axioms, there is precisely one way to allocate the total payoff to the $q$ players of the cooperative game, resulting in the Shapley values $\psi _j$ assigned to players $j \in \mathcal{Q}$ , given by

(4.3) \begin{equation} \psi _j = \frac {1}{q} \sum _{\mathcal{C} \subseteq \mathcal{Q} \setminus \{j\}} {\begin{pmatrix}q-1\\|\mathcal{C}|\end{pmatrix}}^{-1} (\nu (\mathcal{C} \cup \{j\}) - \nu (\mathcal{C}) ), \qquad \text{ for all }j \in {\mathcal{Q}}, \end{equation}

for a given value function $\mathcal{C}\mapsto \nu (\mathcal{C})$ . In fact, starting from the given value function $\nu$ , the Shapley values $(\psi _j)_{j=1}^q$ are the unique solution to a set of four axioms called Shapley’s fairness axioms; see Shapley (Reference Shapley1953) and Lundberg and Lee (Reference Lundberg and Lee2017).

There are two equivalent reformulations of the Shapley values (4.3), namely, there is the kernel SHAP formulation of Lundberg and Lee (Reference Lundberg and Lee2017, Theorem 2), and there is the permutation SHAP version of Štrumbelj and Kononenko (Reference Štrumbelj and Kononenko2010, Reference Štrumbelj and Kononenko2014). In this work, we focus on the permutation SHAP formulation. Denote by $\pi =(\pi _1,\ldots , \pi _q )$ a permutation of the ordered set $(1,\ldots , q)$ . Let $\kappa (j) \in \mathcal{Q}$ be the index defined by $\pi _{\kappa (j)}=j$ and set

(4.4) \begin{equation} \mathcal{C}_{\pi , j} = \{ \pi _1,\ldots , \pi _{\kappa (j)-1} \}\,\subset \,\mathcal{Q}. \end{equation}

These are all components of $\pi$ preceding $\pi _{\kappa (j)}=j$ ; for $\kappa (j)=1$ , this is the empty set $\emptyset$ . The permutation SHAP version of Štrumbelj and Kononenko (Reference Štrumbelj and Kononenko2010, Reference Štrumbelj and Kononenko2014) proves that the Shapley values (4.3) can equivalently be computed by

(4.5) \begin{equation} \psi _j = \frac {1}{q!} \sum _{\pi }\nu (\mathcal{C}_{\pi , j}\cup \{j\} ) -\nu (\mathcal{C}_{\pi , j} ) , \end{equation}

where the summation runs over all permutations $\pi$ . These are $q!$ permutations $\pi$ , which can be an extremely large number making Shapley values computationally unattainable for large $q$ . However, there is a nice result which says that if the value function $\nu$ only includes interactions of maximal order two, then an arbitrary permutation $\pi =(\pi _1,\ldots , \pi _q)$ , and its reverted version $\rho (\pi )=(\pi _q,\ldots , \pi _1)$ is sufficient to compute the Shapley values

(4.6) \begin{equation} \psi _j = \frac {1}{2} (\nu (\mathcal{C}_{\pi , j}\cup \{j\} ) -\nu (\mathcal{C}_{\pi , j} ) +\nu (\mathcal{C}_{\rho (\pi ), j}\cup \{j\} ) -\nu (\mathcal{C}_{\rho (\pi ), j} ) ), \end{equation}

if the value function $\nu$ is a bilinear form; see Mayer and Wüthrich (Reference Mayer and Wüthrich2025, Proposition 3.1). For example, the permutation $\pi =(1,\ldots , q)$ and its reverted version $\rho (\pi )=(q,\ldots , 1)$ do the job.

For explaining machine learning predictions, one typically works on the link scale, i.e., one drops the link function $g$ in (2.9). For a given feature $\boldsymbol{x}$ , this motivates the total value of the grand coalition $\mathcal{Q}$

\begin{equation*} \nu (\mathcal{Q}):=\nu _{\boldsymbol{x}}(\mathcal{Q}):= g^{-1} (f_{\mathrm {PIN}}(\boldsymbol{x}) )= \sum _{1\le j \le k \le q} w_{j,k} \, h_{j,k}(\boldsymbol{x})+b. \end{equation*}

This is the prediction made by $\boldsymbol{x}=(x_1,\ldots , x_q)^\top$ , and we want to explain how the feature components $x_j$ contribute to this prediction. This requires the choice of the value function $\nu (\mathcal{C})=\nu _{\boldsymbol{x}}(\mathcal{C})$ for any coalition $\mathcal{C}\subset \mathcal{Q}$ if the elements in $\mathcal{Q}\setminus \mathcal{C}$ are masked (not available for prediction). The most popular choice for describing the masking is the so-called interventional SHAP version, in our case given by

(4.7) \begin{equation} \nu (\mathcal{C}):=\nu _{\boldsymbol{x}}(\mathcal{C}) := {\mathbb E}\left [\sum _{1\le j \le k \le q} w_{j,k} \, h_{j,k} (\boldsymbol{x}_\mathcal{C}, \boldsymbol{X}_{\mathcal{Q}\setminus \mathcal{C}} )+b\right ], \end{equation}

where the components $\boldsymbol{x}_\mathcal{C}=(x_j)_{j \in \mathcal{C}}$ are kept fixed and the components $\boldsymbol{X}_{\mathcal{Q}\setminus \mathcal{C}}=(X_j)_{j \in \mathcal{Q}\setminus \mathcal{C}}$ are resampled from the population distribution $\boldsymbol{X} \sim {\mathbb P}$ . An important point in the subsequent computation of the Shapley values is that each term $h_{j,k}$ in (4.7) only depends on the components $(j,k)$ , i.e., we only have pairwise interactions, see (2.6). That is, by a slight abuse of notation,

\begin{equation*} h_{j,k} (\boldsymbol{x}_\mathcal{C}, \boldsymbol{X}_{\mathcal{Q}\setminus \mathcal{C}} ) = \left \{ \begin{array}{ll} h_{j,k}(x_j, x_k) & \text{ if $j,k \in \mathcal{C}$,}\\[2pt] h_{j,k}(X_j, x_k) & \text{ if $j \in \mathcal{Q} \setminus \mathcal{C}$ and $k \in \mathcal{C}$,}\\[2pt] h_{j,k}(x_j, X_k) & \text{ if $j \in \mathcal{C}$ and $k\in \mathcal{Q} \setminus \mathcal{C}$,}\\[2pt] h_{j,k}(X_j, X_k) & \text{ if $j,k \in \mathcal{Q} \setminus \mathcal{C}$.}\\ \end{array} \right . \end{equation*}

This interventional SHAP formula (4.7) is usually estimated empirically. For this, we select a background set $\mathcal{B}$ from the learning sample $\mathcal{L}$ . This allows one to determine the value function empirically by

(4.8) \begin{equation} \widehat {\nu }(\mathcal{C})=\widehat {\nu }_{\boldsymbol{x}}(\mathcal{C}) = \frac {1}{|\mathcal{B}|}\sum _{i \in \mathcal{B}} \left [\sum _{1\le j \le k \le q} w_{j,k} \, h_{j,k}\left (\boldsymbol{x}_\mathcal{C}, \boldsymbol{X}^{(i)}_{\mathcal{Q}\setminus \mathcal{C}}\right )+b\right ], \end{equation}

where $\boldsymbol{X}^{(i)}$ is the $i$ -th instance of the background dataset.

The crucial point now is that the value function $\mathcal{C} \mapsto \widehat {\nu }(\mathcal{C})=\widehat {\nu }_{\boldsymbol{x}}(\mathcal{C})$ only involves interactions of maximal order two. Therefore, the paired permutation SHAP formula (4.6) applies (and it is exact for the selected value function $\widehat {\nu }=\widehat {\nu }_{\boldsymbol{x}}$ ). This requires $2(q+1)$ evaluations of the value function  $\widehat {\nu }$ , and each value function evaluation requires $|\mathcal{B}|$ computations of PIN predictions, see (4.8). Thus, altogether, the SHAP decomposition of the prediction $g^{-1}(f_{\mathrm {PIN}}(\boldsymbol{x}) )$ requires $2(q+1)|\mathcal{B}|$ PIN network computations for computing the exact Shapley values $(\psi _j)_{j=1}^q$ of the value function $\widehat {\nu }$ . Any other network architecture that involves higher-order interactions would need $q!\, |\mathcal{B}|$ network computations for receiving the exact Shapley values of the corresponding estimated interventional SHAP value function; thus, the order of computations grows as $q!$ instead of $2(q+1)$ as in our tree-like PIN case.

We revisit the PIN architecture depicted in Table 1 having the best predictive performance of all models presented in Table 2. We emphasize that this network architecture only considers interactions of order two on the link scale; thus, the paired-sampling permutation SHAP formula (4.6) applies, and it gives exact results, for instance, we can use the permutation pair $\pi =(1,\ldots , q)$ and $\rho (\pi )=(q,\ldots , 1)$ . As background dataset, we select at random $|\mathcal{B}|=2,000$ instances from the learning set. This allows us to decompose the predictions $g^{-1}(f_{\mathrm {PIN}}(\boldsymbol{x}) )$ of selected instances $\boldsymbol{x}$ . Decomposing 100 predictions takes 23 s on an ordinary laptop (see footnote on page 36).

Figure 9

(lhs) SHAP waterfall graph $(\psi _j)_{j=1}^q$ of a single instance $\boldsymbol{x}$ and (rhs) SHAP variable importance.

The left-hand side of Figure 9 shows the Shapley values $(\psi _j)_{j=1}^q$ of a single instance $\boldsymbol{x}$ . We have global mean of the background dataset of $\widehat {\nu }_{\boldsymbol{x}}(\emptyset )=-2,776$ , and the prediction of feature $\boldsymbol{x}$ is given by $\widehat {\nu }_{\boldsymbol{x}}(\mathcal{Q})=g^{-1}(f_{\mathrm {PIN}}(\boldsymbol{x}) )=-3.239$ . Thus, feature $\boldsymbol{x}$ has a prediction being $-0.463$ below average. The Shapley values $(\psi _j)_{j=1}^q$ illustrated in Figure 9 (lhs) show how this lower prediction can be explained. The biggest decrease is explained by the bonus-malus level and the density features, whereas the vehicle gas feature leads to the biggest increase in the explanation of this difference from $\widehat {\nu }_{\boldsymbol{x}}(\emptyset )$ to $\widehat {\nu }_{\boldsymbol{x}}(\mathcal{Q})$ .

Going from local to global SHAP explanations, we can decompose the predictions $g^{-1}(f_{\mathrm {PIN}}(\boldsymbol{x}_i))$ of many instances $(\boldsymbol{x}_i)_{i \in \mathcal{I}}$ , for a randomly selected subset $\mathcal{I}$ of the entire data. This provides us with SHAP decompositions $(\psi ^{(i)}_j)_{j=1}^q$ , which allow us to consider different statistics. The following analysis is based on $|\mathcal{I}|=1,000$ instances $\boldsymbol{x}_i$ that have been decomposed by the paired-sampling permutation SHAP method. SHAP variable importance is then obtained by considering

\begin{equation*} \overline {\psi }_j = \frac {1}{|\mathcal{I}|} \sum _{i \in \mathcal{I}} |\psi ^{(i)}_j|. \end{equation*}

Figure 9 (rhs) shows these importance measures $\overline {\psi }_j$ for all feature components $j \in \mathcal{Q}$ . The most significant role is played by the bonus-malus level (past claims history) followed by the driver’s age variable and at the other end is the area code that seems to be the least important variable.

Figure 10

(lhs) Dependence between driver’s age and the bonus-malus variable, this graph is taken from Wüthrich and Merz (Reference Wüthrich and Merz2023, Figure 13.12), and (rhs) Shapley values $\psi _j^{(i)}$ of the selected 100 instances $i\in \mathcal{I}$ for the covariate $j=\texttt {DrivAge}$ ; the coloring is the bonus-malus level of these instances.

Figure 10 Long description

Panel A: A box plot shows the relationship between driver age (DrivAge) and bonus-malus levels (BonusMalus). The x-axis represents driver age in years, grouped into categories, while the y-axis represents the bonus-malus levels ranging from 50 to 120. The box plot displays the distribution of bonus-malus levels for each age group, with visible trends showing higher bonus-malus levels for younger drivers and a decrease in levels as driver age increases. Panel B: A scatter plot illustrates the Shapley values for the covariate driver age (DrivAge) against the driver age. The x-axis represents driver age in years, and the y-axis represents Shapley values ranging from -0.4 to 0.2. The points are colored according to the bonus-malus levels, with a color gradient from red to green indicating lower to higher bonus-malus levels. The plot shows a pattern where Shapley values vary with driver age, with notable clustering and trends visible across different age groups.

In Figure 10 (rhs), we analyze the 1,000 decomposed instances $\boldsymbol{x}_i$ , $i\in \mathcal{I}$ . The plot shows the Shapley values $\psi ^{(i)}_j$ of these 1,000 instances for the variable $j =\texttt {DrivAge}$ . This driver age variable has its strongest interaction with the bonus-malus level, see Figure 8, and the coloring of the Shapley values in Figure 10 (rhs) is selected according to these bonus-malus levels. We observe that the driver’s age variable has a decreasing effect on the prediction for drivers’ ages below 40 (negative Shapley values) and an increasing effect for ages between 40 and 60 (positive Shapley values). To fully understand this picture, one needs to understand the interaction with the bonus-malus level, and we observe for drivers’ ages below 40, that this positive age effect especially holds for low bonus-malus levels. In fact, there is a strong dependence between the driver’s age variable and the bonus-malus level for drivers’ ages below 40, see boxplot of Figure 10 (lhs), because new drivers start at a bonus-malus level of 100 (yellow color in Figure 10, rhs). This level decreases in every year that the driver is free of accidents, and otherwise it increases. Therefore, young drivers typically have a higher bonus-malus level than mature drivers (in a low-frequency environment), and the Shapley values in Figure 10 (rhs) tell us that young drivers with a lower bonus-malus level (green colors) get a decrease compared to their older peers in the same bonus-malus class. Of course, this interaction makes perfect sense. This completes the French MTPL example.