The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
80,941 characters
Conditional nonparametric variable screening by neural factor regression
\def\spacingset#1{\renewcommand{\baselinestretch}
{#1}\small\normalsize} \spacingset{1.0}
\def\r#1{\textcolor{red}{\bf #1}}
\def\b#1{\textcolor{blue}{\bf #1}}
\if00
{
\title{\bf
Conditional nonparametric variable screening by neural factor regression}
\author{Jianqing Fan\\
Operations Research and Financial Engineering, Princeton University\\
Weining Wang \\
Department of Economics, Econometrics and Finance, University of Groningen\\
Yue Zhao\\
Department of Mathematics, University of York}
\maketitle
} \fi
\if10
{
\bigskip
\begin{center}
{\LARGE\bf Conditional nonparametric variable screening via neural network factor regression}
\end{center}
\medskip
} \fi
\bigskip
\bigskip
\begin{abstract}
High-dimensional covariates often admit linear factor structure. To effectively screen correlated covariates in high-dimension, we propose a conditional variable screening test based on non-parametric regression using neural networks due to their representation power. We ask the question whether individual covariates have additional contributions given the latent factors or more generally a set of variables.
Our test statistics are based on the estimated partial derivative of the regression function of the candidate variable for screening and a observable proxy for the latent factors. Hence, our test reveals how much predictors contribute additionally to the non-parametric regression after accounting for the latent factors. Our derivative estimator is the convolution of a deep neural network regression estimator and a smoothing kernel. We demonstrate that when the neural network size diverges with the sample size, unlike estimating the regression function itself, it is necessary to smooth the partial derivative of the neural network estimator to recover the desired convergence rate for the derivative. Moreover, our screening test achieves asymptotic normality under the null after finely centering our test statistics that makes the biases negligible, as well as consistency for local alternatives under mild conditions. We demonstrate the performance of our test in a simulation study and two real world applications.
\end{abstract}
\noindent
{\it Keywords:} Neural networks, factor model, non-parametric regression, non-parametric tests, functional of derivatives, high-dimensionality.
\vfill
\newpage
\spacingset{1.2}
\section{Introduction}
\label{sec:intro}
\subsection{Background}
\label{sec:background}
Variable screening is a powerful tool to expeditiously identify the set of predictors that potentially affect the regression outcome \citep{FanLv2008}. It can reduce a very large number of predictors to a smaller and more manageable set. Then, on this reduced set, one can apply some more refined but computationally demanding variable selection methods such as Lasso, SCAD, Danzig selector \citep{Tibshirani1996,FanLi2001,CandesTao2007,FanLiZhangZou2020} and their non-parametric counterpart FAST-NN \citep{FanGu2023factor}. Conditional marginal screening in parametric regression \citep{barut2016conditional} further augments the screening by conditioning on a \textit{known} set of useful predictors to reduce the impact of the correlations among the predictors, thus making the important {predictors more visible and reducing the false positive and false negative rates in the vanilla screening method}. Despite the aforementioned advances, conducting variable screening or selection in non-parametric regression with high-dimensional inputs remains a challenge due to curse of dimensionality, and
deep learning offers a promising solution thanks to its ability to adapt to unknown low-dimensional structures in multivariate non-parametric regression.
Deep learning has achieved tremendous empirical successes in numerous applications \citep{LeCunBengioHinton2015deep, GoodfellowBengioCourville2016deep}, for instance, in high-dimensional problems such as image recognition \citep{SimonyanZisserman2015, KrizhevskySutskeverHinton2017imagenet}, deep reinforcement learning \citep{mnih2015human}, and large language models \citep{kasneci2023chatgpt,thirunavukarasu2023large}. There is now also a growing literature justifying theoretically the benefit of depth in deep neural networks \citep{Telgarsky2016benefits, Yarotsky2017, ElbrachterPerekrestenkoGrohsBolcskei2021} and their power in alleviating the curse of dimensionality in non-parametric regression via algorithmic learning of unknown low-dimensional structures within complex functions \citep{schmidt2020nonparametric, KohlerLanger2021}. Now being a component of the standard toolbox for statisticians, deep neural networks may nevertheless not be efficient if the dimension of the predictors is very high due to the fundamental limit of multivariate non-parametric regression. Recently \cite{FanGu2023factor} proposed a factor-augmented sparse throughput regression model that simultaneously leveraged the aforementioned adaptivity of deep neural networks and a factor model on the predictors to facilitate variable selection. To complement the variable selection effort in \cite{FanGu2023factor}, \cite{DinhHo2020}, and \cite{HoRichardsonTran2023}, in this paper, we investigate the issue of conditional variable screening with deep neural networks when facing potentially very high-dimensional inputs.
We will assess conditional contribution of a candidate variable for screening by examining its partial derivative in the multivariate non-parametric regression function with a given set of variables and construct our screening test statistics on the moment generating function (MGF) of the smoothed partial derivative of the regression function estimator. We chose a deep neural network as our regression estimator due to its aforementioned algorithmic adaptation to the low-dimensional structure.
Thus, quantifying and improving the performance of derivatives of deep neural networks are integral to our study, and also form an interesting topic of its own right given the importance of derivative estimation in non-parametric regression across diverse research domains and practical applications \citep{GijbelsGoderniaux2005, RondonottiMarronPark2007, horel2020significance}. Traditional non-parametric derivative estimation in general follows one of the following three methods: empirical derivative-type estimation \citep{MullerStadtmullerSchmitt1987, DeBrabanterDeBrabanterGijbelsDeMoor2013, LiuDeBrabanter2020}, kernel/local polynomial-type estimation \citep{GasserMuller1984,fan1996local}, and series/spline-type estimation \citep{Stone1985,ZhouWolfe2000}. However, the first method heavily relies on the existence of an order among the predictor samples and hence naturally applies when the predictor dimension is just one, and the last two methods, just like their original regression estimation counterparts, suffer from the curse of dimensionality problem {when facing predictors of a moderate dimension}.
Despite their deteriorated performances posed by high dimensionality, the traditional derivative estimation methods are relatively amenable to theoretical analysis due to their closed-form solutions. Excluding the empirical derivative method, which avoids explicitly fitting a regression function, the closed-form solutions for the last two methods can be attributed to the {close} connection between estimating the original regression functions and their accompanying derivatives. For instance, in the kernel and the spline methods, a {closed-form} derivative estimator can simply be obtained as the derivative of the original regression function estimator (see, for instance, the discussion between Eqs.~(2) and (3) in \cite{ZhouWolfe2000}). Thus, unsurprisingly, in these methods the quality of the derivative estimation closely follows the quality of the original regression function estimation. However, for deep neural networks, estimating the \textit{derivative} of a regression function can be quite different from the task on the regression function itself due to its smoothness.
Take for example any candidate regression function estimator $m$ within the canonical neural network class, precisely defined in \eqref{eq:base_NN_class} later, built from the popular ReLU activation function. The first order partial derivatives of $m$ are necessarily piecewise constant due to the piecewise linear nature of the ReLU function, and the second order partial derivatives of $m$ are necessarily \textit{zero} almost everywhere, \textit{irrespective} of the underlying truth that the function $m$ may attempt to recover. Numerically, this is easily demonstrable through standard software packages such as PyTorch. In addition, the first order partial derivatives of $m$ could exhibit convergence behaviors qualitatively different from the convergence behavior of $m$ itself as we will explain in Section~\ref{sec:test_stat}. Such phenomena are in sharp contrast to the traditional non-parametric estimators, and can in part be attributed to the features of deep neural networks: they are highly non-linear and lack easily interpretable closed-form expressions. {Moreover, existing asymptotic results on functionals acting on deep neural networks almost invariably assume some continuity of the functionals with respect to their inputs. However, when such functionals in effect act on the derivatives of neural networks, the assumed continuities {could} break down due to the aforementioned different convergence behaviors of the said derivatives.} We will explain and address this discrepancy as we progress through the paper.
\subsection{Our method and contribution}
\label{sec:our_method}
Our main contribution in this paper is a non-parametric conditional variable screening test using deep neural networks when the ambient dimension of the inputs is potentially very high. In addition to complementing the variable selection methods, our contribution is also an advance over the aforementioned paper by \cite{barut2016conditional} in that we work with non-parametric regression and our conditioning variables are not known in advance but instead are extracted from the inputs with the help of the factor model as in \cite{fan2022learning} and \cite{FanGu2023factor}. We note that \cite{horel2020significance} has also conducted screening test based on partial derivatives (in low dimensions). However their study focused on single-layer neural networks, and hence does not benefit from the power nor reveal the intricacies of deep neural networks in derivative estimation.
In addition to the general procedures outlined above, we also make the following contributions which could be of independent interest:
\begin{enumerate}[wide, labelwidth=!, labelindent=0pt, label=(\alph*)]
\item[1)]
We rigorously derive the size and power of our proposed {test statistics employing non-linear functionals of} truly deep neural networks whose architecture can become more complex as the sample size increases. For computational ease, we propose a simple variance estimator whose properties we properly characterize. Our resulting test statistics are straightforward to compute and sufficiently precise to accommodate local alternatives. {Moreover, we address the aforementioned continuity issue of functionals acting on the derivatives of neural networks by proposing improved estimation of the said derivatives which naturally leads to our next contribution.}
\item[2)]
We exploit the regression function algorithmically learned by deep neural networks in order to recover high-quality derivative estimation. To achieve this goal, we employ a smoothing technique on deep neural network estimators to regularize their derivatives. Moreover, out method is applicable to
regularize the derivatives of alternative machine learning techniques, thus paving the way for their application in derivative estimation. Last but not least, although we focus on (conditional) marginal screening in the present paper, our method can be generalized easily to test higher-order effect and variable interaction using higher-order and mixed derivatives respectively.
\end{enumerate}
Through simulation studies, we illustrate the favorable size and power performance of our test statistics, and the benefit of the smoothing operation in generating accurate derivative estimators. We also highlight the potential of our test statistics as a viable tool for model specification in nonlinear factor models via two empirical applications.
\subsection{Notations, conventions, and manuscript organization}
\label{sec:notation}
Let $\mathbb{N}$ denote the set of positive integers. For a vector $\mathbf{v} = (v_1, \ldots, v_d)^\top \in \mathbb{R}^d$, we let $\|\mathbf{v}\|_p = (\sum_{k=1}^d |v_i|^p)^{1/p}$ be the $\ell_p$ norm of $\mathbf{v}$. For a matrix $\mathbf{A} = (a_{j,k})_{1\le j\le l, 1\le k \le m}\in\mathbb{R}^{l\times m}$, we define the operator norm $\|\mathbf{A}\|_{\textup{op}} = \max_{\mathbf{v}\in\mathbb{R}^m: \|\mathbf{v}\|_2 = 1} \|\mathbf{A}\mathbf{v}\|_2$. We let $C$ denote an (absolute) constant that may change for each occurrence, and let ``$\lesssim$'' denote an inequality that holds up to such a multiplicative factor $C$; moreover, let $c$, $C$ and $M$ with super/subscripts denote constants with particular (though often non-specified) values.
Limits are taken as $n\to\infty$ unless otherwise stated. For positive number sequences $(a_n: n\ge 1)$ and $(b_n: n\ge 1)$, we denote $a_n\lesssim b_n$ if there exists a positive constant $C$ such that $a_n/b_n\le C$ (for all $n$), and denote $a_n={\scriptstyle{\mathcal{O}}}(b_n)$ (resp. $a_n\sim b_n$) if $a_n/b_n\to 0$ (resp. $a_n\lesssim b_n$ and $b_n\lesssim a_n$). We use ``$\rightarrow_{d}$'' to denote convergence in distribution. Finally, let $\|\cdot\|_{L_2}=\|\cdot\|_{L_2(\mathbb{P})}$ denote the $L_2$ norm of the argument function with respect to the measure $\mathbb{P}$ to be formally introduced in Section~\ref{sec:reg_model}, so $\|f\|_{L_2}= \{ \int f^2 {\rm d}\mathbb{P} \}^{1/2}$. Sections in the supplement are labelled alphabetically.
We organize the remainder of our paper as follows. Section~\ref{sec:testhypothesis} specifies our regression model and screening test, and develop our derivative estimators and test statistics. Section~\ref{sec:regression_function} depicts the accompanying theoretical properties of the derivative estimators and test statistics from Section~\ref{sec:testhypothesis}. Section~\ref{sec:simulations} presents a simulation study. Section~\ref{sec:emp_app} applies our test to an empirical example in asset pricing and another in macroeconomics time series. Section~\ref{sec:conclusion} concludes and suggests several extensions. Additional results for the empirical application, proofs and supporting details are deferred to the supplementary materials.
\section{Derivative estimator and test statistics}
\label{sec:testhypothesis}
\subsection{A conditional screening test through latent factors}
\label{sec:reg_model}
Our screening test aims to tackle the high-dimensional regime where the ambient dimension $d$ of our observed predictors ${\bm{X}}\in\mathbb{R}^d$ can increase with the sample size $n$. Even with the remarkable capacity of deep neural networks to represent complex functions, our task is still infeasible if $d$ grows too fast.
To address this issue, we assume that ${\bm{X}}$ admits the following factor model \citep{FanGu2023factor}:
\begin{align}
\label{eq:factor_observation}
{\bm{X}}= \mathbf{B} {\bm{F}}+{\bm{u}},\quad \mathbb{E}({\bm{u}}|{\bm{F}})=\mathbf{0} ,
\end{align}
where ${\bm{F}} \in \mathbb{R}^r$ is a vector of latent factors, $\mathbf{B}\in\mathbb{R}^{d\times r}$ (usually $d\gg r$ for the high-dimensional ${\bm{X}}$) is the factor loading matrix, and ${\bm{u}}$ is the vector of the idiosyncratic noises. We let $(Y,{\bm{F}},{\bm{X}})$ have joint distribution $\mathbb{P}$ which also determines the distribution of ${\bm{u}}$ by \eqref{eq:factor_observation}.
The conditional marginal screening is to see whether a component $X_j$ of ${\bm{X}}$ has additional contributions to the response variable $Y$ given ${\bm{F}}$. Then, the problem involves the working regression function $m_0(\cdot):\mathbb{R}^{r+1}\rightarrow\mathbb{R}$ defined through
\begin{equation}
\begin{gathered}
Y = m_0({\bm{F}},X_j) + \epsilon , \quad \text{where}~m_0({\bm{F}},X_j) \equiv \mathbb{E}[Y|{\bm{F}},X_j] .
\end{gathered}
\label{eq:reg_model}
\end{equation}
The additional contribution of variable $X_j$ is measured through the partial derivative
\begin{align*}
\textstyle m_{0,j} = m_{0,j}({\bm{f}},x_j) \equiv \frac{\partial}{\partial x_j} m_0({\bm{f}},x_j) .
\end{align*}
We impose a blanket notational convention that a function with subscript $j$ denotes the partial derivative with respect to the last argument, which will almost always be $x_j$, while keeping the other argument (in this case, ${\bm{f}}$) fixed. Then, to examine whether $X_j$ has additional contribution, we propose the null and alternative hypotheses
\begin{gather}\label{hypo_null}
H_0: \text{for all}~{\bm{f}}, x_j\in\mathbb{R}^{r+1},~m_{0,j}({\bm{f}},x_j) = 0,~\text{against} \\
\label{hypo_alt}
H_A: \text{for some}~{\bm{f}}, x_j\in\mathbb{R}^{r+1},~\text{we have}~m_{0,j}({\bm{f}},x_j) \neq 0 .
\end{gather}
We will assume that $(Y,{\bm{F}},{\bm{X}})$ is a $\mathbb{R}\times[-b,b]^{r+d}$ valued random vector. The function $m_0$ is akin to $\mathbb{E}[Y|X_j]$ in condition~C in \cite{FanFengSong2011} on \textit{un}conditional non-parametric screening, but here we exercise finer control over the potential role of $X_j$ by introducing ${\bm{f}}$ in $m_0=m_0({\bm{f}},x_j)$. Note that {working with} $m_0$ involves estimating the latent factors ${\bm{F}}$, which will be discussed in Section~\ref{sec:highd}. By Remark~\ref{rmk:specification_test} in the supplement, the screening hypotheses for $X_j$ in \eqref{hypo_null} and \eqref{hypo_alt} are also equivalent to a significance test for the $j$-th idiosyncratic term $u_j$. We deliberately avoid working with the full conditional expectation of $Y$, namely $m_0^*({\bm{F}},{\bm{X}}) \equiv \mathbb{E}[Y|{\bm{F}},{\bm{X}}]$, because this is likely infeasible if the dimension of ${\bm{X}}$ is too high {without} additional structure. In contrast, the total number of variables $r+1$ in our working regression function $m_0$ and its derivative $m_{0,j}$ is much less than $d$, which allows us to circumvent the problem of testing a potentially very large number of coordinates of ${\bm{X}}$ simultaneously \citep{FanGu2023factor}. Nevertheless, our model retains its validity even when the regression outcome $Y$ involves more coordinates of ${\bm{X}}$, as we will explain in Remark~\ref{rmk:misspecification}.
Although we shall focus on the high-dimensional ${\bm{X}}$ case with growing $d$, our screening test also easily accommodates the low-dimensional ${\bm{X}}$ case, as we will comment in Section~\ref{sec:lowd}. At the other end of the spectrum, our findings can be generalized to testing ${\bm{X}}$ over a fixed or expanding set of coordinates that form a subset of $\{1,\cdots,d\}$, as we will briefly outline in Section~\ref{sec:conclusion}.
\subsection{Initial regression estimator through diversified projection}
\label{sec:highd}
We focus on multi-layer feed-forward neural networks that are fully connected between adjacent layers \cite[p.~75]{AnthonyBartlett1999}. In this paper we exclusively consider the Rectified Linear Unit (ReLU) activation function defined as $\sigma=\sigma_{\text{ReLU}}(x)=\max \{x, 0\}$ due to its widespread popularity \citep[Section~6.1]{GoodfellowBengioCourville2016deep}.
For now, we consider neural network functions with a generic input dimension ${\overline{d}}$. We finalize the structure of our neural networks by specifying the tuple $(L, \mathbf{k})$ where $L\in\mathbb{N}$ represents the number of hidden layers and the width vector $\mathbf{k}=({\overline{d}}, k_1, \ldots, k_L, k_{L+1}=1) \in \mathbb{N}^{L+2}$ specifies the number of nodes (i.e., neurons) in each hidden layer and the input/output dimensions. More precisely, such a neural network function $f:\mathbb{R}^{\overline{d}}\rightarrow\mathbb{R}$ with architecture $(L, \mathbf{k})$ is given by
\begin{align}
\label{eq:NN_function_output}
\textstyle f(\bm{x})= \mathcal{L}_{L+1} \circ {\overline{\sigma}}_L \circ \mathcal{L}_L \circ {\overline{\sigma}}_{L-1} \circ \dots \circ \mathcal{L}_2 \circ {\overline{\sigma}}_1 \circ \mathcal{L}_1(\bm{x}) ;
\end{align}
here $L_\ell(\bm{z})=\mathbf{w}_\ell \bm{z} +\bm{b}_\ell$ is an affine map with weight matrix {$\mathbf{w}_{\ell}\in\mathbb{R}^{k_\ell \times k_{\ell-1}}$ if $\ell\geq 2$, $\mathbf{w}_{1}\in\mathbb{R}^{k_1 \times \overline{d}}$} and bias vector $\bm{b}_\ell\in\mathbb{R}^{k_\ell}$, and the function ${\overline{\sigma}}_\ell:\mathbb{R}^{k_\ell}\rightarrow\mathbb{R}^{k_\ell}$ applies the ReLU activation function $\sigma$ entry-wise. Additionally, we will truncate the (univarite) output of $f$ at a pre-specified constant level $M>0$ with the truncation operator $T_M(x) \equiv \textup{sgn}(x) \min\{|x|, M\}$. We denote such a collection of truncated $f$ by $\mathcal{F}_n({\overline{d}})$ where for brevity of notation we only retain the dependence on the input dimension ${\overline{d}}$:
\begin{align}
\mathcal{F}_n({\overline{d}}) = \big\{ T_M(f): f\text{ is of the form \eqref{eq:NN_function_output}}~
\text{with $L$ hidden layers and width vector $\mathbf{k}$} \big\},
\label{eq:base_NN_class}
\end{align}
where both $L$ and $\mathbf{k}$ can scale with the sample size $n$.
Next, we address the latency issue of the unobserved factor ${\bm{F}}$. For a given diversified weight $\mathbf{w} \in\mathbb{R}^d$, \cite{fan2022learning} note from \eqref{eq:factor_observation} that within the projection $\mathbf{w}^T {\bm{X}}$, the projected idiosyncratic terms $\mathbf{w}^T {\bm{u}}$ will be negligible in high-dimension due to the law of the averages, so $\mathbf{w}^T {\bm{X}}$ yields an approximate linear combination of ${\bm{F}}$. To yield $r$ factors, we need at least $r$ projections. Since $r$ is {unknown}, one specifies an upper bound ${\overline{r}}$ {of} $r$.
Let $\mathbf{W}\in\mathbb{R}^{d\times {\overline{r}}}$ be a pre-trained \textit{diversified projection matrix} as termed by Definition~3 in \cite{FanGu2023factor}, and let $\mathbf{H}\equiv d^{-1}\mathbf{W}^\top \mathbf{B}\in\mathbb{R}^{{\overline{r}}\times r}$.
Then, by \eqref{eq:factor_observation}, we have
$$
\widetilde{\bm{F}} \equiv d^{-1}\mathbf{W}^\top {\bm{X}} = \mathbf{H} {\bm{F}} + d^{-1}\mathbf{W}^\top {\bm{u}} \approx \mathbf{H} {\bm{F}} \in\mathbb{R}^{{\overline{r}}},
$$
by the law of the averages. Hence, $\mathbf{H}^\dagger \widetilde{\bm{F}} \approx {\bm{F}}$ under some appropriate conditions, where $\mathbf{H}^\dagger$ is the pseudo-inverse of $\mathbf{H}$. We call $\widetilde{\bm{F}}$ as the \textit{diversified factor}, which is observable and a proxy of the latent factor ${\bm{F}}$.
In practice, acquiring $\mathbf{W}$ in advance is necessary, either through domain expertise or data-driven methods. For instance, Proposition~1 in \cite{FanGu2023factor} proposes estimating $\mathbf{W}$ via pretraining, where a tiny portion of size approximately $\log(n)$ of the samples is reserved to {extract the top ${\overline{r}}$ principal components in order to} construct a $\mathbf{W}$ that will satisfy Assumption~\ref{ass:W} with high probability. Henceforth we shall assume that $\mathbf{W}$ is pre-determined, exogenous and fixed.
With the proxies of the latent factors, for each given $j$, the coordinate of the candidate variable for screening, we can compute $\{Y_i, \widetilde{\bm{F}}_i, X_{i,j}\}_{i=1}^n$ based on the observed data, where $\widetilde{\bm{F}}_i= d^{-1}\mathbf{W}^\top {\bm{X}}_i$. Then we fit the neural network regression
\begin{align}
\label{eq:def_wh_m_n_highd}
\textstyle \widehat g_n = \operatorname*{arg\!\min}_{g\in\mathcal{F}_n({\overline{r}}+1)} \frac{1}{n}\sum_{i=1}^n \{Y_i-g(\widetilde{\bm{F}}_i,X_{i,j})\}^2 = \operatorname*{arg\!\min}_{g\in\mathcal{F}_n({\overline{r}}+1)} \mathbb{P}_n \ell(\cdot;g)
\end{align}
where $\ell(y,\widetilde{\bm{f}},x_j;g)=\{y-g(\widetilde{\bm{f}},x_j)\}^2$ is the square loss and $\mathbb{P}_n$ is the empirical distribution.
Intuitively, $\widehat g_n$ is a neural network estimation of
\begin{align}
\label{eq:def_g0}
g_0 = g_0(\widetilde{\bm{f}},x_j) \equiv m_0(\mathbf{H}^\dagger {\bm{f}}, x_j).
\end{align}
Indeed, Lemma~\ref{lemma:factor} shows that $g_0$ approximates $m_0$ well.
\subsection{Initial test statistic and its improvements}
\label{sec:test_stat}
To test our screening hypotheses \eqref{hypo_null} and \eqref{hypo_alt}, we start from $\widehat g_{n,j}$, which by the convention in Section~\ref{sec:reg_model} is the partial derivative of $\widehat g_n$ in \eqref{eq:def_wh_m_n_highd} with respect to $x_j$: $\widehat g_{n,j}(\widetilde{\bm{f}},x_j)=\partial\widehat g_n(\widetilde{\bm{f}},x_j) / \partial x_j$. Now, $\widehat g_{n,j}$ can be regarded as an estimator of the derivative $m_{0,j}$. We then consider the following \textit{initial} MGF/exponentially tilted test statistic:
\begin{align}
\label{eq:eta_t_naive}
\eta_t(\widehat g_n)=\mathbb{P}_n\{ \exp(t\,\widehat g_{n,j}) - 1 \}.
\end{align}
Partly owing to the fact that a MGF uniquely characterizes the distribution of the underlying random variable (specifically when considering the MGF over an interval including zero), MGF-based tests are popular in the literature \citep{EppsSingletonPulley1982,BaringhausEbnerHenze2017}. For instance, they have been extensively employed in testing (multivariate) normality \citep{EbnerHenze2020}. Under the null hypothesis $H_0: m_{0,j}=0$, we expect $\widehat g_{n,j}$ to be close to zero and so $\eta_t(\widehat g_n)$ is also centered around zero.
To enhance the quality of both the derivative estimator $\widehat g_n$ and the test statistic $\eta_t(\widehat g_n)$, we will further conduct the following refinements sequentially:
\begin{enumerate}[wide, labelwidth=!, labelindent=0pt, label=(\alph*), topsep=0pt,itemsep=-1ex,partopsep=1ex,parsep=1ex,leftmargin =0.2 in]
\item {\bf Smoothing}.
\label{improve:smooth}
As already alluded to in the introduction, in derivative estimation, the straightforward plugin estimator $\widehat g_{n,j}$ derived from $\widehat g_n$ in \eqref{eq:def_wh_m_n_highd} may perform poorly due to irregularities of the neural network functions. For example, for the $L$-time iterated sawtooth function $\zeta_L$ (see Lemma~2.4 in \cite{Telgarsky2015} and Figure~\ref{fig:sawtooth}) which is implementable by a neural network of no more than $L$ hidden layers and three nodes per layer, we have {$\|\dot\zeta_L\|_{L_2}\ge C 2^L \|\zeta_L\|_{L_2}$}, which implies that the quality of derivative estimation can become significantly worse than that of regression function estimation as the depth $L$ increases, a regime precisely of interest for deep neural networks. This phenomenon results from the increasingly oscillatory behaviour of the neural network functions as the depth increases, and hence is not tied specifically to the ReLU activation function. To address this issue, we will refine the initial estimator $\widehat g_{n,j}$ to obtain a smoothed derivative estimator $\widehat g_{n,j}^\textup{s}$, which will allow us to recover a faster convergence rate; see Theorem~\ref{thm:m_hat_est_master} and Theorem~\ref{thm:m_check_est_master}; this will in turn address the continuity issue {arising from second-order terms in our MGF test statistics} acting on the derivatives of neural networks. Throughout this paper, we adhere to the convention that the superscript ``$\textup{s}$'' in upright font denotes the smoothed variant of the preceding function. We will provide more detailed calculation for our observation on the sawtooth function and more comprehensive rationale behind the smoothing operation in Section~\ref{sec:sawtooth}.
\item{\bf Centering}.
\label{improve:tweaking}
We must further refine the initial $\widehat g_{n,j}^\textup{s}$ {to arrive at $\widecheck g_{n,j}^\textup{s}$ in \eqref{eq:f_check} in order} to center our test statistics. This debiasing step is analogous to the concept of targeted machine learning in the literature \citep{van2011targeted, VanderLaanRose2018}, and also relates to other research on debiased machine learning \citep{quintas2021riesznet, kennedy2022semiparametric}. However, we will justify this refinement step independently without directly referencing the targeted/debiased machine learning literature.
\item{\bf Truncation}.
\label{improve:truncate}
To treat the technical possibility of the unboundedness of $\widecheck g_{n,j}^\textup{s}$, we will truncate $\widecheck g_{n,j}^\textup{s}$ appropriately in our final test statistics.
\item{\bf Uniformity over $t$}.
\label{improve:uniform}
We will extend the fixed-$t$ statistic to the statistics aggregated over a range of $t$ to enhance the power of the test \citep{fan2001generalized, fan2015power}.
\end{enumerate}
We focus solely on the smoothing step~\ref{improve:smooth} in this section, and will introduce step~\ref{improve:tweaking} in Section~\ref{sec:tweaking}, and both \ref{improve:truncate} and \ref{improve:uniform} in Section~\ref{sec:uniform}.
We first describe the smoothing operation applied to a generic function $g:\mathbb{R}^{\bar{r}+1}\rightarrow\mathbb{R}$ (or analogously $m:\mathbb{R}^{r+1}\rightarrow\mathbb{R}$). Let $K$ be a univariate continuously differentiable kernel supported on $[-1,1]$ with derivative $\dot K$.
Then, for the generic function $g$, denote its smoothed version in the variable $x_j$ by $g^\textup{s}=g^\textup{s}(\widetilde{\bm{f}}, x_j)
= \int g(\widetilde {\bm{f}}, z) K_h(x_j - z)dz$, where $K_h(\cdot)=K(\cdot/h) / h$ with a bandwidth parameter $h$. Now, the smoothed function $g^\textup{s}$ becomes differentiable in $x_j$ everywhere with the partial derivative given by
\begin{align}
\label{eq:smooth_master}
\textstyle g_j^\textup{s}(\widetilde{\bm{f}}, x_j) = \frac{\partial}{\partial x_j} g^\textup{s}(\widetilde{\bm{f}}, x_j) = \int_{x_j-h}^{x_j+h} g(\widetilde{\bm{f}},z) \dot K_h(x_j-z){\rm d} z = \frac{1}{h} \int_{-1}^1 g(\widetilde{\bm{f}},x_j-ah) \dot K(a) {\rm d} a.
\end{align}
Accordingly, we let our initial smoothed derivative estimator based on $\widehat g_n$ be $\widehat g_{n,j}^\textup{s}$. In practice, we can select the bandwidth $h$ through cross-validation; see our Remark~\ref{rmk:bandwidth}.
\subsection{Estimating the score function}
\label{sec:Riesz_est}
To test the screening hypotheses \eqref{hypo_null} and \eqref{hypo_alt}, we proceed to estimate the score function corresponding to the statistic $\eta_t^\textup{s}(\widehat g_n)\equiv\mathbb{P}_n\{ \exp(t\,\widehat g_{n,j}^\textup{s}) - 1 \}$ now refined over \eqref{eq:eta_t_naive}. This score function $\alpha_{t,n}^*:\mathbb{R}^{{\overline{r}}+1}\rightarrow\mathbb{R}$ satisfies, under the null and for any $L_2(\mathbb{P})$-integrable $\alpha$, the equality $\int \alpha\alpha_{t,n}^* {\rm d}\mathbb{P} = \int_{\Omega_h} \alpha_j^\textup{s}(\widetilde{\bm{f}},x_j) {\rm d}\mathbb{P}$, where $\alpha_j^\textup{s}$ is the smoothed derivative of $\alpha$ and $\Omega_h=\{\omega\in\Omega: X_j(\omega)\in\mathcal{B}_h\}$ {is the interior sample space corresponding to} the interior set $\mathcal{B}_h=[-b+h,b-h]$. (Such an $\alpha_{t,n}^*$ exists because it is the Riesz representer of the directional derivative functional of our test statistic, as we will formally explain in Section~\ref{sec:variance_prep}.) Therefore, the mean-square error (MSE) can be expressed as
\begin{align}
&\mathbb{E}(\alpha-\alpha_{t,n}^*)^2(\widetilde{\bm{F}},X_j) = \mathbb{E} \alpha^2(\widetilde{\bm{F}},X_j) - 2 \mathbb{E} (\alpha \alpha_{t,n}^*)(\widetilde{\bm{F}},X_j) + \mathbb{E} \alpha_{t,n}^{*2}(\widetilde{\bm{F}},X_j) \nonumber \\
& \textstyle {\stackrel{\text{under}~H_{0}}{=}} \textstyle \mathbb{E} \alpha^2(\widetilde{\bm{F}},X_j) - 2 \int_{\Omega_h} \alpha_j^\textup{s}(\widetilde{\bm{f}},x_j) {\rm d}\mathbb{P} + \mathbb{E} \alpha_{t,n}^{*2}(\widetilde{\bm{F}},X_j).
\label{eq:Rnull}
\end{align}
Discarding the last term $\mathbb{E}[\alpha_{t,n}^{*2}(\widetilde{\mathbf{F}},X_j)]$, which does not depend on $\alpha$, the MSE can be estimated by the empirical loss
\begin{align}
\label{eq:Rhat}
\textstyle \widehat{R}_n^{\textup{null}}(\alpha) = \frac{1}{n} \sum_{i=1}^n \alpha^2(\widetilde{\bm{F}}_i, X_{i,j}) - 2 \frac{1}{n} \sum_{i\in\mathcal{I}_h} \alpha_j^\textup{s}(\widetilde{\bm{F}}_i, X_{i,j}),
\end{align}
where $\mathcal{I}_h=\{i\in\{1,\dots,n\}: X_{i,j}\in\mathcal{B}_h\}$ is the interior index set.
Hence, we estimate $\alpha_{t,n}^*$ by
\begin{align}
\label{eq:alpha_hat}
\textstyle \widehat\alpha_n = \operatorname*{arg\!\min}_{\alpha\in\mathcal{F}_n({\overline{r}}+1)} \widehat{R}_n^{\textup{null}}(\alpha),
\end{align}
which is a neural network estimator of the score function.
We tailor the loss function $\widehat{R}_n^{\textup{null}}$ and the estimator $\widehat\alpha_n$ to the null hypothesis \eqref{hypo_null}. One notable computational advantage is that $\widehat{R}_n^{\textup{null}}$ and hence $\widehat\alpha_n$ are independent of both the response and the value of $t$: A single optimization suffices to yield $\widehat\alpha_n$ for all $t$. However, the potential trade-off is a bias of $\widehat\alpha_n$ with respect to $\alpha_{t,n}^*$ induced under the alternative hypothesis, which will be characterized in Theorem~\ref{thm:Riesz_est}. Numerically, computing the right-hand side in \eqref{eq:Rhat} is straightforward, as detailed in Lemma~\ref{lemma:alpha_j_numerical} in the supplement.
\subsection{Centering test statistics}
\label{sec:tweaking}
Our refined statistic $\eta_t^\textup{s}(\widehat g_n)$ needs to be further centered to ensure asymptotic normality. To achieve this, ideally we adjust $\widehat g_n$ in the direction of the score function $\alpha_{t,n}^*$ to attain a smaller loss than in \eqref{eq:def_wh_m_n_highd}. With the estimate $\widehat\alpha_n$ of $\alpha_{t,n}^*$ given by \eqref{eq:alpha_hat}, one naturally uses
\begin{align}
\label{eq:f_check}
\widecheck g_n = \widehat g_n + \widehat{\delta}_\textup{t} \widehat \alpha_n
\end{align}
to arrive at a debiased statistic $\eta_t^\textup{s}(\widecheck g_n)$, where
\begin{align}
\textstyle \widehat{\delta}_\textup{t} & \textstyle = \frac{1}{\sum_{i=1}^n \widehat\alpha_n^2(\widetilde{\bm{F}}_i,X_{i,j})} \sum_{i=1}^n \left\{ Y_i - \widehat g_n(\widetilde{\bm{F}}_i,X_{i,j}) \right\} \widehat\alpha_n(\widetilde{\bm{F}}_i,X_{i,j})
\label{eq:k_check}
\end{align}
is chosen to minimize $\mathbb{P}_n \ell(\cdot;\widehat g_n + \delta_\textup{t} \widehat\alpha_n)$ in $\delta_\textup{t}$, as derived in Section~\ref{sec:proof_prop:tweaking}. We can interpret this update step as a form of the targeted machine learning \citep{VanderLaanRose2018} specialized to our deep neural network context.
In Section~\ref{sec:tweaking_proof} we will show that $\widecheck g_n$ indeed satisfies an approximate minimization condition: for an infinitesimal stepsize $r_{\textup{inf},n}$,
\begin{align}
\label{eq:min_con}
\mathbb{P}_n\ell(\cdot; \widecheck g_n)-\mathbb{P}_n\ell(\cdot; \widecheck g_n \pmr_{\textup{inf},n}\alpha_{t,n}^*) \le r_{\textup{inf},n} b_n
\end{align}
for a small tolerance $b_n$. This will in turn center our test statistics.
\subsection{Final test statistic and its uniform extensions}
\label{sec:uniform}
To alleviate the technical possibility of the unboundnessed of our derivative estimator, we let $\Psi: \mathbb{R} \rightarrow \mathbb{R}$ be a bounded truncation function that also satisfies Assumption~\ref{ass:truncation}; we provide an example of such a truncation function below the assumption. We then define the final (fixed-$t$) test statistic to be $\widecheck\eta_t^\textup{s}(\widecheck g_n)$ where the operator $\widecheck\eta_t^\textup{s}(\cdot)$ acts on functions $g:\mathbb{R}^{{\overline{r}}+1}\rightarrow\mathbb{R}$ as
\begin{gather}
\textstyle \widecheck\eta_t^\textup{s}(g) = \frac{1}{n}\sum_{i\in\mathcal{I}_h} \left[ \exp\{t\,\Psi(g_j^\textup{s}(\widetilde{\bm{F}}_i, X_{i,j})) \} - 1 \right] .
\label{eta_t_s_wc}
\end{gather}
Its standard error is estimated as $t \|\widehat\epsilon \widehat\alpha_n\|_{L_2(\mathbb{P}_n)}= t [ \sum_{i=1}^n \{ \widehat\epsilon_i \widehat\alpha_n(\widetilde{\bm{F}}_i,X_{i,j}) \}^2 / n ]^{1/2}$, where $\widehat\epsilon_i=Y_i-\widecheck g_n(\widetilde{\bm{F}}_i,X_i)$ is the $i$-th residual.
This leads to a level-$\alpha$ test of the null hypothesis~\eqref{hypo_null} with the critical region
\begin{align}
\label{test:fixed_t}
\left| \frac{ \sqrt{n} }{t} \frac{\widecheck\eta_t^\textup{s}(\widecheck g_n) }{ \|\widehat\epsilon \widehat\alpha_n\|_{L_2(\mathbb{P}_n)} } \right| > z_{1-\alpha/2}
\end{align}
where $z_{1-\alpha/2}$ is the $1-\alpha/2$ quantile of a standard normal distribution; the choice of this critical value is supported by the first part of \eqref{eq:H0_normality_Ln} in Theorem~\ref{thm:main_null_master}.
To increase the power of our test, we aggregate the ensemble of test statistics $\widecheck\eta_t^\textup{s}(\widecheck g_n)$ over a range of $t$. We consider two types of aggregation: a sup-statistic and an integrated weighted square-statistic (or simply a square-statistic). In both cases, let $\delta > 0$ be a small but fixed constant, let $T > \delta$ be a potentially large but also fixed constant, and define $\mathcal{T}_\delta= [-T, -\delta] \cup [\delta, T]$. The two aforementioned aggregations of test statistics and their associated critical regions are as follows.
\begin{enumerate}[wide, labelwidth=!, labelindent=0pt, label=(\alph*), leftmargin =0.2 in]
\item[1)] { \bf Sup-statistic}:
\begin{equation}\label{sup}
\widehat Z = \sqrt{n} \sup_{t\in\mathcal{T}_\delta} \left| \dfrac{ \widecheck\eta_t^\textup{s}(\widecheck g_n ) }{t} \right|, \quad \text{with critical region:}~\dfrac{ \widehat Z }{ \|\widehat\epsilon \widehat\alpha_n\|_{L_2(\mathbb{P}_n)} } > z_{1-\alpha/2} .
\end{equation}
\item[2)] { \bf Square statistic}:
\begin{equation}\label{square}
\widehat\chi^2 = \dfrac{n}{ \int_{\mathcal{T}_\delta} w(t) {\rm d} t } \int_{\mathcal{T}_\delta} \dfrac{1}{t^2} \left\{ \widecheck\eta_t^\textup{s}(\widecheck g_n) \right\}^2 w(t) {\rm d} t, \quad \text{with critical region}~\dfrac{ \widehat\chi^2 }{ \|\widehat\epsilon \widehat\alpha_n\|_{L_2(\mathbb{P}_n)}^2 } > \chi^2_{1,1-\alpha}
\end{equation}
where $w$ is an integrable weight function such as $w(t)=\exp(-\beta t^2)$ for a given $\beta>0$, and
$\chi^2_{1,1-\alpha}$ is the $1-\alpha$ quantile of a chi-square distribution with one degree of freedom.
\end{enumerate}
The critical values are supported by the last two parts of \eqref{eq:H0_normality_Ln} in Theorem~\ref{thm:main_null_master}.
\subsection{Summary of proposed methods}
The major implementation steps of our partial derivative estimation and conditional screening test procedures in Section~\ref{sec:testhypothesis} are summarized in Algorithm~\ref{algorithmhighd}.
{\small \spacingset{1.5}
\begin{algorithm}
\hrulefill
\caption{Partial derivative estimation \& conditional screening test by deep neural network}
\SetKwInput{KwData}{Input}
\KwData{Observed samples $Y_i, {\bm{X}}_i \in \mathbb{R}^d$, $i\in \{1,\cdots, n\}$, and an exogenous diversified projection matrix $\mathbf{W}\in\mathbb{R}^{d\times{\overline{r}}}$; specification of the coordinate $j$ for screening.}
\KwResult{Preliminary and refined regression estimators $\widehat g_n(\cdot)$ and $\widecheck g_n(\cdot)$; their smoothed derivative estimators $\widehat g_{n,j}^\textup{s}(\cdot)$ and $\widecheck g_{n,j}^\textup{s}(\cdot)$; $\widehat\alpha_n(\cdot)$, and test results.}
\begin{enumerate} \itemsep -0.15in
\item
Compute $\widetilde{\bm{F}}_i = d^{-1} \mathbf{W}^\top {\bm{X}}_i$, for all $i \in \{1,\dots,n\}$.
\item
Estimate the regression function $m_0$ in \eqref{eq:reg_model} by the neural network estimator $\widehat g_n$ from \eqref{eq:def_wh_m_n_highd}. {Conduct cross validation (see Remark~\ref{rmk:bandwidth}), or otherwise, to find a smoothing bandwidth $h$.} \\
\item
Define the loss function $\widehat{R}_n^{\textup{null}}(\cdot)$ as
\begin{align*}
\textstyle \widehat{R}_n^{\textup{null}}(\alpha) = \frac{1}{n} \sum_{i=1}^n \alpha^2(\widetilde{\bm{F}}_i,X_{i,j}) - 2 \frac{1}{n} \sum_{i\in\mathcal{I}_h} \alpha_j^\textup{s}(\widetilde{\bm{F}}_i,X_{i,j}) ,
\end{align*}
which can be approximated via Lemma~\ref{lemma:alpha_j_numerical}, and let the estimator $\widehat\alpha_n$ of $\alpha_{t,n}^*$ be
\begin{align*}
\textstyle \widehat\alpha_n = \operatorname*{arg\!\min}_{\alpha\in\mathcal{F}_n({\overline{r}}+1)} \widehat{R}_n^{\textup{null}}(\alpha).
\end{align*}
\item
{Return $ \widecheck g_n= \widehat g_n + \widehat{\delta}_\textup{t} {\widehat\alpha_n}$, where}
\begin{align*}
\textstyle \widehat{\delta}_\textup{t} = \frac{1}{\frac{1}{n} \sum_{i=1}^n {\widehat\alpha_n}^2(\widetilde{\bm{F}}_i,X_{i,j})} \frac{1}{n} \sum_{i=1}^n \left\{ Y_i - \widehat g_n(\widetilde{\bm{F}}_i,X_{i,j}) \right\} {\widehat\alpha_n}(\widetilde{\bm{F}}_i,X_{i,j}) .
\end{align*}
\item
Let $\|\widehat\epsilon \widehat\alpha_n\|_{L_2(\mathbb{P}_n)}$ be the standard deviation estimator of our test statistic $\widecheck\eta_t^\textup{s}(\widecheck g_n)$ (with $\widecheck\eta_t^\textup{s}$ from \eqref{eta_t_s_wc}); specifically, $\|\widehat\epsilon \widehat\alpha_n\|_{L_2(\mathbb{P}_n)}= [ \frac{1}{n} \sum_{i=1}^n \{ \widehat\epsilon_i \widehat\alpha_n(\widetilde{\bm{F}}_i,X_{i,j}) \}^2 ]^{1/2}$ where $\widehat\epsilon_i=Y_i-\widecheck g_n(\widetilde{\bm{F}}_i,X_i)$ is the $i$-th residual. Then, for a prescribed significance level $0<\alpha<1$, reject the null hypothesis \eqref{hypo_null} in our conditional screening hypotheses: \vspace*{-0.15in}
\begin{itemize}\itemsep -0.15in
\item[a)]
For the fixed-$t$ test, if \eqref{test:fixed_t} holds;
\item[b)]
For the sup test, if the decision rule in \eqref{sup} holds;
\item[c)]
For the square test, if the decision rule in \eqref{square} hold. \vspace*{-0.2in}
\end{itemize}
\end{enumerate}
\label{algorithmhighd}
\hrulefill
\end{algorithm}
}
\section{Theoretical properties of estimators and tests}
\label{sec:regression_function}
\subsection{Initial regression and smoothed derivative estimators}
\label{sec:init_est_theory}
We first collect the assumptions most relevant to the construction and rates of our estimators $\widehat g_n$ and $\widehat g_{n,j}^\textup{s}$; we defer the more technical assumptions to Section~\ref{sec:additional_ass_sec_3}.
In particular, Assumption~\ref{ass:DGP} imposes mild conditions on the data generating processes, and Assumption~\ref{ass:NN_scaling} specifies the function class for $m_0$ and its associated neural networks.
\begin{assumption}[Data generating processes]\label{ass:DGP}
\begin{enumerate*}[label=(\roman*)]
\item\label{ass:DGP:con_1}
$({\bm{F}}, X_j)$ take value on the bounded support $[-b,b]^r\times[-b,b]$;
\item\label{ass:DGP:con_2}
{The dimension ${\overline{r}}$ of the diversified factor $\widetilde{\bm{F}}$ is a fixed constant and satisfies ${\overline{r}}\ge r$;}
\item\label{ass:DGP:con_3}
Regression function $m_0^*({\bm{f}},\bm{x})$, defined below \eqref{hypo_alt}, is bounded in magnitude by $M_\infty$ that satisfies $M_\infty\le M$ for $M$ the truncation level in \eqref{eq:base_NN_class}.
\end{enumerate*}
\end{assumption}
\begin{assumption}[Function class for $m_0$ and neural network scaling]
\label{ass:NN_scaling}
Regression function $m_0$ belongs to the hierarchical composition model $\mathcal{H}(r+1,l,\mathcal{P},C_\mathcal{H})$ in Definition~\ref{def:Hier_comp} \citep{schmidt2020nonparametric,KohlerLanger2021}. The tuple $(L, \mathbf{k}=({\overline{r}}+1, k_1, \ldots, k_L,1))$ for the structure of $\mathcal{F}_n({\overline{r}}+1)$ in \eqref{eq:base_NN_class} satisfies $k_1=\dots=k_L\equiv k_{\textup{0}}$ for some $k_{\textup{0}}\in\mathbb{N}$ and $L \cdot k_{\textup{0}} \sim n^{\frac{1}{4\kappa+2}} \log^{\frac{4\kappa-1}{2\kappa+1}}(n)$ where $\kappa$ is the dimension-adjusted degree of smoothness defined above Definition~\ref{def:Hier_comp}.
\end{assumption}
\begin{assumption}[Idiosyncratic terms]\label{asumidio}
The idiosyncratic terms ${\bm{u}}\in\mathbb{R}^d$ satisfy $\sum_{j=1}^d\mathbb{E}[u_j^2]\lesssim d_u$ for $d_u\leq d$ and a weak dependence condition $\sum_{j\neq{j}^{\prime}}|\mathbb{E}\left[u_{j}u_{j^{\prime}}\right]|\lesssim d_u$.
\end{assumption}
In Theorem~\ref{thm:m_hat_est_master} we establish the convergence of $\widehat g_n-m_0$ for regression function estimation and of $\widehat g_{n,j}^\textup{s} - m_{0,j}$ for derivative estimation. The following rates will appear:
\begin{align}
\label{eq:nu_n_concrete}
p_n = c_{\textup{p}} n^{\frac{1}{2\kappa+1}} \log^{2\frac{4\kappa-1}{2\kappa+1}+{2}}(n), \quad \nu_n = n^{-\frac{\kappa}{2\kappa+1}} \log^{\frac{6\kappa}{2\kappa+1}}(n) , \quad \delta_{\textup{f}}= \{{\overline{r}} d_u/d^2\}^{1/2}
\end{align}
where $c_{\textup{p}}$ is the constant appearing in \eqref{eq:V_calF_bound}, $d_u$ appears in Assumption~\ref{asumidio}, and $\delta_{\textup{f}} \to 0$ as $d\to \infty$. Moreover, because the target of the \textit{smoothed} derivative estimator $\widehat g_{n,j}^\textup{s}$ is the \textit{smoothed} derivative $m_{0,j}^\textup{s}$, not $m_{0,j}$, naturally a bias
\begin{align*}
\textstyle r_{\textup{b},m,j} \equiv \{ \int_{\Omega_h} ( m_{0,j}^\textup{s}-m_{0,j})^2\,{\rm d}\mathbb{P}\}^{1/2} = \left\{ \int_{\Omega_h} [ \int_{-h}^h \{ m_{0,j}({\bm{f}},x_j-z) - m_{0,j}({\bm{f}}, x_j)\} K_h(z) {\rm d} z ]^2{\rm d}\mathbb{P} \right\}^{1/2}
\end{align*}
occurs. However, this bias vanishes under the null in \eqref{hypo_null} because there $m_{0,j}=0$ everywhere.
\begin{theorem}
\label{thm:m_hat_est_master}
Under Assumptions~\ref{ass:DGP} to \ref{asumidio} and \ref{assum_fun_class} to \ref{ass:kernel}, for $p_n$, $\nu_n$ and $\delta_{\textup{f}}$ in \eqref{eq:nu_n_concrete}, on an event $\mathcal{A}_m'$ with $\mathbb{P}(\mathcal{A}_m') \ge 1 - C {\exp(- p_n )}$, and for a constant $c_{m,1}' $, the initial regression function estimator $\widehat g_n$ from \eqref{eq:def_wh_m_n_highd} satisfies
\begin{align}
\label{eq:rn1}
\textstyle \left[ \int \{ \widehat g_n(\widetilde{\bm{f}},x_j ) - m_0({\bm{f}}, x_j) \}^2 {\rm d}\mathbb{P} \right]^{1/2} \le c_{m,1}' ( \nu_n + \delta_{\textup{f}} ) \equiv r_{m,n} .
\end{align}
In addition, its associated smoothed derivative estimator $\widehat g_{n,j}^\textup{s}$ satisfies, on the same event $\mathcal{A}_m'$ and with a possibly different constant $c_{m,2}'$,
\begin{align}
\textstyle \left[ \int_{\Omega_h} \{ \widehat g_{n,j}^\textup{s}(\widetilde{\bm{f}}, x_j)- m_{0,j}({\bm{f}}, x_j) \}^2 {\rm d}\mathbb{P} \right]^{1/2} \le c_{m,2}' h^{-1} ( \nu_n + \delta_{\textup{f}} ) + r_{\textup{b},m,j} .
\label{eq:rn1_deriv}
\end{align}
\end{theorem}
The proof of Theorem~\ref{thm:m_hat_est_master} is deferred to Section~\ref{sec:proof_thm:m_hat_est_master}. With the rate $r_{m,n}$ in \eqref{eq:rn1} for the regression estimator, the rate \eqref{eq:rn1_deriv} for the smoothed \textit{derivative} estimator follows from Lemma~\ref{lemma:derivative_bound_via_original} that applies to the smoothed derivatives of all estimators. Next, the rate $r_{m,n}$ consists of two parts: $\nu_n$ represents the combined effect of stochastic error and neural network approximation bias, while $\delta_{\textup{f}}$ represents the error induced by relying on the diversified factors in place of the latent factors. Theorem~\ref{thm:m_hat_est_master} also covers as a special example the low-dimensional ${\bm{X}}$ case in Section~\ref{sec:lowd} by simply setting $\delta_{\textup{f}}=0$.
\subsection{Estimating the asymptotic variance}
\label{sec:variance}
\subsubsection{Definition of the score function}
\label{sec:variance_prep}
In this section, to complete our discussion in Section~\ref{sec:Riesz_est}, we formally introduce the score function $\alpha_{t,n}^*$ of our test statistics and the bias of $\widehat\alpha_n$ relative to $\alpha_{t,n}^*$. Define the population version of the operator $\widecheck\eta_t^\textup{s}(g)$ in \eqref{eta_t_s_wc} as
\begin{align}
\label{labeltildeg}
\textstyle \widetilde\eta_t^\textup{s}(g) = \int_{\Omega_h} [ \exp\{ t\,\Psi(g_j^\textup{s}(\widetilde{\bm{f}}, x_j)) \} - 1 ] {\rm d}\mathbb{P} ,
\end{align}
and its associated directional derivative (in the direction $v=v(\widetilde{\bm{f}},x_j)$)
\begin{align*}
\textstyle \frac{\partial \widetilde\eta_t^\textup{s}(g_0)}{\partial g}[v] = \frac{\partial \widetilde\eta_t^\textup{s}(g_0+\tau v)}{\partial\tau}\Big|_{\tau=0} .
\end{align*}
By \eqref{eq:target_functional_derivative_master} in the proof of Lemma~\ref{lem:In2_2_master}, $\frac{\partial \widetilde\eta_t^\textup{s}(g_0)}{\partial g}[v] / t$ is a bounded linear functional of $v$. Hence, there exists a \textit{Riesz representer} that becomes our $\alpha_{t,n}^*$ and that satisfies $\frac{\partial \widetilde\eta_t\left(g_0\right)}{\partial g}[v] / t = {\int v \alpha_{t,n}^* {\rm d}\mathbb{P}}$ for all $v$ \cite[Theorem~3.4, Chapter~1]{Conway1990}. In particular, by \eqref{eq:target_functional_derivative_master}, the effect of $\alpha_{t,n}^*$ is captured explicitly as
\begin{align}
& \textstyle \int v \alpha_{t,n}^* {\rm d}\mathbb{P} = \frac{1}{t} \frac{\partial \widetilde\eta_t^\textup{s}(g_0)}{\partial g}[v] = \int_{\Omega_h} \exp\{t \,g_{0,j}^\textup{s}(\widetilde{\bm{f}},x_j) \} v_j^\textup{s}(\widetilde{\bm{f}},x_j) {\rm d}\mathbb{P} \textstyle {\stackrel{\text{under}~H_{0}}{=}} \int_{\Omega_h} v_j^\textup{s}(\widetilde{\bm{f}},x_j) {\rm d}\mathbb{P}.
\label{eq:Riesz_H0_master}
\end{align}
Because the loss $\widehat{R}_n^{\textup{null}}$ relies on the last step of \eqref{eq:Riesz_H0_master}, $\widehat\alpha_n$ consistently estimates $\alpha_{t,n}^*$ under the null, as confirmed by Theorem~\ref{thm:Riesz_est}. However $\widehat\alpha_n$ may not be consistent for $\alpha_{t,n}^*$ under the alternative. Instead, it is not hard to show that
it is consistent for a population limit $\alpha_n^{\textup{null}}$ which exists and satisfies $\int v \alpha_n^{\textup{null}} {\rm d}\mathbb{P} = \int_{\Omega_h} v_j^\textup{s}(\widetilde{\bm{f}},x_j) {\rm d}\mathbb{P}$ for all $v$. See the remark below Eq.~\eqref{eq:target_functional_derivative_master} for details.
To obtain an intuitive idea of what the Riesz representer $\alpha_{t,n}^*$ could look like, assume that $\widetilde{\bm{F}}, X_j$ admit a joint density {$p(\widetilde{\bm{f}},x_j)$} that is differentiable in $x_j$ with the derivative being {$p_j(\widetilde{\bm{f}},x_j)$}. Then, by Lemma~\ref{lem:Riesz_analytic}, under the null hypothesis and in the limit $h\rightarrow 0$, $\alpha_{t,n}^*(\widetilde{\bm{f}},x_j)= - p_j(\widetilde{\bm{f}},x_j) / p(\widetilde{\bm{f}},x_j)$; in this case, if further $\widetilde{\bm{F}}, X_j$ are jointly Gaussian {on their support}, then $\alpha_{t,n}^*(\widetilde{\bm{f}},x_j)$ is simply a linear function in $\widetilde{\bm{f}}$ and $x_j$.
\subsubsection{Convergence rate for the score function estimator}
\label{sec:variance_thm}
To ensure the rate of $\widehat\alpha_n$, our Assumption~\ref{ass:u_t_n} mirrors our earlier Assumptions~\ref{ass:DGP} and \ref{ass:NN_scaling} for estimating the regression function $m_0$ with deep neural networks and in particular imposes a hierarchical composition model on $\alpha_n^{\textup{null}}$ that takes the diversified predictors $(\widetilde{\bm{F}},X_j)$ as arguments. Then, Assumption~\ref{ass:alpha_h} places a mild condition on the smoothing bandwidth.
As in Section~\ref{sec:init_est_theory}, we defer the more technical assumptions to Section~\ref{sec:additional_ass_sec_3}.
\begin{assumption}[Function class and neural network scaling for estimating $\alpha_n^{\textup{null}}$ and $\alpha_{t,n}^*$]\label{ass:u_t_n}
The function $\alpha_n^{\textup{null}}$ is bounded in magnitude by $M_\infty$, and when we restrict the support to $[-c_b b, c_b b]^{{\overline{r}}}\times[-b,b]$ for a constant $c_b>0$, belongs to the hierarchical composition model $\mathcal{H}({\overline{r}}+1,l,\mathcal{P},C_\mathcal{H})$ on the same support and (without loss of generality) with the same $l,\mathcal{P},C_\mathcal{H}$ as in Assumption~\ref{ass:NN_scaling}. The class $\mathcal{F}_n({\overline{r}}+1)$ in \eqref{eq:alpha_hat} satisfies the same structural scaling as in Assumption~\ref{ass:NN_scaling}.
\end{assumption}
\begin{assumption}[Rate of bandwidth]
\label{ass:alpha_h}
The bandwidth $h$ satisfies $h\ge 1/(\sqrt{n}\nu_n)$.
\end{assumption}
To characterize the bias of $\alpha_n^{\textup{null}}$ with respect to $\alpha_{t,n}^*$ under the alternative, define the random variable that represents the signal of the alternative as
\begin{align}
\label{eq:Z_t}
Z_{t,j,h} = [ \exp\{t\,g_{0,j}^\textup{s}(\widetilde{\bm{F}},X_j)\} - 1 ] \mathbbm{1}\{X_j\in\mathcal{B}_h\} .
\end{align}
Note that $Z_{t,j,h}$ always equals zero under the null hypothesis.
\begin{theorem}
\label{thm:Riesz_est}
Under Assumptions~\ref{ass:DGP} to \ref{ass:alpha_h} and \ref{assum_fun_class} to \ref{ass:truncation}, for $p_n$, $\nu_n$ and $\delta_{\textup{f}}$ in \eqref{eq:nu_n_concrete}, on an event $\mathcal{A}_\alpha$ with $\mathbb{P}(\mathcal{A}_\alpha) \ge 1 - C \exp(-p_n)$, and for a constant $c_{\alpha,1}$,
\begin{align}
&\textstyle \left[ \int (\widehat\alpha_n - \alpha_n^{\textup{null}})^2(\widetilde{\bm{f}},x_j ) {\rm d}\mathbb{P} \right]^{1/2} \le c_{\alpha,1} \{ ( h^{-1} \nu_n + \delta_{\textup{f}} ) \wedge M \} \equiv r_{\alpha,n}^{\textup{null}} .
\label{eq:r_n2_null}
\end{align}
Moreover, for a constant $c_{\alpha,2}$ not dependent on $t$, the bias of $\alpha_n^{\textup{null}}$ relative to $\alpha_{t,n}^*$ is
\begin{align}
\textstyle \forall t\in\mathbb{R}, \quad \left[ \int (\alpha_n^{\textup{null}} - \alpha_{t,n}^*)^2(\widetilde{\bm{f}},x_j ) {\rm d}\mathbb{P} \right]^{1/2} \le c_{\alpha,2} ( h^{-1} \| Z_{t,j,h} \|_{L_2} \wedge M ) \equiv r_{\alpha,t,\textup{b}} .
\label{eq:r_alpha}
\end{align}
Consequently, on the event $\mathcal{A}_\alpha$,
\begin{align}
&\textstyle \forall t\in\mathbb{R}, \quad \left[ \int (\widehat\alpha_n - \alpha_{t,n}^*)^2(\widetilde{\bm{f}},x_j ) {\rm d}\mathbb{P} \right]^{1/2} \le r_{\alpha,n}^{\textup{null}} + r_{\alpha,t,\textup{b}} \equiv r_{\alpha,t,n}.
\label{eq:r_n2}
\end{align}
\end{theorem}
The proof of Theorem~\ref{thm:Riesz_est} is deferred to Section~\ref{sec:proof_thm:Riesz_est}. Compared with the rate \eqref{eq:rn1} in Theorem~\ref{thm:m_hat_est_master} for the regression estimator $\widehat g_n$, the rates for the score function estimator $\widehat\alpha_n$ mainly differ in two aspects: first, a factor $h^{-1}$ precedes $\nu_n$ which is the consequence of the smoothing operation in the loss function $\widehat{R}_n^{\textup{null}}$ in \eqref{eq:Rhat}; second, under the alternative hypothesis a bias $h^{-1} \| Z_{t,j,h} \|_{L_2}$ in \eqref{eq:r_alpha} is induced relative to $\alpha_{t,n}^*$.
\subsection{Centering and adjusted estimator}
\label{sec:tweaking_proof}
Recall from Section~\ref{sec:tweaking} that the adjusted estimator $\widecheck g_n$ in \eqref{eq:f_check} was introduced to center our test statistics through a suitable minimization condition \eqref{eq:min_con}. In this section Proposition~\ref{prop:tweaking} first shows that a small $\widehat{\delta}_\textup{t}$ is sufficient to arrive at $\widecheck g_n$ that satisfies \eqref{eq:min_con} with a small $b_n$.
\begin{assumption}[Centering test statistics]
\label{ratexx}
$n$ is large enough such that
\begin{enumerate*}[label=(\roman*)]\label{ratexx:whole}
\item\label{ratexx:con_1}
$c_{\textup{t},1} (\delta_{\textup{f}} + \nu_n/h) \le 1$ for a large enough constant $c_{\textup{t},1}$;
\item\label{ratexx:con_2}
$c_{\textup{t},2} r_{m,n} \le 1$ for the constant $c_{\textup{t},2}$ in Proposition~\ref{prop:tweaking}; also, let {the} infinitesimal positive sequence $r_{\textup{inf},n}$ satisfy $r_{\textup{inf},n}={\scriptstyle{\mathcal{O}}}(r_{m,n} \nu_n)$.
\end{enumerate*}
\end{assumption}
\begin{proposition}
\label{prop:tweaking}
Under Assumptions~\ref{ass:DGP} to \ref{ratexx}\ref{ratexx:con_1} and \ref{assum_fun_class} to \ref{ass:truncation},
on an event $\mathcal{A}_{\textup{t},1}$ with $\mathbb{P}(\mathcal{A}_{\textup{t},1}) \ge 1 - C \exp(-p_n)$, for a constant $c_{\textup{t},2}$, $\widehat{\delta}_\textup{t}$ given by \eqref{eq:k_check} satisfies $|\widehat{\delta}_\textup{t}| \le c_{\textup{t},2} r_{m,n}$. If furthermore Assumption~\ref{ratexx}\ref{ratexx:con_2} holds, then on an event $\mathcal{A}_{\textup{t},2}$ with $\mathbb{P}(\mathcal{A}_{\textup{t},2}) \ge 1 - C \exp(-p_n)$,
for a constant $c_{\textup{t},3}$, condition \eqref{eq:min_con} holds uniformly at all $t\in\mathcal{T}_\delta$ and $b_n=c_{\textup{t},3} r_{m,n} (r_{\alpha,n}^{\textup{null}} + r_{\alpha,t,\textup{b}})$.
\end{proposition}
The proof of Proposition~\ref{prop:tweaking} is deferred to Section~\ref{sec:proof_prop:tweaking}. In the proposition, under the null the tolerance $b_n=r_{m,n} r_{\alpha,n}^{\textup{null}}$, which is faster than $n^{-1/2}$ under appropriate conditions and will ensure the asymptotic normality of our test statistics. Under the alternative, $b_n$ is not necessarily faster than $n^{-1/2}$, but is still fast enough to ensure {consistency under the local alternatives}.
We will start working mostly with the adjusted estimator $\widecheck g_n$ from now on. The next theorem is the counterpart of Theorem~\ref{thm:m_hat_est_master} and shows that $\widecheck g_n$ and its smoothed derivative estimator $\widecheck g_{n,j}^\textup{s}$ maintain convergence rates similar to the unadjusted ones.
\begin{theorem}
\label{thm:m_check_est_master}
Under Assumptions~\ref{ass:DGP} to \ref{ratexx} and \ref{assum_fun_class} to \ref{ass:truncation}, for $p_n$, $\nu_n$ and $\delta_{\textup{f}}$ in \eqref{eq:nu_n_concrete}, on an event $\mathcal{A}_m$ with $\mathbb{P}(\mathcal{A}_m) \ge 1 - C \exp(-p_n)$, and for a constant $c_{m,1}$,
\begin{align}
\label{eq:m_check_rate_master}
\textstyle \left[ \int \{ \widecheck g_n(\widetilde{\bm{f}},x_j ) - m_0({\bm{f}},\bm{x} ) \}^2 {\rm d}\mathbb{P} \right]^{1/2} \le c_{m,1} ( \nu_n + \delta_{\textup{f}} ) .
\end{align}
In addition, its smoothed derivative estimator $\widecheck g_{n,j}^\textup{s}$ satisfies, on the same event $\mathcal{A}_m$ and with a possibly different constant $c_{m,2}$,
\begin{align}
\textstyle \left[ \int_{\Omega_h} \{ \widecheck g_{n,j}^\textup{s}(\widetilde{\bm{f}}, x_j)- m_{0,j}({\bm{f}}, x_j) \}^2 {\rm d}\mathbb{P} \right]^{1/2} \le c_{m,2} h^{-1} ( \nu_n + \delta_{\textup{f}} ) + r_{\textup{b},m,j} .
\label{eq:m_check_rate_master_deriv}
\end{align}
\end{theorem}
The proof of Theorem~\ref{thm:m_check_est_master} is deferred to Section~\ref{sec:proof_thm:f_check_est}.
\subsection{Properties of the conditional screening tests}
\label{sec:test}
This section gives the result on the asymptotic null distributions of the test statistics. The results for consistency against the local alternatives are given in Section~\ref{sec:alter_supp}.
\begin{theorem}
\label{thm:main_null_master}
Suppose that Assumptions~\ref{ass:DGP} to \ref{ratexx} and \ref{assum_fun_class} to \ref{ass:truncation} hold, and in addition condition
\begin{enumerate*}[label=($\ast$)]
\item\label{thm:main_null_master:con_1}
holds: $\delta_{\textup{f}} + h^{-2} (\nu_n^2+\delta_{\textup{f}}^2) = {\scriptstyle{\mathcal{O}}}(n^{-1/2})$.
\end{enumerate*}
Then, with $\mathcal{T}_\delta$ and tests given in Section~\ref{sec:uniform}, we have
\begin{align}
\label{eq:H0_normality_L2}
\text{for all fixed}~t\in\mathcal{T}_\delta, \frac{ \sqrt{n} }{t} \frac{\widecheck\eta_t^\textup{s}(\widecheck g_n) }{\| \epsilon \alpha_n^{\textup{null}}\|_{L_2} }\rightarrow_{d}\mathcal{Z} , \quad \frac{\widehat Z}{\| \epsilon \alpha_n^{\textup{null}}\|_{L_2}} \rightarrow_{d} |\mathcal{Z}|, \quad \frac{\widehat\chi^2}{\| \epsilon \alpha_n^{\textup{null}}\|_{L_2}^2} \rightarrow_{d} \chi_1^2 ,
\end{align}
where $\mathcal{Z}$ stands for a standard normal random variable and $\chi_1^2$ is a chi-square random variable with one degree of freedom. Moreover, we are free to replace $\| \epsilon \alpha_n^{\textup{null}}\|_{L_2}$ in \eqref{eq:H0_normality_L2} above by its empirical counterpart $\|\widehat\epsilon \widehat\alpha_n\|_{L_2(\mathbb{P}_n)}$ from Section~\ref{sec:uniform}, to conclude that
\begin{align}
\label{eq:H0_normality_Ln}
\text{for all fixed}~t\in\mathcal{T}_\delta, \frac{ \sqrt{n} }{t} \frac{ \widecheck\eta_t^\textup{s}(\widecheck g_n) }{\| \widehat\epsilon \widehat\alpha_n \|_{L_2(\mathbb{P}_n)} }\rightarrow_{d} \mathcal{Z} ,~~\frac{\widehat Z}{ \| \widehat\epsilon \widehat\alpha_n \|_{L_2(\mathbb{P}_n)} } \rightarrow_{d} |\mathcal{Z}|,~~\frac{\widehat\chi^2}{ \| \widehat\epsilon \widehat\alpha_n \|_{L_2(\mathbb{P}_n)}^2 } \rightarrow_{d} \chi_1^2 .
\end{align}
\end{theorem}
The proof of Theorem~\ref{thm:main_null_master} is deferred to Section~\ref{sec:Proof_Thm_thm:main_null_master}. The extra condition~\ref{thm:main_null_master:con_1} in Theorem~\ref{thm:main_null_master} is natural due to the presence of the $\sqrt{n}$ scaling factor in the our tests and is mild. Note that if the bandwidth $h$ is held as a constant, then condition~\ref{thm:main_null_master:con_1} simply reduces to $\delta_{\textup{f}}={\scriptstyle{\mathcal{O}}}(n^{-1/2})$, and $\nu_n={\scriptstyle{\mathcal{O}}}(n^{-1/4})$ which is in turn implied by $\kappa>1/2$ for $\kappa$ in Assumption~\ref{ass:NN_scaling} on the hierarchical composition model.
\section{Simulation studies}\label{sec:simulations}
In this section, we conduct simulations to assess the size and power of our conditional screening tests across different scenarios. Recall that we summarize the implementations of these tests in Algorithm~\ref{algorithmhighd}. We consistently employ the Quartic/biweight kernel.
We choose $t=1$ for the fixed-$t$ test statistic, and simply set $\mathcal{T}_\delta = \{-1.25,-0.5,0.5,1.25\}$ in the sup statistic and the square statistic in \eqref{sup} and \eqref{square} respectively. Then, we implement the sup statistic in \eqref{sup} as
$\mbox{sup}_{t \in T_{\delta}}| {\widehat\eta_t^\textup{s}(\widecheck g_n)} |/t$, and the square statistic as $\widehat\chi^2_{\delta} = n \sum_{t\in \mathcal{T}_\delta} \left\{ {\widehat\eta_{t}^\textup{s}(\widecheck g_n)} \right\}^2 / t^2$ (we take $w(t)=1$ in \eqref{square}).
We set $d=200$ and $d=400$ when $n=256$ and $n=512$ respectively, and the full regression model as $Y=m_0^*({\bm{F}},X_1,X_3)+\epsilon^*$ (see $m_0^*$ defined below \eqref{hypo_alt}). To illustrate the performance of our conditional screening test, we have designed $m_0^*$ so that a significant portion of its variation is accounted for by the factors.
We further consider a nonlinear and a linear model of $m_0^*$. Specifically, under the null hypothesis, we set:
\begin{equation}
\begin{gathered}
\mbox{nonlinear}: m_0^*(\cdot) = m_0^{*(null)}(\cdot) = \mbox{sin}(f_1+u_1) + \log(8+f_2) \times \log(8+f_3) + \exp(-f_4^2/2), \\
\mbox{linear}: m_0^*(\cdot) = m_0^{*(null)}(\cdot) = f_1-f_2+f_3+f_4-f_5 .
\end{gathered}
\label{eq:sim_null}
\end{equation}
For our screening test, we always select $X_3$ as the variable of interest. Accordingly, under the alternative hypothesis we add signals in $X_3$ to $m_0^{*(null)}$ above:
\begin{equation}
\begin{gathered}
\mbox{nonlinear}: m_0^*(\cdot) = m_0^{*(null)}(\cdot) + X_3^2 /4, \quad
\mbox{linear}: m_0^*(\cdot) = m_0^{*(null)}(\cdot)+ X_3/16.
\end{gathered}
\label{eq:sim_alt}
\end{equation}
The signal under the linear model is very weak and allows us to discern different test settings.
While the signal under the nonlinear model may not appear weak at first, here the \textit{average} derivative $\mathbb{E} m_{0,j}=0$, so detecting departure from the null critically depends on the higher-order, nonlinear effect of our MGF test statistics. Under the nonlinear model $r=4$ while under the linear model $r=5$. In both models we set ${\overline{r}}=r$.
Under the nonlinear model, $u_1$ is incorporated for a richer structure that specifically leads to the working regression function $m_0({\bm{F}},X_3)$ being different from the full regression function $m_0^*({\bm{F}},X_1,X_3)$ by design; see Remark~\ref{rmk:misspecification}. We incorporate $u_1 = X_1 - \mathbf{B}_{1\cdot}{\bm{F}}$ instead of $X_1$ directly because $u_1$ will be enforced to be independent of $X_3$ under screening. Moreover, under the null hypothesis, the selection of $X_3$ is purely for clarity: under the nonlinear model, we could conduct our test on any $X_j$ for any $j\in\{2,\dots,d\}$, and under the linear model, for any $j\in\{1,\dots,d\}$.
For the high-dimensional predictor ${\bm{X}}$ in the factor model~\eqref{eq:factor_observation}, the factor loading matrix $\mathbf{B}$ is generated with i.i.d.\,Unif$[-\sqrt{3},\sqrt{3}]$ entries; the factor ${\bm{F}}$ and the idiosyncratic terms ${\bm{u}}$ both have i.i.d.~$\mathcal{N}(0,0.6)$ entries. We further draw the noise as $\epsilon^*\sim\mathcal{N}(0,0.3)$. The quantities $\mathbf{B}$, ${\bm{F}}$, ${\bm{u}}$ and $\epsilon^*$ are all drawn independently. We pre-train $\mathbf{W}$ with samples of size $100$.
We set the hyper-parameters for the neural network fitting as follows: for the regression estimator $\widehat g_n$ in \eqref{eq:def_wh_m_n_highd}, we employ neural networks with $L=5$ hidden layers, a common width of $k_{\textup{0}}=16$ per layer, and the ReLU activation function. Training proceeds over $800$ epochs, utilizing a batch size of $256$ and a constant learning rate of $0.005$. For the score function estimator $\widehat\alpha_n$ in \eqref{eq:alpha_hat}, we maintain the same neural network fitting parameters, except that we set $L=2$, batch size to $64$ and the number of epochs to $400$. We employ early stopping as the only regularization technique, and terminate training if there's no improvement after 20 epochs ($\text{patience}=20$) on a validation set.
We conduct our tests at significance levels of either $5\%$ or $10\%$ under both the nonlinear and the linear models, resulting in a total of four combinations summarized in Tables~\ref{tab:farlog1} to \ref{tab:farlinear2}. In each table we present in alternating rows the sizes and powers (the latters in parentheses) of the tests in \eqref{test:fixed_t}, \eqref{sup} and \eqref{square}, and of the same tests but without centering (that is, tests employing the non-adjusted estimator $\widehat g_n$), under different sample sizes. The columns, arranged from left to right, correspond to the four different bandwidths. Each entry is calculated based on $500$ Monte Carlo repetitions.
\begin{table}[htbp]
\centering
\begin{tabular}{c| c|c|c|c|c|c}
\hline\hline
\multirow{9}{*}{$n=256$}& & $h$ & 0.1 & 1.0 & 1.5 & 2.0 \\ \hline
& \multirow{2}{*}{fixed-$t$ test} & Non-centered & 0.16 (0.67) & 0.14 (0.57) & 0.13 (0.56) & 0.12 (0.51) \\
& & Centered & 0.11 (0.69) & 0.10 (0.60) & 0.08 (0.55) & 0.07 (0.50) \\ \cline{2-7}
& \multirow{2}{*}{sup test} & Non-centered & 0.26 (0.97) & 0.21 (0.93) & 0.19 (0.92) & 0.18 (0.87) \\
& &Centered & 0.19 (0.97) & 0.16 (0.93) & 0.13 (0.91) & 0.13 (0.85) \\ \cline{2-7}
& \multirow{2}{*}{square test} & Non-centered & 0.16 (0.86) & 0.14 (0.79) & 0.14 (0.74) & 0.12 (0.65) \\
& & Centered & 0.11 (0.82) & 0.09 (0.78) & 0.09 (0.70) & 0.09 (0.63) \\
\hline
\end{tabular}
\begin{tabular}{c|c|c|c|c|c|c}
\multirow{7}{*}{$n=512$} & \multirow{2}{*}{fixed-$t$ test} & Non-centered & 0.11 (0.88) & 0.10 (0.80) & 0.10 (0.77) & 0.08 (0.70) \\
& & Centered & 0.09 (0.88) & 0.07 (0.83) & 0.07 (0.76) & 0.05 (0.70) \\ \cline{2-7}
& \multirow{2}{*}{sup test} & Non-centered & 0.19 (1.00) & 0.16 (0.98) & 0.15 (0.97) & 0.11 (0.96) \\
& & Centered & 0.15 (1.00) & 0.12 (0.97) & 0.11 (0.97) & 0.08 (0.95) \\ \cline{2-7}
& \multirow{2}{*}{square test} & Non-centered & 0.12 (0.96) & 0.10 (0.91) & 0.09 (0.87) & 0.08 (0.85) \\
& & Centered & 0.08 (0.96) & 0.07 (0.91) & 0.05 (0.87) & 0.05 (0.86) \\ \hline \hline
\end{tabular}
\caption{Performance summary of our test statistics under the nonlinear model in \eqref{eq:sim_null} (for size under the null) and \eqref{eq:sim_alt} (for power under the alternative) at the significance level $5\%$. Specifically, we provide the size and power (the latter displayed in parentheses) of our tests under various combinations of the test statistic (fixed-$t$, sup, or squared statistic), sample size, bandwidth, and the use of either the non-centered (employing the non-adjusted estimator $\widehat g_n$) or the centered (employing the adjusted $\widecheck g_n$) test statistics. Each value in the table represents the average over $500$ Monte Carlo repetitions.}
\label{tab:farlog1}
\end{table}
\begin{table}[htbp]
\centering
\begin{tabular}{c|c|c|c|c|c|c}\hline\hline
\multirow{9}{*}{$n=512$}& & $h$ & 0.1 & 1.0 & 1.5 & 2.0 \\ \hline
& \multirow{2}{*}{fixed-$t$ test} & Non-centered & 0.16 (0.84) & 0.14 (0.78) & 0.14 (0.76) & 0.14 (0.76) \\
& & Centered & 0.08 (0.87) & 0.09 (0.82) & 0.08 (0.80) & 0.07 (0.79) \\
\cline{2-7}
& \multirow{2}{*}{sup test} & Non-centered & 0.25 (0.85) & 0.21 (0.79) & 0.18 (0.77) & 0.18 (0.76) \\
& & Centered & 0.18 (0.88) & 0.13 (0.83) & 0.11 (0.80) & 0.11 (0.80) \\
\cline{2-7}
& \multirow{2}{*}{square test} & Non-centered & 0.14 (0.74) & 0.14 (0.72) & 0.15 (0.72) & 0.14 (0.72) \\
& & Centered & 0.09 (0.77) & 0.09 (0.75) & 0.08 (0.73) & 0.07 (0.75) \\
\hline\hline
\end{tabular}
\caption{ Performance summary of our test statistics under the linear model in \eqref{eq:sim_null} (for size under the null) and \eqref{eq:sim_alt} (for power under the alternative) at the significance level $5\%$.}
\label{tab:farlinear1}
\end{table}
\begin{table}[htbp]
\centering
\begin{tabular}{c| c|c|c|c|c|c}
\hline\hline
\multirow{9}{*}{$n=512$}& & $h$ & 0.1 & 1.0 & 1.5 & 2.0 \\ \hline
& \multirow{2}{*}{fixed-$t$ test} & Non-centered & 0.17 (0.90) & 0.18 (0.84) & 0.18 (0.81) & 0.14 (0.76) \\
& & Centered & 0.13 (0.91) & 0.13 (0.86) & 0.13 (0.81) & 0.13 (0.77) \\ \cline{2-7}
& \multirow{2}{*}{sup test} & Non-centered & 0.27 (1.00) & 0.24 (0.99) & 0.24 (0.98) & 0.20 (0.97) \\
& & Centered & 0.24 (1.00) & 0.19 (0.99) & 0.19 (0.98) & 0.17 (0.98) \\
\cline{2-7}
& \multirow{2}{*}{square test} & Non-centered & 0.19 (0.99) & 0.18 (0.96) & 0.16 (0.93) & 0.13 (0.90) \\
&& Centered & 0.13 (0.98) & 0.12 (0.95) & 0.12 (0.92) & 0.11 (0.90) \\
\hline\hline
\end{tabular}
\caption{ The same caption as Table~\ref{tab:farlog1} except for significant level 10\%.}
\label{tab:farlog2}
\end{table}
\begin{table}[htbp]
\centering
\begin{tabular}{c|c|c|c|c|c|c}\hline\hline
\multirow{9}{*}{$n=512$}& & $h$ & 0.1 & 1.0 & 1.5 & 2.0 \\ \hline
& \multirow{2}{*}{fixed-$t$ test} &
Non-centered & 0.23 (0.89) & 0.21 (0.85) & 0.20 (0.82) & 0.21 (0.82) \\
& & Centered & 0.14 (0.93) & 0.16 (0.88) & 0.13 (0.87) & 0.12 (0.85) \\ \cline{2-7}
& \multirow{2}{*}{sup test} & Non-centered & 0.34 (0.90) & 0.26 (0.86) & 0.26 (0.82) & 0.26 (0.83) \\
&& Centered & 0.26 (0.93) & 0.21 (0.89) & 0.17 (0.88) & 0.16 (0.86) \\\cline{2-7}
& \multirow{2}{*}{square test} & Non-centered & 0.21 (0.83) & 0.21 (0.80) & 0.20 (0.78) & 0.21 (0.79) \\
& & Centered & 0.15 (0.86) & 0.15 (0.84) & 0.12 (0.83) & 0.13 (0.82) \\
\hline\hline
\end{tabular}
\caption{ The same caption as Table~\ref{tab:farlinear1} except for significant level 10\%.}
\label{tab:farlinear2}
\end{table}
We make a few observations from the results in the tables. First, centering generally improves the size of our tests without affecting their power. Moreover, under the null, sizes arrive at the nominal level if we further increase the bandwidth; this is especially evident for the sup test, which naturally tends to reject more often (under both the null and alternative). This shows the improvement of the power of the centered, smoothed test compared to the initial test \eqref{eq:eta_t_naive} if both tests are calibrated to have the same size. Of course, if we increase $h$ more than $1.5$, we expect the biases will arrive {under the alternative}. In the nonlinear case (Tables~\ref{tab:farlog1} and \ref{tab:farlog2}), the ensemble tests by aggregating various $t$ also outperform the fixed-$t$ test in power.
\section{Empirical Applications}
\label{sec:emp_app}
\subsection{Asset Pricing}
\label{sec:asset_pricing}
In our first empirical analysis, we examine a comprehensive dataset containing both returns and specific characteristics of firms. Twelve monthly returns in $2021$ for the 100 largest financial institutions are drawn from the CRSP database as our response variable, along with a set of $d=48$ firm-specific characteristics as predictors; the resulting sample size is {$n=1200$}. The dataset's foundation is credited to \cite{chen2021open}.
We select the smoothing bandwidth from the sequence $\{0.4, 0.6, \dots, 1.8, 2.0\}$ using the bandwidth selection algorithm suggested in Remark~\ref{rmk:bandwidth}. We set ${\overline{r}}=5$ by the conventional practice in finance, the hidden size $k_{\textup{0}}=24$, and the learning rate $\tau=10^{-3}$; we keep the other neural network fitting parameters identical to Section~\ref{sec:simulations}.
In Figure~\ref{assetpricing} we plot for each of the $44$ continuous predictors the intervals whose right ends are the test statistics normalized by their corresponding critical values in \eqref{test:fixed_t}, \eqref{sup} and \eqref{square}. Specifically, we plot the intervals
$$ [ 0, n^{1/2} |\widecheck\eta_t^\textup{s}(\widecheck g_n)| / ( t z_{1-\alpha/2} \|\widehat\epsilon \widehat\alpha_n\|_{L_2(\mathbb{P}_n)} ) ], \quad [ 0, \widehat Z/(z_{1-\alpha/2}\|\widehat\epsilon \widehat\alpha_n\|_{L_2(\mathbb{P}_n)})], \quad [ 0, \sqrt{\widehat\chi^2}/(\sqrt{\chi^2_{1,1-\alpha}}\|\widehat\epsilon \widehat\alpha_n\|_{L_2(\mathbb{P}_n)})],$$
for the fixed-$t$, sup, and square tests respectively. We also plot the same intervals but for the non-centered versions of these tests. In all cases, an interval covering one indicates statistical significance. The results suggest that the presence of some idiosyncratic terms related to tail risk, such as ReturnSkewCAPM and DownsideBeta, exert additional influence on returns and complement the explanatory power of the five factors.
\begin{figure}
\includegraphics[height=20cm, width = 1.0\textwidth]{Figs/CI_RET.pdf}
\caption{Asset Pricing Dataset: Significance of {predictors/idiosyncratic terms} at the $5\%$ significance level. Intervals covering one (indicated by the vertical black line) correspond to the significant {predictors/idiosyncratic terms}. For each variable, six intervals are plotted in the order of: fix-$t$ test, sup test and square test, and within each test the {\color{black}non-centered version is plotted before the centered version}.\color{black}}
\label{assetpricing}
\end{figure}
\subsection{Macroeconomics Time Series}
\label{sec:applications:FRED-MD}
In this section, to substantiate Section~\ref{sec:asset_pricing}, we illustrate our conditional screening test with another empirical application, this time on the macroeconomics dataset FRED-MD introduced in \cite{mccracken2016fred} and also later studied in \cite{FanGu2023factor}. The dataset collects $d=127$ monthly U.S.\,macroeconomic variables, such as unemployment rate and real personal income, starting from 1959/01. It is shown in \cite{mccracken2016fred} that these variables can be explained well by several latent factors.
Our analysis setup is in general similar to Section~C in \cite{FanGu2023factor}. Our target variables are UEMP15T26, TB3SMFFM or TB6SMFFM, and we aim to identify which variables contribute to predicting the target variables beyond the latent factors. The variable UEMP15T26 represents the civilians unemployed for $15-26$ weeks. The variable TB3SMFFM (TB6SMFFM) measures the 3-month (6-month) treasury bill rate minus the effective federal funds rate. For each target response variable $y_{t+1}$ in \{UEMP15T26, TB3SMFFM, TB6SMFFM\}, we regress $y_{t+1}$ on $\bm{x}_t\in\mathbb{R}^{127}$ where $\bm{x}_t$ is the vector of all variables at the previous month. We choose the $n=330$ valid sample pairs $\{\bm{x}_t, y_{t+1}\}$ between January 1980 and July 2022. We employ the same neural network fitting parameters as Section~\ref{sec:asset_pricing} except that we revert back to a hidden size $k_0=16$ as in our simulation studies.
For each variable, we plot the same six intervals derived from our tests as in Section~\ref{sec:asset_pricing} and we recall that an interval covering one indicates statistical significance.
Our analysis, based on the variables categorized by \cite{mccracken2016fred}, underscores the predominant impact of a limited set of variables on the three response variables. As there are 127 variables, it takes four sub-figures to display a figure that is similar Figure~\ref{assetpricing} for each given response variable.
Figures~\ref{FRED10} {to} \ref{FRED13} in the appendix display the intervals, with the response variable being UEMP15T26 (civilian unemployment).
Notably, most variables fail to reach the $5\%$ significance level, except for some variables such as RETAILx, NDMANEMP, and COMPAPFFx. Additional significant variables for the response variables TB3SMFFM and TB6SMFFM are illustrated from Figures~\ref{FRED1} to \ref{FRED8}. As for TB3SMFFM, there are many significant variables. For example, the variables such as IPDMAT in the output and income variable group and DTCTHFNM in the prices group contribute additionally to the factors. TB6SMFFM demonstrates slightly higher susceptibility to variables compared to TB3SMFFM, which, in turn, is notably influenced by variables related to the stock market, employment, interest rates, and money and credit. Our overall findings validate our regression model based on a latent factors plus sparse idiosyncratic terms structure.
\section{Conclusion and further work}
\label{sec:conclusion}
We have introduced a conditional variable screening test for non-parametric regression using deep neural networks; the inputs to the networks are obtained with the help of a factor model that further enables us to handle very high-dimensional predictors. In our test statistics, we employ high-quality estimators of the partial derivatives of the non-parametric regression function, which could be of independent interest. To demonstrate the versatility of our test, we apply it to assess the adequacy of non-parametric factor regression. An intriguing avenue for further exploration involves extending this framework to dependent data, other statistical machine learning losses, and simultaneous testing for multiple variables. Relevant examples for the latter direction include the $\ell_{\infty}$ statistics proposed by \cite{chen2022inference} and the $\ell_2-\ell_{\infty}$ type statistics discussed by \cite{li2022ell}.
This paper focuses on the popular feed-forward networks with the ReLU activation function. The derivative irregularity is not limited to the ReLU activation function, and hence, addressing this issue can benefit various neural network classes. Moreover, in our approach, we apply the smoothing procedure after we have obtained a neural network estimator. This is different from other derivative regularizations, for instance, Sobolev training \citep{Czarneck2017sobolev}, that employ roughness penalties during optimization, which could form a potential future topic. Last but not least, we could also potentially conduct different tests for testing higher- or mixed-derivatives. For instance, testing the monotonicity in a variable in the regression is equivalent to testing the sign of the associated partial derivative, a task already studied in the context of an one-dimensional predictor by, for instance, \cite{BowmanJonesGijbels1998testing}, \cite{GhosalArusharka2000testing} and \cite{HallHuang2001nonparametric}.
{\footnotesize \spacingset{1.0}
\bibliographystyle{apalike}
\bibliography{biball}
}