ANOVA-NODE: An identifiable neural network for the functional ANOVA model for better interpretability

Weakness 2. Hyperparameter Tuning.

As suggested by the reviewer, we additionally conducted an experiments of prediction performance for NODE and XGB by tuning the comprehensive hyper-parameters. We turned NODE by grid search over the following values.

$\bullet$ num layers = [2,4,8]

$\bullet$ tree depth = [6,8]

$\bullet$ # of trees in each layer = [256,512]

Similarly, for XGB, we used below values for the grid search.

$\bullet$ learning rate = [0.0001, 0.005, 0.01, 0.05 , 0.1]

$\bullet$ max depth = [3, 5, 7]

$\bullet$ # of trees = [50,100,200,300,400,500,600,700,800,900,1000]

Table D.2 Results of prediction performance and hyper-parameters in XGB.

Table D.3 Results of prediction performance and hyper-parameters in NODE.

Tables D.2 and D.3 present the best prediction performance obtained through grid search for the XGB and NODE models, along with the corresponding hyper-parameters. Based on the results of the above experiments, we have revised the experimental results in the paper.

Weakness 3.Code Availability.

We have included the source code of ANOVA-NODE in the Supplementary Material.

Weakness 4. Baseline Comparison.

As suggested by the reviewer, we conducted additional experiment of prediction performance in simple decision tree by using sklearn python package. Table D.4 presents the results of prediction performance in simple decision tree. We add these results in Appendix I of the paper.

Table D.4 Results of prediction performance in decision tree.

Weakness 4. Related Work.

ANOVA-NODE estimates each component using the ANOVA-NODT. It is a modified version of NODT introduced in Section 2.3, which satisfies the sum-to-zero condition. On the other hand, NAM estimates each component using a deep neural network which does not satisfy sum-to-zero condition.

Weakness 5. Complexity of Presentation.

As pointed out by the reviewer, we revised the content corresponding to page 5 to make it clearer, more concise, and easier to understand in the revision.

We appreciate your feedback on our paper and have made every effort to address your insightful questions. We have uploaded the revised version of the paper, reflecting the reviewer's feedback.

Weakness 1. Motivation and Relevance.

In Section 3.1, the example which addresses the importance of identifiability was discussed. Moreover, the plots of the components estimated by ANOVA-NODE, which satisfies the sum-to-zero condition, and by NAM and NBM, which do not, are presented in Appendix F.1. While ANOVA-NODE estimates the components stably, NAM and NBM do not.

Furthermore, in NA$^{2}$M, some main effects are estimated as a constant function, which would be partly because the main effects are absorbed into the second order interactions. However, such an issue does not occur in ANOVA-NODE.

The advantages of ANOVA-NODE over SHAP are as follows.

1.1. Functional relationship.

Unlike SHAP, ANOVA-NODE can deterministically identify the relationship between the components and the output. Therefore, in ANOVA-NODE, we can accurately assess the contribution of each component to the output.

1.2. Decomposition of components.

While SHAP provides an interpretation of the feature by considering all higher-order interactions that include the feature, ANOVA-NODE interprets each component individually, offering a more specific analysis. For example, consider $f(X_{1},X_{2}) =X_{1} + X_{2} + X_{1}X_{2}$.

Then, for given data points X$^{1}=(1,1)^{\top}$ and X$^{2}=(2,100)^{\top}$, then the SHAP value $\phi_{1}(f,$X$^{1})$ and $\phi_{1}(f,$X$^{2})$ is calculated as

$\phi_{1}(f,$X$^{1}) = {3\over 2} , \quad \phi_{1}(f,$X$^{2}) = 102$

In this case, while SHAP only provides an interpretation that combines both the main effect and interactions, ANOVA-NODE can provide interpretations for both the main effect and interactions via estimated component functions $(f_{1}, f_{1,2} )$.

1.3. Stability

We compare the stability of the local interpretation from SHAP and ANOVA-NODE on Wine dataset. As in the experiments of the paper, we randomly split the data, trained a XGB model with 800 trees, a learning rate of 0.005, and a max depth of 5, and calculated the SHAP values via Tree-SHAP, repeating this process a total of 10 times. We calculated the stability score, as defined in our paper, for the SHAP values from the 10 trials.

Table D.1 Results of stability score in ANOVA-N$^{1}$ODE + ANOVA-SHAP and XGB + SHAP.

Table D.1 presents the results of stability scores for local interpretations obtained using ANOVA-SHAP in ANOVA-NODE and Tree-SHAP in XGB. We observed that the local interpretation obtained using ANOVA-SHAP in the ANOVA-N$^{1}$ODE is approximately seven times more stable than that of SHAP in XGB.

Weakness 6. Runtime and Efficiency.

We conducted additional experiments on the runtime of NAM, NBM, ANOVA-NODE and NBM-NODE where NBM-NODE is an extended model of ANOVA-NODE, proposed in Appendix G.2 of this paper, which improves the scalability of ANOVA-N$^{1}$ODE by using basis functions, similar to the approach used in NBM. Note that in NBM-NODE, the number of ANOVA-NODT for basis does not depend on the number of components

We conducted additional experiments to assess the improvement in scalability. We consider NA$^{1}$M, which has 3 hidden layers with 16, 16, and 8 units; 10 basis DNNs for NB$^{1}$M, which has 3 hidden layers with 32, 16, and 16 units; 10 trees for each component in ANOVA-N$^{1}$ODE; and 10 basis functions in NBM-N$^{1}$ODE.

Table D.5 Results of computaitonal complexity

Table D.5 presents the runtime for each model when trained on a Abalone dataset and 100 epochs. Our computational environment consists of RTX 3090 and RTX 4090. We observed that our model has fewer weight parameters and takes less runtime compared to existing models (NAM and NBM).

Weakness 7. Performance Reporting.

As suggested by the reviewer, we report the unnormalized stability scores in Appendix G.2 of paper. The results reported in Table 2 are the AUROCs for component selection, not the ranks.

Question 1.

In particular, the black-box models XGB and NODE show good prediction performance, which is due to the presence of important higher-order interactions of degree 3 or higher in these two datasets (Calhousing and Wine dataset).

Question 2.

In Appendix E.1.1, ANOVA-NODE identified the important main effects in the following order: 7, 8, 1, 6, 3, 2, 4, 5. ANOVA-N$^{2}$ODE interpreted main effect 6 as the most important, because the effects of features 1, 7 and 8 in ANOVA-N$^{1}$ODE are decomposed into interactions such as (7,8), (1,7), (3,8), (1,8), (4,8), (2,7) and (1,5), whose importance scores are significant in ANOVA-N$^{2}$ODE. Therefore, the importance scores of main effects 1, 7, and 8 are relatively lower in ANOVA-N$^{2}$ODE.

Question 3.

If a linear model is used instead of ANOVA-NODE within the CBM framework, the interpretation of features for each image is limited to linear relationships, whereas using ANOVA-NODE allows for more flexible interpretations.

The prediction performance results when using a linear model instead of ANOVA-N$^{1}$ODE in CBM are shown in Table D.6. The target variable is set as 'gender'.

Table D.6 Results of prediction performance in CBM + Linear and CBM + ANOVA-N$^{1}$ODE.

We observed a decrease in prediction performance when using a linear model instead of ANOVA-N$^{1}$ODE in CBM. Additionally, the details of how the estimated function becomes monotone due to the monotone constraint are described in Appendix B.1.

W4.

The results in Table D.4 are from a Decision Tree model trained with the default hyper-parameters set in the sklearn package. Additionally, we turned hyper-parameter f the Decision Tree using the Optuna package in Python. We applied Optuna to the following hyper-parameter range.

$\bullet$ Range of max depth = [2 ,12]

$\bullet$ Range of min_samples_leaf = [2,10]

$\bullet$ Range of min_samples_split = [2,10]

$\bullet$ Range of max_leaf_nodes = [2,10]

Table D.9 presents the results of prediction performance in turned decision tree by Optuna package in Python. Additionally, the results obtained from the Decision Tree could not be included in the table of results for other baseline models in the main text, as doing so would exceed the paper's page limit. We have replaced the results in Appendix I with the tuned results.

Table D.9. Results of prediction performance in decision tree on the optimized hyper-parameters by the Optuna python package.

W6.

As suggested by the reviewer, we conducted runtime experiments not only on the Abalone dataset, but also on the Calhousing and Online datasets, which can be found in Table D.10. We observed that the runtime of NBM-NODE, which improves the scalability of ANOVA-NODE, is the shortest. The experimental results for runtime have been added in Appendix J of the revision.

Table D.10. Results of runtime in each model on Abalone, Calhousing, and Online dataset.

W7.

Unlike other models (such as Deep Neural Networks), in explainable models based on the Functional ANOVA model, which estimate components and provide interpretations using the estimated components, stability is crucially important. [1] discussed the stability of component estimation in Section 4.6 and compared the stability of NA$^{1}$M and NB$^{1}$M on Calhousing dataset.

References

[1] Radenovic, Filip, Abhimanyu Dubey, and Dhruv Mahajan. "Neural basis models for interpretability." Advances in Neural Information Processing Systems 35 (2022): 8414-8426.

We appreciate your valuable comments.

W1.

As suggested by the reviewer, a figure 1 corresponding to the motivation has been added to the introduction to enhance the clarity of the paper's contributions.

W2.

As suggested by the reviewer, we conducted additional experiments on the hyper-parameter optimization of XGB. As the first approach, we fixed the number of trees to 1000 and performed a grid search for the learning rate, max_depth, alpha, and lambda using the values below.

$\bullet$ learning rate = [0.001, 0.005, 0.01, 0.05, 0.1]

$\bullet$ max depth = [3, 5, 7]

$\bullet$ alpha = [1e-6, 1e-4, 1e-2, 1.0]

$\bullet$ lambda= [1e-6, 1e-4, 1e-2, 1.0]

Table D.7 presents the results of prediction performance in XGB on the optimized hyper-parameters by the grid search.

Table D.7. Results of prediction performance in XGB on the optimized hyper-parameters by the grid search

As the second approach, we used the Optuna package in Python to tune the hyper-parameters of XGB. The ranges of hyper-parameters for Optuna are as follow.

$\bullet$ The range of # of tree = [100,1000]

$\bullet$ The range of max depth = [2,12]

$\bullet$ The range of learning rate = [0.001, 1.0]

$\bullet$ The range of lambda= [1e-8, 1.0]

$\bullet$ The range of alpha= [1e-8, 1.0]

Table D.8 presents the results of prediction performance in XGB on the optimized hyper-parameters by the Optuna python package.

Table D.8. Results of prediction performance in XGB on the optimized hyper-parameters by the Optuna python package.

We observed that The results in Table D.2, D.7, and D.8 showed little difference in view of prediction performance. Moreover, for some datasets, we found that setting lambda and alpha to their default values of 0 and fine-tuning the number of trees, learning rate, and max depth through grid search yielded better results.

W3.

We have not only added the code as supplementary material, but also included explanations of the data preprocessing process for NAM, NBM, ANOVA-NODE, and NODE-GAM in Section 3 and Appendix C.2. Additionally, we have added an explanation of the challenges encountered during the optimization of ANOVA-NODE in Section 3.

We are grateful for your feedbacks, which have greatly enriched our work. As suggested by the reviewer, we have been provided more detailed information on the evaluation protocol and data preprocessing in Appendix C.2.

Overall impact of paper

Existing explainable models based on the Functional ANOVA model (e.g., NAM(Agarwal et al., 2021), NBM(Radenovic et al., 2022), etc.) provide interpretations using the estimated components. However, if the components are not identifiable, there are multiple interpretation for a given function. For example,

$f(x_{1},x_{2}) = f_{1}(x_{1}) + f_{2}(x_{2}) + f_{1,2}(x_{1},x_{2})$ $\quad \cdots (1)$

where $f_{1}(x_{1}) = x_{1}, f_{2}(x_{2}) = x_{2} , f_{1,2}(x_{1},x_{2}) = x_{1}x_{2}$ can be expressed as

$f(x_{1},x_{2}) = f_{1}^{\*}(x_{1}) + f_{2}^{\*}(x_{2}) + f_{1,2}^{\*}(x_{1},x_{2})$ $\quad \cdots (2)$

where $f_{1}^{\*}(x_{1}) = -x_{1}, f_{2}^{\*}(x_{2}) = x_{2} , f_{1,2}^{\*}(x_{1},x_{2}) = x_{1}(x_{2}+2)$.

In this case, the interpretations of main effect provided by case (1) and case (2) differs. Specifically, in case (1), it is interpreted that the main effect of $x_{1}$ provides a positive contribution to the output as $x_{1}$ increases, whereas in case (2), it is interpreted as providing a negative contribution. In other words, if the components are not identifiable, the interpretations provided by XAI models based on the Functional ANOVA model cannot be trusted.

In this paper, the main contributions of our work are as follows.

1. We proposed ANOVA-NODE which used a specially modified version of NODT. ANOVA-NODE satisfies specific conditions to ensure that the Functional ANOVA model we estimate is identifiable. As a result, the interpretations provided by our model are more reliable compared to those of NAM(Agarwal et al., 2021) and NBM(Radenovic et al., 2022).

2. Second, we proved a universal approximation theorem for our proposed model (ANOVA-NODE).

References

[1] Agarwal, Rishabh, et al. "Neural additive models: Interpretable machine learning with neural nets." Advances in neural information processing systems 34 (2021): 4699-4711.

[2] Radenovic, Filip, Abhimanyu Dubey, and Dhruv Mahajan. "Neural basis models for interpretability." Advances in Neural Information Processing Systems 35 (2022): 8414-8426.

We appreciate your feedback on our paper and have made every effort to address your insightful questions. We have uploaded the revised version of the paper, reflecting the reviewer's feedback.

Weakness N1.

We agree with the reviewer’s comment, and therefore, Theorem 3.1 and Theorem 3.2 have been changed to proposition 3.1 and 3.2. Additionally, references have been added. Note that the main contribution of our paper is not Theorem 3.1 and Theorem 3.2, but rather the proposal of interpretation XAI model that can be trained using gradient-based optimization methods while satisfying the sum-to-zero condition for the functioanl ANOVA model, along with the associated universal approximation theorem.

Weakness N2.

[1] proposed a post-processing method to satisfy the sum-to-zero condition for piecewise-constant function, whereas ANOVA-NODE is a model that is learned while inherently satisfying the sum-to-zero condition.

When applied to a given dataset of size n, the computation order of performing post-processing for a given point x is $\mathcal{O}(dn^{d-1})$, which is very demanding. Furthermore, performing post-processing requires storing a dataset, which causes memory efficiency issues. Therefore, interpretation, either locally or globally, is practically infeasible.

To be specific, for a given point x, the computation order of calculating ANOVA-SHAP $\phi_{j}(f,$x$)$ is $\mathcal{O}(p^{d-1}dn^{d-1})$. Therefore, not only is it impossible to draw dependence plots, but performing local interpretation using ANOVA-SHAP is also infeasible. The method for post-processing in GA$^{d}$M and the process for calculating its computational order have been added in the Appendix L.

Weakness NC3.1.

A universal approximation theorem also may exists for NODE-GAM, but the most important aspect of Theorem 3.3 is that the universal approximation theorem holds for ANOVA-NODE using ANOVA-NODTs that satisfy the sum-to-zero condition. Additionally, since the function space of ANOVA-NODE is smaller than that of NODE-GAM, proving an approximation theorem for ANOVA-NODE is a much more difficult problem than for NODE-GAM.

Weakness NC3.2.

The $entmax({**x** - b \over \gamma})$ in ANOVA-NODE corresponds to the basis function $B(**x**)$ in spline basis representation of [2]. While AHOFM ([2]) does not train the basis function, ANOVA-NODE can be considered to learn the basis by training $b , \gamma$.

Also, we acknowledge that if all higher-order interactions are considered, the method proposed by [2] is more efficient in terms of scalability compared to ANOVA-NODE. However, the experiments in the paper show that considering interactions up to the second order is sufficient for real data. Moreover, when all higher-order interactions are considered in the functional ANOVA model, it is likely to be overfitted.

To examine the relationship between considering all higher-order interactions and model overfitting, we trained ANOVA-N$^{3}$ODE and ANOVA-N$^{4}$ODE on the Calhousing dataset, and the results are shown in Table 3. We observed that considering all interactions of third-order or higher led to overfitting, resulting in a decrease in prediction performance.

Table C.1. Results of prediction performance in ANOVA-N$^{2}$ODE, ANOVA-N$^{3}$ODE, and ANOVA-N$^{4}$ODE.

Therefore, rather than considering all higher-order interactions, we should focus on the important ones. As demonstrated in the high-dimensional data experiments in our paper, this can be done by pre-selecting interactions using an interaction selection algorithm like Neural Interaction Detection ([3]). When these pre-selected interactions are applied to ANOVA-NODE, there is little difference in scalability compared to method in [2].

Weakness NC3.3.

We thank the reviewer for pointing out the errors. The errors have been corrected.

Weakness E.1.

In Appendix I, we confirmed that the prediction performance is better when sum-to-zero conditions are not imposed, and this was discussed.

Weakness E.2.

As suggested by the reviewer, we conducted additional experiments GAMs by using pygam python package, denote it as py-GAM. py-GA$^{d}M estimates each component of GA$^{d}$ by a cubic spline model. Table C.2 shows the average and standard deviation of the RMSE from 10 trials on the Abalone dataset.

Table C.2. Results of prediction performance in ANOVA-N$^{1}$ODE, ANOVA-N$^{2}$ODE, and ANOVA-N$^{3}$ODE.

Weakness E.4.

The sum-to-zero condition is based on the product measure, so $Var(f)= \sum_{S}Var(f_{S})$ does not hold in ANOVA-NODE. Therefore, in the paper, instead of conducting experiments of the analysis of variance, we conducted experiments to analysis the importance components using $l1$ norm, and the experiments results on Calhousing dataset are presented in Section E.1.1.

Weakness P.1.

As suggested by the reviewer, we restructured the paper and uploaded it.

Weakness P.2.

``is neither a unique nor optimal identifiability condition" means that there can be various types of constraints, other than the sum-to-zero condition, that ensure the identifiability of the functional ANOVA decomposition. Because the sum-to-zero condition also ensures the identifiability of components, the initially stated problem is solved.

Weakness P.3.

As suggested by the reviewer, we added axis labels in Figure 2 and changed the term 'similarity measure' to 'stability score'.

Minor l.161.

As suggested by the reviewer, we explicitly express the $f$ into $f_{1},f_{2},f_{1,2}$.

Minor l.285.

We have added a more detailed description of the challenges in optimization to the data preprocessing paragraph in the revision.

Minor l.302.

As suggested by the reviewer, we have replaced all $\mathcal{S}$ symbol with the $\mathcal{SC}$ symbol.

Minor Table 1: Both Table 1 and Figure 2 discuss stability but show different value ranges.

We add the results of unnormalized stability scores in Appendix G.2 of the revision.

Minor Figures 3, 4, and those in the Appendix.

Following the reviewer's suggestion, we have increased the text size in all the figures.

Minor Figures 3, 4, and those in the Appendix.

As suggested by the reviewer, we have added a more detailed comparison with NODE-GAM and included it in Appendix K of the revision.

Minor Manuscript:

The typo have been corrected.

Question 1.

It has already been addressed in the response to the Weakness section.

Question 2. Would it make sense to use an ANOVA kernel as in Blondel et al., 2016 to improve the scaling of ANOVA-N ODES? Question 3. Could a factorization approach as in Rügamer, 2024 incorporated to improve scalability?

Since ANOVA-NODT is also a model represented by the spline basis representation of equation (2) in [2](entmax can be viewed as basis) factorization method can be used to improve the scalability of ANOVA-NODE. Exploring the use of the factorization approach to improve the scalability of ANOVA-NODE seems interesting.

Note that, in our paper, we proposed NBM-NODE which is an extended model of ANOVA-NODE, proposed in Appendix G.2 of this paper, which improves the scalability of ANOVA-NODE by utilizing basis functions, similar to the approach used in NBM. Note that in NBM-NODE, the number of ANOVA-NODT for basis does not depend on the number of components

We conducted additional experiments to assess the improvement in scalability. We consider NA$^{1}$M, which has 3 hidden layers with 16, 16, and 8 units; 10 basis DNNs for NB$^{1}$M, which has 3 hidden layers with 32, 16, and 16 units; 10 trees for each component in ANOVA-N$^{1}$ODE; and 10 basis functions in NBM-N$^{1}$ODE.

Table C.3 Results of computaitonal complexity

Table C.3 presents the runtime for each model when trained on a Abalone dataset and 100 epochs. Our computational environment consists of RTX 3090 and RTX 4090. We observed an improvement in scalability when we extended ANOVA-NODE using basis functions.

References

[1] Lengerich, Benjamin, et al. "Purifying interaction effects with the functional anova: An efficient algorithm for recovering identifiable additive models." International Conference on Artificial Intelligence and Statistics. PMLR, 2020.

[2] Rügamer, David. "Scalable Higher-Order Tensor Product Spline Models." International Conference on Artificial Intelligence and Statistics. PMLR, 2024.

[3] Tsang, Michael, Dehua Cheng, and Yan Liu. "Detecting statistical interactions from neural network weights." arXiv preprint arXiv:1705.04977 (2017).

Dear Authors,

Thank you for providing such detailed feedback. Due to time constraints, I have only skimmed parts of your response. I will follow up on some of your answers hopefully by the end of the day (AoE). However, what I have read so far already addresses many of my concerns. I will raise my score at least to 5, potentially higher after reading your response more carefully.

One small request: Could the authors maybe fix the formatting of their post (for example, the \textbf is not rendering correctly in markdown, I suggest using **[TEXT]**) to improve its readability? Edit: Thanks, looks good now.

Thank you for your valuable feedback. It has greatly enriched our work. As suggested by the reviewer, we will add the discussion of related literature to Appendix and upload it by tomorrow.

Question 4. what the paper's novelty

ANOVA-NODE uses ANOVA-NODT, a modified version of Neural Oblivious Decision Tree (NODT) that satisfies the sum-to-zero condition. Of course, there may be other architectures that can be modified to satisfy the sum-to-zero condition, but it cannot be guaranteed that these modified architectures will increase the stability of component estimation practically. In other words, while the modified model ensures the uniqueness of the functional ANOVA decomposition theoretically, the components estimated by the modified model that satisfies the sum-to-zero condition may vary for each randomly sampled training dataset.

For example, we discussed the post-processing for the sum-to-zero condition in Appendix M of the paper (revision version). We observed that while the neural network-based models, NAM(Agarwal et al., 2021) and NBM(Radenovic et al., 2022), did not show improved stability in component estimation even after post-processing, the stability was improved in the case of GAM-NODE, which estimates components using NODT.

Therefore, we believe that the stability of component estimation in ANOVA-NODE is due to the sum-to-zero condition and the characteristics of the NODT model structure. The relationship between the sum-to-zero condition and NODT will be explored in future work.

Question 5. how it differs from existing literature

The model (Rügamer, 2024 ) mentioned by the reviewer is a highly efficient algorithm for estimating the functional ANOVA model, but it does not satisfy the sum-to-zero condition or other conditions for identifiability.

Therefore, this model may provide multiple interpretations for a given data point x. For example,

$f(x_{1},x_{2}) = f_{1}(x_{1}) + f_{2}(x_{2}) + f_{1,2}(x_{1},x_{2})$ $\quad \cdots (1)$

where $f_{1}(x_{1}) = x_{1}, f_{2}(x_{2}) = x_{2} , f_{1,2}(x_{1},x_{2}) = x_{1}x_{2}$ can be expressed as

$f(x_{1},x_{2}) = f_{1}^{\*}(x_{1}) + f_{2}^{\*}(x_{2}) + f_{1,2}^{\*}(x_{1},x_{2})$ $\quad \cdots (2)$

where $f_{1}^{\*}(x_{1}) = -x_{1}, f_{2}^{\*}(x_{2}) = x_{2} , f_{1,2}^{\*}(x_{1},x_{2}) = x_{1}(x_{2}+2)$.

In this case, the interpretations of main effect provided by case (1) and case (2) differs. Specifically, in case (1), it is interpreted that the main effect of $x_{1}$ provides a positive contribution to the output as $x_{1}$ increases, whereas in case (2), it is interpreted as providing a negative contribution. In other words, if the components are not identifiable, the interpretations provided by models based on the Functional ANOVA model cannot be trusted.

To solve this issue, post-processing would need to be performed, but due to the high computational complexity, it is practically infeasible.

Furthermore, as mentioned in the response to Question 4, even if these models (Blondel et al., 2016 and Rügamer, 2024) are modified to satisfy the sum-to-zero condition, there is no guarantee that component estimation will be stable for each randomly sampled training dataset.

References

[1] Agarwal, Rishabh, et al. "Neural additive models: Interpretable machine learning with neural nets." Advances in neural information processing systems 34 (2021): 4699-4711.

[2] Radenovic, Filip, Abhimanyu Dubey, and Dhruv Mahajan. "Neural basis models for interpretability." Advances in Neural Information Processing Systems 35 (2022): 8414-8426.

[3] Rügamer, David. "Scalable Higher-Order Tensor Product Spline Models." International Conference on Artificial Intelligence and Statistics. PMLR, 2024.

[4] Blondel, Mathieu, et al. "Higher-order factorization machines." Advances in neural information processing systems 29 (2016).

Thank you for your additional questions. We have done our best to respond to your thoughtful questions

Question 6. Do you mean that the model estimation exhibits high variance?

We explain the points we mentioned earlier in more detail. For a randomly sampled training dataset, we can obtain the components $f_{S}$ estimated by the model, and by repeating this process $k$ times, we obtain $(\hat{f}^{i}_{S} , i=1,...,k )$ for $S \subseteq [p]$.

​In this case, for any $S \subseteq [p]$, if the components $(\hat{f}^{i}_{S} , i=1,...,k )$ are estimated similarly, we say that the model is estimating the component stably.

Question 7. quick follow-up question.

Section 5.6.2 of Wood (2017) describes the SSANOVA (Gu, Chong, 2014). SSANOVA estimates the components in functional ANOVA model by using spline models which satisfy the sum-to-zero condition which we used. (They referred to the sum-to-zero condition as the averaging operator.)

As far as I know, py-GAM estimates the components of a Functional ANOVA model using a spline model without identifiability condition (e.g., sum-to-zero condition).

Using the averaging operator of SSANOVA, it seems possible to apply it to AHOFM( Rügamer, 2024), but further research is needed on this.

References

[1] Gu, Chong. "Smoothing spline ANOVA models: R package gss." Journal of Statistical Software 58 (2014): 1-25.

We appreciate your feedback on our paper and have made every effort to address your insightful questions. We have uploaded the revised version of the paper, reflecting the reviewer's feedback.

Weakenss 1. Wrong definitions for (generalized) additive models.

In our paper $f$ can be a function of any quantity. For the regression model, we could have $f($x$)=\mathbb{E}(Y|$X$=$x$),$ and for the logistic regression, it could be $f($x$)= {\rm logit} \Pr(Y=1|$X$=$x$)).$ More generally, in the generalized linear model, we could have $f($x$)=g(\mathbb{E}(Y|$X=x$))$ for a certain link function $g.$ In the revision, we added this explanation right after introducing $f.$

Weakness 2. Novelty in modeling and 3. Missing crucial mathematical background.

As far as I know, KAN is an extension of projection pursuit regression, takes the sum of univariate functions as input to a activation function, and is an interpretable model.

The differences between ANOVA-NODE and KAN are as follows.

2.1. KAN captures interactions through the activation function, whereas the functional ANOVA model explicitly represents interactions with specific interaction terms.

2.2. The $entmax({**x** - b \over \gamma })$ in ANOVA-NODE corresponds to the basis function $B(**x**)$ in the activation function of KAN. While KAN does not train the basis function, ANOVA-NODE can be considered to learn the basis by learning $b , \gamma$.

2.3. In KAN, the closed-form functional relationship between the output and components is not available, whereas in ANOVA-NODE, it can be determined.

2.4. While the KAN model does not guarantee identifiability, ANOVA-NODE satisfies identifiability, ensuring the stability of the interpretations it provides.

Therefore, these two models are different, and both are interpretable and useful. We left the comparison between these two models as a future work.

Weakness 3. Theory.

As mentioned in Line 178 (page 4), the sum-to-zero condition ensure the uniqueness of the functional ANOVA decomposition, thereby providing a stable interpretation. Note that, there are multiple conditions that ensure such uniqueness, and the sum-to-zero condition is just one of them. Also, theoretically, the global minimum does not change even with the sum-to-zero condition. However, since the model is trained by using a local minimum, the performance of the model can vary depending on the identifiability condition. Therefore, we do not know which condition provides the optimal results.

As suggested by the reviewer, we summarize the challenges in the proof of Theorem 3.3 and add the specific assumptions to the theorem. The proof for the case when the interaction order $d \geq 2$ is very complex to express using mathematical formulas, so it was omitted. The most important and challenging part of the proof of Theorem 3.3 is decomposing the ensemble function $f_{\mathcal{E}}$ into ANOVA-NODTs. Therefore, we added how the decomposition is performed using a toy example for $d=2$ in the proof of Theorem 3.3.

Weakness 5. Optimization & Experiments.

As pointed out by the reviewer, individual modeling with respect to each input feature may reduce prediction accuracy. However, we demonstrated through the experiments of prediction performance in the paper that considering up to second-order interactions is sufficient to achieve comparable prediction performance. Moreover, as far as I know, the target loss of many other models such as deep neural networks is also non-convex. The smoothness of ANOVA-NODE can be controlled by imposing constraints on the $\gamma$ parameter in Section 3.3 of the paper.

Weakness 6.

We did determine the threshold $K_{S}$ by empirical search. The experimental results of the prediction performance in ANOVA-NODE as $K_{S}$ increases are provided in Appendix D.1.

Theoretically, as $K_{S}$ increases, ANOVA-NODE can approximate any smooth function well, which can be found in Theorem 3.3. Also, choosing the optimal $K_{S}$ value requires research on the convergence rate, which will be addressed as a future work.

Finally, since we estimate each component of the function ANOVA model using an ensemble of $K_{S}$ NODTs, increasing $K_{S}$ does not affect the interpretability of ANOVA-NODE.

We appreciate your feedback on our paper and have made every effort to address your insightful questions. We have uploaded the revised version of the paper, reflecting the reviewer's feedback

Weakness 1. Reproducibility issue.

As suggested by the reviewers, we have included the source code of ANOVA-NODE in the Supplementary Material and we add the description of the data preprocessing in Appendix C.2 of the revision.

Weakness 2. Inconsistency between subsection title and Content.

There was originally an explanation of NODE, but due to the page limit of the paper, it was removed. However, the title was not updated accordingly. In the revision, we add the description of NODE in subsection 2.3.

Weakness 3. Need for Enhancements in Interpretability.

The results of importance score for second-order interactions and main effects on Calhousing dataset are provided in Appendix E.1.1 of the revision. As suggested by the reviewer, contour plots for the second-order interactions have been added to Appendix E.1.1 of the revision.

Question 1.

As suggested by the reviewer, we conducted additional experiments to compare the computational complexity of NAM, NBM, ANOVA-NODE, and NBM-NODE based on the runtime of each model and the number of weight parameters. Note that NBM-NODE is an extended model of ANOVA-NODE, proposed in Appendix G.2 of this paper, which improves the scalability of ANOVA-NODE by utilizing basis functions, similar to the approach used in NBM.

We consider NA$^{1}$M, which has 3 hidden layers with 16, 16, and 8 units; 10 basis DNNs for NB$^{1}$M, which has 3 hidden layers with 32, 16, and 16 units; 10 trees for each component in ANOVA-N$^{1}$ODE; and 10 basis functions in NBM-N$^{1}$ODE.

Table A.1. Results of computational complexity.

Table A.1 presents the number of weight parameters and runtime for each model when trained on a Abalone dataset with 100 epochs. Our computational environment consists of RTX 3090 and RTX 4090. We observed that our model has fewer weight parameters and takes less runtime compared to existing models (NAM and NBM).

Question 2.

As the reviewer understood, ANOVA-NODE learns the effects of all interaction orders simultaneously. When lower-order and higher-order interactions are trained sequentially, multicollinearity can obscure estimating the optimal model. That is, sequential learning would not be optimal which can be seen even in a linear model. For example, consider the model $Y=\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+\epsilon$, where there is mulicollinearity between $X_{1}$ and $X_{2}$ (i.e. ${\rm corr}(X_1,X_2) \neq 0$).

We estimate $\beta_{0},\beta_{1}$ and $\beta_{2}$ simultaneously by minimizing the loss, whose estimated coefficients are $\hat{\beta_{0}^{\*}}, \hat{\beta_{1}^{\*} }, \hat{\beta_{2}^{\*} }$. Altenatively, we estimate the $\beta_{0},\beta_{1},\beta_{2}$ sequentially by estimating $\beta_{0},\beta_{1}$ based on the model $Y=\beta_{0}+\beta_{1}X_{1}+\epsilon_{1}$ to have $\hat{\beta_{0}},\hat{\beta_{1}}$, and then estimating $\beta_{2}$ based on the model $Y-\hat{\beta_{0}}-\hat{\beta_{1}}X_{1} = \beta_{2}X_{2} + \epsilon_{2}$ to have $\hat{\beta_{2}}$.

In this case, $\{\hat{\beta_{0}^{\*}}, \hat{\beta_{1}^{\*}}, \hat{\beta_{2}^{\*}} \}$ and $\{\hat{\beta_{0}}, \hat{\beta_{1}}, \hat{\beta_{2}} \}$ are different and thus $\{\hat{\beta_{0}}, \hat{\beta_{1}}, \hat{\beta_{2}} \}$ is a non-optimal solution. As far as we understand, sequential learning is used in practice when simultaneous learning is computationally demanding.

Question 3.

We analyzed a wine dataset with 11 feature dimensions and a size of 4k, and applied NID to select the top 10 most important second-order interactions. Table A.2 shows the performance and stability scores of ANOVA-N$^{2}$ODE with and without NID. We observed an increase in stability when applying NID before ANOVA-NODE.

Table A.2. Results of prediction performance and stability score with and without NID.