A Novel Test for Additivity \\ in Supervised Ensemble Learners
Lucas Mentch and Giles Hooker \vspace{7mm
Additivity_arXiv.tex
\documentclass[12pt]{article}

\usepackage[margin=1.25in]{geometry}	
\usepackage{setspace}
\usepackage{natbib}
\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{url}
\usepackage{etoolbox}
\usepackage{bm}
\usepackage{amsmath,amssymb,amsthm}
\usepackage[linesnumbered,ruled,vlined]{algorithm2e}
\usepackage{float}
\usepackage{multirow}
\usepackage{colortbl}
\usepackage{enumerate}
 


\newcommand{\dataset}{{\cal D}}
\newcommand{\fracpartial}[2]{\frac{\partial #1}{\partial  #2}}
\newcommand{\tR}{\mathbb{R}}
\newcommand{\indist}{\overset{d}{\rightarrow}}     
\newcommand{\inprob}{\overset{P}{\rightarrow}}   
\newcommand{\as}{\overset{as}{\rightarrow}} 	    
\newcommand{\iid}{\stackrel{iid}{\sim}}		    
\newcommand{\um}{U_{n,k_n,m_n}}
\newcommand{\uk}{U_{n,k_n}}
\newcommand{\hatukm}{\hat{U}_{n,k_n,m_n}}
\newcommand{\honekn}{h_{1,k_n}}
\newcommand{\xred}{\boldsymbol{X}_R}
\newcommand{\xadd}{\boldsymbol{X}_A}
\newcommand{\xstar}{\mathbb{X}^*}
\newcommand{\predpoint}{\boldsymbol{x}^*}
\newcommand{\treefn}{T}
\newcommand{\featvect}{\boldsymbol{x}_i}
\newcommand{\predpoints}{\boldsymbol{x}_{\text{\tiny TEST}}}
\newcommand{\randparam}{\omega}
\newcommand{\randparamnum}{\mathcal{R}}
\newcommand{\featspace}{\mathcal{X}}
\newcommand{\uinc}{U_{n}^{(m)}}
\newcommand{\gred}{g^{(R)}}
\newcommand{\gredhat}{\hat{g}^{(R)}}
\newcommand{\nxone}{n_{\tilde{\boldsymbol{x}}}}
\newcommand{\nmc}{n_{MC}}
\newcommand{\rf}{r_{n,k_n,m_n}}
\newcommand{\rfstar}{r_{n,k_n,m_n}^{*}}
\newcommand{\randkern}{h_{k_n}^{(\omega_i)}}
\newcommand{\sigone}{\Sigma_{1,k_n}}
\newcommand{\sigk}{\Sigma_{k_n,k_n}}

\newtheorem{theorem}{Theorem}

\begin{document}



\title{A Novel Test for Additivity  \\ in Supervised Ensemble Learners}
\author{Lucas Mentch and Giles Hooker \vspace{7mm} \\
Department of Statistical Science \\ Cornell University \\ 
Ithaca, NY  14850   USA \\
email: \texttt{lkm54@cornell.edu} \\ \hspace{15.5mm} \texttt{gjh27@cornell.edu} \vspace{7mm} }

\maketitle

\begin{center}
\textbf{Abstract}
\end{center}

Additive models remain popular statistical tools due to their ease of interpretation and as a result, hypothesis tests for additivity have been developed to assess the appropriateness of these models.  However, as data grows in size and complexity, learning algorithms continue to gain popularity due to their exceptional predictive performance.  Due to the black-box nature of these learning methods, the increase in predictive power is assumed to come at the cost of interpretability and inference.  However, recent work suggests that many popular learning techniques, such as bagged trees and random forests, have desirable asymptotic properties which allow for formal statistical inference when base learners are built with proper subsamples.  This work extends hypothesis tests previously developed and demonstrates that by enforcing a grid structure on an appropriate test set, we may perform formal hypothesis tests for additivity among features.  We develop notions of total and partial additivity and demonstrate that both tests can be carried out at no additional computational cost.  We also suggest a new testing procedure based on random projections that allows for testing on larger grids, even when the grid size is larger than that of the training set.  Simulations and demonstrations on real data are provided.

\vspace*{.3in}

\noindent\textsc{Keywords}: {Bagging, Random Forest, Random Projection, Interaction, U-statistic}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%  INTRODUCTION  %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\doublespace
\section{Introduction}

Additive models were suggested by \cite{Friedman1981} and later generalized and made popular by \cite{HastieAdditiveModels}.  An underlying regression function $F\colon \mathcal{X} \mapsto \mathbb{R}$ is said to be additive if

\[
F(X_1, ..., X_d) = \sum_{i=1}^{d} F_{i}(X_i)
\]

for some functions $F_1, ..., F_d$.  If the regression function cannot be written as, or at least well-approximated by, a sum of univariate functions, then an interaction exists between at least two features.  Many methods have been developed to estimate the additive functions $F_1, ..., F_d$ including a method based on marginal integration by \cite{Linton1995}, a wavelet method suggested by \cite{Amato2001}, and a tree-based method by \cite{Lou2013}, but the most popular class of estimation methods are based on backfitting algorithms as found in \cite{Buja1989}, \cite{Opsomer1998,Opsomer1999}, and \cite{Mammen1999}.

The popularity of additive models lead naturally to hypothesis tests for additivity to asses whether an additive model is appropriate.  Versions of these lack-of-fit tests have been proposed by \cite{Barry1993}, \cite{Eubank1995}, \cite{Dette2001ANOVA}, \cite{Dette2001additivity}, \cite{Derbort2002}, and \cite{De2004}.  In some cases, such as in the work of \cite{FanAdditivity}, these procedures can also evaluate whether the additive functions belong to a particular parametric class.  Even when additive models are not used as the primary analytical tool, scientists often utilize tests of additivity to determine whether features (covariates) contribute additively to the response.  If no interactions are detected, then levels of one feature may be changed without affecting the contribution to the response of the others.  

Despite their usefulness, simple and intuitive tools like additive models can often fail to fully capture the signal hidden within modern complex data.  Learning algorithms like bagged trees and random forests, both introduced by \cite{bagging, randomforests}, are robust to a variety of regression functions and typically generate significantly more accurate predictions.  However, until recently, little was understood about the asymptotic behavior of these ensemble learners.  Recently, \cite{MentchTrees} showed that by using subsamples instead of full bootstrap samples to build individual trees, the predictions generated by these ensemble learners can be viewed as incomplete, infinite-order U-statistics and as such, are asymptotically normal.  This asymptotic normality allows for confidence intervals to accompany predictions and hypothesis tests to formally assess the significance of features.  Furthermore, by imposing structure on the subsamples used in the ensemble, variance estimates may be computed at no additional cost to the original ensemble.  \cite{WagerIJ} develop an infinitesimal jackknife estimate of variance and \cite{Wager2014} later showed that the conditions on subsample size imposed by \cite{MentchTrees} may be relaxed in the context of random forests for a specific class of tree-building methods.

This paper continues the trend of inference procedures for learning algorithms by developing formal hypothesis tests for additivity within the context of ensemble learners like bagged trees and random forests.  These tests allow practitioners to formally investigate the manner in which features contribute to the response when simpler tools like additive models are insufficient.  In Section \ref{sec:2} we propose a revised test for feature significance that imposes a grid structure on the test set and in Section \ref{sec:3} we demonstrate how this additional structure allows for tests of additivity.  In Section \ref{sec:4} we extend our procedure to the case where a large test grid is needed, so as to cover the case where $p>n$.  In Section \ref{sec:5} we provide simulations to investigate the power and $\alpha$-level of our tests and in Section \ref{sec:6} we apply our procedure to real data.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%  NEW TEST OF SIGNIFICANCE  %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{An alternative test for feature significance}
\label{sec:2}

We begin by recalling the key result from \cite{MentchTrees} which paves the way for statistical inference to be carried out for ensemble learners like bagged trees and random forests.  The authors demonstrate that the predictions from these ensembles may be viewed as incomplete, infinite-order U-statistics when the base learners, most often trees, are built with subsamples instead of full bootstrap samples.  The kernels of these U-statistics are functions that map from the subsample to the prediction via the trees and since the result holds in the infinite order case, the size of each subsample may grow with the size of the training set, thereby allowing trees to be grown to greater depths when more data is available.  Thus, given a set of subbagged trees or a subsampled random forest based on a training set of size $n$ in which $m$ trees are built, each with a subsample of size $k$, the prediction at a given point $\predpoint$ is asymptotically normal and the limiting variance depends only on the number of base learners relative to the size of the training set.  This central limit theorem and more discussion is provided in Appendix \ref{app:A}.

Given this asymptotic normality, \cite{MentchTrees} propose a hypothesis test to evaluate feature significance, which we briefly summarize here.  Suppose that the data consists of only two covariates $X_1$ and $X_2$ and that the response is observed according to $Y = F(X_1, X_2) + \epsilon$.  To test the significance of $X_2$, we begin with a test set $\predpoints$ consisting of $N$ points and build two subsampled ensembles $\hat{F}$ and $\hat{F}_1$ using the same subsamples, but where $\hat{F}$ is built using both $X_1$ and $X_2$ and $\hat{F}_1$ is built using only $X_1$.  Predictions at each point in $\predpoints$ are then made with each ensemble and by extension of the central limit theorem referenced above, the vector of differences in predictions has a multivariate normal distribution with mean $\mu$ and variance $\Sigma$.  Given consistent estimators of these parameters, $\hat{\mu}^{T} \hat{\Sigma}^{-1} \hat{\mu} \sim \chi_{N}^2$ can be used as a test statistic to formally test the hypotheses
\begin{align}
\label{hyp:Sig2feat} &H_0:  F(x_1, x_2) = F_1(x_1) \; \; \forall (x_1, x_2) \in \predpoints \\
&H_1:  F(x_1, x_2) \neq F_1(x_1) \; \; \mbox{for some } (x_1, x_2) \in \predpoints \; \; \mbox{for any } F_1. \notag
\end{align}

\noindent The testing procedure we propose here follows similarly, but we impose additional structure on the test set which will allow us to perform tests for additivity and avoid training an additional set of trees. 

\begin{figure}
  \centering
  \includegraphics[scale=0.4]{JASAlabels.png}
  \caption{\label{fig:grid} A grid of test points shown in the $X_1$$X_2$ plane below the response surface.}
\end{figure}

Define a grid consisting of $N$ total test points as in Figure \ref{fig:grid} with $N_1$ levels of $X_1$ and $N_2$ levels of $X_2$ so that the $(i,j)^{th}$ point in the grid has true value $F_{ij} = F(X_{1_i},X_{2_j})$ and predicted value $\hat{F}_{ij}$.  Further, let $V_F$ and $V_{\hat{F}}$ represent the vectorized versions of these true and predicted values so that $V_{F} = (F_{1,1}, ..., F_{1,N_2}, ..., F_{N_1,1}, ..., F_{N_1, N_2})$ and define

\[
\hat{f}_{i \cdot } = \frac{1}{N_2} \sum_{i=1}^{N_2} \hat{F}_{ij}
\]

\noindent as the average response at the $i^{th}$ level of $X_1$ across all grid levels of $X_2$.  For each point in the grid, define the difference in predictions $\hat{F}_{ij} - \hat{f}_{i \cdot}$ which we can write in vectorized form as $DV_{\hat{F}}$ for some $N \times N$ difference matrix $D$ of rank $N-N_1$.  Let $\Sigma$ denote the covariance of $V_F$ so that we can write $\Sigma_D = cov(DV_F) = D \Sigma D^T$.  We then estimate $\Sigma$ using the internal covariance estimation procedure proposed in \cite{MentchTrees} to calculate $\hat{\Sigma}$ so that $\hat{\Sigma}_{D} = D \hat{\Sigma} D^{T}$ forms a consistent estimate of $\Sigma_{D}$.  Then, by the theoretical results established in \cite{MentchTrees}, $(DV_{\hat{F}})^{T} \hat{\Sigma}^{-1}_{D} DV_{\hat{F}} \sim \chi_{N-N_1}^{2}$ and since we can equivalently write the hypotheses in (\ref{hyp:Sig2feat}) as 
\begin{align}
&H_0:  F_{ij} - f_{i \cdot} = 0 \; \; \forall (x_1,x_2) \in \predpoints \notag \\
&H_1:  F_{ij} - f_{i \cdot} \neq 0 \; \; \mbox{for some } (x_1, x_2) \in \predpoints \notag
\end{align}

\noindent we can use this as a test statistic.  This testing procedure is laid out in Algorithm \ref{algo:1}.  The variables $\nxone$ and $\nmc$ are covariance estimation parameters selected by the user as part of the internal covariance estimation procedure; details can be found in \cite{MentchTrees} and a summary is provided in Appendix \ref{app:A}.

\begin{algorithm}
\hspace{1mm} Select internal covariance estimation parameters $\nxone$ and $\nmc$ \\
\hspace{1mm} Perform internal covariance estimation procedure to build ensemble and estimate $V_F$ and $\Sigma$ \\
\hspace{1mm} Compute $\hat{f}_{\cdot j}$ and calculate $D$ and $\hat{\Sigma}_{D} = D \hat{\Sigma} D^{T}$ \\ \label{algostep:D}
\hspace{1mm} Compute test statistic $(DV_{\hat{F}})^{T} \hat{\Sigma}^{-1}_{D} DV_{\hat{F}}$ \label{algostep:tstat}

\caption{Test for Feature Significance}
\label{algo:1}
\end{algorithm}

Asymptotically, this test statistic has a $\chi_{N-N_1}^{2}$ distribution and thus can be compared to the $1-\alpha$ quantile as the critical value to achieve a test with type 1 error rate $\alpha$.  If the test statistic is larger than the critical value, we reject the null hypothesis and conclude that $X_2$ is significant.  

Note that this testing procedure readily extends to the more general case of $d$ features $X_1, ..., X_d$.  Let $\xred$ and $\xadd$ form a partion of the feature set $\{ X_1, ..., X_d \}$ so that $\xred$ and $\xadd$ are disjoint and $\xred \cup \xadd = \{ X_1, ..., X_d \}$; the set $\xred$ denotes the \emph{reduced} set of features and $\xadd$ represents the \emph{additional} features that we want to test for significance.  To test the hypotheses
\begin{align*}
&H_0:  F(\boldsymbol{x}_{R_i}, \boldsymbol{x}_{A_i}) = F_R(\boldsymbol{x}_{R_i}) \; \; \forall (\boldsymbol{x}_{R_i}, \boldsymbol{x}_{A_i}) \in \predpoints \\
&H_1:  F(\boldsymbol{x}_{R_i}, \boldsymbol{x}_{A_i}) \neq F_R(\boldsymbol{x}_{R_i}) \mbox{ for some } (\boldsymbol{x}_{R_i}, \boldsymbol{x}_{A_i}) \in \predpoints \mbox{ for any } F_R
\end{align*}

\noindent we simply replace $X_1$ with $\xred$ and $X_2$ with $\xadd$ in the above work by appropriately redefining levels of each feature set to form the grid of test points and repeat the procedure in Algorithm \ref{algo:1}.  Note that in this case, each grid point now corresponds to a level of a vector of features. 

It is also worth noting that \cite{MentchTrees} suggest comparing predictions generated with the full training set not only to those produced with the reduced set $\xred$, but also to those generated with $\xred$ and a permuted version of $\xadd$ in order to rule out the possibility that the ensemble is simply making use of additional noise.  The procedure we propose above avoids this problem by utilizing the projection matrix $D$.  



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%  INTERACTIONS PART 1  %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Tests for additivity}
\label{sec:3}

We now demonstrate how the grid structure imposed on the test set allows for tests of additivity.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%  Additivity 1  %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection*{Tests for total additivity}
Again assume that our training set consists of only two features $X_1$ and $X_2$ and that the response is observed according to $Y = F(X_1, X_2) + \epsilon$.  Tests for total additivity asses whether the underlying regression function $F$ is equal to, or at least well-approximated by, a sum of univariate functions so that in this 2-dimensional case, the hypotheses of interest are 
\begin{align}
\label{hyp:TotAdd1} &H_0:  \exists \; F_1, F_2 \mbox{  such that } F(x_1, x_2) = F_1(x_1) + F_2(x_2) \; \; \forall (x_1, x_2) \in \predpoints \\
&H_1:  F(x_1, x_2) \neq F_1(x_1) + F_2(x_2) \; \; \mbox{for some } (x_1, x_2) \in \predpoints \mbox{ for any } F_1, F_2. \notag
\end{align}

Again define a 2-dimensional grid of test points as in Figure \ref{fig:grid} so that each point in the grid has true value $F_{ij}$, predicted value $\hat{F}_{ij}$, and vectorized versions $V_F$ and $V_{\hat{F}}$.  Define $\bar{F}$ to be the mean of all predictions in the grid and define 

\[
\hat{f}_{i \cdot} = \frac{1}{N_2} \sum_{j=1}^{N_2} \hat{F}_{ij}   \hspace{25mm} \mbox{ and } \hspace{25mm}  \hat{f}_{\cdot j} = \frac{1}{N_1} \sum_{i=1}^{N_1} \hat{F}_{ij}
\]

\noindent as the mean prediction at the $i^{th}$ level of $X_1$ across all levels of $X_2$, and the mean prediction at the $j^{th}$ level of $X_2$ across all levels of $X_1$, respectively.  If the features are additive, i.e. under the null hypothesis, all points $(x_{1_i},x_{2_j})$ in the grid can be written as $F_{ij} = f_{i \cdot} + f_{\cdot j} - \mu$ where $\mu = \mathbb{E}\bar{F}$ is the true mean expected prediction across all points in the grid.  Thus, we may equivalently write the hypotheses in (\ref{hyp:TotAdd1}) as
\begin{align}
&H_0:  F_{ij} - f_{i \cdot} - f_{\cdot j} + \mu = 0 \; \; \mbox{for all } (x_i,x_j) \in \predpoints \notag \\
&H_1:  F_{ij} - f_{i \cdot} - f_{\cdot j} + \mu \neq 0 \; \; \mbox{for some } (x_1, x_2) \in \predpoints. \notag
\end{align}

The natural test statistic that arises is $\hat{F}_{ij} - \hat{f}_{i \cdot} - \hat{f}_{\cdot j} + \bar{F}$ which we can write in vectorized form as $D_2 V_{\hat{F}}$ for a difference matrix $D_2$.  Again, let $\Sigma$ denote the covariance of $V_{F}$ so that we can write $\Sigma_{D_2} = cov(D_2 V_{F}) = D_2 \Sigma D^{T}_{2}$.  Given $N_1$ and $N_2$ levels of $X_1$ and $X_2$ respectively, we estimate $P = 1 + (N_1 - 1) + (N_2 - 1)$ parameters and we can use  $(D_2 V_{\hat{F}})^{T} \hat{\Sigma}^{-1}_{D_2} D_2 V_{\hat{F}}\sim \chi_{N-P}^{2}$ as our test statistic.  This testing procedure for total additivity is identical to the procedure in Algorithm \ref{algo:1} but in Steps \ref{algostep:D} and \ref{algostep:tstat} we calculate the appropriate difference matrix and test statistic.


Note that this procedure naturally extends to the case of $d$ features $X_1, ..., X_d$.  To test hypotheses of the form 
\begin{align}
\label{hyp:FeatAdd} &H_0:  \exists \; F_1, ..., F_d \mbox{  s.t. } F(x_1, ..., x_d) = F_1(x_1) + \cdots + F_d(x_d) \; \; \forall (x_1, ...,x_d) \in \predpoints \\
&H_1:  F(x_1, ..., x_d) \neq F_1(x_1) + \cdots + F_d(x_d) \; \; \mbox{for some } (x_1, ..., x_d) \in \predpoints \mbox{ for any } F_1, ..., F_d \notag 
\end{align}

\noindent we require a $d$-dimensional grid of test points so that given $N_i$ levels of each feature $X_i$, our grid contains a total of $N = \prod_{i=1}^{d} N_i$ test points.  Further, define

\[
\hat{f}_{\cdots j \cdots} = \frac{1}{N_1 \cdots N_{p-1} N_{p+1} \cdots N_{d}} \sum_{i_1 = 1}^{N_1} \cdots \sum_{i_{p-1} = 1}^{N_{p-1}} \sum_{i_{p+1} = 1}^{N_{p+1}} \cdots \sum_{i_d = 1}^{N_d} \hat{F}_{i_1 \cdots j \cdots i_d}
\]

\noindent to be the average prediction over all points in the grid at the $j^{th}$ level of $X_p$.  As in the 2-dimensional case, we can rewrite the hypotheses in (\ref{hyp:FeatAdd}) as 
\begin{align}
&H_0:  F_{i_1 ... i_d} - f_{i_1 \cdots} - f_{\cdot i_2 \cdots} - \cdots - f_{\cdots i_d} + (d-1) \mu = 0 \; \; \mbox{for all } (x_{i_1}, ..., x_{i_d}) \in \predpoints \notag \\
&H_1:  F_{i_1 ... i_d} - f_{i_1 \cdots} - f_{\cdot i_2 \cdots} - \cdots - f_{\cdots i_d} + (d-1) \mu \neq 0 \; \; \mbox{for some } (x_{i_1}, ..., x_{i_d}) \in \predpoints \notag
\end{align}

\noindent and write $\hat{F}_{i_1 ... i_d} - \hat{f}_{i_1 \cdots} - \cdots - \hat{f}_{\cdots i_d} + (d-1) \bar{F}$ as $D_d V_{\hat{F}}$.  Again, we define $\Sigma$ to be the covariance of $V_{F}$ so that $\Sigma_{D_d} = cov(D_d V_{F}) = D_d \Sigma D^{T}_{d}$ and we can use  $(D_d V_{\hat{F}})^{T} \hat{\Sigma}^{-1}_{D_d} D_d V_{\hat{F}}\sim \chi_{N-P}^{2}$ as our test statistic, where $P = 1+(N_1-1) + \cdots + (N_d-1)$.

Furthermore, the additive functions need not be univariate.  Define a (disjoint) partition of the feature space $\boldsymbol{S}_1, ..., \boldsymbol{S}_m$ so that $\cup_{i=1}^{m}\boldsymbol{S}_i = \{X_1, ..., X_d \}$.  We can test hypotheses of the form
\begin{align}
 &H_0:  \exists \; F_1, ..., F_m \mbox{  such that } F(\boldsymbol{s}_1, ..., \boldsymbol{s}_m) = F_1(\boldsymbol{s}_1) + \cdots + F_m(\boldsymbol{s}_m) \; \; \forall (\boldsymbol{s}_1, ...,\boldsymbol{s}_m) \in \predpoints \notag \\
&H_1:  F(\boldsymbol{s}_1, ..., \boldsymbol{s}_m) \neq F_1(\boldsymbol{s}_1) + \cdots + F_m(\boldsymbol{s}_m) \; \; \mbox{for some } (\boldsymbol{s}_1, ..., \boldsymbol{s}_m) \in \predpoints \mbox{ for any } F_1, ..., F_m \notag
\end{align}

\noindent in exactly the same fashion by appropriately defining levels of an $m$-dimensional grid of test points.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%  Additivity 2  %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection*{Tests for partial additivity}

We now handle the case where we are interested in testing only whether a proper subset of features contribute additively to the response.  Suppose that our training set consists of three features $X_1, X_2$ and $X_3$ and we are interested in testing
\begin{align}
\label{hyp:PartAdd} &H_0:  \exists \; F_1, F_2 \mbox{  s.t. } F(x_1, x_2, x_3) = F_1(x_1,x_3) + F_2(x_2,x_3) \; \; \forall (x_1, x_2,x_3) \in \predpoints \\
&H_1:  F(x_1, x_2, x_3) \neq F_1(x_1,x_3) + F_2(x_2,x_3) \mbox{ for some } (x_1, x_2,x_3) \in \predpoints \mbox{ for any } F_1, F_2. \notag
\end{align}

\noindent Rejecting this null hypothesis means that an interaction exists between $X_1$ and $X_2$ but implies nothing about potential interactions between $X_1$ and $X_3$ or between $X_2$ and $X_3$.  \cite{Hooker2004} uses the size of the deviation of $F$ from partial additivity as a means of identifying the bivariate and higher-order interactions required to reconstruct some percentage of the variation in the values of $F$.  This is also referred to as the Sobol index for the $X_1$, $X_2$ interaction (\cite{Sobol2001}).  Define a 3-dimensional grid of test points with $N_1, N_2$ and $N_3$ levels of $X_1, X_2$ and $X_3$, respectively and continuing with the dot notation, define

\[
\hat{f}_{i \cdot k} = \frac{1}{N_2} \sum_{j=1}^{N_2} \hat{F}_{ijk}   \hspace{25mm} \mbox{ and } \hspace{25mm}  \hat{f}_{\cdot j k} = \frac{1}{N_1} \sum_{i=1}^{N_1} \hat{F}_{ijk}
\]

\noindent to be the average prediction over all levels of $X_2$ in the grid at the $i^{th}$ and $k^{th}$ levels of $X_1$ and $X_3$, and the average prediction over all levels of $X_1$ in the grid at the $j^{th}$ and $k^{th}$ levels of $X_2$ and $X_3$, respectively.  If there is no interaction between $X_1$ and $X_2$, then $F_{ijk} - f_{i \cdot k} - f_{\cdot j k} + f_{\cdot \cdot k} = 0$ at all levels $(x_{1_i},x_{2_j},x_{3_k})$ in the grid.  Thus, we can rewrite the hypotheses in (\ref{hyp:PartAdd}) as 
\begin{align*}
&H_0:  F_{ijk} - f_{i \cdot k} - f_{\cdot j k} + f_{\cdot \cdot k} = 0 \; \; \forall (x_1, x_2, x_3) \in \predpoints \\
&H_1: F_{ijk} - f_{i \cdot k} - f_{\cdot j k} + f_{\cdot \cdot k} \neq 0 \; \; \mbox{for some } (x_1, x_2, x_3) \in \predpoints 
\end{align*}

\noindent and use the empirical analogues of these parameters to conduct the testing procedure.  Once again, write $\hat{F}_{ijk} - \hat{f}_{i \cdot k} - \hat{f}_{\cdot j k} + \hat{f}_{\cdot \cdot k}$ as $D_3 V_{\hat{F}}$ and define $\Sigma$ to be the covariance of $V_{F}$ so that we can write $\Sigma_{D_3} = cov(D_3 V_{F}) = D_3 \Sigma D^{T}_{3}$ and use  $(D_3 V_{\hat{F}})^{T} \hat{\Sigma}^{-1}_{D_3} D_3 V_{\hat{F}}\sim \chi_{N-P}^{2}$ as our test statistic, where $N=N_1N_2N_3$ and since we must now account for two-way interactions, $P = 1+(N_1-1)+(N_2-1)+(N_3-1)+(N_1-1)(N_3-1)+(N_2-1)(N_3-1)$.  As was the case in testing for total additivity, the testing procedure is the same as in Algorithm \ref{algo:1} with the appropriate difference matrix and test statistic calculated in Steps \ref{algostep:D} and \ref{algostep:tstat}.

Note that we can perform this same test when our training set consists of $d$ features and we are interested in determining whether an interaction exists between $X_i$ and $X_j$.  Denote the set of all features except $X_i$ and $X_j$ as $\boldsymbol{X}_{-i,j}$ so that our hypotheses become
\begin{align*}
&H_0:  \exists \; F_i, F_j \mbox{  such that } F(X_1, ..., X_d) = F_i(X_i,\boldsymbol{X}_{-i,j}) + F_j(X_j,\boldsymbol{X}_{-i,j}) \; \; \forall (x_1, ..., x_d) \in \predpoints \\
&H_1:  F(X_1, ..., X_d) \neq F_i(X_i,\boldsymbol{X}_{-i,j}) + F_j(X_j,\boldsymbol{X}_{-i,j}) \mbox{ for some } (x_1, ..., x_d) \in \predpoints \mbox{ for any } F_i, F_j.
\end{align*}

By simply redefining the grid levels of $X_3$ to be levels of $\boldsymbol{X}_{-i,j}$, the testing procedure remains identical.  The same holds for the case where $X_i$ and $X_j$ are treated as vectors.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%  Additivity 2  %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection*{A general approach}

The above testing procedures were derived by choosing the model parameters that minimized the sum of squared error (SSE) with equal weight placed on each point in the test grid.  We could have instead selected $\hat{F}_1$ and $\hat{F}_2$ to minimize the weighted SSE

\[
WSSE = \sum_{i,j,k} w_{i,j,k} \left(F(x_{1_i},x_{2_j},x_{3_k}) - F_{1}(x_{1_i},x_{3_k}) - F_{2}(x_{2_j},x_{3_k}) \right)^2
\]

\noindent where the $w_{ijk}$ are specified by the user.  \cite{Hooker2007}, for example, recommends basing these weights on an approximation to the density of observations near the point $(x_{1_i},x_{2_j},x_{3_k})$.  This procedure takes the form of a weighted ANOVA.  In particular, define $\vec{F}$ to be the $N_1 N_3 + N_2 N_3$ vector concatenating the $\hat{F}_{1}(x_1,x_3)$ and $\hat{F}_{2}(x_2,x_3)$ and as in the previous sections let $V_{\hat{F}}$ be the vector containing the $\hat{F}_{ijk}$.  Further, let $Z$ be the $N \times N_1 N_3 + N_2 N_3$ matrix defined so that $Z \vec{F}$ produces the corresponding $\hat{F}_{1}(x_1,x_3) + \hat{F}_{2}(x_2,x_3)$ and let $W$ be a diagonal matrix containing the weights.  Then we can write

\[
WSSE= (V_{\hat{F}} - Z \vec{F})^T W (V_{\hat{F}} - Z \vec{F})
\]

\noindent and we know that the solution $\vec{F}$ that minimizes this weighted SSE is given by

\[
\vec{F} = (Z^T W Z)^{-1} Z^T W V_{\hat{F}}
\]

\noindent so that under the null hypothesis

\[
V_{\hat{F}} - Z \vec{F} = (I - Z(Z^T W Z)^{-1} Z^T W) V_{\hat{F}}
\]

\noindent has mean 0.  Further, letting $\Sigma$ denote the covariance of $V_{F}$, the variance of $V_{F} - Z \vec{F}$ is given by 

\[
C = (I - Z(Z^T W Z)^{-1} Z^T W) \Sigma (I - Z(Z^T W Z)^{-1} Z^T W)^T
\]

\noindent so that

\[
\left[ (I - Z(Z^T W Z)^{-1} Z^T W) V_{\hat{F}} \right]^T \hat{C}^{-1} \left[ (I - Z(Z^T W Z)^{-1} Z^T W) V_{\hat{F}} \right]
\]

\noindent has a $\chi^2_{N-P}$ distribution.  For equal weighting ($W$ given by the identity matrix), these calculations reduce to the averages employed above, and for the sake of simplicity we have restricted ourselves to this choice in the examples below.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%  Random Projections  %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Extension to Large Grids}
\label{sec:4}
Note that the above procedures require estimating a covariance matrix of size proportional to the number of points in the test grid.  In our experience, using relatively few trees to estimate the covariance parameters usually leads to an overestimate, which can thereby reduce the power in the testing procedures.  Thus, in situations where many grid points are required, either because certain features contain many levels or because we are interested in evaluating a more complex additive structure, it may be computationally infeasible to directly obtain a reasonable estimate.  To account for this potential large grid problem, we borrow from recently developed methods utilizing random projections.  

\subsection{Random Projections}
We present here a brief overview of random projection theory; for more background we refer the reader to \cite{Bingham2001}.  Let $X_{n \times p}$ represent a data matrix with sample size $n$ and presumably large number of covariates $p$, perhaps even $p>n$.  The idea is to project the data into a lower dimensional subspace via a random projection matrix $R_{p \times r}$ where $1 \leq r < n$ to form the projected data matrix $XR_{n \times r}$.  The Lemma provided by \cite{Johnson1984} states that an orthogonal projection of elements in a higher dimensional space into a lower dimensional space will approximately preserve the distances between projected elements.  Thus, selecting an orthogonal random projection matrix $R$ with relatively small projected dimension $r$ allows us to take initially high dimensional data and work in a more friendly lower dimensional space without a significant distortion of information.

Recently, \cite{pingli2014} proposed a new high dimensional testing procedure called RAPTT ({\bf Ra}ndom {\bf P}rojection {\bf T}-{\bf T}est) that represents an extension of Hotelling's classic $T^2$ test to the $p>n$ case.  The authors utilize random projections to test equality of multivariate means and demonstrate that this procedure performs well with reasonably high power even in a $p \gg n$ setting.  Given two multivariate samples $X_{n_1 \times p}$ and $Y_{n_2 \times p}$, the test statistic for the projected data is given by

\[
T_{R}^{2} = \frac{1}{n_{1}^{-1}+n_{2}^{-1}} (\bar{X}-\bar{Y})' R (R' S R)^{-1} R' (\bar{X}-\bar{Y})
\]

\noindent and the authors obtain a randomized test by defining

\[ \phi(T^{2}_{R}) = \left\{ 
  \begin{array}{l l}
    1 & \quad \text{if $\frac{n-r+1}{r} \frac{T^{2}_{R}}{n} > c_{\alpha}$}\\
    0 & \quad \text{otherwise}
  \end{array} \right.\]

\noindent where $n = n_1 + n_2 - 2$ and $c_{\alpha}$ is taken as the $1-\alpha$ quantile of the $F_{r,n - r +1}$ distribution.  Though this test is unbiased and consistent of exact size $\alpha$, it may have little power and lead to different conclusions based on the choice of $R$.  RAPTT is thus defined by averaging the p-values across $M$ such tests, utilizing $M$ random projections $R_1, ..., R_M$.  The p-value of the $i^{th}$ such test is given by

\[
\theta_i = 1 - F_{r,n-r+1} \left( \frac{n-r+1}{r} \frac{T^{2}_{R}}{n} \right)
\]

\noindent and RAPTT is defined by

\[ \bar{\phi} = \left\{ 
  \begin{array}{l l}
    1 & \quad \text{if $\bar{\theta} < u_{\alpha}$}\\
    0 & \quad \text{otherwise}
  \end{array} \right.\]
  
\noindent where $\bar{\theta} = \frac{1}{M} \sum_{i=1}^{M} \theta_i$ and $u_{\alpha}$ is chosen such that $P \left[ \bar{\theta} < u_{\alpha} \Big| H_0 \right] = \alpha$.  



%%%%%%%%%%%%%%%%%%%%%
%%%%  Application to Ensembles  %%%
%%%%%%%%%%%%%%%%%%%%%
\subsection{Application to Ensemble Tests}
We now demonstrate how this idea of testing using random projections can be applied to our procedures.  Suppose we have a training set of size $n$ and we build an ensemble consisting of $m$ trees, each of which is built with a subsample of size $k$, and that we are interested in predicting at $N$ total test points.  The test statistic for the original test of feature significance proposed in \cite{MentchTrees} can be written as $\hat{\mu}_{N}^{T} \hat{\Sigma}^{-1} \hat{\mu}_N$ where $\hat{\mu}$ is the vector of ensemble predictions, and $\hat{\Sigma}$ is a consistent estimate of the corresponding covariance matrix.  Based on the theory in \cite{MentchTrees}, this test statistic follows a $\chi_{N}^{2}$ distribution and we can employ a testing procedure based on random projections similar to the RAPTT procedure outlined above.  Given $m$ predictions at each of $N$ locations, we can think of our data as an $m \times N$ matrix so that for a set of random projection matrices $R_1, ..., R_M$ and reduced dimension $r < m,N$, we can write the projected test statistic as 

\begin{equation}
\label{tstat:origproj}
T_R = \hat{\mu}_{N}^{T} R (R^{T} \hat{\Sigma} R)^{-1} R^{T} \hat{\mu}_N \sim \chi^{2}_{r}.
\end{equation}

In the case where we perform the internal variance estimation procedure to simultaneously build the ensemble and estimate the variance parameters, we can think of our data matrix as having dimension $\nxone \times N$, in which case we choose $r < \nxone,N$.  

\noindent In the context of testing for additivity via the grid approach outlined above, the procedure remains largely the same.  We are already utilizing a difference matrix $D$ to project into the space of additive models, but so long as the elements of $R$ are independently generated continuous random variables, the overall projection has rank $r$ with probability 1.  Recall that the original test statistic is given by $(D V_{\hat{F}})^{T} \hat{\Sigma}^{-1}_{D} D V_{\hat{F}}$  where $\Sigma_{D} = cov(D V_{F}) = D \Sigma D^{T}$ and so our random projection test statistic becomes

\[
T_R = (D V_{\hat{F}})^{T} R (R^T \hat{\Sigma}_{D} R)^{-1} R^{T} (D V_{\hat{F}}) \sim \chi^{2}_{r}
\]

\noindent where now, $r < N-P$.  In both cases, the p-value for a single implementation of our random projection test based on the projection matrix $R_i$ is given by 

\[
\theta_i = 1 - \Phi^{2}_{r}(T_R)
\]

\noindent where $\Phi^{2}_{r}$ denotes the cdf of the $\chi^{2}_{r}$.  Thus, for $M$ repititions of this test, we can define our final test in exactly the same fashion as RAPTT, so that  

\[ \bar{\phi} = \left\{ 
  \begin{array}{l l}
    1 & \quad \text{if $\bar{\theta} < u_{\alpha}$}\\
    0 & \quad \text{otherwise}
  \end{array} \right.\]
  
where $\bar{\theta} = \frac{1}{M} \sum_{i=1}^{M} \theta_i$ and $u_{\alpha}$ is chosen such that $P \left[ \bar{\theta} < u_{\alpha} \Big| H_0 \right] = \alpha$.


\subsection{Choice of R and r}
\cite{pingli2014} note that the method of sampling projection matrices $R_i$ and the choice of reduced dimension $r$ are important considerations in this procedure as the power of the test can depend on both.  In our case, the covariance matrix parameters $\sigone$ and $\sigk$ are usually difficult to accurately estimate when a large number of grid points is required.  Thus, though the procedure is well defined for $1 \leq r < m$, this practical restriction necessitates a relatively small projected dimension $r$.  In our experience, we see a significant drop in power whenever our grids consist of more than approximately 30 points, so choosing $5 \leq r \leq 15$ should be reasonable and computationally feasible in most situations.  

Note also that because $r$ is small, little dependence remains between the resulting p-values $\theta_i$.  Under the null hypothesis, each p-value is uniformly distributed on $[0,1]$.  Furthermore, the mean of independent standard uniform random variables follows a Bates distribution, so the final cutoff $u_{\alpha}$ can be well approximated by the $\alpha$ quantile of this distribution when $r$ is small.  \cite{pingli2014} suggest simulating under the null hypothesis to estimate $u_{\alpha}$ but in our case, these simulations are difficult to perform since we cannot assume that the rows in our data matrix come from a multivariate normal distribution.  

The most appropriate choice of random projection matrices is less clear.  \cite{pingli2014} demonstrate that any semi-orthogonal matrix $R$ (i.e. $R^TR = I_r$, where $I_r$ is the $r \times r$ identity matrix) whose elements are continuous random variables with finite second moment will satisfy the necessary conditions to perform the tests based on random projections but consider both a uniform sample from the Haar distribution and an idea referred to as \emph{one permutation + one random projection} in developing RAPTT.  \cite{Bingham2001} note that the elements of $R$ are typically sampled from a normal distribution, but that they may be selected from simpler distributions, such as that suggested by \cite{Achlioptas2001}.  

\begin{algorithm}

\hspace{1mm} Compute difference matrix $D$ \\
\hspace{1mm} Select reduced dimension $r$ \\
\hspace{1mm} Generate random projection matrices $R_1, ..., R_M$ \\
\hspace{1mm} {\bf For} $i$ in 1 to $n_{\tilde{\bm{x}}}$

\hspace{8mm} Select initial fixed point $\tilde{\bm{x}}^{(i)}$

\hspace{8mm} {\bf For} $j$ in 1 to $n_{MC}$

\hspace{12mm} Select subsample $\mathcal{S}_{\tilde{\bm{x}}^{(i)},j}$ of size $k_n$ from training set that includes $\tilde{\bm{x}}^{(i)}$ \\
\hspace{12mm} Build tree using subsample $\mathcal{S}_{\tilde{\bm{x}}^{(i)},j}$ \\
\hspace{12mm} Use tree to predict at each of the $N$ grid points to obtain $\hat{V}_j$ \\
\hspace{12mm} Apply each projection to $\hat{V}$ to obtain $\hat{W}_{j,c} = (D\hat{V}_j)^{T}R_c$ \\
	
\hspace{8mm} {\bf End}

\hspace{8mm} Record average of $\hat{W}_{j,c}$ over $j$ for each projection 

\hspace{1mm} {\bf End}

\hspace{1mm} Compute variance of each of the $n_{\tilde{\bm{x}}}$ averages to estimate each $\zeta_{1,k_n}$ \\
\hspace{1mm} Compute variance of all predictions from each projection to estimate each $\zeta_{k_n,k_n}$ \\
\hspace{1mm} Compute mean of all predictions from each projection to estimate each $\theta_{k_n}$ \\
\hspace{1mm} Compute each p-value $\theta_1, ..., \theta_M$ by comparing to $\chi_{r}^{2}$ \\
\hspace{1mm} Record average p-value $\bar{\theta}$ and compare to Bates $\alpha$ quantile 
\caption{Internal Estimation Procedure}
\label{algo:largep}
\end{algorithm}

For our situation, we recommend simply generating random matrices whose elements are sampled from a standard normal, orthogonalizing via a process such as Gram-Schmidt, and selecting the appropriate submatrix.  Though such a method is perhaps not the most computationally efficient -- as noted by \cite{Bingham2001} -- it is straightforward and can be easily implemented in most software packages.


Our testing procedure based on random projections is laid out in Algorithm \ref{algo:largep}.  It is worth emphasizing again that though the original ensemble does not need to be built according to the internal estimation procedure, building in this fashion allows us to more easily use a small projected dimension $r$ and also allows for simultaneous estimation of the variance parameters.  It is also worth noting that the classic sample covariance $S$ is consistent and could potentially be used in place of the variance estimation procedures suggested by \cite{MentchTrees}, but typically does not produce an accurate estimate.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%  SIMULATIONS  %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Simulations}
\label{sec:5}

We now provide a short simulation study to investigate both the $\alpha$-level and power of our proposed testing procedures.  Suppose first that we have two features $X_1$ and $X_2$ and that our responses are generated according to 

\[
Y = X_1 + X_2 + \beta X_1X_2 + \epsilon
\]

\noindent where $\beta=0$ in assessing the $\alpha$-level and $\beta=1$ to evaluate power and $\epsilon \sim \mathcal{N}(0,0.05^2)$.  We repeat the test for total additivity on 1000 datasets when $\beta=0$ and 1000 where $\beta=1$, taking our empirical $\alpha$-level as the proportion of tests that incorrectly reject the null hypothesis (when $\beta=0$) and our estimate of power as the proportion of tests that correctly reject the null hypothesis (when $\beta=1$).  For reference, we also built 1000 linear regression models and used the traditional t-test to determine whether the interaction coefficient is significant on each dataset and recorded the empirical $\alpha$-level and power of this testing procedure.  This was repeated for training set sizes of 250, 500, and 1000 using subsample sizes of 30, 50, and 75 respectively and the results are shown in Table \ref{table:1}.   The test grid was selected as a $4 \times 4$ grid with levels 0.2, 0.4, 0.6, and 0.8.  In each case, our test for total additivity using a subbagged ensemble performed nearly exactly as well as the traditional t-test.


\begin{table}[t]
\label{table:1}
\begin{center}
\begin{tabular}{|r|c|c|c|}
\hline
\multicolumn{1}{|c|}{\textbf{Method}} & \textbf{$\boldsymbol{n}$}            & \textbf{$\boldsymbol{\alpha}$-level} & \textbf{Power} \\ \hline
Linear Model                          & \multirow{2}{*}{250}  & 0.056            & 1.000          \\ \cline{1-1} \cline{3-4} 
Subbagged Ensemble                            &                       & 0.065            & 0.954          \\ \hline
Linear Model                          & \multirow{2}{*}{500}  & 0.048            & 1.000          \\ \cline{1-1} \cline{3-4} 
Subbagged Ensemble                            &                       & 0.047            & 0.998          \\ \hline
Linear Model                          & \multirow{2}{*}{1000} & 0.046            & 1.000          \\ \cline{1-1} \cline{3-4} 
Subbagged Ensemble                            &                       & 0.020            & 0.999          \\ \hline
\end{tabular}
\end{center}
\caption{\label{table:1} Empirical $\alpha$-levels and power for the linear model example.}
\end{table}

We also selected a number of more complex regression functions to further investigate the $\alpha$-level and power.  Many of these underlying regression functions were studied in previous publications relating to tests for additivity, such as \cite{De2004} and \cite{Barry1993}.  Each estimate is the result of 1000 simulations with a sample size of 500, subsample size of 50, and internal covariance estimation parameters $\nxone=50$ and $\nmc=250$ with a $4 \times 4$ test grid (with levels 0.2, 0.4, 0.6, and 0.8) in the 2-dimensional tests for total additivity and a $3 \times 3 \times 3$ grid (with levels 0.3, 0.5, and 0.7) in the 3-dimensional tests for total and partial additivity.  In each case the features were selected uniformly at random from $[0,1]$ and the responses were generated according to $Y = F(\bm{X}) + \epsilon$ where $\epsilon \sim \mathcal{N}(0,0.05^2)$.  The results are shown in Table \ref{table:2}.  Note that even though the response in the first two models does not depend on $X_2$, this additional feature was still included in the training sets and the same test for total additivity was performed.  In each case, we see that our false rejection rate is very conservative and we also maintain high power.


%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%  Table  3  %%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{table}[t]
\label{table:2}
\begin{center}
\footnotesize
\begin{tabular}{|l|c|c|l|c|c|}
\hline
\multicolumn{1}{|c|}{\textbf{Model}} 		 & \textbf{Test} & \textbf{$\boldsymbol{\alpha}$-level} 		    & \multicolumn{1}{c|}{\textbf{Model}}        		& \textbf{Test} & \textbf{Power} \\ \hline
$y=x_1$                                			 & Total         & 0.000            						        	    & $y=x_1x_2$                                    				& Total         & 1.000          \\ \hline
$y=e^{x_1}$                          			 & Total         & 0.000        								    & $y=x_1x_2x_3$                                 				& Partial       & 0.605          \\ \hline
$y = e^{x_1} + \sin(\pi x_2)$            		 & Total         & 0.004        								    & $y =e^{5(x_1 + x_2)} / (1 + e^{5(x_1 + x_2)}) - 1$       & Total         & 0.948          \\ \hline
$y=x_1+x_2+x_3$                           		 & Total         & 0.002         							    & $y = 0.5(1 + \sin(2 \pi (x_1 + x_2)))$                            & Total         & 1.000          \\ \hline
$y = e^{x_1} + e^{x_2} + e^{x_3}$              & Total          & 0.004       								    & $y = 0.5(1 + \sin(2 \pi (x_1 + x_2 + x_3)))$                  & Partial       & 0.778          \\ \hline
$y=x_1x_3+x_2x_3$                        		 & Partial       & 0.000           							    & $y = 64 (x_1 x_2)^3 (1 - x_1 x_2)^3$ 				& Total          & 1.000          \\ \hline
$y=e^{x_1x_3}+e^{x_2x_3}$          		 & Partial       & 0.000           							    & $y = 64 (x_1 x_2 x_3)^3 (1 - x_1 x_2 x_3)^3$ 		& Partial       & 0.956          \\ \hline
\end{tabular}
\end{center}
\caption{\label{table:2} Empirical $\alpha$-level and power for a variety of underlying regression functions.}												
\end{table}


Finally, we repeated these simulations on the same regression functions, this time employing our testing procedures based on random projections.  In the 2-dimensional tests for total additivity, we use a $10 \times 10$ test grid so that $N=100$ and in the 3-dimensional tests for total additivity and the tests for partial additivity, we use a $5 \times 5 \times 5$ test grid for a total of 125 test points.  The results are shown in Table \ref{table:3}.  Note that in these tests, we maintain a reasonable type 1 error rate but have significantly more power due to the finer resolution of the test grid.  


\begin{table}[h]
\label{table:3}
\begin{center}
\footnotesize
\begin{tabular}{|l|c|c|l|c|c|}
\hline
\multicolumn{1}{|c|}{\textbf{Model}} 		 & \textbf{Test} & \textbf{$\boldsymbol{\alpha}$-level} 		    & \multicolumn{1}{c|}{\textbf{Model}}        		& \textbf{Test} & \textbf{Power} \\ \hline
$y=x_1$                                			 & Total         & 0.000            						        	    & $y=x_1x_2$                                    				& Total         & 1.000          \\ \hline
$y=e^{x_1}$                          			 & Total         & 0.000        								    & $y=x_1x_2x_3$                                 				& Partial       & 0.998          \\ \hline
$y = e^{x_1} + \sin(\pi x_2)$            		 & Total         & 0.059        								    & $y =e^{5(x_1 + x_2)} / (1 + e^{5(x_1 + x_2)}) - 1$       & Total         & 0.999          \\ \hline
$y=x_1+x_2+x_3$                           		 & Total         & 0.001         							    & $y = 0.5(1 + \sin(2 \pi (x_1 + x_2)))$                            & Total         & 1.000          \\ \hline
$y = e^{x_1} + e^{x_2} + e^{x_3}$              & Total          & 0.012       								    & $y = 0.5(1 + \sin(2 \pi (x_1 + x_2 + x_3)))$                  & Partial       & 0.959          \\ \hline
$y=x_1x_3+x_2x_3$                        		 & Partial       & 0.008           							    & $y = 64 (x_1 x_2)^3 (1 - x_1 x_2)^3$ 				& Total          & 1.000          \\ \hline
$y=e^{x_1x_3}+e^{x_2x_3}$          		 & Partial       & 0.066           							    & $y = 64 (x_1 x_2 x_3)^3 (1 - x_1 x_2 x_3)^3$ 		& Partial       & 1.000          \\ \hline
\end{tabular}
\end{center}
\caption{\label{table:3} Empirical $\alpha$-level and power for the testing procedure based on random projections.}												
\end{table}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%  REAL DATA  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Real data}
\label{sec:6}

We now demonstrate our testing procedures on a dataset provided by a team of ornithologists at the Cornell University Lab of Ornithology as part of the \emph{ebird} project described in \cite{ebird}.  This ongoing project relies on citizens to submit reports of bird observations and from these reports, researchers are able to monitor things like migration patterns and species abundance.  The dataset we were provided was originally compiled in order to determine how pollution levels affect the change in Wood Thursh population.  The data consists of 3 pollutant features, mercury deposition ($md$), acid deposition ($ad$), and soil PH level ($sph$) as well as 2 non-pollutant features, elevation ($elev$) and abundance ($ab$).  We begin our analysis by testing whether the pollutant and non-pollutant features are additive:
\begin{align}
&H_0:  F(md,ad,sph,elev,ab) = F_P(md,ad,sph) + F_{NP}(elev,ab). \label{hyp:ebird1}
\end{align}

We performed a test for total additivity using 4 levels of each feature set, the 0.20, 0.40, 0.60, and 0.80 quantiles of each feature, for a total of 16 test points.  Our test statistic was 52.30, larger than the critical value, the 0.95 quantile of the $\chi^{2}_{9}$, of 16.92 so we reject the null hypothesis in (\ref{hyp:ebird1}) and conclude that an interaction exists between the pollutant and non-pollutant features.  This result was confirmed by our random projection test, which consisted of 1000 random projections to a dimension of $r=5$ using a $10 \times 10$ test grid.  In this case, the final averaged p-value was only 0.0043, far below the critical value of 0.485.  

Next, we investigated how the pollutants contributed to the response.  Based on preliminary investigations, ebird researchers suspected an interaction between mercury and acid deposition ($md$ and $ad$) but were unsure of the relationship between soil PH ($sph$) and $md$ and $ad$.  In performing these tests for partial additivity, our test grid consisted of 3 points for each feature set, the 0.30, 0.50, and 0.70 quantiles of each feature for a total of 27 test points and a critical value, the 0.95 quantile of the $\chi_{12}^{2}$, of 21.03.  Our test for partial additivity between $md$ and $ad$,
\begin{align}
&H_0:  F(md,ad,sph,elev,ab) = F_1(md,sph,elev,ab) + F_2(ad,sph,elev,ab). \notag
\end{align}

\noindent yielded a significant result with a test statistic of 41.00 so our test supports the belief that an interaction exists between $md$ and $ad$.   Again, this result was supported by our random projection test, which consisted of 1000 random projections to a dimension of $r=5$ using a $5 \times 5 \times 5$ test grid, for a total of 125 test points.  The final averaged p-value was only 0.0064, far below the critical value of 0.485.  
 
 
Our test for partial additivity between $sph$ and the vector $(md,ad)$ 
 \begin{align}
&H_0:  F(md,ad,sph,elev,ab) = F_1(md,ad,elev,ab) + F_2(sph,elev,ab). \notag
\end{align}
 
\noindent yielded a test statistic of 36.43, above the critical value of 21.03, so once again we reject the null hypothesis and conclude that an interaction exists between $sph$ and $(md,ad)$.  This result was again supported by the random projection test based on 1000 random projections to a dimension of $r=5$ using a $5 \times 5 \times 5$ test grid.  We find a final averaged p-value of 0.225, which, though larger than in the previous tests, is still far below the critical value of 0.485.  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%  DISCUSSION  %%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Discussion}

This work continues the recent trend of developing formal statistical inference procedures within the context of learning algorithms.  The desirable asymptotic properties of subsampled ensembles explored by \cite{MentchTrees} and \cite{Wager2014} have begun to shed light on the black-box learning algorithms and the tests for additivity developed here further demonstrate that traditional scientific and statistical questions need not be seen as a sacrifice of less interpretable machine learning procedures. 

An important aspect of our testing procedures is the definition of the test grid.  In our applications we chose a relatively small grid in order achieve reasonably accurate covariance estimates while still maintaining an acceptable amount of power.  When larger testing grids are required, we suggest utilizing random projections in order to maintain high power.  Selecting levels of each feature away from the boundary of the feature space is also important as trees are known to exhibit bias near the edge of the feature space.

In the preceding work, all of our tests were carried out using subbagged ensembles, but the theory established in \cite{MentchTrees} and \cite{Wager2014} also allows for subsampled random forests or any other ensemble learner conforming to the regularity conditions.  However, the predictive improvement often seen with random forests is generally attributed to the increased independence between trees in the ensemble.  Since our base learners are built with subsamples instead of full bootstrap samples and thus are less dependent, we expect the predictive improvement often seen with random forests to be less significant.  

It is also worth noting that the particular additive forms for which we developed testing procedures were selected only because of their scientific utility.  Other similar testing procedures could also be developed in the same manner by establishing the appropriate model parameters from an ANOVA set-up and defining the difference matrix $D$ accordingly.  These methods can also be extended to provide formal statistical guarantees for the screening procedures described in \cite{Hooker2004}.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%  BIBLIOGRAPHY  %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliographystyle{apalike}
\begin{thebibliography}{}

\bibitem[Achlioptas, 2001]{Achlioptas2001}
Achlioptas, D. (2001).
\newblock Database-friendly random projections.
\newblock In {\em Proceedings of the twentieth ACM SIGMOD-SIGACT-SIGART
  symposium on Principles of database systems}, pages 274--281. ACM.

\bibitem[Amato and Antoniadis, 2001]{Amato2001}
Amato, U. and Antoniadis, A. (2001).
\newblock Adaptive wavelet series estimation in separable nonparametric
  regression models.
\newblock {\em Statistics and Computing}, 11(4):373--394.

\bibitem[Barry, 1993]{Barry1993}
Barry, D. (1993).
\newblock Testing for additivity of a regression function.
\newblock {\em The Annals of Statistics}, pages 235--254.

\bibitem[Bingham and Mannila, 2001]{Bingham2001}
Bingham, E. and Mannila, H. (2001).
\newblock Random projection in dimensionality reduction: applications to image
  and text data.
\newblock In {\em Proceedings of the seventh ACM SIGKDD international
  conference on Knowledge discovery and data mining}, pages 245--250. ACM.

\bibitem[Breiman, 1996]{bagging}
Breiman, L. (1996).
\newblock {B}agging predictors.
\newblock {\em Machine Learning}, 24:123--140.

\bibitem[Breiman, 2001]{randomforests}
Breiman, L. (2001).
\newblock {R}andom {F}orests.
\newblock {\em Machine Learning}, 45:5--32.

\bibitem[Buja et~al., 1989]{Buja1989}
Buja, A., Hastie, T., and Tibshirani, R. (1989).
\newblock Linear smoothers and additive models.
\newblock {\em The Annals of Statistics}, pages 453--510.

\bibitem[De~Canditiis and Sapatinas, 2004]{De2004}
De~Canditiis, D. and Sapatinas, T. (2004).
\newblock Testing for additivity and joint effects in multivariate
  nonparametric regression using fourier and wavelet methods.
\newblock {\em Statistics and Computing}, 14(3):235--249.

\bibitem[Derbort et~al., 2002]{Derbort2002}
Derbort, S., Dette, H., and Munk, A. (2002).
\newblock A test for additivity in nonparametric regression.
\newblock {\em Annals of the Institute of Statistical Mathematics},
  54(1):60--82.

\bibitem[Dette and Derbort, 2001]{Dette2001ANOVA}
Dette, H. and Derbort, S. (2001).
\newblock Analysis of variance in nonparametric regression models.
\newblock {\em Journal of multivariate analysis}, 76(1):110--137.

\bibitem[Dette et~al., 2001]{Dette2001additivity}
Dette, H., Wilkau, C. V. L.~U., et~al. (2001).
\newblock Testing additivity by kernel-based methods-what is a reasonable test?
\newblock {\em Bernoulli}, 7(4):669--697.

\bibitem[Eubank et~al., 1995]{Eubank1995}
Eubank, R., Hart, J.~D., Simpson, D., Stefanski, L.~A., et~al. (1995).
\newblock Testing for additivity in nonparametric regression.
\newblock {\em The Annals of Statistics}, 23(6):1896--1920.

\bibitem[Fan and Jiang, 2005]{FanAdditivity}
Fan, J. and Jiang, J. (2005).
\newblock Nonparametric inferences for additive models.
\newblock {\em Journal of the {A}merican {S}tatistical {A}ssociation},
  100(471):890--907.

\bibitem[Friedman and Stuetzle, 1981]{Friedman1981}
Friedman, J.~H. and Stuetzle, W. (1981).
\newblock Projection pursuit regression.
\newblock {\em Journal of the American statistical Association},
  76(376):817--823.

\bibitem[Hastie and Tibshirani, 1990]{HastieAdditiveModels}
Hastie, T.~J. and Tibshirani, R.~J. (1990).
\newblock {\em Generalized {A}dditive {M}odels}, volume~43.
\newblock CRC Press.

\bibitem[Hooker, 2004]{Hooker2004}
Hooker, G. (2004).
\newblock Discovering additive structure in black box functions.
\newblock In {\em Proceedings of the tenth ACM SIGKDD international conference
  on Knowledge discovery and data mining}, pages 575--580. ACM.

\bibitem[Hooker, 2007]{Hooker2007}
Hooker, G. (2007).
\newblock Generalized functional anova diagnostics for high-dimensional
  functions of dependent variables.
\newblock {\em Journal of Computational and Graphical Statistics}, 16(3).

\bibitem[Johnson and Lindenstrauss, 1984]{Johnson1984}
Johnson, W.~B. and Lindenstrauss, J. (1984).
\newblock Extensions of lipschitz mappings into a hilbert space.
\newblock {\em Contemporary mathematics}, 26(189-206):1.

\bibitem[Linton and Nielsen, 1995]{Linton1995}
Linton, O. and Nielsen, J.~P. (1995).
\newblock A kernel method of estimating structured nonparametric regression
  based on marginal integration.
\newblock {\em Biometrika}, pages 93--100.

\bibitem[Lou et~al., 2013]{Lou2013}
Lou, Y., Caruana, R., Gehrke, J., and Hooker, G. (2013).
\newblock Accurate intelligible models with pairwise interactions.
\newblock In {\em Proceedings of the 19th ACM SIGKDD international conference
  on Knowledge discovery and data mining}, pages 623--631. ACM.

\bibitem[Mammen et~al., 1999]{Mammen1999}
Mammen, E., Linton, O., Nielsen, J., et~al. (1999).
\newblock The existence and asymptotic properties of a backfitting projection
  algorithm under weak conditions.
\newblock {\em The Annals of Statistics}, 27(5):1443--1490.

\bibitem[Mentch and Hooker, 2014]{MentchTrees}
Mentch, L. and Hooker, G. (2014).
\newblock {E}nsemble {T}rees and {CLT}s: {S}tatistical {I}nference for
  {S}upervised {L}earning.
\newblock arXiv:1404.6473 [stat.ML].

\bibitem[Opsomer and Ruppert, 1998]{Opsomer1998}
Opsomer, J.~D. and Ruppert, D. (1998).
\newblock A fully automated bandwidth selection method for fitting additive
  models.
\newblock {\em Journal of the American Statistical Association},
  93(442):605--619.

\bibitem[Opsomer and Ruppert, 1999]{Opsomer1999}
Opsomer, J.~D. and Ruppert, D. (1999).
\newblock A root-n consistent backfitting estimator for semiparametric additive
  modeling.
\newblock {\em Journal of Computational and Graphical Statistics},
  8(4):715--732.

\bibitem[Sobol, 2001]{Sobol2001}
Sobol, I.~M. (2001).
\newblock Global sensitivity indices for nonlinear mathematical models and
  their monte carlo estimates.
\newblock {\em Mathematics and computers in simulation}, 55(1-3):271--280.

\bibitem[Srivastava et~al., 2014]{pingli2014}
Srivastava, R., Li, P., and Ruppert, D. (2014).
\newblock {RAPTT}: {A}n {E}xact {T}wo-{S}ample {T}est in {H}igh {D}imensions
  {U}sing {R}andom {P}rojections.
\newblock arXiv:1405.1792 [stat.ME].

\bibitem[Sullivan et~al., 2009]{ebird}
Sullivan, B.~L., Wood, C.~L., Iliff, M.~J., Bonney, R.~E., Fink, D., and
  Kelling, S. (2009).
\newblock ebird: A citizen-based bird observation network in the biological
  sciences.
\newblock {\em Biological Conservation}, 142(10):2282--2292.

\bibitem[Wager, 2014]{Wager2014}
Wager, S. (2014).
\newblock {A}symptotic {T}heory for {R}andom {F}orests.
\newblock arXiv:1405.0352 [math.ST].

\bibitem[Wager et~al., 2014]{WagerIJ}
Wager, S., Hastie, T., and Efron, B. (2014).
\newblock Confidence intervals for random forests: The jackknife and the
  infinitesimal jackknife.
\newblock {\em Journal of Machine Learning Research}, 15:1625--1651.

\end{thebibliography}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%  APPENDIX  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Appendix}
\begin{appendix}
\section{Internal Estimation Procedure}
\label{app:A}

We provide here a review of the internal estimation procedure suggested in \cite{MentchTrees}.  We begin by providing one of the main results of this work which establishes asymptotic normality and allows for formalized inference.
\begin{theorem}
\label{subbaggingthm}
Let $X_1, X_2, ... \iid F_X$ and let $U_{n,k_n,m_n}$ be an incomplete, infinite order U-statistic with kernel $h_{k_n}$ and $\theta_{k_n} = \mathbb{E}h_{k_n}(X_1, ..., X_{k_n})$ such that $\mathbb{E}h_{k_n}^{2}(X_1, ..., X_{k_n}) < \infty$ for all $n$.  Suppose that $\lim \frac{k_n}{\sqrt{n}} = 0$ and let $\lim \frac{n}{m_n} = \alpha$.  Then as long as $\lim k_{n}^{2} \zeta_{1,k_n} \neq 0$,
\begin{enumerate}[(i)]
\item if $\alpha = 0$, then $\frac{\sqrt{n}(U_{n,k_n,m_n} - \theta_{k_n})}{\sqrt{k_{n}^{2} \zeta_{1,k_n}}} \indist \mathcal{N}(0,1)$.
\item if $0 < \alpha < \infty$, then $\frac{\sqrt{m_n}(U_{n,k_n,m_n} - \theta_{k_n})}{\sqrt{\frac{k_{n}^{2}}{\alpha} \zeta_{1,k_n} + \zeta_{k_n,k_n}}} \indist \mathcal{N}(0,1)$.
\item if $\alpha = \infty$, then $\frac{\sqrt{m_n}(U_{n,k_n,m_n} - \theta_{k_n})}{\sqrt{\zeta_{k_n,k_n}}} \indist \mathcal{N}(0,1)$.
\end{enumerate}
\label{subbaggingthm}
\end{theorem}

Thus, provided that the subsample size $k$ is on the order of $\sqrt{n}$, the predictions follow one of the three limiting distributions given above, and the limiting variance depends only on the ratio of the training set size to the number of base learners in the ensemble.  The center of the distribution is estimated with the standard (point) prediction generated by the ensemble and the parameter $\alpha$ can be approximated by simply plugging in $n/{m_n}$ but the variance parameters $\zeta_{1,k_n}$ and $\zeta_{k_n,k_n}$ (or the multivariate analogues $\Sigma_{1,k_n}$ and $\Sigma_{k_n,k_n}$ used in the testing procedures) still need to be estimated.  For a general $c$, these variance parameters are defined as 
\[
\zeta_{c,k_n} = var \Big( \mathbb{E} \big( h_{k_n}(X_1, ..., X_{k_n})\; | \; X_1 = x_1, ..., X_c=x_c \big) \Big)
\]

\noindent so that 

\[
\zeta_{1,k_n} = var \Big( \mathbb{E} \big( h_{k_n}(X_1, ..., X_{k_n})\; | \; X_1 = x_1 \big) \Big)
\]

\noindent and 

\[
\zeta_{k_n,k_n} = var \big( h_{k_n}(X_1, ..., X_{k_n}) \big).
\]

\noindent The procedure to estimate $\zeta_{1,k_n}$ is provided below.  The procedure to estimate $\Sigma_{1,k_n}$ is nearly identical, but use the ensemble to predict at all points in step 6, and instead of taking the variance of the vector of $\nxone$ averages, we record the covariance of $\nxone$ rows of averages.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%  Algorithm to estimate Zeta_1  %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{algorithm}

\hspace{1mm} {\bf For} $i$ in 1 to $n_{\tilde{\bm{x}}}$

\hspace{8mm} Select initial fixed point $\tilde{\bm{x}}^{(i)}$

\hspace{8mm} {\bf For} $j$ in 1 to $n_{MC}$

\hspace{12mm} Select subsample $\mathcal{S}_{\tilde{\bm{x}}^{(i)},j}$ of size $k_n$ from training set that includes $\tilde{\bm{x}}^{(i)}$ \\
\hspace{12mm} Build tree using subsample $\mathcal{S}_{\tilde{\bm{x}}^{(i)},j}$ \\
\hspace{12mm} Use tree to predict at $\predpoint$ \\

\hspace{8mm} {\bf End}

\hspace{8mm} Record average of the $n_{MC}$ predictions \\

\hspace{1mm} {\bf End}

\hspace{1mm} Compute the variance of the $n_{\tilde{\bm{x}}}$ averages
\caption{$\zeta_{1,k_n}$ Estimation Procedure}
\label{algo:zeta1}
\end{algorithm}

To estimate $\zeta_{k_n,k_n}$ we can simply generate a large number of subsamples and record the variance (or covariance matrix) of the predictions.  Alternatively, to be more computationally efficient, we can record all predictions generated in estimating $\zeta_{1,k_n}$ and record the variance (covariance) amongst these as our estimate.

\begin{algorithm}

\hspace{1mm} {\bf For} $i$ in 1 to $n_{\tilde{\bm{x}}}$

\hspace{8mm} Select initial fixed point $\tilde{\bm{x}}^{(i)}$

\hspace{8mm} {\bf For} $j$ in 1 to $n_{MC}$

\hspace{12mm} Select subsample $\mathcal{S}_{\tilde{\bm{x}}^{(i)},j}$ of size $k_n$ from training set that includes $\tilde{\bm{x}}^{(i)}$ \\
\hspace{12mm} Build tree using subsample $\mathcal{S}_{\tilde{\bm{x}}^{(i)},j}$ \\
\hspace{12mm} Use tree to predict at $\predpoint$ and record prediction \\
	
\hspace{8mm} {\bf End}

\hspace{8mm} Record average of the $n_{MC}$ predictions 

\hspace{1mm} {\bf End}

\hspace{1mm} Compute the variance of the $n_{\tilde{\bm{x}}}$ averages to estimate $\zeta_{1,k_n}$ \\
\hspace{1mm} Compute the variance of all predictions to estimate $\zeta_{k_n,k_n}$ \\
\hspace{1mm} Compute the mean of all predictions to estimate $\theta_{k_n}$ \\
\caption{Internal Estimation Procedure}
\label{algo:internal}
\end{algorithm}


These variance estimation procedures can be performed after the original ensemble is computed, which is referred to as \emph{external} estimation.  However, this means that all of the predictions generated in an effort to estimate $\zeta_{1,k_n}$ simply go to waste.  To perform the \emph{internal} variance estimation procedure, we record all predictions used in Algorithm \ref{algo:zeta1}, using the average of all predictions as our point estimate and the variance across all predictions as our estimate of $\zeta_{k_n,k_n}$.  The algorithm to perform this internal estimation is provided in Algorithm \ref{algo:internal}.

\end{appendix}

\end{document}
1406.1845.extracted_bib
[Achlioptas, 2001]{Achlioptas2001 Achlioptas, D. (2001). \newblock Database-friendly random projections. \newblock In {\em Proceedings of the twentieth ACM SIGMOD-SIGACT-SIGART symposium on Principles of database systems, pages 274--281. ACM.
1406.1845.extracted_bib
[Amato and Antoniadis, 2001]{Amato2001 Amato, U. and Antoniadis, A. (2001). \newblock Adaptive wavelet series estimation in separable nonparametric regression models. \newblock {\em Statistics and Computing, 11(4):373--394.
1406.1845.extracted_bib
[Barry, 1993]{Barry1993 Barry, D. (1993). \newblock Testing for additivity of a regression function. \newblock {\em The Annals of Statistics, pages 235--254.
1406.1845.extracted_bib
[Bingham and Mannila, 2001]{Bingham2001 Bingham, E. and Mannila, H. (2001). \newblock Random projection in dimensionality reduction: applications to image and text data. \newblock In {\em Proceedings of the seventh ACM SIGKDD international conference on Knowledge discovery and data mining, pages 245--250. ACM.
1406.1845.extracted_bib
[Breiman, 1996]{bagging Breiman, L. (1996). \newblock {Bagging predictors. \newblock {\em Machine Learning, 24:123--140.
1406.1845.extracted_bib
[Breiman, 2001]{randomforests Breiman, L. (2001). \newblock {Random {Forests. \newblock {\em Machine Learning, 45:5--32.
1406.1845.extracted_bib
[Buja et~al., 1989]{Buja1989 Buja, A., Hastie, T., and Tibshirani, R. (1989). \newblock Linear smoothers and additive models. \newblock {\em The Annals of Statistics, pages 453--510.
1406.1845.extracted_bib
[De~Canditiis and Sapatinas, 2004]{De2004 De~Canditiis, D. and Sapatinas, T. (2004). \newblock Testing for additivity and joint effects in multivariate nonparametric regression using fourier and wavelet methods. \newblock {\em Statistics and Computing, 14(3):235--249.
1406.1845.extracted_bib
[Derbort et~al., 2002]{Derbort2002 Derbort, S., Dette, H., and Munk, A. (2002). \newblock A test for additivity in nonparametric regression. \newblock {\em Annals of the Institute of Statistical Mathematics, 54(1):60--82.
1406.1845.extracted_bib
[Dette and Derbort, 2001]{Dette2001ANOVA Dette, H. and Derbort, S. (2001). \newblock Analysis of variance in nonparametric regression models. \newblock {\em Journal of multivariate analysis, 76(1):110--137.
1406.1845.extracted_bib
[Dette et~al., 2001]{Dette2001additivity Dette, H., Wilkau, C. V. L.~U., et~al. (2001). \newblock Testing additivity by kernel-based methods-what is a reasonable test? \newblock {\em Bernoulli, 7(4):669--697.
1406.1845.extracted_bib
[Eubank et~al., 1995]{Eubank1995 Eubank, R., Hart, J.~D., Simpson, D., Stefanski, L.~A., et~al. (1995). \newblock Testing for additivity in nonparametric regression. \newblock {\em The Annals of Statistics, 23(6):1896--1920.
1406.1845.extracted_bib
[Fan and Jiang, 2005]{FanAdditivity Fan, J. and Jiang, J. (2005). \newblock Nonparametric inferences for additive models. \newblock {\em Journal of the {American {Statistical {Association, 100(471):890--907.
1406.1845.extracted_bib
[Friedman and Stuetzle, 1981]{Friedman1981 Friedman, J.~H. and Stuetzle, W. (1981). \newblock Projection pursuit regression. \newblock {\em Journal of the American statistical Association, 76(376):817--823.
1406.1845.extracted_bib
[Hastie and Tibshirani, 1990]{HastieAdditiveModels Hastie, T.~J. and Tibshirani, R.~J. (1990). \newblock {\em Generalized {Additive {Models, volume~43. \newblock CRC Press.
1406.1845.extracted_bib
[Hooker, 2004]{Hooker2004 Hooker, G. (2004). \newblock Discovering additive structure in black box functions. \newblock In {\em Proceedings of the tenth ACM SIGKDD international conference on Knowledge discovery and data mining, pages 575--580. ACM.
1406.1845.extracted_bib
[Hooker, 2007]{Hooker2007 Hooker, G. (2007). \newblock Generalized functional anova diagnostics for high-dimensional functions of dependent variables. \newblock {\em Journal of Computational and Graphical Statistics, 16(3).
1406.1845.extracted_bib
[Johnson and Lindenstrauss, 1984]{Johnson1984 Johnson, W.~B. and Lindenstrauss, J. (1984). \newblock Extensions of lipschitz mappings into a hilbert space. \newblock {\em Contemporary mathematics, 26(189-206):1.
1406.1845.extracted_bib
[Linton and Nielsen, 1995]{Linton1995 Linton, O. and Nielsen, J.~P. (1995). \newblock A kernel method of estimating structured nonparametric regression based on marginal integration. \newblock {\em Biometrika, pages 93--100.
1406.1845.extracted_bib
[Lou et~al., 2013]{Lou2013 Lou, Y., Caruana, R., Gehrke, J., and Hooker, G. (2013). \newblock Accurate intelligible models with pairwise interactions. \newblock In {\em Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 623--631. ACM.
1406.1845.extracted_bib
[Mammen et~al., 1999]{Mammen1999 Mammen, E., Linton, O., Nielsen, J., et~al. (1999). \newblock The existence and asymptotic properties of a backfitting projection algorithm under weak conditions. \newblock {\em The Annals of Statistics, 27(5):1443--1490.
1406.1845.extracted_bib
[Mentch and Hooker, 2014]{MentchTrees Mentch, L. and Hooker, G. (2014). \newblock {Ensemble {Trees and {CLTs: {Statistical {Inference for {Supervised {Learning. \newblock arXiv:1404.6473 [stat.ML].
1406.1845.extracted_bib
[Opsomer and Ruppert, 1998]{Opsomer1998 Opsomer, J.~D. and Ruppert, D. (1998). \newblock A fully automated bandwidth selection method for fitting additive models. \newblock {\em Journal of the American Statistical Association, 93(442):605--619.
1406.1845.extracted_bib
[Opsomer and Ruppert, 1999]{Opsomer1999 Opsomer, J.~D. and Ruppert, D. (1999). \newblock A root-n consistent backfitting estimator for semiparametric additive modeling. \newblock {\em Journal of Computational and Graphical Statistics, 8(4):715--732.
1406.1845.extracted_bib
[Sobol, 2001]{Sobol2001 Sobol, I.~M. (2001). \newblock Global sensitivity indices for nonlinear mathematical models and their monte carlo estimates. \newblock {\em Mathematics and computers in simulation, 55(1-3):271--280.
1406.1845.extracted_bib
[Srivastava et~al., 2014]{pingli2014 Srivastava, R., Li, P., and Ruppert, D. (2014). \newblock {RAPTT: {An {Exact {Two-{Sample {Test in {High {Dimensions {Using {Random {Projections. \newblock arXiv:1405.1792 [stat.ME].
1406.1845.extracted_bib
[Sullivan et~al., 2009]{ebird Sullivan, B.~L., Wood, C.~L., Iliff, M.~J., Bonney, R.~E., Fink, D., and Kelling, S. (2009). \newblock ebird: A citizen-based bird observation network in the biological sciences. \newblock {\em Biological Conservation, 142(10):2282--2292.
1406.1845.extracted_bib
[Wager, 2014]{Wager2014 Wager, S. (2014). \newblock {Asymptotic {Theory for {Random {Forests. \newblock arXiv:1405.0352 [math.ST].
1406.1845.extracted_bib
[Wager et~al., 2014]{WagerIJ Wager, S., Hastie, T., and Efron, B. (2014). \newblock Confidence intervals for random forests: The jackknife and the infinitesimal jackknife. \newblock {\em Journal of Machine Learning Research, 15:1625--1651.