Variable Screening via Conditional Martingale

Difference Divergence

Lei Fang1, Qingcong Yuan2, Xiangrong Yin1 and Chenglong Ye1

1Dr. Bing Zhang Department of Statistics, University of Kentucky, 2Sanofi

> Abstract: Variable screening has been a useful research area that deals with
> ultrahigh-dimensional data. When there exist both marginally and jointly
> dependent predictors to the response, existing methods such as conditional
> screening or iterative screening often suffer from instability against the
> selection of the conditional set or the computational burden, respectively.
> In this article, we propose a new independence measure, named conditional
> martingale difference divergence ($\textrm{CMD}_{\mathcal{H}}$), that can be
> treated as either a conditional or a marginal independence measure. Under
> regularity conditions, we show that the sure screening property of
> $\textrm{CMD}_{\mathcal{H}}$ holds for both marginally and jointly active
> variables. Based on this measure, we propose a kernel-based model-free
> variable screening method, which is efficient, flexible, and stable against
> high correlation among predictors and heterogeneity of the response. In
> addition, we provide a data-driven method to select the conditional set. In
> simulations and real data applications, we demonstrate the superior
> performance of the proposed method.
>
> Key words and phrases: Independence measure; Reproducing kernel Hilbert
> space; Sure screening property; Conditional screening; Dimension reduction.

##  1  Introduction

Variable screening has been a research area that deals with ultrahigh-
dimensional data, where high-dimensional methods may fail due to the curse of
dimensionality, as Fan et al., (2009) suggested. Fan and Lv, (2008)’s seminal
work suggests to screen the ultrahigh-dimensional data before conducting
variable selection. They proposed a sure independent screening (SIS) method
for linear models to screen out inactive variables based on Pearson
correlation. After that, variable screening receives more attention since it
only requires that the selected set of variables covers the set of active
variables, which is referred to as the sure screening property (Fan and Lv,,
2008). Screening methods with this property suffer less from instability (Yu,,
2013) that is seen in many variable selection methods.

Various model-based screening methods have been developed. For linear
regression models, screening methods have been proposed based on different
measures, including marginal Pearson correlation (Fan and Lv,, 2008), forward
regression (Wang,, 2009), marginal empirical likelihood ratio (Chang et al.,,
2013), and Kendall’s rank correlation (Li et al., 2012a, ). For linear
quantile regression, screening methods based on quantile partial correlation
(Ma et al.,, 2017) and conditional quantile correlation (Zhang and Zhou,,
2018) have been proposed to handle heterogeneous data. In the context of
generalized linear models, screening methods based on the maximum marginal
likelihood or its estimate (Fan and Song,, 2010), the sparsity-restricted
maximum likelihood estimator (Xu and Chen,, 2014), and Kolmogorov-Smirnov
statistic (Mai and Zou,, 2013) have also been proposed. Other screening
methods include model settings such as linear regression models with
interactions (Hao and Zhang,, 2014; Fan et al.,, 2016; Kong et al.,, 2017),
Cox models (Zhao and Li,, 2012), varying coefficient models (Cheng et al.,,
2014; Song et al.,, 2014), and additive models (Fan et al.,, 2011).

An alternative approach is the model-free screening method, which recently has
gained popularity due to its less stringent assumptions. Zhu et al., (2011)
proposed a sure independence ranking and screening approach (SIRS) for index
models. Li et al., 2012b  proposed a screening method (DC-SIS) based on
distance correlation, which can be applied to grouped variables. Cui et al.,
(2015) proposed a model-free screening method (MV-SIS) based on empirical
conditional distribution function for discriminant analysis. Shao and Zhang,
(2014) proposed the use of martingale difference correlation (MDC), which can
be applied to mean and quantile screening. Mai and Zou, (2015) proposed the
fused Kolmogorov filter that works with different types of response variables
and high covariate correlation. Nandy et al., (2021) proposed a screening
method based on covariate information number (CIN) motivated by Fisher
information. Han, (2019) developed a screening framework from the perspective
of loss functions and proposed a screening method based on conditional
strictly convex losses. Based on ball correlation, Pan et al., (2019) proposed
a generic screening method for biomedical discovery.

As pointed out by Li et al., 2012b ; Shao and Zhang, (2014); Tian and Feng,
(2021), screening methods based on marginal measures (e.g., the marginal
correlation between the response and each predictor) will possibly miss the
marginally but not jointly independent predictors. Two types of approaches
have been developed to handle this issue. One approach is conditional
screening, a screening procedure based on a given set of variables. For
example, Barut et al., (2016) proposed conditional sure independence screening
(CSIS) for assessing the conditional (on a given conditional set) contribution
of a predictor to the response in generalized linear models. Based on
conditional distance correlation, Wen et al., (2018) proposed a method that
adjusts for confounding variables. Tong et al., (2022) proposed a conditional
independence measure and its corresponding screening method (CIS) with false
discovery rate (FDR) control, which also works for heavy-tailed
predictors/responses. Another approach is screening via iterative procedures.
For example, the aforementioned forward regression iteratively selects the
variables. Xu and Chen, (2014) considered a method based on sparse MLE, where
the algorithm iteratively updates the coefficients in the link function. Zhou
et al., (2020) proposed a model-free forward screening method that iteratively
updates the conditional set and is robust to outliers. Tian and Feng, (2021)
proposed an iterative variable screening method based on random subspace
ensembles (RaSE) with a theoretical guarantee for iterative screening
procedures.

However, two challenges remain. Model-based iterative methods (e.g., iterative
SIS) may rely on a specific variable selection method, which makes the
procedure less stable. On the other hand, the conditional screening method
requires prior knowledge of the conditional set and its performance becomes
unstable if an unreasonable conditional set is selected. It motivates us to
develop a stable model-free screening method that identifies both marginally
and jointly dependent variables to the response. We propose a kernel-based
measure that captures both conditional and marginal mean-independent
relationships. In particular, via Bochner’s theorem (Wendland,, 2004), we
transform the problem of choosing weights, a key element in our independence
measure, to the problem of choosing kernels and their bandwidths in
reproducing kernel Hilbert space (RKHS). This flexible kernel-based fashion
allows our method to perform well in various settings, as illustrated in the
synthetic and real data analysis.

The advantages of our method are as follows. First, we propose a kernel-based
independence measure ($\textrm{CMD}_{\mathcal{H}}$) that is able to
characterize both conditional and marginal mean independence. Thus, we propose
a $\textrm{CMD}_{\mathcal{H}}$-based screening method that can detect both
marginally and jointly dependent/active variables. Second, the proposed model-
free screening method is stable against outliers, data heterogeneity, and high
covariate correlation. Third, we show the sure screening property holds for
screening both marginally and jointly dependent variables under mild
regularity conditions. We also suggest selecting a data-driven conditional set
for conducting conditional screening when no prior information is available.

The rest of the article is organized as follows. Section 2 introduces the
proposed independence measure and its theoretical properties. In Section 3, we
propose a model-free variable screening procedure and present its sure
screening property. The simulation results and two real data examples are
reported in Section 4, followed by the conclusion in Section 5. Additional
theorems are presented in the appendix. Auxiliary simulation results and
technical proofs are included in the supplementary material.

##  2 General methodology

##### Notations.

Throughout the article, we use upper case (e.g., $V$) to denote a random
variable and use bold font to denote a random vector (e.g., $\boldsymbol{U}$).
We use $(\boldsymbol{U}_{1}^{\prime},\boldsymbol{U}_{2}^{\prime},V^{\prime})$
and
$(\boldsymbol{U}_{1}^{\prime\prime},\boldsymbol{U}_{2}^{\prime\prime},V^{\prime\prime})$
to denote $i.i.d$ copies of $(\boldsymbol{U}_{1},\boldsymbol{U}_{2},V)$. For a
complex function $f(\boldsymbol{s}):\mathbb{R}^{q}\rightarrow\mathbb{C}^{p}$,
denote its RKHS norm as
$||f(\boldsymbol{s})||^{2}_{\mathcal{H}_{k}}=\int_{\mathbb{R}^{q}}|f(\boldsymbol{s})|^{2}w(\boldsymbol{t})d\boldsymbol{t}$,
where $w(\boldsymbol{t})$ corresponds to the kernel function $k$ in RKHS.

For a sequence $\\{t_{ij}\\}$ with double indices $i,j=1,...,n$, we define

$t_{ij}^{*}=t_{ij}-\bar{t}_{i.}-\bar{t}_{.j}+\bar{t}_{..},$ (2.1)

where $\bar{t}_{.j}=\frac{1}{(n-2)}\sum_{i=1}^{n}t_{ij}$,
$\bar{t}_{i.}=\frac{1}{(n-2)}\sum_{j=1}^{n}t_{ij}$, and
$\bar{t}_{..}=\frac{1}{(n-1)(n-2)}\sum_{i=1}^{n}\sum_{j=1}^{n}t_{ij}$. Denote
$\langle\boldsymbol{a},\boldsymbol{b}\rangle$ as the inner product of any two
vectors $\boldsymbol{a},\boldsymbol{b}$ of the same dimension.

###  2.1 Conditional Martingale Difference Divergence

##### A motivating example.

Assume two random variables $X_{1}$ and $X_{2}$ are independent with
$E(X_{2})=0$, and let $Y=X_{1}X_{2}$. We have
$E(Y|X_{1})-E(Y)=E(X_{2})\cdot(X_{1}-E(X_{1}))=0$ as long as $E(X_{2})=0$. The
condition “$E(X_{2})=0$” is mild as we can standardize predictors in the
dataset to have mean 0 in practice. This example indicates that mean
independence measures based on the relationship $E(Y|X_{1})-E(Y)$ will
misleadingly suggest that $X_{1}$ is independent from mean of $Y$. It was
pointed out by Li et al., 2012b  and Shao and Zhang, (2014) that the marginal
measures/methods such as DC and MDC will possibly miss the variables that only
jointly contribute to the response variable. It motivates us to consider the
following:

$E(Y|X_{1},X_{2})-E(Y|X_{2}),$

which equals to $X_{2}\cdot(X_{1}-E(X_{1}))$ in this example. More generally,
we consider the equality
$E(V|\boldsymbol{U}_{1},\boldsymbol{U}_{2})=E(V|\boldsymbol{U}_{1})$, i.e.,
the response variable $V$ and covariate vector $\boldsymbol{U}_{2}$ given the
covariate vector $\boldsymbol{U}_{1}$. If in addition, $\boldsymbol{U}_{1}$
and $\boldsymbol{U}_{2}$ are independent, we have
$E(V|\boldsymbol{U}_{1},\boldsymbol{U}_{2})=E(V|\boldsymbol{U}_{1})$ if and
only if
$E(Ve^{i\langle\boldsymbol{t}_{1},\boldsymbol{U}_{1}\rangle}|\boldsymbol{U}_{2})=E(Ve^{i\langle\boldsymbol{t}_{1},\boldsymbol{U}_{1}\rangle})$
for any $\boldsymbol{t}_{1}$ (the proof is presented in supplementary material
S2 (b)). This motivates us to propose the following independence measure in
Definition 1, which can be treated as either a conditional or a marginal (see
Remark 2) mean independence measure.

######  Definition 1.

Given a random vector $\boldsymbol{U}_{1}\in\mathbb{R}^{p}$, the conditional
martingale difference divergence of a random variable $V$ and a random vector
$\boldsymbol{U}_{2}\in\mathbb{R}^{q}$ is defined as

$\displaystyle\textnormal{$\textrm{CMD}_{\mathcal{H}}$}^{2}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})$
$\displaystyle=\iint|E(Ve^{i(\langle\boldsymbol{t}_{1},\boldsymbol{U}_{1}\rangle+\langle\boldsymbol{t}_{2},\boldsymbol{U}_{2}\rangle)})-E(Ve^{i\langle\boldsymbol{t}_{1},\boldsymbol{U}_{1}\rangle})E(e^{i\langle\boldsymbol{t}_{2},\boldsymbol{U}_{2}\rangle})|^{2}w_{1}(\boldsymbol{t}_{1})w_{2}(\boldsymbol{t}_{2})d\boldsymbol{t}_{1}d\boldsymbol{t}_{2},$

where $w_{1}(\boldsymbol{t}_{1})$ and $w_{2}(\boldsymbol{t}_{2})$ are weight
functions.

The measure $\textrm{CMD}_{\mathcal{H}}$ depends on two ingredients: a mean
independence measure of a random vector and a random variable, and an
adjusting method of the effect of a third vector. It provides valuable
information on the conditional contribution of $\boldsymbol{U}_{2}$ to the
mean of $V$ given $\boldsymbol{U}_{1}$.

######  Remark 1.

We choose the weight function $w_{1}(\boldsymbol{t}_{1})$ and
$w_{2}(\boldsymbol{t}_{2})$ to be integrable, relaxing the strong assumption
of the boundedness of $\boldsymbol{U}_{1}$ and $\boldsymbol{U}_{2}$ as in the
literature. The choice of an integrable weight function makes the proposed
independence measure more flexible. In particular, we can rewrite
$\textnormal{$\textrm{CMD}_{\mathcal{H}}$}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})$
as a functional of kernel functions in RKHS. See more details in Theorem 1(a)
and Remark 3.

We now define a scale-invariant version of the proposed measure.

######  Definition 2.

Let $k_{1}$ and $k_{2}$ be the two kernel functions that correspond to the
weight function $w_{1}(\boldsymbol{t}_{1})$ and $w_{2}(\boldsymbol{t}_{2})$,
as illustrated in Theorem 1(a). We define the conditional martingale
difference correlation

$\textnormal{$\textrm{CMC}_{\mathcal{H}}$}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})=\begin{cases}\dfrac{\textnormal{$\textrm{CMD}_{\mathcal{H}}$}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})}{\sqrt{\textrm{v}(k_{2},\boldsymbol{U}_{2}){v}(k_{1_{V}},\boldsymbol{U}_{1})}}&\text{if
$\textrm{v}(k_{2},\boldsymbol{U}_{2}){v}(k_{1_{V}},\boldsymbol{U}_{1})>0$ }\\\
0&\text{ otherwise,}\end{cases}$

where
$\textrm{v}(k,\boldsymbol{U}):=E{[}k^{2}(\boldsymbol{U},\boldsymbol{U}^{\prime}){]}+E^{2}{[}k(\boldsymbol{U},\boldsymbol{U}^{\prime}){]}-2E[k(\boldsymbol{U},\boldsymbol{U}^{{}^{\prime}})\cdot
k(\boldsymbol{U},\boldsymbol{U}^{{}^{\prime\prime}})]$ and
$k_{V}(\boldsymbol{U},\boldsymbol{U}^{\prime}):=VV^{{}^{\prime}}k(\boldsymbol{U},\boldsymbol{U}^{\prime})$
for any kernel function $k$.

######  Remark 2.

When $\boldsymbol{U}_{1}$ contains no useful information
($\boldsymbol{U}_{1}=\emptyset$, $\boldsymbol{U}_{1}\equiv\mathbf{c}$, or
$\boldsymbol{U}_{1}$ is independent from $(V,\boldsymbol{U}_{2})$), the
definition of $\textrm{CMD}_{\mathcal{H}}$ reduces to a marginal mean
independence measure. That is,

$\displaystyle\textnormal{$\textrm{MD}_{\mathcal{H}}$}^{2}(V,\boldsymbol{U}_{2}):=\int|E(Ve^{i\langle\boldsymbol{t}_{2},\boldsymbol{U}_{2}\rangle})-E(V)E(e^{i\langle\boldsymbol{t}_{2},\boldsymbol{U}_{2}\rangle})|^{2}w(\boldsymbol{t}_{2})d\boldsymbol{t}_{2}.$

Note that $\textnormal{$\textrm{MD}_{\mathcal{H}}$}(V,\boldsymbol{U}_{2})$ is
a generalized version of MDD (Shao and Zhang,, 2014) by kernerlizing the
$L_{2}$ distance of $\boldsymbol{U}_{2}$ and its i.i.d. copy
$\boldsymbol{U}^{\prime}_{2}$. The standardized version
$\textrm{MC}_{\mathcal{H}}$ is defined similarly as
$\textrm{CMC}_{\mathcal{H}}$.

As will be seen in Theorem 1, the definition of the
$\textrm{CMC}_{\mathcal{H}}$ is more convenient for variable screening purpose
since it takes values in $[0,1]$. More detailed discussion of
$\textrm{CMC}_{\mathcal{H}}$ is included in Section 3. Now we show the
theoretical properties of the proposed conditional independence measure.

######  Theorem 1.

Assume $E(V^{2})<\infty$, we have the following properties:

  1. (a).

We can rewrite
$\textnormal{$\textrm{CMD}_{\mathcal{H}}$}^{2}{(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})}$
as

$\displaystyle\textnormal{$\textrm{CMD}_{\mathcal{H}}$}^{2}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})$
$\displaystyle=E(VV^{\prime}k_{1}(\boldsymbol{U}_{1},\boldsymbol{U}_{1}^{{}^{\prime}})k_{2}(\boldsymbol{U}_{2},\boldsymbol{U}_{2}^{{}^{\prime}}))+E(VV^{\prime}k_{1}(\boldsymbol{U}_{1},\boldsymbol{U}_{1}^{{}^{\prime}}))E(k_{2}(\boldsymbol{U}_{2},\boldsymbol{U}_{2}^{{}^{\prime}}))$
$\displaystyle-2E(VV^{\prime}k_{1}(\boldsymbol{U}_{1},\boldsymbol{U}_{1}^{{}^{\prime}})k_{2}(\boldsymbol{U}_{2},\boldsymbol{U}_{2}^{{}^{\prime\prime}})),$

where $k_{1}$ and $k_{2}$ are RHKS kernel functions determined by
$w_{1}(\boldsymbol{t}_{1})$ and $w_{2}(\boldsymbol{t}_{2})$ defined in
Definition 1, respectively.

  2. (b).

$0\leq\textnormal{$\textrm{CMC}_{\mathcal{H}}$}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})\leq
1,$ and
$\textnormal{$\textrm{CMC}_{\mathcal{H}}$}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})=0\Leftrightarrow
E(V|\boldsymbol{U}_{1},\boldsymbol{U}_{2})=E(V|\boldsymbol{U}_{1})$ a.s. if
$\boldsymbol{U}_{1}\perp\boldsymbol{U}_{2}$.

  3. (c).

Given two constants $d\in\mathbb{R}$ and $e\in\mathbb{R}$,
$\textnormal{$\textrm{CMC}_{\mathcal{H}}$}(a+bV,\boldsymbol{c}+d\boldsymbol{U}_{2}|e\boldsymbol{U}_{1})=\textnormal{$\textrm{CMC}_{\mathcal{H}}$}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})$
for any scalars $a,b\in\mathbb{R}$ and $\boldsymbol{c}\in\mathbb{R}^{q}$. If
the kernels $k_{1}$ and $k_{2}$ in (a) are scale-invariant, the above equality
holds for any scalars $d$ and $e$ as well.

  4. (d).

If the random variables $U,V\in\mathbb{R}$ are independent, then

$\textnormal{$\textrm{MC}_{\mathcal{H}}$}^{2}(VU,U)=\frac{E^{2}(V)}{Var(V)+E^{2}(V)+E^{2}(U)\frac{Var(V)}{Var({U})}}\textnormal{$\textrm{MC}_{\mathcal{H}}$}^{2}(U,U).$

Furthermore, if $E(V)=0$, then
$\textnormal{$\textrm{MC}_{\mathcal{H}}$}(VU,U)=0$.

######  Remark 3.

If we take non-integrable weight functions $w_{1}(\boldsymbol{t}_{1})$ and
$w_{2}(\boldsymbol{t}_{2})$ in Definition 1, then $k_{1}$ and $k_{2}$ in
Theorem 1(a) may not be translation-invariant kernels in RKHS (e.g., the
Euclidean distance function). See dCov (Székely et al.,, 2007) for an example
that adopts a non-integrable weight in its definition.

######  Remark 4.

Property (b) shows the equivalence between the conditional mean independence
and $\textrm{CMC}_{\mathcal{H}}$ being 0, which suggests
$\textrm{CMC}_{\mathcal{H}}$ is a suitable tool for conducting variable
screening. Note that the independence of $\boldsymbol{U}_{1}$ and
$\boldsymbol{U}_{2}$ in Property (b) is to ease the proof. Indeed, if we
define a new independence measure
$\textnormal{$\textrm{CMD}_{\mathcal{H}}$}_{,new}^{2}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})=\iint|E(Ve^{i(\langle\boldsymbol{t}_{1},\boldsymbol{U}_{1}\rangle+\langle\boldsymbol{t}_{2},\boldsymbol{U}_{2}\rangle)}|\boldsymbol{U}_{1})-E(Ve^{i\langle\boldsymbol{t}_{1},\boldsymbol{U}_{1}\rangle}|\boldsymbol{U}_{1})E(e^{i\langle\boldsymbol{t}_{2},\boldsymbol{U}_{2}\rangle}|\boldsymbol{U}_{1})|^{2}w_{1}(\boldsymbol{t}_{1})w_{2}(\boldsymbol{t}_{2})d\boldsymbol{t}_{1}d\boldsymbol{t}_{2},$
then

$\textnormal{$\textrm{CMD}_{\mathcal{H}}$}_{,new}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})=0\textit{
a.s.}\Leftrightarrow
E(V|\boldsymbol{U}_{1},\boldsymbol{U}_{2})=E(V|\boldsymbol{U}_{1})\textit{
a.s..}$

This removes the independence condition of $\boldsymbol{U}_{1}$ and
$\boldsymbol{U}_{2}$. Then we need to replace our $U$-statistics estimator
with the conditional $U$-statistics to estimate the new independence measure.
Note that this new measure $\textnormal{$\textrm{CMD}_{\mathcal{H}}$}_{,new}$
is a function of the random vector $\boldsymbol{U}_{1}$. Such a conditional
measure and its associated screening method is left as future work of
interest. In this article, we stick to our original proposed measure
$\textrm{CMD}_{\mathcal{H}}$. In the simulations, as we see, even
$\boldsymbol{U}_{1}$ and $\boldsymbol{U}_{2}$ are not independent (e.g., high
variable correlations $\rho=0.5,0.8,0.9$ in Example 1, and nonlinearly
associated predictors in Example 4), our variable screening method still
performs well, or even outperforms other methods in almost all the simulation
settings.

######  Remark 5.

Property (c) shows that the proposed $\textrm{CMC}_{\mathcal{H}}$ is scale-
invariant. Property (d) directly shows the deficiency of marginal-type mean
independence measure $(\textnormal{$\textrm{MC}_{\mathcal{H}}$})$ in
interaction screening. Thus we propose the variable screening procedure based
on $\textrm{CMC}_{\mathcal{H}}$.

###  2.2 Empirical Estimators and Asymptotic Properties

Based on property (a) in Theorem 1, we construct the $U$-statistic to estimate
$\textrm{CMC}_{\mathcal{H}}$.

######  Definition 3.

Let $(\boldsymbol{U}_{1i},\boldsymbol{U}_{2i},V_{i})_{i=1}^{n}$ be $i.i.d.$
observations of $(\boldsymbol{U}_{1},\boldsymbol{U}_{2},V)$. Denote
$a_{ij}=V_{i}V_{j}k_{1}(\boldsymbol{U}_{1i},\boldsymbol{U}_{1j})$ and
$b_{ij}=k_{2}(\boldsymbol{U}_{2i},\boldsymbol{U}_{2j})$ for $i,j=1,...,n$.
Define the corresponding $a_{ij}^{*}$ and $b_{ij}^{*}$ as in Equation (2.1).
The $U$-statistic estimator of $\textrm{CMD}_{\mathcal{H}}$ is

$\textnormal{$\widehat{\textrm{CMD}}_{\mathcal{H}}$}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})=\frac{1}{n(n-3)}\sum_{1\leq
i\neq j\leq n}a_{ij}^{*}b_{ij}^{*},$

and the corresponding estimator of $\textrm{CMC}_{\mathcal{H}}$ is:

$\textnormal{$\widehat{\textrm{CMC}}_{\mathcal{H}}$}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})=\dfrac{\sum_{1\leq
i\neq j\leq n}a_{ij}^{*}b_{ij}^{*}}{\sqrt{\sum_{1\leq i\neq j\leq
n}a_{ij}^{*2}\sum_{1\leq i\neq j\leq n}b_{ij}^{*2}}}.$

######  Remark 6.

Compared to the adoption of $V$-statistic estimator, we choose the
$U$-statistic because it is unbiased and less computationally expensive.

We now show the strong consistency of the proposed estimators.

######  Theorem 2.

(Consistency) If E $(V^{2})<\infty$, then

$lim_{n\rightarrow\infty}\textnormal{$\widehat{\textrm{CMD}}_{\mathcal{H}}$}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})=\textnormal{$\textrm{CMD}_{\mathcal{H}}$}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})\textrm{
a.s.,}$ (2.2)

and

$lim_{n\rightarrow\infty}\textnormal{$\widehat{\textrm{CMC}}_{\mathcal{H}}$}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})=\textnormal{$\textrm{CMC}_{\mathcal{H}}$}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})\textrm{
a.s.}.$ (2.3)

In the next theorem, we derive the asymptotic distribution for
$\textnormal{$\textrm{CMD}_{\mathcal{H}}$}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})$.
Denote the following functions:
$\textsl{g}_{\boldsymbol{U}_{2}}(\boldsymbol{t}_{2}):=E(e^{i\langle\boldsymbol{t}_{2},\boldsymbol{U}_{2}\rangle})$,
$\textsl{g}_{V,\boldsymbol{U}_{1}}(\boldsymbol{t}_{1}):=E(Ve^{i\langle\boldsymbol{t}_{1},\boldsymbol{\boldsymbol{U}_{1}}\rangle})$,
and
$F(\boldsymbol{t}_{1},\boldsymbol{t}_{2}):=E(V^{2}e^{i\langle\boldsymbol{t}_{1},\boldsymbol{U}_{1}\rangle}e^{i\langle\boldsymbol{t}_{2},\boldsymbol{U}_{2}\rangle})$.
Define the covariance function
$\textnormal{cov}_{\Gamma}((\boldsymbol{t}{}_{1},\boldsymbol{t}{}_{2}),(\boldsymbol{t}^{\prime}_{1},\boldsymbol{t}^{\prime}_{2})):=F(\boldsymbol{t}_{1}-\boldsymbol{t}^{\prime}_{1},\boldsymbol{t}_{2}-\boldsymbol{t}^{\prime}_{2})+(F(\boldsymbol{t}_{1}-\boldsymbol{t}^{\prime}_{1},0)+\textsl{g}_{V,\boldsymbol{U}_{1}}(\boldsymbol{t}{}_{1})\overline{\textsl{g}_{V,\boldsymbol{U}_{1}}(\boldsymbol{t}^{\prime}_{1})})\\{\textsl{g}_{\boldsymbol{U}_{2}}(\boldsymbol{t}{}_{2})\overline{\textsl{g}_{\boldsymbol{U}_{2}}(\boldsymbol{t}^{\prime}_{2})}-\textsl{g}_{\boldsymbol{U}_{2}}(\boldsymbol{t}_{2}-\boldsymbol{t}^{\prime}_{2})\\}$$-F(\boldsymbol{t}_{1}-\boldsymbol{t}^{\prime}_{1},\boldsymbol{t}_{2})\overline{\textsl{g}_{\boldsymbol{U}_{2}}(\boldsymbol{t}^{\prime}_{2})}-F(\boldsymbol{t}_{1}-\boldsymbol{t}^{\prime}_{1},-\boldsymbol{t}^{\prime}_{2})\textsl{g}_{\boldsymbol{U}_{2}}(\boldsymbol{t}{}_{2})$.

######  Theorem 3.

Assume $E(V^{2})<\infty$, we have the following:

  1. a.

If
$\textnormal{$\textrm{CMD}_{\mathcal{H}}$}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})=\textnormal{0}$,
then

$n\textnormal{$\widehat{\textrm{CMD}}_{\mathcal{H}}$}^{2}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})\overset{d}{\rightarrow}||\Gamma(s)||_{\mathcal{H}_{k}}^{2}$
(2.4)

as $n\rightarrow\infty$, where $\Gamma(\cdot)$ is a complex-valued zero-mean
Gaussian random process with covariance function
$\textnormal{cov}_{\Gamma}((\boldsymbol{t}{}_{1},\boldsymbol{t}{}_{2}),(\boldsymbol{t}^{\prime}_{1},\boldsymbol{t}^{\prime}_{2}))$.

  2. b.

If
$\textnormal{CMD${}_{\mathbf{\mathcal{H}}}$}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})=0$
and $E(V^{2}|\boldsymbol{U}_{2})=E(V^{2})$, then

$n\textnormal{$\widehat{\textrm{CMD}}_{\mathcal{H}}$}^{2}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})/S_{n}\overset{d}{\rightarrow}\sum_{j=1}^{\infty}\lambda_{j}Z_{j}$

as $n\rightarrow\infty$, where
$S_{n}=(\frac{1}{n}\sum_{i}V_{i}^{2}-\frac{1}{n(n-1)}\sum_{i\neq
j}a_{ij})(1-\frac{1}{n(n-1)}\sum_{i\neq j}b_{ij})$,
$Z_{j}\overset{i.i.d.}{\sim}\chi_{1}^{2}$, and
$\\{\lambda_{j}\\}_{j=1}^{\infty}$ are nonnegative constants such that
$E(\sum_{j=1}^{\infty}\lambda_{j}Z_{j})=1$.

  3. c.

If
$\textnormal{$\textrm{CMD}_{\mathcal{H}}$}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})>0$,
then
$n\cdot\textnormal{$\widehat{\textrm{CMD}}_{\mathcal{H}}$}^{2}(V,\boldsymbol{U}_{2}|\boldsymbol{U}_{1})/S_{n}\overset{p}{\rightarrow}\infty$
as $n\rightarrow\infty$.

The properties stated in the theorems of this section motivate us to propose
variable screening algorithm based on $\textrm{CMC}_{\mathcal{H}}$ and its
estimate $\widehat{\textrm{CMC}}_{\mathcal{H}}$.

##  3 $\textrm{CMC}_{\mathcal{H}}$-based Variable Screening

In this section, we show the sure screening property of
$\textrm{CMC}_{\mathcal{H}}$ in Section 3.1. In Section 3.2, we introduce a
variable screening algorithm named $\textrm{S-CMC}_{\mathcal{H}}$ to
accommodate the dependence among predictors. The algortihm works reasonably
well when the data suffer from outliers, high correlation, and heterogeneity
as seen in the numerical studies.

###  3.1 Sure Screening Property

Without loss of generality, let $Y$ be a univariate continuous response
variable and $\boldsymbol{X}=(X_{1},...,X_{p})^{T}$ be the predictor vector.
Denote the sample as $(X_{1k},...,X_{pk},Y_{k})_{k=1}^{n}$, where $p\gg n$.
For any index set $S\subseteq\\{1,...,p\\}$, denote
$\boldsymbol{X}_{S}:=\\{\boldsymbol{X}_{j}:j\in S\\}$. Given a conditional set
$\boldsymbol{X}_{S}$ with cardinality $d_{1}$, we define

$\mathcal{D}_{S}=\\{j:\mathbb{E}(Y|(\boldsymbol{X}_{S},X_{j}))\ \text{depends
on}\ X_{j}\\}$

as the index set of dependent/active predictors conditional on
$\boldsymbol{X}_{S}$, and

$\mathcal{I}_{S}=\\{j:\mathbb{E}(Y|(\boldsymbol{X}_{S},X_{j}))\ \text{is
independent from}\ X_{j}\\}$

as the index set of independent/inactive predictors conditional on
$\boldsymbol{X}_{S}$. Note that $\mathcal{D}_{S}$ is a subset of
$\mathcal{D}:=\\{j:\mathbb{E}(Y|\boldsymbol{X})\ \text{depends on}\ X_{j}\\}$,
the set of all dependent predictors. Suppressing $S$, we denote
$\omega_{j}=\textrm{CMC}_{\mathcal{H}}^{2}(Y,X_{j}|\boldsymbol{X}_{S})$ as the
dependence score of $X_{j}$ given $\boldsymbol{X}_{S}$. Let
$\hat{\omega}_{j}=\widehat{\textrm{CMC}_{\mathcal{H}}}^{2}(Y,X_{j}|\boldsymbol{X}_{S})$
be the estimator of $\omega_{j}$ and

$\hat{\mathcal{D}}_{S}=\\{j:\hat{\omega}_{j}\geq cn^{-\kappa},\ \text{for}\
j\in S^{c}\\}$

be the set of selected variables after screening. Before stating the sure
screening property, we assume the following conditions.

(A1)

There exists a constant $s_{0}>0$ such that $E(\exp(sY^{2}))<\infty$ for all
$0<s\leq 2s_{0}$.

(A2)

For any given $\boldsymbol{X}_{S}$, $\min_{j\in\mathcal{D}_{S}}\omega_{j}\geq
2cn^{-\kappa}$ for some constant $c>0$ and $0\leq\kappa<1/2$.

Condition $(A1)$ puts constraint on the tail distribution of the response
variable and Condition (A2) requires that the conditional active/dependent
variables and inactive/independent variables are well separated.

######  Theorem 4.

Under Condition $(A1)$, for any $0<\gamma<1/2-\kappa$, there exist positive
constants $c_{1}$ and $c_{2}$ such that

$P\Bigg{\\{}\underset{1\leq j\leq
p-d_{1}}{\max}|\widehat{\omega}_{j}-\omega_{j}|\geq cn^{-\kappa}\Bigg{\\}}\leq
O((p-d_{1})[\exp(-c_{1}n^{1-2(\kappa+\gamma)})+n\exp(-c_{2}n^{\gamma})]).$
(3.5)

If Conditions $(A2)$ also holds, we have

$P(\mathcal{D}_{S}\subset\hat{\mathcal{D}}_{S})\geq
1-O(s_{n}[\exp(-c_{1}n^{1-2(\kappa+\gamma)})+n\exp(-c_{2}n^{\gamma})])$ (3.6)

for any conditional set $\boldsymbol{X}_{S}$, where $s_{n}$ is the cardinality
of $\mathcal{D}_{S}$. In particular, let
$\delta={\min}_{j\in\mathcal{D}_{S}}\omega_{j}-\textnormal{max}_{j\in\mathcal{I}_{S}}\omega_{j}$,
we have the ranking consistency:

$P\Bigg{\\{}\underset{j\in\mathcal{I}_{S}}{\max}\widehat{\omega}_{j}<\underset{j\in\mathcal{D}_{S}}{\min}\widehat{\omega}_{j}\Bigg{\\}}\geq
1-2O((p-d_{1})[\exp(-c^{\prime}_{1}\delta^{2}n^{1-2\gamma})+n\exp(-c_{2}n^{\gamma})])$
(3.7)

for a positive constant $c_{1}^{\prime}$.

######  Remark 7.

The above theorem shows the sure screening property holds for any given
conditional set $\boldsymbol{X}_{S}$. Define
$\mathcal{M}:=\\{j:\mathbb{E}(Y|X_{j})\ \text{depends on}\ X_{j}\\}$ as the
set of marginally dependent/active predictors. For the special case where the
conditional set $\boldsymbol{X}_{S}=\boldsymbol{X}$, we have
$\mathcal{D}_{S}=\mathcal{D}$, which is the common sure screening property in
the literature. In Appendix B.1, we also show that, similarly to Theorem 4,
the sure screening property holds when the conditional set is empty. In that
case, $\textrm{CMC}_{\mathcal{H}}$ reduces to $\textrm{MC}_{\mathcal{H}}$, and
the sure screening property holds for selecting the set $\mathcal{D}$ as well
as $\mathcal{M}$. We discuss more details of selecting the conditional set
$\boldsymbol{X}_{S}$ in Section 3.2.

######  Remark 8.

The error terms $\exp(-c_{1}n^{1-2(\kappa+\gamma)})$ and
$n\exp(-c_{2}n^{\gamma})$ in (3.5) comes from estimating the three terms in
CMDH as in Theorem 1(a). In the proof, take the first term
$E(VV^{\prime}k_{1}(\boldsymbol{U}_{1},\boldsymbol{U}_{1}^{{}^{\prime}})k_{2}(\boldsymbol{U}_{2},\boldsymbol{U}_{2}^{{}^{\prime}})):=E(h)$
for example, we decompose it into a bounded term $E[hI(h<M)]$ plus an
unbounded term $E[hI(h>M)]$ for some large enough $M>0$. Similar decomposition
is done to the other two terms. Setting $M=n^{\gamma}$ for some
$0<\gamma<1/2-\kappa$, we obtain the two error terms
$\exp(-c_{1}n^{1-2(\kappa+\gamma)})$ and $n\exp(-c_{2}n^{\gamma})$ for
estimating the sum of the bounded terms and that of the unbounded terms,
respectively. The role of the parameter $\gamma$ is a trade-off of estimating
the bounded and unbounded terms. By setting $\gamma=\frac{1-2\kappa}{3}$, we
achieve a balance and obtain the optimal convergence rate. As mentioned in
Shao and Zhang, (2014), their bound (3.5) can be further improved by assuming
a stronger moment condition on $Y$, i.e., $E(\exp(sY^{4}))<\infty$ for all
$s\in(0,2s_{0}]$. Their improved bound is the same as our bound. It is also
worth mentioning that we do not impose moment conditions on the variable X as
in Shao and Zhang, (2014). The reason why our method enjoys a better rate
under a weaker condition is that our proposed measure CMCH only computes the
RKHS kernel functions of X (as shown in Theorem 1 (a)). Such kernels are
bounded, which frees us from assuming additional moment conditions on X. In
contrast, the Martingale Difference Correlation requires to calculate the
Euclidean distance (Theorem 1 (1) in Shao and Zhang, (2014)), which is
unbounded. Thus, they require the stronger assumptions.

###  3.2 A Variable Screening Algorithm $\textrm{S-CMC}_{\mathcal{H}}$

The sure screening property holds for any conditional set. In practice, if the
conditional set is not given, we use the top $d_{1}$ predictors suggested by
$\textrm{MC}_{\mathcal{H}}$, which also enjoys the sure screening property, as
shown in Theorems 9 and 10 in the appendix. We propose to use the following
variable screening algorithm $\textrm{S-CMC}_{\mathcal{H}}$ as stated in
Algorithm 1.

Input:

The conditional set $\boldsymbol{X}_{S}$ (optional) and its cardinality
$d_{1}$ (optional), the number of variables $d_{2}$ to select (optional), and
the data $\\{(y_{i},\boldsymbol{x}_{i})\\}_{i=1}^{n}$.

  1. 1.

If $\boldsymbol{X}_{S}$ is not given, calculate
$\widehat{\textrm{MC}_{\mathcal{H}}}(X_{i}):=\widehat{\textrm{MC}_{\mathcal{H}}}(Y,X_{i})$
for each $i=1,...,p$. Let the conditional set $\boldsymbol{X}_{S}$ be the set
of the top $d_{1}$(if not given, $d_{1}=\left\lfloor\sqrt{n/\log
n}\right\rfloor$ where $\lfloor\cdot\rfloor$ is the floor function) predictors
with the largest $\widehat{\textrm{MC}_{\mathcal{H}}}(X_{i})$.

  2. 2.

For each $i\in\\{1,...,p-d_{1}\\}$, calculate
$\widehat{\textrm{CMC}}_{\mathcal{H}}(X_{c_{i}}):=\widehat{\textrm{CMC}}_{\mathcal{H}}(Y,X_{c_{i}}^{\perp}|\boldsymbol{X}_{S})$,
where each $c_{i}$ is from ${S^{c}}$ and
$X_{c_{i}}^{\perp}=X_{c_{i}}-P_{\boldsymbol{X}_{S}}X_{c_{i}}$ with
$P_{\boldsymbol{X}_{S}}$ being the projection matrix onto the column space of
$\boldsymbol{X}_{S}$.

  3. 3.

Calculate the score
$A_{i}:=\max(\frac{\textrm{MC}_{\mathcal{H}}(X_{c_{i}})}{\underset{1\leq i\leq
p-d_{1}}{\max}\textrm{MC}_{\mathcal{H}}(X_{c_{i}})},\frac{\textrm{CMC}_{\mathcal{H}}(X_{c_{i}})}{\underset{1\leq
i\leq p-d_{1}}{\max}\textrm{CMC}_{\mathcal{H}}(X_{c_{i}})})$ for each
$i=1,...,p-d_{1}$. Keep the top $d_{2}-d_{1}$ predictors with the largest
scores. If $d_{2}$ is not given, $d_{2}=\left\lfloor n/\log n\right\rfloor$.

Output:

The index set of the $d_{2}$ selected variables
$\\{i_{1},...,i_{d_{2}}\\}\subseteq\\{1,...,p\\}$ (the $d_{1}$ variables
selected in Step 1 plus the $d_{2}-d_{1}$ variables selected in Step 3).

Algorithm 1 The procedure of the $\textrm{S-CMC}_{\mathcal{H}}$ for variable
screening.

The R code of $\textrm{S-CMC}_{\mathcal{H}}$ is available in the supplement
file.

######  Remark 9.

The parameters $d_{1}$ and $d_{2}$ are predefined values. In general, larger
$d_{1}$ will lead to worse performance if the set $\boldsymbol{X}_{S}$
contains larger proportion of inactive variables. It will not affect the
theoretical performance of $\textrm{CMC}_{\mathcal{H}}$. However,
computationally, the estimation of the expected kernel function of long
vectors inside $\textrm{CMC}_{\mathcal{H}}$ becomes less reliable if the
sample size is limited. Thus we recommend to choose small $d_{1}$(e.g.,
$d_{1}=\left\lfloor\sqrt{n/\log n}\right\rfloor$) in practice.

######  Remark 10.

The adoption of $X_{c_{i}}^{\perp}$ is to handle the scenario when the
independence assumption of $\boldsymbol{X}_{S}$ and $X_{i}$ is violated. Note
that this only removes the linear dependence. See more discussion in Remark 4.

######  Remark 11.

In the last step, an alternative way is to use a linear combination
$w_{1}\frac{\textnormal{$\textrm{MC}_{\mathcal{H}}$}(X_{c_{i}})}{\underset{1\leq
i\leq
p-d_{1}}{\max}\textnormal{$\textrm{MC}_{\mathcal{H}}$}(X_{c_{i}})}+w_{2}\frac{\textnormal{$\textrm{CMC}_{\mathcal{H}}$}(X_{c_{i}})}{\underset{1\leq
i\leq p-d_{1}}{\max}\textnormal{$\textrm{CMC}_{\mathcal{H}}$}(X_{c_{i}})}$ to
rank all the predictors. The weights can be adaptively chosen driven by the
data, which we leave as a potential future work.

###  3.3 Extension to Quantile Screening

In this section, we extend our method to the quantile screening setting. For a
univariate random response $Y$, denote ${Y}_{{\tau}}=\tau-\boldsymbol{1}(Y\leq
q_{\tau})$ as its binary version with $\tau\in(0,1)$, where $q_{\tau}$ is the
$\tau$-th quantile of the distribution of $Y$. Given $i.i.d.$ observations
$\\{y_{k}\\}_{k=1}^{n}$ of $Y$, denote
$\hat{{Y}}_{{\tau}}=\tau-\boldsymbol{1}(Y\leq\hat{q}_{\tau})$ as the estimate
of ${Y}_{{\tau}}$, where $\hat{q}_{\tau}$ is the sample $\tau$-th quantile. So
for each observation $y_{k}$, we denote
$y_{k_{\tau}}=\tau-\boldsymbol{1}(y_{k}\leq\hat{q}_{\tau})$. Let
$\omega_{j}({Y}_{{\tau}})=\textrm{CMC}_{\mathcal{H}}^{2}({Y}_{{\tau}},X_{j}|\boldsymbol{X}_{S})$
and
$\widehat{\omega}_{j}(\hat{Y}_{\tau})=\widehat{\textrm{CMC}}_{\mathcal{H}}^{2}(\hat{{Y}}_{{\tau}},X_{j}|\boldsymbol{X}_{S})$.
Similarly, we denote
$\mathcal{D}_{q_{\tau}}=\\{j:\mathbb{E}(Y_{\tau}|(\boldsymbol{X}_{S},X_{j}))\
\text{depends on}\ X_{j}\\}$ as the quantile active predictors conditional on
$\boldsymbol{X}_{S}$, and denote
$\mathcal{\widehat{D}}_{q_{\tau}}=\\{j:\widehat{\omega}_{j}(\hat{{Y}}_{{\tau}})\geq
cn^{-\kappa},$ for $j\in S^{c}\\}$ as the selected variables.

Next we show the sure screening property for the quantile version of
$\textrm{CMC}_{\mathcal{H}}$.

######  Theorem 5.

Under condition $(C1)$ in the appendix, for any $0<\gamma<1/2-\kappa$ and
$\kappa\in(0,1/2)$, there exist positive constants $c_{1},c_{2}$ such that for
any $c>0$,

$P\Bigg{\\{}\underset{1\leq j\leq
p-d_{1}}{\max}|\widehat{\omega}_{j}(\hat{Y}_{\tau})-\omega_{j}({Y}_{{\tau}})|\geq
cn^{-\kappa}\Bigg{\\}}\leq
O((p-d_{1})[\exp\\{-c_{1}n^{1-2(\kappa+\gamma)}\\}+n\exp(-c_{2}n^{\gamma})]).$
(3.8)

If Condition $(C2)$ in the appendix holds in addition, we have

$P(\mathcal{D}_{q_{\tau}}\subseteq\hat{\mathcal{D}}_{q_{\tau}})\geq
1-O(\widetilde{s}_{n}[\exp\\{-c_{1}n^{1-2(\kappa+\gamma)}\\}+n\exp(-c_{2}n^{\gamma})]),$
(3.9)

where $\widetilde{s_{n}}$ is the cardinality of $\mathcal{D}_{q_{\tau}}$.

##  4 Numerical Studies

In this section, we evaluate the finite-sample performance of the proposed
method $\textrm{CMC}_{\mathcal{H}}$.

##### The choice of kernel functions

Denote the translation-invariant Gaussian kernel as

$\displaystyle
K(\boldsymbol{x}_{1},\boldsymbol{x}_{2}):=\exp(-\frac{1}{h}(\boldsymbol{x}_{1}-\boldsymbol{x}_{2})^{T}(\boldsymbol{x}_{1}-\boldsymbol{x}_{2})),$
(4.10)

where $\boldsymbol{x}_{1},\boldsymbol{x}_{2}\in\mathbb{R}^{t},t\in\mathbb{N}$,
and $h$ is the bandwidth. For the proposed $\textrm{CMC}_{\mathcal{H}}$, we
use $K(\boldsymbol{x}_{1},\boldsymbol{x}_{2})$ for both $k_{1}$ and $k_{2}$.
In our simulations, the performance of $\textrm{CMC}_{\mathcal{H}}$ for
variable screening is robust against the bandwidths of $k_{1}$ and $k_{2}$. So
we set $h=2$ for both $k_{1}$ and $k_{2}$. For $\textrm{MC}_{\mathcal{H}}$, we
adopt Gaussian kernel and conduct a sensitivity analysis of the bandwidth, the
results of which are presented in the supplementary material S1.1. The
performance of $\textrm{MC}_{\mathcal{H}}$ is sensitive to the bandwidth $h$.
In particular, $\textrm{MC}_{\mathcal{H}}$ with smaller $h$ performs better
for selecting covariates that are linearly related to the response variable,
while larger $h$ is more suitable for selecting nonlinearly related
covariates. To select the conditional set $\boldsymbol{X}_{S}$ and avoid
cherry-picking, we first calculate the values of $\textrm{MC}_{\mathcal{H}}$
using two bandwidths $h=2\hat{\sigma}_{X_{i}}^{2}$ and
$h=6\hat{\sigma}_{X_{i}}^{2}$ for each predictor $X_{i}$, where
$\hat{\sigma}_{X_{i}}^{2}$ is sample variance of $X_{i}$ and $i=1,...,p$. Our
experience is that $h=2\hat{\sigma}_{X_{i}}^{2}$ and
$h=6\hat{\sigma}_{X_{i}}^{2}$ generally perform well for the simulations. Then
we take the maximum of the two $\textrm{MC}_{\mathcal{H}}$ values for each
predictor. One can also use Laplacian kernel and Cauchy kernel in practice.
However, in our examples, they yield similar performance to that of the
Gaussian kernel.

##### Criteria of evaluating variable screening performance

Following Li et al., 2012b , we consider three criteria for evaluating the
variable screening performance: 1) $\mathcal{S}_{q}$: the $(100q)$-th quantile
of the minimal model size required to contain all the active predictors. 2)
$\mathcal{P}_{i}$: the proportion that the predictor $X_{i}$ is selected, and
3) $\mathcal{P}_{\textnormal{all}}$: the proportion of all active predictors
being selected. Essentially, an $\mathcal{S}_{q}$ closer to the total number
of active predictors is preferred. The three criteria are connected in a way
that a smaller minimal model size suggests a larger $\mathcal{P}_{all}$ and a
lager $\mathcal{P}_{i}$ for each active variable.

##### Screening thresholds

We compare $\textrm{S-CMC}_{\mathcal{H}}$ with six variable screening methods,
including two marginal screening method (DCSIS2 in Kong et al., (2017) and
MDC), three conditional screening methods (CSIS in Barut et al., (2016),
CDCSIS in Wen et al., (2018) and CIS in Tong et al., (2022)), and one
iterative methods (RaSE1-eBIC in Tian and Feng, (2021). For each screening
method, we keep the top $d_{2}=\lfloor n/log(n)\rfloor$ variables, where $n$
is the sample size. We report the $\mathcal{P}_{i}$, $\mathcal{P}_{all}$ and
$\mathcal{S}_{0.5}$ values based on 100 repetitions for each example. Since
conditional screening methods require a pre-selected conditional set
$\boldsymbol{X}_{S}$, we either artificially set up the conditional set based
on the variables in the true model or select the top
$d_{1}=\lfloor\sqrt{n/log(n)}\rfloor$ variables suggested by
$\textrm{MC}_{\mathcal{H}}$. We also include a sensitivity analysis of using
different methods (e.g. SIS, LASSO and forward regression) to choose the
conditional set in the supplementary material S1.2.

###  4.1 Simulation

##### Example 1 (Marginally inactive but jointly active predictors).

Following the idea of Fan and Lv, (2008), we generate samples
$\\{Y_{i},\boldsymbol{X}_{i}\\}_{i=1}^{n}$ from the linear regression model

$Y=X_{1}+X_{2}+X_{3}+X_{4}+X_{5}-cX_{6}+\epsilon,$

where the coefficient $c$ is designed so that $\textrm{cov}(X_{6},Y)=0$. That
is, the predictor $X_{6}$ is marginally independent from the response $Y$. The
predictor vector $\boldsymbol{X}=(X_{1},...,X_{p})\sim
N(\boldsymbol{0},\boldsymbol{\Sigma}$), where $\boldsymbol{\Sigma}$ has
$(i,j)$-th entry $\sigma_{ij}=\rho^{I\\{i\neq j\\}}$. The error term
$\epsilon\sim N(0,1)$ and is independent from $\boldsymbol{X}$. We set the
sample size $n=200$ and the dimension $p=3000$. We consider three cases:
$(c,\rho)\in\\{(2.5,0.5),(4,0.8),(4.5,0.9)\\}$. Note that $X_{6}$ is dependent
with $Y$ if given one or more predictors from $\\{X_{1},...,X_{5}\\}$.

| $\mathcal{P}_{1}$ | $\mathcal{P}_{2}$ | $\mathcal{P}_{3}$ | $\mathcal{P}_{4}$ | $\mathcal{P}_{5}$ | $\mathcal{P}_{6}$ | $\mathcal{P}_{all}$ | $\mathcal{S}_{0.5}$
---|---|---|---|---|---|---|---|---
$c=2.5$, $\rho=0.5$
MDC | 0.89 | 0.94 | 0.93 | 0.91 | 0.90 | 0.00 | 0.00 | 3000.0
CSIS ($\boldsymbol{X}_{S_{1}}$) | 1.00 | 0.97 | 0.98 | 0.98 | 0.97 | 0.87 | 0.78 | 14.0
CSIS ($\boldsymbol{X}_{S_{2}}$) | 0.71 | 0.75 | 0.75 | 0.72 | 0.72 | 1.00 | 0.17 | 1580.5
CDC-SIS ($\boldsymbol{X}_{S_{1}}$) | 1.00 | 0.95 | 0.95 | 0.95 | 0.93 | 0.16 | 0.14 | 708.5
CDC-SIS ($\boldsymbol{X}_{S_{2}}$) | 0.73 | 0.76 | 0.74 | 0.72 | 0.72 | 1.00 | 0.16 | 1541.5
CIS($\boldsymbol{X}_{S_{1}}$) | 1.00 | 0.83 | 0.88 | 0.87 | 0.87 | 0.71 | 0.34 | 58.0
CIS($\boldsymbol{X}_{S_{2}}$) | 0.97 | 0.91 | 0.94 | 0.94 | 0.97 | 0.01 | 0.01 | 2997.5
$\textrm{S-CMC}_{\mathcal{H}}$ ($\boldsymbol{X}_{S_{1}}$) | 1.00 | 0.94 | 0.94 | 0.92 | 0.92 | 0.53 | 0.39 | 90.5
$\textrm{S-CMC}_{\mathcal{H}}$ ($\boldsymbol{X}_{S_{2}}$) | 0.91 | 0.95 | 0.94 | 0.93 | 0.92 | 1.00 | 0.72 | 17.5
RaSE1-eBIC | 0.99 | 1.00 | 0.99 | 0.99 | 1.00 | 1.00 | 0.97 | 6.0
$c=4$, $\rho=0.8$
MDC | 0.61 | 0.63 | 0.64 | 0.63 | 0.62 | 0.00 | 0.00 | 3000.0
CSIS ($\boldsymbol{X}_{S_{1}}$) | 1.00 | 0.48 | 0.44 | 0.40 | 0.44 | 1.00 | 0.10 | 320.0
CSIS ($\boldsymbol{X}_{S_{2}}$) | 0.39 | 0.34 | 0.36 | 0.41 | 0.33 | 1.00 | 0.00 | 2988.5
CDC-SIS ($\boldsymbol{X}_{S_{1}}$) | 1.00 | 0.56 | 0.53 | 0.54 | 0.54 | 0.93 | 0.06 | 264.0
CDC-SIS ($\boldsymbol{X}_{S_{2}}$) | 0.39 | 0.34 | 0.36 | 0.41 | 0.34 | 1.00 | 0.00 | 2623.5
CIS($\boldsymbol{X}_{S_{1}}$) | 1.00 | 0.39 | 0.31 | 0.34 | 0.42 | 1.00 | 0.03 | 546.0
CIS($\boldsymbol{X}_{S_{2}}$) | 0.70 | 0.59 | 0.67 | 0.64 | 0.69 | 1.00 | 0.09 | 268.0
$\textrm{S-CMC}_{\mathcal{H}}$ ($\boldsymbol{X}_{S_{1}}$) | 1.00 | 0.73 | 0.71 | 0.66 | 0.67 | 1.00 | 0.26 | 127.5
$\textrm{S-CMC}_{\mathcal{H}}$ ($\boldsymbol{X}_{S_{2}}$) | 0.61 | 0.73 | 0.72 | 0.66 | 0.67 | 1.00 | 0.16 | 156.5
RaSE1-eBIC | 0.80 | 0.92 | 0.76 | 0.79 | 0.87 | 1.00 | 0.35 | 2312.5
$c=4.5$, $\rho=0.9$
MDC | 0.48 | 0.40 | 0.39 | 0.48 | 0.40 | 0.00 | 0.00 | 3000.0
CSIS ($\boldsymbol{X}_{S_{1}}$) | 1.00 | 0.14 | 0.09 | 0.13 | 0.08 | 1.00 | 0.00 | 2339.0
CSIS ($\boldsymbol{X}_{S_{2}}$) | 0.28 | 0.24 | 0.23 | 0.24 | 0.20 | 1.00 | 0.00 | 2992.5
CDC-SIS ($\boldsymbol{X}_{S_{1}}$) | 1.00 | 0.29 | 0.28 | 0.29 | 0.28 | 0.98 | 0.01 | 875.0
CDC-SIS ($\boldsymbol{X}_{S_{2}}$) | 0.28 | 0.24 | 0.23 | 0.24 | 0.21 | 1.00 | 0.00 | 2875.0
CIS($\boldsymbol{X}_{S_{1}}$) | 1.00 | 0.11 | 0.08 | 0.07 | 0.16 | 1.00 | 0.00 | 2037.5
CIS($\boldsymbol{X}_{S_{2}}$) | 0.43 | 0.44 | 0.47 | 0.46 | 0.41 | 1.00 | 0.00 | 897.0
$\textrm{S-CMC}_{\mathcal{H}}$ ($\boldsymbol{X}_{S_{1}}$) | 1.00 | 0.48 | 0.47 | 0.51 | 0.45 | 1.00 | 0.07 | 269.5
$\textrm{S-CMC}_{\mathcal{H}}$ ($\boldsymbol{X}_{S_{2}}$) | 0.53 | 0.48 | 0.47 | 0.51 | 0.47 | 1.00 | 0.03 | 338.0
RaSE1-eBIC | 0.62 | 0.89 | 0.73 | 0.65 | 0.68 | 1.00 | 0.19 | 2316.0

Table 1: The $\mathcal{P}_{i}$, $\mathcal{P}_{all}$ and $\mathcal{S}_{0.5}$ in
Example 1. The conditional set is either $\boldsymbol{X}_{S_{1}}=\\{X_{1}\\}$
or $\boldsymbol{X}_{S_{2}}=\\{\textrm{the first $d_{1}=6$ predictors selected
by MDC${}_{\mathcal{H}}$}\\}$.

The simulation results are reported in Table 1. The marginal screening method
MDC fails to detect the marginally independent predictor $X_{6}$ in all cases,
while other methods identify $X_{6}$ as an active predictor. As the
correlation increases from $\rho=0.5$ to $\rho=0.9$, the selection proportion
$\mathcal{P}_{i}$ decreases for $i=1,2,...,5$. Consequently,
$\mathcal{P}_{all}$ decreases as the correlation increases. Note that if the
conditional set $\boldsymbol{X}_{S_{1}}=\\{X_{1}\\}$, the proportion of
selecting $X_{1}$ is set to 1. We set
$d_{2}=\left\lfloor{n/log(n)}\right\rfloor=37$ in this example. Note that the
minimal model size $\mathcal{S}_{0.5}$ is small only in the first case where
$\rho=0.5$, which explains the low values of $\mathcal{P}_{all}$ for all the
methods for the high correlation case ($\rho=0.8$ or $0.9$). The three
conditional methods: CIS, CSIS and CDC-SIS, fail to detect the active
variables $X_{2},X_{3},X_{4}$ and $X_{5}$ when the correlation increases to
$\rho=0.9$.

When the conditional set changes from an oracle set $\boldsymbol{X}_{S_{1}}$
to a data-dependent set $\boldsymbol{X}_{S_{2}}$, the performances of the
three conditional screening methods (CIS, CSIS and CDC-SIS) deteriorate. Our
proposed S-C$\textrm{MC}_{\mathcal{H}}$ instead shows robustness against the
conditional set. The performance of our method decreases due to the extra cost
of selecting $X_{1}$ when the conditional set becomes data-dependent.

In this example, from the perspective of the minimal model size
$\mathcal{S}_{0.5}$, the minimal model size for RaSE1-eBIC changes from $6$ to
$2316$ when $\rho$ increases from $0.5$ to $0.9$. In contrast, the stable
performance of S-$\textrm{CMC}_{\mathcal{H}}$ in $\mathcal{S}_{0.5}$ indicates
that our method is more robust against the correlation $\rho$. In terms of
$\mathcal{P}_{all}$ and $\mathcal{P}_{i}$, RaSE1-eBIC and
S-$\textrm{CMC}_{\mathcal{H}}$ are better than any other method. RaSE1-eBIC is
better than S-$\textrm{CMC}_{\mathcal{H}}$, which may be due to the good
performance of RaSE1-eBIC in the linear case.

##### Example 2 (Interaction terms).

We consider the following model with interaction terms:

$\displaystyle Y=$ $\displaystyle
X_{1}+X_{5}+X_{10}+X_{1}X_{15}+1.5X_{5}X_{20}+2X_{10}X_{25}+\epsilon.$

The predictor vector $\boldsymbol{X}=(X_{1},...,X_{p})\sim
N(\boldsymbol{0},\boldsymbol{\Sigma}$), where $\boldsymbol{\Sigma}$ has
$(i,j)$-th entry $\sigma_{ij}=\rho^{|i-j|}$. We consider two cases:
$\rho\in\\{0,0.9\\}$. The error term $\epsilon\sim N(0,1)$ and is independent
from $\boldsymbol{X}$. We set the sample size $n=200$ and the dimension
$p=3000$.

| $\mathcal{P}_{1}$ | $\mathcal{P}_{5}$ | $\mathcal{P}_{10}$ | $\mathcal{P}_{15}$ | $\mathcal{P}_{20}$ | $\mathcal{P}_{25}$ | $\mathcal{P}_{all}$ | $\mathcal{S}_{0.5}$
---|---|---|---|---|---|---|---|---
$\rho=0$
MDC | 0.94 | 0.96 | 0.90 | 0.01 | 0.00 | 0.03 | 0.00 | 2662.5
CSIS ($\boldsymbol{X}_{S_{1}}$) | 1.00 | 1.00 | 1.00 | 0.04 | 0.03 | 0.07 | 0.00 | 2452.0
CSIS($\boldsymbol{X}_{S_{2}}$) | 0.97 | 0.93 | 0.90 | 0.03 | 0.02 | 0.04 | 0.00 | 2147.0
CDC-SIS ($\boldsymbol{X}_{S_{1}}$) | 1.00 | 1.00 | 1.00 | 0.01 | 0.03 | 0.26 | 0.00 | 1639.0
CDC-SIS ($\boldsymbol{X}_{S_{2}}$) | 0.84 | 0.86 | 0.84 | 0.07 | 0.17 | 0.43 | 0.00 | 1529.0
DCSIS2 | 0.40 | 0.70 | 0.99 | 0.11 | 0.35 | 0.85 | 0.00 | 1372.5
CIS($\boldsymbol{X}_{S_{1}}$) | 1.00 | 1.00 | 1.00 | 0.02 | 0.12 | 0.45 | 0.00 | 566.0
CIS($\boldsymbol{X}_{S_{2}}$) | 0.99 | 0.98 | 0.96 | 0.01 | 0.01 | 0.02 | 0.00 | 2278.0
$\textrm{S-CMC}_{\mathcal{H}}$ ($\boldsymbol{X}_{S_{1}}$) | 1.00 | 1.00 | 1.00 | 0.18 | 0.62 | 0.99 | 0.09 | 196.5
$\textrm{S-CMC}_{\mathcal{H}}$($\boldsymbol{X}_{S_{2}}$) | 0.94 | 0.96 | 0.91 | 0.11 | 0.37 | 0.77 | 0.01 | 851.5
RaSE1-eBIC | 0.87 | 0.82 | 0.79 | 0.01 | 0.01 | 0.04 | 0.00 | 2295.5
$\rho=0.9$
MDC | 1.00 | 1.00 | 1.00 | 0.97 | 0.52 | 0.14 | 0.13 | 350.5
CSIS ($\boldsymbol{X}_{S_{1}}$) | 1.00 | 1.00 | 1.00 | 0.08 | 0.15 | 0.13 | 0.02 | 1887.0
CSIS ($\boldsymbol{X}_{S_{2}}$) | 0.78 | 0.97 | 0.74 | 0.08 | 0.12 | 0.13 | 0.02 | 1870.0
CDCSIS ($\boldsymbol{X}_{S_{1}}$) | 1.00 | 1.00 | 1.00 | 0.03 | 0.04 | 0.06 | 0.00 | 2024.5
CDCSIS($\boldsymbol{X}_{S_{2}}$) | 0.52 | 0.96 | 0.80 | 0.86 | 1.00 | 1.00 | 0.40 | 66.0
DCSIS2 | 0.93 | 1.00 | 1.00 | 0.87 | 0.91 | 0.81 | 0.62 | 30.5
CIS($\boldsymbol{X}_{S_{1}}$) | 1.00 | 1.00 | 1.00 | 0.28 | 0.93 | 1.00 | 0.27 | 62.0
CIS($\boldsymbol{X}_{S_{2}}$) | 0.78 | 0.99 | 0.83 | 0.43 | 0.92 | 0.85 | 0.2 | 88.5
$\textrm{S-CMC}_{\mathcal{H}}$($\boldsymbol{X}_{S_{1}}$) | 1.00 | 1.00 | 1.00 | 0.90 | 1.00 | 1.00 | 0.90 | 25.0
$\textrm{S-CMC}_{\mathcal{H}}$($\boldsymbol{X}_{S_{2}}$) | 1.00 | 1.00 | 1.00 | 0.90 | 1.00 | 0.99 | 0.89 | 26.0
RaSE1-eBIC | 0.41 | 0.42 | 0.37 | 0.02 | 0.02 | 0.01 | 0.01 | 2403.0

Table 2: The $\mathcal{P}_{i}$, $\mathcal{P}_{all}$ and $\mathcal{S}_{0.5}$ in
Example 2. We set the conditional set to be
$\boldsymbol{X}_{S_{1}}$=$\\{X_{1},X_{5},X_{10}\\}$ or
$\boldsymbol{X}_{S_{2}}=${the first $d_{1}=6$ predictors selected by
$\textrm{MC}_{\mathcal{H}}$}.

The results are reported in Table 2. The three variables
$X_{15},X_{20},X_{25}$ all jointly contribute to the mean of $Y$, but are
marginally independent of the mean of $Y$. It is difficult for the marginal
screening methods to detect these three terms. Note that $X_{25}$ has a larger
coefficient in its interaction term than that of $X_{20}$ or $X_{15}$. As the
signal/coefficient of the interaction term increases, its effect is easier to
be detected ($\mathcal{P}_{25}>\mathcal{P}_{20}>\mathcal{P}_{15}$) for all the
methods. RaSE1-eBIC fails in detecting those three variables in this example.
This is possibly due to the fact that RaSE1-eBIC targets additive models. In
comparison to the interaction screening method DCSIS2, the proposed method
S-C$\textrm{MC}_{\mathcal{H}}$ has a comparable performance in selecting
$X_{15}$, $X_{20}$ and $X_{25}$ for both $\rho=0$ and $\rho=0.9$. But
S-C$\textrm{MC}_{\mathcal{H}}$ performs better than DCSIS2 in selecting the
marginally active variables $X_{1}$, $X_{5}$ and $X_{10}$. What’s more,
S-$\textrm{CMC}_{\mathcal{H}}$ has the smallest $\mathcal{S}_{0.5}$ among all
the methods. Similar to Example 1, the performance of
S-$\textrm{CMC}_{\mathcal{H}}$ is stable against the conditional set. We also
did a sensitivity analysis against the conditional set. We consider three more
cases: the conditional set is selected by Lasso, SIS, and forward regression.
The results are reported in the supplementary material S1.2, demonstrating a
better and more stable performance of $\textrm{S-CMC}_{\mathcal{H}}$ compared
to other conditional methods (CSIS, CDC-SIS). Finally, we include a block
structure correlation setting where the correlation among active predictors
are 0.2 and 0.1 otherwise. We report the result in the supplementary material
S1.3, which clearly shows the advantage of S-$\textrm{CMC}_{\mathcal{H}}$.

##### Example 3 (Heteroscedasticity & Quantile screening).

In this example, we demonstrate that our screening method can help in
heteroskedastic model specification. We consider the following model:

$\displaystyle Y=$ $\displaystyle
X_{1}+X_{5}+X_{1}X_{10}+1.5X_{5}X_{15}+\epsilon\cdot\exp(X_{35}+X_{40}).$

The predictor vector $\boldsymbol{X}=(X_{1},...,X_{p})\sim
N(\boldsymbol{0},\boldsymbol{\Sigma}$), where $\boldsymbol{\Sigma}$ has
$(i,j)$-th entry $\sigma_{ij}=\rho^{|i-j|}$. We consider two cases:
$\rho\in\\{0,0.9\\}$. The error term $\epsilon\sim N(0,1)$ is independent from
$\boldsymbol{X}$. We set the sample size $n=400$ and the dimension $p=3000$.
For the purpose of quantile screening, we change the continuous response $Y$
to a binary response
$Y_{\textrm{$\tau$}}=\tau-\boldsymbol{1}(Y\leq\hat{q}_{\tau})$, where
$\hat{q}_{\tau}$ is the $\tau$-th sample quantile of the response. Then we
apply S-$\textrm{CMC}_{\mathcal{H}}$ and MDC on the data
$(Y_{\tau},\boldsymbol{X})$. Note that in fixed design, the population
quantile $q_{\tau}$ does not depend on $(X_{20},X_{25})$ if and only if
$\tau=0.5$. We consider two choices of the quantile: $\tau=0.5$ and
$\tau=0.75$.

The results are presented in Table 3 with conditional set selected by
$\textrm{MC}_{\mathcal{H}}$. In this example, we also evaluate the method
QaSIS (He et al.,, 2013), a quantile-adaptive model-free variable screening
method for heterogeneous data. Overall, our proposed method
S-C$\textrm{MC}_{\mathcal{H}}$ performs the best across all four combinations
of $(\rho,\tau)$. In particular, S-C$\textrm{MC}_{\mathcal{H}}$ is dominantly
better in ${\cal S}_{0.5}$ than any other method even under high correlation
setting. When there is no correlation among predictors ($\rho=0$), from the
perspective of quantile screening, MDC, CDC-SIS, QaSIS, and
S-C$\textrm{MC}_{\mathcal{H}}$ correctly differentiate the different roles of
$X_{35}$ and $X_{40}$ under two values of $\tau$. However, these three methods
(MDC, CDC-SIS, QaSIS) fail to identify $X_{15}$ (compared to $X_{10}$) when
$\tau=0.5$. In contrast, S-C$\textrm{MC}_{\mathcal{H}}$ has a better
performance in separating the active variables from the inactive variables.
When high correlation exists among predictors, all methods receive improved
performance and S-C$\textrm{MC}_{\mathcal{H}}$ still remains competitive. This
is because the marginal relationship between $X_{10}$ ($X_{15}$) and the
response $Y$ is strengthened by the high correlation among predictors. We also
include the scenario when conditional set is $\\{X_{1},X_{5},X_{10}\\}$ in the
supplementary material S1.4.

Method | $\tau$ | $\rho=0$ |  | $\rho=0.9$
---|---|---|---|---
$\mathcal{P}_{1}$ | $\mathcal{P}_{5}$ | $\mathcal{P}_{10}$ | $\mathcal{P}_{15}$ | $\mathcal{P}_{35}$ | $\mathcal{P}_{40}$ | $\mathcal{P}_{all}$ | $\mathcal{S}_{0.5}$ |  | $\mathcal{P}_{1}$ | $\mathcal{P}_{5}$ | $\mathcal{P}_{10}$ | $\mathcal{P}_{15}$ | $\mathcal{P}_{35}$ | $\mathcal{P}_{40}$ | $\mathcal{P}_{all}$ | $\mathcal{S}_{0.5}$
MDC | 0.5 | 1.00 | 1.00 | 0.06 | 0.03 | 0.04 | 0.02 | 0.00 | 1910.0 |  | 1.00 | 1.00 | 0.97 | 0.49 | 0.06 | 0.05 | 0.49 | 71.0
0.75 | 0.99 | 1.00 | 0.24 | 0.38 | 0.52 | 0.55 | 0.03 | 986.5 |  | 1.00 | 1.00 | 1.00 | 1.00 | 0.91 | 0.92 | 0.86 | 29.5
CSIS | 0.5 | 1.00 | 0.99 | 0.01 | 0.03 | 0.20 | 0.16 | 0.00 | 2025.5 |  | 1.00 | 1.00 | 0.19 | 0.07 | 0.31 | 0.24 | 0.00 | 1854.5
0.75 | 0.97 | 0.95 | 0.06 | 0.14 | 0.39 | 0.42 | 0.00 | 2207.0 |  | 0.86 | 0.96 | 0.39 | 0.11 | 0.31 | 0.30 | 0.00 | 2138.5
CDC-SIS | 0.5 | 1.00 | 0.99 | 0.13 | 0.28 | 0.23 | 0.30 | 0.04 | 757.0 |  | 1.00 | 1.00 | 0.28 | 0.94 | 1.00 | 1.00 | 0.27 | 172.5
0.75 | 0.96 | 0.97 | 0.13 | 0.51 | 0.42 | 0.47 | 0.04 | 1159.5 |  | 0.88 | 0.97 | 0.41 | 0.87 | 0.95 | 0.92 | 0.25 | 214.5
QaSIS | 0.5 | 1.00 | 1.00 | 0.16 | 0.16 | 0.28 | 0.27 | 0.03 | 1085.0 |  | 1.00 | 1.00 | 1.00 | 0.90 | 0.69 | 0.55 | 0.90 | 20.5
0.75 | 0.92 | 0.99 | 0.14 | 0.38 | 0.69 | 0.72 | 0.03 | 600.0 |  | 0.99 | 1.00 | 1.00 | 0.97 | 0.99 | 0.97 | 0.92 | 35.5
DCSIS2 | 0.5 | 0.18 | 0.41 | 0.04 | 0.08 | 0.93 | 0.90 | 0.00 | 1443.5 |  | 0.07 | 0.15 | 0.06 | 0.05 | 1.00 | 1.00 | 0.00 | 918.0
0.75 | 0.18 | 0.41 | 0.04 | 0.08 | 0.93 | 0.90 | 0.00 | 1443.5 |  | 0.07 | 0.15 | 0.06 | 0.05 | 1.00 | 1.00 | 0.00 | 918.0
CIS | 0.5 | 1.00 | 0.99 | 0.07 | 0.04 | 0.01 | 0.02 | 0.01 | 1982.5 |  | 1.00 | 1.00 | 0.31 | 0.89 | 0.95 | 0.98 | 0.29 | 121.5
0.75 | 1.00 | 0.99 | 0.06 | 0.14 | 0.21 | 0.3 | 0.00 | 2537.5 |  | 1.00 | 0.99 | 0.47 | 0.76 | 0.9 | 0.91 | 0.33 | 129.0
$\textrm{S-CMC}_{\mathcal{H}}$ | 0.5 | 1.00 | 1.00 | 0.18 | 0.80 | 0.05 | 0.04 | 0.13 | 240.0 |  | 1.00 | 1.00 | 0.97 | 0.93 | 0.34 | 0.38 | 0.90 | 16.0
0.75 | 0.99 | 1.00 | 0.33 | 0.67 | 0.49 | 0.58 | 0.05 | 574.5 |  | 1.00 | 1.00 | 1.00 | 0.95 | 0.97 | 0.93 | 0.86 | 28.0
RaSE1-eBIC | 0.5 | 0.49 | 0.31 | 0.00 | 0.00 | 0.03 | 0.01 | 0.00 | 1816.0 |  | 0.25 | 0.07 | 0.00 | 0.00 | 0.04 | 0.03 | 0.00 | 2425.5
0.75 | 0.49 | 0.31 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 2427.0 |  | 0.25 | 0.07 | 0.00 | 0.00 | 0.01 | 0.02 | 0.00 | 2427.0

Table 3: The $\mathcal{P}_{i}$, $\mathcal{P}_{all}$ and $\mathcal{S}_{0.5}$ in
Example 3 with $\boldsymbol{X}_{S}$ selected by $\textrm{MC}_{\mathcal{H}}$.

##### Example 4 (Nonlinear case).

We consider the following nonlinear case.

$Y=3I(X_{1}>0.5)X_{2}+3sin^{2}(2\pi
X_{1})X_{3}+3(X_{1}^{2}-1)X_{4}+\exp(X_{1})X_{5}+\epsilon,$

where we first generate $U_{1}$, $U_{2}$
$\overset{i.i.d}{\sim}\textrm{Unif}[0,1]$ and then let $X_{1}=(U_{1}+U_{2})/2$
and $X_{k}=(Z_{k}+2U_{1})/4$ for $k=2,3,...,p$. We consider two scenarios,
symmetric distribution $Z_{k}\stackrel{{\scriptstyle i.i.d}}{{\sim}}N(0,1)$
and asymmetric distribution $Z_{k}\stackrel{{\scriptstyle
i.i.d}}{{\sim}}\chi^{2}_{(1)}$. The error term $\epsilon$ is independently
drawn from standard normal distribution. We set the sample size $n=200$ and
the dimension $p=3000$. The results are reported in Table 4.

$Z_{k}$ $\stackrel{{\scriptstyle i.i.d}}{{\sim}}N(0,1)$ |  |  $Z_{k}$ $\stackrel{{\scriptstyle i.i.d}}{{\sim}}\chi^{2}_{(1)}$ |
---|---|---|---
| $\mathcal{P}_{1}$ | $\mathcal{P}_{2}$ | $\mathcal{P}_{3}$ | $\mathcal{P}_{4}$ | $\mathcal{P}_{5}$ | $\mathcal{P}_{all}$ | $\mathcal{S}_{0.5}$ |  | $\mathcal{P}_{1}$ | $\mathcal{P}_{2}$ | $\mathcal{P}_{3}$ | $\mathcal{P}_{4}$ | $\mathcal{P}_{5}$ | $\mathcal{P}_{all}$ | $\mathcal{S}_{0.5}$
MDC | 1.00 | 0.79 | 0.84 | 0.00 | 0.85 | 0.00 | 2974.5 |  | 1.00 | 0.64 | 0.80 | 0.00 | 0.85 | 0.00 | 2725.0
CSIS | 0.99 | 0.56 | 0.55 | 1.00 | 0.69 | 0.17 | 554.5 |  | 1.00 | 0.61 | 0.70 | 1.00 | 0.81 | 0.26 | 100.5
CDC-SIS | 0.99 | 0.62 | 0.53 | 0.34 | 0.68 | 0.11 | 1101.5 |  | 1.00 | 0.45 | 0.60 | 0.65 | 0.60 | 0.11 | 962.5
DCSIS2 | 1.00 | 0.84 | 0.57 | 0.00 | 0.52 | 0.00 | 2235.0 |  | 1.00 | 0.90 | 0.71 | 0.05 | 0.70 | 0.01 | 1284.5
CIS | 1.00 | 0.95 | 0.90 | 0.31 | 0.94 | 0.24 | 153.0 |  | 1.00 | 0.76 | 0.78 | 0.42 | 0.78 | 0.17 | 155.0
$\textrm{S-CMC}_{\mathcal{H}}$ | 1.00 | 0.79 | 0.84 | 0.95 | 0.85 | 0.52 | 33.5 |  | 1.00 | 0.63 | 0.78 | 0.90 | 0.82 | 0.35 | 77.0
RaSE1-eBIC | 0.28 | 0.60 | 0.60 | 0.58 | 0.74 | 0.06 | 2250.0 |  | 0.98 | 0.47 | 0.63 | 0.95 | 0.73 | 0.18 | 2253.5

Table 4: The $\mathcal{P}_{i}$, $\mathcal{P}_{all}$ and $\mathcal{S}_{0.5}$ in
Example 4 with $\boldsymbol{X}_{S}$ selected by $\textrm{MC}_{\mathcal{H}}$.

In this example, the active predictors are nonlinearly associated with each
other. Each element in this model is an interaction effect of two variables.
The proposed S-CMCH demonstrates its competitive performance for selecting all
the active predictors with relatively high $\mathcal{P}_{i}$ and
$\mathcal{P}_{all}$, and smallest $\mathcal{S}_{0.5}$, in both scenarios. The
methods DCSIS2 and MDC fail to detect $X_{4}$. It is worth pointing out that
the contribution of variable $X_{4}$ is also underestimated by the methods
CDC-SIS and CIS. This may be because that $X_{4}$’s conditional contribution
to $Y$ is diffused by the exponential term exp$(X_{1})X_{5}$, as indicated by
the much larger $\mathcal{P}_{5}$ in the four methods mentioned above.

###  4.2 Data Applications

####  4.2.1 Single Cell Malt Tumor CITE-seq Dataset

The Malt Tumor Cellular Indexing of Transcriptomes and Epitopes by sequencing
(CITE-seq) dataset (http://www.10xgenomics.com) contains single-cell level
sequencing RNA data and as well as the surface protein expression count. We
are interested in identifying the genes that affect the surface protein level.
The dataset contains 33555 genes and proteins from 8412 single cells.
Following the data pre-processing procedure in Tong et al., (2022) (filtering
out cells with more than 90% zero entries and genes that has zero variance),
we obtain a sample of $n=207$ single cells and $p=18702$ genes. And we set
protein CD8 as the response variable and the protein CD3 as conditional
variable. Interested readers are referred to Tong et al., (2022) for detailed
scientific explanations of using CD3 as the conditional variable. To evaluate
the prediction performance of each screening method, we randomly split the
observations into a training set of size 176 and a test set of size 31. We
select $d_{2}=38$ variables by the following four conditional methods: CSIS,
CDCSIS, CIS and SCMCH including CD3. Then a random forest model is fitted with
the selected variables and the response, with log transformations on the
response and the variables. The mean squared errors (MSE) of the methods on
the test set for 100 repetitions are reported in Table 5. In each of the 100
repetitions, we calculate the difference between the MSE of S-CMCH and that of
each competing method. Then we conducted a one-sided paired two-sample t-test
with the alternative hypothesis that our method S-CMCH has a smaller average
mean squared error (MSE). Our method has the smallest MSE and the $p$-values
of the paired $t$-test indicates S-C$\textrm{MC}_{\mathcal{H}}$ outperforms
each method in prediction accuracy.

Method | MSE | paired $t$-test $p$-value
---|---|---
CSIS | 1.364 | $<2.2\times 10^{-16}$
CDC-SIS | 1.162 | $<2.2\times 10^{-16}$
CIS | 1.268 | $<2.2\times 10^{-16}$
$\textrm{S-CMC}_{\mathcal{H}}$ | 0.748 | –
Table 5: The prediction accuracy in MALT data.

####  4.2.2 Riboflavin Production Dataset

The dataset (Lee et al.,, 2001) contains information about riboflavin (vitamin
B2) production by $n=71$ bacillus subtiliswith, where $p$ = 4088 gene
expression levels are recorded. The dataset is provided by Royal DSM
(Switzerland) and is available in the R package $\mathrm{\mathtt{hdi}}$.

Our goal is to find which genes are most related in predicting the riboflavin
production rate. We randomly split the sample into a training set of size 60
and a test set of size 11. To evaluate the prediction performance of each
screening method, we select $d_{2}=16$ variables for each screening method and
train a random forest model on the training set. Then we calculate the mean
squared error (MSE) of each method on the test set. The average MSE’s based on
100 data splittings into training/test sets are reported in Table 6. Similar
to Section 4.2.1, we conducted the same one-sided paired two-sample t-test
with the alternative hypothesis that our method S-CMCH has a smaller averge
mean squared error (MSE). Our method has the smallest MSE and the $p$-values
suggest that it significantly improves the MSE.

Method | MSE | paired $t$-test $p$-value
---|---|---
MDC | 0.301 | $0.017$
CSIS | 0.357 | $2.909\times 10^{-15}$
CDC-SIS | 0.323 | $5.62\times 10^{-5}$
CIS | 0.323 | $5.62\times 10^{-5}$
RaSE1-eBIC | 0.293 | 0.401
$\textrm{S-CMC}_{\mathcal{H}}$ | 0.292 | –
Table 6: The prediction accuracy in Riboflavin production data. The
conditional set $\boldsymbol{X}_{S}$ contains the top $d_{1}=4$ variables
suggested by $\textrm{MC}_{\mathcal{H}}$.

##  5 Conclusion

In this article, we propose the conditional martingale difference divergence
($\textrm{CMD}_{\mathcal{H}}$) to measure the dependence between a response
variable and a predictor vector given a third vector. It is primarily designed
to overcome the limitation of marginal independence measures. Based on
$\textrm{CMD}_{\mathcal{H}}$, we develop a new screening procedure called
$\textrm{S-CMC}_{\mathcal{H}}$ by combining the merits of the
$\textrm{CMC}_{\mathcal{H}}$ and $\textrm{MC}_{\mathcal{H}}$ for selecting
both marginal and jointly active variables. The proposed framework can be
easily extended to quantile screening. The simulations and real data
applications demonstrate that $\textrm{S-CMC}_{\mathcal{H}}$ has a competitive
and stable performance under variety model settings for mean or quantile
screening. We also provide a data-driven method for selecting the conditional
set $\boldsymbol{X}_{S}$. The limitation of this method is that we do need to
predetermine a proper number of variables in conditional set to get a
satisfactory performance. Using $\lfloor\sqrt{n/\log n}\rfloor$, as done in
our numerical study, may not suffice for the cases when the true underlying
model consists jointly only active variables that depends on a large
conditional set. Designing a variable screening method that is free of tuning
the cardinality of the conditional set is a challenging future research topic.

## Supplementary Materials

The supplementary material contains additional simulation results, and the
proofs of theorems and results introduced in the main article. (.pdf file)

## Acknowledgements

The authors gratefully acknowledge the funding from National Science
Foundation (NSF-CIF-1813330).

## References

  * Barut et al., (2016) Barut, E., Fan, J., and Verhasselt, A. (2016).  Conditional sure independence screening.  Journal of the American Statistical Association, 111(515):1266–1277.
  * Chang et al., (2013) Chang, J., Tang, C. Y., and Wu, Y. (2013).  Marginal empirical likelihood and sure independence feature screening.  The Annals of Statistics, 41(4):2123–2148.
  * Cheng et al., (2014) Cheng, M.-Y., Honda, T., Li, J., and Peng, H. (2014).  Nonparametric independence screening and structure identification for ultra-high dimensional longitudinal data.  The Annals of Statistics, 42(5):1819–1849.
  * Cui et al., (2015) Cui, H., Li, R., and Zhong, W. (2015).  Model-free feature screening for ultrahigh dimensional discriminant analysis.  Journal of the American Statistical Association, 110(510):630–641.  PMID: 26392643.
  * Fan et al., (2011) Fan, J., Feng, Y., and Song, R. (2011).  Nonparametric independence screening in sparse ultra-high-dimensional additive models.  Journal of the American Statistical Association, 106(494):544–557.  PMID: 22279246.
  * Fan and Lv, (2008) Fan, J. and Lv, J. (2008).  Sure independence screening for ultrahigh dimensional feature space.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 70(5):849–911.
  * Fan et al., (2009) Fan, J., Samworth, R., and Wu, Y. (2009).  Ultrahigh dimensional feature selection: beyond the linear model.  The Journal of Machine Learning Research, 10:2013–2038.
  * Fan and Song, (2010) Fan, J. and Song, R. (2010).  Sure independence screening in generalized linear models with NP-dimensionality.  The Annals of Statistics, 38(6):3567–3604.
  * Fan et al., (2016) Fan, Y., Kong, Y., Li, D., and Lv, J. (2016).  Interaction pursuit with feature screening and selection.  arXiv preprint arXiv:1605.08933.
  * Han, (2019) Han, X. (2019).  Nonparametric screening under conditional strictly convex loss for ultrahigh dimensional sparse data.  The Annals of Statistics, 47(4):1995–2022.
  * Hao and Zhang, (2014) Hao, N. and Zhang, H. H. (2014).  Interaction screening for ultrahigh-dimensional data.  Journal of the American Statistical Association, 109(507):1285–1301.
  * He et al., (2013) He, X., Wang, L., and Hong, H. G. (2013).  Quantile-adaptive model-free variable screening for high-dimensional heterogeneous data.  The Annals of Statistics, 41(1):342–369.
  * Kong et al., (2017) Kong, Y., Li, D., Fan, Y., and Lv, J. (2017).  Interaction pursuit in high-dimensional multi-response regression via distance correlation.  The Annals of Statistics, 45(2):897–922.
  * Lee et al., (2001) Lee, J. M., Zhang, S., Saha, S., Santa Anna, S., Jiang, C., and Perkins, J. (2001).  Rna expression analysis using an antisense bacillus subtilis genome array.  Journal of bacteriology, 183(24):7371–7380.
  * (15) Li, G., Peng, H., Zhang, J., and Zhu, L. (2012a).  Robust rank correlation based screening.  The Annals of Statistics, 40(3):1846–1877.
  * (16) Li, R., Zhong, W., and Zhu, L. (2012b).  Feature screening via distance correlation learning.  Journal of the American Statistical Association, 107(499):1129–1139.
  * Ma et al., (2017) Ma, S., Li, R., and Tsai, C. L. (2017).  Variable screening via quantile partial correlation.  Journal of the American Statistical Association, 112(518):650–663.  PMID: 28943683.
  * Mai and Zou, (2013) Mai, Q. and Zou, H. (2013).  The kolmogorov filter for variable screening in high-dimensional binary classification.  Biometrika, 100(1):229–234.
  * Mai and Zou, (2015) Mai, Q. and Zou, H. (2015).  The fused Kolmogorov filter: A nonparametric model-free screening method.  The Annals of Statistics, 43(4):1471–1497.
  * Nandy et al., (2021) Nandy, D., Chiaromonte, F., and Li, R. (2021).  Covariate information number for feature screening in ultrahigh-dimensional supervised problems.  Journal of the American Statistical Association, 0(0):1–14.
  * Pan et al., (2019) Pan, W., Wang, X., Xiao, W., and Zhu, H. (2019).  A generic sure independence screening procedure.  Journal of the American Statistical Association, 114(526):928–937.  PMID: 31692981.
  * Shao and Zhang, (2014) Shao, X. and Zhang, J. (2014).  Martingale difference correlation and its use in high-dimensional variable screening.  Journal of the American Statistical Association, 109(507):1302–1318.
  * Song et al., (2014) Song, R., Yi, F., and Zou, H. (2014).  On varying-coefficient independence screening for high-dimensional varying-coefficient models.  Statistica Sinica, 24(4):1735–1752.
  * Székely et al., (2007) Székely, G. J., Rizzo, M. L., and Bakirov, N. K. (2007).  Measuring and testing dependence by correlation of distances.  The Annals of Statistics, 35(6):2769–2794.
  * Tian and Feng, (2021) Tian, Y. and Feng, Y. (2021).  Rase: A variable screening framework via random subspace ensembles.  Journal of the American Statistical Association, 0(0):1–12.
  * Tong et al., (2022) Tong, Z., Cai, Z., Yang, S., and Li, R. (2022).  Model-free conditional feature screening with fdr control.  Journal of the American Statistical Association, pages 1–13.
  * Wang, (2009) Wang, H. (2009).  Forward regression for ultra-high dimensional variable screening.  Journal of the American Statistical Association, 104(488):1512–1524.
  * Wen et al., (2018) Wen, C., Pan, W., Huang, M., and Wang, X. (2018).  Sure independence screening adjusted for confounding covariates with ultrahigh dimensional data.  Statistica Sinica, 28(1):293–317.
  * Wendland, (2004) Wendland, H. (2004).  Scattered data approximation, volume 17.  Cambridge university press.
  * Xu and Chen, (2014) Xu, C. and Chen, J. (2014).  The sparse mle for ultrahigh-dimensional feature screening.  Journal of the American Statistical Association, 109(507):1257–1269.  PMID: 25382886.
  * Yu, (2013) Yu, B. (2013).  Stability.  Bernoulli, 19(4):1484 – 1500.
  * Zhang and Zhou, (2018) Zhang, S. and Zhou, Y. (2018).  Variable screening for ultrahigh dimensional heterogeneous data via conditional quantile correlations.  Journal of Multivariate Analysis, 165:1–13.
  * Zhao and Li, (2012) Zhao, S. D. and Li, Y. (2012).  Principled sure independence screening for cox models with ultra-high-dimensional covariates.  Journal of Multivariate Analysis, 105(1):397–411.
  * Zhou et al., (2020) Zhou, T., Zhu, L., Xu, C., and Li, R. (2020).  Model-free forward screening via cumulative divergence.  Journal of the American Statistical Association, 115(531):1393–1405.
  * Zhu et al., (2011) Zhu, L. P., Li, L., Li, R., and Zhu, L. X. (2011).  Model-free feature screening for ultrahigh-dimensional data.  Journal of the American Statistical Association, 106(496):1464–1475.  PMID: 22754050.

Department of Statistics, University of Kentucky, Lexington KY 40536-0082,
U.S.A. E-mail<EMAIL_ADDRESS>
Department of Statistics, University of Kentucky, Lexington KY 40536-0082,
U.S.A. E-mail<EMAIL_ADDRESS>
##  6 Appendix

##  Appendix A Bochner’s Theorem

######  Lemma 1.

(Wendland,, 2004, Theorem 6.6) A continuous function
$k:\mathbb{R}^{p+q}\rightarrow\mathbb{R}$ is positive semi-definite if and
only if it is the Fourier transform of a finite nonnegative Borel measure
$W(\xi)d\xi$ on $\mathbb{R}_{{}^{p+q}}$, that is,

$\displaystyle
k(z)=\int_{\mathbb{R}^{(p+q)}}e^{-iz^{T}\xi}W(\xi)d\xi,~{}\forall
z\in\mathbb{R}^{p+q}.$

##  Appendix B Properties of $\textrm{MC}_{\mathcal{H}}$

Some important properties of $\textrm{MD}_{\mathcal{H}}$ and
$\textrm{MC}_{\mathcal{H}}$ are presented as follows.

######  Definition 4.

Given i.i.d. observations $(\boldsymbol{U}_{i},V_{i})_{i=1}^{n}$ from the
distribution of $(\boldsymbol{U},V)$. Let $a_{ij}=V_{i}V_{j}$ and
$b_{ij}=k(\boldsymbol{U}_{i}-\boldsymbol{U}_{j})$, where $k$ is a kernel in
RKHS. The unbiased sample RKHS type martingale difference divergence
$\textnormal{$\widehat{\textrm{MD}_{\mathcal{H}}}$}^{2}(V,\boldsymbol{U})$ is
defined as

$\displaystyle\textnormal{$\widehat{\textrm{MD}_{\mathcal{H}}}$}^{2}(V,\boldsymbol{U})=\dfrac{1}{n(n-3)}\sum_{i\neq
j}^{n}a_{ij}^{*}b_{ij}^{*}$ (B.11)

and the unbiased sample RKHS type martingale difference correlation
$\textnormal{$\widehat{\textrm{MC}_{\mathcal{H}}}$}^{2}(V,\boldsymbol{U})$ is
defined by

$\displaystyle\textnormal{$\widehat{\textrm{MC}_{\mathcal{H}}}$}^{2}(V,\boldsymbol{U})=\begin{cases}\dfrac{\textnormal{$\widehat{\textrm{MD}_{\mathcal{H}}}$}^{2}(V,\boldsymbol{U})}{\textnormal{var}_{n}(V)\textnormal{var}_{n\mathcal{H}}(\boldsymbol{U})}&\text{if
{var}${}_{n}(V)${var}${}_{n\mathcal{H}}(\boldsymbol{U})>0$}\\\ 0&\text{
otherwise,}\end{cases}$ (B.12)

where var${}_{n}(V)=(\frac{1}{n(n-3)}\sum_{i\neq
j}^{n}|a_{ij}^{*}|^{2})^{1/2}$, and
var${}_{n\mathcal{H}}(\boldsymbol{U})=(\frac{1}{n(n-3)}\sum_{i\neq
j}^{n}|b_{ij}^{*}|^{2})^{1/2}$.

######  Theorem 6.

The following properties hold if $E(V^{2})<\infty$ :

  1. a.

$\textnormal{MD}_{\mathcal{H}}^{2}(V,\boldsymbol{U})=E[(V-E(V))(V^{\prime}-E(V^{\prime}))k(\boldsymbol{U}-\boldsymbol{U}^{\prime})]$.

  2. b.

$0\leq\textnormal{MC}_{\mathcal{H}}(V,\boldsymbol{U})\leq 1,$ and
$\textnormal{MC}_{\mathcal{H}}(V,\boldsymbol{U})=0\Leftrightarrow
E(V|\boldsymbol{U})=E(V)$ almost surely.

  3. c.

$\textnormal{MC}_{\mathcal{H}}(a+bV,c+\boldsymbol{U})=\textnormal{MC}_{\mathcal{H}}(V,\boldsymbol{U})$
for any scalars $a$, $b\in\mathbb{R}$ and $c\in\mathbb{R}^{q}$.

######  Theorem 7.

If E $(V^{2})<\infty$, then

$lim_{n\rightarrow\infty}\textnormal{$\widehat{\textrm{MD}_{\mathcal{H}}}$}(V,\boldsymbol{U})=\textnormal{$\textrm{MD}_{\mathcal{H}}$}(V,\boldsymbol{U})\textrm{
a.s.,}$ (B.13)

and

$lim_{n\rightarrow\infty}\textnormal{$\widehat{\textrm{MC}_{\mathcal{H}}}$}(V,\boldsymbol{U})=\textnormal{$\textrm{MC}_{\mathcal{H}}$}(V,\boldsymbol{U})\textrm{
a.s.}.$ (B.14)

######  Theorem 8.

Assume $E(V^{2})<\infty$, we have the following:

  1. a.

If $\textnormal{MC}_{\mathcal{H}}(V,\boldsymbol{U})=0$, then

$n\textnormal{$\widehat{\textrm{MD}_{\mathcal{H}}}$}^{2}(V,\boldsymbol{U})\xrightarrow[n\rightarrow\infty]{D}||\Gamma(s)||_{\mathcal{H}_{k}}^{2},$
(B.15)

where $\Gamma(\cdot)$ denotes a complex-valued zero-mean Gaussian random
process with covariance function

$\displaystyle\textnormal{cov}_{\Gamma}(s,s_{0})$
$\displaystyle=F(s-s_{0})-\textsl{g}_{\boldsymbol{U}}(s-s_{0})E^{2}(V)+\\{E(V^{2})+E^{2}(V)\\}\textsl{g}_{\boldsymbol{U}}(s)\overline{\textsl{g}_{\boldsymbol{U}}(s_{0})}$
$\displaystyle-F(s)\overline{\textsl{g}_{\boldsymbol{U}}(s_{0})}-\textsl{g}_{\boldsymbol{U}}(s)\overline{F(s_{0})}$

with $s,s_{0}\in\mathbb{R}^{q}$ and $g_{\boldsymbol{U}}(s)=E(e^{i\langle
s,\boldsymbol{U}\rangle})$,
$F(s)=E[V^{2}\exp(i\langle\boldsymbol{U},s\rangle)]$.

  2. b.

If $\textnormal{$\textrm{MC}_{\mathcal{H}}$}(V,\boldsymbol{U})=0$ and
$E(V^{2}|\boldsymbol{U})=E(V^{2})$, then

$n\textnormal{$\widehat{\textrm{MD}_{\mathcal{H}}}$}^{2}(V,\boldsymbol{U})/S_{n}\xrightarrow[n\rightarrow\infty]{D}Q,$

where $S_{n}=(1-\frac{1}{n(n-1)}\sum_{k\neq
l}k(\boldsymbol{U}_{k}-\boldsymbol{U}_{l}))(\frac{1}{n}\sum_{k}(V_{k}-\bar{V}_{n})^{2})$,
and $Q$ is a nonnegative quadratic form
$Q=\sum_{i=1}^{\infty}\lambda_{i}Z_{i}^{2}$, where $Z_{i}$ are independent
standard normal random variables. $\\{\lambda_{i}\\}$ are nonnegative
constants that depend on the distribution of (U,V) and $E(Q)=1$.

  3. c.

If $\textnormal{$\textrm{MC}_{\mathcal{H}}$}(V,\boldsymbol{U})\neq 0$, then
$n\textnormal{$\widehat{\textrm{MD}_{\mathcal{H}}}$}^{2}(V,\boldsymbol{U})/S_{n}\xrightarrow[n\rightarrow\infty]{P}\infty$.

###  B.1 Sure Screening Property of $\textrm{MC}_{\mathcal{H}}$

Let $\psi_{i}=\textrm{MC}_{\mathcal{H}}(Y,X_{i})$ for predictor $X_{i}$ and
$\widehat{\psi_{i}}=\widehat{\textrm{MC}_{\mathcal{H}}}(Y,X_{i})$. Denote
$\mathcal{\widehat{M}}=\\{j:\widehat{\psi}_{j}\geq cn^{-\kappa},$ for $1\leq
j\leq p\\}$. Similar to $\textrm{CMC}_{\mathcal{H}}$, we need the following
two assumptions.

(B1)

There exists a positive constant $s_{0}$ such that for all $0<s\leq 2s_{0}$,
then $E\\{$exp$(sY^{2})\\}<\infty$.

(B2)

The minimum $\textrm{MC}_{\mathcal{H}}$ value of active predictors is greather
than $2cn^{-\kappa}$, for some constant $c>0$ and $0\leq\kappa<\dfrac{1}{2}$.

######  Theorem 9.

Under Assumption $(B1)$, for any $0<\gamma<1/2-\kappa$, there exist postive
constants $c_{1}$ and $c_{2}$ such that

$P\Bigg{\\{}\underset{1\leq j\leq p}{\max}|\widehat{\psi}_{j}-\psi_{j}|\geq
cn^{-\kappa}\Bigg{\\}}\leq
O(p[\exp\\{-c_{1}n^{1-2(\kappa+\gamma)}\\}+n\exp(-c_{2}n^{\gamma})]).$ (B.16)

Under conditions $(B1)$ and $(B2)$, we have that

$P(\mathcal{M}\subset\widehat{\mathcal{M}})\geq
1-O(s_{n}[\exp\\{-c_{1}n^{1-2(\kappa+\gamma)}\\}+n\exp(-c_{2}n^{\gamma})]),$
(B.17)

where $s_{n}$ is the cardinality of $\mathcal{M}$.

For the sure screening property of quantile screening by
$\textrm{MC}_{\mathcal{H}}$, we require the condition (C1) in the following
Section C. Denote $\mathcal{M}_{q_{\tau}}:=\\{j:\mathbb{E}(Y_{\tau}|X_{j})\
\text{depends on}\ X_{j}\\}$ and
$\mathcal{\widehat{M}}_{q_{\tau}}:=\\{j:\widehat{\psi}_{j}(\hat{{Y}}_{{\tau}})\geq
cn^{-\kappa},$ for $1\leq j\leq p\\}$.

######  Theorem 10.

Under $(C1)$, for any $0<\gamma<1/2-\kappa$ and $\kappa\in(0,1/2)$, there
exists positive constants $c_{1},c_{2}$ such that for any $c>0$,

$P\Bigg{\\{}\underset{1\leq j\leq
p}{\max}|\widehat{\psi}_{j}(\hat{{Y}}_{{\tau}})-\psi_{j}({Y}_{{\tau}})|\geq
cn^{-\kappa}\Bigg{\\}}\leq
O(p[\exp\\{-c_{1}n^{1-2(\kappa+\gamma)}\\}+n\exp(-c_{2}n^{\gamma})]).$ (B.18)

If the minimum $\textrm{MC}_{\mathcal{H}}$ value of active predictors
satisfies $\min_{j\in\mathcal{M}_{q_{\tau}}}\psi_{j}({Y}_{{\tau}})\geq
2cn^{-\kappa}$ for some constant $c>0$ and $0\leq\kappa<1/2$, we can show that

$P(\mathcal{M}_{q_{\tau}}\subseteq\widehat{\mathcal{M}}_{q_{\tau}})\geq
1-O(\widetilde{s}_{n}[\exp\\{-c_{1}n^{1-2(\kappa+\gamma)}\\}+n\exp(-c_{2}n^{\gamma})]),$
(B.19)

where $\widetilde{s_{n}}$ is the cardinality of $\mathcal{M}_{q_{\tau}}$.

##  Appendix C Sure Screening Property of Quantile Screening using
$\textrm{CMC}_{\mathcal{H}}$

We require the two assumptions below:

(C1)

The CDF of $Y$ ($F_{Y}$) is continuously differentiable in a small
neighborhood of $q_{\tau}=q_{\tau}(Y)$, say
$[q_{\tau}-\delta_{0},q_{\tau}+\delta_{0}]$ for $\delta>0$. Let
$G_{1}(\delta_{0})=\inf{}_{y\in[q_{\tau}-\delta_{0},q_{\tau}+\delta_{0}]}f_{Y}(y)$,
and
$G_{2}(\delta_{0})=\sup{}_{y\in[q_{\tau}-\delta_{0},q_{\tau}+\delta_{0}]}f_{Y}(y)$
where $f_{Y}$ is the density function of $Y$. Assume that
$0<G_{1}(\delta_{0})\leq G_{2}(\delta_{0})<\infty$.

(C2)

The minimum C$\textrm{MC}_{\mathcal{H}}$ value of active predictors satisfies
$\min_{j\in\mathcal{D}_{q_{\tau}}}\omega_{j}({Y}_{{\tau}})\geq 2cn^{-\kappa}$
for some constant $c>0$ and $0\leq\kappa<1/2$.

The following proposition from Shao and Zhang, (2014) is necessary for proving
the sure screening property.

######  Proposition 1.

Under condition (C1), there exists $\epsilon_{0}>0$ and $c_{1}>0$, such that
for any $\epsilon\in(0,\epsilon_{0})$,

$P\bigg{(}\dfrac{1}{n}\sum_{l=1}^{n}|\hat{y}_{l_{\tau}}-y_{l_{\tau}}|>\epsilon\bigg{)}\leq
3\exp(-2nc_{1}\epsilon^{2})$ (C.20)