# Test of Significance for High-dimensional Thresholds with Application to
Individualized Minimal Clinically Important Difference

Huijie Feng Yang Ning Jiwei Zhao Department of Statistics and Data Science,
Cornell University, Ithaca, NY 14850, USA; e-mail:
<EMAIL_ADDRESS>of Statistics and Data Science, Cornell
University, Ithaca, NY 14850, USA; e-mail<EMAIL_ADDRESS>of
Biostatistics and Medical Informatics, University of Wisconsin-Madison,
Madison, WI 53726, USA; e-mail<EMAIL_ADDRESS>
###### Abstract

This work is motivated by learning the individualized minimal clinically
important difference, a vital concept to assess clinical importance in various
biomedical studies. We formulate the scientific question into a high-
dimensional statistical problem where the parameter of interest lies in an
individualized linear threshold. The goal of this paper is to develop a
hypothesis testing procedure for the significance of a single element in this
high-dimensional parameter as well as for the significance of a linear
combination of this parameter. The difficulty dues to the high-dimensionality
of the nuisance component in developing such a testing procedure, and also
stems from the fact that this high-dimensional threshold model is nonregular
and the limiting distribution of the corresponding estimator is nonstandard.
To deal with these challenges, we construct a test statistic via a new bias
corrected smoothed decorrelated score approach, and establish its asymptotic
distributions under both the null and local alternative hypotheses. In
addition, we propose a double-smoothing approach to select the optimal
bandwidth parameter in our test statistic and provide theoretical guarantees
for the selected bandwidth. We conduct comprehensive simulation studies to
demonstrate how our proposed procedure can be applied in empirical studies.
Finally, we apply the proposed method to a clinical trial where the scientific
goal is to assess the clinical importance of a surgery procedure.

Keyword: Bandwidth selection, High-dimensional statistical inference, Kernel
method, Nonstandard asymptotics.

##  1 Introduction

###  1.1 Motivation: Individualized Minimal Clinically Important Difference
(iMCID) under High-dimensionality

In clinical studies, the effect of a treatment or intervention is widely
assessed through clinical significance, instead of statistical significance.
By leveraging patient-reported outcomes (PRO) that are directly collected from
the patients without a third party’s interpretation, the aim of assessing
clinical significance is to provide clinicians and policy makers the clinical
effectiveness of the treatment or intervention. For example, in our motivating
study, the ChAMP randomized controlled trial (Bisson et al., 2015), the
interest is to identify _the smallest WOMAC pain score change_ such that the
corresponding improvement and beyond can be claimed as clinically significant.
In Jaeschke et al. (1989), this concept was firstly and formally introduced as
the minimal clinically important difference (MCID), “the smallest difference
in score in the domain of interest which patients perceive as beneficial and
which would mandate a change in the patient’s management”.

There are roughly three approaches to determine the magnitude of MCID (Lassere
et al., 2001; Erdogan et al., 2016; Angst et al., 2017; Jayadevappa et al.,
2017): distribution-based, opinion-based, and anchor-based. Although adopted
in various studies (Wyrwich et al., 1999a, b; Samsa et al., 1999; Norman et
al., 2003; Bellamy et al., 2001), the first two approaches are usually
criticized (McGlothlin and Lewis, 2014) due to, for example, “distribution-
based methods are not derived from individual patients” and “expert opinion
may not be a valid and reliable way to determine what is important to
patients”. The third approach, anchor-based, conceptually determines the MCID
by incorporating both certainty of effective treatment encoded as a continuous
variable and the patient’s satisfaction collected from the anchor question. It
is clinically evident (Wells et al., 2001) that the magnitude of MCID would
depend on various factors such as the demographic variables and the patients’
baseline status. For example, in a shoulder pain reduction study (Heald et
al., 1997), because of the higher expectation for complete recovery, the
healthier patients with mild pain at baseline often deemed greater pain
reduction as “meaningful” than the ones who suffered from chronic disease.
Therefore, it is of scientific interest to generally estimate the
individualized MCID (iMCID) based on each individual patient’s clinical
profile as well as to quantify the uncertainties of those estimates.

Nowadays, there is an increasing use and advancing development of EHR-based
(electronic health records) studies in clinical research. The EHR data are
complex, diverse and high-dimensional (Abdullah et al., 2020). The rich
information contained in the EHR data could facilitate the determination and
quantification of iMCID. Therefore, there is a pressing need to develop
statistical methods that incorporate the high-dimensional data into both
magnitude determination and uncertainty quantification of iMCID.

###  1.2 Problem Formulation

To facilitate the presentation, we first introduce some notation. Let
$X\in\mathbb{R}$ be a continuous variable representing the score change
collected from the PRO, e.g., the WOMAC pain score change from baseline to one
year after surgery in the ChAMP trial. Let $Y=\pm 1$ be a binary variable
derived from the patient’s response to the anchor question, where $Y=1$
represents an improved health condition and $Y=-1$ otherwise. We use a
$d$-dimensional vector $\bm{Z}$ to denote the patient’s clinical profile
including demographic variables, clinical biomarkers, disease histories, among
many others. Suppose the data we observe are $n$ i.i.d. samples
$\\{(x_{i},y_{i},\bm{z}_{i})\\}_{i=1}^{n}$ of $(X,Y,Z)$. We focus on the high-
dimensional setting, i.e., $d\gg n$.

Firstly, if there were no covariate $\bm{Z}$, the MCID can be estimated by
$\mathop{\mathrm{argmax}}_{\tau}\\{\mathbb{P}(X\geq\tau\mid
Y=1)+\mathbb{P}(X<\tau\mid Y=-1)\\}$, which is equivalent to

$\displaystyle\mathop{\mathrm{argmin}}_{\tau}~{}\mathbb{E}[w(Y)L_{01}\\{Y(X-\tau)\\}],$
(1.1)

where $L_{01}(u)=\frac{1}{2}\\{1-\text{sign}(u)\\}$ is the 0-1 loss,
$\text{sign}(u)=1$ if $u\geq 0$ and $-1$ otherwise, $w(1)=1/\pi$,
$w(-1)=1/(1-\pi)$ and $\pi=\mathbb{P}(Y=1)$. When the high-dimensional
covariate $\bm{Z}$ is available, as the focus of this paper, the natural idea
is to consider the iMCID with a functional form of $\bm{Z}$, say
$\tau(\bm{Z})$. In clinical practice, a simple structure, such as linear, is
preferred due to its transparency and convenience for interpretation,
especially for high-dimensional data. Therefore, we focus on the linear
structure $\tau(\bm{Z})=\bm{\beta}^{T}\bm{Z}$ in this paper. The objective
becomes

$\displaystyle\bm{\beta}^{*}=\mathop{\mathrm{argmin}}_{\bm{\beta}}R(\bm{\beta}),~{}~{}\textrm{where}~{}~{}R(\bm{\beta})=\mathbb{E}\big{[}w(Y)L_{01}\\{Y(X-\bm{\beta}^{T}\bm{Z})\\}\big{]}.$
(1.2)

Throughout this paper we assume that $\bm{\beta}^{*}$ exists and is unique.
Denote $\bm{\beta}^{*}=(\theta^{*},\bm{\gamma}^{*T})^{T}$, where $\theta^{*}$
is an arbitrary one-dimensional component of $\bm{\beta}^{*}$ and
$\bm{\gamma}^{*}$ represents the rest of the parameter which is high-
dimensional. In this paper, we start from considering the hypothesis testing
procedure for the parameter $\theta^{*}$. With a simple reparametrization, the
same procedure can be applied to infer the iMCID
$\bm{c}_{0}^{T}\bm{\beta}^{*}$ for some fixed and known vector
$\bm{c}_{0}\in\mathbb{R}^{d}$.

It is also worthwhile to mention that, although the motivation of this paper
is to study iMCID, our formulation of this problem can be similarly applied to
other scenarios as well, such as the covariate-adjusted Youden index (Xu et
al., 2014), one-bit compressed sensing (Boufounos and Baraniuk, 2008), linear
binary response model (Manski, 1975, 1985), and personalized medicine (Wang et
al., 2018). Interested readers could refer to Feng et al. (2019) for those
examples.

###  1.3 From Estimation to Inference

Incorporating high-dimensional data in the objective, i.e., moving forward
from (1.1) to (1.2), is not trivial, even for the purpose of estimation only.
Recently, Mukherjee et al. (2019) established the rate of convergence of the
(penalized) maximum score estimator for (1.2) in growing dimension, where the
dimension $d$ of $\bm{Z}$ is allowed to grow with $n$. In a related work, Feng
et al. (2019) proposed a regularized empirical risk minimization framework
with a smoothed surrogate loss for the estimation of the high-dimensional
parameter $\bm{\beta}^{*}$. They showed that the estimation problem is
nonregular, in the sense that there do not exist estimators of
$\bm{\beta}^{*}$ with root-$n$ convergence rate uniformly over a proper
parameter space.

Under (1.2), developing a valid statistical inference procedure is
challenging, even for fixed dimensional setting. Manski (1975, 1985)
considered the binary response model
$Y=\text{sign}(X-\bm{Z}^{T}\bm{\beta}+\epsilon)$, where $\epsilon$ may depend
on $(X,\bm{Z})$ but with $\text{Median}(\epsilon|X,\bm{Z})=0$. It can be shown
that the true coefficient $\bm{\beta}^{*}$ can be equivalently defined via
(1.2) with $w(-1)=w(1)=1/2$. The maximum score estimator is proposed to
estimate $\bm{\beta}$, and is later shown to have a non-Gaussian limiting
distribution (Kim et al., 1990). To tackle the challenge of nonstandard
limiting distribution of the maximum score estimator, Horowitz (1992) proposed
the smoothed maximum score estimator which is asymptotically normal in fixed
dimension.

On top of the nonregularity of the problem (1.2), the high-dimensionality of
the parameter adds an additional layer of complexity for inference. The reason
is that the estimator that minimizes the penalized loss function does not have
a tractable standard limiting distribution under high dimensionality, due to
the bias induced by the penalty term. For regular models (e.g., generalized
linear models), there is a growing literature on correcting the bias from the
penalty for valid inference, such as Javanmard and Montanari (2014); Zhang and
Zhang (2014); Van de Geer et al. (2014); Belloni et al. (2015); Ning et al.
(2017); Cai and Guo (2017); Fang et al. (2017); Neykov et al. (2018); Feng and
Ning (2019); Fang et al. (2020); Yi and Neykov (2021), among others. Their
main idea is to firstly construct a consistent estimator of the high
dimensional parameter via proper regularization, and then remove the bias (via
debiasing or decorrelation) in order to develop valid inferential statistics.
While these methods enjoy great success under regular models, it remains
unclear whether they can be applied to conduct valid inference in nonregular
models such as the problem we consider in this paper. To the best of our
knowledge, our work is the first one that provides valid inferential tools for
nonregular models in high dimension.

###  1.4 Our Contributions

In this paper, we propose a unified hypothesis testing framework for the one-
dimensional parameter $\theta^{*}$ as well as for the iMCID encoded as a
linear combination of $\bm{\beta}^{*}$. We start from considering the
hypothesis testing problem $H_{0}:\theta^{*}=0$ versus $H_{1}:\theta^{*}\neq
0$, where we treat $\bm{\gamma}^{*}$ as a high-dimensional nuisance parameter.
Built on the smoothed surrogate estimation framework (Feng et al., 2019), we
propose a bias corrected smoothed decorrelated score to form the score test
statistic.

There are several new ingredients in the construction of our score statistic.
First, the score function is derived based on a smoothed surrogate loss to
overcome the nonregularity due to the nonsmoothness of the 0-1 loss. Second,
unlike the existing works on high-dimensional inference, the score function
from the smoothed surrogate loss is biased. By explicitly estimating the bias
term, we derive the bias corrected score. Third, the decorrelation step,
developed by Ning et al. (2017) for regular models, is applied to reduce the
uncertainty of estimating high-dimensional nuisance parameters. Theoretically,
we show that under some conditions, the proposed score test statistic
converges in distribution to a standard Gaussian distribution under the null
hypothesis. We further establish the local asymptotic power of the test
statistic when $\theta^{*}$ deviates from $0$ in a local neighborhood. In
particular, we give the conditions under which the test statistic has
asymptotic power one.

When constructing the bias corrected smoothed decorrelated score, we need to
specify a bandwidth parameter, whose optimal choice depends on the unknown
smoothness of the data distribution. We further propose a double-smoothing
approach to select the optimal bandwidth parameter by minimizing the mean
squared error (MSE) of the score function. To our knowledge, such bandwidth
selection procedures have not been studied for high-dimensional models. We
show that under some extra smoothness assumptions, the ratio of the data-
driven bandwidth to the theoretically optimal bandwidth converges to one in
probability. From a practical perspective, with the selected bandwidth, the
proposed score statistic is fully data-driven and can be easily applied for
empirical research.

###  1.5 Paper Structure and Notation

The organization of this paper is as follows. In Section 2, we first provide
some background on the estimation of iMCID then introduce the bias corrected
smoothed deccorelated score and the associated test statistic. In Section 3,
we discuss the theoretical properties of the score test. The data-driven
bandwidth selection is addressed in Section 4. The corresponding results for
$\bm{c}_{0}^{T}\bm{\beta}^{*}$, the linear combination of parameters, are
briefly summarized in Section 5. Sections 6 and 7 contain simulation studies
and a real data example, respectively. All the technical details and proofs
are contained in the Appendix.

Throughout the paper, we adopt the following notation. For any set
${\mathcal{S}}$, we write $|{\mathcal{S}}|$ for its cardinality. For any
vector $\bm{v}\in\mathbb{R}^{d}$, we use $\bm{v}_{{\mathcal{S}}}$ to denote
the subvector of $\bm{v}$ with entries indexed by the set ${\mathcal{S}}$, and
define its $\ell_{q}$ norm as
$\|\bm{v}\|_{q}=(\sum_{j=1}^{d}|\bm{v}_{j}|^{q})^{1/q}$ for some real number
$q\geq 0$. For any matrix $\bm{M}\in\mathbb{R}^{d_{1}\times d_{2}}$, we denote
$||\bm{M}||_{\max}=\max_{i,j}|M_{ij}|$. For any two sequences $a_{n}$ and
$b_{n}$, we write $a_{n}\lesssim b_{n}$ if there exists some positive constant
$C$ such that $a_{n}\leq Cb_{n}$ for any $n$. We let $a_{n}\asymp b_{n}$ stand
for $a_{n}\lesssim b_{n}$ and $b_{n}\lesssim a_{n}$. Denote $a\vee
b=\max(a,b)$ and $a\wedge b=\min(a,b)$. For some function
$F(\theta,\bm{\gamma})$, we use $\nabla_{\theta}F(\theta,\bm{\gamma})$ and
$\nabla_{\bm{\gamma}}F(\theta,\bm{\gamma})$ to denote the partial derivatives
with respect to the corresponding parameters, and similarly
$\nabla^{2}_{\theta,\theta}F(\theta,\bm{\gamma})$ to denote the second order
derivative.

##  2 Methodology

###  2.1 Review of Penalized Smoothed Surrogate Estimation

Under the condition that $d\gg n$, the estimation of $\bm{\beta}^{*}$ via the
empirical risk minimization following (1.2) induces challenges from both
statistical and computational perspectives. The non-smoothness of $L_{01}(u)$
would cause the estimator to have a nonstandard convergence rate, which
happens even in the fixed low dimensional case (Kim et al., 1990). Moreover,
minimizing the empirical risk function based on the 0-1 loss is
computationally NP-hard and is often very difficult to implement. To tackle
these challenges, Feng et al. (2019) considered the following smoothed
surrogate risk

$\displaystyle
R_{\delta}(\bm{\beta})=\mathbb{E}\bigg{[}w(Y)L_{\delta,K}\big{\\{}Y(X-\bm{\beta}^{T}\bm{Z})\big{\\}}\bigg{]},$
(2.1)

where $L_{\delta,K}(u)=\int_{u/\delta}^{\infty}K(t)dt$ is a smoothed
approximation of $L_{01}(u)$, $K$ is a kernel function defined in Section 3
and $\delta>0$ is a bandwidth parameter. As the bandwidth $\delta$ shrinks to
0, $L_{\delta,K}(u)$ converges pointwisely to $L_{01}(u)$ (for any $u\neq 0$),
from which it can be shown that $\bm{\beta}^{*}$ also minimizes the smoothed
risk $R_{\delta}(\bm{\beta})$ up to a small approximation error. They further
proposed the following penalized smoothed surrogate estimator

$\displaystyle\widehat{\bm{\beta}}:=\mathop{\mathrm{argmin}}_{\bm{\beta}}R_{\delta}^{n}(\bm{\beta})+P_{\lambda}(\bm{\beta}),$
(2.2)

where $P_{\lambda}(\bm{\beta})$ is some sparsity inducing penalty (e.g.,
Lasso) with a tuning parameter $\lambda$, and $R_{\delta}^{n}(\bm{\beta})$ is
the corresponding empirical risk

$\displaystyle
R_{\delta}^{n}(\bm{\beta})=\frac{1}{n}\sum_{i=1}^{n}\bar{R}^{i}_{\delta}(\bm{\beta})=\frac{1}{n}\sum_{i=1}^{n}w(y_{i})L_{\delta,K}\Big{(}y_{i}(x_{i}-\bm{\beta}^{T}\bm{z}_{i})\Big{)}.$
(2.3)

Computationally, the empirical surrogate risk $R_{\delta}^{n}(\bm{\beta})$ is
a smooth function of $\bm{\beta}$, which renders the optimization more
tractable. Statistically, under some conditions, the estimator
$\widehat{\bm{\beta}}$ is shown to be rate-optimal, i.e., the convergence rate
of $\widehat{\bm{\beta}}$ matches the minimax lower bound up to a logarithmic
factor. We refer to Feng et al. (2019) for the detailed results.

###  2.2 Bias Corrected Smoothed Decorrelated Score

While Feng et al. (2019) showed that the penalized smoothed surrogate
estimator $\widehat{\bm{\beta}}$ is consistent, it does not automatically
equip with a practical inferential procedure for $\bm{\beta}^{*}$, or even for
the parameter $\theta^{*}$, mainly because of the sparsity inducing penalty.
In practice, how to draw valid statistical inference is often the ultimate
goal. In our motivating example, it is of critical importance to quantify the
uncertainty of $\bm{c}_{0}^{T}\bm{\beta}^{*}$ where $\bm{c}_{0}$ represents
the realized value of a new patient’s clinical profile. In other words, we
would like to develop a testing procedure for

$\displaystyle H_{0L}:\bm{c}_{0}^{T}\bm{\beta}^{*}=0\mbox{ versus
}H_{1L}:\bm{c}_{0}^{T}\bm{\beta}^{*}\neq 0.$ (2.4)

In this section, we focus on a special case of (2.4), the hypothesis test for
$\theta^{*}$,

$\displaystyle H_{0}:\theta^{*}=0\mbox{ versus }H_{1}:\theta^{*}\neq 0,$ (2.5)

where we treat $\bm{\gamma}$ as the nuisance parameter. Below we mainly
present the results for (2.5). Once they are clear, we can extend the results
to (2.4), which will be presented in Section 5.

In this section, we propose a new bias corrected smoothed decorrelated score
to quantify the uncertainty for testing (2.5). It is well known that the
classical score test is constructed based on the magnitude of the gradient of
the loglikelihood, or more generally, the empirical risk function associated
with $R(\bm{\beta})$ in (1.2). However, this construction breaks down in our
problem due to the following two reasons.

First, to construct the score statistic, one needs to plug in some estimate of
the nuisance parameter $\bm{\gamma}$ such as $\widehat{\bm{\gamma}}$ obtained
by partitioning
$\widehat{\bm{\beta}}=(\widehat{\theta},\widehat{\bm{\gamma}}^{T})^{T}$ in
(2.2). However, since $\bm{\gamma}$ is a high-dimensional parameter, the
estimation error from $\widehat{\bm{\gamma}}$ may become the leading term in
the asymptotic analysis of the score function. To deal with the high-
dimensional nuisance parameter, we use the decorrelated score, where the key
idea is to project the score of the parameter of interest to a high-
dimensional nuisance space (Ning et al., 2017). On the population level, it
takes the form

$\displaystyle\nabla_{\theta}R(\theta,\bm{\gamma})-\bm{\omega}^{*T}\nabla_{\bm{\gamma}}R(\theta,\bm{\gamma}),$
(2.6)

where the decorrelation vector is
$\bm{\omega}^{*}=\big{(}\nabla^{2}_{\bm{\gamma},\bm{\gamma}}R(\bm{\beta}^{*})\big{)}^{-1}\nabla^{2}_{\bm{\gamma},\theta}R(\bm{\beta}^{*})$.
The extra term $\bm{\omega}^{*T}\nabla_{\bm{\gamma}}R(\theta,\bm{\gamma})$ is
the key to reduce the error from estimating $\bm{\gamma}$.

Second, even if the above decorrelated score approach can successfully remove
the effect of the high-dimensional nuisance parameter, one cannot construct
the sample based decorrelated score from (2.6), as the sample version of
$R(\theta,\bm{\gamma})$ is non-differentiable, leading to the so called non-
standard inference. To circumvent this issue, we approximate
$R(\theta,\bm{\gamma})$ in (2.6) by the smoothed surrogate risk
$R_{\delta}(\theta,\bm{\gamma})$ in (2.1), that is

$\displaystyle\nabla_{\theta}R(\theta,\bm{\gamma})-\bm{\omega}^{*T}\nabla_{\bm{\gamma}}R(\theta,\bm{\gamma})=\big{\\{}\nabla_{\theta}R_{\delta}(\theta,\bm{\gamma})-\bm{\omega}^{*T}\nabla_{\bm{\gamma}}R_{\delta}(\theta,\bm{\gamma})\big{\\}}-\textrm{approximation
bias}.$ (2.7)

Since the empirical version of $R_{\delta}(\theta,\bm{\gamma})$ is smooth, we
define the (empirical) smoothed decorrelated score function as

$\displaystyle
S_{\delta}(\theta,\bm{\gamma})=\nabla_{\theta}R^{n}_{\delta}(\theta,\bm{\gamma})-\bm{\omega}^{*T}\nabla_{\bm{\gamma}}R^{n}_{\delta}(\theta,\bm{\gamma}).$

With $\bm{\gamma}$ estimated by $\widehat{\bm{\gamma}}$, the estimated score
function is then naturally defined as

$\displaystyle\widehat{S}_{\delta}(\theta,\widehat{\bm{\gamma}})=\nabla_{\theta}R^{n}_{\delta}(\theta,\widehat{\bm{\gamma}})-\widehat{\bm{\omega}}^{T}\nabla_{\bm{\gamma}}R^{n}_{\delta}(\theta,\widehat{\bm{\gamma}}),$
(2.8)

where $\widehat{\bm{\omega}}$ is an estimator of $\bm{\omega}^{*}$ defined
more precisely in Section 2.3.

In view of (2.7) and (2.8), the sample version of
$\nabla_{\theta}R_{\delta}(\theta,\bm{\gamma})-\bm{\omega}^{*T}\nabla_{\bm{\gamma}}R_{\delta}(\theta,\bm{\gamma})$
is given by $\widehat{S}_{\delta}(\theta,\widehat{\bm{\gamma}})$ and
therefore, to construct a valid score function, it remains to estimate the
approximation bias in (2.7). To proceed, we first analyze this term, which is
simply $\bm{v}^{*T}\nabla R_{\delta}(\bm{\beta}^{*})$ at
$\bm{\beta}=\bm{\beta}^{*}$, where $\bm{v}^{*}=(1,-\bm{\omega}^{*T})^{T}$.
After some analysis, we can show that the magnitude of the approximation bias
depends on the smoothness of $f(x|y,\bm{z})$, the conditional density of $X$
given $Y$ and $\bm{Z}$. To obtain an explicit form of the approximation bias,
we assume that $f(x|y,\bm{z})$ is $\ell$th order differentiable for some
$\ell\geq 2$, which is defined more precisely in Section 3. Under this
assumption, we can show that as the bandwidth parameter $\delta\rightarrow 0$,

$\displaystyle\bm{v}^{*T}\nabla
R_{\delta}(\bm{\beta}^{*})=\delta^{\ell}\mu^{*}(1+o(1)),$ (2.9)

where

$\displaystyle\mu^{*}:=\bm{v}^{*T}\bm{b}^{*}$ $\displaystyle=$
$\displaystyle\bm{v}^{*T}\Big{(}\int
K(u)\frac{u^{\ell}}{\ell!}du\Big{)}\sum_{y\in\\{-1,1\\}}w(y)\int
y\bm{z}f^{(\ell)}(\bm{\beta}^{*T}\bm{z}|y,\bm{z})f(y,\bm{z})d\bm{z},$ (2.10)
$\displaystyle=$ $\displaystyle\underbrace{\Big{(}\int
K(u)\frac{u^{\ell}}{\ell!}du\Big{)}}_{\gamma_{K,\ell}}\bm{v}^{*T}\underbrace{\mathbb{E}\Big{[}w(Y)Y\bm{Z}f^{(\ell)}(\bm{\beta}^{*T}\bm{Z}|Y,\bm{Z})\Big{]}}_{T^{(\ell)}(\bm{\beta}^{*})},$

and $f^{(\ell)}(x|y,\bm{z})$ denotes the $\ell$th order derivative of
$f(x|y,\bm{z})$ with respect to $x$.

To estimate the approximation bias $\bm{v}^{*T}\nabla
R_{\delta}(\bm{\beta}^{*})$, it suffices to estimate $\mu^{*}$. From (2.10),
we can see that, once $f^{(\ell)}(x|y,\bm{z})$ at $x=\bm{\beta}^{*T}\bm{z}$ is
estimated, we can construct a plug-in estimator for $\mu^{*}$. To be specific,
assume that a pilot kernel estimator with some kernel function $U$ and
bandwidth $h$ is available to estimate
$f^{(\ell)}(\bm{\beta}^{*T}\bm{z}|y,\bm{z})$. Then we cam estimate $\mu^{*}$
by

$\displaystyle\widehat{\mu}=\gamma_{K,\ell}\widehat{\bm{v}}^{T}\widehat{T}_{h,U}^{(\ell),n}(\widehat{\bm{\beta}}),$
(2.11)

where
$\widehat{T}_{h,U}^{(\ell),n}(\widehat{\bm{\beta}}):=\frac{1}{n}\sum_{i=1}^{n}w(y_{i})y_{i}\frac{\bm{z}_{i}}{h^{1+\ell}}U^{(\ell)}\big{(}\frac{\widehat{\bm{\beta}}^{T}\bm{z}_{i}-x_{i}}{h}\big{)}$
and $\widehat{\bm{v}}=(1,-\widehat{\bm{\omega}}^{T})^{T}$.

The last step to construct a valid score test is to find the asymptotic
variance of the smoothed decorrelated score $S_{\delta}(\bm{\beta}^{*})$.
Lemma 1 in the next section shows that the asymptotic variance of the
standardized decorrelated score $(n\delta)^{1/2}S_{\delta}(\bm{\beta}^{*})$ is
$\sigma^{*2}=\bm{v}^{*T}\bm{\Sigma}^{*}\bm{v}^{*}$, where

$\displaystyle\bm{\Sigma}^{*}$ $\displaystyle:=$
$\displaystyle\sum_{y\in\\{-1,1\\}}w(y)^{2}\int\bm{z}\bm{z}^{T}\int
K(u)^{2}duf(\bm{\beta}^{*T}\bm{z}|y,\bm{z})f(y,\bm{z})d\bm{z},$ (2.12)
$\displaystyle=$ $\displaystyle\underbrace{\Big{(}\int
K(u)^{2}du\Big{)}}_{\widetilde{\mu}_{K}}\underbrace{\mathbb{E}\Big{[}w(Y)^{2}\bm{Z}\bm{Z}^{T}f(\bm{\beta}^{*T}\bm{Z}|Y,\bm{Z})\Big{]}}_{H(\bm{\beta}^{*})},$

and thus $\sigma^{*}$ can be estimated by

$\displaystyle\widehat{\sigma}=\sqrt{\widetilde{\mu}_{K}\widehat{\bm{v}}^{T}\widehat{H}^{n}_{g,L}(\widehat{\bm{\beta}})\widehat{\bm{v}}},$
(2.13)

where
$\widehat{H}_{g,L}^{n}(\widehat{\bm{\beta}})=\frac{1}{n}\sum_{i=1}^{n}w^{2}(y_{i})\bm{z}_{i}\bm{z}_{i}^{T}\frac{1}{g}L(\frac{x_{i}-\widehat{\bm{\beta}}^{T}\bm{z}_{i}}{g})$
with some kernel function $L$ and bandwidth $g$.

Equipped with the smoothed decorrelated score
$\widehat{S}_{\delta}(\theta,\widehat{\bm{\gamma}})$ in (2.8), the estimate of
the approximation bias $\delta^{\ell}\widehat{\mu}$ in (2.11) and the estimate
of the asymptotic variance $\widehat{\sigma}^{2}$ in (2.13), we define the
bias corrected smoothed decorrelated score statistic as

$\displaystyle\widehat{U}_{n}=\sqrt{n\delta}\Big{(}\frac{\widehat{S}_{\delta}(0,\widehat{\bm{\gamma}})-\delta^{\ell}\widehat{\mu}}{\widehat{\sigma}}\Big{)}.$
(2.14)

######  Remark 1.

Compared to the existing decorrelated score approach (Ning et al., 2017), our
methodological innovation is to develop an explicit bias correction step to
remove the approximation bias in (2.7) induced by the smoothed surrogate risk.
From the theoretical aspect, our test statistic $\widehat{U}_{n}$ is rescaled
by $(n\delta)^{1/2}$ rather than the classical $n^{1/2}$ factor, which leads
to the non-standard rate of the decorrelated score not only under the null but
also under local alternatives; see Section 3.

###  2.3 Detailed Implementation

For numerical implementation, we follow the path-following algorithm presented
in Feng et al. (2019) to compute the initial estimator $\widehat{\bm{\beta}}$.
For the estimator $\widehat{\bm{\omega}}$, recall that $\bm{\omega}^{*}$
satisfies
$\nabla^{2}_{\bm{\gamma},\bm{\gamma}}R(\bm{\beta}^{*})\bm{\omega}^{*}=\nabla^{2}_{\bm{\gamma},\theta}R(\bm{\beta}^{*})$.
Since $\nabla^{2}R(\bm{\beta}^{*})$ can be approximated by the Hessian of the
smoothed surrogate loss $\nabla^{2}R^{n}_{\delta}(\bm{\beta}^{*})$, we
consider the following Dantzig type estimator $\widehat{\bm{\omega}}$, where

$\displaystyle\widehat{\bm{\omega}}=\mathop{\mathrm{argmin}}_{\bm{\omega}}||\bm{\omega}||_{1}~{}~{}~{}~{}~{}~{}~{}s.t.~{}~{}||\nabla^{2}_{\bm{\gamma},\theta}R_{\delta}(\widehat{\bm{\beta}})-\nabla^{2}_{\bm{\gamma},\bm{\gamma}}R_{\delta}(\widehat{\bm{\beta}})\bm{\omega}||_{\infty}\leq\lambda^{\prime},$
(2.15)

for some tuning parameter $\lambda^{\prime}>0$. Lemma 10 presented in the
Appendix shows the convergence rate of this estimator $\widehat{\bm{\omega}}$
as required in Assumption 6.

For implementing $\widehat{U}_{n}$, we note that the analysis of the
asymptotic distribution of $\widehat{U}_{n}$ is complicated by the dependence
between the estimator $\widehat{\bm{\beta}}$ and
$S_{\delta}(\theta,\bm{\gamma})$. To decouple the dependence and ease
theoretical development, we apply the cross-fitting technique to construct the
bias corrected smoothed decorrelated score. Specifically, instead of utilizing
the same set of samples for estimating
$\widehat{\bm{\beta}},\widehat{\bm{\omega}}$ and constructing the score
function $S_{\delta}(\theta,\bm{\gamma})$, we will firstly estimate
$\widehat{\bm{\beta}}$ using one set of samples, and then use the rest of
samples for estimating $\widehat{\bm{\omega}}$ and constructing
$S_{\delta}(\theta,\bm{\gamma})$. We can further switch the samples and
aggregate the decorrelated score. Without loss of generality, assume the
sample size $n$ is even and we divide the samples into two halves with equal
size for this purpose. Formally, denote
$\widehat{\bm{\beta}}^{(i)},\widehat{\bm{\omega}}^{(i)},i=1,2$ as the
estimator based on the $i$th fold of the samples, $\mathcal{N}_{i}$, and
similarly $\nabla
R_{\delta}^{n_{(i)}}(\bm{\beta}),\nabla^{2}R_{\delta}^{n_{(i)}}(\bm{\beta})$
as the corresponding gradient and Hessian. Define

$\widehat{S}^{(1)}_{\delta}(\theta,\widehat{\bm{\gamma}}^{(2)})=\nabla_{\theta}R^{n_{(1)}}_{\delta}(\theta,\widehat{\bm{\gamma}}^{(2)})-\widehat{\bm{\omega}}^{(1)T}\nabla_{\bm{\gamma}}R^{n_{(1)}}_{\delta}(\theta,\widehat{\bm{\gamma}}^{(2)}),$

and $\widehat{S}^{(2)}_{\delta}(\theta,\widehat{\bm{\gamma}}^{(1)})$ in a
similar way. The estimated decorrelated score via cross-fitting is

$\displaystyle\widehat{S}_{\delta}(\theta,\widehat{\bm{\gamma}})=\frac{1}{2}\big{(}\widehat{S}^{(1)}_{\delta}(\theta,\widehat{\bm{\gamma}}^{(2)})+\widehat{S}^{(2)}_{\delta}(\theta,\widehat{\bm{\gamma}}^{(1)})\big{)}.$
(2.16)

Similarly, we define the cross-fitted estimators $\widehat{\mu}$ and
$\widehat{\sigma}$ as

$\displaystyle\widehat{\mu}=$
$\displaystyle\frac{1}{2}\gamma_{K,\ell}(\widehat{\bm{v}}^{(1)T}\widehat{T}_{h,U}^{(\ell),n_{(1)}}(\widehat{\bm{\beta}}^{(2)})+\widehat{\bm{v}}^{(2)T}\widehat{T}^{(\ell),n_{(2)}}_{h,U}(\widehat{\bm{\beta}}^{(1)})),$
(2.17) $\displaystyle\widehat{\sigma}^{2}=$
$\displaystyle\frac{\widetilde{\mu}_{K}}{2}\bigg{[}\widehat{\bm{v}}^{(1)T}\widehat{H}^{n_{(1)}}_{g,K}(\widehat{\bm{\beta}}^{(2)})\widehat{\bm{v}}^{(1)}+\widehat{\bm{v}}^{(2)T}\widehat{H}^{n_{(2)}}_{g,K}(\widehat{\bm{\beta}}^{(1)})\widehat{\bm{v}}^{(2)}\bigg{]},$

where

$\displaystyle\widehat{T}_{h,U}^{(\ell),n_{(1)}}(\widehat{\bm{\beta}}^{(2)})=$
$\displaystyle\frac{1}{|\mathcal{N}_{1}|}\sum_{i\in\mathcal{N}_{1}}w(y_{i})y_{i}\frac{\bm{z}_{i}}{h^{1+\ell}}U^{(\ell)}\bigg{(}\frac{\widehat{\bm{\beta}}^{(2)T}\bm{z}_{i}-x_{i}}{h}\bigg{)},$
(2.18) $\displaystyle\widehat{H}_{g,L}^{n_{(1)}}(\widehat{\bm{\beta}}^{(2)})=$
$\displaystyle\frac{1}{|\mathcal{N}_{1}|}\sum_{i\in\mathcal{N}_{1}}w^{2}(y_{i})\bm{z}_{i}\bm{z}_{i}^{T}\frac{1}{g}L(\frac{x_{i}-\widehat{\bm{\beta}}^{(2)T}\bm{z}_{i}}{g}),$

and similarly for
$\widehat{T}_{h,U}^{(\ell),n_{(2)}}(\widehat{\bm{\beta}}^{(1)}),\widehat{H}_{g,L}^{n_{(2)}}(\widehat{\bm{\beta}}^{(1)})$.
Given $\widehat{S}_{\delta}(\theta,\widehat{\bm{\gamma}})$ in (2.16) and the
above estimators $\widehat{\mu}$ and $\widehat{\sigma}$, we can form the score
test statistic $\widehat{U}_{n}$ in the same way as in (2.14).

##  3 Theory

###  3.1 Assumptions

In this paper, we consider the following definition of function smoothness.

######  Definition 1.

We say the conditional density $f(x|y,\bm{z})$ of $X$ given $Y,\bm{Z}$ is
$\ell$th order smooth, if for any $\bm{z}$ and $y\in\\{-1,1\\}$, the
conditional density $f(x|y,\bm{z})$ is $\ell$-times continuously
differentiable in $x$ with derivatives $f^{(i)}(x|y,\bm{z})$ bounded by a
constant $C$, $|f^{(i)}(x|y,\bm{z})|\leq C$ for $i=1,\dotso,\ell$, and
$f^{(\ell)}(x|y,\bm{z})$ is Hölder continuous with some exponent $0<\zeta\leq
1$, that is, for any $\bm{z},\triangle$ and $y\in\\{-1,1\\}$,

$\displaystyle|f^{(\ell)}(x+\triangle|y,\bm{z})-f^{(\ell)}(x|y,\bm{z})|\leq
L\triangle^{\zeta},$ (3.1)

where $L>0$ is some constant.

######  Assumption 1.

We assume $f(x|y,\bm{z})$ is $\ell$th order smooth with some integer $\ell\geq
2$.

Assumption 1 concerns the smoothness of $f(x|y,\bm{z})$. To see why the
smoothness condition is important, notice that the gradient functions of (2.1)
and (1.2) are

$\displaystyle\nabla R_{\delta}(\bm{\beta})=$
$\displaystyle\sum_{y\in\\{-1,1\\}}w(y)\int
yz\Big{[}\int\frac{1}{\delta}K(\frac{y(x-\bm{\beta}^{T}\bm{z})}{\delta})f(x|y,\bm{z})dx\Big{]}f(y,\bm{z})d\bm{z}$
(3.2) $\displaystyle\nabla R(\bm{\beta})=$
$\displaystyle\sum_{y\in\\{-1,1\\}}w(y)\int
yzf(\bm{\beta}^{T}\bm{z}|y,\bm{z})f(y,\bm{z})d\bm{z},$

from which we can see that $f(\bm{\beta}^{T}\bm{z}|y,\bm{z})$ in $\nabla
R(\bm{\beta})$ is substituted by its kernel approximation
$\int\frac{1}{\delta}K(\frac{y(x-\bm{\beta}^{T}\bm{z})}{\delta})f(x|y,\bm{z})dx$,
and thus the difference between $\nabla R_{\delta}(\bm{\beta})$ and $\nabla
R(\bm{\beta})$ naturally depends on the smoothness of $f(x|y,\bm{z})$.

Notice that our smoothness condition in Definition 1 is slightly stronger than
the standard Hölder smoothness condition in the nonparametric literature
(Tsybakov, 2009). In particular, we require that $f^{(\ell)}(x|y,\bm{z})$ is
Hölder continuous with some exponent $0<\zeta\leq 1$ in (3.1). This additional
assumption is essential to show the rate of the bias estimator $\widehat{\mu}$
in (2.11).

Our next assumption is about the kernel function $K(t)$ that we first
introduced in the surrogate risk $R_{\delta}(\bm{\beta})$ in (2.1).

######  Assumption 2.

We assume $K(t)$ is a kernel function with support $[-1,1]$ that satisfies:
$K(t)=K(-t)$, $|K(t)|\leq K_{\max}<\infty\;\forall\;t\in\mathbb{R}$, $\int
K(t)dt=1$, $\int K^{2}(t)dt<\infty$, and $|K^{\prime}|<\infty$. We also assume
that $K$ degenerates at the boundaries. A kernel is said to be of order
$\ell\geq 1$ if it satisfies $\int
t^{j}K(t)dt=0,\;\forall\;j=1,\dotso,\ell-1$, $\int t^{\ell}K(t)dt\neq 0$, and
$\int|t|^{q}|K(t)|dt$ are bounded by a constant for any $q\in[\ell,\ell+1]$.

We now impose regularity conditions on $(X,Y,\bm{Z})$.

######  Assumption 3.

There exists a constant $c>0$ such that $c\leq\mathbb{P}(Y=1)\leq 1-c$ and the
weight function $w(\cdot)$ is positive and upper bounded by a constant.

######  Assumption 4.

We assume $\forall~{}j=1,\dotso,d,|Z_{j}|\leq M_{n}$ for some $M_{n}$ that
possibly depends on $n$, where $M_{n}^{2}\leq C\sqrt{n\delta/\log(d)}$ for
some constant $C>0$. We also assume that $\mathbb{E}[Z_{j}^{4}|Y=y]$ is
bounded by a constant for $y\in\\{1,-1\\}$.

######  Assumption 5.

We assume $\sigma^{*}=\sqrt{\bm{v}^{*T}\bm{\Sigma}^{*}\bm{v}^{*}}$ and
$|\mu^{*}|=|\bm{v}^{*T}\bm{b}^{*}|$ are bounded away from 0 and infinity by
some constants.

Assumption 4 requires the boundedness of $\bm{Z}$ and the fourth order moment.
Notice that if each component of $\bm{Z}$ is sub-Gaussian with bounded sub-
Gaussian norm, Assumption 4 is satisfied with high probability providing
$(\log d)^{3}/(n\delta)=\mathcal{O}(1)$, which is a mild assumption.
Assumption 5 ensures that the asymptotic variance of the smoothed decorrelated
score $\sigma^{*}$ and the approximation bias $\mu^{*}$ are nonzero and
bounded.

We impose the following assumption on the estimators of $\bm{\beta}^{*}$ and
$\bm{\omega}^{*}$.

######  Assumption 6.

Assume there are estimators $\widehat{\bm{\beta}}$ and
$\widehat{\bm{v}}=(1,-\widehat{\bm{\omega}}^{T})^{T}$ with

$||\widehat{\bm{\beta}}-\bm{\beta}^{*}||_{1}\lesssim\eta_{1}(n)~{}~{}~{}~{}\text{and}~{}~{}~{}~{}||\widehat{\bm{v}}-\bm{v}^{*}||_{1}/||\bm{v}^{*}||_{1}\lesssim\eta_{2}(n),$

for some non-random sequences $\eta_{1}(n),\eta_{2}(n)$ converging to $0$ as
$n\rightarrow\infty$.

It is shown by Feng et al. (2019) that, under some conditions, the estimator
$\widehat{\bm{\beta}}$ in (2.2) achieves the (near) minimax-optimal rate
$\eta_{1}(n)=\sqrt{s}(\frac{s\log(d)}{n})^{\ell/(2\ell+1)}$ by choosing
$\lambda\asymp\sqrt{\log d/(n\delta)}$ and $\delta\asymp(s\log
d/n)^{1/(2\ell+1)}$, where $s=\|\bm{\beta}^{*}\|_{0}$. For
$\widehat{\bm{v}}=(1,-\widehat{\bm{\omega}}^{T})^{T}$, we assume
$||\widehat{\bm{v}}-\bm{v}^{*}||_{1}\lesssim||\bm{v}^{*}||_{1}\eta_{2}(n)$.
Notice that the term $||\bm{v}^{*}||_{1}$ is not absorbed into $\eta_{2}(n)$
only for notational simplicity. In Lemma 10 in the Appendix, we show that a
Dantzig type estimator $\widehat{\bm{v}}$ attains a fast rate $\eta_{2}(n)$.

###  3.2 Theoretical Results

We start from the following lemma which characterizes the asymptotic
distribution of the decorrelated score function evaluated at the true
parameter $\bm{\beta}^{*}$.

######  Lemma 1.

Under Assumptions 1 \- 5, if
$(||\bm{v}^{*}||_{1}M_{n})^{3}/(n\delta)^{1/2}=o(1)$ and $\delta=o(1)$, then

$\displaystyle\sqrt{n\delta}\frac{\bm{v}^{*T}(\nabla
R_{\delta}^{n}(\bm{\beta}^{*})-\nabla
R_{\delta}(\bm{\beta}^{*}))}{\sqrt{\bm{v}^{*T}\bm{\Sigma}^{*}\bm{v}^{*}}}\stackrel{{\scriptstyle
d}}{{\rightarrow}}N(0,1),$ (3.3)

where

$\displaystyle\bm{v}^{*T}\nabla
R_{\delta}(\bm{\beta}^{*})=\delta^{\ell}\bm{v}^{*T}\bm{b}^{*}(1+o(1)).$ (3.4)

This lemma illustrates the asymptotic bias and variance of the smoothed
decorrelated score. The former is essentially the approximation bias term in
(2.7). Since $\mu^{*}=\bm{v}^{*T}\bm{b}^{*}$ and
$\sigma^{*}=\sqrt{\bm{v}^{*T}\bm{\Sigma}^{*}\bm{v}^{*}}$ are bounded by
constants, the asymptotic bias and standard deviation of $\bm{v}^{*T}\nabla
R_{\delta}^{n}(\bm{\beta}^{*})$ are of order $\delta^{\ell}$ and
$(n\delta)^{-1/2}$, respectively. Thus, choosing $\delta=cn^{-1/(2\ell+1)}$
for any constant $c>0$ attains the optimal bias and variance trade-off in
terms of the rate. Indeed, we show in Section 4 that the bandwidth $\delta$
that minimizes an estimate of the MSE of the decorrelated score not only
attains the optimal order $n^{-1/(2\ell+1)}$ but also finds the optimal
constant $c$. Note that in this lemma we require
$(||\bm{v}^{*}||_{1}M_{n})^{3}/(n\delta)^{1/2}=o(1)$ to verify the Lindeberg
condition in the central limit theorem, which holds as long as $\delta$ does
not shrink to 0 too fast.

Our first main theorem characterizes the asymptotic normality of the
decorrelated score under the null hypothesis with nuisance parameters
$\bm{\gamma}^{*}$ and $\bm{\omega}^{*}$ estimated by those in Assumption 6.

######  Theorem 1.

Under Assumptions 1 \- 6, if
$(||\bm{v}^{*}||_{1}M_{n})^{3}/(n\delta)^{1/2}=o(1)$,
$\frac{\log(d)}{n\delta^{3}}=o(1)$, $n\delta^{2\ell+1}=O(1)$, and

$\displaystyle(n\delta)^{1/2}||\bm{v}^{*}||_{1}\bigg{(}\frac{\eta_{1}(n)}{\delta}\vee\eta_{2}(n)\bigg{)}\bigg{(}\sqrt{\frac{\log(d)}{n\delta}}\vee\delta^{\ell}\vee
M_{n}^{2}\eta_{1}(n)\bigg{)}=o(1),$ (3.5)

then under $H_{0}:\theta^{*}=0$, it holds that

$\displaystyle\frac{\sqrt{n\delta}\widehat{S}_{\delta}(0,\widehat{\bm{\gamma}})-\sqrt{n\delta^{2\ell+1}}\mu^{*}}{\sigma^{*}}\stackrel{{\scriptstyle
d}}{{\rightarrow}}N(0,1).$ (3.6)

Theorem 1 implies that the decorrelated score with some high-dimensional plug-
in estimators $\widehat{\bm{\gamma}}$ and $\widehat{\bm{\omega}}$ has the same
asymptotic distribution as in Lemma 1. Several conditions are needed to show
this result. The first condition
$(||\bm{v}^{*}||_{1}M_{n})^{3}/(n\delta)^{1/2}=o(1)$ is from Lemma 1, and the
second condition $\frac{\log(d)}{n\delta^{3}}=o(1)$ is also mild as long as
$\delta$ does not go to zero too fast. The third condition
$n\delta^{2\ell+1}=O(1)$ guarantees that the higher order bias of the
decorrelated score can be ignored and therefore it suffices to only correct
for the leading bias term in (3.4). In particular, if the rate optimal
bandwidth $\delta=cn^{-1/(2\ell+1)}$ is chosen, from (3.6), we can see that
the asymptotic bias of
$\sqrt{n\delta}\widehat{S}_{\delta}(0,\widehat{\bm{\gamma}})$ is a constant
$\sqrt{n\delta^{2\ell+1}}\mu^{*}=c^{\ell+1/2}\mu^{*}$, which is further
removed to attain the mean zero Gaussian limiting distribution.

We now elaborate the condition (3.5). Roughly speaking, the term
$\sqrt{\frac{\log(d)}{n\delta}}\vee\delta^{\ell}\vee M_{n}^{2}\eta_{1}(n)$
comes from the bound for $\|\nabla
R^{n_{(1)}}_{\delta}(\theta,\widehat{\bm{\gamma}}^{(2)})-\nabla
R(\theta,{\bm{\gamma}})\|_{\infty}$, in which the three terms correspond to
the stochastic and approximation error of $\nabla
R^{n_{(1)}}_{\delta}(\theta,\widehat{\bm{\gamma}}^{(2)})$ and the plug-in
error from $\widehat{\bm{\gamma}}^{(2)}$, respectively. Indeed, the cross-
fitting technique guarantees the independence between
$\widehat{\bm{\gamma}}^{(2)}$ and $\nabla
R^{n_{(1)}}_{\delta}(\theta,{\bm{\gamma}})$, which is the key to obtain a
sharp bound for the stochastic and approximation error. Thus, condition (3.5)
simply means that these error terms interacting with the estimation error of
$\widehat{\bm{\gamma}}$ and $\widehat{\bm{\omega}}$ are sufficiently small. We
can further simplify the condition (3.5) by plugging the order of
$\eta_{1}(n)$ and $\eta_{2}(n)$. Recall that Feng et al. (2019) derived
$\eta_{1}(n)=\sqrt{s}(\frac{s\log(d)}{n})^{\ell/(2\ell+1)}$, where
$s=\|\bm{\beta}^{*}\|_{0}$. The Lemma 10 in the Appendix implies that
$\eta_{2}(n)=s^{\prime}(\log(d)/n)^{(\ell-1)/(2\ell+1)}$, where
$s^{\prime}=\|\bm{\omega}^{*}\|_{0}$. Assuming
$||\bm{v}^{*}||_{1},M_{n}=\mathcal{O}(1)$ for now, condition (3.5) becomes

$\displaystyle s^{(4\ell+1)/(4\ell+2)}(s^{(4\ell+1)/(4\ell+2)}\vee
s^{\prime})n^{-(\ell-1)/(2\ell+1)}(\log d)^{(4\ell-1)/(4\ell+2)}=o(1)$ (3.7)

when taking $\delta\asymp n^{-1/(2\ell+1)}$. If we consider the extreme case
with $\ell\rightarrow\infty$, (LABEL:eq_simple_rate) reduces to $s(s\vee
s^{\prime})\log d=o(n^{1/2})$.

Recall that in our score statistic $\widehat{U}_{n}$ in (2.14), we plug in the
estimators $\widehat{\mu}$ and $\widehat{\sigma}$ for $\mu^{*}$ and
$\sigma^{*}$. In Lemmas 11 and 12 in the Appendix, we establish the rate of
convergence of $\widehat{\mu}$ and $\widehat{\sigma}$. Under the assumption
that $|\widehat{\mu}-\mu^{*}|=o_{p}(1)$ and
$|\widehat{\sigma}-\sigma^{*}|=o_{p}(1)$, the Slutsky’s theorem implies that
the bias corrected decorrelated score statistic
$\widehat{U}_{n}\stackrel{{\scriptstyle d}}{{\rightarrow}}N(0,1)$ under the
null hypothesis.

Accordingly, given the desired significance level $\alpha$, we define the test
function as

$T_{DS}=I(|\widehat{U}_{n}|>\Phi^{-1}(1-\alpha/2)),$

where $\Phi^{-1}(\cdot)$ is the inverse function of the cdf of the standard
normal distribution. Thus, our result shows that the Type I error of the test
$T_{DS}$ converges to $\alpha$ asymptotically, i.e.,
$\mathbb{P}(T_{DS}=1|H_{0})\rightarrow\alpha$.

Now denote

$\nabla^{2}_{\theta|\bm{\gamma}}R(\bm{\beta}^{*})=\nabla^{2}_{\theta\theta}R(\bm{\beta}^{*})-\nabla^{2}_{\theta\bm{\gamma}}R(\bm{\beta}^{*})(\nabla^{2}_{\bm{\gamma}\bm{\gamma}}R(\bm{\beta}^{*}))^{-1}\nabla^{2}_{\bm{\gamma}\bm{\theta}}R(\bm{\beta}^{*}).$

Our second main theorem characterizes the limiting behavior of
$\widehat{U}_{n}$ under the local alternative hypothesis
$H_{1}:\theta^{*}=\widetilde{C}n^{-\phi}$ for some constants
$\widetilde{C}\neq 0$ and $\phi>0$.

######  Theorem 2.

Under the conditions in Theorem 1 and Lemmas 11 and 12 in the Appendix, we
further assume

$\displaystyle||\bm{v}^{*}||_{1}^{2}M_{n}^{4}n^{1-4\phi}/\delta=o(1),\quad(n\delta)^{1/2}||\bm{v}^{*}||_{1}(\eta_{1}(n)\vee\eta_{2}(n))M_{n}n^{-\phi}=o(1),$
(3.8)

and $\widehat{\mu},\widehat{\sigma}$ are consistent estimators of
$\mu^{*},\sigma^{*}$. Then the following results hold under the local
alternative hypothesis $H_{1}:\theta^{*}=\widetilde{C}n^{-\phi}$.

  * •

If
$\widetilde{C}(n\delta)^{1/2}n^{-\phi}\nabla^{2}_{\theta|\bm{\gamma}}R(\bm{\beta}^{*})\sigma^{*-1}\rightarrow\xi$
for some constant $\xi$, then it holds that

$\displaystyle\widehat{U}_{n}\stackrel{{\scriptstyle
d}}{{\rightarrow}}N(-\xi,1).$ (3.9)

  * •

If
$\widetilde{C}(n\delta)^{1/2}n^{-\phi}\nabla^{2}_{\theta|\bm{\gamma}}R(\bm{\beta}^{*})\sigma^{*-1}\rightarrow\infty$,
then for any fixed $t$, it holds that

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}(|\widehat{U}_{n}|>t)=1.$
(3.10)

In addition to the conditions imposed in Theorem 1 and Lemmas 11 and 12 in the
Appendix, we further require two additional conditions involving the magnitude
of $\theta^{*}$ in (LABEL:eq_power_condition). The first condition
$||\bm{v}^{*}||_{1}^{2}M_{n}^{4}n^{1-4\phi}/\delta=o(1)$ is imposed to ensure
the local asymptotic normality (LAN) in terms of the parameter $\theta^{*}$.
Assuming $\|\bm{v}^{*}\|_{1}$ and $M_{n}=O(1)$, this condition requires
$\phi>\frac{\ell+1}{2(2\ell+1)}$ with the rate optimal bandwidth $\delta\asymp
n^{-1/(2\ell+1)}$. The second condition in (LABEL:eq_power_condition) is
similar to (3.5), that controls the magnitude of $||\nabla
R(0,\bm{\gamma}^{*})||_{\infty}$ and
$||\nabla^{2}R_{\delta}(\theta^{*},\bm{\gamma}^{*})-\nabla^{2}R_{\delta}(0,\bm{\gamma}^{*})||_{\max}$
under the alternative hypothesis.

The conclusions in Theorem 2 are shown in two settings (for now assume
$\nabla^{2}_{\theta|\bm{\gamma}}R(\bm{\beta}^{*})$ is a constant and
$\delta\asymp n^{-1/(2\ell+1)}$). Firstly, if
$\widetilde{C}(n\delta)^{1/2}n^{-\phi}\nabla^{2}_{\theta|\bm{\gamma}}R(\bm{\beta}^{*})\sigma^{*-1}\rightarrow\xi$,
which holds when $\phi=\ell/(2\ell+1)$, then it implies that the asymptotic
power of the test $T_{DS}$ at level $\alpha$ is

$\mathbb{P}(T_{DS}=1|H_{1})\rightarrow
A_{\alpha}:=1-\Phi(\xi+\Phi^{-1}(1-\alpha/2))+\Phi(\xi-\Phi^{-1}(1-\alpha/2)).$

Notice that $A_{\alpha}$ is greater than the level $\alpha$ as long as
$\xi\neq 0$. This indicates the asymptotic unbiasedness of the test $T_{DS}$.
Second, the power of the test tends to 1 when
$(n\delta)^{1/2}n^{-\phi}\nabla^{2}_{\theta|\bm{\gamma}}R(\bm{\beta}^{*})\sigma^{*-1}\rightarrow\infty$,
which holds when $\phi<\ell/(2\ell+1)$. This implies that the proposed test
can successfully detect the nonzero $\theta^{*}$ whose magnitude exceeds the
order of $n^{-\ell/(2\ell+1)}$. In contrast, for regular models, the local
alternative that is detectable is of the standard parametric rate $n^{-1/2}$.
This is a major theoretical difference between our work and the existing ones
on debiased and decorrelated methods.

##  4 Data-Driven Bandwidth Selection

In the previous section, we establish the theoretical property of the bias
corrected smoothed decorrelated score when the underlying conditional density
$f(x|y,\bm{z})$ is $\ell$th order smooth. However, this smoothness parameter
$\ell$ is typically unknown in practice, leading to the following two
complications. First, in Assumption 2, a kernel function $K$ of the same order
is applied, which implicitly requires the knowledge on the smoothness
parameter $\ell$. In practice, the choice of kernel functions is often
determined by the user’s preference rather than the theory. Since high order
kernels may exacerbate the problem of variability, choosing low order kernels
of 2 or 4 is often recommended (even if the density is more smooth); see
Härdle et al. (1992). Second, the optimal bandwidth $\delta\asymp
n^{-1/(2\ell+1)}$ that balances the asymptotic bias and variance of the
decorrelated score in Lemma 1 also depends on the unknown $\ell$. It is well
known from the nonparametric literature that the choice of bandwidth is an
extremely important problem of both theoretical and practical values
(Silverman, 1986; Bowman, 1984; Sheather and Jones, 1991; Hall et al., 1992;
Jones et al., 1996). In this section, we focus on how to choose the bandwidth
$\delta$ in a data-driven manner. In view of the above discussion on the
kernels, we assume that a low order kernel $K$ is chosen (for simplicity, we
still denote its order by $\ell$) and meanwhile the underlying conditional
density has a higher order smoothness parameter.

######  Assumption 7.

We assume that the kernel $K$ is of order $\ell$ and $f(x|y,\bm{z})$ is
$(\ell+r)$th order smooth for some $\ell\geq 2$ and $r>0$.

We define the optimal bandwidth $\delta^{*}$ as the one that minimizes the MSE
of the smoothed decorrelated score:

$\displaystyle\delta^{*}=\mathop{\mathrm{argmin}}_{\delta}M(\delta),~{}~{}\textrm{where}~{}~{}M(\delta)=\mathbb{E}[(\bm{v}^{*T}\nabla
R_{\delta}^{n}(\bm{\beta}^{*}))^{2}].$ (4.1)

A direct bias-variance decomposition of $M(\delta)$ gives

$\displaystyle
M(\delta)=\frac{1}{n}\mathbb{E}[(\bm{v}^{*T}\nabla\bar{R}_{\delta}^{1}(\bm{\beta}^{*}))^{2}]+\frac{n-1}{n}(\bm{v}^{*T}\nabla
R_{\delta}(\bm{\beta}^{*}))^{2}:=\frac{1}{n}V(\delta)+\frac{n-1}{n}SB(\delta),$
(4.2)

where $\nabla\bar{R}_{\delta}^{1}(\bm{\beta}^{*})$ is defined in (2.3). To
simplify the estimation, we define
$V(\delta)=\mathbb{E}[(\bm{v}^{*T}\nabla\bar{R}_{\delta}^{1}(\bm{\beta}^{*}))^{2}]$
and use it as a proxy for the variance. We use $SB(\delta)=(\bm{v}^{*T}\nabla
R_{\delta}(\bm{\beta}^{*}))^{2}$ to denote the squared bias.

To estimate $\delta^{*}$, the main idea is to construct estimators
$\widehat{V}(\delta)$ and $\widehat{SB}(\delta)$ for $V(\delta)$ and
$SB(\delta)$ and then estimate $\delta^{*}$ by

$\widehat{\delta}=\mathop{\mathrm{argmin}}_{\delta}\widehat{M}(\delta)~{}~{}\textrm{where}~{}~{}\widehat{M}(\delta)=\frac{1}{n}\widehat{V}(\delta)+\frac{n-1}{n}\widehat{SB}(\delta).$

######  Remark 2.

Before we proceed to develop the estimators of $V(\delta)$ and $SB(\delta)$,
we first clarify the difference between the variance and squared bias terms
$V(\delta)$ and $SB(\delta)$ and the asymptotic variance and bias
$\sigma^{*2}$ and $\delta^{\ell}\mu^{*}$ in Lemma 1. From the proof of Lemma
1, we can show that ${SB}(\delta)=(\delta^{\ell}\mu^{*})^{2}(1+o(1))$ and
${V}(\delta)=\delta^{-1}\sigma^{*2}(1+o(1))$ as $\delta\rightarrow 0$, and
thus $\sigma^{*2}/\delta$ and $(\delta^{\ell}\mu^{*})^{2}$ are the asymptotic
versions of $V(\delta)$ and $SB(\delta)$, respectively. As a result, one may
attempt to estimate the optimal bandwidth by minimizing the asymptotic MSE
$\widehat{\sigma}^{2}/\delta+(\delta^{\ell}\widehat{\mu})^{2}$ with the plug-
in estimators $\widehat{\sigma}$ and $\widehat{\mu}$ developed in the previous
section. However, the asymptotic MSE depends on the unknown smoothness $\ell$
and therefore is not appropriate for bandwidth selection in practice. Thus, we
have to use a different strategy to estimate $V(\delta)$ and $SB(\delta)$.

To estimate $V(\delta)$, we consider the following moment estimator

$\displaystyle\widehat{V}(\delta)=\frac{1}{2}\big{(}\widehat{\bm{v}}^{(1)T}\widehat{\bm{\Gamma}}^{(1)}(\delta)\widehat{\bm{v}}^{(1)}+\widehat{\bm{v}}^{(2)T}\widehat{\bm{\Gamma}}^{(2)}(\delta)\widehat{\bm{v}}^{(2)}\big{)},$
(4.3)

where
$\widehat{\bm{\Gamma}}^{(1)}(\delta)=\frac{1}{|\mathcal{N}_{1}|}\sum_{i\in\mathcal{N}_{1}}\nabla\bar{R}_{\delta}^{i}(\widehat{\bm{\beta}}^{(2)})\nabla\bar{R}_{\delta}^{i}(\widehat{\bm{\beta}}^{(2)})^{T}$
and similarly for $\widehat{\bm{\Gamma}}^{(2)}$.

Next, we focus on the problem of estimating the squared bias $SB(\delta)$.
After some algebra, the bias term $B(\delta)$ can be written as

$\displaystyle B(\delta)=\bm{v}^{*T}\underbrace{(\nabla
R_{\delta}(\bm{\beta}^{*})-\nabla
R(\bm{\beta}^{*}))}_{A(\bm{\beta}^{*},\delta)}=\int_{u}K(u)\biggl{\\{}\bm{v}^{*T}(\nabla
R(u\delta,\bm{\beta}^{*})-\nabla R(\bm{\beta}^{*}))\biggr{\\}}du,$ (4.4)

where with a bit abuse of notation, we use $\nabla
R(u\delta,\bm{\beta}^{*})=\sum_{y}w(y)\int_{\bm{z}}\bm{z}yf(u\delta+\bm{\beta}^{*T}\bm{z}|y,\bm{z})f(y,\bm{z})d\bm{z}$
to denote the population gradient with a bias induced by $u\delta$. Consider
the following estimator with cross-fitting

$\displaystyle\widehat{B}(\delta)=\frac{1}{2}\bigg{(}\widehat{\bm{v}}^{(1)T}\frac{1}{|\mathcal{N}_{1}|}\sum_{i\in\mathcal{N}_{1}}A_{i}(\widehat{\bm{\beta}}^{(2)},\delta)+\widehat{\bm{v}}^{(2)T}\frac{1}{|\mathcal{N}_{2}|}\sum_{i\in\mathcal{N}_{2}}A_{i}(\widehat{\bm{\beta}}^{(1)},\delta)\bigg{)},$
(4.5)

where

$\displaystyle A_{i}(\widehat{\bm{\beta}},\delta)=$
$\displaystyle\int_{u}K(u)w(y_{i})\frac{\bm{z}_{i}y_{i}}{b}\Big{[}J(\frac{x_{i}-\widehat{\bm{\beta}}^{T}\bm{z}_{i}-u\delta}{b})-J(\frac{x_{i}-\widehat{\bm{\beta}}^{T}\bm{z}_{i}}{b})\Big{]}du,$
(4.6)

and $J$ is another pilot kernel function of order $r$ with bandwidth $b$.
Essentially, we substitute $\nabla R(u\delta,\bm{\beta}^{*})$ and $\nabla
R(\bm{\beta}^{*})$ in (4.4) with their corresponding kernel estimators. This
approach is motivated by the “double-smoothing” technique in nonparametric
statistics (Härdle et al., 1992; Neumann et al., 1995), and is also related to
the “smoothed cross validation” approach (Hall et al., 1992). In these works,
when the target function has higher order smoothness than the kernel function
applied for estimation, a second kernel smoothing procedure is applied for
estimating the bias. Using a similar strategy, in our case, we apply the
second smoother with kernel $J$ in the estimators of $\nabla
R(u\delta,\bm{\beta}^{*})$ and $\nabla R(\bm{\beta}^{*})$. We now estimate the
squared bias by $\widehat{SB}(\delta)=\widehat{B}(\delta)^{2}$.

To analyze theoretical properties of the estimates, let’s define
$\Delta=[q_{1}n^{-1+\epsilon_{1}},q_{2}n^{-1+\epsilon_{2}}]$ as the range of
bandwidth $\delta$ for some constants $0<q_{1}\leq q_{2}$ and
$0<\epsilon_{1}<\epsilon_{2}<1$. Since the optimal bandwidth $\delta^{*}$ is
of order $n^{-1/(2\ell+1)}$, we can guarantee $\delta^{*}\in\Delta$ for some
suitable $\epsilon_{1}$ and $\epsilon_{2}$. Under some conditions, the uniform
convergence rates of $\widehat{V}(\delta)$ and $\widehat{SB}(\delta)$ are
given by

$|\widehat{V}(\delta)-V(\delta)|\lesssim\frac{\psi_{1}(n,\delta)}{\delta},~{}~{}~{}~{}|\widehat{SB}(\delta)-SB(\delta)|\lesssim\delta^{2\ell}\psi_{2}(n,\delta),$

uniformly over all $\delta\in\Delta$, where

$\psi_{1}(n,\delta)=||\bm{v}^{*}||_{1}^{2}\bigg{(}\eta_{2}(n)\vee\sqrt{\frac{\log(n\vee
d)}{n\delta}}\vee M_{n}\eta_{1}(n)\bigg{)},$

and

$\psi_{2}(n,\delta)=||\bm{v}^{*}||_{1}^{2}\bigg{(}\sqrt{\frac{\log(n\vee
d)}{nb^{2\ell+1}}}\vee(\delta\vee b)^{r}\vee
M_{n}\eta_{1}(n)(1\vee\frac{M_{n}\eta_{1}(n)}{\delta^{\ell}})\vee\eta_{2}(n)\bigg{)}.$

We refer to Lemmas 2 and 3 in Appendix A for the formal statement of the
results and further interpretations of the rates. Now we are ready to present
the following theorem.

######  Theorem 3.

Under the conditions in Lemmas 2 and 3 in Appendix A, if
$\psi_{1}(n,\delta),\psi_{2}(n,\delta)=o(1)$, then with probability tending to
1,

$\frac{\widehat{\delta}-\delta^{*}}{\delta^{*}}\lesssim\psi_{1}(n,\delta^{*})\vee\psi_{2}(n,\delta^{*}).$

Thus, $\widehat{\delta}/\delta^{*}\rightarrow 1$ in probability.

Notice that in Theorem 3, the MSE-optimal bandwidth satisfies
$\delta^{*}\asymp n^{-1/(2\ell+1)}$. To examine the order of
$\psi_{1}(n,\delta^{*})\vee\psi_{2}(n,\delta^{*})$, we consider $d\gtrsim n$
(high-dimensional case) and optimize the bandwidth $b$ in
$\psi_{2}(n,\delta^{*})$, leading to $b\asymp(\log d/n)^{1/(2\ell+2r+1)}$.
Furthermore, take $\eta_{1}(n),\eta_{2}(n)$ as the rate from Feng et al.
(2019) and Lemma 10, and assume that
$||\bm{v}^{*}||_{1},M_{n}=\mathcal{O}(1)$, then
$\psi_{1}(n,\delta^{*})\vee\psi_{2}(n,\delta^{*})$ can be simplified to

$\displaystyle\Big{(}\frac{\log d}{n}\Big{)}^{r/(2\ell+2r+1)}\vee
s^{(4\ell+1)/(2\ell+1)}\Big{(}\frac{\log^{2}d}{n}\Big{)}^{\ell/(2\ell+1)}\vee
s^{\prime}\Big{(}\frac{\log d}{n}\Big{)}^{(\ell-1)/(2\ell+1)}.$ (4.7)

There are in general two scenarios: firstly, if $r$ is small relative to
$\ell$ while $s,s^{\prime}$ do not grow too fast, then $(\log
d/n)^{r/(2\ell+2r+1)}$ from the extra smoothing step for estimating
$SB(\delta)$ will be the dominant term. However, if $r$ is relatively large
such that $r/(2\ell+2r+1)\geq(\ell-1)/(2\ell+1)$, then the last two terms in
(4.7) due to $\widehat{\bm{\beta}}$ and $\widehat{\bm{\omega}}$ dominate. In
this case, even if $f(x|y,\bm{z})$ has a large amount of extra smoothness $r$,
the convergence rate of $\widehat{\delta}$ cannot be further improved, as the
rate is dominated by the error from $\widehat{\bm{\beta}}$ and
$\widehat{\bm{\omega}}$.

##  5 Hypothesis Testing for the Linear Combination

Our results presented thus far are for the hypothesis testing (2.5), where
$\theta$ is simply a single element in $\bm{\beta}$. In applications,
researchers may also be interested in the hypothesis testing (2.4), a linear
combination of the parameter $\bm{\beta}$. For example, as we mentioned at the
beginning of Section 2.2, in the study of inferring iMCID, it is of interest
to test $\bm{c}_{0}^{T}\bm{\beta}^{*}=0$ where $\bm{c}_{0}$ represents the
realized value of a new patient’s clinical profile and we assume $c_{01}\neq
0$ without loss of generality. Below, we show that all the methods and results
developed for (2.5) can be extended to (2.4) in a parallel manner.

To test the hypothesis (2.4), consider the one to one reparametrization
$(\theta,\bm{\gamma})\rightarrow(\xi,\bm{\gamma})$, where
$\xi=\bm{c}_{0}^{T}\bm{\beta}$. Under this new set of parameters, the null
hypothesis can be written as $H_{0L}:\xi^{*}=0$, and the smoothed surrogate
loss reduces to
$R^{n}_{\delta}(\frac{\xi-\bm{c}_{02}^{T}\bm{\gamma}}{c_{01}},\bm{\gamma})$,
where we write $\bm{c}_{0}=(c_{01},\bm{c}_{02}^{T})^{T}$ with
$\bm{c}_{02}\in\mathbb{R}^{d-1}$. Define
$\bm{C}=\begin{bmatrix}\frac{1}{c_{01}}&0\\\\[3.00003pt]
\frac{-\bm{c}_{02}}{c_{01}}&\bm{I}_{d-1}\end{bmatrix}\in\mathbb{R}^{d\times
d}$. From the chain rule, we can show that

$\displaystyle\nabla_{(\xi,\bm{\gamma})}R^{n}_{\delta}(\frac{\xi-\bm{c}_{02}^{T}\bm{\gamma}}{c_{01}},\bm{\gamma})=\bm{C}\nabla
R_{\delta}^{n}(\theta,\bm{\gamma}),~{}~{}\nabla^{2}_{(\xi,\bm{\gamma}),(\xi,\bm{\gamma})}R^{n}_{\delta}(\frac{\xi-\bm{c}_{02}^{T}\bm{\gamma}}{c_{01}},\bm{\gamma})=\bm{C}\nabla^{2}R_{\delta}^{n}(\theta,\bm{\gamma})\bm{C}^{T},$
(5.1)

and similarly for
$R_{\delta}(\frac{\xi-\bm{c}_{02}^{T}\bm{\gamma}}{c_{01}},\bm{\gamma})$ and
$R(\frac{\xi-\bm{c}_{02}^{T}\bm{\gamma}}{c_{01}},\bm{\gamma})$. Therefore,
following the same idea as in Section 2.2, we define the smoothed decorrelated
score as

$\displaystyle
S^{L}_{\delta}(\xi,\bm{\gamma},\bm{\omega}^{*}_{L})=\nabla_{\xi}R^{n}_{\delta}(\frac{\xi-\bm{c}_{02}^{T}\bm{\gamma}}{c_{01}},\bm{\gamma})-\bm{\omega}^{*T}_{L}\nabla_{\bm{\gamma}}R^{n}_{\delta}(\frac{\xi-\bm{c}_{02}^{T}\bm{\gamma}}{c_{01}},\bm{\gamma}),$
(5.2)

where
$\bm{\omega}^{*}_{L}=\Big{[}\nabla^{2}_{\bm{\gamma},\bm{\gamma}}R(\frac{\xi-\bm{c}_{02}^{T}\bm{\gamma}}{c_{01}},\bm{\gamma})\Big{]}^{-1}\nabla^{2}_{\bm{\gamma},\xi}R(\frac{\xi-\bm{c}_{02}^{T}\bm{\gamma}}{c_{01}},\bm{\gamma}).$

Write $\bm{v}^{*}_{L}=(1,\bm{\omega}^{*T}_{L})^{T}$ and denote
$\mu^{*}_{L},\sigma^{*}_{L}$ as the (scaled) asymptotic bias and standard
deviation of the score function
$S^{L}_{\delta}(\xi,\bm{\gamma},\bm{\omega}^{*}_{L})$, i.e.,

$\displaystyle\mu^{*}_{L}=\bm{v}^{*T}_{L}\bm{C}\bm{b}^{*},~{}~{}~{}~{}~{}~{}\sigma^{*}_{L}=\sqrt{\bm{v}^{*T}_{L}\bm{C}\bm{\Sigma}^{*}\bm{C}^{T}\bm{v}^{*}_{L}},$
(5.3)

where $\bm{b}^{*}$ and $\bm{\Sigma}^{*}$ are defined in (2.10) and (2.12),
respectively. From above we can see that the estimation methods for
$\bm{\omega}^{*},\mu^{*}$ and $\sigma^{*}$ proposed in Section 2.2 can be
easily extended to obtain corresponding estimators for
$\bm{\omega}^{*}_{L},\mu^{*}_{L}$ and $\sigma^{*}_{L}$. Given these
estimators, we define the test statistics for $H_{0L}$ as

$\widehat{U}^{L}_{n}=(n\delta)^{1/2}\frac{S^{L}_{\delta}(0,\widehat{\bm{\gamma}},\widehat{\bm{\omega}}_{L})-\delta^{\ell}\widehat{\mu}_{L}}{\widehat{\sigma}_{L}}.$

Accordingly, all the parallel results presented in Section 3 and Section 4 can
be developed. In the interest of space, we only present the following result
that characterizes the asymptotic distribution of $\widehat{U}_{n}^{L}$ under
the null hypothesis $H_{0L}$. All other parallel results are omitted.

######  Theorem 4.

If Assumptions 1 \- 6 hold with
$\mu^{*},\sigma^{*},\bm{v}^{*},\widehat{\bm{v}}$ substituted by
$\mu^{*}_{L},\sigma^{*}_{L},\bm{v}^{*}_{L},\widehat{\bm{v}}_{L}$, and in
addition $\widehat{\mu}_{L}$ and $\widehat{\sigma}_{L}$ are consistent
estimators of $\mu^{*}_{L}$ and $\sigma^{*}_{L}$, respectively, then under the
same conditions as in Theorem 1 and the null hypothesis $H_{0L}:\xi^{*}=0$, it
holds that

$\widehat{U}_{n}^{L}\stackrel{{\scriptstyle d}}{{\rightarrow}}N(0,1).$

##  6 Simulation Studies

In this section, we evaluate the empirical performance of the proposed
methods. Although many models can be formulated as special cases of our
problem (1.2), here we mainly consider the following binary response model

$\displaystyle Y=\text{sign}(X-\bm{\beta}^{*T}\bm{Z}+\epsilon),$ (6.1)

where $\epsilon$ possibly depends on $X$ and $\bm{Z}$ but the median of
$\epsilon$ given $X$ and $\bm{Z}$ is 0.

###  6.1 Experiments with pre-specified bandwidth

In the first set of experiments, we evaluate the performance of the proposed
test statistic with pre-specified bandwidth. We use Gaussian kernel $K$ of
order 2 with bandwidth pre-specified at $\delta=n^{-1/5}$. The choices of
other tuning parameters are detailed in Appendix D. Throughout this
subsection, we consider sample size $n=800$, dimension $d=100,500,1000$ and
generate $\beta^{*}_{2},\dotso,\beta^{*}_{s}$ by sampling from a uniform
distribution within $[1,2]$ for $s=3,10$. The first coordinate $\beta_{1}^{*}$
would vary depending on the purpose of the experiment, and the rest
coordinates of $\bm{\beta}^{*}$ are all set to 0. After that, the coefficient
vector is then normalized such that $||\bm{\beta}^{*}||_{2}=1$. We generate
$X\sim N(0,1)$ and $\bm{Z}\sim N(0,\bm{\Sigma}_{\rho})$, where
$(\bm{\Sigma}_{\rho})_{jk}=\rho^{|j-k|}$ with $\rho=0.2,0.5,0.7$. For all
cases, the simulations are repeated 250 times.

In the first scenario, we let $\epsilon\sim
N(0,0.2^{2}(1+2(X-\bm{\beta}^{*T}\bm{Z})^{2}))$, which is referred to as
Heteroskedastic Gaussian scenario later on. We compare the proposed smoothed
decorrelated score test (SDS) with the decorrelated score test method (DS)
(Ning et al., 2017) and Honest confidence region method (Honest) (Belloni et
al., 2016) from the “hdm” package. We fix the significance level at $0.05$ and
firstly evaluate the performance of the tests under the null hypothesis
$H_{0}:\beta^{*}_{1}=0$. In this case, we set $\beta^{*}_{1}=0$. Note that the
R code for the DS and Honest approaches is tailored for the high-dimensional
logistic regression, which differs from the above data generating process.

Table 6.1 reports the empirical Type I error rate under the first scenario.
The error rate from the SDS method is generally close to the nominal
significance level 0.05, which empirically verifies the theoretical results in
Theorem 1. For both Honest and DS methods, the empirical Type I error rate
seems to be consistently higher or lower than the nominal level. This is
expected as these two methods only work for the logistic regression. By taking
a closer look at the Normal Q-Q plot of the test statistics, we observe that
the distribution of the test statistics from the Honest and DS methods deviate
substantially from Gaussian, as opposed to those yield by the proposed SDS
method. Please see Appendix E.2 for more details.

Table 6.1: The empirical Type I error rate of the tests under the Heteroskedastic Gaussian scenario from SDS, DS and Honest methods. |  |  | $s=3$ |  |  | $s=10$ |
---|---|---|---|---|---|---|---
d | method | $\rho=0.2$ | $\rho=0.5$ | $\rho=0.7$ | $\rho=0.2$ | $\rho=0.5$ | $\rho=0.7$
100 | SDS | 5.6% | 5.0% | 6.4% | 4.8% | 4.8% | 5.2%
| DS | 1.2% | 2.0% | 2.0% | 2.0% | 1.8% | 1.8%
| Honest | 5.2% | 5.6% | 7.6% | 5.4% | 5.2% | 6.8%
500 | SDS | 4.8% | 4.4% | 5.6% | 5.6% | 5.0% | 4.8%
| DS | 0.2% | 0.4% | 0.4% | 0.2% | 0.0% | 0.4%
| Honest | 7.0% | 10.8% | 7.6% | 8.2% | 6.8% | 7.2%
1000 | SDS | 4.4% | 6.0% | 5.6% | 5.0% | 5.4% | 5.0%
| DS | 0.0% | 0.4% | 0.2% | 0.0% | 0.0% | 0.4%
| Honest | 10.0% | 10.4% | 12.6% | 12.4% | 6.4% | 15.2%
Table 6.2: The empirical Type I error rate of the tests under the Heteroskedastic Uniform scenario from SDS, DS and Honest methods. |  |  | $s=3$ |  |  | $s=10$ |
---|---|---|---|---|---|---|---
d | method | $\rho=0.2$ | $\rho=0.5$ | $\rho=0.7$ | $\rho=0.2$ | $\rho=0.5$ | $\rho=0.7$
100 | SDS | 6.0% | 5.6% | 6.8% | 7.2% | 6.8% | 7.2%
| DS | 2.0% | 0.4% | 0.8% | 1.6% | 0.8% | 1.2%
| Honest | 9.6% | 17.6% | 19.2% | 9.6% | 18.4% | 21.6%
500 | SDS | 6.4% | 6.8% | 6.4% | 6.0% | 7.6% | 7.2%
| DS | 1.6% | 0.8% | 0.0% | 0.8% | 1.2% | 0.4%
| Honest | 11.2% | 15.6% | 20.0% | 13.6% | 17.2% | 22.8%
1000 | SDS | 8.4% | 7.6% | 8.8% | 9.2% | 7.6% | 8.0%
| DS | 0.0% | 0.4% | 0.4% | 1.2% | 0.8% | 0.0%
| Honest | 13.6% | 15.2% | 22.4% | 16.0% | 19.2% | 20.8%

In the second scenario, we let $\epsilon\sim
0.2\cdot\text{Unif}(-G(X,\bm{Z}),G(X,\bm{Z}))$, where
$G(x,\bm{z})=\sqrt{1+2(x-\bm{\beta}^{*T}\bm{z})^{2}}$. In other words, the
error $\epsilon$ follows a uniform distribution such that its range depends on
the covariates $X,\bm{Z}$ (we will call it Heteroskedastic Uniform scenario).
Similar to the Heteroskedastic Gaussian case, we compare SDS method with DS
and Honest methods and study the empirical Type I error rate. From Table 6.2
we can see that the proposed method yields Type I error close to the nominal
level as opposed to the other two. The Normal QQ-plots in Figure E.1 in
Appendix E.2 further confirm the asymptotic normality of our SDS test
statistics. The above results suggest that in practice, if the underlying data
generating process is the binary response model, our proposed approach
provides valid inferential results while the existing approaches fail.

Next, we investigate the empirical power of the SDS method. We use the same
data generating processes as in the above two scenarios, but instead of
setting $\beta^{*}_{1}=0$, we vary $\beta^{*}_{1}$ in the grid
$\\{0.02,0.05,0.075,0.10,0.15,0.20,0.25,0.30\\}$ for the Heteroskedastic
Gaussian case, and $\\{0.025,0.05,0.075,0.10,0.125,0.15,0.175\\}$ for the
Heteroskedastic Uniform case. Similarly, we consider $s=3,10$,
$d=100,500,1000$ and $\rho=0.2,0.5,0.7$. Figure 6.1 shows the empirical
rejection rate of the SDS method when $s=10$ (see Appendix E.3 for the results
when $s=3$). Note that we do not compare with the DS and Honest methods for
the empirical power, because these two tests do not maintain the desired Type
I error in our scenarios. We can see that for all considered cases, the
empirical power converges to 1 as the magnitude of the signal $\beta^{*}_{1}$
becomes larger, which agrees with Theorem 2. In addition, we find that the
dimension $d$ has minor effects on the empirical power, which is reasonable as
Theorem 2 only depends on $\log d$ via the condition
(LABEL:eq_power_condition). Finally, we note that the power of the test
deteriorates as the correlation of the design increases.

Figure 6.1: Empirical rejection rate of the proposed test under both scenarios
with $s=10$, $d=100,500,1000$ and $\rho=0.2,0.5,0.7$.

###  6.2 Experiments with data-driven bandwidth

In the next set of experiments, we study the empirical performance of the
data-driven bandwidth selection approach. We firstly study the type I error
and power of our SDS method for testing $H_{0}:\beta^{*}_{1}=0$ versus
$H_{1}:\beta^{*}_{1}\neq 0$ with data-driven bandwidth. We consider the same
data generating processes as in Section 6.1 with $n=800,d=100,s=3,10$ and
$\rho=0.2,0.5,0.7$. We seek for the minimizer of the estimated MSE over
$\delta\in[0.1,1.2]$ and each experiment is repeated 250 times. After
$\widehat{\delta}$ is obtained, we plug-in it into the test statistic and
estimate the bias and variance as discussed in Appendix D. With the same
implementations, we also evaluate the empirical power of the test by varying
$\beta^{*}_{1}$ in the same grid as in Section 6.1.

Table 6.3: The empirical Type I error rate of the tests under the Heteroskedastic Gaussian and Uniform scenarios with data-driven bandwidth $\widehat{\delta}$. |  | $s=3$ |  |  | $s=10$ |
---|---|---|---|---|---|---
Data generating process | $\rho=0.2$ | $\rho=0.5$ | $\rho=0.7$ | $\rho=0.2$ | $\rho=0.5$ | $\rho=0.7$
Heteroskedastic Gaussian | 6.8% | 7.2% | 5.6% | 8.4% | 7.2% | 6.4%
Heteroskedastic Uniform | 8.8% | 8.0% | 7.6% | 8.4% | 6.8% | 5.2%

Figure 6.2: Empirical power of the tests under the Heteroskedastic Gaussian
and Uniform scenarios with data-driven bandwidth $\widehat{\delta}$.

Table 6.3 shows the empirical Type I error rate over 250 repetitions when
$\beta^{*}_{1}=0$, and Figure 6.2 shows the empirical power for different
$\beta^{*}_{1}\neq 0$ in these scenarios. Similar to the case when the
bandwidth $\delta$ is pre-specified, the empirical Type I errors are generally
close to the nominal level 0.05, and the empirical power converges to 1 as
$\beta^{*}_{1}$ becomes larger.

In the Appendix E.4, we show that the data-driven bandwidth is indeed close to
the theoretically optimal bandwidth that minimizes the MSE. From the above
results, we recommend using the data-driven bandwidth selection approach in
practice.

##  7 Analysis of ChAMP Trial

In this section we analyze the ChAMP (Chondral Lesions And Meniscus
Procedures) trial (Bisson et al., 2017), which contains clinical information
about $n=138$ patients undergoing arthroscopic partial meniscectomy (APM), a
knee surgery for meniscal tears. The response variable is $Y=1$ if the patient
is healthy/satisfactory and $-1$ otherwise, obtained from the SF-36 survey.
The continuous measurement $X$ encodes the WOMAC pain score change from the
baseline to one-year after the surgery. The dataset also contains $d=160$
additional variables from the patient’s clinical profile, denoted by $\bm{Z}$.
The scientific question is to determine the iMCID, defined as a linear
combination of the covariates $\bm{\beta}^{T}\bm{Z}$, such that the treatment
of debriding chondral lesions can be claimed as clinically significant by
comparing the WOMAC pain score change with this individualized threshold. As
we can see, this application naturally fits into our formulation (1.2) with
weight function $w(y)=1/\mathbb{P}(Y=y)$. The goal of the analysis is to
address this question by providing valid inferential results for each
component of $\bm{\beta}$.

We apply the proposed SDS test for $H_{0j}:\beta_{j}=0$ versus
$H_{1j}:\beta_{j}\neq 0$, where $1\leq j\leq d$. We use the same tuning
parameter setting for estimating $\bm{\beta}^{*},\omega^{*}$ and the
asymptotic bias and variance of the score function following Appendix D. For
comparison, we also apply the DS and Honest methods discussed in Section 6.1.

Table 7.1: The three significant variables (with p-value $<0.05/d$) identified by the proposed SDS method, and their corresponding p-values obtained from the DS and Honest methods. p-value | KQOL_6wk | flex_inj_pre | KSymp_3mo
---|---|---|---
SDS | 3.583e-07 | 5.575e-05 | 4.482e-05
DS | 0.0168 | 0.0340 | 0.0213
Honest | 0.0591 | 0.0241 | 0.4383

Table 7.1 lists the three significant variables (i.e., those with p-value
smaller than $0.05/d=$3.125e-04) from the proposed SDS approach.
Interestingly, all of them are clinically relevant and can provide meaningful
implications for iMCID. The significance of the variable KQOL_6wk, which
represents the KOOS score for quality of life at 6-week, definitely indicates
how the patients recover at a relatively early stage after the surgery. The
variable flex_inj_pre means the degree of flexion right before the surgery.
Its significance recommends that the baseline disease severity would affect
the magnitude of iMCID—this similar phenomenon was also discovered in the
clinical literature for other types of diseases, such as the shoulder pain
reduction study (Heald et al., 1997). The third variable KSymp_3mo is the KOOS
score for other symptoms at 3-month. In some previous analysis of ChAMP trial
where only estimate is available but without inference results, this variable
has the second largest coefficient (Feng et al., 2019).

The results from the DS and the Honest methods are different. Firstly, the DS
method only yields 1 significant variable and the Honest method yields 13.
From Table 7.1, the three significant variables identified by the proposed SDS
method cannot be identified by either DS or Honest. In general, compared to
the proposed SDS method, DS identifies fewer significant variables while
Honest identifies more. This phenomenon is also evident from the simulation
results in Section 6.1. On the other hand, the results from the three methods
do not completely contradict with each other. For instance, the two
significant variables, KQOL_6wk and KSymp_3mo, identified from the proposed
SDS method, has the fourth and fifth smallest p-values in the DS method.
Please refer to Appendix E.5 for more results of this analysis.

In general, recall that DS and Honest methods are devised for the logistic
regression, while our proposed SDS method can produce valid inference results
under the binary response model (6.1), which is more flexible since the
distribution of $\epsilon$ is left unspecified. Therefore, we expect that the
significant variables identified by SDS are potentially more reliable and
clinically more relevant. Our results presented in this session echo this
rationale.

## Acknowledgment

Ning is supported in part by National Science Foundation (NSF) CAREER award
DMS-1941945 and NSF award DMS-1854637. Zhao is supported in part by NSF award
DMS-2122074 and a start-up grant from the University of Wisconsin-Madison.

## References

  * Abdullah et al. (2020) Abdullah, S. S., Rostamzadeh, N., Sedig, K., Garg, A. X. and McArthur, E. (2020).  Visual analytics for dimension reduction and cluster analysis of high dimensional electronic health records.  In Informatics, vol. 7. Multidisciplinary Digital Publishing Institute.
  * Angst et al. (2017) Angst, F., Aeschlimann, A. and Angst, J. (2017).  The minimal clinically important difference raised the significance of outcome effects above the statistical level, with methodological implications for future studies.  Journal of Clinical Epidemiology 82 128–136.
  * Bellamy et al. (2001) Bellamy, N., Carr, A., Dougados, M., Shea, B. and Wells, G. (2001).  Towards a definition of “difference” in osteoarthritis.  The Journal of Rheumatology 28 427–430.
  * Belloni et al. (2015) Belloni, A., Chernozhukov, V. and Kato, K. (2015).  Uniform post-selection inference for least absolute deviation regression and other z-estimation problems.  Biometrika 102 77–94.
  * Belloni et al. (2016) Belloni, A., Chernozhukov, V. and Wei, Y. (2016).  Post-selection inference for generalized linear models with many controls.  Journal of Business & Economic Statistics 34 606–619.
  * Bisson et al. (2015) Bisson, L. J., Kluczynski, M. A., Wind, W. M., Fineberg, M. S., Bernas, G. A., Rauh, M. A., Marzo, J. M. and Smolinski, R. J. (2015).  Design of a randomized controlled trial to compare debridement to observation of chondral lesions encountered during partial meniscectomy: The champ (chondral lesions and meniscus procedures) trial.  Contemporary clinical trials 45 281–286.
  * Bisson et al. (2017) Bisson, L. J., Kluczynski, M. A., Wind, W. M., Fineberg, M. S., Bernas, G. A., Rauh, M. A., Marzo, J. M., Zhou, Z. and Zhao, J. (2017).  Patient outcomes after observation versus debridement of unstable chondral lesions during partial meniscectomy: the chondral lesions and meniscus procedures (champ) randomized controlled trial.  The Journal of Bone and Joint Surgery 99 1078\.
  * Boufounos and Baraniuk (2008) Boufounos, P. T. and Baraniuk, R. G. (2008).  1-bit compressive sensing.  In 2008 42nd Annual Conference on Information Sciences and Systems. IEEE.
  * Bowman (1984) Bowman, A. W. (1984).  An alternative method of cross-validation for the smoothing of density estimates.  Biometrika 71 353–360.
  * Cai and Guo (2017) Cai, T. T. and Guo, Z. (2017).  Confidence intervals for high-dimensional linear regression: Minimax rates and adaptivity.  The Annals of Statistics 45 615–646.
  * Chernozhukov et al. (2018) Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W. and Robins, J. (2018).  Double/debiased machine learning for treatment and structural parameters.  The Econometrics Journal 21 C1–C68.
  * Erdogan et al. (2016) Erdogan, B. D., Leung, Y. Y., Pohl, C., Tennant, A. and Conaghan, P. G. (2016).  Minimal clinically important difference as applied in rheumatology: an omeract rasch working group systematic review and critique.  The Journal of Rheumatology 43 194–202.
  * Fang et al. (2020) Fang, E. X., Ning, Y. and Li, R. (2020).  Test of significance for high-dimensional longitudinal data.  Annals of Statistics 48 2622–2645.
  * Fang et al. (2017) Fang, E. X., Ning, Y. and Liu, H. (2017).  Testing and confidence intervals for high dimensional proportional hazards models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 79 1415–1437.
  * Feng and Ning (2019) Feng, H. and Ning, Y. (2019).  High-dimensional mixed graphical model with ordinal data: Parameter estimation and statistical inference.  In The 22nd International Conference on Artificial Intelligence and Statistics.
  * Feng et al. (2019) Feng, H., Ning, Y. and Zhao, J. (2019).  Nonregular and minimax estimation of individualized thresholds in high dimension with binary responses.  arXiv preprint arXiv:1905.10888 .
  * Hall et al. (1992) Hall, P., Marron, J. and Park, B. U. (1992).  Smoothed cross-validation.  Probability Theory and Related Fields 92 1–20.
  * Härdle et al. (1992) Härdle, W., Hall, P. and Marron, J. (1992).  Regression smoothing parameters that are not far from their optimum.  Journal of the American Statistical Association 87 227–233.
  * Heald et al. (1997) Heald, S. L., Riddle, D. L. and Lamb, R. L. (1997).  The shoulder pain and disability index: the construct validity and responsiveness of a region-specific disability measure.  Physical Therapy 77 1079–1089.
  * Horowitz (1992) Horowitz, J. L. (1992).  A smoothed maximum score estimator for the binary response model.  Econometrica: journal of the Econometric Society 505–531.
  * Jaeschke et al. (1989) Jaeschke, R., Singer, J. and Guyatt, G. H. (1989).  Measurement of health status: ascertaining the minimal clinically important difference.  Controlled Clinical Trials 10 407–415.
  * Javanmard and Montanari (2014) Javanmard, A. and Montanari, A. (2014).  Confidence intervals and hypothesis testing for high-dimensional regression.  The Journal of Machine Learning Research 15 2869–2909.
  * Jayadevappa et al. (2017) Jayadevappa, R., Cook, R. and Chhatre, S. (2017).  Minimal important difference to infer changes in health-related quality of life—a systematic review.  Journal of Clinical Epidemiology 89 188–198.
  * Jones et al. (1996) Jones, M. C., Marron, J. S. and Sheather, S. J. (1996).  A brief survey of bandwidth selection for density estimation.  Journal of the American Statistical Association 91 401–407.
  * Kim et al. (1990) Kim, J., Pollard, D. et al. (1990).  Cube root asymptotics.  The Annals of Statistics 18 191–219.
  * Lassere et al. (2001) Lassere, M., van der HEIJDE, D. and Johnson, K. R. (2001).  Foundations of the minimal clinically important difference for imaging.  The Journal of Rheumatology 28 890–891.
  * Manski (1975) Manski, C. F. (1975).  Maximum score estimation of the stochastic utility model of choice.  Journal of econometrics 3 205–228.
  * Manski (1985) Manski, C. F. (1985).  Semiparametric analysis of discrete response: Asymptotic properties of the maximum score estimator.  Journal of econometrics 27 313–333.
  * McGlothlin and Lewis (2014) McGlothlin, A. E. and Lewis, R. J. (2014).  Minimal clinically important difference: defining what really matters to patients.  Journal of the American Medical Association 312 1342–1343.
  * Mukherjee et al. (2019) Mukherjee, D., Banerjee, M. and Ritov, Y. (2019).  Non-standard asymptotics in high dimensions: Manski’s maximum score estimator revisited.  arXiv preprint arXiv:1903.10063 .
  * Neumann et al. (1995) Neumann, M. H. et al. (1995).  Automatic bandwidth choice and confidence intervals in nonparametric regression.  The Annals of Statistics 23 1937–1959.
  * Neykov et al. (2018) Neykov, M., Ning, Y., Liu, J. S., Liu, H. et al. (2018).  A unified theory of confidence regions and testing for high-dimensional estimating equations.  Statistical Science 33 427–443.
  * Ning et al. (2017) Ning, Y., Liu, H. et al. (2017).  A general theory of hypothesis tests and confidence regions for sparse high dimensional models.  The Annals of Statistics 45 158–195.
  * Norman et al. (2003) Norman, G. R., Sloan, J. A. and Wyrwich, K. W. (2003).  Interpretation of changes in health-related quality of life: the remarkable universality of half a standard deviation.  Medical Care 582–592.
  * Samsa et al. (1999) Samsa, G., Edelman, D., Rothman, M. L., Williams, G. R., Lipscomb, J. and Matchar, D. (1999).  Determining clinically important differences in health status measures.  Pharmacoeconomics 15 141–155.
  * Sheather and Jones (1991) Sheather, S. J. and Jones, M. C. (1991).  A reliable data-based bandwidth selection method for kernel density estimation.  Journal of the Royal Statistical Society: Series B (Methodological) 53 683–690.
  * Silverman (1986) Silverman, B. W. (1986).  Density Estimation for Statistics and Data Analysis, vol. 26.  CRC press.
  * Tsybakov (2009) Tsybakov, A. B. (2009).  Introduction to Nonparametric Estimation.  Springer.
  * Van de Geer et al. (2014) Van de Geer, S., Bühlmann, P., Ritov, Y., Dezeure, R. et al. (2014).  On asymptotically optimal confidence regions and tests for high-dimensional models.  The Annals of Statistics 42 1166–1202.
  * Wang et al. (2018) Wang, L., Zhou, Y., Song, R. and Sherwood, B. (2018).  Quantile-optimal treatment regimes.  Journal of the American Statistical Association 113 1243–1254.
  * Wells et al. (2001) Wells, G., Beaton, D., Shea, B., Boers, M., Simon, L., Strand, V., Brooks, P. and Tugwell, P. (2001).  Minimal clinically important differences: review of methods.  The Journal of Rheumatology 28 406–412.
  * Wyrwich et al. (1999a) Wyrwich, K. W., Nienaber, N. A., Tierney, W. M. and Wolinsky, F. D. (1999a).  Linking clinical relevance and statistical significance in evaluating intra-individual changes in health-related quality of life.  Medical Care 469–478.
  * Wyrwich et al. (1999b) Wyrwich, K. W., Tierney, W. M. and Wolinsky, F. D. (1999b).  Further evidence supporting an sem-based criterion for identifying meaningful intra-individual changes in health-related quality of life.  Journal of Clinical Epidemiology 52 861–873.
  * Xu et al. (2014) Xu, T., Wang, J. and Fang, Y. (2014).  A model-free estimation for the covariate-adjusted youden index and its associated cut-point.  Statistics in medicine 33 4963–4974.
  * Yi and Neykov (2021) Yi, Y. and Neykov, M. (2021).  A new perspective on debiasing linear regressions.  arXiv preprint arXiv:2104.03464 .
  * Zhang and Zhang (2014) Zhang, C.-H. and Zhang, S. S. (2014).  Confidence intervals for low dimensional parameters in high dimensional linear models.  Journal of the Royal Statistical Society: Series B: Statistical Methodology 217–242.

Supplementary Appendix

##  Appendix A Discussions and Lemmas for the bandwidth selection

######  Lemma 2.

Under Assumptions 1 \- 6, if $\sqrt{\frac{\log(n\vee d)}{n\delta}}\vee
M_{n}\eta_{1}(n)=o(1)$, and $M_{n}^{2}\sqrt{\log(n\vee
d)/n^{\epsilon_{1}}}=\mathcal{O}(1)$, it holds that

$\displaystyle|\widehat{V}(\delta)-V(\delta)|\lesssim\frac{\psi_{1}(n,\delta)}{\delta},$
(A.1)

uniformly over all $\delta\in\Delta$, where

$\displaystyle\psi_{1}(n,\delta)=||\bm{v}^{*}||_{1}^{2}\bigg{(}\eta_{2}(n)\vee\sqrt{\frac{\log(n\vee
d)}{n\delta}}\vee M_{n}\eta_{1}(n)\bigg{)}.$ (A.2)

###### Proof.

See Section B.4 for a detailed proof. ∎

######  Lemma 3.

Suppose Assumptions 1\- 7 hold. Choose $J$ as a proper kernel function of
order $r$ satisfying the conditions of $K$ in Assumption 2. In addition,
assume $J$ is $\ell$ times continuously differentiable and $J^{(i)}$
degenerates at the boundary for $i=1,\dotso,\ell-1$. If $M_{n}\leq
C\sqrt{nb/\log(nd)}$ for some constant $C$ and $\frac{\log(n\vee
d)}{nb^{2\ell+1}}=o(1)$, then it holds that

$\displaystyle|\widehat{SB}(\delta)-SB(\delta)|\lesssim\delta^{2\ell}\psi_{2}(n,\delta)$
(A.3)

uniformly over all $\delta\in\Delta$, where

$\displaystyle\psi_{2}(n,\delta)=||\bm{v}^{*}||_{1}^{2}\bigg{(}\sqrt{\frac{\log(n\vee
d)}{nb^{2\ell+1}}}\vee(\delta\vee b)^{r}\vee
M_{n}\eta_{1}(n)(1\vee\frac{M_{n}\eta_{1}(n)}{\delta^{\ell}})\vee\eta_{2}(n)\bigg{)}.$
(A.4)

###### Proof.

See Section B.5 for a detailed proof. ∎

In view of (A.3) and (A.4), the error terms $\sqrt{\frac{\log(n\vee
d)}{nb^{2\ell+1}}}$ and $(\delta\vee b)^{r}$ correspond to the variance and
bias when estimating $A(\bm{\beta}^{*},\delta)$ in (4.4) with the extra
smoothing step. This rate is similar to the case of estimating the $\ell$th
derivative of a density function using a kernel of order $r$ in the context of
kernel density estimation under the smoothness Assumption 7. The terms
$M_{n}\eta_{1}(n)(1\vee\frac{M_{n}\eta_{1}(n)}{\delta^{\ell}})$ and
$\eta_{2}(n)$ come from the plug-in error of $\widehat{\bm{\beta}}$ and
$\widehat{\bm{v}}$, respectively.

##  Appendix B Proofs of Main Results

###  B.1 Proof of Lemma 1

Recall that

$\displaystyle\nabla\bar{R}_{\delta}^{i}(\bm{\beta}^{*})=w(y_{i})\frac{y_{i}\bm{z}_{i}}{\delta}K\bigg{(}\frac{y_{i}(x_{i}-\bm{\beta}^{*T}\bm{z}_{i})}{\delta}\bigg{)}.$
(B.1)

Let
$T_{i}=\frac{\sqrt{\delta}\bm{v}^{*T}(\nabla\bar{R}_{\delta}^{i}(\bm{\beta}^{*})-\nabla
R_{\delta}(\bm{\beta}^{*}))}{\sqrt{\bm{v}^{*T}\bm{\Sigma}^{*}\bm{v}^{*}}}$, we
know by definition $\mathbb{E}T_{i}=0$. Consider

$\displaystyle Var[\bm{v}^{*T}\nabla\bar{R}_{\delta}^{i}(\bm{\beta}^{*})]=$
$\displaystyle\mathbb{E}[(\bm{v}^{*T}\nabla\bar{R}_{\delta}^{i}(\bm{\beta}^{*}))^{2}]-(\bm{v}^{*T}\nabla
R_{\delta}(\bm{\beta}^{*}))^{2}.$ (B.2)

Here

$\displaystyle\mathbb{E}[(\bm{v}^{*T}\nabla\bar{R}_{\delta}^{i}(\bm{\beta}^{*}))^{2}]=$
$\displaystyle\sum_{y\in\\{-1,1\\}}w(y)^{2}\int\frac{(\bm{v}^{*T}\bm{z})^{2}}{\delta^{2}}\int
K^{2}(\frac{x-\bm{\beta}^{*T}\bm{z}}{\delta})f(x|y,\bm{z})dxf(y,\bm{z})d\bm{z}$
(B.3) $\displaystyle=$
$\displaystyle\sum_{y\in\\{-1,1\\}}w(y)^{2}\int\frac{(\bm{v}^{*T}\bm{z})^{2}}{\delta}\int
K^{2}(u)f(u\delta+\bm{\beta}^{*T}\bm{z}|y,\bm{z})duf(y,\bm{z})d\bm{z}$
$\displaystyle=$
$\displaystyle\frac{1}{\delta}\sum_{y\in\\{-1,1\\}}w(y)^{2}\int(\bm{v}^{*T}\bm{z})^{2}\int
K^{2}(u)(f(\bm{\beta}^{*T}\bm{z}|y,\bm{z})+u\delta
f^{\prime}(\bm{\beta}^{*T}\bm{z}|y,\bm{z})+o(\delta))duf(y,\bm{z})d\bm{z}$
$\displaystyle=$
$\displaystyle\frac{1}{\delta}(\bm{v}^{*T}\bm{\Sigma}^{*}\bm{v}^{*}(1+o(1))),$

where the second equality is due to a change of variable, and the third
equality is due to Assumption 1. Meanwhile, we know

$\displaystyle\bm{v}^{*T}\nabla R_{\delta}(\bm{\beta}^{*})=$
$\displaystyle\sum_{y\in\\{-1,1\\}}w(y)y\int\frac{(\bm{v}^{*T}\bm{z})}{\delta}\int
K(\frac{x-\bm{\beta}^{*T}\bm{z}}{\delta})f(x|y,\bm{z})dxf(y,\bm{z})d\bm{z}$
(B.4) $\displaystyle=$
$\displaystyle\sum_{y\in\\{-1,1\\}}w(y)y\int(\bm{v}^{*T}\bm{z})\int
K(u)f(u\delta+\bm{\beta}^{*T}\bm{z}|y,\bm{z})duf(y,\bm{z})d\bm{z}$
$\displaystyle=$
$\displaystyle\sum_{y\in\\{-1,1\\}}w(y)y\int(\bm{v}^{*T}\bm{z})\int
K(u)\frac{(u\delta)^{\ell}}{\ell!}\bigg{(}f^{(\ell)}(\bm{\beta}^{*T}\bm{z}|y,\bm{z})+\mathcal{O}((u\delta)^{\zeta})\bigg{)}duf(y,\bm{z})d\bm{z}$
$\displaystyle=$ $\displaystyle\delta^{\ell}\bm{v}^{*T}\bm{b}^{*}(1+o(1)).$

This together with Assumption 5 implies that $Var[\bm{v}^{*T}\nabla
R_{\delta}^{i}(\bm{\beta}^{*})]=\frac{1}{\delta}\bm{v}^{*T}\bm{\Sigma}^{*}\bm{v}^{*}(1+o(1))$
and therefore $Var(T_{i})=1+o(1).$ Now we verify the Lyapunov condition

$\displaystyle\frac{1}{n^{3/2}}\sum_{i}^{n}\mathbb{E}|T_{i}|^{3}=$
$\displaystyle\frac{1}{n^{3/2}}\sum_{i}^{n}\mathbb{E}\bigg{|}\frac{\sqrt{\delta}\bm{v}^{*T}(\nabla\bar{R}_{\delta}^{i}(\bm{\beta}^{*})-\nabla
R_{\delta}(\bm{\beta}^{*}))}{\sqrt{\bm{v}^{*T}\bm{\Sigma}^{*}\bm{v}^{*}}}\bigg{|}^{3}.$
(B.5)

By Assumption 5,

$\displaystyle\mathbb{E}\bigg{|}\frac{\sqrt{\delta}\bm{v}^{*T}(\nabla\bar{R}_{\delta}^{i}(\bm{\beta}^{*})-\nabla
R_{\delta}(\bm{\beta}^{*}))}{\sqrt{\bm{v}^{*T}\bm{\Sigma}^{*}\bm{v}^{*}}}\bigg{|}^{3}\lesssim$
$\displaystyle\delta^{3/2}\mathbb{E}\bigg{|}\bm{v}^{*T}(\nabla\bar{R}_{\delta}^{i}(\bm{\beta}^{*})-\nabla
R_{\delta}(\bm{\beta}^{*}))\bigg{|}^{3}$ (B.6) $\displaystyle\lesssim$
$\displaystyle\delta^{3/2}(\mathbb{E}|\bm{v}^{*T}\nabla\bar{R}_{\delta}^{i}(\bm{\beta}^{*})|^{3}+|\bm{v}^{*T}\nabla
R_{\delta}(\bm{\beta}^{*})|^{3}),$

where we can show that

$\displaystyle\mathbb{E}|\bm{v}^{*T}\nabla\bar{R}_{\delta}^{i}(\bm{\beta}^{*})|^{3}\leq$
$\displaystyle||\bm{v}^{*}||_{1}^{3}\sum_{y\in\\{-1,1\\}}\int\delta\max_{j}\bigg{|}\frac{y\bm{z}_{j}}{\delta}K(u)\bigg{|}^{3}f(u\delta+\bm{\beta}^{*T}\bm{z},y,\bm{z})dud\bm{z}$
(B.7) $\displaystyle\lesssim$
$\displaystyle||\bm{v}^{*}||_{1}^{3}M_{n}^{3}/\delta^{2},$

and from Lemma 4 it’s easy to see that the first term on the RHS of (B.6) is
dominant. This implies that

$\displaystyle\frac{1}{n^{3/2}}\sum_{i}^{n}\mathbb{E}|T_{i}|^{3}=\mathcal{O}((||\bm{v}^{*}||_{1}M_{n})^{3}/(n\delta)^{1/2}).$
(B.8)

Therefore, under the condition of this Lemma the Lyapunov condition holds,
which completes the proof by applying Lindeberg Feller Central Limit Theorem.

###  B.2 Proof of Theorem 1

It suffices to show that
$(n\delta)^{1/2}|\widehat{S}_{\delta}(\widehat{\bm{\beta}}_{0})-S_{\delta}(\bm{\beta}^{*})|=o_{\mathbb{P}}(1)$
where $\widehat{\bm{\beta}}_{0}=(0,\widehat{\gamma})$. Here we only show that
$(n\delta)^{1/2}|\widehat{S}_{\delta}^{(1)}(\widehat{\bm{\beta}}_{0}^{(2)})-S_{\delta}(\bm{\beta}^{*})|=o_{\mathbb{P}}(1)$
and the desired result shall follow naturally from Lemma 1.

By definition, we have

$\displaystyle(n\delta)^{1/2}|\widehat{S}_{\delta}^{(1)}(\widehat{\bm{\beta}}_{0}^{(2)})-S_{\delta}(\bm{\beta}^{*})|$
(B.9) $\displaystyle=$
$\displaystyle(n\delta)^{1/2}|\widehat{\bm{v}}^{(1)T}\nabla
R_{\delta}^{n_{(1)}}(\widehat{\bm{\beta}}_{0}^{(2)})-\bm{v}^{*T}\nabla
R_{\delta}^{n_{(1)}}(\bm{\beta}^{*})|$ $\displaystyle\leq$
$\displaystyle(n\delta)^{1/2}|\bm{v}^{*T}(\nabla
R_{\delta}^{n_{(1)}}(\widehat{\bm{\beta}}_{0}^{(2)})-\nabla
R_{\delta}^{n_{(1)}}(\bm{\beta}^{*}))|+(n\delta)^{1/2}|(\widehat{\bm{v}}^{(1)}-\bm{v}^{*})\nabla
R_{\delta}^{n_{(1)}}(\widehat{\bm{\beta}}_{0}^{(2)})|$ $\displaystyle:=$
$\displaystyle I_{1}+I_{2}.$

The fact that $\bm{v}^{*T}\nabla^{2}_{\cdot\bm{\gamma}}R(\bm{\beta}^{*})=0$
and Lemma 7 implies that

$\displaystyle I_{1}=$
$\displaystyle(n\delta)^{1/2}|\bm{v}^{*T}\nabla^{2}_{\cdot\bm{\gamma}}R^{n_{(1)}}_{\delta}(\bm{\beta}^{*})(\widehat{\bm{\beta}}^{(2)}-\bm{\beta}^{*})|+o_{\mathbb{P}}(1)$
(B.10) $\displaystyle\leq$
$\displaystyle(n\delta)^{1/2}||\bm{v}^{*}||_{1}\big{(}||\nabla^{2}_{\cdot\bm{\gamma}}R^{n_{(1)}}_{\delta}(\bm{\beta}^{*})-\nabla^{2}_{\cdot\bm{\gamma}}R_{\delta}(\bm{\beta}^{*})||_{\max}+||\nabla^{2}_{\cdot\bm{\gamma}}R_{\delta}(\bm{\beta}^{*})-\nabla^{2}_{\cdot\bm{\gamma}}R(\bm{\beta}^{*})||_{\max}\big{)}||\widehat{\bm{\beta}}^{(2)}-\bm{\beta}^{*}||_{1}+o_{\mathbb{P}}(1)$
$\displaystyle\lesssim$
$\displaystyle(n\delta)^{1/2}||\bm{v}^{*}||_{1}\frac{\eta_{1}(n)}{\delta}(\sqrt{\frac{\log(d)}{n\delta}}+\delta^{\ell})+o_{\mathbb{P}}(1)$
$\displaystyle=$ $\displaystyle o_{\mathbb{P}}(1).$

Similarly, since $\nabla R(\bm{\beta}^{*})=0$, imply that

$\displaystyle
I_{2}=(n\delta)^{1/2}|(\widehat{\bm{v}}^{(1)}-\bm{v}^{*})(\nabla
R_{\delta}^{n_{(1)}}(\widehat{\bm{\beta}}_{0}^{(2)})-\nabla
R(\bm{\beta}^{*}))|$ (B.11) $\displaystyle\leq$
$\displaystyle(n\delta)^{1/2}||\widehat{\bm{v}}^{(1)}-\bm{v}^{*}||_{1}||\nabla
R_{\delta}^{n_{(1)}}(\widehat{\bm{\beta}}_{0}^{(2)})-\nabla
R(\bm{\beta}^{*})||_{\infty}$ $\displaystyle\leq$
$\displaystyle(n\delta)^{1/2}||\widehat{\bm{v}}^{(1)}-\bm{v}^{*}||_{1}\bigg{(}||\nabla
R_{\delta}^{n_{(1)}}(\widehat{\bm{\beta}}_{0}^{(2)})-\nabla
R_{\delta}(\widehat{\bm{\beta}}_{0}^{(2)})||_{\infty}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+||\nabla
R_{\delta}(\widehat{\bm{\beta}}_{0}^{(2)})-\nabla
R(\widehat{\bm{\beta}}_{0}^{(2)})||_{\infty}+||\nabla
R(\widehat{\bm{\beta}}_{0}^{(2)})-\nabla
R(\bm{\beta}^{*})||_{\infty}\bigg{)}.$

Since $\widehat{\bm{\beta}}_{0}^{(2)}$ depends on the set of samples that is
disjoint with $\mathcal{N}_{1}$, Lemma 5 together with Lemma 13 implies that

$||\nabla R_{\delta}^{n_{(1)}}(\widehat{\bm{\beta}}_{0}^{(2)})-\nabla
R_{\delta}(\widehat{\bm{\beta}}_{0}^{(2)})||_{\infty}=\mathcal{O}_{\mathbb{P}}(\sqrt{\frac{\log(d)}{n\delta}}).$

This in combine with Lemma 4 and 6 further implies that

$I_{2}\lesssim(n\delta)^{1/2}||\bm{v}^{*}||_{1}\eta_{2}(n)\bigg{(}M_{n}\eta_{1}(n)+\sqrt{\frac{\log(d)}{n\delta}}+\delta^{\ell}\bigg{)}=o_{\mathbb{P}}(1).$

Putting all above together with the results from Lemma 11 and Lemma 12 as well
as the Slutsky’s theorem, we obtain that the bias corrected decorrelated score
statistic $\widehat{U}_{n}\stackrel{{\scriptstyle d}}{{\rightarrow}}N(0,1)$
under the null hypothesis.

###  B.3 Proof of Theorem 2

Recall that
$\widehat{S}_{\delta}(0,\widehat{\bm{\gamma}})=\frac{1}{2}\big{(}\widehat{S}^{(1)}_{\delta}(0,\widehat{\bm{\gamma}}^{(2)})+\widehat{S}^{(2)}_{\delta}(0,\widehat{\bm{\gamma}}^{(1)})\big{)}$.
Let’s focus on $\widehat{S}^{(1)}_{\delta}(0,\widehat{\bm{\gamma}}^{(2)})$. By
definition

$\displaystyle\widehat{S}^{(1)}(0,\widehat{\bm{\gamma}}^{(2)})=\widehat{\bm{v}}^{(1)T}\nabla
R_{\delta}^{n_{(1)}}(0,\widehat{\bm{\gamma}}^{(2)})$ (B.12) $\displaystyle=$
$\displaystyle\underbrace{(\widehat{\bm{v}}^{(1)}-\bm{v}^{*})^{T}\nabla
R_{\delta}^{n_{(1)}}(0,\widehat{\bm{\gamma}}^{(2)})}_{I_{1}}+\bm{v}^{*T}\nabla
R_{\delta}^{n_{(1)}}(0,\widehat{\bm{\gamma}}^{(2)})$ $\displaystyle=$
$\displaystyle I_{1}+\underbrace{\bm{v}^{*T}\big{(}\nabla
R_{\delta}^{n_{(1)}}(0,\widehat{\bm{\gamma}}^{(2)})-\nabla
R_{\delta}^{n_{(1)}}(0,\bm{\gamma}^{*})\big{)}}_{I_{2}}+\bm{v}^{*T}\nabla
R_{\delta}^{n_{(1)}}(0,\bm{\gamma}^{*})$ $\displaystyle=$ $\displaystyle
I_{1}+I_{2}+\bm{v}^{*T}\nabla
R_{\delta}^{n_{(1)}}(\bm{\beta}^{*})-\theta^{*}\bm{v}^{*T}\nabla^{2}_{\cdot\theta}R(\bm{\beta}^{*})+\underbrace{\bm{v}^{*T}(\nabla
R_{\delta}^{n_{(1)}}(0,\bm{\gamma}^{*})-\nabla
R_{\delta}^{n_{(1)}}(\bm{\beta}^{*}))+\theta^{*}\bm{v}^{*T}\nabla^{2}_{\cdot\theta}R(\bm{\beta}^{*})}_{I_{3}}$
$\displaystyle=$ $\displaystyle
S^{(1)}(\bm{\beta}^{*})-\theta^{*}\bm{v}^{*T}\nabla^{2}_{\cdot\theta}R(\bm{\beta}^{*})+I_{1}+I_{2}+I_{3}.$

For $I_{1}$, we know

$I_{1}\leq||\widehat{\bm{v}}^{(1)}-\bm{v}^{*}||_{1}||\nabla
R_{\delta}^{n}(0,\widehat{\bm{\gamma}})||_{\infty}=\mathcal{O}_{\mathbb{P}}(||\bm{v}^{*}||_{1}\eta_{2}(n)(\sqrt{\log(d)/(n\delta)}+\delta^{\ell}+M_{n}(\theta^{*}+\eta_{1}(n))))=o_{\mathbb{P}}((n\delta)^{-1/2}).$

Following a similar proof of Lemma 7, we can show that

$I_{2}=\mathcal{O}_{\mathbb{P}}\big{(}||\bm{v}^{*}||_{1}\eta_{1}(n)\big{[}\delta^{\ell-1}+\sqrt{\log(d)/(n\delta^{3})}+M_{n}\theta^{*}\big{]}\big{)}=o_{\mathbb{P}}((n\delta)^{-1/2}).$

Finally, for $I_{3}$, since by definition

$\displaystyle\bm{v}^{*T}\nabla^{2}_{\cdot\theta}R(\bm{\beta}^{*})=\nabla^{2}_{\theta|\bm{\gamma}}R(\bm{\beta}^{*}),$
(B.13)

Lemma 9 implies that $I_{3}=o_{\mathbb{P}}((n\delta)^{-1/2})$. Put all pieces
together, we have shown that

$|\widehat{S}^{(1)}(0,\widehat{\bm{\gamma}}^{(2)})-(S^{(1)}(\bm{\beta}^{*})-\theta^{*}\nabla^{2}_{\theta|\bm{\gamma}}R(\bm{\beta}^{*}))|=o_{\mathbb{P}}((n\delta)^{-1/2}).$

A similar result will also hold for
$\widehat{S}^{(2)}(0,\widehat{\bm{\gamma}}^{(1)})$, and thus we conclude

$\displaystyle|\widehat{S}(0,\widehat{\bm{\gamma}})-(S(\bm{\beta}^{*})-\theta^{*}\nabla^{2}_{\theta|\bm{\gamma}}R(\bm{\beta}^{*}))|=o_{\mathbb{P}}((n\delta)^{-1/2}).$
(B.14)

At this point the conclusion of (3.9) is shown. To show (3.10), the above
formula also implies that

$\displaystyle\mathbb{P}(|(n\delta)^{1/2}\widehat{S}(0,\widehat{\bm{\gamma}})/\sigma^{*}|\leq
t)\leq\mathbb{P}(L(n)\leq(n\delta)^{1/2}S(\bm{\beta}^{*})/\sigma^{*}\leq
U(n))+o(1),$ (B.15)

where
$L(n)=-t-q(n)+\widetilde{C}\sigma^{*-1}\nabla^{2}_{\theta|\bm{\gamma}}R(\bm{\beta}^{*})n^{-\phi}(n\delta)^{1/2}$,
$U(n)=t+q(n)+\widetilde{C}\sigma^{*-1}\nabla^{2}_{\theta|\bm{\gamma}}R(\bm{\beta}^{*})n^{-\phi}(n\delta)^{1/2}$,
and $q(n)=o(1)$ is some deterministic sequence. Since
$n^{-\phi}(n\delta)^{1/2}\nabla^{2}_{\theta|\bm{\gamma}}R(\bm{\beta}^{*})\sigma^{*-1}\rightarrow\infty$,
it is easily seen that
$\mathbb{P}(|(n\delta)^{1/2}\widehat{S}(0,\widehat{\bm{\gamma}})/\sigma^{*}|\leq
t)\rightarrow 0$. Since $\widehat{\sigma}$ is consistent, for $n$ large enough
we will get $|\widehat{\sigma}/\sigma-1|\leq 3$, which finally implies the
desired result. This completes the proof.

###  B.4 Proof of Lemma 2

It suffices to show the rate for
$\widehat{\bm{v}}^{(1)T}\widehat{\bm{\Gamma}}^{(1)}(\delta)\widehat{\bm{v}}^{(1)}$.
Let
$\widetilde{\bm{\Gamma}}(\delta)=\mathbb{E}\bigg{[}\nabla\bar{R}_{\delta}^{1}(\widehat{\bm{\beta}}^{(2)})\nabla\bar{R}_{\delta}^{1}(\widehat{\bm{\beta}}^{(2)})^{T}\bigg{]}$,
and
$\bm{\Gamma}(\delta)=\mathbb{E}[\nabla\bar{R}^{1}_{\delta}(\bm{\beta}^{*})\nabla\bar{R}^{1}_{\delta}(\bm{\beta}^{*})^{T}]$.
We firstly look at

$\displaystyle||\widehat{\bm{\Gamma}}^{(1)}(\delta)-\bm{\Gamma}(\delta)||_{\max}\leq||\widehat{\bm{\Gamma}}^{(1)}(\delta)-\widetilde{\bm{\Gamma}}(\delta)||_{\max}+||\widetilde{\bm{\Gamma}}(\delta)-\bm{\Gamma}(\delta)||_{\max}.$
(B.16)

Now we bound the first term. For any $\delta\in\Delta$, we can always find a
$\widetilde{\delta}\in\widetilde{\Delta}$ such that
$|\delta-\widetilde{\delta}|\lesssim n^{-\rho}$, where $\rho$ is some constant
that can be arbitrarily large and the set $\widetilde{\Delta}$ that is a
subset of $\Delta$ has cardinality of order $\mathcal{O}(n^{\rho})$. For such
pair $(\delta,\widetilde{\delta})$, we can show that

$\displaystyle\sqrt{n\delta^{3}}||\widehat{\bm{\Gamma}}^{(1)}(\delta)-\widetilde{\bm{\Gamma}}(\delta)||_{\max}$
(B.17) $\displaystyle\lesssim$
$\displaystyle\sqrt{n\widetilde{\delta}^{3}}\bigg{(}||\widehat{\bm{\Gamma}}^{(1)}(\widetilde{\delta})-\widetilde{\bm{\Gamma}}(\widetilde{\delta})||_{\max}+||\widehat{\bm{\Gamma}}^{(1)}(\delta)-\widehat{\bm{\Gamma}}^{(1)}(\widetilde{\delta})||_{\max}+||\widetilde{\bm{\Gamma}}(\widetilde{\delta})-\widetilde{\bm{\Gamma}}(\delta)||_{\max}\bigg{)}$
$\displaystyle\lesssim$
$\displaystyle\sqrt{n\widetilde{\delta}^{3}}||\widehat{\bm{\Gamma}}^{(1)}(\widetilde{\delta})-\widetilde{\bm{\Gamma}}(\widetilde{\delta})||_{\max}+\mathcal{O}(n^{-1/2}),$

where the first inequality is by taking $\rho$ large enough. To see last step
of above, recall that under Assumptions 2, 3 and 4

$\displaystyle||\widehat{\bm{\Gamma}}^{(1)}(\delta)-\widehat{\bm{\Gamma}}^{(1)}(\widetilde{\delta})||_{\max}=$
$\displaystyle\max_{j,k}\bigg{|}\frac{1}{|\mathcal{N}_{1}|}\sum_{i\in\mathcal{N}_{1}}w(y_{i})^{2}\bm{z}_{ij}\bm{z}_{ik}\bigg{[}\frac{K^{2}((x_{i}-\widehat{\bm{\beta}}^{(2)T}\bm{z}_{i})/\delta)}{\delta^{2}}-\frac{K^{2}((x_{i}-\widehat{\bm{\beta}}^{(2)T}\bm{z}_{i})/\widetilde{\delta})}{\widetilde{\delta}^{2}}\bigg{]}\bigg{|}$
(B.18) $\displaystyle\lesssim$ $\displaystyle
M_{n}^{2}(\frac{1}{\delta^{2}}-\frac{1}{\widetilde{\delta}^{2}}).$

Thus, taking $\rho$ large enough will ensure
$\sqrt{n\widetilde{\delta}^{3}}(1+o(1))||\widehat{\bm{\Gamma}}^{(1)}(\delta)-\widehat{\bm{\Gamma}}^{(1)}(\widetilde{\delta})||_{\max}=\mathcal{O}(n^{-1/2})$.
With a similar derivation we will also obtain
$\sqrt{n\widetilde{\delta}^{3}}(1+o(1))||\widetilde{\bm{\Gamma}}(\delta)-\widetilde{\bm{\Gamma}}(\widetilde{\delta})||_{\max}=\mathcal{O}(n^{-1/2})$,
and thus (LABEL:eq_var_adapt_1) holds.

Now we start to bound
$\max_{\widetilde{\delta}\in\widetilde{\Delta}}\sqrt{n\widetilde{\delta}^{3}}||\widehat{\bm{\Gamma}}^{(1)}(\widetilde{\delta})-\widetilde{\bm{\Gamma}}(\widetilde{\delta})||_{\max}$.
Since

$|\widehat{\bm{\Gamma}}^{(1)}_{ijk}(\widetilde{\delta})|=|(\nabla\bar{R}_{\widetilde{\delta}}^{i}(\widehat{\bm{\beta}}^{(2)})\nabla\bar{R}_{\widetilde{\delta}}^{i}(\widehat{\bm{\beta}}^{(2)})^{T})_{jk}|\lesssim\frac{M_{n}^{2}}{\widetilde{\delta}^{2}}$

and

$Var[(\nabla\bar{R}_{\widetilde{\delta}}^{i}(\widehat{\bm{\beta}}^{(2)})\nabla\bar{R}_{\widetilde{\delta}}^{i}(\widehat{\bm{\beta}}^{(2)})^{T})_{jk}]\lesssim\frac{1}{\widetilde{\delta}^{3}},$

by applying Bernstein inequality similar to the proof of Lemma 5 , conditioned
on $\widehat{\bm{\beta}}^{(2)}$ we can show that

$\displaystyle\mathbb{P}_{\widehat{\bm{\beta}}^{(2)}}(\max_{\widetilde{\delta}\in\widetilde{\Delta}}\sqrt{n\widetilde{\delta}^{3}}||\widehat{\bm{\Gamma}}^{(1)}(\widetilde{\delta})-\widetilde{\bm{\Gamma}}(\widetilde{\delta})||_{\max}>t)$
(B.19) $\displaystyle\leq$
$\displaystyle\sum_{\widetilde{\delta}\in\widetilde{\Delta}}\sum_{j=1}^{d}\sum_{k=1}^{d}\mathbb{P}_{\widehat{\bm{\beta}}^{(2)}}(\sqrt{n\widetilde{\delta}^{3}}|\widehat{\bm{\Gamma}}^{(1)}(\widetilde{\delta})-\widetilde{\bm{\Gamma}}(\widetilde{\delta})|_{jk}>t)$
$\displaystyle\lesssim$ $\displaystyle
n^{\rho}d^{2}\exp\bigg{(}-\frac{\frac{1}{2}t^{2}/\widetilde{\delta}^{3}}{\frac{c_{0}}{\widetilde{\delta}^{3}}+\frac{c1M_{n}^{2}t}{3\widetilde{\delta}^{3}\sqrt{n\widetilde{\delta}}}}\bigg{)},$

where $c_{0},c_{1}>0$ are some constants. By taking $t=c_{2}\sqrt{\log(nd)}$
for some constant $c_{2}$ large enough, we obtain that conditioned on
$\widehat{\bm{\beta}}^{(2)}$, with probability greater than
$1-\mathcal{O}((nd)^{-1})$

$\displaystyle||\widehat{\bm{\Gamma}}^{(1)}(\widetilde{\delta})-\widetilde{\bm{\Gamma}}(\widetilde{\delta})||_{\max}\lesssim\sqrt{\frac{\log(nd)}{n\widetilde{\delta}^{3}}}$
(B.20)

uniformly over $\widetilde{\delta}\in\widetilde{\Delta}$. This combines with
Lemma 13 further implies that with probability approaching to 1

$\displaystyle||\widehat{\bm{\Gamma}}^{(1)}(\delta)-\widetilde{\bm{\Gamma}}(\delta)||_{\max}\lesssim\sqrt{\frac{\log(nd)}{n\delta^{3}}}$
(B.21)

uniformly over $\delta\in\Delta$.

For the second term on the RHS of (B.16), similar to Lemma 6, we can show that

$\displaystyle||\widetilde{\bm{\Gamma}}(\delta)-\bm{\Gamma}(\delta)||_{\max}\lesssim\frac{M_{n}\eta_{1}(n)}{\delta}.$
(B.22)

This implies that

$\displaystyle||\widehat{\bm{\Gamma}}(\delta)-\bm{\Gamma}(\delta)||_{\max}=\frac{1}{\delta}\mathcal{O}_{\mathbb{P}}(\sqrt{\frac{\log(nd)}{n\delta}}\vee
M_{n}\eta_{1}(n)).$ (B.23)

Meanwhile, since

$\displaystyle||\delta\cdot\bm{\Gamma}(\delta)||_{\max}=$
$\displaystyle\max_{i,j}\bigg{|}\sum_{y}w(y)^{2}\int\frac{z_{i}z_{j}}{\delta}\int
K^{2}(\frac{y(x-\bm{\beta}^{*T}\bm{z})}{\delta})f(x|y,\bm{z})dxf(y,\bm{z})d\bm{z}\bigg{|}$
(B.24) $\displaystyle=$ $\displaystyle\mathcal{O}(1),$

under the condition that $\sqrt{\frac{\log(nd)}{n\delta}}\vee
M_{n}\eta_{1}(n)=o(1)$, we know
$||\delta\widehat{\bm{\Gamma}}(\delta)||_{\infty}=\mathcal{O}_{\mathbb{P}}(1)$.
Finally, by triangle inequality

$\displaystyle|\widehat{\bm{v}}^{(1)T}\widehat{\bm{\Gamma}}^{(1)}(\delta)\widehat{\bm{v}}^{(1)}-\bm{v}^{*T}{\bm{\Gamma}}(\delta)\bm{v}^{*}|$
(B.25) $\displaystyle\leq$
$\displaystyle||\widehat{\bm{v}}^{(1)}-\bm{v}^{*}||_{1}||\widehat{\bm{\Gamma}}^{(1)}(\delta)(\widehat{\bm{v}}^{(1)}-\bm{v})||_{\infty}+2||\bm{v}^{*T}\widehat{\bm{\Gamma}}^{(1)}(\delta)||_{\infty}||\widehat{\bm{v}}^{(1)}-\bm{v}^{*}||_{1}+|\bm{v}^{*T}(\widehat{\bm{\Gamma}}^{(1)}(\delta)-{\bm{\Gamma}}(\delta))\bm{v}^{*}|,$

which further implies that with probability approaching 1

$|\widehat{\bm{v}}^{(1)T}\widehat{\bm{\Gamma}}^{(1)}(\delta)\widehat{\bm{v}}^{(1)}-\bm{v}^{*T}{\bm{\Gamma}}(\delta)\bm{v}^{*}|\lesssim\frac{||\bm{v}^{*}||_{1}^{2}}{\delta}\bigg{(}\eta_{2}(n)+\sqrt{\frac{\log(nd)}{n\delta}}+M_{n}\eta_{1}(n)\bigg{)}$

uniformly over $\delta\in\Delta$. This completes the proof.

###  B.5 Proof of Lemma 3

It suffices to bound
$\bigg{|}\bigg{(}\widehat{\bm{v}}^{(1)T}\frac{1}{|\mathcal{N}_{1}|}\sum_{i\in\mathcal{N}_{1}}A_{i}(\widehat{\bm{\beta}}^{(2)},\delta)\bigg{)}^{2}-SB(\delta)\bigg{|}$.
Define

$\displaystyle\widetilde{A}(\bm{\beta},\delta)=$
$\displaystyle\mathbb{E}_{\bm{\beta}}(A_{1}(\bm{\beta},\delta))$ (B.26)
$\displaystyle=$
$\displaystyle\sum_{y}w(y)y\int_{\bm{z}}\bm{z}\int_{t}J(t)\int_{u}K(u)[f(u\delta+tb+\bm{\beta}^{T}\bm{z}|y,\bm{z})-f(tb+\bm{\beta}^{T}\bm{z}|y,\bm{z})]dudtf(y,\bm{z})d\bm{z}.$

By definition,

$\displaystyle\bigg{|}\bigg{(}\widehat{\bm{v}}^{(1)T}\frac{1}{|\mathcal{N}_{1}|}\sum_{i\in\mathcal{N}_{1}}A_{i}(\widehat{\bm{\beta}}^{(2)},\delta)\bigg{)}^{2}-SB(\delta)\bigg{|}$
(B.27) $\displaystyle=$
$\displaystyle\bigg{|}\bigg{(}\widehat{\bm{v}}^{(1)T}\frac{1}{|\mathcal{N}_{1}|}\sum_{i\in\mathcal{N}_{1}}A_{i}(\widehat{\bm{\beta}}^{(2)},\delta)\bigg{)}^{2}-\big{(}\bm{v}^{*T}A(\bm{\beta}^{*},\delta)\big{)}^{2}\bigg{|}$
$\displaystyle\leq$
$\displaystyle\bigg{|}\bigg{(}\widehat{\bm{v}}^{(1)T}\frac{1}{|\mathcal{N}_{1}|}\sum_{i\in\mathcal{N}_{1}}A_{i}(\widehat{\bm{\beta}}^{(2)},\delta)\bigg{)}^{2}-\big{(}\widehat{\bm{v}}^{(1)T}A(\widehat{\bm{\beta}}^{(2)},\delta)\big{)}^{2}\bigg{|}+\bigg{|}\big{(}\widehat{\bm{v}}^{(1)T}A(\widehat{\bm{\beta}}^{(2)},\delta)\big{)}^{2}-\big{(}\bm{v}^{*T}A(\bm{\beta}^{*},\delta)\big{)}^{2}\bigg{|}$
$\displaystyle\leq$
$\displaystyle\underbrace{\bigg{|}\bigg{(}\widehat{\bm{v}}^{(1)T}\big{(}\frac{1}{|\mathcal{N}_{1}|}\sum_{i\in\mathcal{N}_{1}}A_{i}(\widehat{\bm{\beta}}^{(2)},\delta)-\widetilde{A}(\widehat{\bm{\beta}}^{(2)},\delta)\big{)}\bigg{)}^{2}\bigg{|}}_{I_{1}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}+2\underbrace{\bigg{|}\bigg{(}\widehat{\bm{v}}^{(1)T}\widetilde{A}(\widehat{\bm{\beta}}^{(2)},\delta)\cdot\widehat{\bm{v}}^{(1)T}\big{(}\frac{1}{|\mathcal{N}_{1}|}\sum_{i\in\mathcal{N}_{1}}A_{i}(\widehat{\bm{\beta}}^{(2)},\delta)-\widetilde{A}(\widehat{\bm{\beta}}^{(2)},\delta)\big{)}\bigg{|}}_{I_{2}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}+\underbrace{\bigg{|}\big{(}\widehat{\bm{v}}^{(1)T}\widetilde{A}(\widehat{\bm{\beta}}^{(2)},\delta)\big{)}^{2}-\big{(}\widehat{\bm{v}}^{(1)T}A(\widehat{\bm{\beta}}^{(2)},\delta)\big{)}^{2}\bigg{|}}_{I_{3}}+\underbrace{\bigg{|}\big{(}\widehat{\bm{v}}^{(1)T}A(\widehat{\bm{\beta}}^{(2)},\delta)\big{)}^{2}-\big{(}\bm{v}^{*T}A(\bm{\beta}^{*},\delta)\big{)}^{2}\bigg{|}}_{I_{4}}.$

Now we start to bound each term.

  * •

For $I_{1}$, we firstly study
$\delta^{-\ell}||\frac{1}{|\mathcal{N}_{1}|}\sum_{i\in\mathcal{N}_{1}}A_{i}(\widehat{\bm{\beta}}^{(2)},\delta)-\widetilde{A}(\widehat{\bm{\beta}}^{(2)},\delta)||_{\infty}$
condition on $\widehat{\bm{\beta}}^{(2)}$. Similar to the proof of Lemma 12,
it suffices to bound

$\displaystyle\widetilde{\delta}^{-\ell}||\frac{1}{|\mathcal{N}_{1}|}\sum_{i\in\mathcal{N}_{1}}A_{i}(\widehat{\bm{\beta}}^{(2)},\widetilde{\delta})-\widetilde{A}(\widehat{\bm{\beta}}^{(2)},\widetilde{\delta})||_{\infty}$
(B.28)

uniformly over $\widetilde{\delta}\in\widetilde{\Delta}$, where
$\widetilde{\Delta}$ is a subset of $\Delta$ with cardinality of order
$\mathcal{O}(n^{\rho})$ for some constant $\rho>0$, such that for each
$\delta\in\Delta$, there exist $\widetilde{\delta}\in\widetilde{\Delta}$ with
$|\delta-\widetilde{\delta}|\lesssim n^{-\rho}$.

Let
$T_{ij}(\widetilde{\delta})=\widetilde{\delta}^{-\ell}\big{(}A_{i}(\widehat{\bm{\beta}}^{(2)},\widetilde{\delta})-\widetilde{A}(\widehat{\bm{\beta}}^{(2)},\widetilde{\delta})\big{)}_{j}$,
we know $\mathbb{E}_{\widehat{\bm{\beta}}^{(2)}}T_{ij}(\widetilde{\delta})=0$.
Moreover for all $j$,

$\displaystyle\widetilde{\delta}^{-\ell}|(A_{i}(\widehat{\bm{\beta}}^{(2)},\widetilde{\delta}))_{j}|=$
$\displaystyle\widetilde{\delta}^{-\ell}\bigg{|}w(y_{i})\frac{z_{ij}y_{i}}{b}\int
K(u)\big{[}J(\frac{x_{i}-\widehat{\bm{\beta}}^{(2)T}\bm{z}_{i}-u\widetilde{\delta}}{b})-J(\frac{x_{i}-\widehat{\bm{\beta}}^{(2)T}\bm{z}_{i}}{b})\big{]}\biggr{\\}}du\bigg{|}$
(B.29) $\displaystyle=$
$\displaystyle\widetilde{\delta}^{-\ell}\bigg{|}w(y_{i})\frac{z_{ij}y_{i}}{b}\int
K(u)\frac{\widetilde{\delta}^{\ell}}{b^{\ell}}\frac{u^{\ell}}{\ell!}J^{\ell}(\frac{x_{i}-\widehat{\bm{\beta}}^{(2)T}\bm{z}_{i}-\tau
u\widetilde{\delta}}{b})\big{]}\biggr{\\}}du\bigg{|}$ $\displaystyle\lesssim$
$\displaystyle\frac{M_{n}}{b^{\ell+1}},$

where $\tau\in[0,1]$ and the second equality is because $K$ is of order $\ell$
and mean value theorem. This implies that
$|T_{ij}(\widetilde{\delta})|\lesssim\frac{M_{n}}{b^{\ell+1}}$. Now we look at
$\mathbb{E}_{\widehat{\bm{\beta}}^{(2)}}[T^{2}_{ij}(\widetilde{\delta})]$.

Notice that for all $j$

$\displaystyle\widetilde{A}(\widehat{\bm{\beta}}^{(2)},\widetilde{\delta})_{j}=$
$\displaystyle\sum_{y}w(y)\int_{\bm{z}}\bm{z}_{j}y\int_{t}J(t)\int_{u}K(u)[f(u\widetilde{\delta}+tb+\widehat{\bm{\beta}}^{(2)T}\bm{z}|y,\bm{z})-f(tb+\widehat{\bm{\beta}}^{(2)T}\bm{z}|y,\bm{z})]dudtf(y,\bm{z})d\bm{z}$
(B.30) $\displaystyle=$
$\displaystyle\sum_{y}w(y)\int_{\bm{z}}\bm{z}_{j}y\int_{t}J(t)\int_{u}K(u)\frac{(u\widetilde{\delta})^{\ell}}{\ell!}f^{(\ell)}(\tau
u\widetilde{\delta}+tb+\widehat{\bm{\beta}}^{(2)T}\bm{z}|y,\bm{z})dudtf(y,\bm{z})d\bm{z}$
$\displaystyle\lesssim$ $\displaystyle\widetilde{\delta}^{\ell},$

where the second equality is because $K$ is of order $\ell$ and
$\tau\in[0,1]$. For the second moment, we have

$\displaystyle\mathbb{E}_{\widehat{\bm{\beta}}^{(2)}}[(\widetilde{\delta}^{-\ell}A_{1}(\widehat{\bm{\beta}}^{(2)},\widetilde{\delta}))_{j}^{2}]$
(B.31) $\displaystyle=$
$\displaystyle\widetilde{\delta}^{-2\ell}\mathbb{E}_{\widehat{\bm{\beta}}^{(2)}}\bigg{[}w(y)^{2}\frac{z_{j}^{2}}{b^{2}}\bigg{(}\int
K(u)\big{[}J(\frac{x-\widehat{\bm{\beta}}^{(2)T}\bm{z}-u\widetilde{\delta}}{b})-J(\frac{x-\widehat{\bm{\beta}}^{(2)T}\bm{z}}{b})\big{]}\biggr{\\}}du\bigg{)}^{2}\bigg{]}$
$\displaystyle=$
$\displaystyle\widetilde{\delta}^{-2\ell}\mathbb{E}_{\widehat{\bm{\beta}}^{(2)}}\bigg{[}w(y)^{2}\frac{z_{j}^{2}}{b^{2}}\bigg{(}\int
K(u)\frac{\widetilde{\delta}^{\ell}}{b^{\ell}}\frac{u^{\ell}}{\ell!}J^{\ell}(\frac{x-\widehat{\bm{\beta}}^{(2)T}\bm{z}-\tau
u\widetilde{\delta}}{b})du\bigg{)}^{2}\bigg{]}$ $\displaystyle=$
$\displaystyle\frac{1}{b^{2\ell+1}}\bigg{|}\sum_{y}\int_{z}w(y)^{2}z_{j}^{2}\int_{t}\bigg{(}\int
K(u)\frac{u^{\ell}}{\ell!}J^{\ell}(t-\frac{\tau
u\widetilde{\delta}}{b})du\bigg{)}^{2}f(tb+\widehat{\bm{\beta}}^{(2)T}\bm{z}|y,\bm{z})dtf(y,\bm{z})d\bm{z}\bigg{|}$
$\displaystyle\lesssim$ $\displaystyle\frac{1}{b^{2\ell+1}},$

where the last equality is due to a change of variable
$t=\frac{x-\widehat{\bm{\beta}}^{(2)T}\bm{z}}{b}$. Since $b\rightarrow 0$, we
conclude that
$\mathbb{E}_{\widehat{\bm{\beta}}^{(2)}}[T^{2}_{ij}(\widetilde{\delta})]\lesssim\frac{1}{b^{2\ell+1}}.$
Now applying Bernstein inequality similar to Lemma 12 with
$M_{n}\sqrt{\frac{\log(d)}{nb}}=\mathcal{O}(1)$, second moment of order
$\frac{1}{b^{2\ell+1}}$ and each term is bounded by
$\frac{M_{n}}{b^{\ell+1}}$, we will obtain that with probability greater than
$1-\mathcal{O}((nd)^{-1})$, conditioned on $\widehat{\bm{\beta}}^{(2)}$

$\displaystyle||\frac{1}{|\mathcal{N}_{1}|}\sum_{i\in\mathcal{N}_{1}}A_{i}(\widehat{\bm{\beta}}^{(2)},\widetilde{\delta})-\mathbb{E}A_{i}(\widehat{\bm{\beta}}^{(2)},\widetilde{\delta})||_{\infty}\lesssim\frac{\widetilde{\delta}^{\ell}}{b^{\ell}}\sqrt{\frac{\log(nd)}{nb}},$
(B.32)

which together with Lemma 13 implies that with probability approaching 1

$\displaystyle||\frac{1}{|\mathcal{N}_{1}|}\sum_{i\in\mathcal{N}_{1}}A_{i}(\widehat{\bm{\beta}}^{(2)},\delta)-\mathbb{E}_{\widehat{\bm{\beta}}^{(2)}}A_{i}(\widehat{\bm{\beta}}^{(2)},\delta)||_{\infty}\lesssim\frac{\delta^{\ell}}{b^{\ell}}\sqrt{\frac{\log(nd)}{nb}}$
(B.33)

holds uniformly for all $\delta\in\Delta$.

Therefore we conclude that

$\displaystyle
I_{1}=\mathcal{O}_{\mathbb{P}}(||\bm{v}^{*}||_{1}^{2}\frac{\delta^{2\ell}}{b^{2\ell}}\frac{\log(nd)}{nb}).$
(B.34)

  * •

For $I_{2}$, (B.30) implies that
$||\mathbb{E}_{\widehat{\bm{\beta}}^{(2)}}A_{1}(\widehat{\bm{\beta}}^{(2)},\delta)||_{\infty}\lesssim\delta^{\ell}$
regardless the choice of $\widehat{\bm{\beta}}^{(2)}$. This combined with
(B.33) further implies that

$\displaystyle
I_{2}=\mathcal{O}_{\mathbb{P}}(||\bm{v}^{*}||_{1}^{2}\delta^{\ell}\sqrt{\frac{\delta^{2\ell}}{b^{2\ell}}\frac{\log(nd)}{nb}})=\mathcal{O}_{\mathbb{P}}(||\bm{v}^{*}||_{1}^{2}\frac{\delta^{2\ell}}{b^{\ell}}\sqrt{\frac{\log(nd)}{nb}}).$
(B.35)

  * •

For $I_{3}$, by definition for all $j$ we can write

$\displaystyle\big{(}\widetilde{A}(\widehat{\bm{\beta}}^{(2)},\delta)-A(\widehat{\bm{\beta}}^{(2)},\delta)\big{)}_{j}$
(B.36) $\displaystyle=$
$\displaystyle\sum_{y}w(y)\int_{\bm{z}}\bm{z}_{j}y\int_{t}J(t)\int_{u}K(u)[(f(u\delta+tb+\widehat{\bm{\beta}}^{(2)T}\bm{z}|y,\bm{z})-f(tb+\widehat{\bm{\beta}}^{(2)T}\bm{z}|y,\bm{z}))$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}-(f(u\delta+\widehat{\bm{\beta}}^{(2)T}\bm{z}|y,\bm{z})-f(\widehat{\bm{\beta}}^{(2)T}\bm{z}|y,\bm{z}))]dudtf(y,\bm{z})d\bm{z}$
$\displaystyle=$
$\displaystyle\sum_{y}w(y)\int_{\bm{z}}\bm{z}_{j}y\int_{u}\sum_{j=\ell}^{\ell+r-1}\frac{(u\delta)^{j}}{j!}K(u)\int_{t}J(t)[f^{(j)}(tb+\widehat{\bm{\beta}}^{(2)T}\bm{z}|y,\bm{z})-f^{(j)}(\widehat{\bm{\beta}}^{(2)T}\bm{z}|y,\bm{z})]dtduf(y,\bm{z})d\bm{z}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}+\mathcal{O}(\delta^{\ell+r})$
$\displaystyle=$
$\displaystyle\sum_{y}w(y)\int_{\bm{z}}\bm{z}_{j}y\int_{u}\sum_{j=\ell}^{\ell+r-1}\frac{(u\delta)^{j}}{j!}K(u)\int_{t}J(t)\frac{(tb)^{\ell+r-j}}{(\ell+r-j)!}f^{(\ell+r)}(\tau
tb+\widehat{\bm{\beta}}^{(2)T}\bm{z}|y,\bm{z})dtduf(y,\bm{z})d\bm{z}+\mathcal{O}(\delta^{\ell+r})$
$\displaystyle=$ $\displaystyle\mathcal{O}(\delta^{\ell}(\delta\vee b)^{r}),$

where $\tau\in[0,1]$. Here the first equality is by definition, the second
equality uses the property that $\int K(u)u^{i}du=0~{}\forall~{}i<\ell$, and
the third equality uses the property that $\int J(u)u^{i}du=0~{}\forall~{}i<r$
and the mean value theorem. This implies that

$\displaystyle||\widetilde{A}(\widehat{\bm{\beta}}^{(2)},\delta)-A(\widehat{\bm{\beta}}^{(2)},\delta)||_{\infty}=\mathcal{O}(\delta^{\ell}(\delta\vee
b)^{r}).$ (B.37)

Meanwhile, by (B.30) and a similar derivation, we know both
$||\widetilde{A}(\widehat{\bm{\beta}}^{(2)},\delta)||_{\infty}$ and
$||A(\widehat{\bm{\beta}}^{(2)},\delta)||_{\infty}$ are of order
$\delta^{\ell}$. These together imply that

$\displaystyle
I_{3}=\mathcal{O}(||\bm{v}^{*}||_{1}^{2}\delta^{2\ell}(\delta\vee b)^{r}).$
(B.38)

  * •

For $I_{4}$, again
$||A(\widehat{\bm{\beta}}^{(2)},\delta)||_{\infty}\asymp||A(\bm{\beta}^{*},\delta)||_{\infty}=\mathcal{O}(\delta^{\ell})$.
For the difference, we have

$\displaystyle||A(\widehat{\bm{\beta}}^{(2)},\delta)-A(\bm{\beta}^{*},\delta)||_{\infty}$
(B.39) $\displaystyle=$
$\displaystyle\max_{j}\bigg{|}\int_{u}K(u)\bigg{[}\sum_{y}w(y)\int_{\bm{z}}z_{j}y\big{(}f(u\delta+\widehat{\bm{\beta}}^{(2)T}\bm{z}|y,\bm{z})-f(\widehat{\bm{\beta}}^{(2)T}\bm{z}|y,\bm{z})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-f(u\delta+\bm{\beta}^{*T}\bm{z}|y,\bm{z})+f(\bm{\beta}^{*T}\bm{z}|y,\bm{z})\big{)}f(y,\bm{z})d\bm{z}\bigg{]}du\bigg{|}$
$\displaystyle=$
$\displaystyle\max_{j}\bigg{|}\sum_{y}w(y)\int_{\bm{z}}z_{j}y\int_{u}K(u)\bigg{[}\big{(}f(u\delta+\widehat{\bm{\beta}}^{(2)T}\bm{z}|y,\bm{z})-f(u\delta+\bm{\beta}^{*T}\bm{z}|y,\bm{z})\big{)}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\big{(}f(\widehat{\bm{\beta}}^{(2)T}\bm{z}|y,\bm{z})-f(\bm{\beta}^{*T}\bm{z}|y,\bm{z})\big{)}\bigg{]}f(y,\bm{z})d\bm{z}du\bigg{|}$
$\displaystyle=$
$\displaystyle\max_{j}\bigg{|}\sum_{y}w(y)\int_{\bm{z}}z_{j}y\int_{u}K(u)\bigg{[}(\widehat{\bm{\beta}}-\bm{\beta}^{*})^{T}\bm{z}\big{(}f^{\prime}(u\delta+\bm{\beta}^{*T}\bm{z}|y,\bm{z})-f^{\prime}(\bm{\beta}^{*T}\bm{z}|y,\bm{z})\big{)}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+((\widehat{\bm{\beta}}-\bm{\beta}^{*})^{T}\bm{z})^{2}\big{(}f^{\prime\prime}(u\delta+\widetilde{\bm{\beta}}^{T}\bm{z})-f^{\prime\prime}(\breve{\bm{\beta}}^{T}\bm{z})\big{)}\bigg{]}duf(y,\bm{z})d\bm{z}\bigg{|}$
$\displaystyle\lesssim$ $\displaystyle
M_{n}||\widehat{\bm{\beta}}-\bm{\beta}^{*}||_{1}\delta^{\ell}+M_{n}^{2}||\widehat{\bm{\beta}}-\bm{\beta}^{*}||_{1}^{2},$

where the last step follows from an $\ell$ order Taylor expansion of
$f^{\prime}(u\delta+\bm{\beta}^{*T}\bm{z}|y,\bm{z})-f^{\prime}(\bm{\beta}^{*T}\bm{z}|y,\bm{z})$
for the first term and boundedness of $f^{\prime\prime}$ for the second term.
This implies that

$\displaystyle||A(\widehat{\bm{\beta}}^{(2)},\delta)-A(\bm{\beta}^{*},\delta)||_{\infty}\lesssim
M_{n}\eta_{1}(n)(\delta^{\ell}+M_{n}\eta_{1}(n))$ (B.40)

with probability approaching 1 by Assumption 6. Therefore, we can obtain

$\displaystyle I_{4}\leq$
$\displaystyle\bigg{|}(\bm{v}^{(1)T}A(\widehat{\bm{\beta}}^{(2)},\delta)-\bm{v}^{*T}A(\bm{\beta}^{*},\delta))(\bm{v}^{(1)T}A(\widehat{\bm{\beta}}^{(2)},\delta)+\bm{v}^{*T}A(\bm{\beta}^{*},\delta))\bigg{|}$
(B.41) $\displaystyle\leq$
$\displaystyle\bigg{|}\bm{v}^{(1)T}(A(\widehat{\bm{\beta}}^{(2)},\delta)-A(\bm{\beta}^{*},\delta))+(\bm{v}^{(1)T}-\bm{v}^{*T})A(\bm{\beta}^{*},\delta)\bigg{|}\bigg{|}\bm{v}^{(1)T}A(\widehat{\bm{\beta}}^{(2)},\delta)+\bm{v}^{*T}A(\bm{\beta}^{*},\delta)\bigg{|}$
$\displaystyle=$
$\displaystyle||\bm{v}^{*}||_{1}^{2}\delta^{2\ell}\mathcal{O}_{\mathbb{P}}(M_{n}\eta_{1}(n)(1\vee
M_{n}\eta_{1}(n)/\delta^{\ell})+\eta_{2}(n)).$

Combining all bounds for $I_{1},\dotso,I_{4}$, we know that uniformly over all
$\delta\in\Delta$

$\displaystyle\bigg{|}\bigg{(}\widehat{\bm{v}}^{(1)T}\frac{1}{|\mathcal{N}_{1}|}\sum_{i\in\mathcal{N}_{1}}A_{i}(\widehat{\bm{\beta}}^{(2)},\delta)\bigg{)}^{2}-SB(\delta)\bigg{|}$
(B.42) $\displaystyle\lesssim$
$\displaystyle||\bm{v}^{*}||_{1}^{2}\delta^{2\ell}\bigg{(}\frac{\log(nd)}{nb^{2\ell+1}}+\sqrt{\frac{\log(nd)}{nb^{2\ell+1}}}+(\delta\vee
b)^{r}+M_{n}\eta_{1}(n)(1\vee\frac{M_{n}\eta_{1}(n)}{\delta^{\ell}})+\eta_{2}(n)\bigg{)}.$

This completes the proof.

###  B.6 Proof of Theorem 3

By Lemma 2 and Lemma 3, uniformly over all $\delta\in\Delta$

$\displaystyle|\widehat{M}(\delta)-M(\delta)|$ (B.43) $\displaystyle=$
$\displaystyle||\bm{v}^{*}||_{1}^{2}\frac{1}{n\delta}\mathcal{O}_{\mathbb{P}}\bigg{(}\eta_{2}(n)+\sqrt{\frac{\log(nd)}{n\delta}}+M_{n}\eta_{1}(n)\bigg{)}$
$\displaystyle~{}~{}~{}~{}~{}~{}+||\bm{v}^{*}||_{1}^{2}\delta^{2\ell}\mathcal{O}_{\mathbb{P}}\bigg{(}\frac{\log(nd)}{nb^{2\ell+1}}+\sqrt{\frac{\log(nd)}{nb^{2\ell+1}}}+(\delta\vee
b)^{r}+M_{n}\eta_{1}(n)(1\vee\frac{M_{n}\eta_{1}(n)}{\delta^{\ell}})+\eta_{2}(n)\bigg{)}$
$\displaystyle\lesssim$
$\displaystyle\frac{1}{n\delta}\psi_{1}(n)+\delta^{2\ell}\psi_{2}(n).$

Recall that $M(\delta)=\frac{1}{n}V(\delta)+\frac{n-1}{n}SB(\delta),$ where

$\displaystyle
V(\delta)=\frac{1}{\delta}(\bm{v}^{*T}\bm{\Sigma}^{*}\bm{v}^{*})(1+o(1))~{}~{}~{}~{}~{}\text{and}~{}~{}{SB(\delta)}=\delta^{2\ell}(\bm{v}^{*T}\bm{b}^{*})^{2}(1+o(1)).$
(B.44)

This implies that
$\delta^{*}=(\frac{1}{n}\frac{\bm{v}^{*T}\bm{\Sigma}^{*}\bm{v}^{*}}{2\ell(\bm{v}^{*T}\bm{b}^{*})^{2}})^{1/(2\ell+1)}(1+o(1))$.
By Assumption 5 we know $M(\delta)\gtrsim\frac{1}{n\delta}+\delta^{2\ell}$,
which implies

$\frac{|\widehat{M}(\delta)-M(\delta)|}{M(\delta)}\lesssim\psi_{1}(n)\vee\psi_{2}(n)=o_{\mathbb{P}}(1),$

uniformly over all $\delta\in\Delta$. That means for any $\epsilon>0$,
uniformly over all $\delta\in\Delta$, the event
$(1-\epsilon)M(\delta)\leq\widehat{M}(\delta)\leq(1+\epsilon)M(\delta)$ holds
with probability tending to 1. Under this event, we have that

$\displaystyle M(\widehat{\delta})$
$\displaystyle=\widehat{M}(\widehat{\delta})+[M(\widehat{\delta})-\widehat{M}(\widehat{\delta})]$
(B.45)
$\displaystyle\leq\widehat{M}(\delta^{*})+[M(\widehat{\delta})-\widehat{M}(\widehat{\delta})]$
$\displaystyle\leq(1+\epsilon)M(\delta^{*})+\epsilon M(\widehat{\delta}),$

where the first inequality holds because $\widehat{\delta}$ is the minimizer
of $\widehat{M}(\delta)$ and the second inequality follows from
$\widehat{M}(\delta^{*})\leq(1+\epsilon)M(\delta^{*})$ and
$\widehat{M}(\widehat{\delta})\geq(1-\epsilon)M(\widehat{\delta})$. As a
result, we obtain
$M(\widehat{\delta})\leq\frac{1+\epsilon}{1-\epsilon}M(\delta^{*})$, from
which we can claim that $\widehat{\delta}\leq\frac{3}{2}\delta^{*}$. To see
this, let us consider the complement case
$\widehat{\delta}>\frac{3}{2}\delta^{*}$. When it holds, we know from (B.44)
that

$M(\widehat{\delta})\geq(\frac{3}{2})^{2\ell}SB(\delta^{*})(1+o(1))=(\frac{3}{2})^{2\ell}\frac{1}{2\ell+1}M(\delta^{*})(1+o(1))\geq
1.01M(\delta^{*})(1+o(1)),$

where the second step follows from
$M(\delta^{*})=(2\ell+1)SB(\delta^{*})(1+o(1))$ and last step holds by the
condition $\ell\geq 2$. However, the above inequality contradicts with
$M(\widehat{\delta})\leq\frac{1+\epsilon}{1-\epsilon}M(\delta^{*})$ for some
sufficiently small $\epsilon$. This justifies the statement
$\widehat{\delta}\leq\frac{3}{2}\delta^{*}$. Following a similar argument, we
can show that $\widehat{\delta}\geq\frac{1}{2}\delta^{*}$. Thus, with
probability tending to 1, we have
$|\widehat{\delta}-\delta^{*}|/\delta^{*}\leq 1/2$.

By definition, we have $M^{\prime}(\delta^{*})=0$ and
$\widehat{M}^{\prime}(\widehat{\delta})=0$. This implies that

$\displaystyle\widehat{M}^{\prime}(\widehat{\delta})-M^{\prime}(\widehat{\delta})=-M^{\prime\prime}(\widetilde{\delta})(\widehat{\delta}-\delta^{*}),$
(B.46)

for some intermediate value $\widetilde{\delta}$, which gives

$\displaystyle(\widehat{\delta}-\delta^{*})=\frac{\widehat{M}^{\prime}(\widehat{\delta})-M^{\prime}(\widehat{\delta})}{-M^{\prime\prime}(\widetilde{\delta})}.$
(B.47)

For any ${\delta}\asymp\delta^{*}$, after some algebra, we can show that

$\displaystyle M^{\prime\prime}(\delta)\asymp\frac{1}{\delta^{3}n}.$ (B.48)

Similar to the proof of (B.43), we can show that
$(\widehat{M}^{\prime}-M^{\prime})(\delta)$ has the same order as
$(\widehat{M}-M)(\delta)$, except for an additional factor $\frac{1}{\delta}$.
Since $|\widehat{\delta}-\delta^{*}|/\delta^{*}\leq 1/2$ holds in probability,
it implies $\widehat{\delta}$ and $\widetilde{\delta}$ are of the same order
of $\delta^{*}$. Thus, (B.47) implies

$\displaystyle\frac{\widehat{\delta}-\delta^{*}}{\delta^{*}}\lesssim$
$\displaystyle\psi_{1}(n)\vee\psi_{2}(n).$ (B.49)

This completes the proof.

###  B.7 Proof of Theorem 4

The proof is very similar to the proof of Theorem 1 so we only give a sketch
here. Following the same derivation of Lemma 1, we can show that

$\displaystyle(n\delta)^{1/2}\frac{S_{\delta}^{L}(0,\bm{\gamma}^{*},\bm{\omega}^{*}_{L})-\delta^{\ell}\mu^{*}_{L}}{\sigma^{*}_{L}}\stackrel{{\scriptstyle
d}}{{\rightarrow}}N(0,1).$ (B.50)

Next, following the derivation of Theorem 1, we can show that
$I_{1}=(n\delta)^{1/2}|\bm{v}^{*T}_{L}\bm{C}(\nabla
R_{\delta}^{n_{(1)}}(\widehat{\bm{\beta}}_{0}^{(2)})-\nabla
R_{\delta}^{n_{(1)}}(\bm{\beta}^{*}))|$ and
$I_{2}=(n\delta)^{1/2}|(\widehat{\bm{v}}_{L}^{(1)}-\bm{v}_{L}^{*})\bm{C}\nabla
R_{\delta}^{n_{(1)}}(\widehat{\bm{\beta}}_{0}^{(2)})|$ are
$o_{\mathbb{P}}(1)$, which further implies that
$(n\delta)^{1/2}|S_{\delta}^{L(1)}(0,\widehat{\bm{\gamma}}^{(2)},\widehat{\bm{\omega}}^{(1)}_{L})-S_{\delta}^{L}(0,\bm{\gamma}^{*},\bm{\omega}^{*}_{L})|=o_{\mathbb{P}}(1)$.
This shall hold similarly for
$S_{\delta}^{L(2)}(0,\widehat{\bm{\gamma}}^{(1)},\widehat{\bm{\omega}}^{(2)}_{L})$.
Invoking Slutsky’s theorem completes the proof.

##  Appendix C Proofs of additional technical lemmas

######  Lemma 4.

Under Assumptions 1 \- 4, for any fixed $\bm{\beta}$, we have

$\displaystyle||\nabla R_{\delta}(\bm{\beta})-\nabla
R(\bm{\beta})||_{\infty}\lesssim\delta^{\ell},$ (C.1)
$\displaystyle||\nabla^{2}R_{\delta}(\bm{\beta})-\nabla
R(\bm{\beta})||_{\max}\lesssim\delta^{\ell-1}.$

###### Proof of Lemma 4.

We focus on proving the result for
$||\nabla^{2}R_{\delta}(\bm{\beta})-\nabla^{2}R(\bm{\beta})||_{\max}$ and the
proof for the other one is very similar. By definition

$\displaystyle||\nabla^{2}R_{\delta}(\bm{\beta})-\nabla^{2}R(\bm{\beta})||_{\max}$
(C.2) $\displaystyle=$
$\displaystyle\max_{j,k}|\sum_{y}w(y)\int_{\bm{z}}yz_{j}z_{k}(\int\frac{-1}{\delta^{2}}K^{\prime}(\frac{x-\bm{\beta}^{T}\bm{z}}{\delta})f(x|y,\bm{z})dx-f^{\prime}(\bm{\beta}^{T}\bm{z}|y,\bm{z}))f(y,\bm{z})d\bm{z}|$
$\displaystyle=$
$\displaystyle\max_{j,k}|\sum_{y}w(y)\int_{\bm{z}}yz_{j}z_{k}(\int\frac{-1}{\delta}K^{\prime}(u)f(u\delta+\bm{\beta}^{T}\bm{z}|y,\bm{z})du-f^{\prime}(\bm{\beta}^{T}\bm{z}|y,\bm{z}))f(y,\bm{z})d\bm{z}|$
$\displaystyle=$
$\displaystyle\max_{j,k}|\sum_{y}w(y)\int_{\bm{z}}yz_{j}z_{k}\int
K(u)(f^{\prime}(u\delta+\bm{\beta}^{T}\bm{z}|y,\bm{z})-f^{\prime}(\bm{\beta}^{T}\bm{z}|y,\bm{z}))duf(y,\bm{z})d\bm{z}|$
$\displaystyle=$
$\displaystyle\max_{j,k}|\sum_{y}w(y)\int_{\bm{z}}yz_{j}z_{k}\int
K(u)\frac{(u\delta)^{\ell-1}}{(\ell-1)!}f^{(\ell)}(\tau
u\delta+\bm{\beta}^{T}\bm{z}|y,\bm{z})duf(y,\bm{z})d\bm{z}|$
$\displaystyle\lesssim$ $\displaystyle\mathcal{O}(\delta^{\ell-1}),$

where the first equality is by definition, the second equality follows from a
change of variable, the third equality follows from an integration by parts,
the last equality follows from Assumptions 1, 2 and the last inequality
follows from Assumptions 3, 4. The proof is complete. ∎

######  Lemma 5.

Under Assumptions 1 \- 4, for any fixed $\bm{\beta}$, we have with probability
greater than $1-\mathcal{O}(d^{-1})$

$\displaystyle||\nabla R^{n}_{\delta}(\bm{\beta})-\nabla
R_{\delta}(\bm{\beta})||_{\infty}\lesssim\sqrt{\frac{\log(d)}{n\delta}},$
(C.3)
$\displaystyle||\nabla^{2}R^{n}_{\delta}(\bm{\beta})-\nabla^{2}R_{\delta}(\bm{\beta})||_{\max}\lesssim\sqrt{\frac{\log(d)}{n\delta^{3}}}.$

###### Proof of Lemma 5.

We focus on proving the result for
$||\nabla^{2}R^{n}_{\delta}(\bm{\beta})-\nabla^{2}R_{\delta}(\bm{\beta})||_{\max}$
and the proof for the other one is very similar. Denote

$\displaystyle T_{ijk}=$
$\displaystyle(\nabla^{2}\bar{R}^{i}_{\delta}(\bm{\beta})-\nabla^{2}R_{\delta}(\bm{\beta}))_{jk}$
(C.4) $\displaystyle=$
$\displaystyle-w(y_{i})y_{i}\frac{z_{ij}z_{ik}}{\delta^{2}}K^{\prime}(\frac{x_{i}-\bm{\beta}^{T}\bm{z}_{i}}{\delta})-(\nabla^{2}R_{\delta}(\bm{\beta}))_{jk}.$

By definition, we know $\mathbb{E}[T_{ijk}]=0$. Meanwhile for $Var[T_{ijk}]$,
we know

$\displaystyle\mathbb{E}[(\nabla^{2}\bar{R}^{i}_{\delta}(\bm{\beta}))_{jk}^{2}]=$
$\displaystyle\sum_{y}\int
w(y)^{2}\frac{z_{j}^{2}z_{k}^{2}}{\delta^{4}}K^{\prime
2}(\frac{x-\bm{\beta}^{T}\bm{z}}{\delta})f(x|y,\bm{z})dxf(y,\bm{z})d\bm{z}$
(C.5) $\displaystyle=$ $\displaystyle\sum_{y}\int
w(y)^{2}\frac{z_{j}^{2}z_{k}^{2}}{\delta^{3}}K^{\prime
2}(u)f(u\delta+\bm{\beta}^{T}\bm{z}|y,\bm{z})duf(y,\bm{z})d\bm{z}$
$\displaystyle=$ $\displaystyle\mathcal{O}(\frac{1}{\delta^{3}}),$

and with a similar derivation
$(\nabla^{2}R_{\delta}(\bm{\beta}))^{2}_{jk}=\mathcal{O}(\frac{1}{\delta^{2}})$.
This shows that $Var[T_{ijk}]=\mathcal{O}(\frac{1}{\delta^{3}})$. Since
$|T_{ijk}|\lesssim\frac{M_{n}^{2}}{\delta^{2}}$, by Bernstein inequality we
can show that with probability greater than $1-\mathcal{O}(d^{-1})$

$\displaystyle||\nabla^{2}R^{n}_{\delta}(\bm{\beta})-\nabla^{2}R_{\delta}(\bm{\beta})||_{\max}\lesssim\sqrt{\frac{\log(d)}{n\delta^{3}}}.$
(C.6)

This completes the proof. ∎

######  Lemma 6.

Under Assumptions 1 \- 4, it holds that

$\displaystyle||\nabla R(\widehat{\bm{\beta}})-\nabla
R(\bm{\beta}^{*})||_{\infty}\lesssim
M_{n}||\widehat{\bm{\beta}}-\bm{\beta}^{*}||_{1},$ (C.7)
$\displaystyle||\nabla^{2}R(\widehat{\bm{\beta}})-\nabla^{2}R(\bm{\beta}^{*})||_{\max}\lesssim
M_{n}||\widehat{\bm{\beta}}-\bm{\beta}^{*}||_{1}.$

###### Proof of Lemma 6.

Here we prove the second inequality and the first one should follow similarly.
By definition,

$\displaystyle||\nabla^{2}R(\widehat{\bm{\beta}})-\nabla^{2}R(\bm{\beta}^{*})||_{\max}$
(C.8) $\displaystyle=$