Error-controlled non-additive interaction discovery in machine learning models

Consider a prediction task where we have $n$ independent and identically distributed (i.i.d.) samples $\mathbf{X}=\left\{x_{i}\right\}_{i=1}^{n}\in\mathbb{R}^{n\times p}$ and $\mathbf{Y}=\left\{y_{i}\right\}_{i=1}^{n}\in\mathbb{R}^{n\times 1}$, denoting the data matrix with $p$-dimensional features and the corresponding response, respectively. The response of the prediction task can be categorical labels such as disease status or numerical measurements such as body mass index, and the features could include gene expression data, microbial taxa abundance profiles, and more. The prediction task can be described using an ML model $f:\mathbb{R}^{p}\mapsto\mathbb{R}$ that maps from the input $x\in\mathbb{R}^{p}$ to the response $y\in\mathbb{R}$. When modeling the task, the function $f$ learns the dependence structure of $\mathbf{Y}$ on $\mathbf{X}$, enabling effective prediction with the fitted ML model.

In this study, we focus on interpreting the fitted ML model for data-driven scientific discovery by identifying the non-additive feature interactions that contribute most to the model’s predictions. We say that $\mathcal{I}\subset\{1,\cdots,p\}$ is a non-additive interaction of function $f$ if and only if $f$ cannot be decomposed into an addition of $\left|\mathcal{I}\right|$ subfunctions $f_{i}$, each of which excludes a corresponding interaction feature [49, 52], i.e., $f(x)\neq\sum_{i\in\mathcal{I}}f_{i}\left(x_{\{1,\cdots,p\}\backslash i}\right)$. For example, the multiplication between two features $x_{i}$ and $x_{j}$ is a non-additive interaction because it cannot be decomposed into a sum of univariate functions, i.e., $x_{i}x_{j}\neq f_{i}(x_{j})+f_{j}(x_{i})$. Assume that there exists a group of interactions $\mathcal{S}=\left\{\mathcal{I}_{1},\mathcal{I}_{2},\cdots\right\}$ such that conditional on interactions $\mathcal{S}$, the response $\mathbf{Y}$ is independent of interactions in the complement $\mathcal{S}^{c}=\{1,\cdots,p\}\times\{1,\cdots,p\}\backslash\mathcal{S}$. In this setting, our goal is to accurately discover feature interactions in $\mathcal{S}$ without erroneously reporting too many incorrect interactions in $\mathcal{S}^{c}$.

Diamond is designed to discover feature interactions from fitted ML models while maintaining a controlled FDR [4]. For a set of feature interactions $\widehat{S}\subset\{1,\cdots,p\}\times\{1,\cdots,p\}$ discovered by some interaction detection method, the FDR is defined as:

Commonly used procedures, such as the Benjamini–Hochberg procedure [4], achieve FDR control by working with p-values computed against some null hypothesis. In the interaction discovery setting, for each feature interaction, one tests the significance of the statistical association between the specific interaction and the response, either jointly or marginally, and obtains a p-value under the null hypothesis that the interaction is irrelevant. These p-values are then used to rank the features for FDR control. However, controlling FDR in generic machine learning models, especially deep learning models, is challenging because, to our knowledge, the field lacks a method for producing meaningful p-values that reflect interaction importance. To bypass the use of p-values but still achieve FDR control, we draw inspiration from the model-X knockoffs framework [2, 6]. The core idea of knockoffs is to generate dummy features that perfectly mimic the empirical dependence structure among the original features but are conditionally independent of the response given the original features. These knockoff features can then be used as a control by comparing the feature importance between the original features and their knockoff counterparts to estimate FDR (Fig. 1a).

Diamond trains a generic ML model that takes as input an augmented data matrix $(\mathbf{X},\tilde{\mathbf{X}})\in\mathbb{R}^{n\times 2p}$, created by concatenating the original data matrix $\mathbf{X}\in\mathbb{R}^{n\times p}$ with its knockoffs $\tilde{\mathbf{X}}\in\mathbb{R}^{n\times p}$. After training, Diamond quantifies feature interactions by interpreting the trained model and produces a ranked list of these interactions. Diamond subsequently uses a calibration procedure to distill non-additive or interactive effects from the reported interaction importance measures from existing methods, thereby maintaining FDR control at the target level (Fig. 1b). The calibration process is designed to remove both prediction-dependent marginal effects and prediction-independent feature biases from the reported feature interactions, leaving only the prediction-dependent non-additive interactive effects (Fig. 1c).

We started by evaluating the performance of Diamond on simulated datasets, assessing its ability to identify important non-additive interactions while controlling the FDR. We benchmarked Diamond on a test suite of 10 simulated datasets generated by different simulation functions proposed by [52]. These datasets contain a mixture of univariate functions and multivariate interactions, exhibiting varied order, strength, and nonlinearity (Fig. 2a). Since our goal is to detect pairwise interactions, high-order interaction functions (e.g., $F(x_{1},x_{2},x_{3})=x_{1}x_{2}x_{3}$) are decomposed into pairwise interactions (e.g., $(x_{1},x_{2})$, $(x_{1},x_{3})$, and $(x_{2},x_{3})$) to serve as the ground truth.

Following the settings used in [52], we employed a sample size of $n=20,000$, equally divided into training and test sets. In addition, the number of features is set at $p=30$, and all features are sampled randomly from a continuous uniform distribution, $U(0,1)$. Only the first $10$ out of $30$ features contribute to the corresponding response, while the remaining features serve as noise to increase the task’s complexity. For robustness, we repeated the experiment 20 times for each simulated dataset using different random seeds. Each repetition involved data generation, knockoff generation using KnockoffsDiagnostics [5], ML model training, and interaction-wise FDR estimation. For all simulation settings, we reported the mean performance with 95% confidence intervals, fixing the target FDR level at $q=0.2$.

Our analysis shows that Diamond consistently identifies important non-additive interactions with controlled FDR across various ML models (Fig. 2b). Under controlled FDR, the multi-layer perceptron (MLP) model demonstrates better performance in identifying important interactions than other ML models, as measured by the statistical power and the area under the receiver operating characteristic curve (AUROC). This superior performance is attributed to the inclusion of a knockoff-tailored plugin pairwise-coupling layer design [33], which has been shown to maximize statistical power (See details in Sec. S.1). We discovered that the proposed calibration procedure (Sec 4.2) is critical; without it, the FDR cannot be controlled by naively using reported interaction importance values from existing methods (Fig. 2b and Fig. S.2).

To gain insight into the FDR control failure in the absence of calibration, we conducted a qualitative comparison assessing interaction importance before and after calibration using the simulation function $F_{1}$ (Fig. 2c and Fig. S.3). The primary cause of the FDR control failure appears to lie in the distributional disparity between original-only interactions (i.e., interactions involving two features from the original feature set rather than knockoffs) and interactions involving knockoffs. This observation suggests a violation of the knockoff filter’s assumption in controlling the FDR (See discussion in Sec 4.2). The proposed calibration procedure mitigates the disparity by extracting non-additive interactive effects from the reported interaction importance measures, thereby enhancing the utility of knockoff-involving interactions as a negative control for FDR estimation.

Finally, we verified whether alternative baseline methods can accurately identify important non-additive interactions with controlled FDR. We compared two baseline methods for FDR estimation: one based on permutation-based interaction-wise p-values coupled with the Benjamini-Hochberg procedure, and the other representing interaction-wise FDR as the aggregation of feature-wise FDR (See details in Sec. 4.6). Our analysis shows that neither method correctly controls the FDR. This greatly reduces the utility of these methods, despite their reported high power and AUROC.

The design of Diamond incorporates several key components, including knockoff generation and interaction importance measures. To demonstrate the robustness of Diamond in FDR control, we conducted a control study where we modified Diamond by replacing these two components with alternative solutions. Specifically, we considered two variants of Diamond using the MLP model. In the first setting, we replaced KnockoffsDiagnostics with three alternatives: KnockoffGAN [23], Deep knockoffs [42] or VAE knockoffs [31], while keeping the MLP model and the interaction importance measure (i.e., Expected Hessian [14]) unchanged. In the second setting, we replaced the interaction importance measure, Expected Hessian, with two alternatives: Integrated Hessian [22] and a model-specific interaction importance measure derived from the MLP weights (See details in Sec. S.2), while keeping the MLP model and the knockoff generation unchanged. For each variant, we applied Diamond to the test suite of 10 simulated datasets using the same settings.

The results indicate that Diamond is robust to various knockoff designs and importance measures. In the first study, our analysis shows that Diamond consistently discovers important non-additive interactions with controlled FDR across various knockoff designs (Fig. 3a). Despite the methodological differences in knockoff generation between KnockoffsDiagnostics, KnockoffGAN, Deep knockoffs, and VAE knockoffs, Diamond achieves comparable statistical power and AUROC in each case. This result demonstrates the Diamond’s stability in identifying important interactions in knockoff-based FDR estimation.

In the second study, our analysis shows that Diamond discovers important non-additive interactions with controlled FDR across various interaction importance measures (Fig. 3b). Despite the methodological differences between Expected Hessian, Integrated Hessian, and model-specific measures, which involve mechanistic disparities in elucidating the relationships between features and responses, Diamond achieves valid FDR control while maintaining reasonably good statistical power and AUROC. This result demonstrates Diamond’s flexibility in identifying important interactions based on interaction-wise importance ranking.

Finally, it is worth mentioning that regardless of the knockoff generation and interaction importance measures we used, there is room for improvement in achieving better statistical power while maintaining controlled FDR. Diamond tends to be conservative by overestimating the FDR, suggesting that an improved FDR control procedure could potentially boost statistical power.

We then evaluated the performance of Diamond on a real dataset [13], assessing its ability to identify important interactions in the progression of disease. We used a quantitative study with $n=442$ diabetes patients, each characterized by $10$ standardized baseline features, including age, sex, body mass index, average blood pressure, and six blood serum measurements (Fig. 4b). The task is to construct an ML model to predict the response of interest, which is a quantitative measure of disease progression one year after the baseline. We considered four different types of ML models: MLP, LightGBM, XGBoost, and factorization machines (FM). For robustness, we repeated each experiment 20 times using different random seeds. Each repetition involved knockoff generation, ML model training, and interaction-wise FDR estimation.

Our analysis shows that Diamond consistently identifies the important interaction between body mass index (BMI) and serum triglycerides level (STL) across three of the four different ML models (Fig. 4a). The only exception is the FM method, for which Diamond fails to identify any important interactions. We investigated the identified interaction via literature search and found that high BMI and high STL are associated with an increased risk of diabetes, as reported by multiple studies [17, 59, 20]. The literature evidence is further supported by the qualitative evaluation (Fig. 4c). Specifically, we examined the identified interaction in four ways: the contribution of feature values to response prediction, the marginal importance measure, the interaction importance measure, and the contribution of the marginal importance measures to the interaction importance measure. This analysis showed that high BMI and high STL contribute to more severe progression of diabetes. Meanwhile, higher STL further reinforces the contribution of BMI to the progression of diabetes synergistically.

Finally, we find that Diamond does not necessarily report a lower FDR for interactions composed of two marginally important features. All $10$ baseline features from diabetes patients are effective in predicting disease progression, as measured by the Expected Gradient scores in the MLP model (Fig. 4b). Among these features, in addition to BMI and STL, blood pressure, blood sugar level, and high-density lipoproteins are more predictive than others. It is important to note that Diamond still associates most interactions composed of these highly predictive features with a high FDR threshold.

We then evaluated the performance of Diamond on a Drosophila enhancer dataset [3] to investigate the relationship between enhancer activity and DNA occupancy for transcription factor (TF) binding and histone modifications in Drosophila embryos. We used a quantitative study of DNA occupancy for $23$ TFs and $13$ histone modifications with the enhancer status labeled for $7,809$ genomic sequence samples from blastoderm Drosophila embryos (Fig. S.4). The enhancer status for each sequence is binarized as the response, depending on whether the sequence drives patterned expression in blastoderm embryos. To predict the enhancer status, the maximum value of normalized fold-enrichment [27] of ChIP-seq and ChIP-chip assays for each TF or histone modification served as our features. For robustness, we repeated the experiment 20 times using different random seeds.

We started by comparing the identified interactions against a list of well-characterized interactions in the early Drosophila embryos, summarized by [3]. If Diamond is successful in identifying important interactions while controlling the FDR, then the identified interactions should considerably overlap with the ground truth list. Our results show that three of the top five interactions reported by Diamond are present in the ground truth list for MLP and tree-based models, while four out of five interactions are present in the ground truth list for FM (Fig. 5a). In comparison, the baseline methods exhibit significantly different behaviors. One method, which is based on permutation-based interaction-wise p-values and the Benjamini-Hochberg procedure, reported no ground truth interactions among its top five identifications. Another method, which aggregates feature-wise FDR, reported both ground truth along with a large number of non-ground truth interactions.

Next, we investigated the identified interactions that are not included in the ground truth list. For example, we investigated the interaction between the TFs Snail (UniProt ID: P08044) and Twist (UniProt ID: P10627) identified by Diamond+MLP through databases containing assorted experimentally verified interactions [30] and literature [21]. Research indicates that Snail represses Twist targets because their target sequences are very similar, and their binding is mutually exclusive. This literature evidence is further supported by a qualitative evaluation (Fig. 5b-d). Specifically, we examined the Snail–Twist interaction in terms of the contribution of feature values to the enhancer status prediction, observing that high Twist expression suppresses the contribution of Snail to enhancer activation.

Furthermore, we investigated the remaining identifications that lack supporting ground truth or literature evidence. We discovered that these interactions can be logically explained through transitive effects. Specifically, even without a direct interaction between TF1 and TF3, we might still expect a strong interaction between them if solid interactions exist between TF1 and TF2 and between TF2 and TF3. On this basis, for example, the interaction between the TFs Twist (UniProt ID: P10627) and Zeste (UniProt ID: P09956), identified by Diamond+MLP, can be classified as a transitive interaction, supported by experimentally validated interactions [57]: (1) the interaction between Twist and Ultrabithorax (UniProt ID: P83949), Decapentaplegic (UniProt ID: P07713), Pho (UniProt ID: Q8ST83), Raf (UniProt ID: P11346), Psc (UniProt ID: P35820), Trr (UniProt ID: Q8IRW8) and (2) the interaction between these proteins and Zeste. Analogously, the interaction between the TFs Twist (UniProt ID: P10627) and Dichaete (UniProt ID: Q24533), identified by Diamond together with two tree-based models, can also be supported by transitive interactions [57]: (1) the interaction between Twist and Pho, Raf, Psc, Trr, Myc (UniProt ID: Q9W4S7), Trx (UniProt ID: P20659 ) and (2) the interaction between these proteins and Dichaete.

Finally, Diamond is able to identify non-additive interactions with various interactive patterns (Fig. 5b-d). Among the top five interactions reported by Diamond with MLP, the Zelda–Twist interaction, the Bicoid–Twist interaction, and the Krueppel–Twist interaction exhibit synergistic effects. The high expression of both TFs in these interactions reinforces the contribution of each individual factor to enhancer activation. In contrast, the Snail–Twist interaction, which was mentioned earlier, exhibits repressive effects, where high expression of one TF suppresses the contribution of the other to enhancer activation. Lastly, the Twist–Zeste interaction exhibits alternative effects, where high expression of either TF enhances the contribution of each to enhancer activation.

We lastly evaluated the performance of Diamond on a real mortality risk dataset [10] to investigate the relationship between mortality risk factors and long-term health outcomes within the US population. We used a mortality dataset from the National Health and Nutrition Examination Survey (NHANES I) and NHANES I Epidemiologic Follow-up Study (NHEFS). The dataset incorporates $35$ clinical and laboratory measurements from $14,407$ US participants between 1971 and 1974 (Fig. 6b). Note that some of the features are presented categorically, which makes model’s trained from this data challenging to interpret. We addressed this problem by converting each categorical feature into a list of binary features using one-hot encoding. The dataset also reports the mortality status of participants as of 1992, with $4,785$ individuals known to have passed away before 1992. The task is to construct an MLP-based model to predict mortality status using a survival loss function. For robustness, we repeated the experiment 20 times using different random seeds.

Given the absence of ground truth interactions for this task, we began by evaluating the interactions identified by Diamond through literature support. Our analysis shows that three of the top ten selected interactions are directly supported by existing literature (Fig. 6a). For example, we investigated the interaction between sex and sedimentation rate (SR) and found that, although a high SR is strongly associated with a higher risk of overall mortality, this association can be attenuated by sex [15]. Additionally, we investigated the interaction between creatinine and the blood urea nitrogen (BUN) identified by Diamond. Research indicates that the BUN/creatinine ratio is recognized to have a nonlinear association with all-cause mortality and a linear association with cancer mortality [46]. The literature evidence is further supported by qualitative evaluation (Fig. 5c-e). Specifically, we examined the BUN-creatinine interaction in terms of the contribution of feature values to mortality status prediction, observing that when the creatinine level exceeds a certain level, a high BUN level further increases the mortality risk.

Finally, we probed into the remaining identifications that lack supporting literature evidence. As in the Drosophila analysis, we discovered that these interactions can be logically explained through transitive effects. For example, the interaction between BUN and potassium can be justified by combining the following two established facts: (1) BUN level is indicative of chronic kidney disease development [9], and (2) patients with chronic kidney diseases show a progressively increasing mortality rate with abnormal potassium levels [44]. Additionally, the interaction between BUN and SR reflects the fact that (1) the creatinine level and SR are strongly associated to clinical pathology types and prognosis of patients [28], and (2) creatinine and BUN have a nonlinear association [46].

Diamond achieves FDR control by leveraging the model-X knockoffs framework [2, 6], which was proposed in the setting of error-controlled feature selection. The core idea of knockoffs is to generate dummy features that perfectly mimic the empirical dependence structure among the original features but are conditionally independent of the response given the original features. Briefly speaking, the knockoff filter achieves FDR control in two steps: (1) construction of knockoff features and (2) filtering using knockoff statistics.

For the first step, the knockoff features are defined as follows:

According to Definition 1, the construction of the knockoffs must be independent of the response $\mathbf{Y}$. Thus, if we can construct a set $\tilde{X}$ of model-X knockoff features properly, then by comparing the original features with these control features, FDR can be controlled at target level $q$. In the Gaussian setting, i.e., $\mathbf{X}\sim\mathcal{N}(0,\mathbf{\Sigma})$ with covariance matrix $\mathbf{\Sigma}\in\mathbb{R}^{p\times p}$, the model-X knockoff features can be constructed easily:

where $\mathrm{diag}\{\mathbf{s}\}$ is a diagonal matrix with all components of $\mathbf{s}$ being positive such that the conditional covariance matrix in Equation 1 is positive definite. As a result, the original features and the model-X knockoff features constructed by Equation 1 have the following joint distribution:

With the constructed knockoff $\tilde{\mathbf{X}}$, feature importances are quantified by computing the knockoff statistics $W_{j}=g_{j}(Z_{j},\tilde{Z}_{j})$ for $1\leq j\leq p$, where $Z_{j}$ and $\tilde{Z}_{j}$ represent feature importance measures for the $j$-th feature $X_{j}$ and its knockoff counterpart $\tilde{X}_{j}$, respectively, and $g_{j}(\cdot,\cdot)$ is an antisymmetric function satisfying $g_{j}(Z_{j},\tilde{Z}_{j})=-g_{j}(\tilde{Z}_{j},Z_{j})$. The knockoff statistics $W_{j}$ should satisfy a coin-flip property such that swapping an arbitrary pair $X_{j}$ and its knockoff counterpart $\widetilde{X}_{j}$ only changes the sign of $W_{j}$ but keeps the signs of other $W_{k}$ ($k\neq j$) unchanged [6]. A desirable property for knockoff statistics $W_{j}$’s is that important features are expected to have large absolute values, whereas unimportant ones should have small symmetric values around 0.

Finally, the absolute values of the knockoff statistics $|W_{j}|$’s are sorted in decreasing order, and FDR-controlled features are selected whose $W_{j}$’s exceed some threshold $T$. In particular, the choice of threshold $T$ follows $T=\min\left\{t\in\mathcal{W},\frac{1+|\left\{j:W_{j}\leq-t\right\}|}{|\left\{j:W_{j}\geq t\right\}|}\leq q\right\}$ where $\mathcal{W}=\left\{|W_{j}|:1\leq j\leq p\right\}\char 92\relax\left\{0\right\}$ is the set of unique nonzero values from $|W_{j}|$’s and $q\in(0,1)$ is the desired FDR level specified by the user.

As a key precursor to FDR estimation, Diamond quantifies feature interactions from trained ML models and produces a ranked list of these interactions, with higher-ranked interactions indicating greater importance. For notational simplicity, we use indices for both original features and knockoffs as $\left\{1,2,\cdots,2p\right\}$, with $\left\{1,\cdots,p\right\}$ and $\left\{p+1,\cdots,2p\right\}$ corresponding to the original features and their respective knockoff counterparts. Here, we define $\mathbf{E}^{\text{2D}}=\left[e_{ij}\right]_{i,j=1}^{2p}\in\mathbb{R}^{2p\times 2p}$ as a reported interaction importance measure from existing methods. There are many feature interaction importance measures available for $\mathbf{E}^{\text{2D}}$, each attributing the prediction influence to feature pairs in different ways. However, it is important to note that such measures favor pairs where both features are simultaneously important for a model’s prediction, rather than capturing the true non-additive or interactive effects between the two features [54]. Further supported by simulation studies (Fig. 2c), we observed that many feature interaction importance measures tend to assign higher interaction scores to two marginally important but non-interacting features compared to two random ones, even though neither pair has a real interaction, leading to the failure of FDR control.

The direct reason for the interaction importance measure falling short in controlling FDR is its violation of the knockoff filter’s assumption. Specifically, the knockoff filter requires that the importance scores of knockoff-involving interactions and false interactions have a similar distribution. To resolve this issue, we introduce a calibration procedure to be applied on top of existing interaction importance measures. Specifically, we consider that a reported interaction importance measure from existing methods comprises a mixture of several factors: prediction-dependent marginal effects for individual features, prediction-independent feature biases, independent random noise, and potential non-additive interactive effects between feature pairs. Thus, the reported interaction between features $i$ and $j$ is represented as:

where $\varepsilon_{ij}\in\mathbb{R}$ is random noise independent of both features and predictions, $e_{ij}\in\mathbb{R}$ and $e_{i},e_{j}\in\mathbb{R}$ are reported pairwise and univariate feature importance measures that are dependent on the model’s predictions, respectively. The functions $g_{i},g_{j}:\mathbb{R}\mapsto\mathbb{R}$ adapt univariate feature importance to be compatible with feature interaction importance. The function $b:\mathbb{R}^{2p}\mapsto\mathbb{R}$ models the feature-specific biases that are independent of the model’s predictions, where $I_{ij}\in\left\{0,1\right\}^{2p}$ indicates the presence of feature $i$ and $j$. Our goal is to identify $s_{ij}$, the potential non-additive interactive effects between features $i$ and $j$.

We formulate the identification of interactive effects $s_{ij}$ as the residuals of a regression task:

where $w_{ij}>0$ is the conditional probability of being either original-only (i.e., $i,j\leq p$) or knockoff-involving (i.e., $i>p$ or $j>p$) feature pairs given two univariate feature importance measures $e_{i}$ and $e_{j}$, estimated by a logistic regression model [16]. The rationale is based on the important observation that most feature pairs do not exhibit non-additive interactions, especially those involving knockoff features. Therefore, we want to focus more on potential non-additive interactions that have large univariate feature importance, as important interactions naturally consist of significant marginal features. In this study, we parameterize the functions $b:\mathbb{R}^{2p}\mapsto\mathbb{R}$ and $g_{i},g_{j}:\mathbb{R}\mapsto\mathbb{R}$ using generalized additive models and optimize Eq. 4 using the pyGAM library [45].

After calculating the non-additive interactive effects using Eq. 4, we denote the resultant set of interactive effects as $\Gamma=\left\{s_{ij}|i<j,i\neq j-p\right\}$. We arrange $\Gamma$ in decreasing order and select interactions for which the interactive effect, $\Gamma_{j}$, exceeds some threshold, $T$. This selection ensures that the chosen interactions adhere to a desired FDR level $q\in(0,1)$.

However, a point of complexity is introduced due to the heterogeneous interactions, which include original-only and knockoff-involving interactions. The latter further comprises original-knockoff, knockoff-original, and knockoff-knockoff interactions. Following the strategy outlined by [55], the threshold $T$ is determined by:

where $\mathcal{D}$ and $\mathcal{DD}$ respectively denote the sets of interactions that include at least one knockoff feature and both knockoff features, while $\mathcal{T}$ refers to the set of unique nonzero values present in $\Gamma$.

Diamond is designed to be compatible with a wide range of ML models. To demonstrate the broad applicability of Diamond, we use it in conjunction with machine learning models from methodologically distinct categories, including DNN models, tree-based models, and factorization-based models.

For DNNs, we configured a pyramid-shaped multi-layer perceptron (MLP) model with the exponential linear unit activation and four hidden layers, each having neuron sizes of 2p, p, p/2, and p/4 respectively, where p denotes the input dimensionality. The model follows the guidance of [33] and includes a plugin pairwise-coupling layer, connecting each original feature with its knockoff counterpart in a pairwise fashion to maximize statistical power (See details in Sec. S.1). Using an MLP as the ML model, we employed three representative interaction importance measures to clarify the relationship between features and responses. We began with a model-specific interaction importance measure derived from the model weights, specifically designed for the MLP architecture (See details in Sec. S.2). To demonstrate Diamond’s flexibility, we also used Integrated Hessian [22] and Expected Hessian [14], two state-of-the-art, model-agnostic interaction importance measures, to elucidate the relationships between features and responses without making assumptions about any specific model architecture (See details in Sec. S.3). Both model-agnostic interaction importance measures are calculated using the Path-Explain library [22].

For tree-based models, we employed two widely used gradient-boosted decision tree representatives—XGBoost [8] and LightGBM [24]. We used the implementations provided by the two libraries with default settings. Using either XGBoost or LightGBM as the ML model, we utilized TreeSHAP [37], a Shapley value-based interpretation method tailored for tree-based ML models, to elucidate the relationships between features and responses. We used the implementation provided by the SHAP library [36] to calculate both univariate and pairwise feature importance.

For factorization-based models, we used a widely adopted representative method—Factorization machines (FM) [41]. We used 2-Way FM implemented in the xLearn library [39] with default settings, which models up to second-order feature interactions. Since 2-Way FM is explicitly modelled as the weighted sum of univariate and pairwise feature interactions, we used the learned coefficients of these interactions as the corresponding importance measures.

Diamond relies on knockoffs to control the FDR. It is worth mentioning that conventional knockoffs are limited to Gaussian settings, which may not be applicable in many practical scenarios. In this study, we focus on the state-of-the-art knockoff design without assuming any specific feature distribution. In particular, we consider commonly-used non-Gaussian knockoff generation methods such as KnockoffsDiagnostics [5], KnockoffGAN [23], Deep knockoffs [42], and VAE knockoffs [31].

We evaluate the performance of Diamond in comparison to other baseline methods. For the first baseline method, we employ a permutation-based approach to calculate the interaction-wise FDR. Specifically, this involves using a previously described permutation procedure tailored for neural networks to assess the significance of interactions and calculate permutation p-values [12], followed by the Benjamini–Hochberg procedure [4] to estimate the FDR.

For the second baseline method, we consider an ensemble-based approach that represents interaction-wise FDR as the aggregation of feature-wise FDR. This approach follows the intuition that an important feature interaction is composed of important univariate features. Specifically, we use a previously described knockoff-based procedure tailored for neural networks to estimate the feature-wise FDR of each univariate feature [33]. We then approximate the interaction-wise FDR as the maximum of the two comprising univariate feature-wise FDRs.

The software implementation and all models described in this study are available at https://github.com/batmen-lab/diamond. All public datasets used for evaluating the models are available at https://github.com/batmen-lab/diamond/data.