# Screening methods for linear errors-in-variables models
in high dimensions

Linh H. Nghiem  Research School of Finance, Actuarial Studies and Statistics,
Australian National University, Australia  School of Mathematics and
Statistics, University of Sydney, Australia  Francis K.C. Hui  Research School
of Finance, Actuarial Studies and Statistics, Australian National University,
Australia  Samuel Müller  Department of Mathematics and Statistics, Macquarie
University, Australia  School of Mathematics and Statistics, University of
Sydney, Australia  A.H. Welsh  Research School of Finance, Actuarial Studies
and Statistics, Australian National University, Australia

###### Abstract

Microarray studies, in order to identify genes associated with an outcome of
interest, usually produce noisy measurements for a large number of gene
expression features from a small number of subjects. One common approach to
analyzing such high-dimensional data is to use linear errors-in-variables
models; however, current methods for fitting such models are computationally
expensive. In this paper, we present two efficient screening procedures,
namely corrected penalized marginal screening and corrected sure independence
screening, to reduce the number of variables for final model building. Both
screening procedures are based on fitting corrected marginal regression models
relating the outcome to each contaminated covariate separately, which can be
computed efficiently even with a large number of features. Under mild
conditions, we show that these procedures achieve screening consistency and
reduce the number of features considerably, even when the number of covariates
grows exponentially with the sample size. Additionally, if the true covariates
are weakly correlated, corrected penalized marginal screening can achieve full
variable selection consistency. Through simulation studies and an analysis of
gene expression data for bone mineral density of Norwegian women, we
demonstrate that the two new screening procedures make estimation of linear
errors-in-variables models computationally scalable in high dimensional
settings, and improve finite sample estimation and selection performance
compared with estimators that do not employ a screening stage.

Keywords: dimension reduction; forward regression; measurement error;
penalized regression; regularization; sure independence screening

##  1 Introduction

In microarray studies, to identify genes that are associated with an outcome
of interest, a large number of gene expressions (potentially tens of
thousands) are measured from typically a much smaller number of subjects
(often in the tens to hundreds). The gene expression measurements tend to be
noisy, where measurement errors come from many sources such as sample
preparation, labeling, and hybridization; for example, see Rocke and Durbin,
(2001) and Zakharkin et al., (2005). The gene measurements are also often
analyzed on the log scale, making the assumption of additive measurement
errors more plausible (Nghiem and Potgieter,, 2019). Furthermore, as in common
genome wide association studies (Do et al.,, 2011; Zhou et al.,, 2018, among
others), it is usually assumed that only a few genes are related to the
outcome of interest, i.e a sparsity assumption on the statistical model. As a
specific motivating example, after some preprocessing steps, our Bone Mineral
Density data in Section 4 contains noisy measurements of $p=993$ features
(genes) from $n=84$ observations (Norwegian women), and we are interested in
identifying genes that are associated with the total hip T-score.

For analyzing such data, a commonly used approach which we focus on in this
paper is the classical linear errors-in-variables (EIV) model

$\mathbf{y}=\operatorname{\mathbf{X}}\bm{\beta}_{0}+\bm{\varepsilon},\quad\operatorname{\mathbf{W}}=\operatorname{\mathbf{X}}+{\operatorname{\mathbf{U}}},$
(1)

where $\mathbf{y}\in\mathbb{R}^{n}$ is a random vector of outcomes from $n$
independent and identically distributed (iid) observations,
$\mathbf{X}\in\mathbb{R}^{n\times p}$ is the deterministic true covariate
matrix (typically the true gene expressions),
$\operatorname{\bm{\beta}}_{0}=(\beta_{01},\ldots,\beta_{0p})^{\top}\in\mathbb{R}^{p}$
is the true coefficient vector, and $\bm{\varepsilon}\in\mathbb{R}^{n}$ is the
model error vector term whose components are assumed to be iid with zero mean
and variance $\sigma^{2}$. Due to the existence of measurement error, the true
covariate matrix $\operatorname{\mathbf{X}}$ is not observed; instead, we
observe the random matrix $\mathbf{W}\in\mathbb{R}^{n\times p}$, which is a
noisy version of $\mathbf{X}$, contaminated by an additive random measurement
error matrix $\mathbf{U}\in\mathbb{R}^{n\times p}$ independent of
$\bm{\varepsilon}$. The rows of $\mathbf{U}$ are assumed to be iid random
vectors with zero mean and covariance matrix $\bm{\Sigma}_{u}$. We focus on
the model (1) in high dimensional settings, where the number of covariates $p$
can be bigger than and grow with the sample size $n$, potentially at an
exponential rate. We also assume the true coefficient vector
$\operatorname{\bm{\beta}}_{0}$ to be sparse, meaning that only a few
components of $\operatorname{\bm{\beta}}_{0}$ are non-zero.

When the true covariate matrix $\operatorname{\mathbf{X}}$ is observed,
penalized regression methods (Tibshirani,, 1996; Fan et al.,, 2004; Huang et
al.,, 2008; Simon et al.,, 2013; Piironen et al.,, 2017; Ida et al.,, 2019,
among others) are widely used to estimate and perform variable selection on
$\operatorname{\bm{\beta}}_{0}$. However, when the true covariate matrix
$\mathbf{X}$ is not observed, replacing $\mathbf{X}$ with $\mathbf{W}$ leads
to a naive estimator that is inconsistent in both estimation and variable
selection of $\operatorname{\bm{\beta}}_{0}$ (Sørensen et al.,, 2015). To
address this challenge, several corrections for measurement error in high
dimensional linear EIV models have been proposed. For example, Rosenbaum et
al., (2010) and Rosenbaum et al., (2013) proposed the matrix uncertainty (MU)
selector and its improved version, respectively, while Belloni et al., (2017)
proved its near-optimal minimax properties and developed a conic programming
estimator that can achieve the minimax bound. Both the MU selector and conic
estimator require appropriate choice of multiple tuning parameters, which is
typically very challenging in practice, especially in high dimensional
settings. Another approach for handling measurement error is to modify the
loss function or the conditional score functions commonly seen in the error-
free penalized regressions; examples include the corrected lasso method of Loh
and Wainwright, (2012) and Sørensen et al., (2015), and the convex conditioned
lasso of Datta et al., (2017). More recently, Romeo and Thoresen, (2019)
presented a simulation study to compare the performance of the MU, corrected
lasso, and convex conditioned lasso against the naive estimator in both
estimation and variable selection; they concluded that the relative
performance of those estimators depend on the structure of $\bm{\Sigma}_{u}$.
Brown et al., (2019) introduced a boosting algorithm based on the estimating
equation of the corrected lasso, but the theoretical properties of the final
estimates were not examined. In another line of research, Nghiem and
Potgieter, (2019) proposed a SIMSELEX estimator that first uses simulation to
evaluate the effect of measurement error on estimated coefficients, then
selects important covariates based on these simulated effects, and finally
extrapolates the simulation to the scenario with no measurement error present.
Furthermore, Byrd and McGee, (2019) presented an EM-type correction method,
where they iteratively sample the true covariate from the conditional
distribution of $\operatorname{\mathbf{X}}$ given $\mathbf{y}$ and
$\operatorname{\mathbf{W}}$, and fitted penalized regressions of $\mathbf{y}$
on these sampled covariates. A major disadvantage of all these methods is that
they are not computationally efficient when the number of covariates $p$ is
very large. Specifically, while the corrected lasso is defined to be a global
minimum of a non-convex optimization problem, the convex conditioned lasso
requires computation of the nearest semi-positive definite matrix measured in
element-wise max norm in high dimensions; we will elaborate these two
estimators in Section 3.2. In addition, the SIMSELEX procedure requires
running the lasso on a large number of simulated datasets, while the method of
Byrd and McGee, (2019) requires sampling from a large $p$-dimensional
multivariate distribution, which is slow in high dimensional settings.

In this paper, motivated by the Bone Mineral Density data, we address the
large $p$ problem for linear EIV models by proposing two efficient corrected
marginal screening methods, namely _corrected penalized marginal screening_
and _corrected sure independence screening_ , respectively. These screening
methods aim to quickly identify a screening index set
$\hat{\mathcal{Q}}\subset\\{1,\ldots,p\\}$ from the observed data, such that
$\hat{\mathcal{Q}}$ has much smaller cardinality than $p$ but still retains
all the important covariates, a property known as screening consistency. When
no measurement error is present, screening methods in high dimensions are
usually carried out based on marginal regressions, with the first two proposed
methods being the penalized marginal bridge regression of Huang et al., (2008)
and sure independence screening of Fan and Lv, (2008) for the linear model,
followed by a vast literature that improved these two methods and applied them
to more complex models, see Fan and Lv, (2008, among others); Fan et al.,
(2010, among others); Li et al., (2012, among others); Barut et al., (2016,
among others); Wen et al., (2018, among others). Compared to penalized methods
when no measurement error is present, screening methods have received far less
attention for EIV models. A screening method for linear EIV models was briefly
mentioned in Kaul et al., (2016). As detailed in Section 2, this method is a
special case of our proposed corrected sure independence screening method when
all the covariates have the same measurement error variance. In our new
methods, after screening, we only compute penalized estimators using the
variables indexed by $\hat{\mathcal{Q}}$; since the cardinality of the set
$\hat{\mathcal{Q}}$ is much smaller than $p$, the total estimation times of
two-stage estimators are much reduced compared to those of one-stage
estimators that do not employ a screening. We demonstrate that the benefit of
the proposed screening procedures is so substantial that it can make many
estimators (such as the convex conditioned lasso) computationally feasible in
high dimensional settings. Theoretically, we show that under mild conditions,
even when the number of covariates $p$ grows at an exponential rate with the
sample size $n$, our proposed screening methods achieve screening consistency
while reducing the number of variables to below the sample size. Moreover,
under a partial orthogonality condition for the true covariates, we
demonstrate that the corrected penalized marginal screening approach can
achieve full variable selection consistency. Our simulation studies and an
analysis of the motivating Bone Mineral Density data both verify our
theoretical results, and demonstrate that our proposed screening procedures
lead to remarkable gains in both computational cost and finite sample
performance.

The remainder of this paper is organized as follows. Section 2 introduces the
proposed screening procedures and establishes their theoretical properties. In
Section 3, we present simulation studies to demonstrate the strong empirical
performance of screening procedures and several two-stage estimators. Section
4 applies the methodologies to analyze the motivating Bone Mineral Density
data, and Section 5 offers some concluding remarks.

The following notation is used throughout the paper. For a generic matrix
$\mathbf{A}$, let $A_{ij}$ denote the $(i,j)$ element of $\mathbf{A}$ and let
$\lVert\mathbf{A}\rVert_{2}$, $\lVert\mathbf{A}\rVert_{F}$, and
$\lVert\mathbf{A}\rVert_{\text{max}}$ denote the $\ell_{2}$ norm, Frobenius
norm, and element-wise max norm, respectively. For a square matrix
$\mathbf{A}$, let $\lambda_{\text{max}}(\mathbf{A})$ denote its maximum
eigenvalue. For any vector $\mathbf{v}$, let $v_{j}$ denote its $j$th
component, $\lVert\mathbf{v}\rVert_{2}$ and $\lVert\mathbf{v}\rVert_{1}$
denote its $\ell_{2}$ norm and $\ell_{1}$ norm, respectively. For any set
$\mathcal{S}$, let $\mathcal{S}^{c}$ denote its complement, and
$|\mathcal{S}|$ denote its cardinality. Finally, for any sequence $a_{n}$ and
$b_{n}$, we write $a_{n}\sim b_{n}$ if there exist positive constants $c_{1}$
and $c_{2}$ such that $c_{1}a_{n}\leq b_{n}\leq c_{2}a_{n}$.

##  2 Corrected marginal screening procedures

Consider the linear EIV model (1), with the observed data consisting of the
outcome vector $\mathbf{y}\in\mathbb{R}^{n}$ and the surrogate matrix
$\operatorname{\mathbf{W}}\in\mathbb{R}^{n\times p}$. Let
$\sigma_{1}^{2},\ldots,\sigma_{p}^{2}$ denote the diagonal elements of
$\bm{\Sigma}_{u}$; if any true covariate $X_{ij}$ is measured without error,
then the corresponding element $\sigma_{j}^{2}$ is set to zero. We assume only
$s$ components of $\operatorname{\bm{\beta}}_{0}$ are non-zero, where
$s\ll\min(n,p)$, and without loss of generality, let
$\mathcal{S}=\left\\{1,\ldots,s\right\\}$ and
$\mathcal{S}^{c}=\left\\{s+1,\ldots,p\right\\}$ denote the set of indices
corresponding to non-zero and zero components of $\bm{\beta}_{0}$,
respectively. For the remainder of this paper, unless otherwise stated, we
assume the covariance matrix $\bm{\Sigma}_{u}$ is known; in practice,
$\bm{\Sigma}_{u}$ is usually estimated from replicate data (Carroll et al.,,
2006). For example, in our microarray data analysis in Section 4, replications
are available in the form of multiple probes for each gene expression, and we
follow a common procedure to estimate $\bm{\Sigma}_{u}$ from these data
(further details are presented in Appendix B).

We develop screening procedures that can reduce the number of covariates for
model (1) but still maintain all the important variables. Moreover, these
procedures are designed to be computationally scalable to high dimensional
settings with $p\gg n$. Specifically, we propose to screen variables by
minimizing

$L(\bm{\beta})=\sum_{j=1}^{p}L_{j}(\beta_{j})=\frac{1}{n}\sum_{j=1}^{p}\left\\{\sum_{i=1}^{n}\left(y_{i}-W_{ij}\beta_{j}\right)^{2}-{\sigma}_{j}^{2}\beta_{j}^{2}+\lambda_{n}p_{n}(|\beta_{j}|)\right\\},$
(2)

with $\lambda_{n}$ being a non-negative tuning parameter and
$p_{n}(|\beta_{j}|)$ a penalty function on $|\beta_{j}|$. Let
$\hat{\bm{\beta}}=\operatorname*{argmin}_{\bm{\beta}}L(\operatorname{\bm{\beta}})$.
The function $L(\bm{\beta})$ to be minimized in (2) consists of two parts. The
first part,
$\sum_{j=1}^{p}n^{-1}\sum_{i=1}^{n}\left(y_{i}-W_{ij}\beta_{j}\right)^{2}-\sum_{j=1}^{p}{\sigma}_{j}^{2}\beta_{j}^{2}$,
is the sum of all the $\ell_{2}$ losses for the regression of the outcome on
each surrogate predictor (i.e column of $\operatorname{\mathbf{W}}$)
separately, with measurement error accounted for via a (negative) $\ell_{2}$
penalty term, $-\sum_{j=1}^{p}\sigma_{j}^{2}\beta_{j}^{2}$. When $p$ is fixed
and $n\to\infty$, this loss part converges to
$n^{-1}\sum_{j=1}^{p}\sum_{i=1}^{n}\left(y_{i}-X_{ij}\beta_{j}\right)^{2}$,
which is the loss that was used by both Huang et al., (2008) and Fan et al.,
(2010) for their corresponding marginal screening procedures when no
measurement error is present. In the theoretical analysis below, we allow $p$
to diverge to $\infty$; in this case, this loss part still enables us to
achieve screening consistency. For the second part of (2), the penalty
function $\lambda_{n}\sum_{j=1}^{p}p_{n}(|\beta_{j}|)$ regularizes the
estimates, i.e with appropriate choices for $\lambda_{n}$, many components of
the estimated $\hat{\bm{\beta}}$ are set to zero. Common choices of the
penalty function $p_{n}(|\beta_{j}|)$ include the lasso $|\beta_{j}|$, the
bridge penalty $|\beta_{j}|^{\alpha}$ with $0<\alpha<1$ (Frank and Friedman,,
1993), and the SCAD penalty (Fan et al.,, 2004), among others. In this paper,
we use the bridge penalty, which makes the solution of (2) relatively fast to
compute and easy to analyze theoretically. However, we note that other
penalties can be used in practice, although pursuing this is beyond the scope
of this paper. Computationally, for any value of $\lambda_{n}$, each component
$\hat{\beta}_{j}$ can be obtained by minimizing $L_{j}(\beta_{j})$ separately.
When $n$ is large enough, the coefficient
$n^{-1}\sum_{i=1}^{n}W_{ij}^{2}-\sigma_{j}^{2}$ associated with
$\beta_{j}^{2}$ is positive, so $L_{j}(\beta_{j})$ is a convex function of
$\beta_{j}$ and can be minimized quickly using univariate convex optimization
routines. Therefore, the minimizer of (2) can be computed efficiently even
with very large $p$, making it and subsequent screening methods scalable in
high dimensions.

We consider two screening approaches arising from (2). For the first approach,
which will be referred to as _corrected penalized marginal screening_ (PMSc),
we consider the case when the tuning parameter $\lambda_{n}$ is strictly
positive. In this case, although the minimizer of (2) does not generally have
a closed form, some components of $\hat{\bm{\beta}}$ can be set exactly to
zero, so the corresponding screening set is defined to be
$\hat{\mathcal{Q}}_{\text{PMSc}}=\\{j:\hat{\beta}_{j}\neq 0\\}$. For the
second approach, we consider the case when $\lambda_{n}=0$, so the minimizer
of (2) has components

$\tilde{\beta}_{j}=\dfrac{\sum_{i=1}^{n}W_{ij}y_{i}}{\sum_{i=1}^{n}W_{ij}^{2}-n\sigma_{j}^{2}},\quad
j=1,\ldots,p,$ (3)

each of which is non-zero with probability one. Note that each component
$\tilde{\beta}_{j}$ is a consistent estimator of the slope $\gamma_{j}$ in the
marginal univariate model between the outcome and the true (unobserved)
covariate $y_{i}=\gamma_{j}X_{ij}+\varepsilon_{i},~{}i=1,\ldots,n$, which is
the quantity that is used by Fan and Lv, (2008) and Fan et al., (2010) in sure
independence screening (SIS) when no measurement error exists. In order to
reduce the number of dimensions in this unpenalized approach, similar to SIS,
we keep the $d$ components with largest magnitude
$|\tilde{\beta}_{j}|,j=1,\ldots,p$, and refer to this approach as _corrected
sure independence screening_ (SISc). The corresponding screening set is then
defined to be

$\hat{\mathcal{Q}}_{\text{SISc}}=\\{1\leq j\leq p:|\tilde{\beta}_{j}|\text{ is
among the first}~{}d~{}\text{largest of all
components}~{}|\tilde{\beta}_{1}|,\ldots,|\tilde{\beta}_{p}|\\}.$ (4)

We remark that the SISc approach in our paper is more general than the
corrected screening method introduced by Kaul et al., (2016) in the first step
of their two-step estimation procedure; specifically, Kaul et al., (2016)
screened the variables based on
$\zeta_{j}=\sum_{i=1}^{n}W_{ij}y_{i},~{}j=1,\ldots,p$, which is the numerator
of (3). The rank of $|\zeta_{j}|$ is asymptotically the same as the rank of
$|\tilde{\beta}_{j}|$ if all the covariates have the same measurement error
variance and are measured on the same scale. Therefore, the method of Kaul et
al., (2016) may be considered as a special case of SISc when all the
covariates have the same measurement error variances.

In the next two subsections, we study the theoretical properties of the two
screening procedures PMSc and SISc separately; all the proofs are provided in
Appendix A. Throughout the development below, we note that the true covariate
matrix $\operatorname{\mathbf{X}}$ is assumed to be deterministic; this
assumption could be relaxed to allow $\operatorname{\mathbf{X}}$ to be random,
although we do not explore this in the paper.

###  2.1 Corrected penalized marginal bridge screening

First, we study the properties of the PMSc approach with the penalty function
$p_{n}(|\beta_{j}|)=|\beta_{j}|^{\alpha}$ for $0<\alpha<1$. Since the
screening step aims only to reduce the dimension of the feature space, we are
interested only in variable selection properties, noting that
$\hat{\operatorname{\bm{\beta}}}$ is generally not estimation consistent for
the true vector $\operatorname{\bm{\beta}}_{0}$. For the theoretical
development below, we define

$\tilde{\xi}_{nj}=\frac{1}{n}\mathbb{E}\left(\sum_{i=1}^{n}y_{i}X_{ij}\right)=\frac{1}{n}\sum_{i=1}^{n}\left(\sum_{k=1}^{s}X_{ik}{\beta}_{0k}\right)X_{ij},$
$\xi_{nj}=\frac{1}{n}\mathbb{E}\left[\left(\sum_{i=1}^{n}y_{i}W_{ij}\right)\bigg{\rvert}~{}U_{ij}\right]=\frac{1}{n}\sum_{i=1}^{n}\left(\sum_{k=1}^{s}X_{ik}{\beta}_{0k}\right)W_{ij}=\tilde{\xi}_{nj}+\frac{1}{n}\sum_{i=1}^{n}\left(\sum_{k=1}^{s}X_{ik}{\beta}_{0k}\right)U_{ij}.$

Here, $\tilde{\xi}_{nj}$ can be considered as the marginal covariance between
the outcome and the $j$th true covariate. As in the no measurement error case,
this quantity $\tilde{\xi}_{nj}$ with $j\in\mathcal{S}$ will play an essential
role in determining whether the component $\hat{\beta}_{j}$ is asymptotically
non-zero. However, because the $X_{ij}$’s are not observed, the quantities
$\tilde{\xi}_{nj}$ are not directly computable, so the random variables
$\xi_{nj}$ are used as surrogates for $\tilde{\xi}_{nj}$, $j=1,\ldots,p.$
Next, we assume the following conditions.

  1. (C1)

The model error terms $\varepsilon_{1},\varepsilon_{2},\ldots\varepsilon_{n}$
are iid sub-Gaussian random variables with mean zero and finite variance
$\sigma^{2}$,

  2. (C2)

The measurement errors $U_{ij},~{}i=1,\ldots,n$ are independent sub-Gaussian
random variables with finite variance $\sigma_{j}^{2}$, $j=1,\ldots,p$.

  3. (C3)

There exist constants $b_{0}$ and $b_{1}$ such that $\displaystyle
0<b_{0}\leq\min_{k\in\mathcal{S}}\left|\beta_{0k}\right|\leq\max_{k\in\mathcal{S}}\left|\beta_{0k}\right|\leq
b_{1}<\infty.$

  4. (C4)

There exist constants $C_{1}$ and $C_{2}$ such that $\displaystyle
0<C_{1}\leq\min_{j}n^{-1}\sum_{i=1}^{n}X_{ij}^{2}<\max_{j}n^{-1}\sum_{i=1}^{n}X_{ij}^{2}\leq
C_{2}<\infty.$

  5. (C5)

There exists a constant $\xi_{0}$ such that
$\min_{k\in\mathcal{S}}|\tilde{\xi}_{nj}|\geq\xi_{0}>0$ for all $n$.

Conditions (C(C11))-(C(C11)) are standard in the high dimensional measurement
error literature; e.g these conditions were used in Rosenbaum et al., (2013)
and Belloni et al., (2017). Condition (C(C11)) requires that all the true non-
zero components of $\operatorname{\bm{\beta}}_{0}$ are bounded away from zero
and infinity. Condition (C(C11)) assumes that each true covariate is well-
controlled; in particular, the lower bound implies that the “signal” part on
each covariate is not too small (for example, the true covariate matrix
$\operatorname{\mathbf{X}}$ should not be too sparse). With the above
conditions, the PMSc procedure can achieve screening consistency, as made
precise in Theorem 1.

######  Theorem 1.

Consider the linear EIV model (1) under conditions (C(C11))-(C(C11)), and the
screening set $\hat{\mathcal{Q}}_{\text{PMSc}}=\\{j:\hat{\beta}_{j}\neq 0\\}$,
where $\hat{\bm{\beta}}$ minimizes (2) with the tuning parameter
$\lambda_{n}\rightarrow 0$ and the number of important variables $s=o(n)$.
Then, we have
$\mathbb{P}(\hat{\mathcal{Q}}_{\text{PMSc}}\supseteq\mathcal{S})\rightarrow
1,$ meaning all the important variables are included in the screening set with
probability tending to one as $n\to\infty$.

The conditions for screening consistency as stated in Theorem 1 are relatively
weak. Indeed, it is easy to see that setting the tuning parameter
$\lambda_{n}=0$, such as in the case of SISc, always guarantees screening
consistency, although no dimension reduction is achieved in this case. If the
true important and non-important covariates are weakly correlated, as
formalized in the condition (C(C11)) below, then the PMSc procedure can
achieve full selection consistency, meaning that the screening set contains
important variables only. Specifically, we assume the following additional
conditions.

  1. (C6)

There exists constants $d_{0}>0$ and $0\leq\theta\leq 1/2$ such that
$\displaystyle\left|n^{-1}\sum_{i=1}^{n}X_{ij}X_{ik}\right|\leq
d_{0}n^{-\theta},~{}~{}j\in\mathcal{S}^{c},~{}k\in\mathcal{S}.$

  2. (C7)

Assume the tuning parameter $\lambda_{n}$ satisfies:

    1. (a)

$\lambda_{n}\to 0$ and
$\lambda_{n}n^{\theta(2-\alpha)}s^{\alpha-2}\rightarrow\infty$,

    2. (b)

$\log\left(2m\right)=o(1)\times\left\\{\lambda_{n}n^{(2-\alpha)/2}\right\\}^{1/(2-\alpha)},$
where $m=p-s$.

Condition (C(C11)) is similar to the partial orthogonality condition of Huang
et al., (2008), which is necessary to establish the full variable selection
consistency of the penalized marginal bridge estimator when no measurement
error is present. Such a condition is also necessary when measurement errors
are present. Condition (C(C11)C(C11)a) implies that the number of important
variables $s$ is of order $o(n^{\theta})$ (hence is of order $o(n^{1/2})$ at
most), while condition (C(C11)C(C11)b) implies the number of zero coefficients
$m=p-s$ can be of order $o(\exp(n^{1/2}))$. As a result, the number of
variables $p$ is allowed to grow at an exponential rate with $n^{1/2}$, while
the number of important variables $s$ grows at a comparably slower rate. These
two conditions also control the behavior of the tuning parameter $\lambda_{n}$
that should be used for PMSc to achieve full variable selection consistency.

######  Theorem 2.

Consider the linear EIV model (1) under conditions (C(C11))-(C(C5)), and the
screening set $\hat{\mathcal{Q}}_{\text{PMSc}}=\\{j:\hat{\beta}_{j}\neq 0\\}$.
Then we have
$\mathbb{P}(\hat{\mathcal{Q}}_{\text{PMSc}}=\mathcal{S})\rightarrow 1,$
meaning the screening set $\hat{\mathcal{Q}}_{\text{PMSc}}$ contains all and
only the important variables as $n\to\infty$.

Note that in both Theorem 1 and Theorem 2, we do not make any assumption on
the structure of the covariance matrix $\bm{\Sigma}_{u}$, except that the
measurement error variance on each covariate is bounded for all
$j=1,\ldots,p$, as formalized in condition (C(C11)). In Appendix C, we
demonstrate empirically that the conclusion from Theorem 2 holds even when the
measurement error among the covariates is highly correlated.

In practice, since the true covariate matrix $\operatorname{\mathbf{X}}$ is
not observed, condition (C(C11)) is hard to verify, and even when it is
assumed to hold, we stress that the theoretical choice of $\lambda_{n}$ in
condition (C(C5)) depends on unknown quantities. Therefore, in practice, for
the proposed PMSc procedure, we use cross-validation to select the tuning
parameter $\lambda_{n}$, where the loss on test data is computed based on the
loss part of (2). This procedure will be referred to as PMSc cross-validation
(PMSc${}_{\text{CV}}$) and is further elaborated upon in Subsection 3.2.
Alternatively, we can start with a sufficiently large value for $\lambda_{n}$
such that all the coefficients are set to zero, then incrementally decrease
$\lambda_{n}$ until a certain number of covariates, $M$, is included. Note
that with this version of PMSc, rather than having to choose the tuning
parameter $\lambda_{n}$, we instead choose $M$, the desired number of
covariates to keep from the screening procedure. This strategy will be
referred to as PMSc forward stepwise (PMSc${}_{\text{FS}}$). As demonstrated
in Appendix C, when the true model is sparse as is assumed in our high
dimensional settings, in order to retain all the important variables, $M$ can
be substantially smaller than the number of covariates $d$ to keep in the SISc
approach, and hence also be substantially smaller than the sample size $n$.

###  2.2 Corrected sure independence screening

Because SISc corresponds to the minimizer of (2) with $\lambda_{n}=0$, it
screens out important variables by keeping only the $d$ components with
largest magnitude among $\\{|\tilde{\beta}_{1}|,\ldots,|\tilde{\beta}_{p}|\\}$
as formalized by $\hat{\mathcal{Q}}_{\text{SISc}}$ in (4). In this section, we
prove that under certain conditions, the scalar $d$ can be chosen such that we
can reduce the number of dimensions to below the sample size $n$ and still
retain all the important variables. To achieve this, we follow the theoretical
development of SIS in the linear model with no measurement error by Fan and
Lv, (2008), and aim to theoretically set $d=\lfloor\gamma n\rfloor$, where
$\gamma$ can be of order $n^{-\tilde{\theta}}$ for some $0<\tilde{\theta}<1$.
In addition to the conditions (C(C11))-(C(C11)) for screening consistency as
in Theorem 1, we impose the following additional assumptions.

  1. (C8)

$p>n$, $\log(p)/n\to C_{1}$ with $C_{1}$ being a constant as both $n$ and $p$
diverge to $\infty$.

  2. (C9)

Let $\mathbf{Z}=\operatorname{\mathbf{U}}\bm{\Sigma}_{u}^{-1/2}$. Then,

    1. (a)

For some constants $c_{1},c_{2}>1$ and $D>0$, the matrix $\mathbf{Z}$ follows
a spherical distribution satisfying
$\mathbb{P}\left\\{\lambda_{\max}\left(\tilde{p}^{-1}\tilde{\mathbf{Z}}\widetilde{\mathbf{Z}}^{\top}\right)>c_{2}\text{
and
}\lambda_{\min}\left(\widetilde{p}^{-1}\tilde{\mathbf{Z}}\tilde{\mathbf{Z}}^{\top}\right)<1/c_{2}\right\\}\leq
e^{-Dn},~{}$for any $n\times\tilde{p}$ submatrix $\tilde{\mathbf{Z}}$ of
$\mathbf{Z}$ with $c_{1}n\leq\tilde{p}\leq p$.

    2. (b)

There exists positive constants $\tau_{1}>0$ and $c_{3}>0$, such that
$\lambda_{\max}(\bm{\Sigma}_{u})\leq c_{3}n^{\tau_{1}}.$

  3. (C10)

There exists positive constants $c_{4}>0$ and $\tau_{2}>0$, such that
$\lambda_{\text{max}}\left(n^{-1}\operatorname{\mathbf{X}}^{\top}\operatorname{\mathbf{X}}\right)\leq
c_{4}n^{\tau_{2}}.$

  4. (C11)

$\tau_{1}+\tau_{2}+\log_{n}(s)<1$, where $\log_{n}(s)$ denotes the logarithm
base $n$ of $s$.

Condition (C(C11)) implies that we can allow the number of covariates $p$ to
grow exponentially with the sample size $n$. Condition (C(C11)C(C11)a) is
referred to as the Concentration Property by Fan and Lv, (2008), meaning that
with large probability, the $n$ nonzero singular values of the
$n\times\tilde{p}$ submatrix $\tilde{\mathbf{Z}}$ of $\mathbf{Z}$ are of the
same order; Fan and Lv, (2008) suggested that a sufficient condition for
(C(C11)C(C11)a) is that each row $\mathbf{U}_{i}$ follows a $p$-variate
Gaussian distribution. Conditions (C(C11)C(C11)b) and (C(C11)) indicate that
the maximum eigenvalues of the covariance measurement error matrix
$\bm{\Sigma}_{u}$ and of the scaled Gram matrix
$n^{-1}\operatorname{\mathbf{X}}^{\top}\operatorname{\mathbf{X}}$ can only
grow polynominally with the sample size $n$. Furthermore, condition (C(C11))
restricts the degree of that polynomial growth to be less than one, which
rules out the case when either true covariates or measurement errors are
highly correlated. Condition (C(C11)) also implies that the number of
important covariates $s$ is smaller than the sample size $n$, so that reducing
the dimension to below $n$ is reasonable. Compared with the theoretical
conditions used to establish screening and full variable selection consistency
of PMSc (Theorems 1 and 2), conditions (C(C11))-(C(C11)) impose stricter
conditions on the distribution and covariance matrix of measurement error, but
allow the number of covariates $p$ to grow at a faster rate.

With the conditions above, we establish the following theorem regarding
screening consistency of the proposed SISc procedure.

######  Theorem 3.

Consider the linear EIV model (1) under conditions (C(C11))-(C(C11)) and
(C(C11))-(C(C11)), and the screening set $\hat{\mathcal{Q}}_{\text{SISc}}$
defined in (4) with $d=\lfloor\gamma n\rfloor$ and $\tilde{\beta}_{j}$ given
by (3). Then there exists some
$\tilde{\theta}<1-\tau_{1}-\tau_{2}-\log_{n}(s)$ such that when $\gamma\sim
cn^{-\tilde{\theta}}$ with $c>0$, we have
$\mathbb{P}\left(\hat{\mathcal{Q}}_{\text{SISc}}\supseteq\mathcal{S}\right)=1-O\left(p\exp\left(-Cn\right)\right)$
for some positive constant $C$.

Theorem 3 implies that the SISc procedure can reduce exponentially high
dimension to a relatively small dimension $d=O(n^{1-\tilde{\theta}})<n$ while
retaining all the important covariates. In practice, as suggested by Fan and
Lv, (2008), common choices for $d$ include $d=\lfloor n/\log(n)\rfloor$ or
$d=n-1$. Finally, unlike PMSc, a screening procedure based on ranking
components such as SISc can only achieve full selection consistency if we know
the true number of important covariates $s$ in advance. Since that is rarely
the case in practice, then we do not study the theoretical conditions under
which we can choose $d=s$.

##  3 Simulation studies

###  3.1 Simulation setup

We conducted simulation studies to demonstrate the benefit of the proposed
screening procedures by comparing the performance of two-stage estimators,
which screen variables in the first stage and compute a corrected penalized
estimator on the retained variables in the second stage, against one-stage
estimators, which do not employ any screening method. We simulated data from
model (1), where each row of the true covariate matrix
$\operatorname{\mathbf{X}}_{i}$ was generated from a $p$-variate Gaussian
distribution with a zero mean vector and two choices for the covariance matrix
$\bm{\Sigma}_{x}$. In the first scenario, $\bm{\Sigma}_{x}$ had an
autoregressive AR(1) structure with elements $\sigma_{ij}=\rho_{x}^{|i-j|}$.
In the second scenario, $\bm{\Sigma}_{x}$ had an homogeneous (exchangeable)
structure where all diagonal and off-diagonal elements were set to 1 and
$\rho_{x}$ correspondingly. For both scenarios, we varied
$\rho_{x}\in\\{0.3,0.5\\}$. The true $p$-dimensional vector $\bm{\beta}_{0}$
was constructed such that the first $s=5$ non-zero components were generated
from the uniform distribution $U(1,1.5)$. Turning to the measurement errors,
we generated $\operatorname{\mathbf{U}}_{i}$ independently of
$\operatorname{\mathbf{X}}_{i}$ from another $p$-variate Gaussian distribution
with a zero mean vector and two choices of the covariance matrix
$\bm{\Sigma}_{u}$. For the first choice, $\bm{\Sigma}_{u}$ was a diagonal
matrix with elements randomly generated from the uniform $U(0.1,0.5)$
distribution; hence, the measurement errors on each covariate were independent
of those of the other covariates with the noise-to-signal ratio ranging from
10% to 50%. For the second choice, $\bm{\Sigma}_{u}$ had a block diagonal
structure, where the $p$ covariates were divided into $p/4$ non-overlapping
groups of size $4$ with the correlation between any pair in the same group
equal to $0.2$, and the diagonal elements of $\bm{\Sigma}_{u}$ equal to $0.4$;
hence, the measurement errors on one covariate were positively correlated with
those on three other covariates with the noise-to-signal ratio being $40\%$.
We set the sample size to $n=500$, and varied the number of covariates
$p\in\\{1000,2000\\}$. Finally, the elements of the model error term
$\bm{\varepsilon}$ were independently generated from a Gaussian distribution
with mean zero and variance $\sigma^{2}=0.25.$

###  3.2 Estimators and performance metrics

We computed several one-stage and two-stage estimators on each simulated
dataset. For one-stage estimators, we implemented the corrected lasso
estimator of Loh and Wainwright, (2012) and Sørensen et al., (2015), and the
convex conditioned (CoCo) lasso of Datta et al., (2017). Without using any
screening procedure, the one-stage corrected lasso is given by

$\hat{\bm{\beta}}^{\text{CL}}=\operatorname*{argmin}_{\operatorname{\bm{\beta}}}\left\\{\dfrac{1}{n}\sum_{i=1}^{n}(y_{i}-\operatorname{\mathbf{W}}_{i}^{\top}\bm{\beta})^{2}-\bm{\beta}^{\top}\bm{\Sigma}_{u}\bm{\beta}+\mu\lVert\bm{\beta}\rVert_{1}\right\\},~{}\text{subject
to ~{} }\lVert\bm{\beta}\rVert_{1}\leq R,$ (5)

with $\mu$ and $R$ being two positive tuning parameters. As noted in Loh and
Wainwright, (2012), the problem (5) is generally non-convex when $p>n$,
because the matrix
$n^{-1}\sum_{i=1}^{n}\operatorname{\mathbf{W}}_{i}^{\top}\operatorname{\mathbf{W}}_{i}-\bm{\Sigma}_{u}$
has a large number of negative eigenvalues and is not positive definite.
Therefore, classic gradient descent algorithms, as outlined in Loh and
Wainwright, (2012) and Sørensen et al., (2015), are only guaranteed to
converge to a local minimum under a careful choice of the tuning parameters
$\mu$ and $R$. Sørensen et al., (2015) suggested that the tuning parameter
$\mu$ be chosen as the minimum of the ten-fold cross-validation curve of the
naive lasso

$\hat{\bm{\beta}}^{\text{naive}}=\operatorname*{argmin}_{\bm{\beta}}\left\\{\dfrac{1}{n}\sum_{i=1}^{n}(y_{i}-\operatorname{\mathbf{W}}_{i}^{\top}\bm{\beta})^{2}+\mu\lVert\bm{\beta}\rVert_{1}\right\\},$

and $R$ be chosen by another ten-fold cross-validation from a grid of equally
spaced values between $10^{-3}\kappa$ and $\kappa$, with
$\kappa=2\lVert\hat{\bm{\beta}}^{\text{naive}}\rVert_{1}$. The one-stage CoCo
lasso estimator is given by

$\hat{\bm{\beta}}^{\text{CoCo}}=\underset{\operatorname{\bm{\beta}}}{\operatorname*{argmin}}\left\\{\dfrac{1}{2}\operatorname{\bm{\beta}}^{\prime}\widetilde{\bm{\Sigma}}\operatorname{\bm{\beta}}-\tilde{\rho}^{\top}\operatorname{\bm{\beta}}+\tilde{\mu}\|\operatorname{\bm{\beta}}\|_{1}\right\\},$

where $\tilde{\rho}=\operatorname{\mathbf{W}}^{\top}\bm{y}$ and
$\widetilde{\bm{\Sigma}}=\underset{\bm{\Sigma}\in\mathbb{R}^{p\times
p}_{+}}{\operatorname*{argmin}}\lVert\bm{\Sigma}-\widehat{\bm{\Sigma}}\rVert_{\text{max}},$
with $\widehat{\bm{\Sigma}}=n^{-1}\mathbf{W}^{\top}\mathbf{W}-\bm{\Sigma}_{u}$
and $\mathbb{R}^{p\times p}_{+}$ the set of $p\times p$ positive semi-definite
matrices. In other words, $\widetilde{\bm{\Sigma}}$ is the nearest positive
semi-definite matrix to $\hat{\bm{\Sigma}}$ measured by the element-wise max
norm, which was computed by an alternating direction method of multipliers
(ADMM) algorithm. We point out that this element-wise max norm makes
computation of $\widetilde{\bm{\Sigma}}$ so expensive when $p$ is large that
in our simulation, we only computed the CoCo lasso estimator if the number of
variables was smaller than $1000$. The non-negative tuning parameter
$\tilde{\mu}$ for the CoCo lasso estimator was selected via a corrected cross-
validation procedure; see Datta et al., (2017) and Datta and Zou, (2020) for
more details.

For two-stage estimators, we first implemented either PMSc${}_{\text{CV}}$,
PMSc${}_{\text{FS}}$, or SISc, and then applied either the corrected lasso or
CoCo lasso to the covariates selected from the first stage. For the two
versions of PMSc, we set $\alpha=0.5$ as is commonly done when the bridge
penalty is used in practice (Huang et al.,, 2008, 2009; Polson et al.,, 2014).
For PMSc${}_{\text{CV}}$, we selected the tuning parameter $\lambda_{n}$ via
the following five-fold cross-validation procedure: on each simulated dataset,
the data were randomly split into 5 folds
$\mathcal{F}_{1},\ldots,\mathcal{F}_{5}$ of equal sizes; for the $k$th
iteration, the fold $\mathcal{F}_{k}$ was left out to be the test set and PMSc
was applied on the remaining 4 folds on a grid consisting of $40$ equally
spaced values of the tuning parameter $\lambda_{n}$. Let
$\hat{\operatorname{\bm{\beta}}}_{\lambda}^{(-k)}=(\hat{\beta}_{1,\lambda}^{(-k)},\ldots,\hat{\beta}_{p,\lambda}^{(-k)})^{\top}$
denote the solution of (2) when $\mathcal{F}_{k}$ was left out. Then the final
$\lambda_{n}$ was selected to minimize

$\sum_{k=1}^{5}\left\\{\sum_{i\in\mathcal{F}_{k}}\sum_{j=1}^{p}\left(y_{i}-\hat{\beta}_{j,\lambda}^{(-k)}\tilde{W}_{ij}\right)^{2}-|\mathcal{F}_{k}|\sum_{j=1}^{p}\left(\hat{\beta}_{j,\lambda}^{(-k)}\right)^{2}\sigma_{j}^{2}\right\\},$

where $|\mathcal{F}_{k}|$ denotes the cardinality of the fold
$\mathcal{F}_{k}$. For the PMSc${}_{\text{FS}}$ and SISc, the number of
covariates to be included and kept was chosen to be $M=d=\left\lfloor
n/\log(n)\right\rfloor$, respectively. This choice of $d$ and $M$ follows a
common practice for screening procedures in high dimensional settings (Fan and
Lv,, 2008; Zhu et al.,, 2011; Cui et al.,, 2015). The two-stage estimator when
PMSc${}_{\text{CV}}$ and corrected lasso are used in the first and the second
stage respectively is referred to as the PMSc${}_{\text{CV}}$-Corrected
estimator. Similar definitions hold for other estimators resulting from other
combinations of the methods used in the first and second stage. Also, we refer
to all the estimators that used SISc in the first stage as SISc-based
estimators; similar definitions are applied to the estimators that used
another screening procedure in the first stage.

We report false positive rate (FPR) and false negative rate (FNR) of the
screening procedures, and compare all the one-stage and two-stage estimators
based on FPR, FNR, and $\ell_{2}$ estimation error (i.e the $\ell_{2}$ norm of
the difference between an estimate and the true coefficient vector). All
metrics were averaged across 500 simulations. We also report the mean total
computation time for each estimator. All the reported times are based on
implementation on Artemis, a high performance computing cluster at the
University of Sydney where every sample was computed on a single core
(Intel(R) Xeon(R) CPU E5-2697A v4 2.60GHz and $4$GB RAM).

###  3.3 Simulation results

Tables 1 and 2 present the summary results for the simulations with
$\rho_{x}=0.5$. The results for other settings show similar trend; complete
results, including the standard errors and separation of total computation
time into the first and second stages for all the settings, can be found in
Appendix D.

When only a few true covariates are correlated with one another (i.e
$\bm{\Sigma}_{x}$ had an $\text{AR}(1)$ structure), Table 1 shows that all the
screening procedures (first stage) were able to keep all the important
variables and reduce the number of dimensions effectively. The
PMSc${}_{\text{CV}}$ method had FPR and FNR closest to zero; as a result,
after the second stage, the PMSc${}_{\text{CV}}$-based estimators had the
lowest estimation error. The PMSc${}_{\text{FS}}$ and SISc screening
procedures had zero FNR but some false positives, which were expected from
further investigation of Figure 1 in Appendix C, where we demonstrate that the
choice of $d$ and $M$ in this simulation was greater than the minimum number
of variables that PMSc${}_{\text{FS}}$ and SISc needed to include in order to
retain all the important variables. However, after the second stage, the FPRs
of the two-stage PMSc${}_{\text{FS}}$ and SISc-based estimators were reduced
to close to zero, and notably lower than the FPRs of the corresponding one-
stage estimators; the gain was more pronounced when the CoCo lasso was used in
the second stage. In turn, the estimation errors of the two-stage estimators
were also lower than those of the one-stage estimators, with the improvement
increasing when $p$ increased from $1000$ to $2000$. When $\bm{\Sigma}_{u}$
was diagonal, there was little difference in the $\ell_{2}$ error of the two-
stage estimators; however, when $\bm{\Sigma}_{u}$ was block-diagonal, the
SISc-based estimators tended to perform worse than the PMSc${}_{\text{FS}}$
and PMSc${}_{\text{CV}}$-based estimators. Regarding computation time, among
the three screening methods, SISc was the fastest to compute, followed by
PMSc${}_{\text{FS}}$ and PMSc${}_{\text{CV}}$. As a whole, the total
computation times of the two-stage estimators were considerably lower than
those of the one-stage estimators; further investigation of Table C.1 in
Appendix C shows that the screening stage substantially decreased the time
needed for the second stage estimation. Finally, when $p$ increased from
$p=1000$ to $p=2000$, there was a remarkable increase in the computation time
of the one-stage estimators; however, the corresponding increases for the two-
stage estimators were small, reflecting their computational scalability.

Table 1: Performance of the one-stage and two-stage estimators in the
simulation study based on mean false positive rate (FPR, in percentage), false
negative rate (FNR, in percentage), $\ell_{2}$ error and computation time (in
seconds) when $\bm{\Sigma}_{x}$ has an AR(1) structure with autocorrelation
$\rho_{x}=0.5$. The CoCo estimator was only computed when the number of
variables was no more than $1000$ (either in 1st or 2nd step).

$\bm{\Sigma}_{u}$ | $p$ | Estimator | 1st step | 2nd step | Time
---|---|---|---|---|---
|  |  | FPR | FNR | FPR | FNR | $\ell_{2}$ |
Diagonal | 1000 | One-stage Corrected | - | - | 2.2 | 0.0 | 0.35 | 49.7
|  | PMSc${}_{\text{CV}}$-Corrected | 0.1 | 0.0 | 0.1 | 0.0 | 0.26 | 23.8
|  | PMSc${}_{\text{FS}}$-Corrected | 7.5 | 0.0 | 0.0 | 0.0 | 0.28 | 7.0
|  | SISc-Corrected | 7.5 | 0.0 | 1.2 | 0.0 | 0.33 | 3.5
|  | One-stage CoCo | - | - | 7.6 | 0.0 | 0.47 | 516.8
|  | PMSc${}_{\text{CV}}$-CoCo | 0.1 | 0.0 | 0.1 | 0.0 | 0.26 | 23.4
|  | PMSc${}_{\text{FS}}$-CoCo | 7.5 | 0.0 | 0.0 | 0.0 | 0.27 | 6.1
|  | SISc-CoCo | 7.5 | 0.0 | 0.1 | 0.0 | 0.31 | 2.2
| 2000 | One-stage Corrected | - | - | 1.3 | 0.0 | 0.36 | 136.5
|  | PMSc${}_{\text{CV}}$-Corrected | 0.0 | 0.0 | 0.0 | 0.0 | 0.26 | 38.2
|  | PMSc${}_{\text{FS}}$-Corrected | 3.8 | 0.0 | 0.6 | 0.0 | 0.29 | 20.0
|  | SISc-Corrected | 3.8 | 0.0 | 0.6 | 0.0 | 0.33 | 4.5
|  | PMSc${}_{\text{CV}}$-CoCo | 0.0 | 0.0 | 0.0 | 0.0 | 0.26 | 38.0
|  | PMSc${}_{\text{FS}}$-CoCo | 3.8 | 0.0 | 0.1 | 0.0 | 0.27 | 18.0
|  | SISc-CoCo | 3.8 | 0.0 | 0.1 | 0.0 | 0.31 | 3.5
Block diagonal | 1000 | One-stage Corrected | - | - | 1.8 | 0.0 | 0.62 | 27.8
| PMSc${}_{\text{CV}}$-Corrected | 0.0 | 0.0 | 0.0 | 0.0 | 0.43 | 23.8
|  | PMSc${}_{\text{FS}}$-Corrected | 7.5 | 0.0 | 0.0 | 0.0 | 0.48 | 7.2
|  | SISc-Corrected | 7.5 | 0.0 | 1.0 | 0.0 | 0.57 | 4.1
|  | One-stage CoCo | - | - | 7.4 | 0.0 | 0.74 | 417.4
|  | PMSc${}_{\text{CV}}$-CoCo | 0.0 | 0.0 | 0.0 | 0.0 | 0.43 | 23.4
|  | PMSc${}_{\text{FS}}$-CoCo | 7.5 | 0.0 | 0.0 | 0.0 | 0.46 | 6.2
|  | SISc-CoCo | 7.5 | 0.0 | 0.1 | 0.0 | 0.54 | 3.2
| 2000 | One-stage Corrected | - | - | 1.0 | 0.0 | 0.64 | 68.0
|  | PMSc${}_{\text{CV}}$-Corrected | 0.0 | 0.0 | 0.0 | 0.0 | 0.44 | 38.2
|  | PMSc${}_{\text{FS}}$-Corrected | 3.8 | 0.0 | 0.5 | 0.1 | 0.46 | 19.5
|  | SISc-Corrected | 3.8 | 0.0 | 0.5 | 0.0 | 0.57 | 5.1
|  | PMSc${}_{\text{CV}}$-CoCo | 0.0 | 0.0 | 0.0 | 0.0 | 0.45 | 37.9
|  | PMSc${}_{\text{FS}}$-CoCo | 3.8 | 0.0 | 0.0 | 0.0 | 0.46 | 18.4
|  | SISc-CoCo | 3.8 | 0.0 | 0.1 | 0.0 | 0.54 | 4.2

When all the true covariates were highly correlated with each other (i.e
$\bm{\Sigma}_{x}$ had a homogeneous structure), Table 2 shows that
PMSc${}_{\text{CV}}$ was not helpful in either improving performance or
reducing computation time. For the first step, PMSc${}_{\text{CV}}$ kept all
the variables, meaning no dimension reduction was achieved. As a result, the
two-stage PMSc${}_{\text{CV}}$-based estimators had essentially the same
performance as the corresponding one-stage estimators. Additionally,
PMSc${}_{\text{FS}}$ and SISc had similar and much lower FPRs than
PMSc${}_{\text{CV}}$, with the former tending to have lower FNR than the
latter. As further demonstrated in Web Figure 1 of Appendix C, these FNRs were
expected because the minimum number of variables SISc had to keep in order to
retain all the important variables was much higher than that of
PMSc${}_{\text{FS}}$. After the second stage, the FPRs of all the PMSc-based
(both versions) and SISc-based estimators were much reduced, while the
corresponding FNRs were unchanged. Regarding estimation error, the two-stage
PMSc${}_{\text{FS}}$ and SISc-based estimators had noticeably lower estimation
errors than the corresponding one-stage estimators, with the improvement being
greater when $p$ increased from $1000$ to $2000$. When $\bm{\Sigma}_{u}$ was
diagonal, the improvement was considerable for both the corrected and CoCo
lasso. In contrast, when $\bm{\Sigma}_{u}$ was block-diagonal, the most
remarkable improvement was made for the CoCo lasso only. In all the settings,
the PMSc${}_{\text{FS}}$-CoCo estimators had the smallest $\ell_{2}$
estimation error. Finally, regarding computation time, it was not surprising
that the SISc was still computationally the fastest screening method. The two-
stage PMSc${}_{\text{FS}}$-based estimators were also fast to compute, and
when $p$ increased to $2000$, the increase in computation time for PMSc and
SISc was small compared to that of the one-stage estimators.

Table 2: Performance of the one-stage and two-stage estimators in the
simulation study based on false positive rate (FPR, in percentage), false
negative rate (FNR, in percentage) $\ell_{2}$ error, and computation time (in
seconds) when $\bm{\Sigma}_{x}$ has an homogeneous structure with
$\rho_{x}=0.5$. The CoCo estimator was only computed when the number of
variables was no more than $1000$ (either in 1st or 2nd step).

$\bm{\Sigma}_{u}$ | $p$ | Estimator | 1st step | 2nd step | Time
---|---|---|---|---|---
|  |  | FPR | FNR | FPR | FNR | $\ell_{2}$ |
Diagonal | 1000 | One-stage Corrected | - | - | 20.0 | 0.1 | 0.96 | 18.9
|  | PMSc${}_{\text{CV}}$-Corrected | 100.0 | 0.0 | 19.7 | 0.2 | 0.95 | 31.4
|  | PMSc${}_{\text{FS}}$-Corrected | 7.6 | 3.2 | 1.6 | 3.2 | 0.56 | 6.8
|  | SISc-Corrected | 7.6 | 6.7 | 1.0 | 6.7 | 0.54 | 2.2
|  | One-stage CoCo | - | - | 11.5 | 0.0 | 0.77 | 481.5
|  | PMSc${}_{\text{CV}}$-CoCo | 100.0 | 0.0 | 11.5 | 0.0 | 0.77 | 489.4
|  | PMSc${}_{\text{FS}}$-CoCo | 7.6 | 3.2 | 1.0 | 3.2 | 0.45 | 7.2
|  | SISc-CoCo | 7.6 | 6.7 | 1.0 | 6.7 | 0.52 | 2.9
| 2000 | One-stage Corrected | - | - | 29.1 | 1.1 | 1.57 | 49.9
|  | PMSc${}_{\text{CV}}$-Corrected | 100.0 | 0.0 | 28.3 | 1.0 | 1.53 | 74.5
|  | PMSc${}_{\text{FS}}$-Corrected | 3.8 | 4.0 | 0.7 | 4.0 | 0.62 | 9.6
|  | SISc-Corrected | 3.8 | 11.0 | 0.5 | 11.0 | 0.61 | 3.6
|  | PMSc${}_{\text{CV}}$-CoCo | 100.0 | 0.0 | - | - | - | -
|  | PMSc${}_{\text{FS}}$-CoCo | 3.8 | 4.0 | 0.5 | 4.0 | 0.45 | 10.2
|  | SISc-CoCo | 3.8 | 11.0 | 0.5 | 11.0 | 0.60 | 4.6
Block diagonal | 1000 | One-stage Corrected | - | - | 23.5 | 1.1 | 1.25 | 22.7
| PMSc${}_{\text{CV}}$-Corrected | 100.0 | 0.0 | 25.6 | 0.8 | 1.30 | 37.0
|  | PMSc${}_{\text{FS}}$-Corrected | 7.6 | 3.4 | 3.0 | 3.6 | 1.15 | 6.7
|  | SISc-Corrected | 7.6 | 14.2 | 1.5 | 14.2 | 0.97 | 2.4
|  | One-stage CoCo | - | - | 14.3 | 0.4 | 1.23 | 417.2
|  | PMSc${}_{\text{CV}}$-CoCo | 100.0 | 0.0 | 14.3 | 0.5 | 1.23 | 427.9
|  | PMSc${}_{\text{FS}}$-CoCo | 7.6 | 3.4 | 0.9 | 3.6 | 0.68 | 8.2
|  | SISc-CoCo | 7.6 | 14.2 | 0.9 | 14.2 | 0.84 | 3.6
| 2000 | One-stage Corrected | - | - | 32.4 | 3.0 | 1.86 | 52.5
|  | PMSc${}_{\text{CV}}$-Corrected | 100.0 | 0.0 | 31.8 | 3.3 | 1.84 | 79.8
|  | PMSc${}_{\text{FS}}$-Corrected | 3.8 | 6.6 | 1.6 | 6.6 | 1.21 | 9.4
|  | SISc-Corrected | 3.8 | 21.4 | 0.8 | 21.4 | 1.09 | 3.8
|  | PMSc${}_{\text{CV}}$-CoCo | 100.0 | 0.0 | - | - | - | -
|  | PMSc${}_{\text{FS}}$-CoCo | 3.8 | 6.6 | 0.5 | 6.7 | 0.73 | 10.8
|  | SISc-CoCo | 3.8 | 21.4 | 0.5 | 21.4 | 0.97 | 5.3

In summary, the simulation study both confirms the theoretical results for the
screening procedures when the true covariates are not highly correlated, and
demonstrates the superior finite sample performance, computational efficiency
and scalability of the proposed two-stage estimators compared to the one-stage
estimators. Among the three screening methods implemented in the simulations,
the PMSc${}_{\text{FS}}$ was the most reliable and efficient screening method,
taking account of both performance metrics and computation time when the true
covariates were either moderately or highly correlated.

##  4 Analysis of microarray data

We applied the proposed methodology to analyze the motivating Bone Mineral
Density Data consisting of gene expression measurements of $54,675$ probe sets
and bone mineral density (BMD) for $n=84$ Norwegian women. The dataset is
publicly available at the European Bioinformatics Institute ArrayExpress
repository under the access number E-MEXP-1618. Microarray measurements are
known to be noisy (Rocke and Durbin,, 2001); although biological variation in
the data is usually of primary interest to investigators, it can be obscured
by measurement errors coming from many sources (Zakharkin et al.,, 2005).
Furthermore, a distinctive feature of the Affymetrix microarray dataset is
that multiple probes are used to measure each gene expression; hence, these
replicated measurements can be used to estimate the covariance matrix of
measurement errors $\bm{\Sigma}_{u}$.

We followed the procedures described in Sørensen et al., (2015), Nghiem and
Potgieter, (2019) and Romeo and Thoresen, (2019) to process the raw data using
the BGX package of Hein et al., (2005), and assumed the measurement error on
each gene was mutually independent from that on the other. As a result, the
measurement error covariance matrix $\bm{\Sigma}_{u}$ was set to be diagonal.
Similar to Sørensen et al., (2015), we also kept $p=993$ genes that had
between-patient variability greater than measurement error variance, i.e the
noise-to-signal ratio was smaller than $1$. After processing, we obtained the
surrogate matrix $\operatorname{\mathbf{W}}$ along with $\bm{\Sigma}_{u}$.
More details on the raw data processing and the estimated covariance matrix of
measurement error can be found in Appendix B. The response variable was chosen
to be the (centered) total hip T-score.

We computed the eight estimators as in the simulation study. For the screening
step, the PMSc${}_{\text{CV}}$ procedure selected $610$ genes, so it did not
reduce the number of dimensions by any great amount. This was likely due to
the high correlation among the gene expressions, similar to the case when
$\bm{\Sigma}_{x}$ had a homogeneous structure in the simulation. Table 3 shows
the number of genes selected by one-stage and two-stage estimators, as well as
the $\ell_{2}$ norm of the estimated coefficients.

Table 3: The number of selected genes and the corresponding $\ell_{2}$ norm of the estimated coefficients obtained from each estimator in the microarray data analysis. Computation time is based on implementation on a laptop with one Dual-Core Intel Core i5 2.7GHz processor and 8GB RAM Estimator | # of selected genes |  $\ell_{2}$ norm | Time (in seconds)
---|---|---|---
| 1st step | 2nd step |  | 1st step | 2nd step
One-stage Corrected | - | 108 | 0.33 | 0.00 | 85.97
PMSc${}_{\text{CV}}$-Corrected | 798 | 104 | 0.34 | 74.02 | 43.48
PMSc${}_{\text{FS}}$-Corrected | 19 | 4 | 2.10 | 6.54 | 2.49
SISc-Corrected | 19 | 4 | 2.00 | 0.73 | 2.49
One-stage CoCo | - | 18 | 0.76 | 0.00 | 2702.81
PMSc${}_{\text{CV}}$-CoCo | 798 | 17 | 0.78 | 74.02 | 1648.86
PMSc${}_{\text{FS}}$-CoCo | 19 | 6 | 1.05 | 6.54 | 0.21
SISc-CoCo | 19 | 5 | 1.12 | 0.73 | 0.30

It can be seen that while the one-stage corrected lasso and
PMSc${}_{\text{CV}}$-Corrected estimators select many more genes than all the
other estimators, the $\ell_{2}$ norms of the corresponding estimated
coefficients were smaller or equivalent to that of the other estimators;
hence, the one-stage corrected lasso and PMSc${}_{\text{CV}}$-Corrected
estimators likely contained many false positives. Among the $p=993$ genes,
there were 15 genes selected by one estimator, 113 genes selected by two
estimators, 1 gene selected by 3 estimators, 1 gene selected by 4 estimators,
and notably 3 genes selected by 6 estimators. Table 4 demonstrates that the
estimated coefficients associated with these three genes had larger magnitude
than those of other selected genes. Furthermore, compared to the one-stage
CoCo lasso and PMSc${}_{\text{CV}}$-CoCo estimators, the two-stage
PMSc${}_{\text{FS}}$-CoCo and SISc-CoCo estimators magnified the effects of
these three genes. These estimators might contain some false negatives, but
they help practitioners obtain stronger signals of the important variables
that are present, as has been noted in the literature of measurement error
correction in high dimensional settings; see Sørensen et al., (2015) and
Nghiem and Potgieter, (2019). Finally, regarding computation time, the two-
stage estimators were generally faster to compute, with the exception of the
PMSc${}_{\text{CV}}$-based estimators; the gain in computation time was the
largest when the CoCo estimator was used in the second stage.

Table 4: The estimated coefficients and their corresponding rank in terms of magnitude for the three genes most frequently selected by most estimators. Estimator | Estimate | Rank (1 = largest)
---|---|---
| Gene 1 | Gene 2 | Gene 3 | Gene 1 | Gene 2 | Gene 3
PMSc${}_{\text{FS}}$-Corrected | 0.96 | -1.54 | 0.62 | 2 | 1 | 4
SISc-Corrected | 0.86 | -1.43 | 0.63 | 2 | 1 | 4
One-stage CoCo | 0.61 | -0.25 | 0.27 | 1 | 3 | 2
PMSc${}_{\text{CV}}$-CoCo | 0.62 | -0.28 | 0.29 | 1 | 3 | 2
PMSc${}_{\text{FS}}$-CoCo | 0.78 | -0.47 | 0.29 | 1 | 2 | 4
SISc-CoCo | 0.49 | -0.94 | 0.18 | 2 | 1 | 4

##  5 Conclusion

This paper proposes two screening procedures for linear errors-in-variables
models in high dimensional settings, namely the corrected penalized marginal
regression and the corrected sure independence screening procedures. Both
procedures are based on fitting corrected marginal regressions of the outcome
on each contaminated covariate, which could be computed efficiently in high
dimensions. Under mild technical conditions, these procedures are shown to
achieve screening consistency, meaning all the important variables are fully
retained. Under a stronger condition of partial orthogonality among the true
covariates, we further illustrate that the corrected penalized marginal
screening approach (using the bridge penalty) can achieve full selection
consistency. We demonstrated the advantages of these screening procedures in
practice through a simulation study and an analysis of a microarray data
concerning gene expressions of bone mineral density for Norweigian women, both
of which showed that the proposed screening procedures reduced computation
time considerably and/or improved performance metrics for both estimation and
variable selection.

Future research could aim at reducing the false negative rates for screening
procedures when the true covariates are highly collinear; for such a setting,
an iterative corrected marginal screening procedure similar to the iterative
sure independent screening of Fan and Lv, (2008) may be considered.
Furthermore, the screening procedures presented in this paper can be extended
to more complicated errors-in-variables models, such as generalized linear
models and non-parametric regression settings, although it may be more
challenging to correct for measurement errors in marginal regressions of the
outcome on each contaminated covariate in these models.

## 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.
  * Belloni et al., (2017) Belloni, A., Rosenbaum, M., and Tsybakov, A. B. (2017).  Linear and conic programming estimators in high dimensional errors-in-variables models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79(3):939–956.
  * Brown et al., (2019) Brown, B., Weaver, T., and Wolfson, J. (2019).  Meboost: Variable selection in the presence of measurement error.  Statistics in Medicine, 38(15):2705–2718.
  * Byrd and McGee, (2019) Byrd, M. and McGee, M. (2019).  A simple correction procedure for high-dimensional general linear models with measurement error.  arXiv preprint arXiv:1912.11740.
  * Carroll et al., (2006) Carroll, R. J., Ruppert, D., Stefanski, L. A., and Crainiceanu, C. M. (2006).  Measurement Error in Nonlinear Models: A Modern Perspective.  CRC press.
  * 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.
  * Datta and Zou, (2020) Datta, A. and Zou, H. (2020).  A note on cross-validation for lasso under measurement errors.  Technometrics, 62(4):549–556.
  * Datta et al., (2017) Datta, A., Zou, H., et al. (2017).  Cocolasso for high-dimensional error-in-variables regression.  The Annals of Statistics, 45(6):2400–2426.
  * Do et al., (2011) Do, C. B., Tung, J. Y., Dorfman, E., Kiefer, A. K., Drabant, E. M., Francke, U., Mountain, J. L., Goldman, S. M., Tanner, C. M., Langston, J. W., et al. (2011).  Web-based genome-wide association study identifies two novel loci and a substantial genetic component for parkinson’s disease.  PLoS Genet, 7(6):e1002141.
  * 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., (2004) Fan, J., Peng, H., et al. (2004).  Nonconcave penalized likelihood with a diverging number of parameters.  The Annals of Statistics, 32(3):928–961.
  * Fan et al., (2010) Fan, J., Song, R., et al. (2010).  Sure independence screening in generalized linear models with np-dimensionality.  The Annals of Statistics, 38(6):3567–3604.
  * Frank and Friedman, (1993) Frank, L. E. and Friedman, J. H. (1993).  A statistical view of some chemometrics regression tools.  Technometrics, 35(2):109–135.
  * Hein et al., (2005) Hein, A.-M. K., Richardson, S., Causton, H. C., Ambler, G. K., and Green, P. J. (2005).  Bgx: A fully Bayesian integrated approach to the analysis of Affymetrix Genechip data.  Biostatistics, 6(3):349–373.
  * Huang et al., (2008) Huang, J., Horowitz, J. L., Ma, S., et al. (2008).  Asymptotic properties of bridge estimators in sparse high-dimensional regression models.  The Annals of Statistics, 36(2):587–613.
  * Huang et al., (2009) Huang, J., Ma, S., Xie, H., and Zhang, C.-H. (2009).  A group bridge approach for variable selection.  Biometrika, 96(2):339–355.
  * Ida et al., (2019) Ida, Y., Fujiwara, Y., and Kashima, H. (2019).  Fast sparse group lasso.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 32, pages 1702–1710.
  * Kaul et al., (2016) Kaul, A., Koul, H. L., Chawla, A., and Lahiri, S. N. (2016).  Two stage non-penalized corrected least squares for high dimensional linear models with measurement error or missing covariates.  arXiv preprint arXiv:1605.03154.
  * Li et al., (2012) Li, G., Peng, H., Zhang, J., Zhu, L., et al. (2012).  Robust rank correlation based screening.  The Annals of Statistics, 40(3):1846–1877.
  * Loh and Wainwright, (2012) Loh, P.-L. and Wainwright, M. J. (2012).  High-dimensional regression with noisy and missing data: Provable guarantees with nonconvexity.  The Annals of Statistics, 40(3):1637–1664.
  * Nghiem and Potgieter, (2019) Nghiem, L. and Potgieter, C. (2019).  Simulation-selection-extrapolation: Estimation in high-dimensional errors-in-variables models.  Biometrics, 75(0):1133–1144.
  * Piironen et al., (2017) Piironen, J., Vehtari, A., et al. (2017).  Sparsity information and regularization in the horseshoe and other shrinkage priors.  Electronic Journal of Statistics, 11(2):5018–5051.
  * Polson et al., (2014) Polson, N. G., Scott, J. G., and Windle, J. (2014).  The Bayesian bridge.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 76:713–733.
  * Rocke and Durbin, (2001) Rocke, D. M. and Durbin, B. (2001).  A model for measurement error for gene expression arrays.  Journal of Computational Biology, 8(6):557–569.
  * Romeo and Thoresen, (2019) Romeo, G. and Thoresen, M. (2019).  Model selection in high-dimensional noisy data: a simulation study.  Journal of Statistical Computation and Simulation, 89(11):2031–2050.
  * Rosenbaum et al., (2010) Rosenbaum, M., Tsybakov, A. B., et al. (2010).  Sparse recovery under matrix uncertainty.  The Annals of Statistics, 38(5):2620–2651.
  * Rosenbaum et al., (2013) Rosenbaum, M., Tsybakov, A. B., et al. (2013).  Improved matrix uncertainty selector.  In From Probability to Statistics and Back: High-Dimensional Models and Processes–A Festschrift in Honor of Jon A. Wellner, pages 276–290. Institute of Mathematical Statistics.
  * Simon et al., (2013) Simon, N., Friedman, J., Hastie, T., and Tibshirani, R. (2013).  A sparse-group lasso.  Journal of Computational and Graphical Statistics, 22(2):231–245.
  * Sørensen et al., (2015) Sørensen, Ø., Frigessi, A., and Thoresen, M. (2015).  Measurement error in Lasso: Impact and likelihood bias correction.  Statistica Sinica, 25:809–829.
  * Tibshirani, (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 58(1):267–288.
  * Wainwright, (2019) Wainwright, M. J. (2019).  High-dimensional Statistics: A Non-Asymptotic Viewpoint, volume 48.  Cambridge University Press.
  * 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:293–317.
  * Zakharkin et al., (2005) Zakharkin, S. O., Kim, K., Mehta, T., Chen, L., Barnes, S., Scheirer, K. E., Parrish, R. S., Allison, D. B., and Page, G. P. (2005).  Sources of variation in Affymetrix microarray experiments.  BMC bioinformatics, 6(1):1–11.
  * Zhou et al., (2018) Zhou, T., Thung, K.-H., Liu, M., and Shen, D. (2018).  Brain-wide genome-wide association study for alzheimer’s disease via joint projection learning and sparse regression model.  IEEE Transactions on Biomedical Engineering, 66(1):165–175.
  * 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.

Appendix

##  A Technical proofs

We first recall the technical conditions that are given in the main paper.

  1. (C1)

The model error terms $\varepsilon_{1},\varepsilon_{2},\ldots\varepsilon_{n}$
are iid sub-Gaussian random variables with mean zero and finite variance
$\sigma^{2}$, i.e there exists a finite constant $\sigma_{*}>0$ such that for
all $t\in\mathbb{R}$,

$\mathbb{E}\left\\{\exp\left(t\varepsilon_{i}\right)\right\\}\leq\exp\left(\frac{\sigma_{*}^{2}t^{2}}{2}\right),~{}i=1,\ldots,n,$

where the constant $\sigma_{*}^{2}$ is referred to as the “variance proxy“ for
$\varepsilon_{i}$.

  2. (C2)

The measurement errors $U_{ij},~{}i=1,\ldots,n$ are independent sub-Gaussian
random variables with variance proxy $\sigma_{*}^{2}$, $j=1,\ldots,p$.

  3. (C3)

There exist constants $b_{0}$ and $b_{1}$ such that $\displaystyle
0<b_{0}\leq\min_{k\in\mathcal{S}}\left|\beta_{0k}\right|\leq\max_{k\in\mathcal{S}}\left|\beta_{0k}\right|\leq
b_{1}\leq\infty.$

  4. (C4)

There exist constants $C_{1}$ and $C_{2}$ such that $\displaystyle
0<C_{1}\leq\min_{j}n^{-1}\sum_{i=1}^{n}X_{ij}^{2}<\max_{j}n^{-1}\sum_{i=1}^{n}X_{ij}^{2}\leq
C_{2}<\infty.$

  5. (C5)

There exists a constant $\xi_{0}$ such that
$\displaystyle\min_{k\in\mathcal{S}}|\tilde{\xi}_{nj}|>\xi_{0}>0$ for all $n$.

  6. (C6)

There exists constants $d_{0}>0$ and $0\leq\theta\leq 1/2$ such that

$\left|n^{-1}\sum_{i=1}^{n}X_{ij}X_{ik}\right|\leq d_{0}n^{-\theta},\quad
j\in\mathcal{S}^{c},~{}k\in\mathcal{S},$

  7. (C7)

Assume the tuning parameter $\lambda_{n}$ for PMSc satisfies:

    1. (a)

$\lambda_{n}\to 0$ and
$\lambda_{n}n^{\theta(2-\alpha)}s^{\alpha-2}\rightarrow\infty$,

    2. (b)

$\log\left(2m\right)=o(1)\times\left\\{\lambda_{n}n^{(2-\alpha)/2}\right\\}^{1/(2-\alpha)},$
where $m=p-s$.

  8. (C8)

$p>n$, $\log(p)/n\to C_{1}$ with $C_{1}$ being a constant as both $n$ and $p$
diverge to $\infty$.

  9. (C9)

Let $\mathbf{Z}=\operatorname{\mathbf{U}}\bm{\Sigma}_{u}^{-1/2}$. Then

    1. (a)

For some constants $c_{1},c_{2}>1$ and $D>0$, the matrix $\mathbf{Z}$ follows
a spherical distribution satisfying

$\mathbb{P}\left(\lambda_{\max}\left(\tilde{p}^{-1}\tilde{\mathbf{Z}}\widetilde{\mathbf{Z}}^{\top}\right)>c_{2}\text{
and
}\lambda_{\min}\left(\widetilde{p}^{-1}\tilde{\mathbf{Z}}\tilde{\mathbf{Z}}^{\top}\right)<1/c_{2}\right)\leq
e^{-Dn},$ (6)

for any $n\times\tilde{p}$ submatrix $\tilde{\mathbf{Z}}$ of $\mathbf{Z}$ with
$c_{1}n\leq\tilde{p}\leq p$.

    2. (b)

There exist some constants $\tau_{1}>0$ and $c_{3}>0$, such that
$\lambda_{\max}(\bm{\Sigma}_{u})\leq c_{3}n^{\tau_{1}}.$

  10. (C10)

There exists positive constants $c_{4}>0$ and $\tau_{2}>0$, such that
$\lambda_{\text{max}}\left(n^{-1}\operatorname{\mathbf{X}}^{\top}\operatorname{\mathbf{X}}\right)\leq
c_{4}n^{\tau_{2}}.$

  11. (C11)

$\tau_{1}+\tau_{2}+\log_{n}(s)<1$, where $\log_{n}(s)$ denotes the logarithm
base $n$ of $s$.

###  A.1 Proof of Theorem 1

First, we state the following lemma from Huang et al. (2008) that will be used
in the proof.

######  Lemma 1.

Let $g(u)=u^{2}-2au+\lambda|u|^{\alpha}$ where $a\neq 0$, $\lambda\geq 0$ and
$0<\alpha<1$. Define

$c_{\alpha}=\left(\frac{2}{2-\alpha}\right)\left\\{\frac{2(1-\alpha)}{2-\alpha}\right\\}^{1-\alpha}.$

Then $\arg\min(g)=0$ if and only if $\lambda>c_{\alpha}|a|^{2-\alpha}$.

Additionally, we state the following properties of sub-Gaussian random
variables from Wainwright, (2019):

  1. (B1)

Let $S_{1},\ldots,S_{n}$ be zero-mean independent sub-Gaussian random
variables with variance proxy $\sigma_{0}^{2}$. Denote
$\bm{v}=(v_{1},\ldots,v_{n})\in{R}^{n}$. Then $\sum_{i=1}^{n}v_{i}S_{i}$ is
sub-Gaussian with variance bounded by
$\left(\sum_{i=1}^{n}v_{i}^{2}\right)\sigma_{0}^{2}$.

  2. (B2)

Let $S_{1},\ldots,S_{n}$ be zero-mean sub-Gaussian random variables (not
necessarily independent) with common variance proxy $\sigma_{0}^{2}$. Then

$\mathbb{E}\left(\max_{i=1,\ldots,n}|S_{i}|\right)\leq\sigma_{0}\sqrt{2\log(2n)}.$

  3. (B3)

If $S_{1}$ and $S_{2}$ are zero-mean independent sub-Gaussian random variables
with finite variance proxies $\sigma_{01}^{2}$ and $\sigma_{02}^{2}$
respectively, then $V=S_{1}S_{2}$ is sub-exponential with (finite) parameter
$(\nu_{0},t_{0})$, meaning that
$\mathbb{E}(\exp(sV))\leq\exp(\frac{1}{2}\nu_{0}^{2}s^{2})$ for all
$|s|>t_{0}.$

  4. (B4)

Let $S_{1},\ldots,S_{n}$ be sub-exponential random variables (not necessarily
independent) with common parameter $(\nu_{0},t_{0})$. Then

$\mathbb{E}\left(\max_{i=1,\ldots,n}|S_{i}|\right)\leq\nu_{0}\sqrt{2\log(2n)}+t_{0}\log(2n).$

###### Proof.

Recall that,

$\tilde{\xi}_{nj}=n^{-1}\sum_{i=1}^{n}\left(\sum_{k=1}^{s}X_{ik}\bm{\beta}_{0k}\right)X_{ij},\quad\xi_{nj}=\tilde{\xi}_{nj}+n^{-1}\sum_{i=1}^{n}\left(\sum_{k=1}^{s}X_{ik}\bm{\beta}_{0k}\right)U_{ij},$

so we have

$\displaystyle L(\bm{\beta})$
$\displaystyle=\dfrac{1}{n}\sum_{j=1}^{p}\sum_{i=1}^{n}\left(y_{i}-W_{ij}\beta_{j}\right)^{2}-\sum_{j=1}^{p}\sigma_{j}^{2}\beta_{j}^{2}+\lambda_{n}\sum_{j=1}^{p}\left|\beta_{j}\right|^{\alpha}$
$\displaystyle=\sum_{j=1}^{p}\sum_{i=1}^{n}\dfrac{1}{n}\left(\sum_{k=1}^{s}X_{ik}\beta_{0k}+\varepsilon_{i}-W_{ij}\beta_{j}\right)^{2}-\sum_{j=1}^{p}\sigma_{j}^{2}\beta_{j}^{2}+\lambda_{n}\sum_{j=1}^{p}\left|\beta_{j}\right|^{\alpha}$
$\displaystyle=\sum_{j=1}^{p}\left[\dfrac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}^{2}+\left(\dfrac{1}{n}\sum_{i=1}^{n}W_{ij}^{2}-\sigma_{j}^{2}\right)\beta_{j}^{2}-\frac{2}{n}\left(\bm{\varepsilon}^{\top}\mathbf{a}_{j}+n\xi_{nj}\right)\beta_{j}+\frac{1}{n}\sum_{i=1}^{n}\left(\sum_{k=1}^{s}X_{ik}\beta_{0k}\right)^{2}+\lambda_{n}\left|\beta_{j}\right|^{\alpha}\right]$

where
$\textbf{a}_{j}=(W_{1j},\ldots,W_{nj})^{\top}=(X_{1j},\ldots,X_{nj})^{\top}+(U_{1j},\ldots,U_{nj})^{\top}=\operatorname{\mathbf{X}}_{j}+\operatorname{\mathbf{U}}_{j}$.
We can ignore the terms that do not contain $\bm{\beta}$, so minimizing
$L(\bm{\beta})$ is equivalent to minimizing $\sum_{j=1}^{p}h_{j}(\beta_{j})$,
where

$h_{j}(\beta_{j})=\left(\dfrac{1}{n}\sum_{i=1}^{n}W_{ij}^{2}-\sigma_{j}^{2}\right)\beta_{j}^{2}-\frac{2}{n}\left(\bm{\varepsilon}^{\top}\mathbf{a}_{j}+n\xi_{nj}\right)\beta_{j}+\lambda_{n}\left|\beta_{j}\right|^{\alpha},\quad
j=1,\ldots,p.$

To simplify the notation, let
$V_{j}=n^{-1}\sum_{i=1}^{n}W_{ij}^{2}-\sigma_{j}^{2}$. By Lemma 1,
$\hat{\beta}_{j}=0$ is the only solution of
$\arg\min_{\beta_{j}}h_{j}(\beta_{j})$ if and only if

$V_{j}^{-1}\lambda_{n}>c_{\alpha}\left\\{\left(nV_{j}\right)^{-1}\left|\bm{\varepsilon}^{\top}\mathbf{a}_{j}+n\xi_{nj}\right|\right\\}^{2-\alpha}.$

By some algebra, letting
$\phi_{n}=c_{\alpha}^{-1/(2-\alpha)}\lambda_{n}^{1/(2-\alpha)}n^{1/2}$, the
above inequality is equivalent to

$\phi_{n}>n^{-1/2}V_{j}^{(-1+\alpha)/(2-\alpha)}\left|\bm{\varepsilon}^{\top}\mathbf{a}_{j}+n\xi_{nj}\right|.$

Hence, Theorem 1 will follow if we can prove that

$\lim_{n\to\infty}\mathbb{P}\left\\{\phi_{n}>n^{-1/2}\min_{j\in\mathcal{S}}V_{j}^{(-1+\alpha)/(2-\alpha)}\left|\bm{\varepsilon}^{\top}\mathbf{a}_{j}+n\xi_{nj}\right|\right\\}\rightarrow
0.$ (7)

In order to prove (7), note that as $n\to\infty$, for every $j=1,\ldots,p$, we
have

$\displaystyle\left|\left(\frac{1}{n}\sum_{i=1}^{n}W_{ij}^{2}-\sigma_{j}^{2}\right)-\dfrac{1}{n}\sum_{i=1}^{n}X_{ij}^{2}\right|$
$\displaystyle=\left|\frac{1}{n}\sum_{i=1}^{n}\left(X_{ij}+U_{ij}\right)^{2}-\sigma_{j}^{2}-\dfrac{1}{n}\sum_{i=1}^{n}X_{ij}^{2}\right|$
$\displaystyle\leq\left|\dfrac{1}{n}\left(\sum_{j=1}^{n}U_{ij}^{2}\right)-\sigma_{j}^{2}\right|+\left|\dfrac{2}{n}\sum_{i=1}^{n}X_{ij}U_{ij}\right|=o_{p}(1),$

by the law of large numbers and that
$E(U_{ij})=0,E(U_{ij}^{2})=\sigma_{j}^{2}$ and $E(U_{ij}^{4})<\infty$. Let
$\tau_{j}=(n^{-1}\sum_{i=1}^{n}X_{ij}^{2})^{(-1+\alpha)/(2-\alpha)}$. By
condition (C(C11)) and because $0<\alpha<1$, there exist constants $C_{3}$ and
$C_{4}$ such that $C_{3}\leq\min_{j}\tau_{j}\leq\max_{j}\tau_{j}\leq
C_{4}<\infty$. It suffices to show that

$\mathbb{P}\left\\{\phi_{n}>n^{-1/2}\min_{j\in\mathcal{S}}\tau_{j}\left|\bm{\varepsilon}^{\top}\mathbf{a}_{j}+n\xi_{nj}\right|\right\\}\rightarrow
0.$ (8)

Now, to prove (8), we have

$\displaystyle\mathbb{P}\left(\phi_{n}>\min_{j\in\mathcal{S}}\tau_{j}\left|n^{-1/2}\bm{\varepsilon}^{\top}\mathbf{a}_{j}+n^{1/2}\xi_{nj}\right|\right)$
$\displaystyle=\mathbb{P}\left(\bigcup_{j\in\mathcal{S}}\left\\{\tau_{j}\left|n^{-1/2}\bm{\varepsilon}^{\top}\mathbf{a}_{j}+n^{1/2}\xi_{nj}\right|<\phi_{n}\right\\}\right)$
(9)
$\displaystyle\leq\sum_{j\in\mathcal{S}}\mathbb{P}\left(\tau_{j}\left|n^{-1/2}\bm{\varepsilon}^{\top}\mathbf{a}_{j}+n^{1/2}\xi_{nj}\right|<\phi_{n}\right)$
$\displaystyle\leq\sum_{j\in\mathcal{S}}\mathbb{P}\left(\left|n^{-1/2}\bm{\varepsilon}^{\top}\mathbf{a}_{j}+n^{1/2}\xi_{nj}\right|<C_{3}^{-1}\phi_{n}\right)$

where the last inequality follows from $\tau_{j}\geq C_{3}$. Then noting that
for all $j\in\mathcal{S}$, we have

$\mathbb{P}\left(\left|n^{-1/2}\bm{\varepsilon}^{\top}\mathbf{a}_{j}+n^{1/2}\xi_{nj}\right|<C_{3}^{-1}\phi_{n}\right)=1-\mathbb{P}\left(\left|n^{-1/2}\bm{\varepsilon}^{\top}\mathbf{a}_{j}+n^{1/2}\xi_{nj}\right|\geq
C_{3}^{-1}\phi_{n}\right),$ (10)

and hence

$\displaystyle\mathbb{P}\left(\left|n^{-1/2}\bm{\varepsilon}^{\top}\mathbf{a}_{j}+n^{1/2}\xi_{nj}\right|\geq
C_{3}^{-1}\phi_{n}\right)\geq\mathbb{P}\left(n^{1/2}\left|\xi_{nj}\right|-n^{-1/2}\left|\bm{\varepsilon}^{\top}\mathbf{a}_{j}\right|\geq
C_{3}^{-1}\phi_{n}\right)$
$\displaystyle=1-\mathbb{P}\left(n^{-1/2}\left|\bm{\varepsilon}^{\top}\mathbf{a}_{j}\right|>n^{1/2}\left|\xi_{nj}\right|-C_{3}^{-1}\phi_{n}\right)$
$\displaystyle\geq
1-\mathbb{P}\left(n^{-1/2}\left|\bm{\varepsilon}^{\top}\mathbf{a}_{j}\right|>n^{1/2}|\tilde{\xi}_{nj}|-n^{-1/2}\left|\sum_{i=1}^{n}\left(\sum_{k=1}^{s}X_{ik}\beta_{0k}U_{ij}\right)\right|-C_{3}^{-1}\phi_{n}\right)$
$\displaystyle\geq
1-\mathbb{P}\left(n^{-1/2}\left|\bm{\varepsilon}^{\top}\mathbf{X}_{j}\right|+n^{-1/2}\left|\bm{\varepsilon}^{\top}\mathbf{U}_{j}\right|+n^{-1/2}\left|\sum_{i=1}^{n}\left(\sum_{k=1}^{s}X_{ik}\beta_{0k}U_{ij}\right)\right|>n^{1/2}|\xi_{0}|-C_{3}^{-1}\phi_{n}\right)$
$\displaystyle\geq
1-\mathbb{P}\left(\left|n^{-1/2}\bm{\varepsilon}^{\top}\mathbf{X}_{j}+n^{-1/2}\bm{\varepsilon}^{\top}\mathbf{U}_{j}+n^{-1/2}\sum_{i=1}^{n}\left(\sum_{k=1}^{s}X_{ik}\beta_{0k}U_{ij}\right)\right|>n^{1/2}|\xi_{0}|-C_{3}^{-1}\phi_{n}\right)$
$\displaystyle\geq
1-\frac{\mathbb{E}\left\\{n^{-1/2}\bm{\varepsilon}^{\top}\mathbf{X}_{j}+n^{-1/2}\bm{\varepsilon}^{\top}\mathbf{U}_{j}+n^{-1/2}b_{1}\sum_{i=1}^{n}\left(\sum_{k=1}^{s}X_{ik}U_{ij}\right)\right\\}^{2}}{\left\\{n^{1/2}\xi_{0}-C_{3}^{-1}\phi_{n}\right\\}^{2}}$
$\displaystyle\geq
1-\frac{n^{-1}\text{Var}\left[\bm{\varepsilon}^{\top}\mathbf{X}_{j}+\bm{\varepsilon}^{\top}\mathbf{U}_{j}+b_{1}\sum_{i=1}^{n}\left(\sum_{k=1}^{s}X_{ik}U_{ij}\right)\right]}{\left\\{n^{1/2}\xi_{0}-C_{3}^{-1}\phi_{n}\right\\}^{2}}.$

Now, consider the numerator of the last expression. Note that,

$\text{Cov}(\bm{\varepsilon}^{\top}\operatorname{\mathbf{X}}_{j},\bm{\varepsilon}^{\top}\operatorname{\mathbf{U}}_{j})=\sum_{i=1}^{n}X_{ij}\text{Cov}\left(\epsilon_{i},\epsilon_{i}U_{ji}\right)=\sum_{i=1}^{n}X_{ij}\left[E(\epsilon_{i}^{2})E(U_{ij})-E(\epsilon_{i})E(\epsilon_{i}U_{ij})\right]=0$

and also

$\text{Cov}(X_{ik}U_{ij},\bm{\varepsilon}^{\top}\operatorname{\mathbf{U}}_{j})=\text{Cov}(X_{ik}U_{ij},\varepsilon_{i}U_{ij})=X_{ik}E(\varepsilon_{i})E(U_{ij}^{2})=0.$

Hence,

$\displaystyle\mathrm{Var}\left[\bm{\varepsilon}^{\top}\mathbf{X}_{j}+\bm{\varepsilon}^{\top}\mathbf{U}_{j}+b_{1}\sum_{i=1}^{n}\left(\sum_{k=1}^{s}X_{ik}U_{ij}\right)\right]$
$\displaystyle=\textrm{Var}(\bm{\varepsilon}^{\top}\mathbf{X}_{j})+\textrm{Var}(\bm{\varepsilon}^{\top}\mathbf{U}_{j})+\textrm{Var}\left[b_{1}\sum_{i=1}^{n}\left(\sum_{k=1}^{s}X_{ik}U_{ij}\right)\right]$
$\displaystyle=\sigma^{2}\sum_{i=1}^{n}X_{ij}^{2}+\sum_{i=1}^{n}E(\epsilon_{i}^{2})E(U_{ij}^{2})+b_{1}^{2}\sum_{i=1}^{n}\left(\sum_{k=1}^{s}X_{ik}\right)^{2}E(U_{ij}^{2})$
$\displaystyle\leq
nC_{2}\sigma^{2}+n\sigma^{2}\sigma_{*}^{2}+nb_{1}^{2}C_{2}s\sigma_{*}^{2}=nO(1)+nsO(1).$

It follows that

$\mathbb{P}\left(\left|n^{-1/2}\bm{\varepsilon}^{\top}\mathbf{a}_{j}+n^{1/2}\xi_{nj}\right|<C_{3}^{-1}\phi_{n}\right)\leq\dfrac{O(1)+sO(1)}{\left\\{n^{1/2}\xi_{0}-C_{3}^{-1}\phi_{n}\right\\}^{2}}.$

Since $\lambda_{n}=o(1)$, we then have

$\frac{\phi_{n}}{n^{1/2}}=O(1)\lambda_{n}^{1/(2-\alpha)}=o(1),$

and $n^{-1}s=o(1)$. Therefore, we have

$\mathbb{P}\left(\phi_{n}>\min_{j\in\mathcal{S}}\left|n^{-1/2}\bm{\varepsilon}^{\top}\mathbf{a}_{j}+n^{1/2}\xi_{nj}\right|\right)=O(1)n^{-1}s^{2}\rightarrow
0$

and (8) follows. ∎

###  A.2 Proof of Theorem 2

Using the same notation as in the proof of Theorem 1, it suffices to show that

$\mathbb{P}\left\\{\phi_{n}>n^{-1/2}\max_{j\in\mathcal{S}^{c}}\tau_{j}\left|\bm{\varepsilon}^{\top}\mathbf{a}_{j}+n\xi_{nj}\right|\right\\}\rightarrow
1.$ (11)

To prove (11), by condition (C(C11)), for all $j\in\mathcal{S}^{c}$, we obtain

$n^{1/2}\left|\tilde{\xi}_{nj}\right|=n^{-1/2}\left|\sum_{k=1}^{s}\sum_{i=1}^{n}X_{ik}X_{ij}\beta_{0k}\right|\leq
n^{-1/2}b_{1}\left(\sum_{l=1}^{s}\left|\sum_{i=1}^{n}X_{ik}X_{ij}\right|\right)\leq
b_{1}d_{0}n^{1/2-\theta}s=c_{1}n^{1/2-\theta}s$

with $c_{1}=b_{1}d_{0}$. Next,

$\displaystyle\mathbb{P}\left(\phi_{n}>n^{-1/2}\max_{j\in\mathcal{S}^{c}}\tau_{j}\left|\bm{\varepsilon}^{\top}\mathbf{a}_{j}+n\xi_{nj}\right|\right)\geq\mathbb{P}\left(C_{4}^{-1}\phi_{n}>n^{-1/2}\max_{j\in\mathcal{S}^{c}}\left|\bm{\varepsilon}^{\top}\mathbf{a}_{j}\right|+n^{1/2}\max_{j\in\mathcal{S}^{c}}\left|\xi_{nj}\right|\right)$
$\displaystyle\geq\mathbb{P}\left(C_{4}^{-1}\phi_{n}>n^{-1/2}\max_{j\in\mathcal{S}^{c}}\left|\bm{\varepsilon}^{\top}\mathbf{X}_{j}\right|+n^{-1/2}\max_{j\in\mathcal{S}^{c}}\left|\bm{\varepsilon}^{\top}\operatorname{\mathbf{U}}_{j}\right|+n^{1/2}\max_{j\in\mathcal{S}^{c}}|\tilde{\xi}_{nj}|\right.$
$\displaystyle~{}\qquad\qquad\qquad\qquad\qquad\qquad\qquad\left.+\max_{j\in\mathcal{S}^{c}}\left|n^{-1/2}\sum_{i=1}^{n}\left(\sum_{k=1}^{s}X_{ik}\beta_{0k}\right)U_{ij}\right|\right)$
$\displaystyle\geq\mathbb{P}\left(C_{4}^{-1}\phi_{n}>n^{-1/2}\max_{j\in\mathcal{S}^{c}}\left|\bm{\varepsilon}^{\top}\mathbf{X}_{j}\right|+n^{-1/2}\max_{j\in\mathcal{S}^{c}}\left|\bm{\varepsilon}^{\top}\operatorname{\mathbf{U}}_{j}\right|+c_{1}n^{1/2-\theta}s+b_{1}\max_{j\in\mathcal{S}^{c}}\left|n^{-1/2}\sum_{i=1}^{n}\left(\sum_{k=1}^{s}X_{ik}\right)U_{ij}\right|\right)$
$\displaystyle=1-$
$\displaystyle\mathbb{P}\left(n^{-1/2}\max_{j\in\mathcal{S}^{c}}\left|\bm{\varepsilon}^{\top}\mathbf{X}_{j}\right|+n^{-1/2}\max_{j\in\mathcal{S}^{c}}\left|\bm{\varepsilon}^{\top}\operatorname{\mathbf{U}}_{j}\right|+b_{1}\max_{j\in\mathcal{S}^{c}}\left|n^{-1/2}\sum_{i=1}^{n}\left(\sum_{k=1}^{s}X_{ik}\right)U_{ij}\right|\geq
C_{4}^{-1}\phi_{n}-c_{1}n^{1/2-\theta}s\right)$ $\displaystyle\geq
1-\frac{\mathbb{E}\left(n^{-1/2}\max_{j\in\mathcal{S}^{c}}\left|\bm{\varepsilon}^{\top}\operatorname{\mathbf{X}}_{j}\right|+n^{-1/2}\max_{j\in\mathcal{S}^{c}}\left|\bm{\varepsilon}^{\top}\operatorname{\mathbf{U}}_{j}\right|+b_{1}\max_{j\in\mathcal{S}^{c}}\left|n^{-1/2}\sum_{i=1}^{n}\left(\sum_{k=1}^{s}X_{ik}\right)U_{ij}\right|\right)}{C_{4}^{-1}\phi_{n}-c_{1}n^{1/2-\theta}s}.$

Consider each term in the numerator. For the first term, by property (B1),
because $\varepsilon_{i}$ is sub-Gaussian with common variance
$\sigma^{2},~{}i=1,\ldots,n$, then
$\bm{\varepsilon}^{\top}\operatorname{\mathbf{X}}_{j}$ is sub-Gaussian with
variance $\sigma^{2}\sum_{i=1}^{n}X_{ij}^{2}$ and variance proxy
$C_{2}n\sigma^{2}$. Therefore, by property (B3), we have

$\mathbb{E}\left\\{n^{-1/2}\max_{j\in\mathcal{S}^{c}}\left|\bm{\varepsilon}^{\top}\operatorname{\mathbf{X}}_{j}\right|\right\\}\leq
C_{2}^{1/2}\sigma\sqrt{\log(2m)}=O(1)\left\\{\log(2m)\right\\}^{1/2}.$

For the second term, by property (B3), each variable $\varepsilon_{i}{U}_{ij}$
is sub-exponential with parameter $(\nu_{0},t_{0})$, so
$\bm{\varepsilon}^{\top}\operatorname{\mathbf{U}}_{j}=\sum_{i=1}^{n}\epsilon_{i}U_{ij}$
is sub-exponential with parameter $(n\nu_{0},t_{0})$. Hence, by property (B4),
we have

$\mathbb{E}\left\\{n^{-1/2}\max_{j\in\mathcal{S}^{c}}\left|\bm{\varepsilon}^{\top}\operatorname{\mathbf{U}}_{j}\right|\right\\}\leq\nu_{0}^{1/2}\sqrt{\log(2m)}+n^{-1/2}t_{0}\log(2m).$

For the third term, by a similar argument, each variable
$\sum_{i=1}^{n}\left(\sum_{k=1}^{s}X_{ik}\right)U_{ij}$ is sub-Gaussian with
variance proxy $C_{2}\sigma_{*}^{2}ns$. Hence,

$\mathbb{E}\left\\{\max_{j\in\mathcal{S}^{c}}\left|n^{-1/2}\sum_{i=1}^{n}\left(\sum_{k=1}^{s}X_{ik}\right)U_{ij}\right|\right\\}\leq
C_{2}^{1/2}\sigma_{*}\sqrt{s\log(2m)}=O(1)\sqrt{s\log(2m)}.$

Putting the results together, we have

$\mathbb{P}\left(\phi_{n}>n^{-1/2}\max_{j\in\mathcal{S}^{c}}\tau_{j}\left|\bm{\varepsilon}^{\top}\mathbf{a}_{j}+n\xi_{nj}\right|\right)\geq
1-\frac{\log(2m)^{1/2}\left\\{O(1)+\sqrt{s}+n^{-1/2}\sqrt{\log(2m)}\right\\}}{{C_{4}^{-1}\phi_{n}-c_{1}n^{1/2-\theta}s}}.$
(12)

Next, by condition (C(C11)C(C11)a),

$\frac{n^{1/2-\theta}s}{C_{4}^{-1}\phi_{n}}=O(1)\left(\frac{s^{(2-\alpha)}}{\lambda_{n}n^{\theta(2-\alpha)}}\right)^{1/(2-\alpha)}=o(1),$

and by condition (C(C11)C(C11)b),

$\frac{\log(2m)}{C_{4}^{-1}\phi_{n}}=O(1)\left(\frac{\log(2m)^{(2-\alpha)}}{\lambda_{n}n^{(2-\alpha)/2}}\right)^{1/(2-\alpha)}=o(1).$

Therefore, $\sqrt{\log(2m)}/\left(C_{4}^{-1}\phi_{n}\right)=o(1)$ and also
$\sqrt{s\log(2m)}/\left(C_{4}^{-1}\phi_{n}\right)=o(1)$. Finally, we have

$\mathbb{P}\left(\phi_{n}>n^{-1/2}\max_{j\in\mathcal{S}^{c}}\tau_{j}\left|\bm{\varepsilon}^{\top}\mathbf{a}_{j}+n\xi_{nj}\right|\right)\rightarrow
1,$

as required.

###  A.3 Proof of Theorem 3

Recall that the screening set for SISc is defined as

$\hat{\mathcal{Q}}_{\text{SISc}}=\left\\{1\leq i\leq
p:\lvert\tilde{\beta}_{i}\rvert\text{ is among the first }d=\lfloor\gamma
n\rfloor\text{ largest of all }\right\\}.$

In this proof, we will call this set $\hat{\mathcal{Q}}_{\gamma}$ to emphasize
the dependence of the screening set on $\gamma$, while the subscript SISc is
omitted to ease the notation. Following Fan & Lv (2007), the proof consists of
two main steps. In the first step, we define the following set

$\widetilde{\mathcal{M}}_{\delta}^{1}=\left\\{1\leq i\leq
p:\lvert\tilde{\beta}_{i}\rvert\text{ is among the first }[\delta p]\text{
largest of all }\right\\},$

and prove that for some positive constant $C$, we have

$\mathbb{P}\left(\widetilde{\mathcal{M}}_{\delta}^{1}\supseteq\mathcal{S}\right)=O(p\exp(-Cn)).$
(13)

In the second step, we will apply the dimensional reduction procedure above
successively until the number of covariates to be kept is smaller than $n$.

#### Step 1.

We will prove (13) by bounding $\|\tilde{\operatorname{\bm{\beta}}}\|_{2}^{2}$
from above and $|\tilde{\beta}_{j}|,~{}j\in\mathcal{S}$ from below. Note that

$\tilde{\beta}_{j}=\dfrac{n^{-1}\sum_{i=1}^{n}W_{ij}y_{i}}{n^{-1}\sum_{i=1}^{n}W_{ij}^{2}-\sigma_{j}^{2}},\quad
j=1,\ldots,p,$ (14)

so for the denominator of (14), we have

$\mathbb{E}\left(\dfrac{1}{n}\sum_{i=1}^{n}W_{ij}^{2}-\sigma_{j}^{2}\right)=\dfrac{1}{n}\sum_{i=1}^{n}X_{ij}^{2},$

and

$\dfrac{1}{n}\sum_{i=1}^{n}W_{ij}^{2}-\sigma_{j}^{2}-\dfrac{1}{n}\sum_{i=1}^{n}X_{ij}^{2}=\dfrac{2}{n}\sum_{i=1}^{n}X_{ij}U_{ij}+\dfrac{1}{n}\sum_{i=1}^{n}U_{ij}^{2}-\sigma_{j}^{2}.$

By condition (C(C11)) and properties of sub-Gaussian random variables, the
variable $\sum_{i=1}^{n}X_{ij}U_{ij}$ is zero-mean sub-Gaussian with variance
proxy $(\sum_{i=1}^{n}X_{ij}^{2})\sigma_{*}^{2}$, which is of order $O(n)$ by
condition (C(C11)). Therefore, by the Hoeffding inequality, there exists
positive constants $d_{1}$ and $K_{1}$ such that

$\mathbb{P}\left(\left|\dfrac{1}{n}\sum_{i=1}^{n}X_{ij}U_{ij}\right|\geq
d_{1}\right)=O\left(\exp\left(-K_{1}n\right)\right).$

Similarly, by condition (C(C11)), the variables $U_{ij}^{2}$’s are sub-
exponential with (finite) parameter ($\nu_{0},\alpha_{0})$, hence
$\sum_{i=1}^{n}U_{ij}^{2}$ is also sub-exponential with parameter
$(\nu_{0}\sqrt{n},\alpha_{0})$. Because
$\mathbb{E}(U_{ij}^{2})=\sigma_{j}^{2}$, there exist constants
$d_{2}\in[0,\nu_{0}^{2}/\alpha_{0}]$ and $K_{2}$ such that

$\mathbb{P}\left(\left|\dfrac{1}{n}\sum_{i=1}^{n}U_{ij}^{2}-\sigma_{j}^{2}\right|\geq
d_{2}\right)=O\left(\exp\left(-K_{2}n\right)\right).$

Combining these results, for some positive constants $d_{3}$ and $K_{3}$, we
have

$P\left(\left|\dfrac{1}{n}\sum_{i=1}^{n}W_{ij}^{2}-\sigma_{j}^{2}-\dfrac{1}{n}\sum_{i=1}^{n}X_{ij}^{2}\right|\geq
d_{3}\right)=O(\exp(-K_{3}n)).$ (15)

By condition (C(C11)), the term $n^{-1}\sum_{i=1}^{n}X_{ij}^{2}=O(1)$, so
inequality (15) implies that $n^{-1}\sum_{i=1}^{n}W_{ij}^{2}-\sigma_{j}^{2}$
is both bounded above and bounded away from zero with probability tending to
one as $n\to\infty$. Now for the numerator of $\eqref{eq:tildebetaj}$, we have

$\mathbb{E}\left(\dfrac{1}{n}\sum_{i=1}^{n}W_{ij}y_{i}\right)=\mathbb{E}\left(\dfrac{1}{n}\sum_{i=1}^{n}X_{ij}y_{i}\right)=\xi_{nj},$

and

$\dfrac{1}{n}\sum_{i=1}^{n}W_{ij}y_{i}-\xi_{nj}=\dfrac{1}{n}\sum_{i=1}^{n}X_{ij}\varepsilon_{i}+\dfrac{1}{n}\sum_{j=1}^{n}U_{ij}\varepsilon_{i}.$

By conditions (C(C11)) and (C(C11)), using similar arguments, we obtain

$\mathbb{P}\left(\left|\dfrac{1}{n}\sum_{i=1}^{n}W_{ij}y_{i}-\xi_{nj}\right|\geq
d_{4}\right)=O(\exp(-K_{4}n)),$

for some positive constant $d_{4}$ and $K_{4}$. Hence if $j\in\mathcal{S}$,
along with the condition (C(C11)), the above inequality implies that
$|n^{-1}\sum_{i=1}^{n}W_{ij}y_{i}|$ is bounded away from zero with probability
at least $1-O(\exp(-Cn))$. Using the union bound, there exist positive
constants $d_{7}=d_{6}/d_{5}$ with
$d_{5}=n^{-1}\sum_{i=1}^{n}X_{ij}^{2}+d_{3}$ and $d_{6}=\xi_{nj}-d_{4}$ such
that for $j\in\mathcal{S}$,

$\displaystyle\mathbb{P}(|\tilde{\beta_{j}}|\geq d_{7})$
$\displaystyle\geq\mathbb{P}\left(\left|\dfrac{1}{n}\sum_{i=1}^{n}W_{ij}y_{i}\right|\geq
d_{6}\text{ and
}\left|\dfrac{1}{n}\sum_{i=1}^{n}W_{ij}^{2}-\sigma_{j}^{2}\right|\leq
d_{5}\right)$
$\displaystyle=1-\mathbb{P}\left(\left|\dfrac{1}{n}\sum_{i=1}^{n}W_{ij}y_{i}\right|\leq
d_{6}\text{ or
}\left|\dfrac{1}{n}\sum_{i=1}^{n}W_{ij}^{2}-\sigma_{j}^{2}\right|\geq
d_{5}\right)$ $\displaystyle\geq
1-\left\\{\mathbb{P}\left(\left|\dfrac{1}{n}\sum_{i=1}^{n}W_{ij}y_{i}\right|\leq
d_{6}\right)+\mathbb{P}\left(\left|\dfrac{1}{n}\sum_{i=1}^{n}W_{ij}^{2}-\sigma_{j}^{2}\right|\geq
d_{5}\right)\right\\}$
$\displaystyle=1-(O\exp(-K_{4}n)+O(\exp(-K_{3}n)))=1-O(\exp(-K_{5}n)),$ (16)

for some constant $C$. The final inequality states that for all
$j\in\mathcal{S}$ we have $|\hat{\beta}_{j}|$ is bounded away from 0 with
probability tending to one as $n\to\infty$.

Next, we will bound the norm $\lVert\tilde{\bm{\beta}}\rVert_{2}^{2}$ from
above. Let $\bm{\Lambda}$ denote the $p\times p$ diagonal matrix with elements
$n^{-1}\sum_{i=1}^{n}W_{ij}^{2}-\sigma_{j}^{2}~{}j=1,\ldots,p$, then we can
write $\tilde{\bm{\beta}}$ as

$\tilde{\bm{\beta}}=\bm{\Lambda}^{-1}\left(\dfrac{1}{n}\operatorname{\mathbf{W}}^{\top}\mathbf{y}\right)=\bm{\Lambda}^{-1}\left(\dfrac{1}{n}\operatorname{\mathbf{X}}^{\top}\operatorname{\mathbf{X}}\bm{\beta}_{0}+\dfrac{1}{n}\operatorname{\mathbf{U}}^{\top}\operatorname{\mathbf{X}}\bm{\beta}_{0}+\dfrac{1}{n}\operatorname{\mathbf{X}}^{\top}\bm{\varepsilon}+\dfrac{1}{n}\operatorname{\mathbf{U}}^{\top}\bm{\varepsilon}\right).$

Hence, by the submultiplicity of matrix norm and the triangle inequality, we
have

$\lVert\tilde{\bm{\beta}}\rVert_{2}\leq\|\bm{\Lambda}^{-1}\|_{2}\left(\left\|\dfrac{1}{n}\operatorname{\mathbf{X}}^{\top}\operatorname{\mathbf{X}}\bm{\beta}_{0}\right\|_{2}+\left\|\dfrac{1}{n}\operatorname{\mathbf{U}}^{\top}\operatorname{\mathbf{X}}\bm{\beta}_{0}\right\|_{2}+\left\|\dfrac{1}{n}\operatorname{\mathbf{X}}^{\top}\bm{\varepsilon}\right\|_{2}+\left\|\dfrac{1}{n}\operatorname{\mathbf{U}}^{\top}\bm{\varepsilon}\right\|_{2}\right)$
(17)

where
$\left\|\bm{\Lambda}^{-1}\right\|_{2}=\max_{j=1,\ldots,p}\left\\{1/(|n^{-1}\sum_{i=1}^{n}W_{ij}^{2}-\sigma_{j}^{2}|)\right\\}=\left\\{\min_{j=1,\ldots,p}|n^{-1}\sum_{i=1}^{n}W_{ij}^{2}-\sigma_{j}^{2}|\right\\}^{-1}$
is the magnitude of the largest eigenvalue of $\bm{\Lambda}^{-1}$ . From
equation (15) and condition (C(C11)), there exist sufficiently large constants
$d_{8}$ and $K_{8}$ such that

$\displaystyle\mathbb{P}\left(\left\|\bm{\Lambda}^{-1}\right\|_{2}\leq
d_{8}\right)$
$\displaystyle=\mathbb{P}\left(\min_{j=1,\ldots,p}\left\\{\left|\dfrac{1}{n}\sum_{i=1}^{n}W_{ij}^{2}-\sigma_{j}^{2}\right|\right\\}\geq\dfrac{1}{d_{8}}\right)$
$\displaystyle=\mathbb{P}\left(\bigcap\limits_{j=1,\ldots,p}\left|\dfrac{1}{n}\sum_{i=1}^{n}W_{ij}^{2}-\sigma_{j}^{2}\right|\geq\dfrac{1}{d_{8}}\right)$
$\displaystyle=1-\mathbb{P}\left(\bigcup\limits_{j=1,\ldots,p}\left|\dfrac{1}{n}\sum_{i=1}^{n}W_{ij}^{2}-\sigma_{j}^{2}\right|\leq\dfrac{1}{d_{8}}\right)$
$\displaystyle\geq
1-\sum_{j=1}^{p}\mathbb{P}\left(\left|\dfrac{1}{n}\sum_{i=1}^{n}W_{ij}^{2}-\sigma_{j}^{2}\right|\leq\dfrac{1}{d_{8}}\right)$
$\displaystyle=1-O(p\exp(-K_{8}n)).$ (18)

Since $\log(p)/n\to C_{1}$ by condition (C(C11)), as long as $K_{8}>C_{1}$,
the operator norm $\left\|\bm{\Lambda}^{-1}\right\|_{2}$ is bounded above by a
constant with probability tending to one $n\to\infty$.

Next, we will bound each term in the bracket of (17). First,

$\left\|\dfrac{1}{n}\operatorname{\mathbf{X}}^{\top}\operatorname{\mathbf{X}}\bm{\beta}_{0}\right\|_{2}^{2}=\operatorname{\bm{\beta}}_{0}^{\top}\left(\dfrac{1}{n}\operatorname{\mathbf{X}}^{\top}\operatorname{\mathbf{X}}\right)^{2}\operatorname{\bm{\beta}}_{0}\leq\lambda_{\max}^{2}\left(\dfrac{1}{n}\operatorname{\mathbf{X}}^{\top}\operatorname{\mathbf{X}}\right)\|\operatorname{\bm{\beta}}_{0}\|_{2}^{2}=O(sn^{2\tau_{2}})$
(19)

by conditions (C(C11)) and (C(C11)). For the second term,

$\displaystyle\left\|\dfrac{1}{n}\operatorname{\mathbf{U}}^{\top}\operatorname{\mathbf{X}}\bm{\beta}_{0}\right\|_{2}^{2}$
$\displaystyle=\dfrac{1}{n^{2}}\operatorname{\bm{\beta}}_{0}^{\top}\operatorname{\mathbf{X}}^{\top}\operatorname{\mathbf{U}}\operatorname{\mathbf{U}}^{\top}\operatorname{\mathbf{X}}\operatorname{\bm{\beta}}_{0}=\dfrac{1}{n^{2}}\operatorname{\bm{\beta}}_{0}^{\top}\operatorname{\mathbf{X}}^{\top}\operatorname{\mathbf{Z}}\bm{\Sigma}_{u}\operatorname{\mathbf{Z}}^{\top}\operatorname{\mathbf{X}}\operatorname{\bm{\beta}}_{0}$
$\displaystyle\leq\dfrac{1}{n^{2}}\lambda_{\max}(\bm{\Sigma}_{u}){\operatorname{\bm{\beta}}}_{0}^{\top}\operatorname{\mathbf{X}}^{\top}\operatorname{\mathbf{Z}}\operatorname{\mathbf{Z}}^{\top}\operatorname{\mathbf{X}}\operatorname{\bm{\beta}}_{0}$
$\displaystyle\leq\dfrac{1}{n^{2}}p\lambda_{\max}(\bm{\Sigma}_{u})\lambda_{\max}(p^{-1}\operatorname{\mathbf{Z}}\operatorname{\mathbf{Z}}^{\top})\operatorname{\bm{\beta}}_{0}^{\top}\operatorname{\mathbf{X}}^{\top}\operatorname{\mathbf{X}}\operatorname{\bm{\beta}}_{0}$
$\displaystyle\leq
p\lambda_{\max}(\bm{\Sigma}_{u})\lambda_{\max}(p^{-1}\operatorname{\mathbf{Z}}\operatorname{\mathbf{Z}}^{\top})\lambda_{\max}\left(\dfrac{1}{n}\operatorname{\mathbf{X}}^{\top}\operatorname{\mathbf{X}}\right)\|\operatorname{\bm{\beta}}_{0}\|_{2}^{2}$
$\displaystyle\stackrel{{\scriptstyle(i)}}{{\leq}}pn^{\tau_{1}+\tau_{2}-1}\lambda_{\max}(p^{-1}\operatorname{\mathbf{Z}}\operatorname{\mathbf{Z}}^{\top})\|\operatorname{\bm{\beta}}_{0}\|_{2}^{2},$

where step $(i)$ follows conditions (C(C11)C(C11)b) and (C(C11)). Therefore,
by conditions (C(C11)C(C11)a) and (C(C11)) there exist positive constants
$d_{9}$ and $K_{9}$ such that

$\mathbb{P}\left(\left\|\dfrac{1}{n}\operatorname{\mathbf{U}}^{\top}\operatorname{\mathbf{X}}\bm{\beta}_{0}\right\|_{2}^{2}\geq
d_{9}psn^{\tau_{1}+\tau_{2}-1})\right)\leq O(\exp(-K_{9}n)).$

For the third term, we have

$\left\|\dfrac{1}{n}\operatorname{\mathbf{X}}^{\top}\bm{\varepsilon}\right\|_{2}^{2}=\dfrac{1}{n}\bm{\varepsilon}^{\top}\left(\dfrac{1}{n}\operatorname{\mathbf{X}}^{\top}\operatorname{\mathbf{X}}\right)\bm{\varepsilon}\leq\lambda_{\max}\left(\dfrac{1}{n}\operatorname{\mathbf{X}}^{\top}\operatorname{\mathbf{X}}\right)\dfrac{1}{n}\|\bm{\varepsilon}\|_{2}^{2}.$

By condition (C(C11)), each term $\varepsilon_{i}^{2}$ is sub-exponential with
finite parameters; furthermore, since $\varepsilon_{1},\ldots,\varepsilon_{n}$
are mutually independent, there exist positive constants $d_{10}$ and $K_{10}$
such that

$\mathbb{P}\left(\dfrac{1}{n}\|\bm{\varepsilon}\|_{2}^{2}\geq
d_{10}\right)=\mathbb{P}\left(\dfrac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}^{2}\geq
d_{10}\right)=O(\exp(-K_{10}n)).$ (20)

As a result, we obtain

$\mathbb{P}\left(\left\|\dfrac{1}{n}\operatorname{\mathbf{X}}^{\top}\bm{\varepsilon}\right\|_{2}^{2}\geq
d_{11}n^{\tau_{2}}\right)=O(\exp(-K_{11}n))$

for some positive constants $d_{11}$ and $K_{11}$. Finally, for the last term
in the bracket of (17), we have

$\displaystyle\left\|\dfrac{1}{n}\operatorname{\mathbf{U}}^{\top}\bm{\varepsilon}\right\|_{2}^{2}$
$\displaystyle=\dfrac{1}{n^{2}}\bm{\varepsilon}^{\top}\operatorname{\mathbf{U}}\operatorname{\mathbf{U}}^{\top}\bm{\varepsilon}=\dfrac{1}{n^{2}}\bm{\varepsilon}^{\top}\operatorname{\mathbf{Z}}\bm{\Sigma}_{u}\operatorname{\mathbf{Z}}^{\top}\bm{\varepsilon}$
$\displaystyle\leq\dfrac{1}{n^{2}}\lambda_{\max}(\bm{\Sigma}_{u})\bm{\varepsilon}^{\top}\operatorname{\mathbf{Z}}\operatorname{\mathbf{Z}}^{\top}\bm{\varepsilon}$
$\displaystyle\leq\dfrac{1}{n}p\lambda_{\max}(\bm{\Sigma}_{u})\lambda_{\max}(p^{-1}\operatorname{\mathbf{Z}}\operatorname{\mathbf{Z}}^{\top})\left(\dfrac{1}{n}\|\bm{\varepsilon}\|_{2}^{2}\right)$
$\displaystyle\leq
c_{3}pn^{\tau_{1}-1}\lambda_{\max}(p^{-1}\operatorname{\mathbf{Z}}\operatorname{\mathbf{Z}}^{\top})\left(\dfrac{1}{n}\|\bm{\varepsilon}\|_{2}^{2}\right).$
(21)

Let $d_{12}=c_{2}c_{3}d_{10}$, then by condition (C(C11)C(C11)a) and (20), we
obtain

$\displaystyle\mathbb{P}\left(\left\|\dfrac{1}{n}\operatorname{\mathbf{U}}^{\top}\bm{\varepsilon}\right\|_{2}^{2}\geq
d_{12}pn^{\tau_{1}-1}\right)$
$\displaystyle\leq\mathbb{P}\left\\{\lambda_{\max}(p^{-1}\operatorname{\mathbf{Z}}\operatorname{\mathbf{Z}}^{\top})\left(\dfrac{1}{n}\|\bm{\varepsilon}\|_{2}^{2}\right)\geq
c_{2}d_{10}\right\\}$
$\displaystyle=1-\mathbb{P}\left\\{\lambda_{\max}(p^{-1}\operatorname{\mathbf{Z}}\operatorname{\mathbf{Z}}^{\top})\left(\dfrac{1}{n}\|\bm{\varepsilon}\|_{2}^{2}\right)\leq
c_{2}d_{10}\right\\}$ $\displaystyle\leq
1-\mathbb{P}\left\\{\lambda_{\max}(p^{-1}\operatorname{\mathbf{Z}}\operatorname{\mathbf{Z}}^{\top})\leq
c_{2}\text{ and }\dfrac{1}{n}\|\bm{\varepsilon}\|_{2}^{2}\leq d_{10}\right\\}$
$\displaystyle=\mathbb{P}\left\\{\lambda_{\max}(p^{-1}\operatorname{\mathbf{Z}}\operatorname{\mathbf{Z}}^{\top})\geq
c_{2}\text{ or }\dfrac{1}{n}\|\bm{\varepsilon}\|_{2}^{2}\geq d_{10}\right\\}$
$\displaystyle\leq\mathbb{P}\left(\lambda_{\max}(p^{-1}\operatorname{\mathbf{Z}}\operatorname{\mathbf{Z}}^{\top})\geq
c_{2}\right)+\mathbb{P}\left(\dfrac{1}{n}\|\bm{\varepsilon}\|_{2}^{2}\geq
d_{10}\right)$ $\displaystyle=O(\exp(-K_{12}n))$ (22)

for some positive constant $K_{12}$. Substituting (18)-(22) to (17), for some
sufficiently large constants $d_{13}$ and $K_{13}>C_{1}$ and under the
condition that $p\geq n^{\tau_{2}-\tau_{1}+1}$, we have

$\mathbb{P}\left(\|\tilde{\operatorname{\bm{\beta}}}\|_{2}^{2}\geq
d_{13}psn^{\tau_{1}+\tau_{2}+-1})\right)\leq O(p\exp(-K_{13}n))$ (23)

Finally, by Bonferonni’s inequality, it follows from (16) and (23) that

$\mathbb{P}\left(\min_{i\in\mathcal{S}}|\tilde{\beta}_{i}|<d_{7}\text{ or
}\lVert\tilde{\bm{\beta}}\rVert_{2}^{2}>d_{13}psn^{\tau_{1}+\tau_{2}-1}\right)\leq
O\left(p\exp\left(-K_{14}n\right)\right),$ (24)

for some constants $K_{14}>C_{1}$. Hence, for sufficiently large constants
$C>C_{1}$, with probability $1-O\left(p\exp\left(-Cn\right)\right)$, the
magnitudes of $\tilde{\beta}_{i},i\in\mathcal{S}$ are bounded away from zero
and for some $d>0$,

$\text{card}\left\\{1\leq k\leq
p:|\tilde{\beta}_{k}|\geq\min_{i\in\mathcal{S}}|\tilde{\beta}_{i}|\right\\}\leq
dpn^{\tau_{1}+\tau_{2}+\log_{n}(s)-1}.$

Therefore, if $\delta\to 0$ and satisfies $\delta
n^{1-\tau_{1}-\tau_{2}-\log_{n}(s)}\to\infty$ as $n\to\infty$, then (13) holds
with a constant $C>0$ larger than $K_{14}$ in (24).

#### Step 2.

This step follows Step 2 in the proof of Theorem 1 of Fan & Lv (2007). In this
step, we will use $C$ to denote generic constants that are larger than $C_{1}$
in condition (C(C11)). Fix an arbitrary $r\in(0,1)$ and choose a shrinking
factor $\delta$ of the form $\delta=(\frac{n}{p})^{1/(k-r)}$ for some integer
$k\geq 1$. We successively perform dimension reduction until the number of
remaining variables drops to below the sample size $n$. In other words we
obtain a sequence of nested sets

$\widetilde{\mathcal{M}}_{\delta}^{k}\subset\widetilde{\mathcal{M}}_{\delta}^{k-1}\subset\ldots\subset\widetilde{\mathcal{M}}_{\delta}^{1},$

where each set
$\widetilde{\mathcal{M}}_{\delta}=\widetilde{\mathcal{M}}_{\delta}^{j}$ has
cardinality $\lceil\delta^{j}p\rceil$ and
$d=\lceil\delta^{k}p\rceil=\lceil\delta^{r}n\rceil<n$ but
$\lceil\delta^{k-1}p\rceil=\lceil\delta^{r-1}n\rceil>n$. Hence we see that
$\widetilde{\mathcal{M}}_{\delta}=\hat{\mathcal{Q}}_{\gamma}$, with
$\gamma=\delta^{r}<1$.

Next, fix an arbitrary $\theta_{1}\in(0,1-\tau_{1}-\tau_{2}-\log_{n}(s))$ and
pick some $r<1$ very close to 1 such that
$\theta_{0}=\theta_{1}/r<1-\tau_{1}-\tau_{2}-\log_{n}(s)$, and choose a
sequence of integers $k\geq 1$ in a way such that when $n\to\infty$,

$\delta n^{1-\tau_{1}-\tau_{2}-\log_{n}(s)}\to\infty\text{ and }\delta
n^{\theta_{0}}\to 0,$ (25)

where $\delta=(\frac{n}{p})^{1/(k-r)}$. Therefore, with the above dimension
reduction process, we can raise both sides of (25) to the $r$th power, and
hence the set $\widetilde{\mathcal{M}}_{\delta}=\hat{\mathcal{Q}}_{\gamma}$
with $\gamma=\delta^{r}$ satisfies

$\gamma n^{r(1-\tau_{1}-\tau_{2}-\log_{n}(s))}\to\infty\text{ and }\gamma
n^{\theta_{1}}\to 0.$

Since for any principal submatrix $\bm{\Sigma}_{u}^{0}$ of $\bm{\Sigma}_{u}$,
we have
$\lambda_{\max}(\bm{\Sigma}_{u}^{0})<\lambda_{\max}(\bm{\Sigma}_{u})\leq
c_{2}n^{\tau_{1}}$, and that property (6) in condition (C(C11)C(C11)a) holds
for any $n\times\tilde{p}$ submatrix $\widetilde{\mathbf{Z}}$ of $\mathbf{Z}$
with $c_{1}n<\tilde{p}\leq p$, for some constant $C>0$, in each step $1\leq
i\leq k$ in the above dimension reduction framework, we have

$\mathbb{P}\left(\widetilde{\mathcal{M}}_{\delta}^{i}\supseteq\mathcal{S}\mid\widetilde{\mathcal{M}}_{\delta}^{i-1}\supseteq\mathcal{S}\right)=1-O\left(p\exp\left(-Cn\right)\right),$

and hence by Bonferroni’s inequality we have

$\mathbb{P}\left(\hat{\mathcal{Q}}_{\gamma}\supseteq\mathcal{S}\right)=1-O\left(kp\exp\left(-Cn\right)\right).$
(26)

It follows from (25) that we require $\delta\to 0$ for all $r$ when both
$n\to\infty$ and $p\ \to\infty$ (i.e $\log(n)/k-\log(p)/k\to-\infty)$. Since
$p>n$ implies $k=O(\log(p)/\log(n))$, a suitable increase of the constant
$C>0$ in (26) gives

$\mathbb{P}\left(\hat{\mathcal{Q}}_{\gamma}\supseteq\mathcal{S}\right)=1-O\left(p\exp\left(-Cn\right)\right).$

This probability bound holds for any $\gamma\sim cn^{-\tilde{\theta}}$, with
$\tilde{\theta}<1-\tau_{1}-\tau_{2}-\log_{n}(s)$ and $c>0$, completing the
proof.

##  B Preprocessing for the microarray data

In this subsection, we provide more detail on the steps that we used to
preprocess the microarray data in Section 4 of the main paper.

As noted in Hein et al., (2005), on Affymetrix “GeneChips” oligonucleotide
arrays, each gene is represented by a probe set, consisting of a number of
probe pairs. A probe pair further contains a perfect match (PM) probe and a
mismatch probe (MM). On the one hand, the intensity observed for the PM
measurement for probe $r$ at gene $j$ was assumed due partly to the binding of
the cRNA that perfectly matches the sequence on the array, denoted as $S_{jr}$
(true signal), and partly to hybridization of the cRNA that did not perfectly
match the sequence, denoted as $H_{jr}$ (nonspecific hybridization). On the
other hand, the intensity observed for the MM measurement for probe $r$ at
gene $j$ was due partly to binding of a fraction $\phi\in(0,1)$ of the true
signal and partly due to nonspecific hybridization. Furthermore, it was
assumed that

$\displaystyle\text{PM}_{jr}$ $\displaystyle\sim N(S_{jr}+H_{jr},\tau^{2}),$
$\displaystyle\text{MM}_{jr}$ $\displaystyle\sim N(\phi
S_{jr}+H_{jr},\tau^{2}),$

where $N(\mu,\sigma^{2})$ denotes the normal distribution with mean $\mu$ and
variance $\sigma^{2}$. Next, Hein et al., (2005) modeled the true signal
$S_{gj}$ and nonspecific hybridization $H_{gj}$ on the log scale. Since they
allowed them to be zero, they assumed that

$\displaystyle\log(S_{jr}+1)\sim\text{TN}(X_{j},\sigma_{j}^{2}),$
$\displaystyle\log(H_{jr}+1)\sim\text{TN}(\lambda,\eta^{2}),$

where $\text{TN}(\mu,\sigma^{2})$ denotes the normal distribution with mean
$\mu$ and variance $\sigma^{2}$ left truncated at 0. Furthermore, to account
for outlying probes and stabilize the gene-specific variance parameters, all
the variances $\sigma_{g}^{2}$ were assumed to be exchangeable,

$\log(\sigma_{j}^{2})\stackrel{{\scriptstyle iid}}{{\sim}}N(a,b^{2}).$

Finally, the Bayesian model was fully specified by assuming non-informative
priors on $X_{j}$, $\phi$, $\lambda$, $\tau^{-2}$ and $\eta^{-2}$, while $a$
and $b$ were fixed at values obtained by an empirical procedure. The primary
parameter of interest was $X_{j}$, the mean gene expression level for gene $j$
on the log scale.

We fitted the above Bayesian model for each patient $i=1,\ldots,n$ separately
by the BGX package of Hein et al., (2005) in R and obtained the posterior
distribution for $X_{ij}$. Let $\widehat{W}_{ij}$ and
$\text{var}(\widehat{X}_{ij})$ denote the posterior mean and variance of
$X_{ij}$ respectively; so we considered $\widehat{W}_{ij}$ as a surrogate for
the true $X_{ij}$ contaminated by some measurement error $U_{ij}$ with
variance $\text{Var}(\widehat{X}_{ij})$. Next, for each gene $j$, we
standardized the posterior mean $\tilde{W}_{ij}$ to obtain
$W_{ij}={(\widehat{W}_{ij}-\overline{W}_{j})}/{s_{j}},\ i=1,\ldots,n$, where
$\overline{W}_{j}=n^{-1}\sum_{i=1}^{n}\widehat{W}_{ij}$ and
$s_{j}^{2}=n^{-1}\sum_{i=1}^{n}(\widehat{W}_{ij}-\overline{W}_{j})^{2}$. To
simplify the analysis, we assumed that measurement error was independent of
the patient’s true gene expression levels and that the associated variance was
constant across patients for a given gene; hence, the matrix $\bm{\Sigma}_{u}$
was set to be diagonal. We averaged the posterior variance,
$\widehat{\sigma}^{2}_{u,j}=n^{-1}\sum_{i=1}^{n}\text{Var}(\widehat{X}_{ij})$,
and the diagonal elements of the measurement error covariance matrix
$\bm{\Sigma}_{u}$ were estimated as
$(\widehat{\bm{\Sigma}}_{u})_{j,j}={\widehat{\sigma}^{2}_{uj}}/s^{2}_{j},j=1,\ldots,p$.
As noted by Sørensen et al., (2015), if the measurement error variance is too
large compared to the between-sample variance of the between-patient
variability, then little can be done. Hence, only the $p=993$ genes with
$\widehat{\sigma}^{2}_{u,j}<0.5s_{j}^{2}$, i.e. estimated noise-to-signal
ratio less than $1$, were retained for further analysis.

##  C Efficiency of the screening procedures

In this section, we present additional simulation results to investigate the
efficiency of the proposed screening procedures. We generated data in the same
way as in the subsection 3.1 of the main paper, with the exception that we
considered one more scenario for $\bm{\Sigma}_{u}$, the covariance matrix of
the measurement error and set $\rho_{x}=0.5$. In this new scenario,
$\bm{\Sigma}_{u}$ was set to have a homogeneous structure, with all off-
diagonal elements equal to $0.2$ and diagonal elements equal to $0.4$; hence
the measurement errors on all covariates were correlated with one another.

We report the mean false positive rate (FPR) and mean false negative rate
(FNR) across $500$ samples for PMSc${}_{\text{CV}}$ (Table C.1), and the
empirical cumulative distribution of the minimum number of variables that had
to be included for PMSc forward stepwise and SISc in order to retain all the
important variables (Figure 1). Note that the latter is always between $s$ and
$p$. On the one hand, a low minimum number of variables to be included
indicates a high screening efficiency; on the other hand, if the minimum
number of variables to be included is close to $p$, then dimension reduction
is not achievable via screening.

Table C.1: Average false positive rates (FPR, in percentage) and false
negative rates (FNR, in percentage) of PMSc${}_{\text{CV}}$ with tuning
parameters selected via five fold cross-validation. Standard errors are
included in parentheses.

$\bm{\Sigma}_{x}$ | $p$ |  $\bm{\Sigma}_{u}$ diagonal |  $\bm{\Sigma}_{u}$ block diagonal |  $\bm{\Sigma}_{u}$ homogeneous
---|---|---|---|---
|  | FPR | FNR | FPR | FNR | FPR | FNR
AR(1) | 1000 | 0.1 (0.1) | 0.0 (0.0) | 0.0 (0.1) | 0.0 (0.0) | 0.1 (0.1) | 0.0 (0.0)
| 2000 | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0)
Homogen | 1000 | 100.0 (0.0) | 0.0 (0.0) | 100.0 (0.0) | 0.0 (0.0) | 100.0 (0.0) | 0.0 (0.0)
2000 | 100.0 (0.0) | 0.0 (0.0) | 100.0 (0.0) | 0.0 (0.0) | 100.0 (0.0) | 0.0 (0.0)

(a) $\bm{\Sigma}_{x}$ has an homogeneous structure.

(b) $\bm{\Sigma}_{x}$ has an AR(1) structure.

Figure 1: Minimum number of variables that has to be kept for the SISc and
PMSc${}_{\text{FS}}$ to retain all the important variables in the simulation
study.

Table C.1 and Figure 1 indicate that when $\bm{\Sigma}_{x}$ had an AR(1)
structure, all the screening methods had high efficiency in variable
selection. Specifically, regardless of the structure of measurement error
variance $\bm{\Sigma}_{u}$, the PMSc${}_{\text{CV}}$ method had both false
positive and false negative rates very close to zero. For the SISc and
PMSc${}_{\text{FS}}$, the number of variables each had to keep was also close
to $5$, the true number of non-zero coefficients in
$\operatorname{\bm{\beta}}_{0}$. These results confirm both the variable
selection consistency of PMSc and the ability of SISc to reduce the number of
dimensions dramatically when all the true covariates are not highly collinear.

By contrast, when $\bm{\Sigma}_{x}$ had a homogeneous structure, the
efficiency of screening methods decreased considerably. The
PMSc${}_{\text{CV}}$ method was essentially unable to reduce the number of
variables in that case, as evidenced by the false positive rates close to
100%. Similarly, to ensure all the important variables were retained, the SISc
method had to keep 75-80% of the variables, which was usually more than the
sample size $n$. Among the three screening methods, the PMSc${}_{\text{FS}}$
was the most effective; in all the considered settings, it could effectively
reduce the dimension to below the sample size $n=500$. Looking at the
distribution of the minimum number of variables that must be kept for
PMSc${}_{\text{FS}}$ more closely, it is interesting that the efficiency of
PMSc${}_{\text{FS}}$ seemed to increase when measurement errors were more
correlated (i.e $\bm{\Sigma}_{u}$ is homogeneous).

##  D Detailed simulation results for one-stage and two-stage estimators

Table D.1: Performance of the one-stage and two-stage estimators in the
simulation study based on false positive rate (FPR, in percentage), false
negative rate (FNR, in percentage), $\ell_{2}$ estimation error, and
computation time (in seconds) when $\bm{\Sigma}_{x}$ has an AR(1) structure
with autocorrelation $\rho_{x}=0.5$. Standard error are included in
parentheses. The CoCo estimator was only computed when the number of variables
was no more than $1000$ (either in the 1st or 2nd step).

$\bm{\Sigma}_{u}$ | $p$ | Estimator | 1st step | 2nd steps | Time
---|---|---|---|---|---
|  |  | FPR | FNR | FPR | FNR | $\ell_{2}$ | 1st step | 2nd step
Diagonal | 1000 | One-stage Corrected | - | - | 2.2 (1.1) | 0.0 (0.0) | 0.35 (0.10) | 0.0 (0.0) | 49.7 (19.5)
|  | PMSc${}_{\text{CV}}$-Corrected | 0.1 (0.1) | 0.0 (0.0) | 0.1 (0.1) | 0.0 (0.0) | 0.26 (0.11) | 23.4 (0.3) | 0.4 (0.1)
|  | PMSc${}_{\text{FS}}$-Corrected | 7.5 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.28 (0.10) | 5.4 (1.1) | 1.6 (0.5)
|  | SISc-Corrected | 7.5 (0.0) | 0.0 (0.0) | 1.2 (0.6) | 0.0 (0.0) | 0.33 (0.10) | 1.5 (0.0) | 2.0 (0.6)
|  | One-stage CoCo | - | - | 7.6 (2.0) | 0.0 (0.0) | 0.47 (0.11) | 0.0 (0.0) | 516.8 (32.1)
|  | PMSc${}_{\text{CV}}$-CoCo | 0.1 (0.1) | 0.0 (0.0) | 0.1 (0.0) | 0.0 (0.0) | 0.26 (0.11) | 23.4 (0.3) | 0.0 (0.0)
|  | PMSc${}_{\text{FS}}$-CoCo | 7.5 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.27 (0.09) | 5.4 (1.1) | 0.7 (0.2)
|  | SISc-CoCo | 7.5 (0.0) | 0.0 (0.0) | 0.1 (0.1) | 0.0 (0.0) | 0.31 (0.09) | 1.5 (0.0) | 0.8 (0.2)
| 2000 | One-stage Corrected | - | - | 1.3 (0.6) | 0.0 (0.0) | 0.36 (0.10) | 0.0 (0.0) | 136.5 (55.7)
|  | PMSc${}_{\text{CV}}$-Corrected | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.26 (0.10) | 38.0 (2.0) | 0.3 (0.1)
|  | PMSc${}_{\text{FS}}$-Corrected | 3.8 (0.0) | 0.0 (0.0) | 0.6 (0.3) | 0.0 (0.0) | 0.29 (0.10) | 17.2 ( 3.1) | 2.4 ( 0.7)
|  | SISc-Corrected | 3.8 (0.0) | 0.0 (0.0) | 0.6 (0.3) | 0.0 (0.0) | 0.33 (0.10) | 2.8 (0.1) | 1.8 (0.5)
|  | PMSc${}_{\text{CV}}$-CoCo | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.26 (0.11) | 38.0 (2.0) | 0.0 (0.0)
|  | PMSc${}_{\text{FS}}$-CoCo | 3.8 (0.0) | 0.0 (0.0) | 0.1 (0.0) | 0.0 (0.0) | 0.27 (0.09) | 17.2 ( 3.1) | 0.8 ( 0.2)
|  | SISc-CoCo | 3.8 (0.0) | 0.0 (0.0) | 0.1 (0.0) | 0.0 (0.0) | 0.31 (0.09) | 2.8 (0.1) | 0.7 (0.2)
Block diagonal | 1000 | One-stage Corrected | - | - | 1.8 (0.9) | 0.0 (0.0) | 0.62 (0.17) | 0.0 (0.0) | 27.8 (14.9)
| PMSc${}_{\text{CV}}$-Corrected | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.43 (0.18) | 23.4 (0.4) | 0.4 (0.1)
|  | PMSc${}_{\text{FS}}$-Corrected | 7.5 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.4) | 0.48 (0.16) | 5.3 (1.0) | 1.9 (0.5)
|  | SISc-Corrected | 7.5 (0.0) | 0.0 (0.0) | 1.0 (0.5) | 0.0 (0.0) | 0.57 (0.18) | 1.5 (0.0) | 2.6 (0.8)
|  | One-stage CoCo | - | - | 7.4 (1.9) | 0.0 (0.0) | 0.74 (0.13) | 0.0 (0.0) | 417.4 (22.2)
|  | PMSc${}_{\text{CV}}$-CoCo | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.43 (0.18) | 23.4 (0.4) | 0.0 (0.0)
|  | PMSc${}_{\text{FS}}$-CoCo | 7.5 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.4) | 0.46 (0.13) | 4.9 (1.6) | 1.3 (0.2)
|  | SISc-CoCo | 7.5 (0.0) | 0.0 (0.0) | 0.1 (0.0) | 0.0 (0.0) | 0.54 (0.15) | 1.5 (0.0) | 1.7 (0.2)
| 2000 | One-stage Corrected | - | - | 1.0 (0.4) | 0.0 (0.9) | 0.64 (0.17) | 0.0 (0.0) | 68.0 (36.1)
|  | PMSc${}_{\text{CV}}$-Corrected | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.44 (0.17) | 37.9 (1.6) | 0.3 (0.1)
|  | PMSc${}_{\text{FS}}$-Corrected | 3.8 (0.0) | 0.0 (0.0) | 0.5 (0.2) | 0.1 (1.3) | 0.46 (0.14) | 16.7 (3.6) | 2.8 (0.8)
|  | SISc-Corrected | 3.8 (0.0) | 0.0 (0.0) | 0.5 (0.2) | 0.0 (0.0) | 0.57 (0.17) | 2.8 (0.1) | 2.3 (0.6)
|  | PMSc${}_{\text{CV}}$-CoCo | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.45 (0.18) | 37.9 (1.6) | 0.0 (0.0)
|  | PMSc${}_{\text{FS}}$-CoCo | 3.8 (0.0) | 0.0 (0.0) | 0.1 (0.0) | 0.0 (0.9) | 0.44 (0.13) | 16.7 (3.6) | 1.7 (0.2)
|  | SISc-CoCo | 3.8 (0.0) | 0.0 (0.0) | 0.1 (0.0) | 0.0 (0.0) | 0.54 (0.15) | 2.8 (0.1) | 1.4 (0.2)

Table D.2: Performance of the one-stage and two-stage estimators in the
simulation study based on false positive rate (FPR, in percentage), false
negative rate (FNR, in percentage), $\ell_{2}$ estimation error, and
computation time (in seconds) when $\bm{\Sigma}_{x}$ has an homogeneous
structure with $\rho_{x}=0.5$. Standard error are included in parentheses. The
CoCo estimator was only computed when the number of variables was no more than
$1000$ (either in the 1st or 2nd step).

$\bm{\Sigma}_{u}$ | $p$ | Estimator | 1st step | 2nd step | Time
---|---|---|---|---|---
|  |  | FPR | FNR | FPR | FNR | $\ell_{2}$ | 1st step | 2nd step
Diagonal | 1000 | One-stage Corrected | - | - | 20.0 (29.0) | 0.1 (1.5) | 0.96 (0.77) | 0.0 (0.0) | 18.9 (3.3)
|  | PMSc${}_{\text{CV}}$-Corrected | 100.0 (0.0) | 0.0 (0.0) | 19.7 (28.8) | 0.2 (2.0) | 0.95 (0.77) | 12.7 (0.7) | 18.7 (3.3)
|  | PMSc${}_{\text{FS}}$-Corrected | 7.6 (0.0) | 3.2 (7.8) | 1.2 (1.6) | 3.2 (7.8) | 0.56 (0.41) | 5.4 (0.7) | 1.4 (0.2)
|  | SISc-Corrected | 7.6 (0.1) | 6.7 (10.5) | 1.0 (0.9) | 6.7 (10.5) | 0.54 (0.30) | 1.3 (0.1) | 1.0 (0.2)
|  | One-stage CoCo | - | - | 11.5 (2.2) | 0.0 (0.0) | 0.77 (0.13) | 0.0 (0.0) | 481.5 (27.5)
|  | PMSc${}_{\text{CV}}$-CoCo | 100.0 (0.0) | 0.0 (0.0) | 11.5 (2.2) | 0.0 (0.0) | 0.77 (0.13) | 12.7 (0.7) | 476.7 (26.8)
|  | PMSc${}_{\text{FS}}$-CoCo | 7.6 (0.0) | 3.2 (7.8) | 1.0 (0.4) | 3.2 (7.8) | 0.45 (0.18) | 5.4 (0.7) | 1.8 (0.3)
|  | SISc-CoCo | 7.6 (0.1) | 6.7 (10.5) | 1.0 (0.4) | 6.7 (10.5) | 0.52 (0.23) | 1.3 (0.1) | 1.6 (0.2)
| 2000 | One-stage Corrected | - | - | 29.1 (21.0) | 1.1 (5.1) | 1.57 (0.84) | 0.0 (0.0) | 49.9 (7.0)
|  | PMSc${}_{\text{CV}}$-Corrected | 100.0 (0.0) | 0.0 (0.0) | 28.3 (21.3) | 1.0 (5.1) | 1.53 (0.84) | 25.2 (1.4) | 49.3 (6.9)
|  | PMSc${}_{\text{FS}}$-Corrected | 3.8 (0.0) | 4.0 (8.5) | 0.7 (0.9) | 4.0 (8.5) | 0.62 (0.51) | 8.5 (1.2) | 1.1 (0.2)
|  | SISc-Corrected | 3.8 (0.0) | 11.0 (13.1) | 0.5 (0.4) | 11.0 (13.1) | 0.61 (0.31) | 2.7 (0.1) | 0.9 (0.2)
|  | PMSc${}_{\text{CV}}$-CoCo | 100.0 (0.0) | 0.0 (0.0) | - | - | - | 25.2 (1.4) | -
|  | PMSc${}_{\text{FS}}$-CoCo | 3.8 (0.0) | 4.0 (8.5) | 0.5 (0.2) | 4.0 (8.5) | 0.45 (0.18) | 8.5 (1.2) | 1.7 (0.2)
|  | SISc-CoCo | 3.8 (0.0) | 11.0 (13.1) | 0.5 (0.3) | 11.0 (13.1) | 0.60 (0.29) | 2.7 (0.1) | 1.9 (0.2)
Block diagonal | 1000 | One-stage Corrected | - | - | 23.5 (29.0) | 1.1 (5.2) | 1.25 (0.73) | 0.0 (0.0) | 22.7 (4.5)
| PMSc${}_{\text{CV}}$-Corrected | 100.0 (0.0) | 0.0 (0.0) | 25.6 (29.4) | 0.8 (4.2) | 1.30 (0.74) | 14.5 (0.7) | 22.5 (4.6)
|  | PMSc${}_{\text{FS}}$-Corrected | 7.6 (0.0) | 3.4 (8.1) | 3.0 (3.2) | 3.6 (8.2) | 1.15 (0.69) | 5.3 (0.7) | 1.4 (0.3)
|  | SISc-Corrected | 7.6 (0.1) | 14.2 (14.8) | 1.5 (2.1) | 14.2 (14.9) | 0.97 (0.48) | 1.3 (0.1) | 1.1 (0.3)
|  | One-stage CoCo | - | - | 14.3 (2.3) | 0.4 (3.2) | 1.23 (0.17) | 0.0 (0.0) | 417.2 (28.2)
|  | PMSc${}_{\text{CV}}$-CoCo | 100.0 (0.0) | 0.0 (0.0) | 14.3 (2.2) | 0.5 (3.3) | 1.23 (0.17) | 14.5 (0.7) | 413.4 (28.5)
|  | PMSc${}_{\text{FS}}$-CoCo | 7.6 (0.0) | 3.4 (8.1) | 0.9 (0.4) | 3.6 (8.2) | 0.68 (0.18) | 5.3 (0.7) | 2.9 (0.3)
|  | SISc-CoCo | 7.6 (0.1) | 14.2 (14.8) | 0.9 (0.4) | 14.2 (14.8) | 0.84 (0.31) | 1.3 (0.1) | 2.3 (0.2)
| 2000 | One-stage Corrected | - | - | 32.4 (17.1) | 3.0 (10.1) | 1.86 (0.65) | 0.0 (0.0) | 52.5 (9.7)
|  | PMSc${}_{\text{CV}}$-Corrected | 100.0 (0.0) | 0.0 (0.0) | 31.8 (17.4) | 3.3 (10.1) | 1.84 (0.65) | 27.7 (1.5) | 52.1 (9.6)
|  | PMSc${}_{\text{FS}}$-Corrected | 3.8 (0.0) | 6.6 (10.8) | 1.6 (1.6) | 6.6 (10.9) | 1.21 (0.69) | 8.3 (1.1) | 1.1 (0.2)
|  | SISc-Corrected | 3.8 (0.0) | 21.4 (16.6) | 0.8 (1.0) | 21.4 (16.6) | 1.09 (0.48) | 2.8 (0.1) | 1.0 (0.3)
|  | PMSc${}_{\text{CV}}$-CoCo | 100.0 (0.0) | 0.0 (0.0) | - | - | - | 27.7 (1.5) | -
|  | PMSc${}_{\text{FS}}$-CoCo | 3.8 (0.0) | 6.6 (10.8) | 0.5 (0.2) | 6.7 (10.9) | 0.73 (0.22) | 8.3 (1.1) | 2.5 (0.2)
|  | SISc-CoCo | 3.8 (0.0) | 21.4 (16.6) | 0.5 (0.3) | 21.4 (16.6) | 0.97 (0.36) | 2.8 (0.1) | 2.5 (0.2)

Table D.3: Performance of the one-stage and two-stage estimators in the
simulation study based on false positive rate (FPR, in percentage), false
negative rate (FNR, in percentage), $\ell_{2}$ estimation error, and
computation time (in seconds) when $\bm{\Sigma}_{x}$ has an AR(1) structure
with autocorrelation $\rho_{x}=0.3$. Standard error are included in
parentheses. The CoCo estimator was only computed when the number of variables
was no more than $1000$ (either in the 1st or 2nd step).

$\bm{\Sigma}_{u}$ | $p$ | Estimator | 1st step | 2nd step | Time
---|---|---|---|---|---
|  |  | FPR | FNR | FPR | FNR | $\ell_{2}$ | 1st step | 2nd step
Diagonal | 1000 | One-stage Corrected | - | - | 2.5 (1.3) | 0.0 (0.0) | 0.35 (0.07) | 0.0 (0.0) | 31.5 (10.9)
|  | PMSc${}_{\text{CV}}$-Corrected | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.19 (0.07) | 23.4 (0.3) | 0.3 (0.0)
|  | PMSc${}_{\text{FS}}$-Corrected | 7.5 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.25 (0.07) | 5.3 (1.0) | 1.3 (0.3)
|  | SISc-Corrected | 7.5 (0.0) | 0.0 (0.0) | 1.2 (0.6) | 0.0 (0.0) | 0.31 (0.07) | 1.5 (0.0) | 1.8 (0.3)
|  | One-stage CoCo | - | - | 8.6 (1.9) | 0.0 (0.0) | 0.49 (0.08) | 0.0 (0.0) | 519.1 (33.7)
|  | PMSc${}_{\text{CV}}$-CoCo | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.19 (0.07) | 23.4 (0.3) | 0.0 (0.0)
|  | PMSc${}_{\text{FS}}$-CoCo | 7.5 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.26 (0.07) | 5.3 (1.0) | 0.6 (0.1)
|  | SISc-CoCo | 7.5 (0.0) | 0.0 (0.0) | 0.2 (0.1) | 0.0 (0.0) | 0.32 (0.07) | 1.5 (0.0) | 0.6 (0.1)
| 2000 | One-stage Corrected | - | - | 1.5 (0.7) | 0.0 (0.0) | 0.36 (0.08) | 0.0 (0.0) | 80.6 (30.8)
|  | PMSc${}_{\text{CV}}$-Corrected | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.19 (0.07) | 38.2 (1.6) | 0.2 (0.1)
|  | PMSc${}_{\text{FS}}$-Corrected | 3.8 (0.0) | 0.0 (0.0) | 0.6 (0.3) | 0.0 (0.0) | 0.26 (0.07) | 16.9 (2.9) | 2.0 (0.4)
|  | SISc-Corrected | 3.8 (0.0) | 0.0 (0.0) | 0.6 (0.3) | 0.0 (0.0) | 0.31 (0.08) | 2.8 (0.1) | 1.5 (0.3)
|  | PMSc${}_{\text{CV}}$-CoCo | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.19 (0.07) | 38.2 (1.6) | 0.0 (0.0)
|  | PMSc${}_{\text{FS}}$-CoCo | 3.8 (0.0) | 0.0 (0.0) | 0.1 (0.1) | 0.0 (0.0) | 0.27 (0.07) | 16.9 (2.9) | 0.6 (0.1)
|  | SISc-CoCo | 3.8 (0.0) | 0.0 (0.0) | 0.1 (0.1) | 0.0 (0.0) | 0.31 (0.07) | 2.8 (0.1) | 0.6 (0.1)
Block diagonal | 1000 | One-stage Corrected | - | - | 2.0 (0.8) | 0.0 (0.0) | 0.59 (0.12) | 0.0 (0.0) | 18.9 (7.4)
| PMSc${}_{\text{CV}}$-Corrected | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.36 (0.13) | 23.4 (0.4) | 0.3 (0.1)
|  | PMSc${}_{\text{FS}}$-Corrected | 7.5 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.2) | 0.45 (0.13) | 5.4 (1.1) | 1.6 (0.3)
|  | SISc-Corrected | 7.5 (0.0) | 0.0 (0.0) | 1.0 (0.5) | 0.0 (0.0) | 0.53 (0.14) | 1.5 (0.0) | 2.2 (0.4)
|  | PMSc${}_{\text{CV}}$-CoCo | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.36 (0.13) | 23.4 (0.4) | 0.0 (0.0)
|  | PMSc${}_{\text{FS}}$-CoCo | 7.5 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.3) | 0.46 (0.11) | 5.4 (1.1) | 1.3 (0.1)
|  | SISc-CoCo | 7.5 (0.0) | 0.0 (0.0) | 0.1 (0.0) | 0.0 (0.0) | 0.55 (0.12) | 1.5 (0.0) | 1.6 (0.2)
| 2000 | One-stage Corrected | - | - | 1.2 (0.5) | 0.0 (0.0) | 0.61 (0.13) | 0.0 (0.0) | 45.7 (23.4)
|  | PMSc${}_{\text{CV}}$-Corrected | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.37 (0.13) | 37.6 (1.6) | 0.3 (0.1)
|  | PMSc${}_{\text{FS}}$-Corrected | 3.8 (0.0) | 0.0 (0.0) | 0.5 (0.3) | 0.0 (0.0) | 0.45 (0.12) | 17.2 (3.0) | 2.4 (0.4)
|  | SISc-Corrected | 3.8 (0.0) | 0.0 (0.0) | 0.5 (0.3) | 0.0 (0.0) | 0.55 (0.13) | 2.8 (0.1) | 1.9 (0.3)
|  | PMSc${}_{\text{CV}}$-CoCo | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.0) | 0.0 (0.9) | 0.37 (0.17) | 37.6 (1.6) | 0.0 (0.0)
|  | PMSc${}_{\text{FS}}$-CoCo | 3.8 (0.0) | 0.0 (0.0) | 0.1 (0.0) | 0.1 (1.3) | 0.46 (0.12) | 17.2 (3.0) | 1.7 (0.2)
|  | SISc-CoCo | 3.8 (0.0) | 0.0 (0.0) | 0.1 (0.0) | 0.0 (0.0) | 0.56 (0.13) | 2.8 (0.1) | 1.4 (0.1)

Table D.4: Performance of the one-stage and two-stage estimators in the
simulation study based on false positive rate (FPR, in percentage), false
negative rate (FNR, in percentage), $\ell_{2}$ estimation error, and
computation time (in seconds) when $\bm{\Sigma}_{x}$ has an homogeneous
structure with $\rho_{x}=0.3$. Standard error are included in parentheses. The
CoCo estimator was only computed when the number of variables was no more than
$1000$ (either in the 1st or 2nd step).

$\bm{\Sigma}_{u}$ | $p$ | Estimator | 1st step | 2nd step | Time
---|---|---|---|---|---
|  |  | FPR | FNR | FPR | FNR | $\ell_{2}$ | 1st step | 2nd step
Diagonal | 1000 | One-stage Corrected | - | - | 20.0 (29.0) | 0.1 (1.5) | 0.96 (0.77) | 0.0 (0.0) | 18.9 (3.3)
|  | PMSc${}_{\text{CV}}$-Corrected | 100.0 (0.0) | 0.0 (0.0) | 19.7 (28.8) | 0.2 (2.0) | 0.95 (0.77) | 12.7 (0.7) | 18.7 (3.3)
|  | PMSc${}_{\text{FS}}$-Corrected | 7.5 (0.0) | 1.2 (4.8) | 3.8 (3.4) | 1.2 (4.8) | 1.12 (0.85) | 6.0 (1.1) | 0.9 (0.1)
|  | SISc-Corrected | 7.6 (0.1) | 6.7 (10.5) | 1.0 (0.9) | 6.7 (10.5) | 0.54 (0.30) | 1.3 (0.1) | 1.0 (0.2)
|  | One-stage CoCo | - | - | 11.5 (2.2) | 0.0 (0.0) | 0.77 (0.13) | 0.0 (0.0) | 481.5 (27.5)
|  | PMSc${}_{\text{CV}}$-CoCo | 100.0 (0.0) | 0.0 (0.0) | 11.5 (2.2) | 0.0 (0.0) | 0.77 (0.13) | 12.7 (0.7) | 476.7 (26.8)
|  | PMSc${}_{\text{FS}}$-CoCo | 7.5 (0.0) | 1.2 (4.8) | 0.8 (0.3) | 1.2 (4.8) | 0.36 (0.12) | 6.0 (1.1) | 0.9 (0.2)
|  | SISc-CoCo | 7.6 (0.1) | 6.7 (10.5) | 1.0 (0.4) | 6.7 (10.5) | 0.52 (0.23) | 1.3 (0.1) | 1.6 (0.2)
| 2000 | One-stage Corrected | - | - | 29.1 (21.0) | 1.1 (5.1) | 1.57 (0.84) | 0.0 (0.0) | 49.9 (7.0)
|  | PMSc${}_{\text{CV}}$-Corrected | 100.0 (0.0) | 0.0 (0.0) | 28.3 (21.3) | 1.0 (5.1) | 1.53 (0.84) | 25.2 (1.4) | 49.3 (6.9)
|  | PMSc${}_{\text{FS}}$-Corrected | 3.8 (0.0) | 1.1 (4.5) | 1.8 (1.7) | 1.1 (4.5) | 1.09 (0.86) | 10.1 (1.7) | 0.7 (0.1)
|  | SISc-Corrected | 3.8 (0.0) | 11.0 (13.1) | 0.5 (0.4) | 11.0 (13.1) | 0.61 (0.31) | 2.7 (0.1) | 0.9 (0.2)
|  | PMSc${}_{\text{CV}}$-CoCo | 100.0 (0.0) | 0.0 (0.0) | - | - | - | 25.2 (1.4) | -
|  | PMSc${}_{\text{FS}}$-CoCo | 3.8 (0.0) | 1.1 (4.5) | 0.4 (0.2) | 1.1 (4.5) | 0.33 (0.11) | 10.1 (1.7) | 0.9 (0.2)
|  | SISc-CoCo | 3.8 (0.0) | 11.0 (13.1) | 0.5 (0.3) | 11.0 (13.1) | 0.60 (0.29) | 2.7 (0.1) | 1.9 (0.2)
Block diagonal | 1000 | One-stage Corrected | - | - | 23.5 (29.0) | 1.1 (5.2) | 1.25 (0.73) | 0.0 (0.0) | 22.7 (4.5)
| PMSc${}_{\text{CV}}$-Corrected | 100.0 (0.0) | 0.0 (0.0) | 25.6 (29.4) | 0.8 (4.2) | 1.30 (0.74) | 14.5 (0.7) | 22.5 (4.6)
|  | PMSc${}_{\text{FS}}$-Corrected | 7.5 (0.0) | 0.9 (4.4) | 5.3 (3.2) | 0.9 (4.4) | 1.56 (0.74) | 5.9 (1.1) | 1.3 (0.4)
|  | SISc-Corrected | 7.6 (0.1) | 14.2 (14.8) | 1.5 (2.1) | 14.2 (14.9) | 0.97 (0.48) | 1.3 (0.1) | 1.1 (0.3)
|  | One-stage CoCo | - | - | 14.3 (2.3) | 0.4 (3.2) | 1.23 (0.17) | 0.0 (0.0) | 417.2 (28.2)
|  | PMSc${}_{\text{CV}}$-CoCo | 100.0 (0.0) | 0.0 (0.0) | 14.3 (2.2) | 0.5 (3.3) | 1.23 (0.17) | 14.5 (0.7) | 413.4 (28.5)
|  | PMSc${}_{\text{FS}}$-CoCo | 7.5 (0.0) | 0.9 (4.4) | 0.7 (0.4) | 0.9 (4.4) | 0.57 (0.15) | 5.9 (1.1) | 2.2 (0.2)
|  | SISc-CoCo | 7.6 (0.1) | 14.2 (14.8) | 0.9 (0.4) | 14.2 (14.8) | 0.84 (0.31) | 1.3 (0.1) | 2.3 (0.2)
| 2000 | One-stage Corrected | - | - | 32.4 (17.1) | 3.0 (10.1) | 1.86 (0.65) | 0.0 (0.0) | 52.5 (9.7)
|  | PMSc${}_{\text{CV}}$-Corrected | 100.0 (0.0) | 0.0 (0.0) | 31.8 (17.4) | 3.3 (10.1) | 1.84 (0.65) | 27.7 (1.5) | 52.1 (9.6)
|  | PMSc${}_{\text{FS}}$-Corrected | 3.8 (0.0) | 2.0 (6.3) | 2.8 (1.5) | 2.0 (6.3) | 1.63 (0.71) | 10.1 (1.7) | 0.9 (0.3)
|  | SISc-Corrected | 3.8 (0.0) | 21.4 (16.6) | 0.8 (1.0) | 21.4 (16.6) | 1.09 (0.48) | 2.8 (0.1) | 1.0 (0.3)
|  | One-stage CoCo | - | - | - | - | - | 0.0 (0.0) | -
|  | PMSc${}_{\text{CV}}$-CoCo | 100.0 (0.0) | 0.0 (0.0) | - | - | - | 27.7 (1.5) | -
|  | PMSc${}_{\text{FS}}$-CoCo | 3.8 (0.0) | 2.0 (6.3) | 0.4 (0.2) | 2.0 (6.3) | 0.59 (0.16) | 10.1 (1.7) | 1.9 (0.2)
|  | SISc-CoCo | 3.8 (0.0) | 21.4 (16.6) | 0.5 (0.3) | 21.4 (16.6) | 0.97 (0.36) | 2.8 (0.1) | 2.5 (0.2)