# Inference in High-Dimensional Linear Regression

H. S. Battey and N. Reid 111Heather Battey is Lecturer, Department of
Mathematics, Imperial College London, 180 Queen’s Gate, London SW7
(h.battey@imperial.ac.uk). Nancy Reid is University Professor of Statistical
Sciences, Department of Statistical Sciences, University of Toronto, 700
University Ave, 9th Floor, Toronto, Ontario M5G 1Z5, Canada
(nancym.reid@utoronto.ca).

(August 29, 2024)

Summary

We develop an approach to inference in a linear regression model when the
number of potential explanatory variables is larger than the sample size. Our
approach treats each regression coefficient in turn as the interest parameter,
the remaining coefficients being nuisance parameters, and seeks an optimal
interest-respecting transformation. The role of this transformation is to
allow a marginal least squares analysis for each variable, as in a factorial
experiment. One parameterization of the problem is found to be particularly
convenient, both computationally and mathematically. In particular, it permits
an analytic solution to the optimal transformation problem, facilitating
comparison to other work. In contrast to regularized regression such as the
lasso (Tibshirani, 1996) and its extensions, neither adjustment for selection,
nor rescaling of the explanatory variables is needed, ensuring the physical
interpretation of regression coefficients is retained. We discuss the use of
such confidence intervals as part of a broader set of inferential statements,
so as to reflect uncertainty over the model as well as over the parameters.
The considerations involved in extending the work to other regression models
are briefly discussed.

Keywords: confidence sets of models; nuisance parameters; orthogonality; post-
selection inference; high-dimensional regression; sparsity

##  1 Introduction

We consider inference for the coefficient parameters of a high-dimensional
linear regression model, with the number of potential explanatory variables
$p$ larger than the sample size $n$. Our approach is to treat each coefficient
in turn as a parameter of interest, and use an orthogonalization argument to
isolate this parameter in a simple linear regression. An advantage of our
approach is that inference statements about the parameter of interest are
straightforward, and do not require adjustment for selection. Regularized
regression is not used, freeing us from the need to standardize the columns of
the covariate matrix to have unit length, an operation that removes any
physical interpretation of the regression coefficients. Our approach is
inspired in part by the construction of supersaturated designs for factorial
experiments (Booth and Cox, 1962) and by parameter orthogonalization (Cox and
Reid, 1987).

The resulting confidence intervals are presented in the context of the
confidence sets of models described by Cox and Battey (2017). That work
emphasized the exploratory nature of the construction of plausible sets of
models but did not consider inferential statements about the individual
parameters. The use of these intervals as an adjunct to the identification of
confidence sets of models is illustrated on a set of gene-expression data from
Bühlmann et al. (2014), while simulations of high-dimensional models are
presented to assess the coverage properties of the proposed intervals and to
confirm the theoretical analysis.

Our approach differs from those in the literature on inference after
selection. Lockhart _et al._ (2014), Lee _et al._ (2014) and Tian and Taylor
(2018) correct inferential statements after variable selection by conditioning
on the region of the sample space that led to the relevant selection being
made. In contrast Zhang and Zhang (2014), van de Geer et al. (2014) and
Javanmard and Montanari (2014) construct hypothesis tests for all parameters
based on regularized regression and recommend corrections for multiple
testing. An overview of this approach is included in Bühlmann et al. (2014). A
third approach, advocated by Berk _et al._ (2013) and further discussed by
Leeb _et al._ (2015), aims for inferential statements that are universally
valid under all possible model selection procedures. We comment further on
related literature in §4.

##  2 Approximate orthogonalization

Let $x_{i}$ be a vector of measurements on $p$ covariates for individual $i$.
We consider the linear regression model

$Y_{i}=x_{i}^{T}\beta+\varepsilon_{i},\quad i=1,\ldots n,$

where $(\varepsilon_{i})_{i=1}^{n}$ are independent with mean zero and
variance $\tau>0$, and $\beta=(\beta_{1},\ldots,\beta_{p})^{T}$ is of
dimension $p\gg n$, but satisfying the sparsity condition
$\|\beta\|_{0}=\sum_{v=1}^{p}{\rm I1}\\{\beta_{v}\neq 0\\}=s\ll p$. The set of
signal indices is $\mathcal{S}=\\{v:\beta_{v}\neq 0\\}$. For notational
convenience, the arbitrary designated interest parameter will be denoted by
$\psi$, so that the vector $\lambda$ of nuisance parameters has elements
$\\{\beta_{1},\ldots,\beta_{p}\\}\backslash\psi$, an arbitrary component of
which we will denote by $\theta$. Let $Y=(Y_{1},\ldots,Y_{n})^{T}$ and let $X$
be the matrix whose $i$th row is $x_{i}^{T}$. The column of $X$ corresponding
to the interest parameter $\psi$ is written
$x_{\psi}=(x_{1,\psi},\ldots,x_{n,\psi})^{T}$ and the remaining columns are
those of the $n\times(p-1)$ dimensional matrix, $X_{-\psi}$, whose $i$th row
is $x^{T}_{i,-\psi}$. We use $x_{\theta}$ to represent an arbitrary column of
$X_{-\psi}$ so that $\sum_{\theta}x_{\theta}$ specifies summation over columns
corresponding to nuisance parameters without need for double indexing, the use
of Greek rather than Roman letters distinguishing the columns from the rows.
We also write, with a slight imprecision in notation,
$\sum_{\theta\in\mathcal{S}}x_{\theta}$ for summation over signal variables.
It is implicit from the specification of $\theta$ that this summation excludes
$x_{\psi}$ even if $x_{\psi}$ is a signal variable. On writing
$\varepsilon=(\varepsilon_{1},\ldots,\varepsilon_{n})^{T}$, the linear
regression model is

$Y=X\beta+\varepsilon=x_{\psi}\psi+X_{-\psi}\lambda+\varepsilon.$ (1)

If the columns in $X_{-\psi}$ are orthogonal to $x_{\psi}$ then a simple
regression of $Y$ on $x_{\psi}$ estimates $\psi$ without bias. In fact only
the columns corresponding to signal variables need be orthogonal to
$x_{\psi}$, but we do not know which columns these may be.

The parameter $\beta$ in (1) is unchanged by premultiplication of both sides
by an $n\times n$ matrix. Thus, for each interest parameter $\psi$, we propose
to seek a matrix $A^{\psi}$ such that the columns
$\widetilde{x}_{\theta}^{\psi}$ of the transformed covariate matrix
$\widetilde{X}^{\psi}=A^{\psi}X$ are as orthogonal as possible to
$\widetilde{x}_{\psi}^{\psi}$, and to estimate $\psi$ by simple linear
regression of $\widetilde{Y}^{\psi}=A^{\psi}Y$ on
$\widetilde{x}_{\psi}^{\psi}$. The associated least squares estimator is

$\widetilde{\psi}=(\widetilde{x}_{\psi}^{\psi
T}\widetilde{x}_{\psi}^{\psi})^{-1}\widetilde{x}_{\psi}^{\psi
T}\widetilde{Y}^{\psi}=(x_{\psi}^{T}A^{\psi
T}A^{\psi}x_{\psi})^{-1}x_{\psi}^{T}A^{\psi T}A^{\psi}Y.$ (2)

Under model (1)

$\displaystyle\mathbb{E}(\widetilde{\psi})$ $\displaystyle=$
$\displaystyle\psi+\textstyle{\sum_{\theta\in\mathcal{S}}}\vartheta_{\theta}\,\theta\stackrel{{\scriptstyle\text{def}}}{{=}}\psi+b_{\psi},$
(3) $\displaystyle\text{var}(\widetilde{\psi})$ $\displaystyle=$
$\displaystyle(\widetilde{x}_{\psi}^{\psi
T}\widetilde{x}_{\psi}^{\psi})^{-1}\widetilde{x}_{\psi}^{\psi
T}A^{\psi}\mathbb{E}(\varepsilon\varepsilon^{T})A^{\psi
T}\widetilde{x}_{\psi}^{\psi}(\widetilde{x}_{\psi}^{\psi
T}\widetilde{x}_{\psi}^{\psi})^{-1}\stackrel{{\scriptstyle\text{def}}}{{=}}\tau
v_{\psi\psi},$

where

$\displaystyle\vartheta_{\theta}=\vartheta_{\theta}(A^{\psi})$
$\displaystyle=$ $\displaystyle(\widetilde{x}_{\psi}^{\psi
T}\widetilde{x}_{\psi}^{\psi})^{-1}\widetilde{x}_{\psi}^{\psi
T}\widetilde{x}_{\theta}^{\psi},$ (4) $\displaystyle
v_{\psi\psi}=v_{\psi\psi}(A^{\psi})$ $\displaystyle=$
$\displaystyle(\widetilde{x}_{\psi}^{\psi
T}\widetilde{x}_{\psi}^{\psi})^{-2}\widetilde{x}_{\psi}^{\psi
T}A^{\psi}A^{\psi T}\widetilde{x}_{\psi}^{\psi},$

and by the Cauchy-Schwarz inequality,
$b_{\psi}^{2}\leq\|\lambda\|_{2}^{2}\sum_{\theta\in\mathcal{S}}\vartheta^{2}_{\theta}$.
Since we have no control over $\|\lambda\|_{2}$, and since $\mathcal{S}$ is
unknown, this suggests choosing $A^{\psi}$ to minimize
$v_{\psi\psi}(A^{\psi})+\sum_{\theta}\vartheta_{\theta}^{2}(A^{\psi})$, an
upper bound on the observable components of the mean squared error. To this
end, we reparameterize in terms of $q_{\psi}=A^{\psi T}A^{\psi}x_{\psi}$ and
write

$\displaystyle\widetilde{\psi}$ $\displaystyle=$
$\displaystyle(q_{\psi}^{T}x_{\psi})^{-1}q_{\psi}^{T}Y,$ (5) $\displaystyle
v_{\psi\psi}$ $\displaystyle=$
$\displaystyle(q_{\psi}^{T}x_{\psi})^{-2}(q_{\psi}^{T}q_{\psi})$
$\displaystyle\textstyle{\sum_{\theta}}\vartheta_{\theta}^{2}$
$\displaystyle=$
$\displaystyle(q_{\psi}^{T}x_{\psi})^{-2}\textstyle{\sum_{\theta}}(q_{\psi}^{T}x_{\theta})^{2}.$

On noting that

$(q_{\psi}^{T}x_{\theta})^{2}=q_{\psi}^{T}x_{\theta}x_{\theta}^{T}q_{\psi},$

the sum of the variance $\tau v_{\psi\psi}(q_{\psi})$ and squared potential
biasing terms $\textstyle{\sum_{\theta}}\vartheta_{\theta}^{2}(q_{\psi})$ is
minimized when $q_{\psi}$ solves the simple unconstrained optimization problem

$\mathop{\mathrm{argmin}}_{q\in\mathbb{R}^{n}}\;(q^{T}x_{\psi})^{-2}q^{T}(I_{n}+\textstyle{\sum_{\theta}}x_{\theta}x_{\theta}^{T})q.$

A generalization of this is helpful to draw parallels with other work. We
therefore define $q_{\psi}$ by

$q_{\psi}\in\mathcal{Q}_{\psi}^{(\delta)}=\\{\mathop{\mathrm{argmin}}_{q\in\mathbb{R}^{n}}\;(q^{T}x_{\psi})^{-2}q^{T}(\delta
I_{n}+\textstyle{\sum_{\theta}}x_{\theta}x_{\theta}^{T})q\\},$ (6)

for some fixed $\delta>0$ to be chosen. In §3, such solutions will be shown to
exist in closed form for any $\delta>0$ for which the matrix

$(\delta I_{n}+X_{-\psi}X_{-\psi}^{T})-\\{x_{\psi}^{T}(\delta
I_{n}+X_{-\psi}X_{-\psi}^{T})^{-1}x_{\psi}\\}^{-1}x_{\psi}x_{\psi}^{T}$

is non-negative definite, which can always be checked. As the notation in (6)
indicates, a solution is not unique; if
$q_{\psi}\in\mathcal{Q}_{\psi}^{(\delta)}$, then any non-zero scalar multiple
of $q_{\psi}$ is also in $\mathcal{Q}_{\psi}^{(\delta)}$. Hereafter, when we
refer to $q_{\psi}$, an arbitrary solution to the minimization problem of
equation (6) is understood unless otherwise specified.

There is a connection to parameter orthogonalization described by Cox and Reid
(1987) in that, for linear regression, both approaches produce zeros in the
relevant row and column of the Fisher information matrix. In the terminology
of Cox and Reid (1987), an interest-respecting orthogonal parameterization of
model (1) is $(\psi,\eta)$, where

$\eta\stackrel{{\scriptstyle\text{def}}}{{=}}X_{-\psi}^{T}x_{\psi}\psi+X_{-\psi}^{T}X_{-\psi}\lambda,$
(7)

or alternatively $(\psi,\phi)$, where

$\phi\stackrel{{\scriptstyle\text{def}}}{{=}}(X_{-\psi}^{T}X_{-\psi})^{-1}\eta=\lambda+(X_{-\psi}^{T}X_{-\psi})^{-1}X^{T}_{-\psi}x_{\psi}\psi,$
(8)

as orthogonal parameters are only defined up to linear transformation.
Expressing (1) in an orthogonal parametrization requires $p<n$.

If $(\varepsilon_{i})_{i=1}^{n}$ are normally distributed, a $(1-\alpha)$
confidence interval for $\psi$ is

$\mathcal{C}_{\psi}(\alpha)=[\widetilde{\psi}-z_{1-\alpha/2}(\tau
v_{\psi\psi})^{1/2},\widetilde{\psi}+z_{1-\alpha/2}(\tau
v_{\psi\psi})^{1/2}],$ (9)

where $\Phi(z_{\alpha})=\alpha$. Without distributional assumptions on the
error terms, the confidence interval $\eqref{eqCSPsi}$ is valid as
$n\rightarrow\infty$, to the extent that $q_{\psi}$ is successful in
orthogonalizing $\widetilde{x}^{\psi}_{\psi}$ to
$\widetilde{x}^{\psi}_{\theta}$ for all $\theta\in\mathcal{S}$. The
qualification here is because the objective function in (6) does not ensure
$\widetilde{x}^{\psi}_{\psi}$ is exactly orthogonal to
$\widetilde{x}^{\psi}_{\theta}$ for all $\theta$, only that the cumulative
non-orthogonality to $\widetilde{x}^{\psi}_{\psi}$ is small. In other words,
if $\sum_{\theta}\vartheta_{\theta}^{2}=m$ after approximate
orthogonalization, we would expect $\vartheta_{\theta}^{2}$ to be of order
$m/(p-1)$ on average over $\theta$. However for some $\psi$ this non-
orthogonality might be concentrated relatively more on the unknown signal
indices. One could in principle attempt to identify positions where the signal
variables are most likely to be, and adjust the objective function
accordingly, as in Zhang and Zhang (2014), which is discussed below in §4. We
do not pursue that route, since this spoils the simple inferential properties
of §3, and since Proposition 5 and our numerical work of §5 suggests such
modifications are unnecessary.

There is a simple transformation of $X$, $\mathring{X}^{\psi}$, say, whose
columns are exactly orthogonal to a given one, $\mathring{x}^{\psi}_{\psi}$
say. This is defined as $\mathring{X}^{\psi}=\mathring{P}^{\psi}X+C^{\psi}$
where
$\mathring{P}^{\psi}=I_{n}-(x_{\psi}^{T}x_{\psi})^{-1}x_{\psi}x_{\psi}^{T}$
projects the columns of $X$ into the subspace of $\mathbb{R}^{n}$ orthogonal
to $x_{\psi}$, and $C^{\psi}$ is a matrix of zeros apart from the column
$c^{\psi}_{\psi}=\mathring{x}^{\psi}_{\psi}=x_{\psi}$, which replaces the
single column of zeros produced by $\mathring{P}^{\psi}X$. Since this
transformation to $\mathring{X}^{\psi}$ is not of the form $A^{\psi}X$, it is
not interest-respecting. An interest-respecting version could be obtained by
projection of $\mathring{X}^{\psi}$, in a suitable matrix norm, onto the space
of matrices of the form $AX$, where $A\in\mathbb{R}^{n\times n}$. This is
however much less convenient than (6), both computationally and theoretically.

For any two parameters $\psi$ and $\chi$, the covariance between
$\widetilde{\psi}$ and $\widetilde{\chi}$ constructed according to equations
(2) and (6), is

$\text{cov}(\widetilde{\psi},\widetilde{\chi})=\tau(q_{\psi}^{T}x_{\psi})^{-1}q_{\psi}^{T}q_{\chi}(q_{\chi}^{T}x_{\chi})^{-1}\stackrel{{\scriptstyle\text{def}}}{{=}}\tau
v_{\psi\chi},$ (10)

which is known up to $\tau$.

We present next some theoretical aspects of the approximate orthogonalization
procedure, particularly the interpretation of $q_{\psi}$ and the associated
matrix $A^{\psi}$, which clarifies the relationship between our proposal and
that of Zhang and Zhang (2014). Proofs of the main results are given in the
Appendix.

##  3 Inferential properties

In the objective function

$m(q)=(q^{T}x_{\psi})^{-2}q^{T}(\delta I_{n}+X_{-\psi}X_{-\psi}^{T})q$ (11)

from equation (6), $\delta>0$ determines the relative importance of the
variance of $\widetilde{\psi}$ to its bias, which in turn determines the
relative importance of the length of each confidence interval to its coverage
probability.

Let $(\delta I_{n}+X_{-\psi}X_{-\psi}^{T})=M_{\delta}$, say. The relevance of
the matrix

$L_{\delta}\stackrel{{\scriptstyle\text{def}}}{{=}}M_{\delta}-\\{x_{\psi}^{T}M_{\delta}^{-1}x_{\psi}\\}^{-1}x_{\psi}x_{\psi}^{T}$

is explained immediately below Proposition 1, which specifies the set of
solutions $\mathcal{Q}_{\psi}^{(\delta)}$ from (6) in closed form.

######  Proposition 1.

Suppose $\delta>0$, bounded away from zero, is such that the eigenvalues of
$L_{\delta}$ are non-negative. Then any $q_{\psi}$ of the form

$q_{\psi}=a(\delta I_{n}+X_{-\psi}X_{-\psi}^{T})^{-1}x_{\psi}$ (12)

is a minimizer of (11), where $a$ is any non-zero real number. If the
eigenvalue condition is violated, then such a $q_{\psi}$ is a saddlepoint of
(11).

Modulo multiplication by a positive scalar, $L_{\delta}$ is the matrix field
$\nabla\nabla^{T}m$ evaluated at any $q_{\psi}$ of the form (12). Let
$u\stackrel{{\scriptstyle\text{def}}}{{=}}\\{x_{\psi}^{T}M_{\delta}^{-1}x_{\psi}\\}^{-1/2}x_{\psi}$
and $v\stackrel{{\scriptstyle\text{def}}}{{=}}-u$ so that
$L_{\delta}=M_{\delta}+uv^{T}$. Then

$\text{det}(L_{\delta})=(1+v^{T}M_{\delta}^{-1}u)\text{det}(M_{\delta})=(1-v^{T}M_{\delta}^{-1}v)\text{det}(M_{\delta})=0,$

where the first equality is the so-called matrix determinant lemma for
determinants of rank one perturbations and the last equality follows as
$v^{T}M_{\delta}^{-1}v=1$. Thus the matrix $L_{\delta}$ has at least one zero
eigenvalue, showing that, if the condition of Proposition 1 is satisfied, $m$
is locally weakly convex at any $q_{\psi}$ of the form (12) but not strictly
convex (e.g. Boyd and Vandenberghe, 2004). This implies that such solutions
minimize $m$ in a “valley” rather than at isolated local points.

It follows from $q_{\psi}=A^{\psi T}A^{\psi}x_{\psi}$ that

$A^{\psi T}A^{\psi}=a(\delta
I_{n}+X_{-\psi}X_{-\psi}^{T})^{-1}\stackrel{{\scriptstyle\text{def}}}{{=}}P^{\psi}(a,\delta).$
(13)

Although any choice of $a$ yields a minimizer of (11), the choice $a=\delta$
improves interpretability in view of equation (14) of Proposition 2. It is
also of theoretical value to consider $a=\delta\searrow 0$.

######  Proposition 2.

Set $a=\delta$ and define $X_{-\psi}^{+}$ to be the Moore-Penrose pseudo-
inverse of $X_{-\psi}$ given in (33) in Appendix A.2. Then

$P^{\psi}(\delta,\delta)=I_{n}-X_{-\psi}(\delta
I_{p-1}+X_{-\psi}^{T}X_{-\psi})^{-1}X_{-\psi}^{T}$ (14)

and

$P^{\psi}\stackrel{{\scriptstyle\text{\emph{def}}}}{{=}}\lim_{\delta\rightarrow
0^{+}}P^{\psi}(\delta,\delta)=I_{n}-X_{-\psi}X_{-\psi}^{+}$ (15)

is an orthogonal projection operator that maps any vector $u\in\mathbb{R}^{n}$
into the null space of $X_{-\psi}^{T}$. Provided that $x_{\psi}$ is not in the
column span of $X_{-\psi}$, it follows that

$\lim_{\delta\rightarrow
0^{+}}b_{\psi}(\delta)=(x_{\psi}^{T}P^{\psi}x_{\psi})^{-1}\lambda^{T}X_{-\psi}^{T}P^{\psi}x_{\psi}=0.$

This limiting case only serves to supply insight and is of low practical
relevance because $x_{\psi}$ will usually be in the column span of $X_{-\psi}$
when $p>n$. When $p<n$,
$X_{-\psi}^{+}=(X_{-\psi}^{T}X_{-\psi})^{-1}X_{-\psi}^{T}$ so that
$P^{\psi}=I_{n}-X_{-\psi}^{T}(X_{-\psi}^{T}X_{-\psi})^{-1}X_{-\psi}^{T}$ and
(14) can be viewed as a ridge regularized version of this projection (Hoerl
and Kennard, 1970). See §4 for further discussion.

The following Berry-Esseen bound for the statistic

$S_{n}(\psi)\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{(q_{\psi}^{T}x_{\psi})(\widetilde{\psi}-\psi-
b_{\psi})}{(q_{\psi}^{T}q_{\psi})^{1/2}}$ (16)

allows us to assess the approximate size and marginal coverage of the
confidence interval (9) when the error terms $(\varepsilon_{i})_{i=1}^{n}$ do
not follow a normal distribution.

######  Proposition 3.

Suppose that $(\varepsilon_{i})_{i=1}^{n}$ are independent mean zero random
variables of variance $\tau$ and bounded moments of order $2+\gamma$ for some
$\gamma>0$. Then

$\sup_{z\in\mathbb{R}}|\text{\emph{pr}}\\{\tau^{-1/2}S_{n}(\psi)\leq
z\\}-\Phi(z)|\rightarrow 0,\text{ as }n\rightarrow\infty.$

If the above condition holds with $\gamma=1$,

$\sup_{z\in\mathbb{R}}|\text{\emph{pr}}\\{\tau^{-1/2}S_{n}(\psi)\leq
z\\}-\Phi(z)|\leq Ce_{n},$ (17)

where $C\in(0.4097,0.5606]$, $(q_{j\psi})_{j=1}^{n}$ are the elements of
$q_{\psi}$, and

$e_{n}=\frac{\sum_{j=1}^{n}|q_{j\psi}|^{3}\mathbb{E}|\varepsilon_{j}|^{3}}{(n\tau)^{3/2}(q_{\psi}^{T}q_{\psi})^{3/4}}.$

With any estimator $V_{n}$ of the error variance that is consistent, in the
sense that $V_{n}/\tau$ converges in probability to one as
$n\rightarrow\infty$, the limiting normal approximation to the Studentized
statistic $V_{n}^{-1/2}S_{n}$ remains valid. This follows immediately from
Proposition 3 by Slutsky’s theorem. Proposition 4 provides more precise
quantification.

######  Proposition 4.

Suppose that the conditions of Proposition 3 hold with $\gamma=1$ and let
$V_{n}$ be any consistent estimator of $\tau$. On writing

$\text{\emph{pr}}(|V_{n}^{-1/2}-\tau^{-1/2}|>\delta_{n})\asymp
r_{n}(\delta_{n}),$ (18)

where $\delta_{n}$ and $r_{n}(\delta_{n})$ are both decreasing in $n$, the
normal approximation to the distribution of $V_{n}^{-1/2}S_{n}(\psi)$ has
uniform error of order

$\sup_{z\in\mathbb{R}}|\text{\emph{pr}}\\{V_{n}^{-1/2}S_{n}(\psi)\leq
z\\}-\Phi(z)|\leq
c\max\\{e_{n},1-\Phi(g_{n}),\delta_{n}g_{n},r_{n}(\delta_{n})\\},$ (19)

where $c$ is a constant that we do not quantify, $g_{n}$ is any nondecreasing
sequence such that $\delta_{n}g_{n}=o(1)$, and $e_{n}$ is defined in
Proposition 3.

While it is unclear which choice of $g_{n}$ delivers the sharpest bound in
(19), a reasonably tight bound is obtained by choosing $g_{n}$ such that
$1-\Phi(g_{n})$ is roughly the same order as $e_{n}$. By well known upper and
lower bounds on the standard normal survival function (e.g Lin and Bai, 2010,
p.10),

$1-\Phi(x)\asymp|x|^{-1}\phi(x)\asymp|x|^{-1}e^{-x^{2}}\lesssim|x|^{-1}e^{-|x|}.$

On setting $g_{n}^{-1}e^{-g_{n}}=e_{n}$ and solving approximately for $g_{n}$
using Lagrange’s inversion formula (e.g. de Bruijn, 1981, p.25–29), we obtain

$g_{n}=[\log(e_{n}^{-1})-\log\\{\log(e_{n}^{-1})\\}]+O[\log\\{\log(e_{n}^{-1})\\}/\log(e_{n}^{-1})],$

which is a very slowly growing function of $n$. Thus, asymptotically, the
confidence interval (9) with $\tau$ replaced by a consistent estimator $V_{n}$
has nominal coverage probability $1-\alpha$, and for any fixed $n$ the error
in this approximation is roughly of order
$\max\\{e_{n},\delta_{n}g_{n},r_{n}(\delta_{n})\\}$.

In the simulations of §5.1 we use (9), i.e. the standard normal approximation
for the distribution of the pivotal quantity

$U_{n}\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{(q_{\psi}^{T}x_{\psi})(\widetilde{\psi}-\psi)}{(V_{n}q_{\psi}^{T}q_{\psi})^{1/2}},$
(20)

ignoring the bias term $(q_{\psi}^{T}x_{\psi})b_{\psi}/(\tau
q_{\psi}^{T}q_{\psi})^{1/2}$. While $q_{\psi}$ is chosen to minimize the
observable part of an upper bound on $b_{\psi}^{2}$, it is inevitable that
$b_{\psi}$ will be larger for some components $\psi$ than for others, due to
lack of knowledge of $\mathcal{S}$. The implication of (20) is that the bias
plays a relatively larger role for components whose associated variances are
smaller. This is visually apparent from Figure 1 of §5. Proposition 5
establishes the behaviour of the bias $b_{\psi}$ as a function of the sample
size.

######  Proposition 5.

Provided that the non-zero elements of $\lambda$ are bounded, the bias
associated with $\widetilde{\psi}$ satisfies $b_{\psi}^{2}=O(s/n)$.

##  4 Related work

Our proposal is closest to that of Zhang and Zhang (2014), who modified the
score equation for the least squares estimator to
$z_{\psi}^{T}(Y-x_{\psi}\psi)=0,$ in our notation, with $z_{\psi}$ to be
chosen. By comparison the estimating equation that defines $\widetilde{\psi}$
is

$\widetilde{x}_{\psi}^{\psi
T}(\widetilde{Y}^{\psi}-\widetilde{x}_{\psi}^{\psi}\psi)=q_{\psi}^{T}(Y-x_{\psi}\psi)=0,\quad\text{with
}q_{\psi}=A^{\psi T}A^{\psi}x_{\psi}.$ (21)

Zhang and Zhang (2014) suggested a bias correction based an a pilot estimator
$\widehat{\theta}^{\text{(init)}}$ of each nuisance parameter $\theta$:

$\frac{x_{\psi}^{T}Y}{z_{\psi}^{T}x_{\psi}}-\sum_{\theta}\frac{\widehat{\theta}^{\text{(init)}}z_{\psi}^{T}x_{\theta}}{z_{\psi}^{T}x_{\psi}}=\psi+\frac{z_{\psi}^{T}\varepsilon}{z_{\psi}^{T}x_{\psi}}.$

Thus the bias is controlled by

$\max_{\theta}|z_{\psi}^{T}x_{\theta}|\cdot\|\widehat{\lambda}^{\text{(init)}}-\lambda\|_{1}$
(22)

and Zhang and Zhang (2014) recommended that
$\widehat{\lambda}^{\text{(init)}}$ be taken as the lasso estimator, so that
$\|\widehat{\lambda}^{\text{(init)}}-\lambda\|_{1}$ is small with high
probability under regularity conditions. They also recommended that $z_{\psi}$
be chosen using the lasso, motivated by ordinary least squares when $p<n$,
where the score equation $z_{\psi}^{T}(Y-x_{\psi}\psi)=0$ holds for the least
squares estimator with

$z_{\psi}=x_{\psi}^{\perp}\stackrel{{\scriptstyle\text{def}}}{{=}}(I_{n}-X_{-\psi}(X_{-\psi}^{T}X_{-\psi})^{-1}X_{-\psi})x_{\psi},$

i.e., the projection of $x_{\psi}$ on the orthogonal complement of the space
spanned by the columns of $X_{-\psi}$. On recalling that $q_{\psi}=A^{\psi
T}A^{\psi}x_{\psi}$, (14) in Proposition 2 shows the parallels between the two
approaches, although they were derived from different perspectives. In
particular, our operations on the sample space, designed to minimize a
feasible upper bound on the mean squared error, recover the ridge regression
residual vector for $z_{\psi}$, rather than the nodewise lasso residuals from
an $\ell_{1}$-constrained linear regression of $x_{\psi}$ on $X_{-\psi}$
recommended by Zhang and Zhang (2014).

An important advantage of our approach is that it does not require rescaling
the columns of $X$ to unit length, as is needed for regularized regression, so
the physical interpretation of the confidence intervals is maintained.

The conditions for valid inference using Zhang and Zhang’s debiased lasso are
more interpretable in the context of van de Geer et al. (2014), who
recommended it in a framework covering generalized linear models with
canonical link. Control of the bias is analogous to (22) but is made explicit
by viewing the matrix $X$ as random. We return to this aspect below. In linear
models the debiased lasso is proved to deliver small bias when a small number
of columns of $X$ are correlated with that corresponding to $\psi$. In more
general models the sparsity assumption is imposed on the Fisher information
matrix. In this respect, our approach seeks to achieve by construction what is
assumed by van de Geer et al. (2014), and is in some ways closer to Cox and
Reid (1987).

Ning and Liu (2017), Fang et al. (2017) and Shah and Bühlmann (2019) also
treat $X$ as random. Ning and Liu (2017) and Fang et al. (2017) consider high-
dimensional versions of score and likelihood ratio tests, assuming sparsity of
the $\psi$ column of the Fisher information matrix. Shah and Bühlmann (2019)
exploit ideas from the double robustness literature to relax stringent scaling
conditions, either on the sparsity of the Fisher information or on sparsity of
$\beta$, with a view to testing conditional independence of the response and a
given explanatory variable, conditional on the remaining variables. See their
extensive numerical study for an illustration of the properties of the
debiased lasso in a range of settings.

Our analysis is based on the conditional distribution of the response. For
inference about $\beta$ in the linear regression model $X$ is an ancillary
statistic, at least when $p<n$ (the definition of ancillarity for $p>n$ has
been questioned by Azriel (2019)). There are also substantive reasons for
preferring a conditional analysis. In observational studies $X$, although
random, is frequently realized in advance of the response $Y$, in which case
it may be possible to interpret the effect of $X$ on $Y$ as causal. For
example $X$ may be information on the genome or epi-genome, and $Y$ a measure
of some condition expected to be caused by changes in $X$. Even in cases where
the values of $(Y_{i},X_{i})_{i=1}^{n}$ are realized at the same time, it is
appropriate to condition on $X$ for studying its association with the outcome
variable $Y$. Otherwise it would be possible that an inference justified on
the basis of its unconditional properties performed poorly for any given
sample of $X$s.

##  5 Performance in simulations

###  5.1 Factorial experiment

In this section we assume the error variance is known; in §6.2 we remove this
assumption. Initially under highly controlled conditions, we calculate
coverage probabilities and lengths of the confidence intervals (9). Our Monte
Carlo experiment was designed as follows.

An $X$ matrix was obtained by generating each row $x_{i}^{T}$ for
$i=1,\ldots,n$ from a normal distribution of mean zero and covariance matrix
$\rho 1_{[p]}1_{[p]}^{T}+(1-\rho)I$. The means of each column were subtracted
so that the column sums of $X$ were zero in sample. Consistent with the
discussion above, $X$ was fixed across Monte Carlo replications. The values of
$n$, $p$ and $\rho$ were taken at all $2^{3}$ combinations of high and low
levels as in Table 1, the value of $n/p$ being the same in two of the four
configurations at high and low values of $\rho$. A sparse vector $\beta$,
consisting of ones in the first five entries and zeros elsewhere, was also
fixed across Monte Carlo replications. In each of 1000 replications, a vector
$\varepsilon$ of independent standard normally distributed errors was
generated anew and a vector of response variables constructed as
$Y=X\beta+\varepsilon$.

Each of the first 1000 entries of $\beta$ was treated in turn as the interest
parameter $\psi$ and the value of $q_{\psi}$ calculated by solving the
optimization problem (6) using formula (12) with $a=\delta=1$. Confidence
intervals at level $\alpha=0.05$ for of the first 1000 entries of $\beta$ were
constructed as in (9). Summary statistics are recorded in Table 1, where
coverage probabilities and lengths are averages over simulations. The median
and modal values reported are summaries across the first 1000 entries of
$\beta$.

|  |  | modal | median | median | median | 95th p.c.
---|---|---|---|---|---|---|---
$\rho$ | $n$ | $p$ | coverage | coverage | length | $\vartheta_{\theta}^{2}$ | $\vartheta_{\theta}^{2}$
0.9 | 70 | 2450 | 0.941 | 0.921 | 1.509 | 0.0065 | 0.056
0.9 | 70 | 1225 | 0.941 | 0.923 | 1.534 | 0.0063 | 0.055
0.9 | 35 | 2450 | 0.941 | 0.909 | 2.127 | 0.0135 | 0.120
0.9 | 35 | 1225 | 0.947 | 0.910 | 2.134 | 0.0133 | 0.117
0.1 | 70 | 2450 | 0.939 | 0.732 | 0.504 | 0.0065 | 0.056
0.1 | 70 | 1225 | 0.942 | 0.745 | 0.511 | 0.0063 | 0.055
0.1 | 35 | 2450 | 0.948 | 0.715 | 0.707 | 0.0134 | 0.118
0.1 | 35 | 1225 | 0.942 | 0.696 | 0.717 | 0.0133 | 0.118
Table 1: Modal and median Monte Carlo coverage probabilities and the median of
the Monte Carlo mean length of confidence intervals across the first 1000
coefficients. The nominal level of the coverage probabilities is
$1-\alpha=0.95$. The median and 95th percentile of the $1000\times(p-1)$
values of $\vartheta_{\theta}^{2}$ for the first 1000 coefficients are also
reported.

The $2^{3}$ factorial experiment just described allows estimation of main
effects and two-way interactions of $\rho$, $n$ and $p$. Table 2 reports
estimates of the main effects, for logistic regression of the coverage
probabilities, and linear regression of the length, using the values in Table
1. For the coverage probabilities, the multiplicative effects are reported,
with estimates close to one indicating a negligible effect on the coverage.
The effects of $\rho$, $n$ and $p$ on the modal coverage probabilities appear
negligible. The largest effect on the median coverage is due to $\rho$, for
which an increase from the low to the high level multiplies the odds of the
median confidence interval covering the corresponding interest parameter by
4.185. Increasing $n$ from its low level to its high level multiplies the odds
by 1.166. The median length of the confidence intervals increases with the
correlation, reflecting expression (5) for the variance of $\widetilde{\psi}$.
Increasing $n$ from its low to high level decreases the median length of the
confidence interval on average.

estimated main effect | $\rho$ | $n$ | $p$
---|---|---|---
modal coverage | $0.995$ | $0.933$ | $0.986$
median coverage | $4.185$ | $1.166$ | $1.005$
median length | $1.216$ | $-0.407$ | $-0.032$
Table 2: Estimated main effects of $\rho$, $n$ and $p$ on the modal and median
coverage probabilities and the median length of the confidence intervals
$\mathcal{C}_{\psi}(0.05)$. The estimated effects of coverage are
multiplicative on the odds scale; the length effect is additive.

More detailed summary information is provided in Figure 1, where we plot the
coverage probabilities in each of the eight settings as a function of the
absolute bias $|b_{\psi}|$. The $z$ axis indicated by the colour bar is the
mean length of the confidence intervals over Monte Carlo replication. As
expected in view of the discussion of §3, the coverage probabilities
associated with the $\rho=0.1$ settings decay more sharply with absolute bias
than those associated with the $\rho=0.9$ settings. The decrease in
$|b_{\psi}|$ for increasing $n$, as quantified in Proposition 5, is also
reflected in these plots.

$\rho=0.9$ $\rho=0.1$

$n$ high, $p$ high

$n$ high, $p$ low

$n$ low, $p$ high

$n$ low, $p$ low

Figure 1: Coverage probabilities (vertical axis) for each $\psi$ plotted
against the associated absolute bias $|b_{\psi}|$ (horizontal axis) and the
mean length of the confidence interval (colour axis). High and low values of
$p$ are $2450$ and $1225$ and those of $n$ are $70$ and $35$.

###  5.2 Simulations using real covariate data

To assess performance beyond the highly controlled conditions of the factorial
experiment, we replace the matrix $X$ from §5.1 with a matrix of covariate
data discussed by Bühlmann et al. (2014), retaining all other aspects of the
simulation study. The purpose of the investigation that generated these data
was to study the association of the response variable, vitamin $B_{2}$
production, with the logarithm of the expression levels of $4088$ genes. The
sample size is $n=71$ and the correlation among columns of $X$ ranges from
$-0.926$ to $0.991$. See Bühlmann et al. (2014) for a complete description of
these data with original references.

The modal and median simulated coverage probabilities over the $4088$
variables and the median of the simulated mean lengths are reported in Table
3. The coverage results are similar to those in Table 1 but the confidence
intervals are longer. This highlights the prospect for signal and noise
variables to be indistinguishable on the basis of such confidence intervals.
To reflect the fact that several models may be equally consistent with the
data, confidence sets of models are constructed in the real example of §6.2.
Our confidence intervals for the components of $\beta$ supplements those sets
of models.

|  |  |  |  | modal | median | median
---|---|---|---|---|---|---|---
min $\rho$ | max $\rho$ | mean $\rho$ | $n$ | $p$ | coverage | coverage | width
-0.926 | 0.991 | 0.029 | 71 | 4088 | 0.942 | 0.902 | 3.32
Table 3: Summary of data and properties of confidence intervals.

##  6 Confidence sets of models

###  6.1 General discussion

Confidence sets of models for high-dimensional regression were proposed by Cox
and Battey (2017). It was argued that if there are several low-dimensional
models that are equally compatible with the data, then one should aim to
report as many as is feasible.

The first phase of the construction is to identify a set of retained
explanatory variables, $\widehat{\mathcal{S}}$, where
$|\widehat{\mathcal{S}}|$ is reasonably large, but less than $n$. The model
that includes all these variables is called the encompassing model. A
confidence set of models, $\mathcal{M}$, is obtained by comparing the fit of
the encompassing model to the fit of any model using five or fewer variables
from $\widehat{\mathcal{S}}$, and retaining the models that predict the
response variable equally well. As these models are by construction on an
equal footing for the data at hand, choice between them requires either more
data or subject-matter expertise. The restriction to models of size five or
smaller is arbitrary and based on an assumption of sparsity.

The confidence intervals in §2 provide a method for refining this set of
models. If the estimated coefficients when fitting a model in $\mathcal{M}$
substantially exceed the limits of $\mathcal{C}_{\psi}(\alpha)$ from equation
(9), this casts doubt on the plausibility of that model. Indeed, the
conclusion could even be that none of the models in $\mathcal{M}$ is
consistent with the collection of confidence intervals
$\mathcal{C}_{\psi}(\alpha)$, suggesting that too severe a reduction has been
used for $\widehat{\mathcal{S}}$, i.e., $|\widehat{\mathcal{S}}|$ is too
small.

###  6.2 Numerical illustration

We use the full data described in §5.2, comprising observations on response
variables as well as covariates. The cross correlations of all potential
explanatory variables were examined and any adjacent variables whose
correlations exceeded 0.95 were treated as a single variable, with the first
of each such pair being retained in the analysis. Among the 4088 variables
initially present, 61 were discarded on these grounds, leaving 4027 variables
for further analysis.

The first analysis is the estimation of each of the 4027 regression
coefficients and their 95% confidence limits using the approximate
orthogonalization approach described in §2, with estimation of error variance
discussed below.

A separate analysis was carried out to choose the comprehensive model
$\widehat{\mathcal{S}}$ and a confidence set of models $\mathcal{M}$. The
variables in $\widehat{\mathcal{S}}$ are identified through a series of model
fits based on partially balanced incomplete block designs (Yates, 1936). Then
each possible low-dimensional sub-model of $\widehat{\mathcal{S}}$ is tested
against the encompassing model using an $F$-test at the $5\%$ level. Models
not rejected by the $F$-test comprise the confidence set $\mathcal{M}$. There
are 21946 models in $\mathcal{M}$, but only 48 of these give estimates of the
regression coefficients that are compatible with our 95% confidence intervals
for the individual parameters in Table 4. These 48 models are reported in
Table 5. Since the decision to report a model in Table 5 depends on the
simultaneous correctness of the associated confidence intervals, in principle
a Bonferroni type correction should be applied. As only 3-5 such intervals are
involved due to sparsity, the effective (Bonferroni-adjusted) levels of the
confidence intervals for the parameters are at most $0.25$. To guarantee
Bonferroni adjusted levels of $0.05$, marginal levels of $0.01$ would suffice.

We implemented a refinement of Cox and Battey’s method, as recent work in
Battey and Cox (2018) demonstrates that using all the available data to
identify first the encompassing model $\widehat{\mathcal{S}}$ and then the
confidence set of models $\mathcal{M}$ leads to lower than nominal coverage
for the confidence set of models. Thus we used a sample splitting procedure,
which serves to give a better calibrated set $\mathcal{M}$ but also enables
unbiased estimation of error variance for construction of the confidence
intervals for components of $\beta$. The sample of size $n=71$ observations
was randomly partitioned into subsets, $\mathcal{I}_{1}$ and
$\mathcal{I}_{2}$, say, of sizes $n_{1}=35$ and $n_{2}=36$. In the terminology
of Cox and Battey (2017), the encompassing model was identified using
$\mathcal{I}_{1}$ for the cube reduction and $\mathcal{I}_{2}$ for the square
reduction. The confidence set of models and estimate of error variance was
constructed using $\mathcal{I}_{1}$.

To improve the stability of this sample-splitting procedure, we repeated it
1000 times, re-randomizing the positions of the variable indices in the cube
each time. The subset $\widehat{\mathcal{S}}$ was taken to be all those
variables that survived the two-stage procedure in at least half of these
randomizations. The confidence intervals constructed as described in §2 for
the 22 variables in $\widehat{\mathcal{S}}$ are reported in Table 4. These
have been ordered according to the proportion of models from $\mathcal{M}$ to
which they belong. For comparison, the lasso and the elastic net, fitted to
all 71 observations and tuned to retain 22 variables, find nine and fourteen
of the same variables. These are indicated in the first column of Table 4.

From the variables constituting the models in Table 5, three were the first in
a string of two or more variables whose empirical correlations exceeded 0.95
(see the first paragraph of this subsection). As far as $\mathcal{M}$ is
concerned, such variables are essentially interchangeable. The affected
variables are 1278 (paired with 1279), 1285 (paired with 1286, 1287, 1288 and
1289) and 4002 (paired with 4003, 4004, 4005 and 4006).

variable index | proportion | $\widetilde{\psi}$ | lower limit | upper limit
---|---|---|---|---
1516L,E | 0.272 |  0.343 |  0.022 |  0.663
2564L,E | 0.272 | -1.481 | -1.801 | -1.160
1503L,E | 0.251 | -0.325 | -0.646 |  -0.0050
2138 | 0.249 | -0.062 | -0.382 |  0.259
4008E | 0.240 | -0.366 | -0.686 | -0.046
4002L,E | 0.240 | -0.505 | -0.825 | -0.185
1639L,E | 0.235 | -0.406 | -0.726 | -0.086
1603 | 0.228 | -1.048 | -1.368 | -0.728
403 | 0.225 |  0.902 |  0.582 |  1.223
3291 | 0.222 | -0.640 | -0.960 | -0.320
978 | 0.222 | -0.259 | -0.580 |  0.061
3808E | 0.222 |  0.677 |  0.356 |  0.997
1069E | 0.221 | -0.398 | -0.718 | -0.078
3514L,E | 0.217 |  1.373 |  1.053 |  1.694
1436E | 0.199 | -0.463 | -0.783 | -0.143
1278L,E | 0.190 |  0.147 | -0.173 |  0.467
1285 | 0.179 |  0.172 | -0.148 |  0.493
1303E | 0.179 |  0.187 | -0.133 |  0.507
1297L,E | 0.179 |  0.219 | -0.102 |  0.539
1423 | 0.176 |  0.043 | -0.277 |  0.363
1290 | 0.171 |  0.189 | -0.131 |  0.510
1312L,E | 0.156 |  0.490 |  0.169 |  0.810
Table 4: Indices for variables in $\widehat{\mathcal{S}}$, ordered by the proportion of models in $\mathcal{M}$ in which they are present. Upper and lower confidence bounds for the associated regression coefficient are calculated using (9) with $\tau$ estimated. The superscripts $L$ and $E$ in the first column indicate that the variable was also found by the lasso and the elastic net fitted to the full sample.  |  1516 |  2564 |  1503 |  2138 |  4008 |  4002 |  1639 |  1603 |  403 |  3291 |  978 |  3808 |  1069 |  3514 |  1436 |  1278 |  1285 |  1303 |  1297 |  1423 |  1290 |  1312
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
1 | - | - | - | - | - | - | * | * | - | - | - | - | * | - | - | - | - | - | - | - | - | -
2 | * | - | * | - | - | - | - | - | - | - | - | - | - | - | - | * | - | - | - | - | - | -
3 | * | - | - | - | - | - | - | - | - | - | - | - | - | - | * | - | * | - | - | - | - | -
4 | * | - | * | - | - | - | - | - | - | - | - | - | - | - | - | - | * | - | - | - | - | -
5 | * | - | - | - | - | - | - | - | - | - | - | - | - | - | * | - | - | - | - | - | * | -
6 | * | - | * | - | - | - | - | - | - | - | - | - | - | - | - | - | - | - | - | - | * | -
7 | * | - | - | - | - | - | - | - | - | - | - | - | - | - | * | - | - | - | * | - | - | -
8 | - | - | - | - | - | - | * | * | - | - | - | - | - | - | * | - | - | - | - | - | - | -
9 | * | * | - | - | - | * | - | - | - | - | - | - | - | - | - | - | - | - | - | - | - | -
10 | * | * | - | - | * | - | - | - | - | - | - | - | - | - | - | - | - | - | - | - | - | -
11 | - | - | - | - | * | - | * | - | - | * | - | - | - | - | - | - | - | - | - | - | - | -
12 | - | - | - | - | - | - | - | - | * | - | * | - | * | - | - | - | * | - | - | - | - | -
13 | - | - | - | - | - | - | - | - | * | - | * | - | * | - | - | - | - | - | * | - | - | -
14 | - | - | - | - | - | - | - | - | * | - | * | - | * | - | - | - | - | * | - | - | - | -
15 | - | - | - | - | - | - | - | - | * | - | * | - | * | - | - | - | - | - | - | - | - | *
16 | - | - | - | - | - | - | - | - | * | - | * | - | * | - | - | - | - | - | - | * | - | -
17 | - | - | - | - | - | - | - | - | * | - | - | - | * | - | - | - | - | - | * | * | - | -
18 | * | - | - | - | - | - | - | - | * | - | - | - | - | - | * | - | - | * | - | - | - | -
19 | - | - | * | - | - | * | - | - | - | - | - | * | * | - | - | - | - | - | - | - | - | -
20 | * | - | * | - | - | - | - | - | - | - | - | - | - | - | * | * | - | - | - | - | - | -
21 | * | - | * | - | - | - | - | - | - | - | - | * | - | - | - | * | - | - | - | - | - | -
22 | * | - | * | - | - | - | - | - | - | - | - | - | - | - | * | - | * | - | - | - | - | -
23 | * | - | * | - | - | - | - | - | - | - | - | * | - | - | - | - | * | - | - | - | - | -
24 | * | - | * | - | - | - | - | - | - | - | - | - | - | - | * | - | - | - | - | - | * | -
25 | * | - | * | - | - | - | - | - | - | - | - | * | - | - | - | - | - | - | - | - | * | -
26 | * | - | * | - | - | - | - | - | - | - | - | - | - | - | * | - | - | - | * | - | - | -
27 | * | - | * | - | - | - | - | - | - | - | - | * | - | - | - | - | - | - | * | - | - | -
28 | * | - | * | - | - | - | - | - | - | - | - | * | - | - | - | - | - | * | - | - | - | -
29 | - | - | - | - | - | * | - | * | - | - | - | * | - | - | * | - | - | - | - | - | - | -
30 | * | - | - | - | - | * | * | - | - | * | - | - | - | - | - | - | - | - | - | - | - | -
31 | * | * | - | - | * | * | - | - | - | - | - | - | - | - | - | - | - | - | - | - | - | -
32 | - | - | - | - | - | - | - | - | * | - | * | - | * | - | - | * | * | - | - | - | - | -
33 | - | - | - | - | - | - | - | - | * | - | * | - | * | - | - | * | - | - | * | - | - | -
34 | - | - | - | - | - | - | - | - | * | - | * | - | * | - | - | * | - | - | - | * | - | -
35 | - | - | - | - | - | - | - | - | * | - | * | - | * | - | - | - | * | - | * | - | - | -
36 | - | - | - | - | - | - | - | - | * | - | * | - | * | - | - | - | * | * | - | - | - | -
37 | - | - | - | - | - | - | - | - | * | - | * | - | * | - | - | - | * | - | - | * | - | -
38 | - | - | - | - | - | - | - | - | * | - | * | - | * | - | - | - | - | - | * | - | * | -
39 | - | - | - | - | - | - | - | - | * | - | * | - | * | - | - | - | - | * | - | - | * | -
40 | - | - | - | - | - | - | - | - | * | - | * | - | * | - | - | - | - | - | - | * | * | -
41 | - | - | - | - | - | - | - | - | * | - | * | - | * | - | - | - | - | * | * | - | - | -
42 | - | - | - | - | - | - | - | - | * | - | * | - | * | - | - | - | - | - | * | * | - | -
43 | - | - | - | - | - | - | - | - | * | - | * | - | * | - | - | - | - | * | - | * | - | -
44 | - | - | - | - | - | - | - | - | * | - | * | - | * | - | - | - | - | - | - | * | - | *
45 | * | - | - | - | * | - | - | - | - | - | - | - | - | - | * | - | - | - | - | * | * | -
46 | - | - | * | - | * | - | - | * | - | - | - | - | - | * | * | - | - | - | - | - | - | -
47 | - | - | * | - | - | * | * | - | - | - | - | * | - | - | * | - | - | - | - | - | - | -
48 | - | - | * | - | * | - | - | - | - | * | - | * | - | - | * | - | - | - | - | - | - | -
Table 5: The 48 models from $\mathcal{M}$ that are compatible with the
confidence intervals in Table 4. Asterisks represent presence of the variables
indicated.

##  7 Discussion and open problems

The logical argument used in §6 seems new and is of the following form.
Suppose for a (possible) contradiction that model $m$ is true. This leads to
standard least squares estimates $\hat{\beta}(m)$ of the low-dimensional
vector $\beta(m)$ associated with model $m$. Under model $m$, $\hat{\beta}(m)$
is unbiased and consistent for $\beta(m)$. Let $\psi$ be an arbitrary entry of
$\beta(m)$. To the extent that the orthogonalization makes $b_{\psi}=0$ our
confidence interval for $\psi$ is unbiased and consistent for any model. Thus
extremity of any entry of $\hat{\beta}(m)$, as calibrated by the associated
limits of $\mathcal{C}_{\psi}(\alpha)$, contradicts validity of $m$.

The extent to which the ideas of §2 apply beyond the linear regression model
remains unclear. We conclude by outlining potential routes to generalization.

In medical applications binary responses are particularly frequently
encountered. With such outcomes encoded arbitrarily as $\\{-1,1\\}$, a widely
used representation is the linear logistic form (Cox, 1958a, b) in which
$\log\\{\text{pr}(Y_{i}=1)/\text{pr}(Y_{i}=-1)\\}$ is assumed to depend
linearly on $x_{i}$. Battey, Cox and Jackson (2019) studied some of the
advantages and disadvantages of using a linear representation of the
probabilities, and showed that the coefficients of such a model should be
fitted by ordinary least squares rather than maximum likelihood, due to
extreme sensitivity of the latter to aberrant observations. The fitted
probabilities potentially exceed the $[0,1]$ interval even when the linear
representation is valid for all observations. Nevertheless, the approach has
several advantages, particularly with regard to interpretation of the
associated parameters.

For $y\in\\{-1,1\\}$ the linear in probability model is (Battey, Cox and
Jackson, 2019)

$\text{pr}(Y_{i}=y)=\frac{1}{2}(1+yx_{i}^{T}\beta),\quad|x_{i}^{T}\beta|<1$
(23)

implying $\mathbb{E}(Y_{i})=x_{i}^{T}\beta$. With $p<n$ the ordinary least
squares estimator $\widehat{\beta}_{\text{OLS}}=(X^{T}X)^{-1}X^{T}Y$ is
computable and unbiased, and remains so if $X$ and $Y$ are replaced by
$\widetilde{X}=AX$ and $\widetilde{Y}=AY$ for some $n\times n$ matrix $A$.
With this linear representation of the probabilities, ideas of the present
paper can be applied almost unchanged, with each parameter calculated as in
(2). A complication is that $\text{var}(Y_{i})=1-(x_{i}^{T}\beta)^{2}$ depends
on the unknown parameter, suggesting that an iterative approach may be needed,
or that the optimization problem (6) needs modification.

A somewhat related idea is to use the sequence of linear approximations to
maximum likelihood estimators provided by the Fisher scoring algorithm
(Fisher, 1925). For generalized linear models the associated derived response
and weights based on the score equation have a convenient canonical form and
the solutions to the iteratively reweighted least squares equations eventually
coincide with the maximum likelihood estimates. See Green (1984) for an
extensive discussion.

While the use of a Wald-type test statistic is reasonable for linear
regression, for other regression models there are compelling reasons, at least
in low dimensions, for preferring a likelihood ratio statistic. Indeed, the
Wald statistic is often justified on the basis that it locally approximates
the logarithm of the likelihood ratio. As noted by Cox and Snell (1989), the
likelihood ratio statistic preserves any natural asymmetry in the confidence
regions due to asymmetry of the log likelihood function. In simple contexts
there are also more sophisticated arguments based on approximate conditioning
(Cox, 1988) that point to its use. This suggests an approach based on direct
generalization of the ideas of Cox and Reid (1987). The possibility for this
is apparent from equation (8), where the factor multiplying $\psi$ is the
least squares estimate of the coefficient vector in a regression of $x_{\psi}$
on $X_{-\psi}$. The introduction of $\delta I_{n}$ would accommodate $p>n$ by
replacing this least squares estimate by a ridge regression estimate (Hoerl
and Kennard, 1970) and would yield a parameter approximately orthogonal to
$\psi$ when $p>n$. Further analysis is needed to understand how far this idea
generalizes. To be useful, any subsequent operation would have to exploit the
resulting sparsity in the relevant column of the Fisher information matrix as
in Ning and Liu (2017), who assumed sparsity of the Fisher information without
the preliminary approximate parameter-orthogonalizing step.

In §6 we used a simple estimator of the error variance, obtained as part of
the construction of calibrated confidence sets of models. A similar approach,
suggested by Fan et al. (2012), is refitted cross-validation based on a
bisection of the data and the use of any preliminary screening procedure that
recovers the true set of signal variables $\mathcal{S}$ with high probability.
Assuming the initial screening is successful, the distribution theory follows
straightforwardly using Cochran’s theorem (Cochran, 1934) for normally
distributed errors. For more general models the asymptotic properties were
established by Fan et al. (2012).

A natural question is whether our approximate orthogonalization from §2 can be
used to consistently estimate the error variance without sample splitting.
Assume for simplicity that the error terms $(\varepsilon_{i})_{i=1}^{n}$ are
normally distributed. If the approximate orthogonalization via equation 6 was
completely effective,

$\frac{(q_{\psi}^{T}x_{\psi})(\widetilde{\psi}-\psi)}{(q_{\psi}^{T}q_{\psi})^{1/2}}$
(24)

would be normally distributed of zero mean and variance $\tau$ for any
interest parameter $\psi\in\\{\beta_{1},\ldots,\beta_{p}\\}$. For most of the
$p$ parameters, it would follow that

$(q_{\psi}^{T}q_{\psi})^{-1/2}(q_{\psi}^{T}x_{\psi})\widetilde{\psi}$ (25)

is normally distributed of mean zero and variance $\tau$. The composite of $p$
such statistics is, however, contaminated by those corresponding to signal
variables and by those for which the optimization problem (6) fails to
orthogonalize $\widetilde{x}_{\psi}^{\psi}$ to $\widetilde{x}_{\theta}^{\psi}$
for one or more signal coefficients $\theta$. One might therefore seek to
construct a robust estimator of variance based on these $p$ statistics, but we
do not pursue this here.

Acknowledgements. We are grateful to Yanbo Tang for important comments on an
earlier draft. The work was partially supported by the UK Engineering and
Physical Sciences Research Council (EP/T01864X/1) and the Natural Sciences and
Engineering Council of Canada.

## References

  * Azriel (2019) Azriel, D. (2019). The conditionality principle in high-dimensional regression. Biometrika, 106, 702–707.
  * Battey and Cox (2018) Battey, H. S. and Cox, D. R. (2018). Large numbers of explanatory variables: a probabilistic assessment. Proc. R. Soc. Lond. A., 474, 20170631
  * Battey, Cox and Jackson (2019) Battey, H. S., Cox, D. R. and Jackson, M. V. (2019). On the linear in probability model for binary data. Royal Society Open Science, 6, 190067.
  * Berk _et al._ (2013) Berk, R., Brown, L., Buja, A., Zhang, K. and Zhao, L. (2013). Valid post-selection inference. Ann. Statist., 41, 802–837.
  * Booth and Cox (1962) Booth, K. H. V. and Cox, D. R. (1962). Some systematic supersaturated designs. Technometrics, 4, 489–495.
  * Boyd and Vandenberghe (2004) Boyd, S. and Vandenberghe, L. (2004). Convex Optimization. Cambridge University Press, Cambridge, UK.
  * de Bruijn (1981) de Bruijn, N. G. (1981). Asymptotic Methods in Analysis. Corrected reprint of the third edition. Dover Publications, New York.
  * Bühlmann et al. (2014) Bühlmann, P., Kalisch, M. and Meir, L. (2014). High-dimensional statistics with a view toward applications in biology. Annu. Rev. Stat. Appl., 1, 255–278.
  * Bühlmann and van de Geer (2011) Bühlmann, P. and van de Geer, S. (2011). Statistics for High-Dimensional Data. Springer-Verlag, Berlin.
  * Cochran (1934) Cochran, W. G. (1934). The distribution of quadratic forms in a normal system, with applications to the analysis of covariance. Proc. Camb. Phil. Soc., 30, 178–191.
  * Cox (1958a) Cox, D. R. (1958a). The regression analysis of binary sequences (with discussion). J. R. Statist. Soc. B, 20, 215–242.
  * Cox (1958b) Cox, D. R. (1958b). Two further applications of a model for binary regression. Biometrika, 45, 562–565.
  * Cox (1988) Cox, D. R. (1988). Some aspects of conditional and asymptotic inference: A review. Sankhya A, 50, 314–337.
  * Cox and Battey (2017) Cox, D. R. and Battey, H. S. (2017). Large numbers of explanatory variables, a semi-descriptive analysis. Proc. Nat. Acad. Sci., 114, 8592–8595.
  * Cox and Reid (1987) Cox, D. R. and Reid, N. (1987). Parameter orthogonality and approximate conditional inference (with discussion). J. R. Statist. Soc. B, 49, 1–39.
  * Cox and Snell (1989) Cox, D. R. and Snell, E. J. (1989). Analysis of Binary Data. Chapman and Hall, London.
  * Fan et al. (2012) Fan, J., Guo, S. and Hao, N. (2012). Variance estimation using refitted cross-validation in ultrahigh dimensional regression. J. R. Statist. Soc. B, 74, 37–65.
  * Fang et al. (2017) Fang, X. E., Ning, Y. and Liu, H. (2017). Testing and confidence intervals for high dimensional proportional hazards models. J. R. Statist. Soc. B, 79, 1415–1437.
  * Fisher (1925) Fisher, R. A. (1925). Theory of statistical estimation. Math. Proc. Cambridge. Phil. Soc., 22, 700–725.
  * Green (1984) Green, P. J. (1984). Iteratively reweighted least squares for maximum likelihood estimation, and some robust and resistant alternatives (with discussion). J. R. Statist. Soc. B, 46, 149–192.
  * Hoerl and Kennard (1970) Hoerl, A. E. and Kennard, R. W. (1970). Ridge regression: biased estimation for nonorthogonal problems. Technometrics, 12, 55–67.
  * Lin and Bai (2010) Lin, Z. and Bai, Z. (2010). Probability Inequalities. Science Press, Beijing.
  * Lee _et al._ (2014) Lee, J. D., Sun, D. L., Sun, Y. and Taylor, J. E. (2016). Exact post-selection inference, with application to the lasso. Ann. Statist., 44, 907–927.
  * Leeb _et al._ (2015) Leeb, H. and Pötscher, B. M. and Ewald, K. (2015). On various confidence intervals post-model-selection. Statist. Sci., 30, 216–227.
  * Lockhart _et al._ (2014) Lockhart, R., Taylor, J. E., Tibshirani, R. J. and Tibshirani, R. (2014). A significance test for the lasso. Ann. Statist., 42, 413–468.
  * Javanmard and Montanari (2014) Javanmard, A. and Montanari, A. (2014). Confidence intervals and hypothesis testing for high-dimensional regression. J. Mach. Learning. Res., 15, 2869–2909.
  * Ning and Liu (2017) Ning, Y. and Liu, H. (2017). A general theory of hypothesis tests and confidence regions for sparse high dimensional models. Ann. Statist., 45, 158–195.
  * Shah and Bühlmann (2019) Shah, R. D. and Bühlmann, P. (2019). Double-estimation-friendly inference for high-dimensional misspecified models. arXiv:1909.10828v1.
  * Tian and Taylor (2018) Tian, X. and Taylor, J. E. (2018). Selective inference with a randomized response. Ann. Statist., 46, 679–710.
  * Tibshirani (1996) Tibshirani, R. (1996). Regression shrinkage and selection via the LASSO J. R. Statist. Soc. B, 58, 267–288.
  * van de Geer et al. (2014) van de Geer, S., Bühlmann, P., Ritov, Y. and Dezeure, R. (2014). On asymptotically optimal confidence regions and tests for high-dimensional models. Ann. Statist., 42, 1166–1202.
  * Yates (1936) Yates, F. (1936). A new method of arranging variety trials involving a large number of varieties. J. Agric. Sci., 26, 424–455.
  * Zhang and Zhang (2014) Zhang, C-H. and Zhang, S. S. (2014). Confidence intervals for low dimensional parameters in high dimensional linear models. J. R. Statist. Soc. B, 76, 217–242.

##  Appendix A Proofs for Section 3

###  A.1 Proof of Proposition 1

###### Proof.

Write the objective function in equation (6) as $m=gh$, where
$g(q)=(q^{T}x_{\psi})^{-2}$ and $h(q)=q^{T}(\delta
I_{n}+X_{-\psi}X_{-\psi}^{T})q$. The associated gradient field evaluated at
$q$ is (e.g. Weatherburn, 1957, §41), $\nabla m=\nabla(gh)=g\nabla h+h\nabla
g$, or more explicitly

$\nabla m=2(q^{T}x_{\psi})^{-2}(\delta
I_{n}+X_{-\psi}X_{-\psi}^{T})q-2(q^{T}x_{\psi})^{-3}\\{q^{T}(\delta
I_{n}+X_{-\psi}X_{-\psi}^{T})q\\}x_{\psi}.$ (26)

Any stationary point, $q_{\psi}^{*}$ say, solves

$(q_{\psi}^{*T}x_{\psi})(\delta
I_{n}+X_{-\psi}X_{-\psi}^{T})q_{\psi}^{*}=\\{q_{\psi}^{*T}(\delta
I_{n}+X_{-\psi}X_{-\psi}^{T})q_{\psi}^{*}\\}x_{\psi}.$ (27)

Since the right hand side belongs to the one-dimensional subspace of
$\mathbb{R}^{n}$ spanned by $x_{\psi}$, we infer by inspection of the left
hand side that $(\delta I_{n}+X_{-\psi}X_{-\psi}^{T})q_{\psi}^{*}$ spans the
same subspace. Equivalently, there exists an $a\in\mathbb{R}\backslash\\{0\\}$
such that

$(\delta I_{n}+X_{-\psi}X_{-\psi}^{T})q_{\psi}^{*}=ax_{\psi}.$

The matrix $(\delta I_{n}+X_{-\psi}X_{-\psi}^{T})$ is invertible for any
values of $n$ and $p$ and any $\delta>0$ bounded away from zero, since the
minimum eigenvalue of the Gramian matrix $X_{-\psi}X_{-\psi}^{T}$ is zero and
for any eigenvalue-eigenvector pair for $X_{-\psi}X_{-\psi}^{T}$,
$(\gamma_{k},\xi_{k})$ say,

$(\delta I_{n}+X_{-\psi}X_{-\psi}^{T})\xi_{k}=(\delta+\gamma_{k})\xi_{k}.$

Thus write

$q_{\psi}^{*}=a(\delta I_{n}+X_{-\psi}X_{-\psi}^{T})^{-1}x_{\psi}.$ (28)

Upon substitution into equation (27) we verify that any such choice of $a$
supplies a valid solution to (27).

This set of solutions consists of local minimizers of $m$ if and only if
$\nabla\nabla^{T}m(q^{*}_{\psi})$ is non-negative definite for any
$q^{*}_{\psi}$ of the form given in equation (28). Equivalently
$v^{T}\\{\nabla\nabla^{T}m(q^{*}_{\psi})\\}v\geq 0$ for any
$v\in\mathbb{R}^{n}$. Otherwise the solutions $q_{\psi}^{*}$ represent other
critical points. Let $v\in\mathbb{R}^{n}$. The corresponding directional
derivative,
$\nabla_{v}m(q)\stackrel{{\scriptstyle\text{def}}}{{=}}v^{T}(\nabla m)$, of
$m$ at $q$ is

$\nabla_{v}m(q)=2(q^{T}x_{\psi})^{-2}v^{T}(\delta
I_{n}+X_{-\psi}X_{-\psi}^{T})q-2(q^{T}x_{\psi})^{-3}\\{q^{T}(\delta
I_{n}+X_{-\psi}X_{-\psi}^{T})q\\}v^{T}x_{\psi}.$ (29)

For $u,v\in\mathbb{R}^{n}$,

$\displaystyle\nabla_{u}\nabla_{v}m(q)=u^{T}\\{\nabla(v^{T}\nabla m)\\}$
$\displaystyle=$ $\displaystyle-4(q^{T}x_{\psi})^{-3}\\{v^{T}(\delta
I_{n}+X_{-\psi}X_{-\psi}^{T})q\\}(u^{T}x_{\psi})$ $\displaystyle+$
$\displaystyle 2(q^{T}x_{\psi})^{-2}u^{T}(\delta
I_{n}+X_{-\psi}X_{-\psi}^{T})v$ $\displaystyle+$ $\displaystyle
6(q^{T}x_{\psi})^{-4}\\{q^{T}(\delta
I_{n}+X_{-\psi}X_{-\psi}^{T})q\\}(v^{T}x_{\psi})(u^{T}x_{\psi})$
$\displaystyle-$ $\displaystyle 4(q^{T}x_{\psi})^{-3}\\{u^{T}(\delta
I_{n}+X_{-\psi}X_{-\psi}^{T})q\\}(v^{T}x_{\psi}),$

so that, on writing $M_{\delta}=(\delta I_{n}+X_{-\psi}X_{-\psi}^{T})$,

$\nabla_{v}\nabla_{v}m(q^{*}_{\psi})=2\\{x_{\psi}^{T}M_{\delta}^{-1}x_{\psi}\\}^{-2}v^{T}\bigl{[}M_{\delta}-\\{x_{\psi}^{T}M_{\delta}^{-1}x_{\psi}\\}^{-1}x_{\psi}x_{\psi}^{T}\bigr{]}v.$
(30)

The conclusion follows because, for any function
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}$,
$v^{T}\\{\nabla\nabla^{T}f(q)\\}v=\nabla_{v}\nabla_{v}f(q)$.

∎

###  A.2 Proof of Proposition 2

###### Proof.

In view of the well known matrix identity

$(A+BCD)^{-1}=A^{-1}-A^{-1}B\\{C^{-1}+DA^{-1}B\\}^{-1}DA^{-1},$

we see that

$A^{\psi T}A^{\psi}=a(\delta
I_{n}+X_{-\psi}X_{-\psi}^{T})^{-1}=a\delta^{-1}[I_{n}-X_{-\psi}(\delta
I_{p-1}+X_{-\psi}^{T}X_{-\psi})^{-1}X_{-\psi}^{T}].$ (31)

On setting $a=\delta$ and taking the limit as $\delta\rightarrow 0^{+}$,

$P^{\psi}=I_{n}-X_{-\psi}\lim_{\delta\rightarrow 0^{+}}(\delta
I_{p-1}+X_{-\psi}^{T}X_{-\psi})^{-1}X_{-\psi}^{T}=I_{n}-X_{-\psi}X_{-\psi}^{+},$
(32)

where $X_{-\psi}^{+}$ is the Moore-Penrose pseudo-inverse of $X_{-\psi}$, the
unique matrix satisfying the equations (Penrose, 1955):

$\displaystyle X_{-\psi}X_{-\psi}^{+}X_{-\psi}$ $\displaystyle=X_{-\psi}\hskip
42.67912pt(X_{-\psi}X_{-\psi}^{+})^{T}$ $\displaystyle=X_{-\psi}X_{-\psi}^{+}$
(33) $\displaystyle X_{-\psi}^{+}X_{-\psi}X_{-\psi}^{+}$
$\displaystyle=X_{-\psi}^{+}\hskip 42.67912pt(X_{-\psi}^{+}X_{-\psi})^{T}$
$\displaystyle=X_{-\psi}^{+}X_{-\psi}.$

It follows from these equations that $I_{n}-X_{-\psi}X_{-\psi}^{+}$ is
symmetric and equal to its square, from which we conclude that it is an
orthogonal projection operator onto the null space of $X_{-\psi}^{T}$.

The bias in (3) is

$\displaystyle b_{\psi}$ $\displaystyle=$ $\displaystyle(x_{\psi}^{T}A^{\psi
T}A^{\psi}x_{\psi})^{-1}\sum_{\theta\in\mathcal{S}}(x_{\psi}^{T}A^{\psi
T}A^{\psi}x_{\theta})\theta$ $\displaystyle=$
$\displaystyle(x_{\psi}^{T}A^{\psi
T}A^{\psi}x_{\psi})^{-1}a\lambda^{T}X_{-\psi}^{T}(\delta
I_{n}+X_{-\psi}X_{-\psi}^{T})^{-1}x_{\psi}.$

Recall that $\lambda=\\{\beta_{1},\ldots,\beta_{p}\\}\backslash\psi$. On
setting $a=\delta$ and taking limits we observe from equation (32) that

$\lim_{\delta\rightarrow 0^{+}}a\lambda^{T}X_{-\psi}^{T}(\delta
I_{n}+X_{-\psi}X_{-\psi}^{T})^{-1}x_{\psi}=\lambda^{T}X_{-\psi}^{T}(I_{n}-X_{-\psi}X_{-\psi}^{+})x_{\psi}=0,$

because $(I_{n}-X_{-\psi}X_{-\psi}^{+})x_{\psi}$ is in the null space of
$X_{-\psi}^{T}$. Since both limits exist, as $x_{\psi}$ is not in the column
span of $X_{-\psi}$ by the statement of the proposition, the limit of the
product is the product of the limits, and the conclusion follows on noting
that $x_{\psi}^{T}(I_{n}-X_{-\psi}X_{-\psi}^{+})x_{\psi}\neq 0$. ∎

Remark. If $x_{\psi}$ is in the column span of $X_{-\psi}$ there exists a
$(p-1)$-dimensional vector $b$ such that $x_{\psi}^{T}=b^{T}X_{-\psi}^{T}$ and
the limit of $b_{\psi}$ is an indeterminate form, perhaps calculable using
L’Hôpital’s rule.

###  A.3 Proof of Proposition 3

###### Proof.

Let $q_{j\psi}$ denote the $j$th entry of $q_{\psi}$. Thus write $S_{n}(\psi)$
as $\sum_{j=1}^{n}\xi_{j}$, where

$\xi_{j}:=(q_{\psi}^{T}q_{\psi})^{-1/2}q_{j\psi}\varepsilon_{j}.$

We have $\mathbb{E}(\xi_{j})=0$ and $\text{var}(\sum_{j}\xi_{j})=\tau$. To
establish a Lindeberg-Feller central limit theorem (e.g. Kallenberg, 1997,
Theorem 4.12) for $\sum_{j}\xi_{j}$, we need only verify the Lindeberg
condition, i.e., for all $\eta>0$,

$\lim_{n\rightarrow\infty}\tau^{-1}\sum_{j=1}^{n}\mathbb{E}(\xi_{j}^{2}\;{\rm
I1}\\{|\xi_{j}|>\eta\tau^{1/2}\\})=0.$

Let $c_{n}$ denote a lower bound on $q_{\psi}^{T}q_{\psi}$. Note that $c_{n}$
is bounded away from zero because otherwise the evaluation of the objective
function in (6) at $q_{\psi}$ is (positive and) unbounded, contradicting the
definition of $q_{\psi}$ as a minimizer. Let
$\vartheta_{n}=\|q_{\psi}\|_{\max}$. Then $|\xi_{j}|\leq
c_{n}^{-1}\vartheta_{n}|\varepsilon_{j}|$ and the event
$\\{|\xi_{j}|>\eta\tau^{1/2}\\}$ is contained in the event
$\\{|\varepsilon_{j}|>\eta\tau^{1/2}c_{n}\vartheta_{n}^{-1}\\}$. It follows
that

$\displaystyle\tau^{-1}\sum_{j=1}^{n}\mathbb{E}(\xi_{j}^{2}\;{\rm
I1}\\{|\xi_{j}|>\eta\tau^{1/2}\\})$ $\displaystyle\leq$
$\displaystyle\tau^{-1}\sum_{j=1}^{n}\mathbb{E}(\xi_{j}^{2}\;{\rm
I1}\\{|\varepsilon_{j}|>\eta\tau^{1/2}c_{n}\vartheta_{n}^{-1}\\})$
$\displaystyle=$
$\displaystyle\tau^{-1}(q_{\psi}^{T}q_{\psi})^{-1}\sum_{j=1}^{n}q_{j\psi}^{2}\mathbb{E}(\varepsilon_{j}^{2}\;{\rm
I1}\\{|\varepsilon_{j}|>\eta\tau^{1/2}c_{n}\vartheta_{n}^{-1}\\})$
$\displaystyle=$ $\displaystyle\tau^{-1}\mathbb{E}(\varepsilon_{j}^{2}\;{\rm
I1}\\{|\varepsilon_{j}|>\eta\tau^{1/2}c_{n}\vartheta_{n}^{-1}\\}).$

Let $\delta=\eta\tau^{1/2}c_{n}\vartheta_{n}^{-1}$. For any $\gamma>0$,

$\mathbb{E}(\varepsilon_{j}^{2}\;{\rm
I1}\\{|\varepsilon_{j}|>\delta\\})\leq\mathbb{E}(\varepsilon_{j}^{2}\delta^{-\gamma}|\varepsilon_{j}|^{\gamma}\;{\rm
I1}\\{|\varepsilon_{j}|>\delta\\})\leq\delta^{-\gamma}\mathbb{E}(|\varepsilon_{j}|^{2+\gamma}).$

Thus, since $\mathbb{E}(|\varepsilon_{j}|^{2+\gamma})<\infty$ and
$\vartheta_{n}^{\gamma}=\|q_{\psi}\|_{\max}^{\gamma}=o(\|q_{\psi}\|_{2}^{\gamma})\asymp
c_{n}^{\gamma}$,

$\lim_{n\rightarrow\infty}\tau^{-1}\sum_{j=1}^{n}\mathbb{E}(\xi_{j}^{2}\;{\rm
I1}\\{|\xi_{j}|>\eta\tau^{1/2}\\})=\lim_{n\rightarrow\infty}(\eta\tau^{1/2}c_{n}\vartheta_{n}^{-1})^{-\gamma}\mathbb{E}(|\varepsilon_{j}|^{2+\gamma})=0$

and the Lindeberg condition is verified. The fact that
$\|q_{\psi}\|_{\max}^{\gamma}=o(\|q_{\psi}\|_{2}^{\gamma})$ for any $\gamma>0$
is evident from the expression (12), which ensures that the number of non-zero
elements of $q_{\psi}$ increases with $n$.

Equation (17) is an application of the Berry-Esseen Theorem for non-
identically distributed random variables (e.g. Petrov, 1995, Chapter 5) with
the constant $C$ improved by Tyurin (2010). ∎

###  A.4 Proof of Proposition 4

###### Proof.

The second term in the decomposition

$\displaystyle\sup_{z\in\mathbb{R}}|\text{pr}\\{V_{n}^{-1/2}S_{n}(\psi)\leq
z\\}-\Phi(z)|$ $\displaystyle\leq$
$\displaystyle\sup_{z\in\mathbb{R}}|\text{pr}\\{V_{n}^{-1/2}S_{n}(\psi)\leq
z\\}-\text{pr}\\{\tau^{-1/2}S_{n}(\psi)\leq z\\}|$ (34) $\displaystyle+$
$\displaystyle\sup_{z\in\mathbb{R}}|\text{pr}\\{\tau^{-1/2}S_{n}(\psi)\leq
z\\}-\Phi(z)|\stackrel{{\scriptstyle\text{def}}}{{=}}\text{I}+\text{II},$

is bounded by $Ce_{n}$ using Proposition 3.

For any $\delta>0$, define the events

$\displaystyle\mathcal{E}_{1}(\delta)$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}$
$\displaystyle\\{(\tau^{-1/2}-V_{n}^{-1/2})S_{n}(\psi)\leq\delta\\},$
$\displaystyle\mathcal{E}_{2}(\delta)$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}$
$\displaystyle\\{|(V_{n}^{-1/2}-\tau^{-1/2})S_{n}(\psi)|\leq\delta\\}.$

On the event $\\{V_{n}^{-1/2}S_{n}(\psi)\leq z\\}\cap\mathcal{E}_{1}(\delta)$,

$\tau^{-1/2}S_{n}(\psi)-\delta\leq V_{n}^{-1/2}S_{n}(\psi)\leq z.$

Thus $\\{V_{n}^{-1/2}S_{n}(\psi)\leq
z\\}\cap\mathcal{E}_{1}(\delta)\subseteq\\{\tau^{-1/2}S_{n}(\psi)\leq
z+\delta\\}$, and it follows on writing

$\\{V_{n}^{-1/2}S_{n}(\psi)\\}=\left\\{\\{V_{n}^{-1/2}S_{n}(\psi)\\}\cap\mathcal{E}_{1}(\delta)\right\\}\cup\left\\{\\{V_{n}^{-1/2}S_{n}(\psi)\\}\cap\mathcal{E}_{1}^{c}(\delta)\right\\}$

that

$\text{pr}\\{V_{n}^{-1/2}S_{n}(\psi)\leq
z\\}\leq\text{pr}\\{\tau^{-1/2}S_{n}(\psi)\leq
z+\delta\\}+\text{pr}\\{\mathcal{E}_{1}^{c}(\delta)\\},$ (35)

where $\mathcal{E}_{1}^{c}(\delta)$ denotes the complement of
$\mathcal{E}_{1}(\delta)$. Later $\delta$ will be replaced by a decreasing
function of $n$, thus control over I in (34) requires an asymptotically
matching lower bound on $\text{pr}\\{V_{n}^{-1/2}S_{n}(\psi)\leq z\\}$ as
$\delta\rightarrow 0$. To this end, consider

$V_{n}^{-1/2}S_{n}(\psi)\leq\tau^{-1/2}S_{n}(\psi)+|(V_{n}^{-1/2}-\tau^{-1/2})S_{n}(\psi)|,$

yielding the containment

$\\{V_{n}^{-1/2}S_{n}(\psi)\leq
z\\}\supseteq\\{\tau^{-1/2}S_{n}(\psi)+|(V_{n}^{-1/2}-\tau^{-1/2})S_{n}(\psi)|\leq
z\\}\stackrel{{\scriptstyle\text{def}}}{{=}}\mathcal{A}.$

Furthermore,

$\mathcal{A}\supseteq\mathcal{A}\cap\mathcal{E}_{2}(\delta)\supseteq\\{\tau^{-1/2}S_{n}(\psi)\leq
z-\delta\\}\cap\mathcal{E}_{2}(\delta),$

so that

$\text{pr}\\{V_{n}^{-1/2}S_{n}(\psi)\\}\geq\text{pr}\\{\tau^{-1/2}S_{n}(\psi)\leq
z-\delta\\}\text{pr}\\{\mathcal{E}_{2}(\delta)\\}$ (36)

and equations (35) and (36) together give

$\displaystyle\sup_{z\in\mathbb{R}}|\text{pr}\\{V_{n}^{-1/2}S_{n}(\psi)\leq
z\\}-\text{pr}\\{\tau^{-1/2}S_{n}(\psi)\leq z\\}|$ $\displaystyle\leq$
$\displaystyle\sup_{z^{\prime}:|z-z^{\prime}|\leq\delta}\max\left\\{|\text{pr}\\{\tau^{-1/2}S_{n}(\psi)\leq
z^{\prime}\\}\text{pr}\\{\mathcal{E}_{2}(\delta)\\}-\right.\text{pr}\\{\tau^{-1/2}S_{n}(\psi)\leq
z\\}|,$ $\displaystyle\hskip
71.13188pt\left.|\text{pr}\\{\tau^{-1/2}S_{n}(\psi)\leq
z^{\prime}\\}-\text{pr}\\{\tau^{-1/2}S_{n}(\psi)\leq
z\\}|+\text{pr}\\{\mathcal{E}_{1}^{c}(\delta)\\}\right\\}$

for any $\delta>0$. Additionally, by Proposition 3 and a Taylor Series
expansion of the standard normal distribution function,

$\displaystyle\sup_{z,z^{\prime}:|z-z^{\prime}|\leq\delta}|\text{pr}\\{\tau^{-1/2}S_{n}(\psi)\leq
z^{\prime}\\}-\text{pr}\\{\tau^{-1/2}S_{n}(\psi)\leq z\\}|$
$\displaystyle\asymp$
$\displaystyle\sup_{z,z^{\prime}:|z-z^{\prime}|\leq\delta}|\Phi(z^{\prime})-\Phi(z)|+e_{n}\asymp\delta+e_{n}.$

As a function of $\delta\rightarrow 0$, we have

$\text{pr}\\{\mathcal{E}_{1}^{c}(\delta)\\}\asymp\text{pr}\\{|(V_{n}^{-1/2}-\tau^{-1/2})S_{n}(\psi)|\leq\delta\\}=1-\text{pr}\\{\mathcal{E}_{2}(\delta)\\}$

so that $\text{pr}\\{\mathcal{E}_{2}(\delta)\\}\rightarrow 1$ at the same rate
in $\delta$ as $\text{pr}\\{\mathcal{E}_{1}^{c}(\delta)\\}\rightarrow 0$. Let
$d_{n}$ be a decreasing function and $g_{n}$ a non-decreasing function, both
to be specified. For two arbitrary constants $c_{1}$ and $c_{2}$, consider the
event

$\mathcal{D}\stackrel{{\scriptstyle\text{def}}}{{=}}\\{|(V_{n}^{-1/2}-\tau^{-1/2})||S_{n}(\psi)|>c_{1}c_{2}d_{n}g_{n}\\}=\\{\mathcal{D}\cap\mathcal{S}\\}\cup\\{\mathcal{S}\cap\mathcal{S}^{c}\\},$

where

$\mathcal{S}=\\{|S_{n}(\psi)|\leq
c_{2}g_{n}\\}=\\{|\tau^{-1/2}S_{n}(\psi)|\leq\tau^{-1/2}c_{2}g_{n}\\}.$

On $\mathcal{S}$, $|S_{n}(\psi)/c_{2}g_{n}|\leq 1$, yielding
$\mathcal{D}\cap\mathcal{S}\subseteq\\{|V_{n}^{-1/2}-\tau^{-1/2}|>c_{1}d_{n}\\}$.
Thus

$\text{pr}\\{|(V_{n}^{-1/2}-\tau^{-1/2})S_{n}(\psi)|>c_{1}c_{2}d_{n}g_{n}\\}\leq\text{pr}(|V_{n}^{-1/2}-\tau^{-1/2}|>c_{1}d_{n})+\text{pr}(\mathcal{S}^{c}),$
(38)

where $\text{pr}(\mathcal{S}^{c})\asymp\\{1-\Phi(g_{n})\\}+e_{n}$ by
Proposition 3. The conclusion follows by combining (34), (A.4) and (38). ∎

###  A.5 Proof of Proposition 5

###### Proof.

On substituting the expression (12) for $q_{\psi}$ in the formula for
$\sum_{\theta}\vartheta_{\theta}^{2}$ in equation (2), we obtain

$\displaystyle\textstyle{\sum_{\theta}}\vartheta_{\theta}^{2}$
$\displaystyle=$ $\displaystyle\\{x_{\psi}^{T}(\delta
I_{n}+X_{-\psi}X_{-\psi}^{T})^{-1}x_{\psi}x_{\psi}^{T}(\delta
I_{n}+X_{-\psi}X_{-\psi}^{T})^{-1}x_{\psi}\\}^{-1}$
$\displaystyle\cdot\;\;\\{x_{\psi}^{T}(\delta
I_{n}+X_{-\psi}X_{-\psi}^{T})^{-1}X_{-\psi}X_{-\psi}^{T}(\delta
I_{n}+X_{-\psi}X_{-\psi}^{T})^{-1}x_{\psi}\\}$ $\displaystyle=$
$\displaystyle\\{x_{\psi}^{T}P^{\psi}(\delta,\delta)x_{\psi}x_{\psi}^{T}P^{\psi}(\delta,\delta)x_{\psi}\\}^{-1}\\{x_{\psi}^{T}P^{\psi}(\delta,\delta)X_{-\psi}X_{-\psi}^{T}P^{\psi}(\delta,\delta)x_{\psi}\\},$

where $P^{\psi}(\delta,\delta)$ is defined in equation (14).

The $(p-1)$-dimensional vector

$u=X_{-\psi}^{T}P^{\psi}(\delta,\delta)x_{\psi}$

has entries close to zero as $P^{\psi}(\delta,\delta)x_{\psi}$ is a
regularized projection of $x_{\psi}$ onto the null space of $X_{-\psi}^{T}$.
More importantly, when $(p-1)>n$, there exists a $(p-1)$-dimensional vector
$b$, $(p-1)-n$ entries of which are zero, such that $x_{\psi}=X_{-\psi}b$.
Thus write

$\textstyle{\sum_{\theta}}\vartheta_{\theta}^{2}=\|u\|_{2}^{2}/(u^{T}bb^{T}u)=O(p/n).$

It follows that $\sum_{\theta\in\mathcal{S}}\vartheta_{\theta}^{2}=O(s/n)$,
while by the Cauchy-Schwarz inequality
$b_{\psi}^{2}\leq\|\lambda\|_{2}^{2}\sum_{\theta\in\mathcal{S}}\vartheta_{\theta}^{2}$.
∎

## References

  * Kallenberg (1997) Kallenberg, O. (1997). Foundations of Modern Probability. Springer, Berlin.
  * Penrose (1955) Penrose, R. (1955). A generalized inverse for matrices. Proc. Cambridge Philos. Soc., 51, 406–413.
  * Petrov (1995) Petrov, V. V. (1995). Limit Theorems of Probability Theory. Oxford University Press, Oxford.
  * Tyurin (2010) Tyurin, I. S. (2010). An improvement of upper estimates of the constants in the Lyapunov theorem. Russian Mathematical Surveys, 65, 201–202.
  * Weatherburn (1957) Weatherburn, C. E. (1957). An Introduction to Riemannian Geometry and the Tensor Calculus. Cambridge University Press, Cambridge, UK.