A Nuisance-Free Inference Procedure Accounting for the Unknown Missingness with Application to Electronic Health Records

Abstract

We study how to conduct statistical inference in a regression model where the outcome variable is prone to missing values and the missingness mechanism is unknown. The model we consider might be a traditional setting or a modern high-dimensional setting where the sparsity assumption is usually imposed and the regularization technique is popularly used. Motivated by the fact that the missingness mechanism, albeit usually treated as a nuisance, is difficult to specify correctly, we adopt the conditional likelihood approach so that the nuisance can be completely ignored throughout our procedure. We establish the asymptotic theory of the proposed estimator and develop an easy-to-implement algorithm via some data manipulation strategy. In particular, under the high-dimensional setting where regularization is needed, we propose a data perturbation method for the post-selection inference. The proposed methodology is especially appealing when the true missingness mechanism tends to be missing not at random, e.g., patient reported outcomes or real world data such as electronic health records. The performance of the proposed method is evaluated by comprehensive simulation experiments as well as a study of the albumin level in the MIMIC-III database.

1. Introduction

A major step towards scientific discovery is to identify useful associations from various features and to quantify their uncertainties. This usually warrants building a regression model for an outcome variable and estimating the coefficient associated with each feature as well as the precision of the estimator. Besides the traditional regression with a small dimensionality, with advances in biotechnology, the modern high-dimensional regression usually posits a sparse parameter in the model, and then applies regularization to select the significant features in order to recover the sparsity. In particular, the post-selection inference could be challenging in a regularized regression framework. In this paper, our main interest is to consider a regression model where the outcome variable is prone to missing values. We study both the traditional setting where regularization is not needed and the modern one with regularization.

The missing data issue is an inevitable concern for statistical analysis in various disciplines ranging from biomedical studies to social sciences. In many applications, the occurrence of missing data is usually not the investigator’s primary interest but complicates the statistical analysis. The validity of any method devised for missing data heavily depends on the assumption of the missingness mechanism [1]. Unfortunately, those assumptions are largely unknown and difficult, if not infeasible, to be empirically tested. Therefore, one prefers to concentrate on analyzing the regression model for the outcome variable, while treating the mechanism model as a nuisance. A flexible assumption imposed at the minimum level on the mechanism would provide protection against model misspecification at this level.

While it is indeed promising to regard the missingness mechanism as a nuisance with a flexible assumption, a potential issue is the model identifiability problem if the mechanism contains missing-not-at-random cases, i.e., allowing the mechanism to depend on the missing values themselves. In the past few years, researchers have made great progress on this topic by introducing a so-called instrument. This instrument could be a shadow variable [2,3,4,5,6,7] or an instrumental variable [8,9]. Both approaches are reasonable and are suitable for different applications. In this paper, we adopt the shadow variable approach as it facilitates the interpretability of the regression model for the outcome. The details of the shadow variable approach will be articulated later throughout the paper.

Therefore, we proceed with a semiparametric framework where our primary interest is a parametric regression, e.g., a linear model, where the statistical task is to estimate the parameter of interest and conduct statistical inference (particularly post-selection inference for the setting with regularization). For the nuisance missingness mechanism, we only impose a nonparametric assumption without specifying a concrete form. We encode the shadow variable as Z, which is one component of the covariate $\mathbf{X}$. In general, a shadow variable with a smaller dimensionality allows more flexibility of the missingness mechanism. Therefore, although it could be multidimensional, we only consider univariate Z throughout the paper. With all of these ingredients, we analyze a conditional likelihood approach which will eventually result in a nuisance-free procedure for parameter estimation and statistical inference.

There are at least two extra highlights of our proposed method that are worth mentioning. The first pertains to the algorithm and computation. Although it looks complicated at first sight, we show that, via some data manipulation strategy, the conditional likelihood function can be analytically written as the likelihood of a conventional logistic regression with some prespecified format. Therefore, our objective function can be readily optimized by many existing software packages. This greatly alleviates the computational burden of our procedure. Second, while the variance estimation under the traditional setting is straightforward following the asymptotic approximation, it is challenging for the setting with regularization. To resolve this problem, we present an easy-to-implement data-driven method to estimate the variance of the regularized estimator via a data perturbation technique. It is noted that the current literature on the inference procedure for regularized estimation in the presence of missing values is very scarce. The authors of [10,11,12] all considered the model selection problem under high dimensionality with missing data; however, none of them studied the post-selection inference in this context.

The remainder of the paper is structured as follows. In Section 2, we first layout our model formulation and introduce the shadow variable and the conditional likelihood. Section 3 details the traditional setting without regularization. We present our algorithm of how to maximize the conditional likelihood function, the theory of how to derive the asymptotic representation of our proposed estimator and how to estimate its variance. In Section 4, we devote ourselves to the modern setting where the sparsity assumption is imposed and the regularization technique is adopted. Both algorithm and theory as well as the variance estimation through the data perturbation technique are presented. In Section 5, we conduct comprehensive simulation studies to examine the finite sample performance of our proposed estimator as well as the comparison to some existing methods. Section 6 is the application of our method to the regression model for the albumin level which suffers from a large amount of missing values in the MIMIC-III study [13]. The paper is concluded with a discussion in Section 7.

2. Methodology

Denote the outcome variable as Y and covariate $\mathbf{X}$. We assume $\mathbf{X}={({\mathbf{U}}^{{}^{\mathrm{T}}},Z)}^{{}^{\mathrm{T}}}$ where $\mathbf{U}$ is p-dimensional and Z univariate, with detailed interpretation later. We consider the linear model

$$\begin{array}{c}\hfill Y=\alpha +{\beta }^{{}^{\mathrm{T}}}\mathbf{U}+\gamma Z+ϵ,\end{array}$$

(1)

where $\beta$ is also p-dimensional, $\alpha$ and $\gamma$ are scalars and the true value of $\gamma$, ${\gamma }_0$, is nonzero, $ϵ\sim N(0,{\sigma }^2)$. We consider the situation that Y has missing values while $\mathbf{X}$ is fully observed. We introduce a binary variable R to indicate missingness: $R=1$ if Y is observed and $R=0$ if missing. To allow the greatest flexibility of the missingness mechanism model, we assume

$$\begin{array}{c}\hfill \mathrm{pr}(R=1|Y,\mathbf{X})=\mathrm{pr}(R=1|Y,\mathbf{U})=s(Y,\mathbf{U}),\end{array}$$

(2)

where $s\left(·\right)$ merely represents an unknown and unspecified function not depending on Z. We reiterate that, as the assumption (2), in a nonparametric flavor, does not specify a concrete form of $s\left(·\right)$, one does not need to be worrisome of the mechanism model misspecification. Moreover, as it allows the dependence on Y, besides missing-completely-at-random (MCAR) and many scenarios of missing-at-random (MAR), the assumption (2) also contains various situations of missing-not-at-random (MNAR).

We term Z the shadow variable following the works in [5,6,7,14]. Its existence depends on whether it is sensible that Z and R are conditionally independent (given Y and $\mathbf{U}$) and that Y heavily relies on Z (as ${\gamma }_0\ne 0$). There are many examples in the literature documenting that the existence of Z is practically reasonable. In application, a surrogate or a proxy of the outcome variable Y, which would not synchronically affect the missingness mechanism, could be a good choice for the shadow variable Z.

We assume independent and identically distributed observations $\{r_i,y_i,{\mathbf{u}}_i,z_i\}$ for $i=1,\dots ,N$ and the first n subjects are free of missing data. Now we present a $s\left(·\right)$-free procedure via the use of the conditional likelihood. Denote $\mathbf{V}={(Y,{\mathbf{U}}^{{}^{\mathrm{T}}})}^{{}^{\mathrm{T}}}$. We start with

$$\begin{array}{c}\hfill \prod _{i=1}^{n}p\left({\mathbf{v}}_i\right|z_i,r_i=1)=\prod _{i=1}^n\frac{s\left({\mathbf{v}}_i\right)}{g\left(z_i\right)}p\left({\mathbf{v}}_i\right|z_i),\end{array}$$

where $g\left(z_i\right)=\mathrm{pr}(r_i=1|z_i)=\int \mathrm{pr}(r_i=1|\mathbf{v})p\left(\mathbf{v}\right|z_i)d\mathbf{v}$ and $p\left(·\right|·)$ is a generic notation for conditional probability density/mass function. If $\mathbf{V}$ were univariate, we denote $\mathcal{A}$ as the rank statistic of $\{v_1,\dots ,v_n\}$, then

$$\begin{array}{ccc}\hfill & & \prod _{i=1}^{n}p\left(v_i\right|z_i,r_i=1)=p(v_1,\dots ,v_n|z_1,\dots ,z_n,r_1=\cdots =r_n=1)\hfill \\ \hfill & =& p\left(\mathcal{A}\right|v_{\left(1\right)},\dots ,v_{\left(n\right)},z_1,\dots ,z_n,r_1=\cdots =r_n=1)p(v_{\left(1\right)},\dots ,v_{\left(n\right)}|z_1,\dots ,z_n,r_1=\cdots =r_n=1).\hfill \end{array}$$

(3)

The conditional likelihood that we use, the first term on the right hand side of (3), is exactly

$$\begin{array}{ccc}& & p\left(\mathcal{A}\right|v_{\left(1\right)},\dots ,v_{\left(n\right)},z_1,\dots ,z_n,r_1=\cdots =r_n=1)=\frac{p(v_1,\dots ,v_n|z_1,\dots ,z_n,r_1=\cdots =r_n=1)}{p(v_{\left(1\right)},\dots ,v_{\left(n\right)}|z_1,\dots ,z_n,r_1=\cdots =r_n=1)}\hfill \\ & =& \frac{{\prod }_{i=1}^{n}p\left(v_i\right|z_i,r_i=1)}{{\Sigma }_{\omega \in \Omega }{\prod }_{i=1}^{n}p\left(v_{\omega \left(i\right)}\right|z_i,r_i=1)}=\frac{{\prod }_{i=1}^{n}p\left(v_i\right|z_i)}{{\Sigma }_{\omega \in \Omega }{\prod }_{i=1}^{n}p\left(v_{\omega \left(i\right)}\right|z_i)},\hfill \end{array}$$

(4)

where $\Omega$ represents the collection of all one-to-one mappings from $\{1,\dots ,n\}$ to $\{1,\dots ,n\}$. Now (4) is nuisance-free and can be used to estimate the unknown parameters in $p\left(v_i\right|z_i)$.

Although $\mathbf{V}$ is multidimensional in our case, the idea presented above can still be applied and it leads to

$$\begin{array}{c}\hfill \frac{{\prod }_{i=1}^{n}p(y_i,{\mathbf{u}}_i|z_i,r_i=1)}{{\Sigma }_{\omega \in \Omega }{\prod }_{i=1}^{n}p(y_{\omega \left(i\right)},{\mathbf{u}}_{\omega \left(i\right)}|z_i,r_i=1)}=\frac{{\prod }_{i=1}^{n}p(y_i,{\mathbf{u}}_i|z_i)}{{\Sigma }_{\omega \in \Omega }{\prod }_{i=1}^{n}p(y_{\omega \left(i\right)},{\mathbf{u}}_{\omega \left(i\right)}|z_i)}.\end{array}$$

(5)

Furthermore, to simplify the computation, we adopt the pairwise fashion of (5) following the previous discussion on pairwise pseudo-likelihood in [15], which results

$$\begin{array}{c}\hfill \prod _{1\le i<j\le n}\frac{p(y_i,{\mathbf{u}}_i|z_i)p(y_j,{\mathbf{u}}_j|z_j)}{p(y_i,{\mathbf{u}}_i|z_i)p(y_j,{\mathbf{u}}_j|z_j)+p(y_i,{\mathbf{u}}_i|z_j)p(y_j,{\mathbf{u}}_j|z_i)}.\end{array}$$

After plugging in model (1) and some algebra, the objective eventually becomes to minimize

$$\begin{array}{c}\hfill L\left(\theta \right)={\left(\genfrac{}{}{0pt}{}{N}{2}\right)}^{-1}\sum _{1\le i<j\le N}{\varphi }_{ij}\left(\theta \right)={\left(\genfrac{}{}{0pt}{}{N}{2}\right)}^{-1}\sum _{1\le i<j\le N}{r}_i{r}_{j}log\{1+W_{ij}exp\left({\theta }^{{}^{\mathrm{T}}}{\mathbf{d}}_{ij}\right)\},\end{array}$$

(6)

where $\theta ={(\tilde{\gamma },{\tilde{\beta }}^{{}^{\mathrm{T}}})}^{{}^{\mathrm{T}}}$, $\tilde{\gamma }=\gamma /{\sigma }^2$, $\tilde{\beta }=\tilde{\gamma }\beta$, ${\mathbf{d}}_{ij}={(-y_{i\setminus j}{z}_{i\setminus j},{\mathbf{u}}_{i\setminus j}^{{}^{\mathrm{T}}}{z}_{i\setminus j})}^{{}^{\mathrm{T}}}$, $y_{i\setminus j}=y_i-y_j$, ${\mathbf{u}}_{i\setminus j}={\mathbf{u}}_i-{\mathbf{u}}_j$, $z_{i\setminus j}=z_i-z_j$ and $W_{ij}=p\left(z_i\right|{\mathbf{u}}_j)p\left(z_j\right|{\mathbf{u}}_i)/\left\{p\left(z_i\right|{\mathbf{u}}_i)p\left(z_j\right|{\mathbf{u}}_j)\right\}$.

Denote the minimizer of (6) as $\widehat{\theta }$. By checking that

$$\frac{{\partial }^2{\varphi }_{ij}\left(\theta \right)}{\partial \theta \partial {\theta }^{{}^{\mathrm{T}}}}=r_i{r}_j{\{1+W_{ij}exp\left({\theta }^{{}^{\mathrm{T}}}{\mathbf{d}}_{ij}\right)\}}^{-2}{W}_{ij}exp\left({\theta }^{{}^{\mathrm{T}}}{\mathbf{d}}_{ij}\right){\mathbf{d}}_{ij}{\mathbf{d}}_{ij}^{{}^{\mathrm{T}}}$$

is positive definite, $\widehat{\theta }$ uniquely exists. To compute $\widehat{\theta }$, one also needs a model for $W_{ij}$. Fortunately, this model only depends on fully observed data ${\mathbf{x}}_i$ and ${\mathbf{x}}_j$. Essentially any existing parametric, semiparametric, or nonparametric modeling technique for $p\left(z\right|\mathbf{u})$ can be used, and $W_{ij}$ can be estimated accordingly. Throughout, we denote ${\widehat{W}}_{ij}$ as an available well-behaved estimator of $W_{ij}$. Although our procedure stems from $p(y,\mathbf{u}|z,r=1)$, which only relies on the data $\{y_i,{\mathbf{x}}_i\}$ with $i=1$, it can be seen that, not only the data $\{y_i,{\mathbf{x}}_i\}$ with $i=1$ are used to compute $\widehat{\theta }$, the data $\left\{{\mathbf{x}}_i\right\}$ with $i=0$ are also used in the process of estimating $W_{ij}$. Therefore, all observed data, both from completely-observed subjects and from partially-observed subjects, are utilized in our procedure.

One can notice that, due to the assumption (2) which allows the greatest flexibility of the mechanism model and the adoption of the conditional likelihood, not all parameters $\alpha$, $\beta$, $\gamma$, and ${\sigma }^2$ are estimable. Nevertheless, the parameter $\beta$, which quantifies the association between Y and $\mathbf{U}$ after adjusting for Z and is of primarily scientific interest, can be fully estimable. The remainder of the paper focuses on the estimation and inference of $\beta$, as well as the variable selection procedure based on $\beta$.

Before moving on, we give some comparison with the existing literature to underline the novel contributions we make in this paper. Based on a slightly different but more restrictive missingness mechanism assumption that $\mathrm{pr}(R=1|Y,\mathbf{X})=a(Y\left)b\right(\mathbf{X})$, Refs. [16,17,18] used the similar idea to analyze non-ignorable missing data for a generalized linear model and a semiparametric proportional likelihood ratio model, respectively. They focused on different aspects of how to use the conditional likelihoods and their consequences such as the partial identifiability issue and the large bias issue. In this paper, we focus on the linear model (1) and we just showed that the parameter $\beta$ is fully identifiable. It can be seen that the method presented in this paper can be applied to different models, but their identifiability problems or some other relevant issues have to be analyzed on a case-by-case basis. For instance, Ref. [19] studied the parameter estimation problem in a logistic regression model with a low dimensionality under assumption (2). They showed that, different from the current paper, all the unknown parameters are identifiable in their context. However, because of the complexity of their objective function, the algorithm studied in [19] is trivial and cannot be extended to a high dimensional setting.

3. Traditional Setting without Regularization

Computation. Directly minimizing $L\left(\theta \right)$ is feasible; however, it is very computationally involved. From rearranging the terms in $L\left(\theta \right)$, we realize that it can be rewritten as the negative log-likelihood function of a standard logistic regression model. To be more specific, let k be the index of pair $(i,j)$ with $k=1,\dots ,K$ and $K=\left(\genfrac{}{}{0pt}{}{n}{2}\right)$. Then,

$$\begin{array}{c}\hfill L\left(\theta \right)=\frac{1}{K}\sum _{k=1}^{K}log\left\{1+exp\left(s_k{\theta }^{{}^{\mathrm{T}}}{\mathbf{t}}_k+log{\widehat{W}}_k\right)\right\},\end{array}$$

(7)

where $s_k=-\mathrm{sign}\left(z_{i\setminus j}\right)$, ${\mathbf{t}}_k=\left(\right|z_{i\setminus j}|y_{i\setminus j},-|z_{i\setminus j}{\left|{\mathbf{u}}_{i\setminus j}^{{}^{\mathrm{T}}}\right)}^{{}^{\mathrm{T}}}$. Denote $g_k=I\{z_{i\setminus j}>0\}$, then one can show that the summand in (7), $log\left\{1+exp\left(s_k{\theta }^{{}^{\mathrm{T}}}{\mathbf{t}}_k+log{\widehat{W}}_k\right)\right\}$, equals,

$$-\left[g_k\left({\theta }^{{}^{\mathrm{T}}}{\mathbf{t}}_k+s_{k}log{\widehat{W}}_k\right)-log\left\{1+exp\left({\theta }^{{}^{\mathrm{T}}}{\mathbf{t}}_k+s_{k}log{\widehat{W}}_k\right)\right\}\right],$$

which is the contribution of the k-th subject to the negative log-likelihood of a logistic regression with $g_k$ as the response, $\theta$ as the coefficient, ${\mathbf{t}}_k$ as the covariate, and $s_{k}log{\widehat{W}}_k$ as the offset term, but without an intercept. Therefore, $\widehat{\theta }$ can be obtained by fitting the aforementioned logistic regression model. Algorithm 1 describes the steps for data manipulation and model fitting to estimate $\theta$ under this traditional setting.

Algorithm 1 Minimization of (6) without penalization

1: Inputs:  $\{y_i,{\mathbf{u}}_i,z_i\},\{y_j,{\mathbf{u}}_j,z_j\},{\widehat{W}}_{ij}$, for $i=1,\dots ,n$ and $j=1,\dots ,n$ 2: Initialize:  $k\leftarrow 0$ 3: for $j\in \{2:n\}$do 4:  for $i\in \{1:(j-1\left)\right\}$ do 5:    $k\leftarrow k+1$ 6:    $y_{i\setminus j}\leftarrow y_i-y_j$, ${\mathbf{u}}_{i\setminus j}\leftarrow {\mathbf{u}}_i-{\mathbf{u}}_j$, $z_{i\setminus j}\leftarrow z_i-z_j$, ${\widehat{W}}_k\leftarrow {\widehat{W}}_{ij}$ 7:    $g_k\leftarrow I\{z_{i\setminus j}>0\}$ 8:    $s_k\leftarrow -\mathrm{sign}\left(z_{i\setminus j}\right)$ 9:    ${\mathbf{t}}_k\leftarrow \left(\right|z_{i\setminus j}|y_{i\setminus j},-|z_{i\setminus j}{\left|{\mathbf{u}}_{i\setminus j}^{{}^{\mathrm{T}}}\right)}^{{}^{\mathrm{T}}}$ 10: Fit logistic regression with response $\mathbf{g}$, covariate $\mathbf{t}$, offset ${\mathbf{s}}^{{}^{\mathrm{T}}}log\widehat{\mathbf{W}}$, and no intercept. 11: Outputs:  $\widehat{\theta }$

Asymptotic Theory. The asymptotic theory of $\widehat{\theta }$ involves a model of $p\left(z\right|\mathbf{u})$, which does not contain any missing values, and therefore any statistical model, either parametric, or semiparametric, or nonparametric, can be used. For simplicity, we only discuss the parametric case here, and any further elaborations will be rendered into Section 7. For a parametric model $p\left(z\right|\mathbf{u};\eta )$, one can apply the standard maximum likelihood estimate $\widehat{\eta }$. Here, we simply assume

$$\begin{array}{c}\hfill \sqrt{N}\left(\widehat{\eta }-{\eta }_0\right)=-{\mathbf{G}}^{-1}\sqrt{N}\frac{1}{N}\sum _{i=1}^N\frac{\partial }{\partial \eta }log\left\{p\left(z_i\right|{\mathbf{u}}_i;{\eta }_0)\right\}+o_p\left(1\right),\end{array}$$

(8)

where $\mathbf{G}=E\left[\frac{{\partial }^2}{\partial \eta \partial {\eta }^{{}^{\mathrm{T}}}}log\left\{p\left(z\right|\mathbf{u};{\eta }_0)\right\}\right]$, $E\parallel \frac{{\partial }^2}{\partial \eta \partial {\eta }^{{}^{\mathrm{T}}}}log\left\{p\left(z\right|\mathbf{u};{\eta }_0)\right\}{\parallel }^2<\infty$, ${\eta }_0$ is the true value of $\eta$, and $\parallel \mathbf{M}\parallel =\sqrt{\mathrm{trace}\left(\mathbf{M}{\mathbf{M}}^{\mathrm{T}}\right)}$ for a matrix $\mathbf{M}$. With this prerequisite, we have the following result for $\widehat{\theta }$, and its proof is provided in Appendix A.

Theorem 1.

Assume (8) as well as $E{∥\frac{{\partial }^2{\varphi }_{ij}({\theta }_0,{\eta }_0)}{\partial \theta \partial {\theta }^{{}^{\mathrm{T}}}}∥}^2<\infty$. Denote ${\theta }_0$ the true value of θ. Then

$$\begin{array}{c}\hfill \sqrt{N}\left(\widehat{\theta }-{\theta }_0\right)\stackrel{d}{\to }N\left(\mathbf{0},{\mathbf{A}}^{-1}\Sigma {\mathbf{A}}^{-1}\right),\end{array}$$

where $\mathbf{A}=E\left\{\frac{{\partial }^2{\varphi }_{ij}({\theta }_0,{\eta }_0)}{\partial \theta \partial {\theta }^{{}^{\mathrm{T}}}}\right\}$, $\Sigma =4E\left\{{\lambda }_{12}({\theta }_0,{\eta }_0){\lambda }_{13}{({\theta }_0,{\eta }_0)}^{{}^{\mathrm{T}}}\right\}$, ${\lambda }_{ij}({\theta }_0,{\eta }_0)=\mathbf{B}{\mathbf{G}}^{-1}{\mathbf{M}}_{ij}\left({\eta }_0\right)-{\mathbf{N}}_{ij}({\theta }_0,{\eta }_0)$, $\mathbf{B}=E\left\{\frac{{\partial }^2{\varphi }_{ij}({\theta }_0,{\eta }_0)}{\partial \theta \partial {\eta }^{{}^{\mathrm{T}}}}\right\}$, ${\mathbf{M}}_{ij}\left({\eta }_0\right)=\frac{1}{2}\left\{\frac{\partial }{\partial \eta }logp\left(z_i\right|{\mathbf{u}}_i;{\eta }_0)+\frac{\partial }{\partial \eta }logp\left(z_j\right|{\mathbf{u}}_j;{\eta }_0)\right\}$, and ${\mathbf{N}}_{ij}({\theta }_0,{\eta }_0)=\frac{\partial {\varphi }_{ij}({\theta }_0,{\eta }_0)}{\partial \theta }$.

If one prefers the asymptotic result of $\widehat{\beta }$, we have

Corollary 1.

Let $\mathbf{C}$ be a $p\times (p+1)$ matrix such that $\mathbf{C}\theta =\beta$, i.e.,

$$\mathbf{C}=\left(\begin{array}{ccccc}0& 1/{\tilde{\gamma }}_0& 0& \cdots & 0\\ 0& 0& 1/{\tilde{\gamma }}_0& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & 1/{\tilde{\gamma }}_0\end{array}\right).$$

Denote ${\beta }_0$ the true value of β. Then, following Theorem 1, we have $\sqrt{N}\left(\widehat{\beta }-{\beta }_0\right)\stackrel{d}{\to }N\left(\mathbf{0},\mathbf{C}{\mathbf{A}}^{-1}\Sigma {\mathbf{A}}^{-1}{\mathbf{C}}^{{}^{\mathrm{T}}}\right)$.

Variance Estimation. With Theorem 1 and Corollary 1, the variance estimation is straightforward using the plugging in strategy. Note that $var\left(\widehat{\theta }\right)=\frac{1}{N}{\mathbf{A}}^{-1}\Sigma {\mathbf{A}}^{-1}$, then one would have the estimate $\widehat{var}\left(\widehat{\theta }\right)=\frac{1}{N}{\widehat{\mathbf{A}}}^{-1}\widehat{\Sigma }{\widehat{\mathbf{A}}}^{-1}$ where $\widehat{\mathbf{A}}={\left(\genfrac{}{}{0pt}{}{N}{2}\right)}^{-1}{\sum }_{1\le i<j\le N}\frac{{\partial }^2{\varphi }_{ij}(\widehat{\theta },\widehat{\eta })}{\partial \theta \partial {\theta }^{{}^{\mathrm{T}}}}$,

$\widehat{\Sigma }=\frac{4}{N-1}{\sum }_{i=1}^N{\left[\frac{1}{N-1}{\sum }_{j=1,j\ne i}^N\left\{\widehat{\mathbf{B}}{\widehat{\mathbf{G}}}^{-1}{\mathbf{M}}_{ij}\left(\widehat{\eta }\right)-{\mathbf{N}}_{ij}(\widehat{\theta },\widehat{\eta })\right\}\right]}^{\otimes 2}$, $\widehat{\mathbf{B}}={\left(\genfrac{}{}{0pt}{}{N}{2}\right)}^{-1}{\sum }_{1\le i<j\le N}\frac{{\partial }^2{\varphi }_{ij}(\widehat{\theta },\widehat{\eta })}{\partial \theta \partial {\eta }^{{}^{\mathrm{T}}}}$, and $\widehat{\mathbf{G}}=\frac{1}{N}{\sum }_{i=1}^N\frac{{\partial }^2}{\partial \eta \partial {\eta }^{{}^{\mathrm{T}}}}log\left\{p\left(z_i\right|{\mathbf{u}}_i;\widehat{\eta })\right\}$.

4. Modern Setting with Regularization

In the past few decades, it has become a standard practice to consider the high-dimensional regression model, where one assumes the parameter $\beta$ is sparse and often uses the regularization technique to recover the sparsity. While it is a prominent problem to analyze this type of model when the data are prone to missing values, the literature is quite scarce primarily because it is cumbersome to rigorously address the missingness under high dimensionality. Therefore, it is valuable to extend the nuisance-free likelihood procedure proposed in Section 3 to the setting with regularization.

Computation. Regularization is a powerful technique to identify the zero elements of a sparse parameter in a regression model. Various penalty functions have been extensively studied, such as LASSO [20], SCAD [21], and MCP [22]. In particular, we study the adaptive LASSO penalty [23] with the objective of minimizing the following function

$$\begin{array}{c}\hfill L_{\lambda }\left(\theta \right)=L\left(\theta \right)+\sum _{j=1}^p\lambda {\left|{\widehat{\tilde{\beta }}}_j\right|}^{-1}\left|{\tilde{\beta }}_j\right|,\end{array}$$

(9)

where $\lambda >0$ is the tuning parameter. Following [23], ${\widehat{\tilde{\beta }}}_j$ is a root-N-consistent estimator of ${\tilde{\beta }}_j$; for example, one can use the estimator via minimizing the unregularized objective Function (6). Obviously, the penalty term in (9) does not alter the numerical characteristic of $L\left(\theta \right)$ that we presented in Section 3. The $L_{\lambda }\left(\theta \right)$ is essentially the regularized log-likelihood of a logistic regression model with the similar format as discussed in (7).

To choose the tuning parameter $\lambda$, one can follow either the cross-validation method or various information-based criteria. Fortunately, all of these approaches have been extensively studied in the literature. In this paper, we follow the Bayesian information criterion (BIC) to determine $\lambda$. Specifically, we choose $\lambda$ to be the minimizer of the following BIC function

$$\begin{array}{c}\hfill \mathrm{BIC}\left(\lambda \right)=2L\left(\theta \right)+p_{\lambda }\frac{log\left(n\right)}{n},\end{array}$$

where $p_{\lambda }$ is the number of nonzero elements in ${\widehat{\tilde{\beta }}}_{\lambda }$ and the minimizer of (9) is encoded as ${\widehat{\theta }}_{\lambda }={({\widehat{\tilde{\gamma }}}_{\lambda },{\widehat{\tilde{\beta }}}_{\lambda }^{\mathrm{T}})}^{\mathrm{T}}$. We summarize the whole computation pipeline as Algorithm 2 below.

Algorithm 2 Minimization of (9) with the ALASSO penalty

1: Inputs:   $\{y_i,{\mathbf{u}}_i,z_i\},\{y_j,{\mathbf{u}}_j,z_j\},{\widehat{W}}_{ij}$, for $i=1,\dots ,n$ and $j=1,\dots ,n$ 2: Initialize:  $k\leftarrow 0$ 3: for$j\in \{2:n\}$ do 4:  for $i\in \{1:(j-1\left)\right\}$ do 5:    $k\leftarrow k+1$ 6:    $y_{i\setminus j}\leftarrow y_i-y_j$, ${\mathbf{u}}_{i\setminus j}\leftarrow {\mathbf{u}}_i-{\mathbf{u}}_j$, $z_{i\setminus j}\leftarrow z_i-z_j$, ${\widehat{W}}_k\leftarrow {\widehat{W}}_{ij}$ 7:    $g_k\leftarrow I\{z_{i\setminus j}>0\}$ 8:    $s_k\leftarrow -\mathrm{sign}\left(z_{i\setminus j}\right)$ 9:    ${\mathbf{t}}_k\leftarrow \left(\right|z_{i\setminus j}|y_{i\setminus j},|z_{i\setminus j}{\left|{\mathbf{u}}_{i\setminus j}^{{}^{\mathrm{T}}}\right)}^{{}^{\mathrm{T}}}$ 10: Fit logistic regression with response $\mathbf{g}$, covariates $\mathbf{t}$, offset ${\mathbf{s}}^{{}^{\mathrm{T}}}log\mathbf{W}$, and no intercept. 11: Obtain $\widehat{\tilde{\theta }}$. 12: Fit logistic regression with ALASSO penalty. 13: Find ${\lambda }^{\star }$ which minimizes the BIC. 14: Outputs:  $\widehat{\theta }\left({\lambda }^{\star }\right)={\widehat{\theta }}_{\lambda }$

Asymptotic Theory. Recall that $\theta ={(\tilde{\gamma },{\tilde{\beta }}^{\mathrm{T}})}^{\mathrm{T}}$. Without loss of generality, we assume the first $p_0$ parameters in $\tilde{\beta }$ are nonzero, where $1\le p_0<p$. For simplicity, we denote ${\theta }_T={(\tilde{\gamma },{\tilde{\beta }}_1,\dots ,{\tilde{\beta }}_{p_0})}^{\mathrm{T}}$ as the vector of nonzero components and ${\theta }_{T^c}={({\tilde{\beta }}_{p_0+1},\dots ,{\tilde{\beta }}_p)}^{\mathrm{T}}$ as the vector of zeros.

In Theorem 1, we defined $\mathbf{A}=E\left\{\frac{{\partial }^2{\varphi }_{ij}({\theta }_0,{\eta }_0)}{\partial \theta \partial {\theta }^{{}^{\mathrm{T}}}}\right\}$, a $(p+1)\times (p+1)$ matrix. Now we assume it can be partitioned as $\mathbf{A}=\left(\begin{array}{cc}{\mathbf{A}}_1& {\mathbf{A}}_2\\ {\mathbf{A}}_2^{{}^{\mathrm{T}}}& {\mathbf{A}}_3\end{array}\right)$, where ${\mathbf{A}}_1$ is a $(p_0+1)\times (p_0+1)$ submatrix corresponding to ${\theta }_T$. Similarly, we defined $\Sigma =4E\left\{{\lambda }_{12}({\theta }_0,{\eta }_0){\lambda }_{13}{({\theta }_0,{\eta }_0)}^{{}^{\mathrm{T}}}\right\}$, and we also assume it can be partitioned as $\Sigma =\left(\begin{array}{cc}{\Sigma }_1& {\Sigma }_2\\ {\Sigma }_2^{{}^{\mathrm{T}}}& {\Sigma }_3\end{array}\right)$, where ${\Sigma }_1$ is a $(p_0+1)\times (p_0+1)$ submatrix corresponding to ${\theta }_T$ as well. We denote the minimizer of (9), ${\widehat{\theta }}_{\lambda }$, as ${\widehat{\theta }}_{\lambda }={({\widehat{\theta }}_{\lambda ,T}^{\mathrm{T}},{\widehat{\theta }}_{\lambda ,T^c}^{\mathrm{T}})}^{\mathrm{T}}$, and its true value ${\theta }_0={({\theta }_{0,T}^{\mathrm{T}},{\theta }_{0,T^c}^{\mathrm{T}})}^{\mathrm{T}}$.

Now, we present the oracle property pertaining to ${\widehat{\theta }}_{\lambda }$, which includes the asymptotic normality for the nonzero components and the variable selection consistency. The proof is provided in Appendix B.

Theorem 2.

Assume (8), ${\mathbf{A}}_1$ is positive definite and $E\parallel \frac{\partial {\varphi }_{ij}({\theta }_0,{\eta }_0)}{\partial \theta }{\parallel }^2<\infty$ for each θ in a neighborhood of ${\theta }_0$. We also assume $\sqrt{N}\lambda \to 0$ and $N\lambda \to \infty$. Then,

$$\sqrt{N}\left({\widehat{\theta }}_{\lambda ,T}-{\theta }_{0,T}\right)\stackrel{d}{\to }N\left(\mathbf{0},{\mathbf{A}}_1^{-1}{\Sigma }_1{\mathbf{A}}_1^{-1}\right).$$

In addition, let $T_N=\{j\in \left\{1,\dots ,p\right\}:{\widehat{\tilde{\beta }}}_{j,\lambda }\ne 0\}$ and $T=\{j\in \left\{1,\dots ,p\right\}:{\tilde{\beta }}_{j,0}\ne 0\}$, then

$$\underset{N\to \infty }{lim}\mathrm{pr}(T_N=T)=1.$$

Variance Estimation. Although the above theory provides a rigorous justification for the asymptotic property of ${\widehat{\theta }}_{\lambda }$, in practice, however, it does not guide the standard error estimation. Here, we propose a data perturbation approach for the variance estimation. Specifically, following [24], we generate a set of independent and identically distributed positive random variables $\Xi =\{{\xi }_i,i=1,\dots ,N\}$ with $E\left({\xi }_i\right)=1$ and $var\left({\xi }_i\right)=1$, e.g., the standard exponential distribution. Since it is based on a U-statistic structure, we perturb our objective function by adding ${\kappa }_{ij}={\xi }_i{\xi }_j$ to each of its pairwise terms. We first obtain the estimator ${\widehat{\theta }}^{\star }$ by minimizing the perturbed version of (6):

$$L^{\star }\left(\theta \right)={\left(\genfrac{}{}{0pt}{}{N}{2}\right)}^{-1}\sum _{1\le i<j\le N}{\kappa }_{ij}{\varphi }_{ij}\left(\theta \right).$$

Then, we obtain the estimator ${\widehat{\theta }}_{\lambda }^{\star }$ by minimizing the perturbed version of (9):

$$L_{\lambda }^{\star }\left(\theta \right)={\left(\genfrac{}{}{0pt}{}{N}{2}\right)}^{-1}\sum _{1\le i<j\le N}{\kappa }_{ij}{\varphi }_{ij}\left(\theta \right)+\sum _{j=1}^p\frac{\lambda }{\left|{\widehat{\tilde{\beta }}}_j^{\star }\right|}\left|{\tilde{\beta }}_j\right|,$$

where the optimal $\lambda$ is also computed by the BIC. We repeat this data perturbation scheme a large number of times, say, M.

Following the theory in [25,26], under some regularity conditions, one can first show that $\sqrt{N}\left({\widehat{\theta }}_{\lambda ,T}^{\star }-{\theta }_{0,T}\right)$ converges in distribution to $N(\mathbf{0},{\mathbf{A}}_1^{-1}{\Sigma }_1{\mathbf{A}}_1^{-1})$, the same limiting distribution of $\sqrt{N}\left({\widehat{\theta }}_{\lambda }-{\theta }_0\right)$. Furthermore, one can also show ${\mathrm{pr}}^{*}\left({\widehat{\theta }}_{\lambda ,T^c}^{\star }=0\right)\to 1$, where ${\mathrm{pr}}^{*}$ is the probability measure generated by the original data $\mathcal{X}$ and the perturbation data $\Xi$. In addition, one can show that the distribution of $\sqrt{N}\left({\widehat{\theta }}_{\lambda ,T}^{\star }-{\widehat{\theta }}_{\lambda ,T}\right)$ conditional on the data can be used to approximate the unconditional distribution of $\sqrt{N}\left({\widehat{\theta }}_{\lambda ,T}-{\theta }_{0,T}\right)$ and that ${\mathrm{pr}}^{*}\left({\widehat{\theta }}_{\lambda ,T^c}^{\star }=0|\mathcal{X}\right)\to 1$.

To achieve a confidence interval for ${\theta }_j$, the j-th coordinate in $\theta$, the lower and upper bounds can be formed by ${\widehat{\theta }}_{\lambda ,j,\alpha /2}^{\star }$ and ${\widehat{\theta }}_{\lambda ,j,1-\alpha /2}^{\star }$, respectively, where ${\widehat{\theta }}_{\lambda ,j,q}^{\star }$ represents the q-th quantile of $\left\{{\widehat{\theta }}_{\lambda ,j,m}^{\star },m=1,\dots ,M\right\}$.

5. Simulation Studies

We conduct comprehensive simulation studies to evaluate the finite sample performance of our proposed estimators and also compare with some currently existing methods. We first present the results under the model without regularization, then with regularization.

5.1. Scenarios without Regularization

For the proposed estimator studied in Section 3, we generate $\{R_i,Y_i,{\mathbf{U}}_i^{\mathrm{T}},Z_i\},i=1,\dots ,N$, independent and identically distributed copies of $(R,Y,{\mathbf{U}}^{\mathrm{T}},Z)$, as follows. We first generate the random vector $\mathbf{U}={(U_1,\dots ,U_p)}^{\mathrm{T}}$ with $U_i\sim N(0.5,1)$ and $p=4$, and then generate $Z={\alpha }_z+{\eta }^{\mathrm{T}}\mathbf{U}+{ϵ}_z$ with ${\alpha }_z=0.5$, $\eta ={(-0.5,1,-1,1.5)}^{\mathrm{T}}$, ${ϵ}_z\sim N(0,1)$. Afterwards, the outcome variable Y is generated following the model (1) with $\alpha =-1$, $\beta ={(-0.5,1,-1,1.5)}^{\mathrm{T}}$, $\gamma =0.5$, and $ϵ\sim N(0,1)$, and the missingness indicator R is generated following $\mathrm{pr}(R=1|Y,\mathbf{U})=I(Y<2.5,U_1<2,U_2<2,U_3<2,U_4<2)$ which results in around 40% missing values. We examine two situations with sample size $N=500$ and $N=1000$ respectively. Besides the estimator studied in Section 3 (Proposed), we also implement the estimator using all simulated data (FullData) and the estimator using completely observed subjects only (CC). Based on 1000 simulation replicates, for each of the three estimators, we summarize the sample bias, sample standard deviation, estimated standard error, and coverage probability of 95% confidence intervals in

Table 1. In Section 5.1, sample bias (Bias), sample standard deviation (SD), estimated standard error (SE), and coverage probability (CP) of 95% confidence interval of the estimator of FullData (using all simulated data), CC (using only completely observed subjects), and of the proposed estimator studied in Section 3.

N	Parameter	Method	Bias	SD	SE	CP
500	$\tilde{\gamma }$	FullData	0.0026	0.0444	0.0450	0.9540
		CC	−0.0329	0.0564	0.0560	0.9100
		Proposed	0.0174	0.0829	0.0789	0.9450
	${\beta }_1$	FullData	0.0022	0.0489	0.0503	0.9510
		CC	0.0376	0.0670	0.0699	0.9300
		Proposed	0.0164	0.1644	0.1607	0.9400
	${\beta }_2$	FullData	−0.0017	0.0657	0.0635	0.9310
		CC	−0.0649	0.0851	0.0835	0.8680
		Proposed	−0.0399	0.2305	0.2239	0.9360
	${\beta }_3$	FullData	0.0022	0.0616	0.0635	0.9540
		CC	0.0778	0.0871	0.0867	0.8430
		Proposed	0.0462	0.2323	0.2298	0.9410
	${\beta }_4$	FullData	−0.0045	0.0792	0.0810	0.9530
		CC	−0.0988	0.1007	0.1043	0.8550
		Proposed	−0.0672	0.3081	0.3047	0.9380
1000	$\tilde{\gamma }$	FullData	−0.0012	0.0317	0.0317	0.9540
		CC	−0.0348	0.0396	0.0393	0.8510
		Proposed	0.0068	0.0573	0.0555	0.9350
	${\beta }_1$	FullData	0.0011	0.0367	0.0355	0.9370
		CC	0.0399	0.0490	0.0494	0.8840
		Proposed	0.0154	0.1154	0.1138	0.9460
	${\beta }_2$	Full Data	0.0020	0.0448	0.0448	0.9500
		CC	−0.0649	0.0577	0.0588	0.8110
		Proposed	−0.0153	0.1531	0.1591	0.9590
	${\beta }_3$	Full Data	−0.0015	0.0458	0.0449	0.9460
		CC	0.0779	0.0605	0.0611	0.7490
		Proposed	0.0135	0.1598	0.1634	0.9480
	${\beta }_4$	Full Data	0.0009	0.0564	0.0571	0.9540
		CC	−0.0949	0.0720	0.0734	0.7550
		Proposed	−0.0242	0.2091	0.2167	0.9430

.

Furthermore, we consider a similar simulation setting where the generation is the same as above except for a logistic missingness mechanism model with $\mathrm{logit}\left\{\mathrm{pr}(R=1|Y,\mathbf{U})\right\}=3-2Y+0.5U_1-U_2+U_3-1.5U_4$, which also results in around 40% missing values. We replicate the results, shown in

Table 2. In Section 5.1, sample bias (Bias), sample standard deviation (SD), estimated standard error (SE), and coverage probability (CP) of 95% confidence interval of the estimator of FullData (using all simulated data), CC (using only completely observed subjects), and of the proposed estimator studied in Section 3, with a logistic missingness mechanism model.

N	Parameter	Method	Bias	SD	SE	CP
500	$\tilde{\gamma }$	FullData	−0.0011	0.0464	0.0451	0.9410
		CC	−0.0306	0.0567	0.0567	0.9200
		Proposed	0.0100	0.0822	0.0787	0.9380
	${\beta }_1$	FullData	−0.0004	0.0509	0.0503	0.9520
		CC	0.0440	0.0636	0.0637	0.8930
		Proposed	0.0146	0.1308	0.1236	0.9420
	${\beta }_2$	FullData	0.0013	0.0639	0.0637	0.9520
		CC	−0.0871	0.0828	0.0821	0.8190
		Proposed	−0.0173	0.1824	0.1753	0.9430
	${\beta }_3$	FullData	−0.0030	0.0655	0.0636	0.9400
		CC	0.0876	0.0847	0.0821	0.8030
		Proposed	0.0214	0.1840	0.1756	0.9440
	${\beta }_4$	FullData	0.0023	0.0845	0.0812	0.9390
		CC	−0.1307	0.1083	0.1061	0.7560
		Proposed	−0.0331	0.2533	0.2384	0.9360
1000	$\tilde{\gamma }$	FullData	0.0004	0.0315	0.0317	0.9490
		CC	−0.0286	0.0396	0.0398	0.8950
		Proposed	0.0060	0.0568	0.0555	0.9390
	${\beta }_1$	FullData	0.0007	0.0362	0.0354	0.9420
		CC	0.0442	0.0451	0.0447	0.8410
		Proposed	0.0079	0.0910	0.0859	0.9290
	${\beta }_2$	FullData	−0.0004	0.0450	0.0448	0.9390
		CC	−0.0879	0.0571	0.0576	0.6640
		Proposed	−0.0044	0.1277	0.1220	0.9420
	${\beta }_3$	FullData	−0.0009	0.0450	0.0448	0.9450
		CC	0.0880	0.0588	0.0577	0.6660
		Proposed	0.0114	0.1309	0.1222	0.9380
	${\beta }_4$	FullData	−0.0005	0.0576	0.0572	0.9510
		CC	−0.1342	0.0755	0.0745	0.5740
		Proposed	−0.0191	0.1757	0.1661	0.9370

.

We can reach the following conclusions from Table 1 and Table 2. For the estimator Proposed, although its bias is slightly larger than the benchmark FullData, it is still very close to zero. The sample standard deviation and the estimated standard error are rather close to each other. The sample coverage probability of the estimated 95% confidence interval is also very close to the nominal level. This observation well matches our theoretical justification in Theorem 1. On the contrary, the estimator CC is clearly biased, resulting in empirical coverage far from the nominal level, and therefore is not recommended to use in practice. It is also clear that, compared to the benchmark FullData, the estimator Proposed has estimation efficiency loss to some extent. This is because the proposed method uses the conditional likelihood approach and it completely eliminates the effect of the nuisance.

5.2. Scenarios with Regularization

For the estimator studied in Section 4, the independent and identically distributed samples are generated as follows. The variable $\mathbf{U}={(U_1,\dots ,U_p)}^{\mathrm{T}}$ is generated from $\mathrm{MVN}(\mathbf{0},{\Sigma }_u)$ with ${\Sigma }_u={(0.5^{|i-j|})}_{1\le i,j\le p}$ and $p=8$. Then, the shadow variable Z is generated following $Z={\alpha }_z+{\eta }^{\mathrm{T}}\mathbf{U}+{ϵ}_z$ with ${\alpha }_z=0$, $\eta ={(-0.5,0.5,-1,1,-0.5,0.5,-1,1)}^{\mathrm{T}}$ and ${ϵ}_z\sim N(0,1)$. The outcome variable Y is generated from model (1) with $\alpha =0$, $\beta ={(3,1.5,0,0,2,0,0,0)}^{\mathrm{T}}$, $\gamma =3$, $ϵ\sim N(0,{\sigma }^2)$ and $\sigma =3$. The distribution of the missingness indicator follows from $\mathrm{logit}\left\{\mathrm{pr}(R=1|Y,\mathbf{U})\right\}=5+5Y+0.2U_1+0.2U_7$, which results in about 45% missing values. Similar to Section 5.1, we also examine two situations with sample size $N=500$ and $N=1000$ respectively, and we implement three estimators FullData, CC, and Proposed. When the estimator Proposed is implemented, we perform $M=500$ perturbations in order to obtain the confidence interval for the unknown parameter. The results summarized below are based on 1000 simulation replicates.

Figure 1 shows the $L_1$, $L_2$, and $L_{\infty }$ norms of the bias for the three different estimators. As sample size increases, there is no doubt that the estimation bias is getting smaller for any method. It is also clear that the bias of the Proposed estimator is larger than the benchmark FullData, but much smaller than the method CC.

We present the statistical inference results in

Table 3. In Section 5.2, with sample size $N=500$, sample bias (Bias), sample standard deviation (SD), estimated standard error (SE), coverage probability (CP), and length (Length) of 95% confidence interval of the estimator of FullData (using all simulated data), CC (using only completely observed subjects) and of the proposed estimator studied in Section 4.

Parameter		Method	Bias	SD	SE	CP	Length
$\tilde{\gamma }$		FullData	0.0001	0.0120	0.0132	0.9480	0.0515
		CC	−0.0729	0.0180	0.0183	0.0370	0.0716
		Proposed	−0.0423	0.0500	0.0498	0.8200	0.1926
True Nonzero	${\beta }_1$	FullData	0.0021	0.1686	0.1649	0.9400	0.6415
		CC	−0.6547	0.2207	0.2114	0.1460	0.8233
		Proposed	0.0354	0.4698	0.4746	0.9320	1.8513
	${\beta }_2$	Full Data	−0.0275	0.1692	0.1791	0.9440	0.6952
		CC	−0.3501	0.2227	0.2174	0.6180	0.8471
		Proposed	−0.2654	0.5843	0.5609	0.8940	1.9237
	${\beta }_5$	Full Data	−0.0172	0.1576	0.1756	0.9650	0.6826
		CC	−0.4478	0.2172	0.2161	0.4370	0.8418
		Proposed	−0.1251	0.4037	0.4611	0.9330	1.8063
True Zero	${\beta }_3$	FullData	0.0085	0.1567	0.1890	0.9960	0.7184
		CC	0.0063	0.2067	0.2304	0.9890	0.8890
		Proposed	0.0109	0.0988	0.1690	1.0000	0.4398
	${\beta }_4$	Full Data	−0.0019	0.1581	0.1900	0.9940	0.7206
		CC	−0.0017	0.2097	0.2307	0.9900	0.8914
		Proposed	0.0126	0.1112	0.1447	1.0000	0.3668
	${\beta }_6$	Full Data	0.0045	0.1212	0.1606	0.9980	0.6146
		CC	−0.0053	0.1749	0.1953	0.9900	0.7560
		Proposed	0.0034	0.0664	0.1160	1.0000	0.2555
	${\beta }_7$	Full Data	0.0014	0.1351	0.1839	0.9980	0.7063
		CC	−0.0055	0.1870	0.2245	0.9950	0.8717
		Proposed	0.0024	0.0386	0.1115	1.0000	0.2538
	${\beta }_8$	Full Data	−0.0072	0.1295	0.1748	0.9990	0.6653
		CC	−0.0062	0.1795	0.2125	0.9940	0.8251
		Proposed	0.0016	0.0741	0.1066	1.0000	0.2284

for $N=500$ and

Table 4. In Section 5.2, with sample size $N=1000$, sample bias (Bias), sample standard derivation (SD), estimated standard error (SE), coverage probability (CP), and length (Length) of 95% confidence interval of the estimator of FullData (using all simulated data), CC (using only completely observed subjects) and of the proposed estimator studied in Section 4.

Parameter		Method	Bias	SD	SE	CP	Length
$\tilde{\gamma }$		FullData	−0.0005	0.0073	0.0088	0.9690	0.0344
		CC	−0.0730	0.0126	0.0130	0.0000	0.0507
		Proposed	−0.0213	0.0311	0.0334	0.8700	0.1293
True Nonzero	${\beta }_1$	FullData	−0.0005	0.1186	0.1170	0.9300	0.4547
		CC	−0.6655	0.1568	0.1507	0.0090	0.5864
		Proposed	0.0211	0.2911	0.2969	0.9300	1.1631
	${\beta }_2$	Full Data	−0.0321	0.1175	0.1249	0.9550	0.4861
		CC	−0.3387	0.1477	0.1534	0.3960	0.5972
		Proposed	−0.0979	0.2907	0.3383	0.9230	1.3115
	${\beta }_5$	Full Data	−0.0225	0.1051	0.1206	0.9590	0.4698
		CC	−0.4485	0.1478	0.1534	0.1770	0.5964
		Proposed	−0.0621	0.2351	0.2526	0.9290	0.9871
True Zero	${\beta }_3$	FullData	−0.0007	0.0621	0.1162	1.0000	0.4253
		CC	0.0023	0.1414	0.1614	0.9920	0.6180
		Proposed	0.0044	0.0581	0.0910	1.0000	0.2091
	${\beta }_4$	Full Data	0.0020	0.0632	0.1170	1.0000	0.4271
		CC	−0.0005	0.1333	0.1608	0.9930	0.6207
		Proposed	0.0063	0.0584	0.0887	1.0000	0.2107
	${\beta }_6$	Full Data	0.0013	0.0571	0.1010	1.0000	0.3670
		CC	−0.0034	0.1159	0.1378	0.9950	0.5313
		Proposed	0.0012	0.0281	0.0688	1.0000	0.1430
	${\beta }_7$	Full Data	−0.0028	0.0599	0.1144	1.0000	0.4231
		CC	−0.0033	0.1243	0.1584	0.9970	0.6131
		Proposed	0.0016	0.0288	0.0698	1.0000	0.1421
	${\beta }_8$	Full Data	0.0039	0.0589	0.1080	1.0000	0.3970
		CC	0.0028	0.1256	0.1497	0.9940	0.5752
		Proposed	0.0000	0.0333	0.0644	1.0000	0.1314

for $N=1000$, respectively, including sample bias, sample standard deviation, estimated standard error, coverage probability, and length of 95% confidence interval for the three different methods. For the nonzero $\beta$’s as well as $\tilde{\gamma }$, similar to Section 5.1, the method CC clearly prompts coverage probability far from the nominal level hence is not reliable. For the method Proposed, its estimation bias is quite close to zero, and its sample standard deviation and estimated standard error are quite close to each other. The coverage probability of the confidence interval converges to the nominal level 95% as the sample size gets larger. For the noisy zero $\beta$’s, the coverage probabilities in the three methods are all close to 1, reflecting the variable selection consistency in the oracle property, even for the CC method. Furthermore, a very nice finite sample property of our proposed estimator is that it produces the confidence interval with the shortest length, which can be clearly seen from both Table 3 and Table 4.

6. Real Data Application

The Medical Information Mart for Intensive Care III (MIMIC-III) is an openly available electronic health records (EHR) database, developed by the MIT Lab for Computational Physiology [13], comprising de-identified health-related data associated with intensive care unit patients with rich information including demographics, vital signs, laboratory test, medications, and more.

Our initial motivation for this data analysis is to understand the missingness mechanism for some laboratory test biomarkers in this EHR system. As for the EHR database, since the data are collected in a non-prescheduled fashion, i.e., only available when the patient seeks care or the physician orders care, the visiting process could be potentially informative about the patients’ risk categories. Therefore, it is very plausible that the data are missing not at random, or a mix of missing not at random and missing at random [27,28]. When we first conducted the data cleaning process briefly, an interesting phenomenon we observe is that, compared to most biomarkers which usually have <3% missing values, the albumin level in the blood sample, a very indicative biomarker associated with different types of diseases [29], has around 30% missingness.

To further understand this phenomenon, we concentrate on a subset of the data with sample size $N=1359$ in which 421 samples have missing values in the albumin level but all other variables are complete. We aim to apply the proposed method to the study of the albumin level (Y). The calcium level in the blood sample, free of missing data, has been shown in the biomedical literature that it has high correlation with the albumin level [30,31,32]; therefore, we adopt the calcium level as the shadow variable Z. Seventeen other variables comprise the vector $\mathbf{U}$, which are either demographics (age and gender), chart events (respiratory rate, glucose, heart rate, systolic blood pressure, diastolic blood pressure, and temperature), other laboratory tests (urea nitrogen, platelets, magnesium, hematocrit, red blood cell, white blood cell, and peripheral capillary oxygen saturation (SpO2)), or aggregated metrics (simplified acute physiology score (SAPS-II) and sequential organ failure assessment score (SOFA)).

We implement the proposed estimator studied in Section 4 to achieve both variable selection and post-selection inference. We also compare it with the CC method which naively fits the regularized linear regression with the ALASSO penalty. For each of the methods, we apply the data perturbation scheme presented in Section 4 with $M=500$ for standard error estimation. The results are summarized in

Table 5. In Section 6, the parameter estimate (Estimate), standard error (SE), and confidence interval (CI) of the estimator of CC (using only completely observed subjects) and of the proposed estimator studied in Section 4 in the MIMIC−III study.

Effect	CC			Proposed
	Estimate	SE	CI	Estimate	SE	CI
Calcium(shadow)	0.7707	0.0691	[0.6532, 0.9153]	1.5271	0.1796	[1.1815, 1.8835]
Red Blood Cell	0.6491	0.0514	[0.5337, 0.7257]	0.7545	0.1631	[0.3594, 1.0109]
Magnesium	0.0000	0.0686	[−0.2073, 0.0000]	0.2731	0.2452	[0.0000, 0.6609]
SOFA	−0.2720	0.0268	[−0.3135, −0.2099]	−0.1852	0.1040	[−0.3467, 0.0000]
Temperature	−0.0360	0.0351	[−0.0883, 0.0659]	0.0000	0.0964	[0.0000, 0.3132]
White Blood Cell	−0.0245	0.0123	[−0.0416, 0.0000]	0.0000	0.0025	[0.0000, 0.0000]
Age	0.0000	0.0008	[0.0000, 0.0000]	0.0000	0.0017	[0.0000. 0.0000]
Gender	0.0000	0.0240	[−0.0477, 0.0662]	0.0000	0.1320	[−0.4025, 0.0000]
Respiratory Rate	0.0000	0.0034	[−0.0141, 0.0000]	0.0000	0.0008	[0.0000, 0.0000]
Glucose	0.0000	0.0000	[0.0000, 0.0000]	0.0000	0.0005	[0.0000, 0.0000]
Heart Rate	0.0000	0.0025	[−0.0091, 0.0000]	0.0000	0.0004	[0.0000, 0.0000]
Systolic BP	0.0000	0.0045	[−0.0139, 0.0000]	0.0000	0.0000	[0.0000, 0.0000]
Diastolic BP	0.0000	0.0072	[0.0000, 0.0223]	0.0000	0.0000	[0.0000, 0.0000]
Urea Nitrogen	0.0000	0.0004	[0.0000, 0.0000]	0.0000	0.0000	[0.0000, 0.0000]
Platelets	0.0000	0.0000	[0.0000, 0.0000]	0.0000	0.0000	[0.0000, 0.0000]
Hematocrit	0.0000	0.0027	[0.0000, 0.0000]	0.0000	0.0000	[0.0000, 0.0000]
SpO2	0.0000	0.0145	[−0.0479, 0.0000]	0.0000	0.0162	[0.0000, 0.0000]
SAPS-II	0.0000	0.0106	[−0.0051, 0.0269]	0.0000	0.0000	[0.0000, 0.0000]

. The solution path of the Proposed method, as the tuning parameter $\lambda$ varies, is also provided in Figure 2.

In general, both methods achieve the goal of variable selection and post-selection inference by leveraging the regularization technique coupled with the data perturbation strategy, and identify many variables as noise with zero coefficients. In particular, the Proposed method provides larger effects for the calcium level (the shadow variable) and the red blood cell count, whereas a smaller effect for the aggregated SOFA score. The Proposed method simplifies the body temperature and the white blood cell count as nonsignificant variables, which are identified as nonzero but with a very small effect using the CC method. It is also worthwhile to mention that the Proposed method signifies the magnesium level with a quite significant coefficient, which was extensively investigated in the scientific literature [33,34,35].

7. Discussion

In this paper, we provide a systematic approach for parameter estimation and statistical inference in both traditional linear model where the regularization is not needed and the modern regularized regression setting, when the outcome variable is prone to missing values and the missingness mechanism can be arbitrarily flexible. A pivotal condition rooted in our procedure is the shadow variable Z, which overcomes the model identifiability problem and enables the nuisance-free conditional likelihood process.

Certainly any method would have its own limitations and could be potentially improved. One needs a model $p\left(z\right|\mathbf{u})$ to implement the proposed estimator in Section 3 and Section 4. As its modeling does not involve any missing data, we simply use the parametric maximum likelihood estimation in our algorithm as well as in the theoretical justification. Indeed, any statistical or machine learning method can be used for modeling $p\left(z\right|\mathbf{u})$. For instance, if one would like to consider a semiparametric model [36], e.g.,

$$\begin{array}{c}\hfill p\left(z\right|\mathbf{u};\eta ,F)=\frac{exp\left({\eta }^{\mathrm{T}}\mathbf{u}z\right)f\left(z\right)}{\int exp\left({\eta }^{\mathrm{T}}\mathbf{u}z\right)dF\left(z\right)},\end{array}$$

where $\eta ={({\eta }_1,\dots ,{\eta }_p)}^{{}^{\mathrm{T}}}$ is a vector of unknown parameters and $f\left(z\right)$ is the density of an unknown baseline distribution function F with respect to some dominating measure $\nu$. With this model fitted, $W_{ij}$ can be simplified to $W_{ij}=exp(-z_{i\setminus j}{\eta }^{{}^{\mathrm{T}}}{\mathbf{u}}_{i\setminus j})$. Therefore, a similar conditional likelihood approach can be used to estimate $\eta$ without estimating the nonparametric component $f\left(z\right)$.