# Efficient Inference on High-Dimensional Linear Models with Missing Outcomes

Yikun Zhang Email<EMAIL_ADDRESS>Alexander Giessing Email<EMAIL_ADDRESS>and Yen-Chi Chen Email<EMAIL_ADDRESS>
( Department of Statistics,
University of Washington

)

###### Abstract

This paper is concerned with inference on the conditional mean of a high-
dimensional linear model when outcomes are missing at random. We propose an
estimator which combines a Lasso pilot estimate of the regression function
with a bias correction term based on the weighted residuals of the Lasso
regression. The weights depend on estimates of the missingness probabilities
(propensity scores) and solve a convex optimization program that trades off
bias and variance optimally. Provided that the propensity scores can be
consistently estimated, the proposed estimator is asymptotically normal and
semi-parametrically efficient among all asymptotically linear estimators. The
rate at which the propensity scores are consistent is essentially irrelevant,
allowing us to estimate them via modern machine learning techniques. We
validate the finite-sample performance of the proposed estimator through
comparative simulation studies and the real-world problem of inferring the
stellar masses of galaxies in the Sloan Digital Sky Survey.

Keywords: High-dimensional inference; Missing data; Semi-parametric
efficiency; Lasso.

##  1 Introduction

In this paper, we develop a novel method to conduct statistical inference on
the conditional mean (or regression function) of a sparse high-dimensional
linear model when outcomes are missing at random. Specifically, let
$Y\in\mathbb{R}$ be an outcome variable and $X\in\mathbb{R}^{d}$ be a high-
dimensional covariate vector. The object of interest is the regression
function $m_{0}(x)=\mathrm{E}(Y|X=x)$.

Valid inference on $m_{0}(x)$ in the presence of high-dimensional covariates
and missing outcomes is of practical importance. For one thing, collecting
high-dimensional data has become common practice across the board in science
and engineering. Examples range from portfolio optimization and cross-
sectional home price analysis in finance (Fan et al., 2011) to the development
of biomarker classifiers in biology (Baek et al., 2009). Moreover, some recent
works in causal inference lean toward generating high-dimensional covariates
in order to control for potential confounders (Wyss et al., 2022). For
another, semi-automatic processing and storing of vast amounts of unstructured
data inevitably entails missingness (Huang and Knowles, 2016). Such
missingness may result from the study dropouts of participants in clinical
trial (Higgins et al., 2008) or noncompliance with the assigned treatments in
a survey (Frumento et al., 2012). The emerging field of semi-supervised
learning in computer science, where additional data samples with the same
distribution are given but without labels, is also a missing-outcome problem
(Chapelle et al., 2006).

Figure 1: Comparison of our debiased estimators under two different choices of
the tuning parameters (“1SE” and “min-feas”) with the conventional Lasso
estimates based on complete-case or oracle data. Here, we adopt the sparse
$\beta_{0}^{sp}$ and dense $x^{(4)}$ with $d=1000,n=900$,
$X\sim\mathcal{N}_{d}\left(0,\Sigma^{\mathrm{cs}}\right)$, and
$\epsilon\sim\mathcal{N}(0,1)$ under the MAR setting (22); see Section 4.1 for
details.

When the outcome variables $Y_{i},i=1,...,n$ are fully observed given a data
sample (which is known as the oracle data setting), statistical estimation on
$m_{0}(x)$ is tractable with regularization and sparsity constraints under
high-dimensional settings (Wainwright, 2019). One of the most well-studied
approaches is Lasso (Tibshirani, 1996), which assumes the linear model and
imposes $L_{1}$-regularization on the regression coefficients. When it comes
to statistical inference, however, the Lasso solution leads to a biased
estimate of $m_{0}(x)$ even when the linear model assumption is correct (van
de Geer et al., 2014; Zhang and Zhang, 2014). Missing outcomes further
exacerbate the bias of the Lasso solution; see the first panel of Figure 1 for
an illustration. While this bias could be partially mitigated via sample-
splitting or re-fitting, doing so would reduce the sample size and increase
the computational cost, leading to less efficient estimates. Consequently, the
goal of this paper is to address one central question:

_“How can we conduct statistically and computationally efficient inference on
$m_{0}(x)$ despite missing outcomes?”_

This is a challenging question; to answer it, we need to impose some
structural assumptions. Throughout the paper, we will assume that the data
sample follows a high-dimensional sparse linear regression model with outcomes
“missing at random” (MAR).

######  Assumption 1 (High-dimensional sparse linear regression model with MAR
outcomes).

  1. (a)

The data sample $\\{(Y_{i},R_{i},X_{i})\\}_{i=1}^{n}$ consists of independent
and identically distributed (i.i.d.) observations drawn from
$(Y,R,X)\in\mathbb{R}\times\\{0,1\\}\times\mathbb{R}^{d}$ generated by the
linear model

$Y=X^{T}\beta_{0}+\epsilon,\quad\quad\mathrm{E}\left(\epsilon\big{|}X\right)=0,\quad\quad\mathrm{E}\left(\epsilon^{2}\big{|}X\right)=\sigma_{\epsilon}^{2},$
(1)

where
$\left|\left|\beta_{0}\right|\right|_{0}=\sum_{k=1}^{d}\mathbbm{1}_{\\{\beta_{0k}\neq
0\\}}=s_{\beta}\ll d$.

  2. (b)

The missingness indicator $R$ is conditionally independent of $Y$ given $X$.

Assumption 1(a) is standard in high-dimensional statistics, and the sparse
linear model can be extended to sparse additive models (Ravikumar et al.,
2009) and partially linear models (Müller and van de Geer, 2015; Belloni et
al., 2019); however, such generalizations are beyond the scope of this paper.
Assumption 1(b) is common in the missing data literature (Tsiatis, 2007;
Little and Rubin, 2019) and is related to the ignorability or unconfoundedness
condition in causal inference (Imbens, 2004). Under the MAR assumption, the
propensity score $\pi(Y,X):=\mathrm{P}(R=1|Y,X)$ depends only on the covariate
vector $X$ so that $\pi(Y,X)\equiv\pi(X):=\mathrm{P}(R=1|X)$, which can thus
be estimated from the fully observed data $\\{(X_{i},R_{i})\\}_{i=1}^{n}$
(Rosenbaum and Rubin, 1983).

Under Assumption 1, our proposed debiasing method directly infers the
regression function $m_{0}(x)=x^{T}\beta_{0}$ by combining a Lasso pilot
estimate of $\beta_{0}$ based on the complete-case data with a bias correction
term based on the weighted residuals of the Lasso regression. The weights
depend on estimates of the propensity scores $\pi(X_{i}),i=1,...,n$ and are
obtained as the solution to a convex debiasing program which trades off bias
and variance in a mean-squared-error-optimal way; see Figure 1(b,c) for a
preview.

A crucial difference between the existing works of debiasing the Lasso
estimate in the literature (Zhang and Zhang, 2014; van de Geer et al., 2014;
Javanmard and Montanari, 2014a, b) and our proposal is that we focus our
inference on the scalar regression function $m_{0}(x)=x^{T}\beta_{0}$ but not
the high-dimensional regression (coefficient) vector
$\beta_{0}\in\mathbb{R}^{d}$. By considering $m_{0}(x)$, we are able to
conduct inference on any individual regression coefficient by setting $x$
equal to a standard basis vector in $\mathbb{R}^{d}$ as well as on joint
effects of several regression coefficients. Additionally, inferring $m_{0}(x)$
itself is more computationally efficient than inferring the regression vector
$\beta_{0}$.

###  1.1 Contribution and Outline of the Paper

1\. Methodology: We describe the detailed procedures of our proposed debiasing
method and reveal its computational feasibility through the dual formulation.
We also provides key interpretations to our debiasing method from the
perspective of bias-variance trade-off and Neyman near-orthogonality; see
Section 2 for details.

2\. Asymptotic Theory: We prove that our proposed debiased estimator is
asymptotically unbiased and normally distributed as long as the propensity
scores are consistently estimated in a mild rate of convergence; see Section
3.4. More essentially, the asymptotic variance of our debiased estimator
attains the semi-parametric efficiency bound among all asymptotically linear
estimators. To establish this asymptotic normality, we derive the consistency
of Lasso pilot estimate with complete-case data and the dual solution of our
debiasing program as building blocks; see Section 3.3.

3\. Simulations and Real-World Applications: To demonstrate the finite-sample
performance of our debiasing method, we compare it with the existing debiasing
methods in the literature through Monte Carlo experiments under comprehensive
simulation settings; see Section 4.1 and Section 4.2. We also apply the
proposed debiasing method to the problem of inferring the stellar masses of
galaxies in the Sloan Digital Sky Survey; see Section 5.

4\. Algorithmic Implementation: We describe the implementation details of our
debiasing method in Section A of the supplement and encapsulate them into both
a Python package “Debias-Infer” and R package “DebiasInfer”. All other codes
for our experiments are available at https://github.com/zhangyk8/Debias-
Infer/tree/main/Paper_Code.

###  1.2 Related Work

Statistical inference for the high-dimensional linear model with complete data
has been studied by Zhang and Zhang (2014); van de Geer et al. (2014);
Javanmard and Montanari (2014a, b, 2018); Bühlmann (2013); Ning and Liu
(2017), where they proposed confidence intervals and multiple testing
adjustments for individual regression coefficients or their projections to
low-dimensional components. The main idea of their approaches is to correct
the bias of Lasso or ridge regression vector (Hoerl and Kennard, 1970) by
estimating a relaxed inverse of the gram matrix or the projection matrix,
which is also known as the one-step update (van der Vaart, 2000) in the
classical study of semi-parametric inference. More recently, Li (2020)
provided a bootstrap procedure for debiasing the Lasso solution, while
Javanmard and Lee (2020) presented a flexible framework for testing general
hypotheses of regression coefficients. Battey and Reid (2023) considered
inferring the linear model by treating each regression coefficient in turn as
the interest parameter and finding an optimal transformation to orthogonalize
the remaining coefficients. Moreover, some of the above ideas have been
extended to parameter estimation and inference for high-dimensional
generalized linear models by Belloni et al. (2016); Guo and Chen (2016); Xia
et al. (2020); Shi et al. (2021); Cai et al. (2021); Guo et al. (2021); Ma et
al. (2021).

The asymptotic normality of the aforementioned debiasing methods requires the
regression vector $\beta_{0}$ to lie in the ultra-sparse regime
$s_{\beta}=o\left(\frac{\sqrt{n}}{\log d}\right)$ as our proposed debiasing
method. Under the known covariance matrix and Gaussian design, Javanmard and
Montanari (2018) alleviated the sparsity requirement to
$s_{\beta}=o\left(\frac{n}{(\log d)^{2}}\right)$. In addition, Cai and Guo
(2017) considered constructing confidence intervals for the regression
function $m_{0}(x)=x^{T}\beta_{0}$ in the moderate-sparse region
$\frac{\sqrt{n}}{\log d}\ll s_{\beta}\lesssim\frac{n}{\log d}$ and establish
the minimaxity and adaptivity for their confidence intervals with a prior
knowledge of the sparsity level $s_{\beta}$; see also Nickl and van de Geer
(2013); Cai and Guo (2018) for the related discussion.

Statistical estimation of the regression function under the MAR mechanism and
its semi-parametric efficiency have also been investigated in the literature
(Robins et al., 1994; Robins and Rotnitzky, 1995; Robins et al., 1995; Graham,
2011; Müller and Keilegom, 2012). Although their studies addressed the
problems from missing outcomes to unobserved covariates, they mainly focused
on low-dimensional data. In the high-dimensional data setting, Loh and
Wainwright (2012) proposed an estimator of the linear regression coefficient
based on projected gradient descent in the presence of MAR covariates. High-
dimensional estimation and inference with missing covariates were also studied
by Wang et al. (2019), but the covariates are assumed to be missing complete
at random (MCAR). Recently, Celentano and Wainwright (2023) considered semi-
parametric estimation of the population mean under MAR outcomes in the high-
dimensional inconsistency regime and proposed their debiasing remedy. The most
closely related work to our paper is Chakrabortty et al. (2019), which also
studied the high-dimensional estimation and inference problems with missing
outcomes. They proposed a general M-estimation framework to conduct
statistical inference via a Lasso-type debiased and doubly robust estimator.
Computationally, their debiased estimator requires the estimation of a
$d\times d$ debiasing matrix, while our debiasing method is more efficient
because we only solve for a $d$-dimensional debiasing vector through a convex
program. Theoretically, both Chakrabortty et al. (2019) and us impose the same
condition on the rate of convergence of the propensity score estimation in
pursuit of the asymptotic normality. Nevertheless, we demonstrate through
simulation studies that our debiased estimator is still asymptotically normal
even when the propensity scores are estimated by nonparametric methods.
Furthermore, the debiased estimator in Chakrabortty et al. (2019)
asymptotically achieves the coordinatewise semi-parametric efficiency bound,
which may not be efficient in the worst possible query direction (Jankova and
van de Geer, 2018), while our debiased estimator is asymptotically efficient
whatever the query direction is.

###  1.3 Notation

The general probability measure and expectation with respect to the
distribution of $(Y,R,X)$ are denoted by $\mathrm{P}$ and $\mathrm{E}$,
respectively. We write $Y{\perp\\!\\!\\!\perp}R$ when the random variables $Y$
and $R$ are independent. We also use the big-$O$ notation $h(x)=O(g(x))$ if
the absolute value of $h(x)$ is upper bounded by a positive constant multiple
of $g(x)$ for all sufficiently large $x$. In contrast, $h(x)=o(g(x))$ when
$\lim_{x\to\infty}|h(x)|/g(x)=0$. For random vectors, the notation $o_{P}(1)$
is short for a sequence of random vectors that converges to zero in
probability. The expression $O_{P}(1)$ denotes a sequence that is bounded in
probability. The norm
$\left|\left|x\right|\right|_{q}=\left(\sum_{k=1}^{d}x_{k}^{q}\right)^{1/q}$
with $q>0$ stands for the $L_{q}$-norm in the Euclidean space
$\mathbb{R}^{d}$, though it is no longer a norm when $0<q<1$. In particular,
$\left|\left|x\right|\right|_{\infty}=\max_{1\leq k\leq d}|x_{k}|$ and
$\left|\left|x\right|\right|_{0}=\sum_{k=1}^{n}\mathbbm{1}_{\\{x_{k}\neq
0\\}}$ indicates the number of nonzero elements in $x\in\mathbb{R}^{d}$.
Furthermore,
$\left|\left|X\right|\right|_{q}=\left(\mathrm{E}|Z|^{q}\right)^{1/q}$ with
$q\geq 1$ is the $L_{q}$-norm for a random variable $Z$. We use the notation
$a_{n}\lesssim b_{n}$ or $b_{n}\gtrsim a_{n}$ when there exists an absolute
constant $A>0$ such that $a_{n}\leq Cb_{n}$ when $n$ is large. If
$a_{n}\gtrsim b_{n}$ and $a_{n}\lesssim b_{n}$, then $a_{n},b_{n}$ are
asymptotically equal and it is denoted by $a_{n}\asymp b_{n}$. Finally, we
denote the unit sphere in $\mathbb{R}^{d}$ by
$\mathbb{S}^{d-1}=\\{x\in\mathbb{R}^{d}:\left|\left|x\right|\right|_{2}=1\\}$
and a ball centered at $x$ with radius $r$ in $\mathbb{R}^{d}$ by
$B_{d}(x,r)=\\{y\in\mathbb{R}^{d}:\left|\left|y-x\right|\right|_{2}\leq r\\}$.

##  2 Methodology

In this section, we outline the detailed procedure of our debiasing inference
method and discuss a feasible solution to the key debiasing program through
its dual form. We further motivate our debiasing method from the perspectives
of bias-variance trade-off and Neyman near-orthogonality.

###  2.1 Debiasing Inference Procedure

Recall our i.i.d. data
$\\{(Y_{i},R_{i},X_{i})\\}_{i=1}^{n}\subset\mathbb{R}\times\\{0,1\\}\times\mathbb{R}^{d}$,
where $Y_{i}$ is the outcome variable with missing indicator $R_{i}$ and
$X_{i}$ is the fully observed high-dimensional covariate vector for
$i=1,...,n$. In order to conduct statistical inference on the regression
function $m_{0}(x)=\mathrm{E}(Y|X=x)=x^{T}\beta_{0}$, we propose a debiasing
method with the following procedure.

  * •

Step 1: Compute the Lasso pilot estimate $\widehat{\beta}$ with complete-case
data

$\widehat{\beta}=\operatorname*{arg\,min}\limits_{\beta\in\mathbb{R}^{d}}\left[\frac{1}{n}\sum\limits_{i=1}^{n}R_{i}(Y_{i}-X_{i}^{T}\beta)^{2}+\lambda\left|\left|\beta\right|\right|_{1}\right],$
(2)

where $\lambda>0$ is a regularization parameter.

  * •

Step 2: Obtain consistent estimates $\widehat{\pi}_{i},i=1,...,n$ of the
propensity scores $\pi_{i}=\pi(X_{i})=\mathrm{P}(R_{i}=1|X_{i}),i=1,...,n$ by
any machine learning method (not necessarily a parametric model) on the data
$\\{(X_{i},R_{i})\\}_{i=1}^{n}\subset\mathbb{R}^{d}\times\\{0,1\\}$.

  * •

Step 3: Solve for the debiasing weight vector
$\widehat{\bm{w}}\equiv\widehat{\bm{w}}(x)=\left(\widehat{w}_{1}(x),...,\widehat{w}_{n}(x)\right)^{T}\in\mathbb{R}^{n}$
through a debiasing program defined as:

$\min_{\bm{w}\in\mathbb{R}^{n}}\left\\{\sum_{i=1}^{n}\widehat{\pi}_{i}w_{i}^{2}:\left|\left|x-\frac{1}{\sqrt{n}}\sum_{i=1}^{n}w_{i}\cdot\widehat{\pi}_{i}\cdot
X_{i}\right|\right|_{\infty}\leq\frac{\gamma}{n}\right\\},$ (3)

where $\gamma>0$ is a tuning parameter.

  * •

Step 4: Define the debiased estimator for $m_{0}(x)$ as:

$\widehat{m}^{\text{debias}}(x;\widehat{\bm{w}})=x^{T}\widehat{\beta}+\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\widehat{w}_{i}(x)R_{i}\left(Y_{i}-X_{i}^{T}\widehat{\beta}\right).$
(4)

  * •

Step 5: Construct the asymptotic $(1-\tau)$-level confidence interval for
$m_{0}(x)$ as:

$\left[\widehat{m}^{\text{debias}}(x;\widehat{\bm{w}})\pm\Phi^{-1}\left(1-\frac{\tau}{2}\right)\cdot\sigma_{\epsilon}\cdot\sqrt{\frac{1}{n}\sum_{i=1}^{n}\widehat{\pi}_{i}\widehat{w}_{i}(x)^{2}}\right],$
(5)

where $\Phi(\cdot)$ denotes the cumulative distribution function (CDF) of
$\mathcal{N}(0,1)$. If $\sigma_{\epsilon}^{2}$ is unknown, then it can be
replaced by any consistent estimator $\widehat{\sigma}_{\epsilon}^{2}$.

In practice, we fit the Lasso pilot estimate $\widehat{\beta}$ in Step 1 by
the scaled Lasso (Sun and Zhang, 2012) so as to automatically select the
regularization parameter $\lambda>0$ and simultaneously produce a consistent
estimator of the noise level $\sigma_{\epsilon}^{2}$. As for the propensity
score estimation in Step 2, the theoretical validity of our proposed method
only requires the estimated propensity scores $\widehat{\pi}_{i},i=1,...,n$ to
be consistent at a mild rate of convergence and thus allows the usage of any
modern machine learning approach; see Assumption 5 with more details in
Section 3.3 and Section 3.4. The formulation of our debiasing program (3) in
Step 3 can be motivated from the perspective of bias-variance trade-off in
Section 2.3, and the subsequent definition of the debiased estimator (4) in
Step 4 is inspired by the original debiased Lasso estimator in Zhang and Zhang
(2014); van de Geer et al. (2014); Javanmard and Montanari (2014a). However,
instead of estimating the (pseudo-)inverse matrix of the complete-case
population gram matrix $\left[\mathrm{E}\left(RXX^{T}\right)\right]^{-1}$ as
the debiased Lasso estimator, we optimize the debiasing weight vector through
our debiasing program (3), which is computationally more efficient. Section
2.2 sheds light on a feasible avenue to solve (3) from its dual formulation
that becomes an unconstrained convex program. Finally, we prove in Section 3.4
that the debiased estimator $\widehat{m}^{\text{debias}}(x;\widehat{\bm{w}})$
is asymptotically normal and that its asymptotic variance can be estimated
consistently via $\sum_{i=1}^{n}\widehat{\pi}_{i}\widehat{w}_{i}^{2}(x)$. This
guarantees the asymptotic validity of the confidence interval in Step 5.

###  2.2 Dual Formulation of the Debiasing Program (3)

The proposed debiasing program (3) is a quadratic programming problem with box
constraints. By standard duality theory in optimization, we obtain the
following result.

######  Proposition 1.

The dual form of the debiasing program (3) is given by

$\min_{\ell\in\mathbb{R}^{d}}\left\\{\frac{1}{4n}\sum_{i=1}^{n}\widehat{\pi}_{i}\left[X_{i}^{T}\ell\right]^{2}+x^{T}\ell+\frac{\gamma}{n}\left|\left|\ell\right|\right|_{1}\right\\}.$
(6)

If strong duality holds, then the solution
$\widehat{\bm{w}}(x)\in\mathbb{R}^{n}$ to the primal debiasing program (3) and
the solution $\widehat{\ell}(x)\in\mathbb{R}^{d}$ to the dual debiasing
program (6) satisfy the following identity:

$\widehat{w}_{i}(x)=-\frac{1}{2\sqrt{n}}\cdot X_{i}^{T}\widehat{\ell}(x),\quad
i=1,...,n.$ (7)

The proof of Proposition 1 is in Section E.2. This result has two important
consequences: First, it is difficult to select the tuning parameter $\gamma>0$
based on the constrained quadratic program (3) alone: if $\gamma>0$ is too
small, the the program is infeasible; if $\gamma>0$ is too large, the estimate
has a non-negligible bias (see Section 2.3). However, since the dual program
(6) is an unconstrained quadratic programming problem, we can select the
tuning parameter $\gamma>0$ via cross-validation on the dual program. Second,
if strong duality holds, the debiased estimator (4) can be expressed via (7)
in terms of the dual solution $\widehat{\ell}(x)$ as:

$\widehat{m}^{\text{debias}}(x;\widehat{\bm{w}})=x^{T}\widehat{\beta}-\frac{1}{2n}\sum_{i=1}^{n}R_{i}X_{i}^{T}\widehat{\ell}(x)\left(Y_{i}-X_{i}^{T}\widehat{\beta}\right),$
(8)

where the second bias-corrected term reduces to a linear combination of the
sample average of $n$ random variables. This representation is the key to
deriving the asymptotic normality of our debiased estimator; see Section 3.1
and Section 3.4 for details. Notice that strong duality holds whenever
$\gamma>0$ is sufficiently large; see Lemma 3 in Section 3.3 for more details.

###  2.3 Bias-Variance Trade-off Perspective

The design of our debiasing program (3) is motivated by controlling the
conditional mean squared error

$\mathrm{E}\left[\left(\sqrt{n}\cdot\widehat{m}^{\text{debias}}(x;\widehat{\bm{w}})-\sqrt{n}\cdot
m_{0}(x)\right)^{2}\big{|}\bm{X}\right]$

of our debiased estimator (4) given $\bm{X}=(X_{1},...,X_{n})$ and can thus be
interpreted from the perspective of a (conditional) bias-variance trade-off.
To explicate this motivation, we consider the generic debiased estimator
$m^{\text{debias}}(x;\bm{w})$ from (4):

$m^{\text{debias}}(x;\bm{w})=x^{T}\beta+\frac{1}{\sqrt{n}}\sum_{i=1}^{n}w_{i}R_{i}\left(Y_{i}-X_{i}^{T}\beta\right),$
(9)

where one would plug in the Lasso pilot estimate $\widehat{\beta}$ from (2)
and the solution $\widehat{\bm{w}}\equiv\widehat{\bm{w}}(x)\in\mathbb{R}^{n}$
from our debiasing program (3) to obtain a concrete debiased estimator (4) in
practice. We decompose the conditional mean squared error of (9) into three
explicit terms as follows.

######  Proposition 2.

Under Assumption 1, the conditional mean squared error of
$\sqrt{n}\,m^{\mathrm{debias}}(x;\bm{w})$ given $\bm{X}=(X_{1},...,X_{n})$ has
the following decomposition as:

$\displaystyle\begin{split}&\mathrm{E}\left[\left(\sqrt{n}\,m^{\mathrm{debias}}(x;\bm{w})-\sqrt{n}\,m_{0}(x)\right)^{2}\Big{|}\bm{X}\right]\\\
&=\underbrace{\sigma_{\epsilon}^{2}\sum_{i=1}^{n}w_{i}^{2}\pi(X_{i})}_{\text{Conditional
variance
I}}+\underbrace{\left(\beta_{0}-\beta\right)^{T}\left[\sum_{i=1}^{n}w_{i}^{2}\pi(X_{i})\left(1-\pi(X_{i})\right)X_{i}X_{i}^{T}\right]\left(\beta_{0}-\beta\right)}_{\text{Conditional
variance II}}\\\
&\quad+\underbrace{\left[\left(\frac{1}{\sqrt{n}}\sum_{i=1}^{n}w_{i}\pi(X_{i})X_{i}-x\right)^{T}\sqrt{n}\left(\beta_{0}-\beta\right)\right]^{2}}_{\text{Conditional
bias}}.\end{split}$ (10)

The derivation for Proposition 2 can be found in Section E.2. By Hölder’s
inequality, we can upper bound the “Conditional bias” term as:

$\displaystyle\left[\left(\frac{1}{\sqrt{n}}\sum_{i=1}^{n}w_{i}\pi(X_{i})X_{i}-x\right)^{T}\sqrt{n}\left(\beta_{0}-\beta\right)\right]^{2}\leq\left[\left|\left|\frac{1}{\sqrt{n}}\sum_{i=1}^{n}w_{i}\pi(X_{i})X_{i}-x\right|\right|_{\infty}\sqrt{n}\left|\left|\beta_{0}-\beta\right|\right|_{1}\right]^{2}.$

The debiasing program (3) is hence designed to optimize the weights
$w_{i},i=1,...,n$ so that the estimated “Conditional variance I” term
$\sum_{i=1}^{n}\widehat{\pi}_{i}w_{i}^{2}$ is minimized as the objective
function while an upper bound of the estimated “Conditional bias” term
$\left|\left|x-\frac{1}{\sqrt{n}}\sum_{i=1}^{n}w_{i}\widehat{\pi}_{i}X_{i}\right|\right|_{\infty}$
is well-controlled by the constraint. We ignore the “Conditional variance II”
term because, by Lemma 4, it is asymptotically negligible and dominated by the
“Conditional variance I” term if $\beta=\widehat{\beta}$ (solution Lasso
program (2)) and $\bm{w}=\widehat{\bm{w}}(x)$ (solution to the debiasing
program (3)). Similar bias-variance trade-offs also appear in Cai et al.
(2021) for generalized linear models and Giessing and Wang (2021) for quantile
regression models. However, the motivation for the procedures in these two
papers is ad-hoc and not grounded in a rigorous decomposition of the
conditional mean squared error. Upon minimizing the (conditional) variance, we
expect our debiased estimator to be asymptotically more efficient than other
estimators; see our discussions in Section 3.4.

###  2.4 Neyman Near-Orthogonalization Viewpoint

Our debiasing method in Section 2.1 can also be interpreted from the
perspective of Neyman near-orthogonalization (Neyman, 1959, 1979; Chernozhukov
et al., 2018). To this end, we regard the regression function $m\equiv
m(x)\in\mathbb{R}$ as the main parameter that we want to infer and the
regression vector $\beta\in\mathbb{R}^{d}$ as the high-dimensional nuisance
parameter. Under Assumption 1, the true values of these two parameters are
$m_{0}=m_{0}(x)=x^{T}\beta_{0}$ and $\beta_{0}$, respectively. We define the
generic score function as:

$\Xi_{x}(Y,R,X;m,\beta)=m-x^{T}\beta-\sqrt{n}\cdot w\cdot R(Y-X^{T}\beta),$

where $w\in\mathbb{R}$ is a symbolic weight. Given the observed data
$\\{(Y_{i},R_{i},X_{i})\\}_{i=1}^{n}$ and $\bm{w}\in\mathbb{R}^{n}$,
$\beta\in\mathbb{R}^{d}$ arbitrary, the generic debiased estimator
$m^{\text{debias}}(x,\bm{w})$ solves the sample-based estimating equation

$\frac{1}{n}\sum_{i=1}^{n}\Xi_{x}(Y_{i},R_{i},X_{i};m^{\text{debias}},\beta)=m^{\text{debias}}(x;\bm{w})-x^{T}\beta-\frac{1}{\sqrt{n}}\sum_{i=1}^{n}w_{i}\cdot
R_{i}\left(Y_{i}-X_{i}^{T}\beta\right)=0.$ (11)

Provided that Assumption 1 holds and the weights $\bm{w}\in\mathbb{R}^{n}$
satisfy the box constraint in (3), the Neyman near-orthogonalization condition
given the data $\bm{X}=(X_{1},...,X_{n})^{T}\in\mathbb{R}^{n\times d}$ at
$(m_{0},\beta_{0})$ requires that

$\displaystyle\begin{split}&\mathrm{E}\left[\frac{1}{n}\sum_{i=1}^{n}\Xi_{x}(Y_{i},R_{i},X_{i};m_{0},\beta_{0})\bigg{|}\bm{X}\right]=0,\\\
&\sup_{\beta\in\mathcal{T}_{n}}\left|\left\\{\frac{\partial}{\partial\beta}\mathrm{E}\left[\frac{1}{n}\sum_{i=1}^{n}\Xi_{x}(Y_{i},R_{i},X_{i};m,\beta)\Big{|}\bm{X}\right]\bigg{|}_{(m_{0},\beta_{0})}\right\\}^{T}(\beta-\beta_{0})\right|\leq\frac{\delta_{n}}{\sqrt{n}},\end{split}$
(12)

where $\mathcal{T}_{n}$ is a properly shrinking neighborhood of $\beta_{0}$
and $\delta_{n}=o(1)$; see Definition 2.2 in Chernozhukov et al. (2018).
Indeed, the first condition obviously holds true, while the second condition
is also satisfied because for any $\beta\in\mathcal{T}_{n}$,

$\displaystyle\left|\left\\{\frac{1}{n}\sum_{i=1}^{n}\frac{\partial}{\partial\beta}\mathrm{E}\left[\Xi_{x}(Y_{i},R_{i},X_{i};m,\beta)|X\right]\big{|}_{(m_{0},\beta_{0})}\right\\}^{T}(\beta-\beta_{0})\right|$
$\displaystyle=\left|\left[x-\frac{1}{\sqrt{n}}\sum_{i=1}^{n}w_{i}\cdot\pi(X_{i})X_{i}\right]^{T}(\beta_{0}-\beta)\right|$
$\displaystyle\stackrel{{\scriptstyle\text{(a)}}}{{``\leq"}}\left|\left|x-\frac{1}{\sqrt{n}}\sum_{i=1}^{n}w_{i}\cdot\widehat{\pi}_{i}\cdot
X_{i}\right|\right|_{\infty}\left|\left|\beta-\beta_{0}\right|\right|_{1}$
$\displaystyle\>\>\stackrel{{\scriptstyle\text{(b)}}}{{\leq}}\frac{\gamma}{n}\left|\left|\beta-\beta_{0}\right|\right|_{1}\leq\frac{\delta_{n}}{\sqrt{n}},$

where the approximate inequality (a) follows from Hölder’s inequality and
substituting the estimated propensity scores $\widehat{\pi}_{i},i=1,...,n$
while (b) holds since $\bm{w}\in\mathbb{R}^{n}$ satisfy the box constraint in
(3). Therefore, the nuisance realization set $\mathcal{T}_{n}$ can be chosen
as
$\\{\beta\in\mathcal{B}\subset\mathbb{R}^{d}:\left|\left|\beta-\beta_{0}\right|\right|_{1}\leq\frac{\sqrt{n}\delta_{n}}{\gamma}\\}$
for some convex set $\mathcal{B}$ containing $\beta_{0}$. We will show in
Theorem 1 that our Lasso pilot estimate $\widehat{\beta}$ satisfies
$\left|\left|\widehat{\beta}-\beta_{0}\right|\right|_{1}=O_{P}\left(s_{\beta}\sqrt{\frac{\log
d}{n}}\right)$ and thus, our final debiased estimator (4) satisfies the Neyman
near-orthogonality condition (12) with a fine-tuned parameter $\gamma>0$.

By Theorem 3.1 in Chernozhukov et al. (2018), the (asymptotic) variance of
$m^{\text{debias}}(x,\bm{w})$ is
$\sigma_{\epsilon}^{2}\,\mathrm{E}\left[\sum_{i=1}^{n}w_{i}^{2}R_{i}^{2}|\bm{X}\right]=\sigma_{\epsilon}^{2}\sum_{i=1}^{n}w_{i}^{2}\pi(X_{i})$.
Therefore, our debiasing program (3) is designed to pursue the debiased
estimator with the smallest (estimated) variance among all the estimators that
(approximately) satisfies the constraints in (12). Consequently, it is
expected that our debiased estimator is more efficient than most other
estimators even without explicitly constructing an efficient score function
satisfying the exact Neyman orthogonality; see Section 2.2.2 in Chernozhukov
et al. (2018) with references therein. Finally, in light of (12), our
debiasing method de-correlates the Lasso pilot regression from the propensity
score estimation and the optimization of weights for the debiased estimator so
that we are able to establish the asymptotic normality of our debiased
estimator without resorting to sample-splitting; see Section 3.4 for more
details.

##  3 Asymptotic Theory

In this section, we study the consistency results of our debiasing inference
method in Section 2.1 and its asymptotic normality property. The high-level
motivation and overview of the technical results are given in Section 3.1 with
more detailed studies in the follow-up subsections.

###  3.1 Heuristics and Overview

In order to derive the asymptotic normality of our debiased estimator (4), we
need to analyze the asymptotic behavior of
$\sqrt{n}\left[\widehat{m}^{\text{debias}}(x;\widehat{\bm{w}})-m_{0}(x)\right]$.
Under Assumption 1 and the definition (4) of
$\widehat{m}^{\text{debias}}(x;\widehat{\bm{w}})$, we deduce that

$\displaystyle\sqrt{n}\left[\widehat{m}^{\text{debias}}(x;\widehat{\bm{w}})-m_{0}(x)\right]$
$\displaystyle=\sum_{i=1}^{n}\widehat{w}_{i}(x)R_{i}\epsilon_{i}+\left[x-\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\widehat{w}_{i}(x)R_{i}X_{i}\right]^{T}\sqrt{n}\left(\widehat{\beta}-\beta_{0}\right),$

where $\epsilon_{i},i=1,...,n$ are i.i.d. noise variables under model (1).
From above expression, it is unclear how one can apply the central limit
theorem or related asymptotic theories to derive the asymptotic distribution
of
$\sqrt{n}\left[\widehat{m}^{\text{debias}}(x;\widehat{\bm{w}})-m_{0}(x)\right]$,
because the limiting distribution of the weight vector
$\widehat{\bm{w}}\equiv\widehat{\bm{w}}(x)=\left(\widehat{w}_{1}(x),...,\widehat{w}_{n}(x)\right)^{T}\in\mathbb{R}^{n}$
is intractable and may not even exist. However, the dual formulation in
Proposition 1 provides a promising direction of approaching this problem. By
the dual representation (8) of the debiased estimator, we have

$\displaystyle\sqrt{n}\left[\widehat{m}^{\text{debias}}(x;\widehat{\bm{w}})-m_{0}(x)\right]$
$\displaystyle=-\frac{1}{2\sqrt{n}}\sum_{i=1}^{n}R_{i}\epsilon_{i}X_{i}^{T}\widehat{\ell}(x)+\left[x+\frac{1}{2n}\sum_{i=1}^{n}R_{i}X_{i}^{T}\widehat{\ell}(x)X_{i}\right]^{T}\sqrt{n}\left(\beta_{0}-\widehat{\beta}\right)$
$\displaystyle=-\frac{1}{2\sqrt{n}}\sum_{i=1}^{n}R_{i}\epsilon_{i}X_{i}^{T}\ell_{0}(x)+\text{``Bias
terms''},$ (13)

where $\ell_{0}(x)\in\mathbb{R}^{d}$ is the deterministic solution to the
population dual program for which the dual solution
$\widehat{\ell}(x)\in\mathbb{R}^{d}$ is consistent; see (17) in Section 3.2
for its definition. The leading term in (3.1) is a sample average of $n$
i.i.d. scalar random variables, which can be handled via the standard central
limit theorem (CLT), while the “Bias terms” in (3.1) are of order $o_{P}(1)$
under mild regularity conditions; see Theorem 5.

The rest of this section corroborates the above heuristic arguments. After
introducing the notations and necessary assumptions, we derive the following
main results:

  * •

Consistency of the Lasso pilot estimate (2): Under a suitable choice of
$\lambda>0$, it holds that
$\left|\left|\widehat{\beta}-\beta_{0}\right|\right|_{2}=O_{P}\left(\frac{1}{\kappa_{R}^{2}}\sqrt{\frac{s_{\beta}\log
d}{n}}\right)$ when $\log d=o(n)$ with $d>n$, where $\kappa_{R}^{2}>0$ is the
minimum eigenvalue of the complete-case gram matrix
$\mathrm{E}\left(RXX^{T}\right)$; see Theorem 1 for details.

  * •

Consistency of the solution to the dual program (6): Under a suitable choice
of $\gamma>0$, we deduce that
$\left|\left|\widehat{\ell}(x)-\ell_{0}(x)\right|\right|_{2}=O_{P}\left(\frac{1}{\kappa_{R}^{3}}\sqrt{\frac{s_{\ell}(x)\log
d}{n}}+\frac{r_{\ell}}{\kappa_{R}^{4}}+\frac{r_{\pi}\sqrt{s_{\ell}(x)}}{\kappa_{R}^{4}}\right)$,
where $\ell_{0}(x)$ is the solution to the population dual program,
$r_{\ell}>0$ is a parameter for the sparse approximation to $\ell_{0}(x)$, and
$r_{\pi}>0$ is the rate of convergence of $\max_{1\leq i\leq
n}\left|\widehat{\pi}_{i}-\pi(X_{i})\right|$; see Theorem 2 and its related
discussions.

  * •

Asymptotic normality of the debiased estimator (4): Combining the above
consistency results with other mild regularity conditions, we prove that

$\sqrt{n}\left[\widehat{m}^{\mathrm{debias}}(x;\widehat{\bm{w}})-m_{0}(x)\right]\stackrel{{\scriptstyle
d}}{{\to}}\mathcal{N}\left(0,\,\sigma_{m}^{2}(x)\right)$

and discuss the semi-parametric efficiency of $\sigma_{m}^{2}(x)$; see Theorem
5 for details.

###  3.2 Notations and Assumptions

We define the $\psi_{\alpha}$-Orlicz norm for any random variable $Z$ and
$\alpha\in(0,2]$ by

$\left|\left|Z\right|\right|_{\psi_{\alpha}}=\inf\left\\{t>0:\mathrm{E}\left[\exp\left(\frac{|Z|^{\alpha}}{t^{\alpha}}\right)\right]\leq
2\right\\}.$ (14)

For $\alpha\in(0,1)$, there is an alternative definition that modifies (14)
into a legitimate norm which can be used interchangeably with (14) up to an
absolute constant that only depends on $\alpha$; see page 266 in van der Vaart
and Wellner (1996).

We enumerate the regularity conditions under the high-dimensional sparse
linear modeling Assumption 1 for our subsequent theoretical analysis.

######  Assumption 2 (Sub-Gaussian covariate).

The covariate vector $X\in\mathbb{R}^{d}$ is sub-Gaussian, _i.e._ ,
$\left|\left|u^{T}X\right|\right|_{\psi_{2}}\lesssim\left[\mathrm{E}\left(u^{T}XX^{T}u\right)\right]^{1/2}$
for any $u\in\mathbb{R}^{d}$.

We introduce the standard sub-Gaussian condition as Assumption 2 on the
covariate vector $X$ in order to control the exponential tail behavior of $X$
by its second moment and facilitate our proofs. Notice that the covariate
vector $X\in\mathbb{R}^{d}$ is not necessarily centered in our definition of
sub-Gaussian covariate vector. We denote the $s$-sparse maximum eigenvalue of
the population gram matrix
$\mathrm{E}\left(XX^{T}\right)\in\mathbb{R}^{d\times d}$ by

$\phi_{s}:=\sup_{u\in\mathbb{S}^{d-1},\left|\left|u\right|\right|_{0}\leq
s}\mathrm{E}\left[(X^{T}u)^{2}\right],$ (15)

where $\mathbb{S}^{d-1}=\\{u\in\mathbb{R}^{d}:\|u\|_{2}=1\\}$. In particular,
$\mathrm{E}(X_{ik}^{2})\leq\max\limits_{1\leq i\leq n}\max\limits_{1\leq k\leq
d}\mathrm{E}(X_{ik}^{2})\leq\phi_{1}$ for each covariate vector
$X_{i}=(X_{i1},...,X_{id})^{T}\in\mathbb{R}^{d}$. It is noteworthy that under
the high-dimensional setting with $s\ll d$, the $s$-sparse maximum eigenvalue
$\phi_{s}$ is substantially smaller than the unrestricted maximum eigenvalue
$\phi_{d}$ of $\mathrm{E}\left(XX^{T}\right)$ and can even be independent of
the dimension $d$.

######  Assumption 3 (Sub-Gaussian noise).

The noise variable $\epsilon\in\mathbb{R}$ is sub-Gaussian, _i.e._ ,
$\left|\left|\epsilon\right|\right|_{\psi_{2}}\lesssim\left[\mathrm{E}(\epsilon^{2})\right]^{1/2}=\sigma_{\epsilon}$
with $\sigma_{\epsilon}^{2}>0$ being the noise level.

Assumption 3 is restrictive and simplifies theoretical arguments. In
principle, it can be relaxed to a moment conditions; see Belloni and
Chernozhukov (2013). We conduct numerical experiments on our debiasing method
in which the sub-Gaussian noise assumption is violated in Section 4.3 and
demonstrate that our proposed method is robust to various heavy-tailed noise
distributions. To establish the consistency of our Lasso pilot estimate (2)
and the dual solution $\widehat{\ell}(x)\in\mathbb{R}^{d}$ to (6), we impose
the following eigenvalue lower bound on the complete-case population gram
matrix $\mathrm{E}\left(RXX^{T}\right)\in\mathbb{R}^{d\times d}$.

######  Assumption 4 (Lower eigenvalue bound on
$\mathrm{E}\left(RXX^{T}\right)$).

The minimum eigenvalue of the complete-case gram matrix
$\mathrm{E}\left(RXX^{T}\right)\in\mathbb{R}^{d\times d}$ is bounded away from
zero, _i.e._ , there exists a constant $\kappa_{R}>0$ such that

$\inf_{v\in\mathbb{S}^{d-1}}\mathrm{E}\left[R(X^{T}v)^{2}\right]\geq\kappa_{R}^{2}.$

Assumption 4 ensures that the regression vector $\beta_{0}\in\mathbb{R}^{d}$
is identifiable. We explicitly introduce the “constant” $\kappa_{R}>0$ to
handle cases in which the minimum eigenvalue of
$\mathrm{E}\left(RXX^{T}\right)$ vanishes (slowly) as $d,n\to\infty$. When
proving the consistency of the Lasso pilot estimate (2), however, we allow the
common relaxation of Assumption 4 that the minimum eigenvalue of
$\mathrm{E}\left(RXX^{T}\right)$ is only bounded away from zero over the cone
$\mathcal{C}_{1}(S_{\beta},3)$:

$\inf_{v\in\mathcal{C}_{1}(S_{\beta},3)\cap\mathbb{S}^{d-1}}\mathrm{E}\left[R(X^{T}v)^{2}\right]\geq\kappa_{R}^{2},$
(16)

where
$\mathcal{C}_{q}(S,\vartheta)=\left\\{v\in\mathbb{R}^{d}:\left|\left|v_{S^{c}}\right|\right|_{q}\leq\vartheta\left|\left|v_{S}\right|\right|_{q}\right\\}$
is a cone of dominant coordinates (Section 7.2 in Koltchinskii 2011) for some
index set $S\subset\\{1,...,d\\}$ and the parameter $\vartheta>0$ with respect
to the norm $\left|\left|\cdot\right|\right|_{q}$ for $q\geq 1$. Here,
$v_{S}\in\mathbb{R}^{|S|}$ denotes the subvector with elements of $v$ indexed
by $S$ and $v_{S^{c}}$ is defined in an analogous manner. This relaxed
assumption (16) is known as the restricted eigenvalue condition in the
literature; see Section 7.3.1 in Wainwright 2019 and Bickel et al. (2009); van
de Geer and Bühlmann (2009).

While it is common in the literature to impose the restricted eigenvalue
conditions on the sample-based complete-case gram matrix
$\frac{1}{n}\sum_{i=1}^{n}R_{i}X_{i}X_{i}^{T}$, Lemma 7 in Section E.9 shows
that under Assumption 2, the population and sample-based gram matrices are
close with high probability, so (16) holds interchangeably on these two gram
matrices in the setting of this paper.

As for the propensity score estimation, we introduce a mild assumption on the
consistency of the estimated propensity scores $\widehat{\pi}_{i}$ to the true
propensity scores $\pi_{i}=\pi(X_{i})$ for $i=1,...,n$.

######  Assumption 5 (Stochastic control of the propensity score estimation).

Given any $n\geq 1$ and $\delta\in(0,1)$, there exists $r_{\pi}\equiv
r_{\pi}(n,\delta)>0$ such that

$\mathrm{P}\left(\max\limits_{1\leq i\leq
n}\left|\widehat{\pi}_{i}-\pi_{i}\right|>r_{\pi}\right)<\delta.$

Assumption 5 is a non-asymptotic version of the stochastic boundedness for the
estimation errors of the propensity scores evaluated at $X_{i},i=1,...n$. In
applications, the exact rate of convergence $r_{\pi}$ depends on the
estimator, which can be parametric, nonparametric, or semiparametric. We
demonstrate in Theorem 2 and Theorem 5 below that as long as
$r_{\pi}=r_{\pi}(n,\delta)\to 0$ for some $\delta\equiv\delta_{n}\to 0$ as
$n\to\infty$, the solution to the dual program (6) is consistent and
consequently, the debiased estimator (4) is asymptotically normal; see Section
3.3 and Section 3.4 for details. The rate at which $r_{\pi}\to 0$ can be very
slow as in (20). These results unveil the possibility of combining our
debiased estimator with any consistent estimates of the propensity scores; see
Section 4.4 for the numerical experiments.

######  Remark 1.

In contrast to the existing literature on regression models with missing
outcomes, incomplete covariates, or potential outcomes (Rosenbaum and Rubin,
1983; Robins et al., 1994; Chakrabortty et al., 2019), we do not require the
positivity condition $\pi(x)=\mathrm{P}(R=1|X=x)>\pi_{\min}>0$ for all
$x\in\mathrm{support}(X)\subset\mathbb{R}^{d}$. As noted by D’Amour et al.
(2021), the positivity condition is too stringent under the context of high-
dimensional covariates. Instead, we merely assume that the complete-case
population gram matrix $\mathrm{E}\left(RXX^{T}\right)$ is positive definite
as in Assumption 4, which is valid even when the positivity condition is
violated. For example, $\mathrm{E}\left(RXX^{T}\right)$ will be positive
definite if $\mathrm{support}(X)$ has a nonzero Lebesgue measure in
$\mathbb{R}^{d}$ and $\pi(x)$ is positive within a subset of
$\mathrm{support}(X)$ with a nonzero Lebesgue measure.

We now state the definition of the population dual solution for which the
solution of our dual formulation (6) is consistent. The population dual
program for any $x\in\mathbb{R}^{d}$ is defined by taking $n\to\infty$ in (6)
under the true propensity scores as:

$\displaystyle\begin{split}\min_{\ell\in\mathbb{R}^{d}}\left\\{\frac{1}{4}\,\mathrm{E}\left[R\left(X^{T}\ell\right)^{2}\right]+x^{T}\ell\right\\}.\end{split}$
(17)

The exact solution to the above population dual program (17) is given by

$\ell_{0}(x):=-2\left[\mathrm{E}\left(RXX^{T}\right)\right]^{-1}x\in\mathbb{R}^{d}$
(18)

accordingly, provided that Assumption 4 holds. Since the population dual
solution $\ell_{0}(x)$ is not necessarily sparse, we consider the following
definition of a $r_{\ell}$-approximation $\widetilde{\ell}(x)$ in order to
establish the consistency of $\widehat{\ell}(x)$; see also Assumption 6.

######  Definition 1 ($r_{\ell}$-approximation to $\ell_{0}(x)$).

For any $x\in\mathbb{R}^{d}$, we define a $r_{\ell}$-approximation to the
population dual solution $\ell_{0}(x)$ as:

$\widetilde{\ell}(x)\equiv\widetilde{\ell}(x;r_{\ell}):=\operatorname*{arg\,min}_{\ell(x)\in\mathbb{R}^{d}}\left\\{\left|\left|\ell(x)\right|\right|_{0}:\left|\left|\ell(x)-\ell_{0}(x)\right|\right|_{2}\leq
r_{\ell}\left|\left|\ell_{0}(x)\right|\right|_{2}\right\\},$ (19)

where $r_{\ell}\in\left[0,\frac{1}{2}\right]$ controls the accuracy of
approximation.

In general, a larger value of $r_{\ell}$ gives rise to a sparser approximation
$\widetilde{\ell}(x)$ to the vector $\ell_{0}(x)$ of interest. On the other
hand, the consistency of $\widehat{\ell}(x)$ requires the approximation
parameter $r_{\ell}$ for $\widetilde{\ell}(x)$ to converge to 0 in a certain
rate as $n\to\infty$; see the discussion after Theorem 2 in Section 3.3 below.
In other words, as the sample size $n$ increases, the exact population dual
solution $\ell_{0}(x)$ must tend to a sparse vector in order for the
consistency result to hold. We furnish several sufficient conditions and
motivating examples under which Assumption 6 below holds in Section B.

######  Assumption 6 (Sparsity of $\widetilde{\ell}(x)$).

The $r_{\ell}$-approximation $\widetilde{\ell}(x)\in\mathbb{R}^{d}$ to the
exact population dual solution to $\ell_{0}(x)\in\mathbb{R}^{d}$ exists and
satisfies

$S_{\ell}(x):=\mathrm{support}\left(\widetilde{\ell}(x)\right)=\left\\{1\leq
k\leq d:\widetilde{\ell}_{k}(x)\neq 0\right\\}\;\text{ and
}\;s_{\ell}(x):=\left|\left|\widetilde{\ell}(x)\right|\right|_{0}\ll\min\\{n,d\\}.$

######  Remark 2.

Under the stronger Assumption 4, $\mathrm{E}(RXX^{T})$ also satisfies the
$\mathcal{C}_{1}(S_{\ell}(x),3)$-restricted eigenvalue condition with some
parameter $\kappa_{R}^{2}>0$ defined in (16).

###  3.3 Consistency Results

In this section, we unroll the consistency results of our Lasso pilot estimate
(2) and the solution to the dual program (6). Additionally, we discuss the
strong duality theory of our debiasing inference method and validate the
asymptotic negligibility of the “Conditional variance II” term in Proposition
2 of Section 2.3.

######  Theorem 1 (Consistency of the Lasso Pilot Estimate with Complete-Case
Data).

Let $\delta\in(0,1)$ and $A_{1},A_{2}>0$ be some absolute constants. Suppose
that Assumptions 1, 2, and 4 (more precisely, Eq.(16)) hold and $n$ is
sufficiently large so that

$\sqrt{\frac{s_{\beta}\log(ed/s_{\beta})}{n}}+\sqrt{\frac{\log(1/\delta)}{n}}<\min\left\\{\frac{\kappa_{R}^{2}}{A_{1}\phi_{s_{\beta}}},1\right\\}.$

  1. (a)

If
$\lambda>\frac{4}{n}\left|\left|\sum_{i=1}^{n}R_{i}X_{i}\epsilon_{i}\right|\right|_{\infty}>0$
and $d\geq s_{\beta}+2$, then with probability at least $1-\delta$,

$\left|\left|\widehat{\beta}-\beta_{0}\right|\right|_{2}\lesssim\frac{\sqrt{s_{\beta}}\lambda}{\kappa_{R}^{2}}.$

  2. (b)

If, in addition, Assumption 3 holds and
$\lambda=A_{2}\sigma_{\epsilon}\sqrt{\frac{\phi_{1}\log(d/\delta)}{n}}$, then
with probability at least $1-\delta$,

$\displaystyle\left|\left|\widehat{\beta}-\beta_{0}\right|\right|_{2}$
$\displaystyle\lesssim\frac{\sigma_{\epsilon}}{\kappa_{R}^{2}}\sqrt{\frac{\phi_{1}s_{\beta}\log(d/\delta)}{n}}.$

######  Remark 3.

If $\log d=o(n)$ with $d>n$ and $\sigma_{\epsilon},\phi_{1}>0$ are independent
of $n$ and $d$, then we obtain the asymptotic consistency result
$\left|\left|\widehat{\beta}-\beta_{0}\right|\right|_{2}=O_{P}\left(\frac{1}{\kappa_{R}^{2}}\sqrt{\frac{s_{\beta}\log
d}{n}}\right)$ by setting $\delta=\frac{1}{d}$.

The proof of Theorem 1 is in Section E.3. It implies that the Lasso pilot
estimate $\widehat{\beta}$ in (2) based on the complete-case data is
consistent. This property is a consequence of our MAR assumption and the lower
eigenvalue bound condition on $\mathrm{E}(RXX^{T})$ (Assumption 4), because
the conditional distribution of $Y$ given $X$ remains unchanged on the
complete-case data ($R=1$) and the loss function in (2) is “strictly convex”
with the cone in which $\widehat{\beta}$ lies.

Next, we present the consistency result for the dual solution
$\widehat{\ell}(x)$ to the sparse $r_{\ell}$-approximation
$\widetilde{\ell}(x)$. Among other things, the rate of consistency for
$\widehat{\ell}(x)$ depends on the rate of consistency
$r_{\pi}=r_{\pi}(n,\delta)$ for estimating propensity scores.

######  Theorem 2 (Consistency of the dual solution $\widehat{\ell}(x)$).

Let $\delta\in(0,1)$, $x\in\mathbb{R}^{d}$ be a fixed covariate vector, and
$A_{1},A_{2},A_{3}>0$ be some large absolute constants. Suppose that
Assumptions 1, 2, 4, 5, and 6 hold as well as

$\sqrt{\frac{s_{\ell}(x)\log(ed/s_{\ell}(x))}{n}}+\sqrt{\frac{\log(1/\delta)}{n}}+r_{\pi}<\min\left\\{\frac{\kappa_{R}^{2}}{A_{1}\phi_{s_{\ell}(x)}},1\right\\}\quad\text{
and }\quad d\geq s_{\ell}(x)+2.$

  1. (a)

Suppose that
$\gamma\geq\left|\left|\frac{5}{2}\sum_{i=1}^{n}\widehat{\pi}_{i}X_{i}X_{i}^{T}\ell_{0}(x)+5nx\right|\right|_{\infty}$
and
$r_{\ell}<\min\left\\{\frac{\kappa_{R}^{2}}{\left|\left|x\right|\right|_{2}}\sqrt{\frac{2\gamma}{13\nu}},\,\frac{1}{2}\right\\}$,
where $\nu>0$ is defined in (58) of Lemma 9. Then, with probability at least
$1-\delta$,

$\left|\left|\widehat{\ell}(x)-\widetilde{\ell}(x)\right|\right|_{2}\lesssim\frac{\gamma\sqrt{s_{\ell}(x)}}{n\cdot\kappa_{R}^{2}}+\frac{\phi_{s_{\ell}(x)}\left|\left|x\right|\right|_{2}}{\kappa_{R}^{4}}\cdot
r_{\ell}.$

  2. (b)

We denote $r_{\gamma}\equiv
r_{\gamma}(n,\delta)=\frac{\sqrt{\phi_{1}}}{\kappa_{R}}\left|\left|x\right|\right|_{2}\sqrt{\frac{\log(d/\delta)}{n}}+\frac{\sqrt{\phi_{1}\phi_{s_{\ell}(x)}}\left|\left|x\right|\right|_{2}}{\kappa_{R}^{2}}\cdot
r_{\pi}(n,\delta)$. If, in addition, $\frac{\gamma}{n}=A_{2}\cdot r_{\gamma}$
and

$r_{\ell}\leq\min\left\\{\frac{1}{2},\,\left[\left(\frac{A_{2}^{2}\,\kappa_{R}^{6}\,\phi_{1}}{A_{3}^{2}\left|\left|x\right|\right|_{2}\phi_{s_{\ell}(x)}}\right)^{\frac{1}{4}}\left(\frac{\log(d/\delta)}{n}\right)^{\frac{1}{4}}+\left(\frac{A_{2}^{2}\phi_{1}\kappa_{R}^{4}}{A_{3}^{2}\phi_{s_{\ell}(x)}\left|\left|x\right|\right|_{2}^{2}}\right)^{\frac{1}{4}}\sqrt{r_{\pi}}\right]\right\\},$

then with probability at least $1-\delta$,

$\displaystyle\left|\left|\widehat{\ell}(x)-\widetilde{\ell}(x)\right|\right|_{2}$
$\displaystyle\lesssim\frac{\left|\left|x\right|\right|_{2}}{\kappa_{R}^{3}}\left[\sqrt{\frac{\phi_{1}s_{\ell}(x)\log(d/\delta)}{n}}+\frac{\phi_{s_{\ell}(x)}}{\kappa_{R}}\cdot
r_{\ell}+\frac{\sqrt{\phi_{1}\phi_{s_{\ell}(x)}s_{\ell}(x)}}{\kappa_{R}}\cdot
r_{\pi}\right].$

######  Remark 4.

If $\log d=o(n)$ with $d>n$, $\left|\left|x\right|\right|_{2}=O(1)$, and
$\phi_{1},\phi_{s_{\ell}(x)}$ are independent of $n$ and $d$, then we obtain
the asymptotic consistency result

$\left|\left|\widehat{\ell}(x)-\widetilde{\ell}(x)\right|\right|_{2}=O_{P}\left(\frac{1}{\kappa_{R}^{3}}\sqrt{\frac{s_{\ell}(x)\log
d}{n}}+\frac{r_{\ell}}{\kappa_{R}^{4}}+\frac{r_{\pi}\sqrt{s_{\ell}(x)}}{\kappa_{R}^{4}}\right)$

by setting $\delta=\frac{1}{d}$. This asymptotic rate of convergence comprises
two parts. The first part is of the order
$O_{P}\left(\frac{1}{\kappa_{R}^{3}}\sqrt{\frac{s_{\ell}(x)\log
d}{n}}+\frac{r_{\ell}}{\kappa_{R}^{4}}\right)$, which is the oracle rate of
convergence for
$\left|\left|\widehat{\ell}(x)-\widetilde{\ell}(x)\right|\right|_{2}$ when the
true propensity scores are known; see Lemma 1 of Section E.4. The second part
is of the order
$O_{P}\left(\frac{r_{\pi}\sqrt{s_{\ell}(x)}}{\kappa_{R}^{2}}\right)$, which is
related to the consistency of the propensity score estimation.

The proof of Theorem 2 is in Section E.4. An immediate consequence of this
result is that the asymptotic rate of convergence of
$\left|\left|\widehat{\ell}(x)-\ell_{0}(x)\right|\right|_{2}$ is the same as
in Remark 4, because by triangle inequality,

$\displaystyle\left|\left|\widehat{\ell}(x)-\ell_{0}(x)\right|\right|_{2}$
$\displaystyle\leq\left|\left|\widehat{\ell}(x)-\widetilde{\ell}(x)\right|\right|_{2}+\left|\left|\widetilde{\ell}(x)-\ell_{0}(x)\right|\right|_{2}\lesssim\left|\left|\widehat{\ell}(x)-\widetilde{\ell}(x)\right|\right|_{2}+\frac{\left|\left|x\right|\right|_{2}}{\kappa_{R}^{2}}\cdot
r_{\ell},$

and the last term on the right-hand side is asymptotically negligible compared
to the second term
$\frac{\phi_{s_{\ell}(x)}\left|\left|x\right|\right|_{2}}{\kappa_{R}^{2}}\cdot
r_{\ell}$ in Theorem 2(b). This implies that the dual solution
$\widehat{\ell}(x)$ is consistent to the population dual solution
$\ell_{0}(x)$ under the setup of Remark 4 if

$\frac{1}{\kappa_{R}^{3}}\sqrt{\frac{s_{\ell}(x)\log d}{n}}\to
0,\quad\frac{r_{\ell}}{\kappa_{R}^{4}}\to
0,\quad\frac{r_{\pi}\sqrt{s_{\ell}(x)}}{\kappa_{R}^{4}}\to 0\quad\text{ as
}\quad n\to\infty.$

In particular, we only need to consistently estimate the propensity scores,
while the rate of consistency $r_{\pi}$ is essentially irrelevant. Hence, we
are completely free to choose whatever consistent estimator suits us best,
including those from modern machine learning methods such as deep neural
networks (Farrell et al., 2021), random forests (Gao et al., 2022), and
support vector machines with universal kernels (such as the Gaussian radial
basis function) (Steinwart, 2001). For concreteness, we present in Proposition
1 of Section C that under comparatively mild regularity conditions, the Lasso-
type generalized linear regression estimator (26) is a feasible option.

The assumptions in Theorem 2 are also sufficient to establish the strong
duality between the primal debiasing program (3) and its dual (6):

######  Lemma 3 (Sufficient conditions for strong duality).

Let $\delta\in(0,1)$. Under the assumptions of Theorem 2(b), there exists
$\gamma>0$ such that $\frac{\gamma}{n}\asymp r_{\gamma}\equiv
r_{\gamma}(n,\delta)$ and the strong duality as well as the relation (7) in
Proposition 1 hold with probability at least $1-\delta$.

The consistency results for the Lasso pilot estimate $\widehat{\beta}$ and the
dual solution $\widehat{\ell}(x)$ also imply that under mild regularity
conditions, the “Conditional variance II” term in Proposition 2 is
asymptotically negligible. The proof of Lemma 4 is in Section E.5. This result
completes the motivation of the debiased program (3) from the perspective of
(conditional) bias-variance trade-off in Section 2.3.

######  Lemma 4 (Negligible conditional variance term).

Suppose that Assumptions 1, 2, 3, 4, 5, and 6 hold as well as

$\left[s_{\beta}\log(d/s_{\beta})+s_{\ell}(x)\log(d/s_{\ell}(x))\right]^{2}=O\left(\sqrt{n}\right),\quad\frac{\left|\left|x\right|\right|_{2}}{\kappa_{R}^{2}}\sqrt{\phi_{2s_{\ell}(x)}\phi_{2s_{\beta}}}=O(1),\quad\text{
and }\quad r_{\ell}=O(1).$

If $\left|\left|\widehat{\beta}-\beta_{0}\right|\right|_{2}=o_{P}(1)$ and
$\left|\left|\widehat{\ell}(x)-\widetilde{\ell}(x)\right|\right|_{2}=o_{P}(1)$,
then the “Conditional variance II” term in Proposition 2 is asymptotically
negligible in the sense that

$\displaystyle\left(\beta_{0}-\widehat{\beta}\right)^{T}\left[\sum_{i=1}^{n}\widehat{w}_{i}^{2}(x)\widehat{\pi}_{i}\left(1-\widehat{\pi}_{i}\right)X_{i}X_{i}^{T}\right]\left(\beta_{0}-\widehat{\beta}\right)=o_{P}(1).$

###  3.4 Asymptotic Normality

In this section, we leverage the consistency results from Section 3.3 to
establish the asymptotic normality of the debiased estimator
$\widehat{m}^{\text{debias}}(x;\widehat{\bm{w}})$ in (4). As mentioned in
Section 3.1, the key observation for applying the CLT is the linear expansion
(3.1) based on the dual solution.

######  Theorem 5 (Asymptotic normality of the debiased estimator).

Let $x\in\mathbb{R}^{d}$ be a fixed query point and
$s_{\max}=\max\left\\{s_{\beta},s_{\ell}(x)\right\\}$. Suppose that
Assumptions 1, 2, 3, 4, 5, and 6 hold and $\lambda,\gamma,r_{\ell}>0$ are
specified as in Theorem 1(b) and Theorem 2(b), respectively. Furthermore, we
assume that
$\sigma_{\epsilon}\phi_{1}\phi_{s_{\max}}^{3/2}\left|\left|x\right|\right|_{2}=O(1)$,

$\displaystyle\frac{(1+\kappa_{R}^{2})s_{\max}\log(nd)}{\kappa_{R}^{4}}=o\left(\sqrt{n}\right),\quad\text{
and
}\quad\frac{(1+\kappa_{R}^{4})\sqrt{s_{\max}\log(nd)}}{\kappa_{R}^{6}}\,(r_{\ell}+r_{\pi})=o(1).$
(20)

Then,

$\sqrt{n}\left[\widehat{m}^{\mathrm{debias}}(x;\widehat{\bm{w}})-m_{0}(x)\right]\stackrel{{\scriptstyle
d}}{{\to}}\mathcal{N}\left(0,\,\sigma_{m}^{2}(x)\right).$

The proof of Theorem 5 can be found in Section E.6, and we make two remarks
for this paramount result.

######  Remark 5 (Semi-parametric efficiency).

Given any fixed dimension $d>0$, the asymptotic variance

$\sigma_{m}^{2}(x)=\sigma_{\epsilon}^{2}\cdot
x^{T}\left[\mathrm{E}\left(RXX^{T}\right)\right]^{-1}x=\sigma_{\epsilon}^{2}\cdot
x^{T}\left\\{\mathrm{E}\left[\pi(X)XX^{T}\right]\right\\}^{-1}x$

for our debiased estimator achieves the semi-parametric efficiency bound among
all asymptotically linear estimators with the MAR outcome (Müller and
Keilegom, 2012). It also attains the efficiency bound for de-sparsified Lasso
in Theorem 3 of Jankova and van de Geer (2018) under Gaussian covariates and
noises. This key property demonstrates the superiority of our debiasing
inference method in terms of statistical efficiency.

######  Remark 6 (Necessary sparsity condition).

Compared with the consistency results in Theorem 2, we impose a stricter
growth condition $s_{\max}=o\left(\frac{\sqrt{n}}{\log d}\right)$ on the
sparsity level when $d>n$ in order to establish the asymptotic normality,
provided that other model parameters
$\sigma_{\epsilon},\phi_{s_{\max}},\kappa_{R}$ are independent of $n$ and $d$.
Such a requirement on the sparsity level is a standard and essentially
necessary condition in the high-dimensional debiased inference literature in
pursuit of asymptotic normality (van de Geer et al., 2014; Javanmard and
Montanari, 2014a, 2018; Cai et al., 2021); see also the discussion in Section
8.6 of Jankova and van de Geer (2018).

In order to utilize the asymptotic normality of Theorem 5 and conduct
inference on $m_{0}(x)$ in practice, we need to estimate the asymptotic
variance $\sigma_{m}^{2}(x)$. The following Proposition 6 verifies that the
objective function $\sum_{i=1}^{n}\widehat{\pi}_{i}\widehat{w}_{i}^{2}(x)$ in
our debiasing program (3) is a consistent estimator of $\sigma_{m}^{2}(x)$
under the known noise level $\sigma_{\epsilon}^{2}$; see Section E.7 for its
proof. This result again motivates the design of our debiasing program as in
Section 2.3 and supports the formulation of our asymptotically valid
confidence interval (5). When the noise level $\sigma_{\epsilon}^{2}$ is
unknown in practice, one can always replace it with a consistent estimator
$\widehat{\sigma}_{\epsilon}^{2}$ without losing any asymptotic efficiency and
obtain the final estimator
$\widehat{\sigma}_{\epsilon}^{2}\sum_{i=1}^{n}\widehat{\pi}_{i}\widehat{w}_{i}^{2}(x)$
for the asymptotic variance. There are various approaches available in the
literature to construct such an estimator for $\sigma_{\epsilon}^{2}$ (Zhang,
2010; Fan et al., 2012; Dicker, 2012; Bayati et al., 2013; Belloni and
Chernozhukov, 2013; Reid et al., 2016). In our simulation studies and real-
world application, we use the scaled Lasso (Sun and Zhang, 2012) to obtain
$\widehat{\sigma}_{\epsilon}^{2}$; see Section A for details.

######  Proposition 6 (Consistent estimate of the asymptotic variance).

Let $x\in\mathbb{R}^{d}$ be a fixed query point. Suppose that Assumptions 1,
2, 4, 5, and 6 hold and $\gamma,r_{\ell}>0$ are specified as in Theorem 2(b).
Furthermore, we assume that
$\left|\left|x\right|\right|_{2}^{2}\phi_{s_{\ell}(x)}^{2}=O(1)$,
$\frac{(1+\kappa_{R}^{3})}{\kappa_{R}^{5}}\sqrt{\frac{s_{\ell}(x)\log(nd)}{n}}=o(1)$,
and
$\frac{(1+\kappa_{R}^{4})}{\kappa_{R}^{6}}\left[r_{\ell}+r_{\pi}\sqrt{s_{\ell}(x)}\right]=o(1)$.
Then,

$\displaystyle\left|\sum_{i=1}^{n}\widehat{\pi}_{i}\widehat{w}_{i}^{2}(x)-x^{T}\left[\mathrm{E}\left(RXX^{T}\right)\right]^{-1}x\right|=o_{P}(1).$

##  4 Experiments

In this section, we evaluate the empirical performance of our proposed
debiasing inference method in Section 2.1 and compare it with some existing
high-dimensional inference methods in the literature through comprehensive
simulation studies.

###  4.1 Basic Simulation Design and Methods to Be Compared

We generate the i.i.d. data
$\\{(Y_{i},R_{i},X_{i})\\}_{i=1}^{n}\subset\mathbb{R}\times\\{0,1\\}\times\mathbb{R}^{d}$
from the following linear model

$Y_{i}=X_{i}^{T}\beta_{0}+\epsilon_{i}\quad\text{ with }\quad
X_{i}{\perp\\!\\!\\!\perp}\epsilon_{i},\quad
Y_{i}{\perp\\!\\!\\!\perp}R_{i}|X_{i},\quad\text{ and }\quad
X_{i}\sim\mathcal{N}_{d}(\bm{0},\Sigma),$ (21)

where $d=1000$ and $n=900$ unless stated otherwise. In order to investigate
the effect of different covariance structures $\Sigma\in\mathbb{R}^{d\times
d}$, sparsity patterns of the true regression vector
$\beta_{0}\in\mathbb{R}^{d}$, and designs of the query point
$x\in\mathbb{R}^{d}$ on the finite-sample performances of the candidate high-
dimensional inference methods, we consider all possible combinations of the
following simulation designs.

Choices of the covariance matrix: We consider two standard designs of
$\Sigma=\left(\Sigma_{jk}\right)_{j,k=1}^{d}$ from the literature:

  1. (i)

The circulant symmetric matrix $\Sigma^{\mathrm{cs}}$ in Javanmard and
Montanari (2014a) defined by $\Sigma_{jj}=1$, $\Sigma_{jk}=0.1\text{ when
}j+1\leq k\leq j+5\text{ or }j+d-5\leq k\leq j+d-1$ with $\Sigma_{jk}=0$
elsewhere for $j\leq k$, and $\Sigma_{jk}=\Sigma_{kj}$;

  2. (ii)

The Toeplitz (or auto-regressive) matrix $\Sigma^{\mathrm{ar}}$ in van de Geer
et al. (2014) defined by $\Sigma_{jk}=0.9^{|j-k|}$.

Designs of the regression vector: For the true value of
$\beta_{0}\in\mathbb{R}^{d}$, we consider three different scenarios as:

  1. (i)

$\beta_{0}^{sp}=\left(\underbrace{\sqrt{5},...,\sqrt{5}}_{5},0,...,0\right)^{T}\in\mathbb{R}^{d}$
is sparse;

  2. (ii)

$\beta_{0}^{de}\propto\left(1,\frac{1}{\sqrt{2}},...,\frac{1}{\sqrt{d}}\right)^{T}\in\mathbb{R}^{d}$
is dense with $\left|\left|\beta_{0}^{de}\right|\right|_{2}=5$;

  3. (iii)

$\beta_{0}^{pd}\propto\left(1,\frac{1}{2},...,\frac{1}{d}\right)^{T}\in\mathbb{R}^{d}$
is pseudo-dense $\left|\left|\beta_{0}^{pd}\right|\right|_{2}=5$.

Designs of the query point: We conduct experiments on the following four
different choices of the query point $x\in\mathbb{R}^{d}$:

  1. (i)

$x^{(1)}=(1,0,...,0)^{T}\in\mathbb{R}^{d}$ to infer the individual effect of
the first (significantly nonzero) component of $\beta_{0}$;

  2. (ii)

$x^{(2)}=\left(1,\frac{1}{2},\frac{1}{4},0,0,0,\frac{1}{2},\frac{1}{8},0,...,0\right)^{T}\in\mathbb{R}^{d}$
to infer the joint effects of a few components of $\beta_{0}$;

  3. (iii)

$x^{(3)}=\left(0,...,0,\underbrace{1}_{100^{th}},0,...,0\right)^{T}\in\mathbb{R}^{d}$
to infer an inactive (_i.e._ , truly zero or close to zero) component of
$\beta_{0}$;

  4. (iv)

$x^{(4)}=\left(1,\frac{1}{2^{2}},...,\frac{1}{d^{2}}\right)^{T}\in\mathbb{R}^{d}$
for the circulant symmetric covariance cases, and
$x^{(4)}=\left(1,\frac{1}{2},...,\frac{1}{d}\right)^{T}\in\mathbb{R}^{d}$ for
the Toeplitz covariance cases. The purpose of choosing a relatively dense
query point $x^{(4)}$ is to study the joint effect of all the components of
$\beta_{0}$.

MAR mechanism: After $Y_{i},i=1,...,n$ are sampled according to (21), we
define the missingness indicators $R_{i},i=1,...,n$ for $Y_{i},i=1,...,n$
through the MAR mechanism by

$\mathrm{P}(R_{i}=1|X_{i})=\frac{1}{1+\exp\left(-1+X_{i7}-X_{i8}\right)}\quad\text{
for }\quad i=1,...,n.$ (22)

In general, the above MAR mechanism yields around 28% of missingness for the
outcome variables $Y_{i},i=1,...,n$, making the complete-case data even more
high-dimensional. Additional simulation results under a simpler MCAR setting
are deferred to Section D.2.

Noise distribution: We first generate the noise variables
$\epsilon_{i},i=1,...,n$ independently from $\mathcal{N}(0,1)$. Other types of
the noise distributions that violate the sub-gaussian noise condition
(Assumption 3) are considered in Section 4.3.

Methods to be compared: To demonstrate the statistical efficiency of our
proposed debiasing method (“Debias”), we compare it with several existing
inference methods in the literature. The detailed implementation of our
proposed method can be found in Section A. There we also explain the three
different criteria “min-CV”, “1SE”, and “min-feas” for selecting the tuning
parameter $\gamma>0$ via cross-validation. We implement all of the following
comparative methods on (i) the complete-case (CC) data
$(X_{i},Y_{i})\in\mathbb{R}^{d}\times\mathbb{R}$ with $R_{i}=1$ for
$i=1,...,n$, (ii) the inverse probability weighted (IPW) data
$\left(\frac{X_{i}}{\sqrt{\widehat{\pi}_{i}}},\frac{Y_{i}}{\sqrt{\widehat{\pi}_{i}}}\right)\in\mathbb{R}^{d}\times\mathbb{R}$
with $R_{i}=1$ for $i=1,...,n$, and (iii) the oracle fully observed data
$(X_{i},Y_{i}),i=1,...,n$. Here, $\widehat{\pi}_{i},i=1,...,n$ are estimated
propensity scores by the Lasso-type logistic regression as described in
Section C.

  1. (i)

“DL-Jav” is the debiased Lasso proposed by Javanmard and Montanari (2014a)
that solves $d$ convex programs for an unbiased estimator for $\beta_{0}$ and
a consistent estimate for the asymptotic covariance matrix. The method is
implemented by the source code sslasso provided in their paper. We select the
tuning parameters as their Section 5.1.

  2. (ii)

“DL-vdG” is the debiased Lasso proposed by van de Geer et al. (2014) that
solves $d$ nodewise regressions for the debiased estimator of $\beta_{0}$ and
the estimate of its asymptotic covariance matrix. This method is implemented
via the function lasso.proj in the R package hdi (Dezeure et al., 2015).

  3. (iii)

“R-Proj” is the ridge projection method proposed by Bühlmann (2013) that
produces a bias-corrected estimator of $\beta_{0}$ via a ridge regression
estimator. This method is implemented by the function ridge.proj in the R
package hdi. We initialize the estimation of $\beta_{0}$ with an estimate
obtained from the scaled Lasso. Since it is difficult to obtain a consistent
estimate of the asymptotic covariance matrix from the function ridge.proj, we
only run “R-Proj” when the query point is $x^{(1)}$ or $x^{(3)}$. This means
that we only conduct inference on single coordinates of $\beta_{0}$ under
these scenarios.

  4. (iv)

“Refit” first runs the scaled Lasso to obtain a pilot estimate
$\widehat{\beta}$ and then computes a final the least-square estimate based on
covariates in the support set of $\widehat{\beta}$ (Belloni and Chernozhukov,
2013).

Confidence intervals for “DL-Jav”, “DL-vdG”, and “R-Proj” are constructed
according to their asymptotic normality results, while the confidence
intervals from “Refit” are based on the assumption that the selected model by
the Lasso pilot estimate $\widehat{\beta}$ contains the true support set.
Contrary to our asymptotic normality theory in Section 3.4, these estimators
need not be semi-parametrically efficient.

###  4.2 Simulation Results Under the Gaussian Noise Distribution

We evaluate the performances of our proposed debiasing method and the existing
inference methods stated above via the absolute bias
$\left|\widehat{m}^{\mathrm{debias}}(x;\widehat{\bm{w}})-m_{0}(x)\right|$, the
average coverage of confidence intervals with the 95% nominal coverage
probability, and the average length of confidence intervals across 1000 Monte
Carlo experiments for each simulation scenario stated in Section 4.1. Since it
is very time-consuming to run “DL-vdG”, we only execute 500 Monte Carlo
experiments for this method under each scenario.

Figure 2: Simulation results with sparse $\beta_{0}^{sp}$ and sparse $x^{(2)}$
for the circulant symmetric covariance matrix $\Sigma^{\mathrm{cs}}$ under the
Gaussian noise $\mathcal{N}(0,1)$ and the MAR setting. Top Left: Boxplots of
the absolute bias. Top Right: Coverage probabilities with standard error bars.
Bottom Left: Average lengths of confidence intervals with standard error bars.
Bottom Right: Four comparative “QQ-plots” of the transformed studentized
debiased estimators
$\Phi\left(\sqrt{n}\cdot\widehat{\sigma}_{n}^{-1}(x)\left[\widehat{m}(x)-m_{0}(x)\right]\right)$
obtained from different debiasing methods, where $\widehat{\sigma}_{n}(x)^{2}$
is the estimated (asymptotic) variance of $\widehat{m}(x)$ by each method and
$\Phi(\cdot)$ is the CDF of $\mathcal{N}(0,1)$. We also highlight the
recommended rules “1SE” and “min-feas” for our proposed debiasing method via
dashed rectangles in the first three panels.

We provide comparative plots in terms of the three evaluation metrics and
normality test on selected simulation settings in Figure 2 and Figure 3 while
deferring the full simulation results for the circulant symmetric and Toeplitz
(or auto-regressive) covariance design to Table 1 and Table 2 in Section D.1,
respectively. In summary, our proposed debiasing method removes roughly the
same amount of bias as the debiased Lasso “DL-Jav” and the Lasso refitting
“Refit” but yields relatively shorter confidence intervals and more accurate
coverage probabilities than these two methods. While the ridge projection
“R-Proj” sometimes produces better coverage probabilities than other methods,
it always has larger biases and more conservative confidence intervals.

The superiority of our debiasing method is better revealed when the query
point $x\in\mathbb{R}^{d}$ has many nonzero entries. Since the competing
methods are asymptotically normal only coordinate-wise, they suffer from
substantial biases and inflated variances if applied to a dense query point.
In contrast, the debiased estimator from our method is asymptotically semi-
parametrically efficient for arbitrary $x$ in Theorem 5 and its asymptotic
variance can be naturally estimated in Proposition 6. Hence, the asymptotic
normality of our debiased estimator can better exemplify in finite samples
than the competing methods. All these simulation results consolidate the
derived theoretical properties of our debiasing method in Section 3 and
demonstrate its effectiveness in conducting valid statistical inference under
missing outcomes.

Figure 3: Simulation results with pseudo-dense $\beta_{0}^{pd}$ and sparse
$x^{(2)}$ for the Toeplitz (or auto-regressive) covariance matrix
$\Sigma^{\mathrm{ar}}$ under the Gaussian noise $\mathcal{N}(0,1)$ and the MAR
setting. Top Left: Boxplots of the absolute bias. Top Right: Coverage
probabilities with standard error bars. Bottom Left: Average lengths of
confidence intervals with standard error bars. Bottom Right: Four comparative
“QQ-plots” of the transformed studentized debiased estimators
$\Phi\left(\sqrt{n}\cdot\widehat{\sigma}_{n}^{-1}(x)\left[\widehat{m}(x)-m_{0}(x)\right]\right)$
obtained from different debiasing methods, where $\widehat{\sigma}_{n}(x)^{2}$
is the estimated (asymptotic) variance of $\widehat{m}(x)$ by each method and
$\Phi(\cdot)$ is the CDF of $\mathcal{N}(0,1)$. We also highlight the
recommended rules “1SE” and “min-feas” for our proposed debiasing method via
dashed rectangles in the first three panels.

###  4.3 Simulation Results Under Heavy-Tailed Noise Distributions

The consistency and asymptotic normality results for our debiasing inference
methods are derived under the sub-gaussian noise condition (Assumption 3), and
our previous simulation results in Section 4.2 are conducted with normally
distributed noises accordingly. We showcase through additional simulations
that the asymptotic normality of our debiased estimator (4) is valid even when
the distribution of the noise $\epsilon$ in model (1) is not sub-gaussian. To
this end, we follow the simulation designs in Section 4.1 but explore two
heavy-tailed noise distributions as:

  1. (i)

$\epsilon$ follows a Laplace distribution with mean 0 and scale parameter
$\frac{1}{\sqrt{2}}$. Note that $\epsilon$ becomes sub-exponentially
distributed and has variance $\mathrm{Var}(\epsilon)=1$.

  2. (ii)

$\epsilon$ follows a mean-zero $t$-distribution with 2 degrees of freedom. In
this case, $\epsilon$ is not even sub-exponential and has infinite variance.

Since “DL-vdG” is time-consuming to run and follows a similar debiasing
approach as “DL-Jav”, we choose not to implement “DL-vdG” on the above two
heavy-tailed distributions of $\epsilon$. We present the selected plots and
asymptotic normality validations in Figure 4 and Figure 5, while again
deferring the full comparative simulation results for the circulant symmetric
and Toeplitz covariance designs to Table 3-4 and Table 5-6 in Section D.1,
respectively. In short, the debiased estimator appears to be approximately
normally distributed in finite samples even when the noise distribution is
heavy-tailed, while other competing methods fail to maintain their asymptotic
normality. Other conclusions from these simulation results are similar to what
we have described in Section 4.2.

Figure 4: Simulation results with dense $\beta_{0}^{de}$ and sparse $x^{(4)}$
for the circulant symmetric covariance matrix $\Sigma^{\mathrm{cs}}$ under the
$\mathrm{Laplace}\left(0,\frac{1}{\sqrt{2}}\right)$ distributed noises and the
MAR setting as in Section 4.1. Top Middle: Coverage probabilities with
standard error bars. Top Right: Average lengths of confidence intervals with
standard error bars. Middle Panels: Three comparative “QQ-plots” of the
transformed studentized debiased estimators
$\Phi\left(\sqrt{n}\cdot\widehat{\sigma}_{n}^{-1}(x)\left[\widehat{m}(x)-m_{0}(x)\right]\right)$
obtained from different debiasing methods, where $\widehat{\sigma}_{n}(x)^{2}$
is the estimated (asymptotic) variance of $\widehat{m}(x)$ by each method and
$\Phi(\cdot)$ is the CDF of $\mathcal{N}(0,1)$.

Figure 5: Simulation results with pseudo-dense $\beta_{0}^{pd}$ and dense
$x^{(4)}$ for the Toeplitz (or auto-regressive) covariance matrix
$\Sigma^{\mathrm{ar}}$ under the $t_{2}$ distributed noises and the MAR
setting as in Section 4.1. Top Left: Boxplots of the absolute bias. Top
Middle: Coverage probabilities with standard error bars. Top Right: Average
lengths of confidence intervals with standard error bars. Middle Panels: Three
comparative “QQ-plots” of the transformed studentized debiased estimators
$\Phi\left(\sqrt{n}\cdot\widehat{\sigma}_{n}^{-1}(x)\left[\widehat{m}(x)-m_{0}(x)\right]\right)$
obtained from different debiasing methods, where $\widehat{\sigma}_{n}(x)^{2}$
is the estimated (asymptotic) variance of $\widehat{m}(x)$ by each method and
$\Phi(\cdot)$ is the CDF of $\mathcal{N}(0,1)$.

###  4.4 Simulation Results for the Proposed Debiasing Method With
Nonparametrically Estimated Propensity Scores

As shown in Theorem 2 and Theorem 5, the consistency and asymptotic normality
of the proposed debiasing method do not require any parametric assumption on
the propensity score $\mathrm{P}(R=1|X)$ under MAR condition. Therefore, we
now conduct additional simulation studies for our debiasing method in which
the propensity scores $\pi(X_{i}),i=1,...,n$ estimated by the following
nonlinear/nonparametric machine learing methods based on the observed data
$\\{(X_{i},R_{i})\\}_{i=1}^{n}$. All of these methods are implemented in the
scikit-learn package (Pedregosa et al., 2011) in Python.

Naive Bayes (“NB”): We adapt the Gaussian naive Bayes method (Zhang, 2004) to
model the propensity scores through

$\pi(X_{i})=\mathrm{P}(R_{i}|X_{i})=\mathrm{P}(R_{i}|X_{i1},...,X_{in})=\frac{\mathrm{P}(R_{i})\cdot\prod_{k=1}^{n}\mathrm{P}(X_{ik}|R_{i})}{\mathrm{P}(X_{i1},...,X_{in})}$

for any $i=1,...,n$, where $\mathrm{P}(X_{ik}|R_{i})$ is the density of
$\mathcal{N}(\mu_{R},\sigma_{R}^{2})$ and $\mu_{R},\sigma_{R}$ are estimated
by the maximum likelihood approach.

Random Forest (“RF”): We implement the random forest method (Breiman, 2001)
with 100 trees, bootstrapping samples, and the Gini impurity to measure the
quality of a split.

Support Vector Machine (“SVM”): We fit a support vector machine (Chen et al.,
2005) with the Gaussian radial basis function to classify the missingness
indicator $R_{i}$ based on the observed covariate vector
$X_{i}\in\mathbb{R}^{d}$ for $i=1,...,n$. The outputs of the trained SVM are
regarded as the surrogates of primitive estimated propensity scores.

Neural Network (“NN”): We train a neural network model (Hinton, 1990) with two
hidden layers of size $80\times 50$ and use the rectified linear unit function
$h(x)=\max\\{x,0\\}$ as the activation function. The learning rate
$\eta_{\text{NN}}$ is initially set to $0.001$ as long as the training loss
keeps decreasing and is adaptively changed to $\eta_{\text{NN}}/5$ when two
consecutive epochs fail to decrease the training loss.

For the above methods, we also consider their calibrated versions through the
Platt’s logistic model (Platt, 1999) that fits an extra logistic regression on
the estimated propensity scores $\widehat{\pi}_{i},i=1,...,n$. However, since
the errors of estimated propensity scores and the final performances of the
proposed debiased estimator worsens after we apply the calibration to the
estimated propensity scores, we choose not to report them here.

The covariate vectors $X_{i}\in\mathbb{R}^{d},i=1,...,n$ are again sampled
independently from $\mathcal{N}_{d}(\bm{0},\Sigma^{\mathrm{cs}})$ with
$d=1000$ and $n=900$, and the noise variables $\epsilon_{i},i=1,...,n$ are
generated independently from $\mathcal{N}(0,1)$ as in Section 4.1. To increase
the complexity of estimating the propensity scores, we generate the missing
indicators $R_{i},i=1,...,n$ for the outcome variables $Y_{i},i=1,...,n$
through a different MAR mechanism than the one in Section 4.1 as:

$\mathrm{P}(R_{i}=1|X_{i})=\Phi\left(-4+\sum_{k=1}^{K}Z_{ik}\right),$ (23)

where $\Phi(\cdot)$ is the CDF of $\mathcal{N}(0,1)$ and the vector
$(Z_{i1},...,Z_{iK})$ contains all polynomial combinations of the first eight
components $X_{i1},...,X_{i8}$ of the covariate vector $X_{i}$ with degrees
less than or equal to two (_i.e._ , including the linear, quadratic, and one-
way interaction terms).

The simulation results for our debiasing method under the oracle propensity
scores (“Oracle”), the Lasso-type logistic regression (27) (“LR”), and the
aforementioned nonparametric methods for the propensity score estimation are
shown in Figure 6 and Figure 7. We provide more comprehensive results with the
evaluation metrics from Section 4.2 and an additional measure for the average
mean absolute error (“Avg-MAE”) for the estimated propensity scores in Table 7
as well as additional “QQ-plots” in Figure 9 of Section D.1.

Figure 6: Simulation results for our debiasing method with dense
$\beta_{0}^{de}$ and sparse $x^{(2)}$ when the propensity scores are estimated
by various machine learning methods under the new MAR setting (23). Top Left:
Boxplots of the absolute bias. Top Right: Coverage probabilities with standard
error bars. Middle Left: Average lengths of confidence intervals with standard
deviation bars. Middle Right: Letter-valued plots of the MAE of estimated
propensity scores. Bottom Two Rows: Histograms of the studentized debiased
estimators under different methods. We also order the performance metrics
according to their means in the first four panels.

Figure 7: Simulation results for our debiasing method with sparse
$\beta_{0}^{sp}$ and (weakly) dense $x^{(4)}$ when the propensity scores are
estimated by various machine learning methods under the new MAR setting (23).
Top Left: Boxplots of the absolute bias. Top Right: Coverage probabilities
with standard error bars. Middle Left: Average lengths of confidence intervals
with standard deviation bars. Middle Right: Letter-valued plots of the MAE of
estimated propensity scores. Bottom Two Rows: Histograms of the studentized
debiased estimators under different methods. We also order the performance
metrics according to their means in the first four panels.

The main conclusion from this simulation study is that the performance of our
debiasing method can be further improved when the propensity scores are better
estimated. Under the MAR mechanism (23) and Gaussian design on covariate
vectors $X_{i},i=1,...,n$, neural network and Gaussian naive Bayes models
outperform other machine learning methods in estimating the propensity scores
and lead to the subsequent debiased estimators with lower biases, better
coverages, and shorter confidence intervals in general. Conversely, the Lasso-
type logistic regression is mis-specified and consequently degrades the
performance of the resulting debiased estimator. Another interesting
phenomenon from these simulation results is that our debiasing method
sometimes has better performances when using the estimated propensity scores
$\widehat{\pi}_{i},i=1,...,n$ than plugging in the oracle propensity scores
$\pi_{i}=\pi(X_{i}),i=1,...,n$. Such a paradox has been analyzed in the
propensity score matching (Rosenbaum, 1987), the IPW estimator for causal
inference (Robins et al., 1992; Su et al., 2023) and other general parameter
estimation problems in the presence of a nuisance parameter (Henmi and Eguchi,
2004; Hitomi et al., 2008; Lok, 2021). A rigorous study of this paradox on our
debiasing method is beyond the scope of this paper.

##  5 Real-World Application

We showcase a real application of our proposed debiasing method to the stellar
mass inference on selected galaxies in the Sloan Digital Sky Survey, Fourth
Phase and Data Release 16 (SDSS-IV DR16; Ahumada et al. 2020)111See
https://www.sdss4.org/dr16..

###  5.1 Background and Study Design

SDSS-IV DR16 contains three main survey data: the Extended Baryon Oscillation
Spectroscopic Survey (eBOSS; Dawson et al. 2016), Mapping Nearby Galaxies at
APO (MaNGA; Bundy et al. 2014), and APO Galactic Evolution Experiment 2
(APOGEE-2; Majewski et al. 2017), among which the eBOSS survey and its
predecessor BOSS cover a broad range of cosmological scales and observes most
galaxies in the Universe. While a couple of value added catalogs for the
estimated stellar masses of observed galaxies based on the BOSS spectra are
available (Blanton et al., 2005; Conroy et al., 2009; Chang et al., 2015), we
focus on the eBOSS Firefly value-added catalog (Comparat et al., 2017) 222See
https://www.sdss4.org/dr17/spectro/galaxy_firefly/. for our stellar mass
inference study here given its timeliness and completeness. In more details,
this catalog fits combinations of single-burst stellar population models to
spectroscopic data that have a positive definite redshift and are classified
as galaxies by SDSS-IV DR16, following an iterative best-fitting process
controlled by the Bayesian Information Criterion in order to produce stellar
masses and other stellar properties of the galaxies (Wilkinson et al., 2017).
There are a fraction of missing data in the estimated stellar masses by the
Firefly catalog due to the limiting usage of the observational run in SDSS-IV
DR16 dedicated to galaxy targets, data contamination, and misclassification of
galaxies as stars (Comparat et al., 2017). Given that we will incorporate
various spectroscopic and photometric properties of the observed galaxies into
the covariate vector for our stellar mass inference task, it is reasonable
that the estimated stellar masses are missing at random.

Since nearly all the existing catalogs only yield point estimates to the
stellar mass, we intend to leverage our debiasing inference method to answer
two scientific questions:

  * •

Question 1: How can we produce uncertainty measures (or confidence intervals)
for the estimated stellar mass of a newly observed galaxy based on its
observed spectroscopic and photometric properties?

  * •

Question 2: Is it statistically significant that the stellar mass is
negatively correlated with its distance to the nearby cosmic filament
structures?

To this end, we fetch the spectroscopic and photometric properties of $n=1185$
observed galaxies within a thin redshift slice $0.4\sim 0.4005$ from SDSS-IV
DR16 database as basic covariates. Such a thin redshift slice helps control
the redshift distortions caused by the Kaiser effect (_i.e._ , galaxies
falling towards the galaxy cluster; Kaiser 1987) and the “Finger-Of-God”
effects (_i.e._ , the elongation of galaxy distributions along the light-of-
sight direction; Jackson 1972). Among the galaxies in the selected redshift
slice, 30.2% of their stellar masses are missing in the Firefly catalog under
the “Chabrier” initial stellar mass function (Chabrier, 2003) and the “MILES”
stellar population model (Maraston and Strömbäck, 2011). In addition, the
covariate vector for each galaxy consists of right ascension, declination,
ugriz bands and their measurement errors, model and cmodel magnitudes for
ugriz bandpasses, galactic extinction corrections for ugriz bandpass, median
signal-to-noise per pixel within each ugriz bandpass, original and best-fit
template spectra projected onto ugriz filters as well as their inverse
variances, sky flux in each of the ugriz imaging filters, and signal-to-noise
squared for spectrograph #1 and #2, at g=20.20, r=20.20, i=20.20 for SDSS
spectrograph spectra. To capture the nonlinear patterns of stellar mass
estimation and handle those extreme values, we generate additional covariates
by applying the logarithmic transformations
$x\mapsto\text{sign}(x)\cdot\log(1+|x|)$ to ugriz bands, galactic extinction
corrections, and those flux features. We also remove those covariates whose
correlation coefficients are higher than 0.95 before generating univariate
B-spline base covariates of polynomial order 3 with 40 knots. The reason why
we adopt the B-spline base covariates instead of the usual polynomial
combinations is that the covariates generated by univariate B-spline bases
tends to be less linearly correlated and can facilitate the subsequent
inference task. To address _Question 2_ above, we compute the angular diameter
distances from the galaxies to the two-dimensional spherical cosmic filaments
constructed by the directional subspace constrained mean shift algorithm in
Zhang and Chen (2023); Zhang et al. (2022) on SDSS-IV DR16 data.333This cosmic
web catalog can be downloaded at https://doi.org/10.5281/zenodo.6244866. We
also control the confounding effects of nearby galaxy clusters by including
the angular diameter distances from galaxies to the estimated local modes and
intersections of filaments as extra covariates. The final dimension of the
covariate vector for each galaxy in the selected redshift slice is $d=1409$,
and we take the usual logarithm of the Firefly stellar mass of each galaxy as
the outcome variable.

###  5.2 Results

Figure 8: Stellar mass inference with our debiasing and related methods. Left
Panel: Estimated propensity scores by “LR”, “NN”, and “NNcal”. Middle Panel:
95% confidence intervals by different methods for the estimated stellar mass
of the new galaxy, addressing _Question 1_. Right Panel: 95% confidence
intervals by different methods for the estimated regression coefficient
associated with the distance to cosmic filaments, addressing _Question 2_.

When approaching _Question 1_ , we randomly select a galaxy in the nearby
redshift slice $0.4005\sim 0.401$ as a new observation and construct its
associated covariate vector as the query point $x^{(Q1)}$ following the
procedures in Section 5.1. As for _Question 2_ , we set the query point as
$x^{(Q2)}=(0,0,1,0,...,0)^{T}\in\mathbb{R}^{d}$, where the only nonzero entry
corresponds to the covariate of the angular diameter distance to cosmic
filaments. We estimate the propensity scores (_i.e._ , the probabilities of
having a non-missing stellar masses) on the galaxies in the redshift slice by
the Lasso-type logistic regression (“LR”) in Section C and the neural network
(“NN”) model described in Section 4.4. Nevertheless, as shown by the left
panel of Figure 8 that the estimated propensity scores by the neural network
model are too extreme, _i.e._ , close to 0 or 1, we also consider calibrating
them through the Platt’s logistic model (Platt, 1999), which is denoted by
“NNcal”.

Figure 8 displays the inference results by our proposed debiasing method and
the classical approaches in Section 4.1 applied to the complete-case data. Due
to the sparsity structure of the design matrix created by B-spline bases, “DL-
Jav” and “DL-vdG” fail to produce their inference results. For _Question 1_ ,
the 95% confidence intervals obtained from our debiasing method cover the
stellar mass of the new galaxy in the Firefly catalog no matter how we
estimate the propensity scores or select the tuning parameter, while the Lasso
refitting method underestimates the stellar mass in a statistically
significant level. For _Question 2_ , the 95% confidence intervals produced by
our debiasing and Lasso refitting methods for the regression coefficient of
the distance to cosmic filaments are both below 0, aligning with the existing
findings that galaxies are less massive when they are farther away from the
cosmic filaments (Alpaslan et al., 2016; Kraljic et al., 2018; Malavasi et
al., 2022). On the contrary, the 95% confidence interval obtained from the
ridge projection method contains zero and cannot draw a statistically
significant conclusion about the relation between stellar masses of galaxies
and their distances to cosmic filaments. These results consolidate the
usefulness of our proposed debiasing method in practice.

##  6 Discussion

This paper proposes a novel debiasing method for conducting valid inference on
high-dimensional linear models with MAR outcomes. We establish the asymptotic
normality of the debiasing method and illuminate its statistical and
computational efficiencies through the dual formulation of the key debiasing
program. Simulation studies and real-world applications of our proposed method
demonstrate its superiority over the existing high-dimensional inference
methods in the literature. There are several potential applications and
extensions that can further advance the impacts of our debiasing method.

1\. Applications in causal inference: As discussed by Ding and Li (2018);
Chakrabortty et al. (2019), the missing outcome setting in this paper
incorporates the classical setup in causal inference under the framework of
potential outcome (Rubin, 1974) as a special case. Specifically, the
observable data in causal inference problems are of the form
$\left(\mathbb{Y},T,X\right)\in\mathbb{R}\times\\{0,1\\}\times\mathbb{R}^{d}$,
where $T\in\\{0,1\\}$ denotes a binary treatment assignment indicator and
$\mathbb{Y}=T\cdot Y(1)+(1-T)\cdot Y(0)$ with $Y(0),Y(1)$ being the potential
outcome of $Y$ within the control group $T=0$ and the treatment group $T=1$
respectively. Given that at most one of the potential outcome $Y(0),Y(1)$ for
each subject is observed, the potential outcome framework can be related to
our problem setting by taking $(Y,R)=(Y(1),T)$ in the treatment group or
$(Y,R)=(Y(0),1-T)$ in the control group. Therefore, our debiasing inference
method can be applied to high-dimensional causal inference problems.

First, our debasing method can be used to conduct valid statistical inference
on the regression function (or the conditional mean outcome) of the treatment
group if we focus on the observational data
$\\{(Y_{i},R_{i},X_{i})\\}_{i=1}^{n}=\\{(Y(1)_{i},T_{i},X_{i})\\}_{i=1}^{n}$
or the control group if the data of interest is
$\\{(Y_{i},R_{i},X_{i})\\}_{i=1}^{n}=\\{(Y(0)_{i},T_{i},X_{i})\\}_{i=1}^{n}$.
Unlike the classical approaches in causal inference, where the observational
data $\\{(\mathbb{Y}_{i},T_{i},X_{i})\\}_{i=1}^{n}$ are split into treatment
and control halves in order to conduct statistical inference on the estimand
of interest, we utilize all the covariate vectors
$X_{i}\in\mathbb{R}^{d},i=1,...,n$ in both treatment and control groups to
address the above inference problems, providing extra efficiency gains.

Second, under Assumption 1, our proposed debiasing program (3) can also be
extended to conduct statistical inference on the linear average conditional
treatment effect (ACTE) $\mathrm{E}[Y(1)-Y(0)|X]$ with no unmeasured
confounding. To this end, borrowing an idea from Giessing and Wang (2021), the
debiasing program (3) can be modified as follows:

$\displaystyle\operatorname*{arg\,min}_{\bm{w}_{(0)},\bm{w}_{(1)}\in\mathbb{R}^{n}}\sum_{i=1}^{n}\left[\widehat{\pi}_{i}w_{i(1)}^{2}+(1-\widehat{\pi}_{i})w_{i(0)}^{2}\right]$
$\displaystyle\text{ subject to
}\left|\left|x-\frac{1}{\sqrt{n}}\sum_{i=1}^{n}w_{i(1)}\cdot\widehat{\pi}_{i}\cdot
X_{i}\right|\right|_{\infty}\leq\frac{\gamma_{1}}{n}\quad\text{ and
}\quad\left|\left|x-\frac{1}{\sqrt{n}}\sum_{i=1}^{n}w_{i(0)}\left(1-\widehat{\pi}_{i}\right)X_{i}\right|\right|_{\infty}\leq\frac{\gamma_{2}}{n},$

where $\gamma_{1},\gamma_{2}>0$ are the tuning parameters and
$\widehat{\pi}_{i},i=1,...,n$ are the estimated treatment assignment
probabilities (_i.e._ , propensity scores). The debiased estimator (4) will
become

$\displaystyle\widehat{m}^{\text{debias}}(x;\widehat{\bm{w}}_{(1)},\widehat{\bm{w}}_{(0)})$
$\displaystyle=x^{T}\left(\widehat{\beta}_{(1)}-\widehat{\beta}_{(0)}\right)+\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[\widehat{w}_{i(1)}\cdot
T_{i}\left(Y(1)_{i}-X_{i}^{T}\widehat{\beta}_{(1)}\right)-\widehat{w}_{i(0)}\cdot\left(1-T_{i}\right)\left(Y(0)_{i}-X_{i}^{T}\widehat{\beta}_{(0)}\right)\right].$

However, this new estimator may no longer embrace the asymptotically semi-
parametric efficiency, and we will leave the thorough theoretical
investigation to the future work. Other related works for inferring the linear
CATE with high-dimensional covariates can be found in Chernozhukov and
Semenova (2018).

2\. Model mis-specification and robustness: Theoretical guarantees of our
debiasing method are validated under the linear model (1) and sub-Gaussianity
assumptions. While we have conducted additional simulations in Section 4.3 to
certify the robustness of our debiasing method against the violations of the
sub-Gaussian noise assumption, it is of research interest to investigate the
performances and theoretical implications of our debiasing method when the
assumed linear model is misspecified. Given that our analysis of the inference
with random design is unconditional on the covariate vector
$X\in\mathbb{R}^{d}$, the modification technique in Bühlmann and van de Geer
(2015) might provide some insights.

3\. Missing covariates: Extending our debiasing method to address issues
related to incomplete covariates remains an open problem. Notably, the problem
is more challenging under high-dimensional settings because there are numerous
covariates, resulting in a large number of potential missing patterns (up to
$2^{d}$ at most). While there have been some efforts to leverage the MCAR
assumption (Wang et al., 2019) to tackle this problem, it is noteworthy that
MCAR is a strong assumption that may not hold in many real-world scenarios. A
possible remedy is to use the multiple imputation method (Carpenter et al.,
2023), but it is unclear how to design a reliable imputation model.

## Acknowledgement

We thank Armeen Taeb for helpful discussions about the computational part of
this paper. YZ is supported in part by YC’s NSF grant DMS-2141808. AG is
supported by NSF grant DMS-2310578. YC is supported by NSF grants DMS-1952781,
2112907, 2141808, and NIH U24-AG07212.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P.
Sloan Foundation, the U.S. Department of Energy Office of Science, and the
Participating Institutions. SDSS-IV acknowledges support and resources from
the Center for High Performance Computing at the University of Utah. The SDSS
website is www.sdss.org.

SDSS-IV is managed by the Astrophysical Research Consortium for the
Participating Institutions of the SDSS Collaboration including the Brazilian
Participation Group, the Carnegie Institution for Science, Carnegie Mellon
University, Center for Astrophysics — Harvard & Smithsonian, the Chilean
Participation Group, the French Participation Group, Instituto de Astrofísica
de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and
Mathematics of the Universe (IPMU) / University of Tokyo, the Korean
Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut
für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA
Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-
Institut für Extraterrestrische Physik (MPE), National Astronomical
Observatories of China, New Mexico State University, New York University,
University of Notre Dame, Observatório Nacional / MCTI, The Ohio State
University, Pennsylvania State University, Shanghai Astronomical Observatory,
United Kingdom Participation Group, Universidad Nacional Autónoma de México,
University of Arizona, University of Colorado Boulder, University of Oxford,
University of Portsmouth, University of Utah, University of Virginia,
University of Washington, University of Wisconsin, Vanderbilt University, and
Yale University.

## References

  * Agrawal et al. (2018) A. Agrawal, R. Verschueren, S. Diamond, and S. Boyd.  A rewriting system for convex optimization problems.  _Journal of Control and Decision_ , 5(1):42–60, 2018.
  * Ahumada et al. (2020) R. Ahumada, C. A. Prieto, A. Almeida, F. Anders, S. F. Anderson, B. H. Andrews, B. Anguiano, R. Arcodia, E. Armengaud, M. Aubert, et al.  The 16th data release of the sloan digital sky surveys: first release from the apogee-2 southern survey and full release of eboss spectra.  _The Astrophysical Journal Supplement Series_ , 249(1):3, 2020.
  * Alpaslan et al. (2016) M. Alpaslan, M. Grootes, P. M. Marcum, C. Popescu, R. Tuffs, J. Bland-Hawthorn, S. Brough, M. J. Brown, L. J. Davies, S. P. Driver, et al.  Galaxy and mass assembly (gama): stellar mass growth of spiral galaxies in the cosmic web.  _Monthly Notices of the Royal Astronomical Society_ , 457(3):2287–2300, 2016.
  * Baek et al. (2009) S. Baek, C.-A. Tsai, and J. J. Chen.  Development of biomarker classifiers from high-dimensional data.  _Briefings in bioinformatics_ , 10(5):537–546, 2009.
  * Battey and Reid (2023) H. S. Battey and N. Reid.  On inference in high-dimensional regression.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 85(1):149–175, 2023.
  * Bayati et al. (2013) M. Bayati, M. A. Erdogdu, and A. Montanari.  Estimating lasso risk and noise level.  _Advances in Neural Information Processing Systems_ , 26, 2013.
  * Belloni and Chernozhukov (2013) A. Belloni and V. Chernozhukov.  Least squares after model selection in high-dimensional sparse models.  _Bernoulli_ , 19(2):521–547, 2013.
  * Belloni et al. (2016) A. Belloni, V. Chernozhukov, and Y. Wei.  Post-selection inference for generalized linear models with many controls.  _Journal of Business & Economic Statistics_, 34(4):606–619, 2016.
  * Belloni et al. (2019) A. Belloni, V. Chernozhukov, and K. Kato.  Valid post-selection inference in high-dimensional approximately sparse quantile regression models.  _Journal of the American Statistical Association_ , 114(526):749–758, 2019.
  * Bickel et al. (2009) P. J. Bickel, Y. Ritov, and A. B. Tsybakov.  Simultaneous analysis of lasso and dantzig selector.  _The Annals of statistics_ , 37(4):1705–1732, 2009.
  * Blanton et al. (2005) M. R. Blanton, D. J. Schlegel, M. A. Strauss, J. Brinkmann, D. Finkbeiner, M. Fukugita, J. E. Gunn, D. W. Hogg, Ž. Ivezić, G. Knapp, et al.  New york university value-added galaxy catalog: a galaxy catalog based on new public surveys.  _The Astronomical Journal_ , 129(6):2562, 2005\.
  * Breiman (2001) L. Breiman.  Random forests.  _Machine learning_ , 45:5–32, 2001.
  * Breiman et al. (1984) L. Breiman, J. Friedman, C. J. Stone, and R. Olshen.  _Classification and Regression Trees_.  Chapman and Hall/CRC, 1984.
  * Bühlmann (2013) P. Bühlmann.  Statistical significance in high-dimensional linear models.  _Bernoulli_ , 19(4):1212 – 1242, 2013.
  * Bühlmann and van De Geer (2011) P. Bühlmann and S. van De Geer.  _Statistics for high-dimensional data: methods, theory and applications_.  Springer Science & Business Media, 2011.
  * Bühlmann and van de Geer (2015) P. Bühlmann and S. van de Geer.  High-dimensional inference in misspecified linear models.  _Electronic Journal of Statistics_ , 9(1):1449 – 1473, 2015.
  * Bundy et al. (2014) K. Bundy, M. A. Bershady, D. R. Law, R. Yan, N. Drory, N. MacDonald, D. A. Wake, B. Cherinka, J. R. Sánchez-Gallego, A.-M. Weijmans, et al.  Overview of the sdss-iv manga survey: mapping nearby galaxies at apache point observatory.  _The Astrophysical Journal_ , 798(1):7, 2014.
  * Cai and Guo (2017) T. T. Cai and Z. Guo.  Confidence intervals for high-dimensional linear regression: Minimax rates and adaptivity.  _The Annals of Statistics_ , 45(2):615 – 646, 2017.
  * Cai and Guo (2018) T. T. Cai and Z. Guo.  Accuracy assessment for high-dimensional linear regression.  _The Annals of Statistics_ , 46(4):1807–1836, 2018.
  * Cai et al. (2016) T. T. Cai, W. Liu, and H. H. Zhou.  Estimating sparse precision matrix: Optimal rates of convergence and adaptive estimation.  _The Annals of Statistics_ , 44(2):455–488, 2016\.
  * Cai et al. (2021) T. T. Cai, Z. Guo, and R. Ma.  Statistical inference for high-dimensional generalized linear models with binary outcomes.  _Journal of the American Statistical Association_ , pages 1–14, 2021\.
  * Carpenter et al. (2023) J. R. Carpenter, J. Bartlett, T. Morris, A. Wood, M. Quartagno, and M. G. Kenward.  _Multiple imputation and its application_.  John Wiley & Sons, 2023.
  * Celentano and Wainwright (2023) M. Celentano and M. J. Wainwright.  Challenges of the inconsistency regime: Novel debiasing methods for missing data models.  _arXiv preprint arXiv:2309.01362_ , 2023.
  * Chabrier (2003) G. Chabrier.  Galactic stellar and substellar initial mass function1.  _Publications of the Astronomical Society of the Pacific_ , 115(809):763, 2003.
  * Chakrabortty et al. (2019) A. Chakrabortty, J. Lu, T. T. Cai, and H. Li.  High dimensional m-estimation with missing outcomes: A semi-parametric framework.  _arXiv preprint arXiv:1911.11345_ , 2019.
  * Chang et al. (2015) Y.-Y. Chang, A. van der Wel, E. da Cunha, and H.-W. Rix.  Stellar masses and star formation rates for 1 m galaxies from sdss+ wise.  _The Astrophysical Journal Supplement Series_ , 219(1):8, 2015.
  * Chapelle et al. (2006) O. Chapelle, B. Schölkopf, and A. Zien.  _Semi-Supervised Learning_.  The MIT Press, 2006.
  * Chen et al. (2005) P.-H. Chen, C.-J. Lin, and B. Schölkopf.  A tutorial on $\nu$-support vector machines.  _Applied Stochastic Models in Business and Industry_ , 21(2):111–136, 2005.
  * Chen and Yang (2021) Y. Chen and Y. Yang.  The one standard error rule for model selection: Does it work?  _Stats_ , 4(4):868–892, 2021.
  * Chernozhukov and Semenova (2018) V. Chernozhukov and V. Semenova.  Simultaneous inference for Best Linear Predictor of the Conditional Average Treatment Effect and other structural functions.  CeMMAP working papers CWP40/18, Centre for Microdata Methods and Practice, Institute for Fiscal Studies, July 2018.
  * Chernozhukov et al. (2018) V. Chernozhukov, D. Chetverikov, M. Demirer, E. Duflo, C. Hansen, W. Newey, and J. Robins.  Double/debiased machine learning for treatment and structural parameters.  _The Econometrics Journal_ , 21(1):C1–C68, 01 2018.
  * Comparat et al. (2017) J. Comparat, C. Maraston, D. Goddard, V. Gonzalez-Perez, J. Lian, S. Meneses-Goytia, D. Thomas, J. R. Brownstein, R. Tojeiro, A. Finoguenov, A. Merloni, F. Prada, M. Salvato, G. B. Zhu, H. Zou, and J. Brinkmann.  Stellar population properties for 2 million galaxies from SDSS DR14 and DEEP2 DR4 from full spectral fitting.  _arXiv preprint arXiv:1711.06575_ , 2017.
  * Conroy et al. (2009) C. Conroy, J. E. Gunn, and M. White.  The propagation of uncertainties in stellar population synthesis modeling. i. the relevance of uncertain aspects of stellar evolution and the initial mass function to the derived physical properties of galaxies.  _The Astrophysical Journal_ , 699(1):486, 2009\.
  * Dawson et al. (2016) K. S. Dawson, J.-P. Kneib, W. J. Percival, S. Alam, F. D. Albareti, S. F. Anderson, E. Armengaud, É. Aubourg, S. Bailey, J. E. Bautista, et al.  The sdss-iv extended baryon oscillation spectroscopic survey: overview and early data.  _The Astronomical Journal_ , 151(2):44, 2016.
  * Dezeure et al. (2015) R. Dezeure, P. Bühlmann, L. Meier, and N. Meinshausen.  High-dimensional inference: confidence intervals, p-values and r-software hdi.  _Statistical science_ , pages 533–558, 2015.
  * Diamond and Boyd (2016) S. Diamond and S. Boyd.  CVXPY: A Python-embedded modeling language for convex optimization.  _Journal of Machine Learning Research_ , 17(83):1–5, 2016.
  * Dicker (2012) L. H. Dicker.  Residual variance and the signal-to-noise ratio in high-dimensional linear models.  _arXiv preprint arXiv:1209.0012_ , 2012.
  * Ding and Li (2018) P. Ding and F. Li.  Causal Inference: A Missing Data Perspective.  _Statistical Science_ , 33(2):214 – 237, 2018\.
  * D’Amour et al. (2021) A. D’Amour, P. Ding, A. Feller, L. Lei, and J. Sekhon.  Overlap in observational studies with high-dimensional covariates.  _Journal of Econometrics_ , 221(2):644–654, 2021\.
  * Fan et al. (2011) J. Fan, J. Lv, and L. Qi.  Sparse high-dimensional models in economics.  _Annual Review of Economics_ , 3(1):291–317, 2011\.
  * Fan et al. (2012) J. Fan, S. Guo, and N. Hao.  Variance estimation using refitted cross-validation in ultrahigh dimensional regression.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 74(1):37–65, 2012.
  * Farrell et al. (2021) M. H. Farrell, T. Liang, and S. Misra.  Deep neural networks for estimation and inference.  _Econometrica_ , 89(1):181–213, 2021.
  * Feller (1968) W. Feller.  _An introduction to probability theory and its application_.  John Wiley & Sons Inc., New York, 3rd edition, 1968.
  * Foucart and Rauhut (2013) S. Foucart and H. Rauhut.  _A Mathematical Introduction to Compressive Sensing_.  Applied and Numerical Harmonic Analysis. Springer Birkhäuser Basel, 2013.
  * Friedman et al. (2008) J. Friedman, T. Hastie, and R. Tibshirani.  Sparse inverse covariance estimation with the graphical lasso.  _Biostatistics_ , 9(3):432–441, 2008.
  * Frumento et al. (2012) P. Frumento, F. Mealli, B. Pacini, and D. B. Rubin.  Evaluating the effect of training on wages in the presence of noncompliance, nonemployment, and missing outcome data.  _Journal of the American Statistical Association_ , 107(498):450–466, 2012.
  * Fu et al. (2020) A. Fu, B. Narasimhan, and S. Boyd.  CVXR: An R package for disciplined convex optimization.  _Journal of Statistical Software_ , 94(14):1–34, 2020.  doi: 10.18637/jss.v094.i14.
  * Gao et al. (2022) W. Gao, F. Xu, and Z.-H. Zhou.  Towards convergence rate analysis of random forests for classification.  _Artificial Intelligence_ , 313:103788, 2022.
  * Geng and Huang (2000) Z. Geng and L. Huang.  Robust stability of the systems with mixed uncertainties under the iqc descriptions.  _International Journal of Control_ , 73(9):776–786, 2000.
  * Giessing (2022) A. Giessing.  Inequalities for suprema of unbounded empirical processes.  _Technical Report_ , 2022.
  * Giessing and Wang (2021) A. Giessing and J. Wang.  Debiased inference on heterogeneous quantile treatment effects with regression rank-scores.  _arXiv preprint arXiv:2102.01753_ , 2021.  URL https://arxiv.org/abs/2102.01753.
  * Giné and Nickl (2021) E. Giné and R. Nickl.  _Mathematical foundations of infinite-dimensional statistical models_.  Cambridge university press, 2021.
  * Gorski et al. (2007) J. Gorski, F. Pfeuffer, and K. Klamroth.  Biconvex sets and optimization with biconvex functions: a survey and extensions.  _Mathematical methods of operations research_ , 66(3):373–407, 2007.
  * Graham (2011) B. S. Graham.  Efficiency bounds for missing data models with semiparametric restrictions.  _Econometrica_ , 79(2):437–452, 2011.
  * Guo and Chen (2016) B. Guo and S. X. Chen.  Tests for high dimensional generalized linear models.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 78(5):1079–1102, 2016.
  * Guo et al. (2021) Z. Guo, P. Rakshit, D. S. Herman, and J. Chen.  Inference for the case probability in high-dimensional logistic regression.  _The Journal of Machine Learning Research_ , 22(1):11480–11533, 2021.
  * Henmi and Eguchi (2004) M. Henmi and S. Eguchi.  A paradox concerning nuisance parameters and projected estimating functions.  _Biometrika_ , 91(4):929–941, 2004.
  * Higgins et al. (2008) J. P. Higgins, I. R. White, and A. M. Wood.  Imputation methods for missing outcome data in meta-analysis of clinical trials.  _Clinical trials_ , 5(3):225–239, 2008.
  * Hinton (1990) G. E. Hinton.  Connectionist learning procedures.  In _Machine learning_ , pages 555–610. Elsevier, 1990.
  * Hitomi et al. (2008) K. Hitomi, Y. Nishiyama, and R. Okui.  A puzzling phenomenon in semiparametric estimation problems with infinite-dimensional nuisance parameters.  _Econometric Theory_ , 24(6):1717–1728, 2008\.
  * Hoerl and Kennard (1970) A. E. Hoerl and R. W. Kennard.  Ridge regression: applications to nonorthogonal problems.  _Technometrics_ , 12(1):69–82, 1970.
  * Huang and Knowles (2016) H. Huang and L. L. Knowles.  Unforeseen Consequences of Excluding Missing Data from Next-Generation Sequences: Simulation Study of RAD Sequences.  _Systematic Biology_ , 65(3):357–365, 07 2016\.
  * Huang and Zhang (2012) J. Huang and C.-H. Zhang.  Estimation and selection via absolute penalized convex minimization and its multistage adaptive applications.  _The Journal of Machine Learning Research_ , 13(1):1839–1864, 2012.
  * Imbens (2004) G. W. Imbens.  Nonparametric estimation of average treatment effects under exogeneity: A review.  _Review of Economics and Statistics_ , 86(1):4–29, 2004.
  * Jackson (1972) J. Jackson.  A critique of rees’s theory of primordial gravitational radiation.  _Monthly Notices of the Royal Astronomical Society_ , 156(1):1P–5P, 1972.
  * Jankova and van de Geer (2018) J. Jankova and S. van de Geer.  Semiparametric efficiency bounds for high-dimensional models.  _The Annals of Statistics_ , 46(5):2336–2359, 2018.
  * Javanmard and Lee (2020) A. Javanmard and J. D. Lee.  A flexible framework for hypothesis testing in high dimensions.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 82(3):685–718, 2020.
  * Javanmard and Montanari (2014a) A. Javanmard and A. Montanari.  Confidence intervals and hypothesis testing for high-dimensional regression.  _The Journal of Machine Learning Research_ , 15(1):2869–2909, 2014a.
  * Javanmard and Montanari (2014b) A. Javanmard and A. Montanari.  Hypothesis testing in high-dimensional regression under the gaussian random design model: Asymptotic theory.  _IEEE Transactions on Information Theory_ , 60(10):6522–6554, 2014b.
  * Javanmard and Montanari (2018) A. Javanmard and A. Montanari.  Debiasing the lasso: Optimal sample size for gaussian designs.  _The Annals of Statistics_ , 46(6A):2593–2622, 2018.
  * Kaiser (1987) N. Kaiser.  Clustering in real space and in redshift space.  _Monthly Notices of the Royal Astronomical Society_ , 227(1):1–21, 1987.
  * Koltchinskii (2011) V. Koltchinskii.  _Oracle inequalities in empirical risk minimization and sparse recovery problems: École D’Été de Probabilités de Saint-Flour XXXVIII-2008_ , volume 2033.  Springer Science & Business Media, 2011.
  * Kraljic et al. (2018) K. Kraljic, S. Arnouts, C. Pichon, C. Laigle, S. de la Torre, D. Vibert, C. Cadiou, Y. Dubois, M. Treyer, C. Schimd, et al.  Galaxy evolution in the metric of the cosmic web.  _Monthly Notices of the Royal Astronomical Society_ , 474(1):547–571, 2018.
  * Li (2020) S. Li.  Debiasing the debiased Lasso with bootstrap.  _Electronic Journal of Statistics_ , 14(1):2298 – 2337, 2020.
  * Lindsey (2004) J. K. Lindsey.  _Statistical analysis of stochastic processes in time_ , volume 14.  Cambridge University Press, 2004.
  * Little and Rubin (2019) R. J. Little and D. B. Rubin.  _Statistical analysis with missing data_.  John Wiley & Sons, third edition, 2019.
  * Loh and Wainwright (2012) P.-L. Loh and M. J. Wainwright.  High-dimensional regression with noisy and missing data: Provable guarantees with nonconvexity.  _The Annals of Statistics_ , 40(3):1637 – 1664, 2012.
  * Lok (2021) J. J. Lok.  How estimating nuisance parameters can reduce the variance (with consistent variance estimation).  _arXiv preprint arXiv:2109.02690_ , 2021.
  * Ma et al. (2021) R. Ma, T. Tony Cai, and H. Li.  Global and simultaneous hypothesis testing for high-dimensional logistic regression models.  _Journal of the American Statistical Association_ , 116(534):984–998, 2021.
  * Majewski et al. (2017) S. R. Majewski, R. P. Schiavon, P. M. Frinchaboy, C. A. Prieto, R. Barkhouser, D. Bizyaev, B. Blank, S. Brunner, A. Burton, R. Carrera, et al.  The apache point observatory galactic evolution experiment (apogee).  _The Astronomical Journal_ , 154(3):94, 2017.
  * Malavasi et al. (2022) N. Malavasi, M. Langer, N. Aghanim, D. Galárraga-Espinosa, and C. Gouin.  Relative effect of nodes and filaments of the cosmic web on the quenching of galaxies and the orientation of their spin.  _Astronomy & Astrophysics_, 658:A113, 2022.
  * Maraston and Strömbäck (2011) C. Maraston and G. Strömbäck.  Stellar population models at high spectral resolution.  _Monthly Notices of the Royal Astronomical Society_ , 418(4):2785–2811, 2011.
  * Müller and van de Geer (2015) P. Müller and S. van de Geer.  The partial linear model in high dimensions.  _Scandinavian Journal of Statistics_ , 42(2):580–608, 2015.
  * Müller and Keilegom (2012) U. U. Müller and I. V. Keilegom.  Efficient parameter estimation in regression with missing responses.  _Electronic Journal of Statistics_ , 6(none):1200 – 1219, 2012.
  * Negahban et al. (2010) S. N. Negahban, P. Ravikumar, M. J. Wainwright, and B. Yu.  A unified framework for high-dimensional analysis of $m$-estimators with decomposable regularizers.  _arXiv preprint arXiv:1010.2731v1_ , 2010.  URL https://arxiv.org/abs/1010.2731v1.
  * Neyman (1959) J. Neyman.  Optimal asymptotic tests of composite hypotheses.  _Probability and statsitics_ , pages 213–234, 1959.
  * Neyman (1979) J. Neyman.  C($\alpha$) tests and their use.  _Sankhyā: The Indian Journal of Statistics, Series A_ , pages 1–21, 1979.
  * Nickl and van de Geer (2013) R. Nickl and S. van de Geer.  Confidence sets in sparse regression.  _The Annals of Statistics_ , 41(6):2852 – 2876, 2013.
  * Ning and Liu (2017) Y. Ning and H. Liu.  A general theory of hypothesis tests and confidence regions for sparse high dimensional models.  _The Annals of Statistics_ , 45(1):158–195, 2017\.
  * Pedregosa et al. (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011\.
  * Platt (1999) J. Platt.  Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods.  _Advances in large margin classifiers_ , 10(3):61–74, 1999.
  * Ravikumar et al. (2009) P. Ravikumar, J. Lafferty, H. Liu, and L. Wasserman.  Sparse additive models.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 71(5):1009–1030, 2009.
  * Reid et al. (2016) S. Reid, R. Tibshirani, and J. Friedman.  A study of error variance estimation in lasso regression.  _Statistica Sinica_ , pages 35–67, 2016.
  * Renegar (2001) J. Renegar.  _A mathematical view of interior-point methods in convex optimization_.  MOS-SIAM Series on Optimization. Society for Industrial and Applied Mathematics, 2001.
  * Robins and Rotnitzky (1995) J. M. Robins and A. Rotnitzky.  Semiparametric efficiency in multivariate regression models with missing data.  _Journal of the American Statistical Association_ , 90(429):122–129, 1995.
  * Robins et al. (1992) J. M. Robins, S. D. Mark, and W. K. Newey.  Estimating exposure effects by modelling the expectation of exposure conditional on confounders.  _Biometrics_ , pages 479–495, 1992.
  * Robins et al. (1994) J. M. Robins, A. Rotnitzky, and L. P. Zhao.  Estimation of regression coefficients when some regressors are not always observed.  _Journal of the American statistical Association_ , 89(427):846–866, 1994.
  * Robins et al. (1995) J. M. Robins, A. Rotnitzky, and L. P. Zhao.  Analysis of semiparametric regression models for repeated outcomes in the presence of missing data.  _Journal of the american statistical association_ , 90(429):106–121, 1995.
  * Rosenbaum (1987) P. R. Rosenbaum.  Model-based direct adjustment.  _Journal of the American statistical Association_ , 82(398):387–394, 1987.
  * Rosenbaum and Rubin (1983) P. R. Rosenbaum and D. B. Rubin.  The central role of the propensity score in observational studies for causal effects.  _Biometrika_ , 70(1):41–55, 1983.
  * Rosenblatt (2015) M. Rosenblatt.  Short range and long range dependence.  _Asymptotic Laws and Methods in Stochastics: A Volume in Honour of Miklós Csörgő_ , pages 283–294, 2015.
  * Rubin (1974) D. B. Rubin.  Estimating causal effects of treatments in randomized and nonrandomized studies.  _Journal of educational Psychology_ , 66(5):688, 1974.
  * Shi et al. (2021) C. Shi, R. Song, W. Lu, and R. Li.  Statistical inference for high-dimensional models via recursive online-score estimation.  _Journal of the American Statistical Association_ , 116(535):1307–1318, 2021.
  * Steinwart (2001) I. Steinwart.  On the influence of the kernel on the consistency of support vector machines.  _Journal of machine learning research_ , 2(Nov):67–93, 2001.
  * Su et al. (2023) F. Su, W. Mou, P. Ding, and M. Wainwright.  When is the estimated propensity score better? high-dimensional analysis and bias correction.  _arXiv preprint arXiv:2303.17102_ , 2023.
  * Sun and Zhang (2012) T. Sun and C.-H. Zhang.  Scaled sparse linear regression.  _Biometrika_ , 99(4):879–898, 2012.
  * Sun and Zhang (2013) T. Sun and C.-H. Zhang.  Sparse matrix inversion with scaled lasso.  _The Journal of Machine Learning Research_ , 14(1):3385–3418, 2013.
  * Tibshirani (1996) R. Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 58(1):267–288, 1996.
  * Tsiatis (2007) A. A. Tsiatis.  _Semiparametric Theory and Missing Data_.  Springer Series in Statistics. Springer New York, 2007.
  * van de Geer et al. (2014) S. van de Geer, P. Bühlmann, Y. Ritov, and R. Dezeure.  On asymptotically optimal confidence regions and tests for high-dimensional models.  _The Annals of Statistics_ , 42(3):1166–1202, 2014.
  * van de Geer (2008) S. A. van de Geer.  High-dimensional generalized linear models and the lasso.  _The Annals of Statistics_ , 36(2):614–645, 2008\.
  * van de Geer and Bühlmann (2009) S. A. van de Geer and P. Bühlmann.  On the conditions used to prove oracle results for the lasso.  _Electronic Journal of Statistics_ , 3:1360–1392, 2009\.
  * van der Vaart (2000) A. W. van der Vaart.  _Asymptotic Statistics_.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2000.
  * van der Vaart and Wellner (1996) A. W. van der Vaart and J. A. Wellner.  _Weak Convergence and Empirical Processes: With Applications to Statistics_.  Springer Series in Statistics. Springer, 1996.
  * Vershynin (2018) R. Vershynin.  _High-dimensional Probability: An introduction with applications in data science_ , volume 47.  Cambridge university press, 2018.
  * Wainwright (2019) M. J. Wainwright.  _High-dimensional statistics: A non-asymptotic viewpoint_ , volume 48.  Cambridge University Press, 2019.
  * Wang et al. (2019) Y. Wang, J. Wang, S. Balakrishnan, and A. Singh.  Rate optimal estimation and confidence intervals for high-dimensional regression with missing covariates.  _Journal of Multivariate Analysis_ , 174:104526, 2019.
  * Wilkinson et al. (2017) D. M. Wilkinson, C. Maraston, D. Goddard, D. Thomas, and T. Parikh.  Firefly (fitting iteratively for likelihood analysis): a full spectral fitting code.  _Monthly Notices of the Royal Astronomical Society_ , 472(4):4297–4326, 2017.
  * Wright (2015) S. J. Wright.  Coordinate descent algorithms.  _Mathematical Programming_ , 151(1):3–34, 2015\.
  * Wyss et al. (2022) R. Wyss, C. Yanover, T. El-Hay, D. Bennett, R. W. Platt, A. R. Zullo, G. Sari, X. Wen, Y. Ye, H. Yuan, et al.  Machine learning for improving high-dimensional proxy confounder adjustment in healthcare database studies: An overview of the current literature.  _Pharmacoepidemiology and Drug Safety_ , 31(9):932–943, 2022.
  * Xia et al. (2020) L. Xia, B. Nan, and Y. Li.  A revisit to de-biased lasso for generalized linear models.  _arXiv preprint arXiv:2006.12778_ , 2020.
  * Zhang (2010) C.-H. Zhang.  Nearly unbiased variable selection under minimax concave penalty.  _The Annals of Statistics_ , 38(2):894 – 942, 2010.
  * Zhang and Zhang (2014) C.-H. Zhang and S. S. Zhang.  Confidence intervals for low dimensional parameters in high dimensional linear models.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 76(1):217–242, 2014.
  * Zhang (2004) H. Zhang.  The optimality of naive bayes.  In V. Barr and Z. Markov, editors, _Proceedings of the Seventeenth International Florida Artificial Intelligence Research Society Conference (FLAIRS 2004)_. AAAI Press, 2004.
  * Zhang and Chen (2023) Y. Zhang and Y.-C. Chen.  Linear convergence of the subspace constrained mean shift algorithm: from euclidean to directional data.  _Information and Inference: A Journal of the IMA_ , 12(1):210–311, 2023.
  * Zhang et al. (2022) Y. Zhang, R. S. de Souza, and Y.-C. Chen.  Sconce: a cosmic web finder for spherical and conic geometries.  _Monthly Notices of the Royal Astronomical Society_ , 517(1):1197–1217, 2022.

Supplementary Materials to “Efficient Inference on High-Dimensional Linear
Models with Missing Outcomes”

Yikun Zhang, Alexander Giessing, and Yen-Chi Chen Department of Statistics,
University of Washington

Contents

[sections] [sections]l1

The supplementary material contains auxiliary technical results, additional
simulation studies, and proofs. We organize its content as follows.

  * •

Section A Implementation of the Proposed Debiasing Method We explain the
implementation steps of our debiasing method in detail.

  * •

Section B Sufficient Conditions for Assumption 6 We provide sufficient
conditions and motivating examples under which Assumption 6 holds.

  * •

Section C Generalized Linear Model with Lasso Regularization for Propensity
Score Estimation We showcase an example of estimating the propensity scores
via the Lasso-type generalized linear regression and derive an explicit
formula for rate of convergence $r_{\pi}$ in Assumption 5.

  * •

Section D Additional Simulation Results We furnish extensive simulation
results that are not covered by Section 4 in the main paper.

  * •

Section E Proofs We provide proofs for all the theorems and technical results
in the main paper.

##  Appendix A Implementation of the Proposed Debiasing Method

Our proposed debiasing method consists of three main estimation steps: (i)
computing the Lasso pilot estimate (2), (ii) estimating the propensity scores
$\pi_{i}=\pi(X_{i})=\mathrm{P}(R_{i}=1|X_{i}),i=1,...,n$, and (iii) solving
the debiasing programs (3) and (6). We explain the implementation details and
strategies of selecting two tuning parameters $\lambda>0$ in (2) and
$\gamma>0$ in (3) as follows.

1\. Lasso pilot estimate: We use the complete-case (or equivalently, observed)
data $(X_{i},Y_{i})\in\mathbb{R}^{d}\times\mathbb{R}$ with $R_{i}=1$ for
$i=1,...,n$ to construct the Lasso pilot estimate $\widehat{\beta}$. In order
to select the regularization parameter $\lambda>0$ in a data adaptive way, we
adopt the scaled Lasso (Sun and Zhang, 2012) with its universal regularization
parameter $\lambda_{0}=\sqrt{\frac{2\log d}{n}}$ as the initialization.
Specifically, it provides an iterative algorithm to obtain the Lasso estimate
$\widehat{\beta}$ and a consistent estimator $\widehat{\sigma}_{\epsilon}$ of
the noise level $\sigma_{\epsilon}$ from the following jointly convex
optimization problem:

$\left(\widehat{\beta}(\widetilde{\lambda}),\widehat{\sigma}_{\epsilon}(\widetilde{\lambda})\right)=\operatorname*{arg\,min}_{\beta\in\mathbb{R}^{d},\sigma_{\epsilon}>0}\left[\frac{1}{2n\sigma_{\epsilon}}\sum_{i=1}^{n}R_{i}\left(Y_{i}-X_{i}^{T}\beta\right)^{2}+\frac{\sigma_{\epsilon}}{2}+\widetilde{\lambda}\left|\left|\beta\right|\right|_{1}\right],$

where the regularization parameter $\widetilde{\lambda}>0$ is updated
iteratively as well. We also utilize this estimated noise level
$\widehat{\sigma}_{\epsilon}$ to construct our asymptotic $1-\tau$ confidence
interval (5).

2\. Propensity score estimation: When comparing our debiasing method with
other existing high-dimensional inference methods, we utilize the Lasso-type
logistic regression (27) in Section C to estimate the propensity scores
$\pi_{i}=\mathrm{P}(R_{i}=1|X_{i}),i=1,...,n$, in which the regularization
parameter $\zeta$ is selected on 40 logarithmically spacing values within the
interval $\left[0.1\sqrt{\frac{\log d}{n}},300\sqrt{\frac{\log d}{n}}\right]$
via 5-fold cross-validations. However, as demonstrated in Section 4.4 through
simulation studies, our debiasing method can perform even better when the
propensity scores are estimated more accurately by other nonparametric
methods, reflecting the theoretical results in Section 3.

3\. Solving the debiasing program: When our primal debiasing program is
strictly feasible, strong duality holds by Slater’s condition (Theorem 3.2.6
in Renegar 2001), and the debiasing weights $\widehat{w}_{i}(x),i=1,...,n$ can
be estimated from either the primal program (3) or the dual program (6) by the
identity (7) in Proposition 1. We implement our debiasing method in both
Python and R. In particular, we solve the primal program (3) by the Python
package “CVXPY” (Diamond and Boyd, 2016; Agrawal et al., 2018) or the R
package “CVXR” (Fu et al., 2020). For the dual program (6), we formulate and
implement the classical coordinate descent algorithm (Wright, 2015) as:

$\left[\widehat{\ell}(x)\right]_{j}\leftarrow\frac{\mathcal{S}_{\frac{\gamma}{n}}\left(-\frac{1}{2n}\sum_{i=1}^{n}\widehat{\pi}_{i}\left(\sum_{k\neq
j}X_{ik}X_{jk}\left[\widehat{\ell}(x)\right]_{k}\right)-x_{j}\right)}{\frac{1}{2n}\sum_{i=1}^{n}\widehat{\pi}_{i}X_{ij}^{2}}\quad\text{
for }\quad j=1,...,d,$ (24)

where
$\mathcal{S}_{\frac{\gamma}{n}}(u)=\mathrm{sign}(u)\cdot\left(u-\frac{\gamma}{n}\right)_{+}$
is the soft-thresholding operator, and the above procedure will iterate
cyclically from $j=1$ to $d$ until convergence.

The tuning parameter $\gamma>0$ of our debiasing program is selected from 41
equally spacing values between 0 and $n\left|\left|x\right|\right|_{\infty}$
via 5-fold cross validations according to the dual problem (6), in which
$x\in\mathbb{R}^{d}$ is the current query point. Here, the upper bound for the
candidate values of $\gamma>0$ is set to
$n\left|\left|x\right|\right|_{\infty}$, because the optimal solution
$\widehat{w}(x)$ for the primal problem (3) becomes 0 when
$\frac{\gamma}{n}\geq\left|\left|x\right|\right|_{\infty}$. There are three
different criteria for selecting the final value of $\gamma>0$ from the cross-
validation results.

  1. (i)

The first criterion selects the value of $\gamma>0$ that minimizes the average
cross-validate risk of the dual debiasing program (6), which is known as the
“min-CV” rule.

  2. (ii)

The second criterion picks the value of $\gamma>0$ whose associated average
cross-validated risk is at least one standard error (calculated across the
5-folded data) away from the smallest cross-validated risk and is also smaller
than the one corresponding to the smallest average cross-validated risk; see
Section 3.4.3 in Breiman et al. (1984) and Chen and Yang (2021) for details.
We name this criterion as the “1SE” rule. The rationale behind the “1SE” rule
is to improve the inference performance of our debiasing method by sacrificing
some statistical efficiency (_i.e._ , potentially yielding a larger value of
the estimated asymptotic variance) for a less biased estimator.

  3. (iii)

The third criterion chooses the minimal value of $\gamma>0$ in which the