Projective Resampling Imputation Mean Estimation Method for Missing Covariates Problem

In this paper, let ${\bm{Y}}=(Y_{1},Y_{2},\cdots,Y_{n})^{\top}$ be the response variable of interest and ${\bm{X}}=\{X_{ij}:i=1,2,\cdots,n,j=1,2,\cdots,p\}$ be the covariate matrix. We assume that $n>p$. We consider the presence of missing covariates by $r_{ij}$, which denotes the missing-data indicator for $X_{ij}$, where $r_{ij}$ is 1 if $X_{ij}$ is missing and is 0 otherwise. For each unit $i$, $A_{i}=\{j:r_{ij}=0,j=1,2,\cdots,p\}$ denotes the available covariates set, e.g., $A_{i}=\{1,2,\cdots,p\}$ for the complete case.

In this study, we focus on a linear regression model. Assume that the random sample $\{(Y_{i},{\bm{X}}_{i}):i=1,2,\cdots,n,\ j=1,2,\cdots,p\}$ is generated by:

where ${\bm{\beta}}=(\beta_{1},\beta_{2},\cdots,\beta_{p})^{\top}$ is the coefficient vector for the covariates and the $\varepsilon_{i}$’s are independently identically distributed random errors. We assume $E(\varepsilon_{i}|{\bm{X}}_{i})=0$ and $E(\varepsilon_{i}^{2}|{\bm{X}}_{i})=\sigma^{2}$.

Covariates can be divided into two parts based on $A_{i}$: ${\bm{X}}_{i,A_{i}}=(X_{ij}:j\in A_{i})^{\top}$ for observed covariates and ${\bm{X}}_{i,\bar{A}_{i}}=(X_{ij}:j\notin A_{i})^{\top}$ for missing covariates. Thus, equation (1) can be expressed as

where ${\bm{\beta}}_{A_{i}}$ and ${\bm{\beta}}_{\bar{A}_{i}}$ denote the regression coefficients for the complete and incomplete covariates, respectively. For ${\bm{X}}_{i,\bar{A}_{i}}$, which is unobserved, refer to Lin et al. [9], we take the expectation of $Y_{i}$ given the observed covariates. Thus, we obtain the following equation:

Hence, by equation (2), we can impute the incomplete part using the information on ${\bm{X}}_{i,A_{i}}$ to obtain an estimator for ${\bm{\beta}}$. We use the following estimator to estimate the missing components of the covariates for unit $i$:

In equation (3), $K_{h}(\cdot)=K(\cdot/h)/h$, where $K(\cdot)$ is a kernel function and $h$ is a bandwidth. In this way, we can make use of the information as fully as possible. To tackle the problem of “the curse of dimensionality”, we further transform the estimator in equation (3) using the projective resampling method. As shown in Figure 1, for subject $i^{\prime}$, we project ${\bm{X}}_{i^{\prime},A_{i}}$ onto random directions $\{{\bm{v}}_{b,A_{i}}\in\mathcal{R}^{1\times|A_{i}|},b=1,2,\cdots,B\}$, where $|A|$ denotes the cardinality of a set A, then obtain $B$ kernel values using the resulting scalars $\{{\bm{X}}_{i^{\prime},A_{i}}^{\top}{\bm{v}}_{b,A_{i}},b=1,2,\cdots,B\}$, and finally integrate the $B$ kernels through the geometric mean.

where ${\bm{v}}_{b,A_{i}}$ is a random vector, with each entry $v_{b,A_{i},j}$ chosen independently from a distribution $\mathcal{D}$ that is symmetric about the origin with $E(v_{b,A_{i},j}^{2})=1$. In practice, we usually generate $v_{b,A_{i},j}$ from $N(0,1)$ or $U(-1,1)$, but random vectors from other possible distributions can also be used.

Applying the above imputation strategy, we can obtain ${\bm{Z}}_{i}$ by using observed data for part $A_{i}$ and imputed data for part $\bar{A}_{i}$.

Therefore, we propose the following estimation equation:

Thus, the estimator of the regression coefficient can be solved by

Specifically, we propose the following algorithm:

Because the number of disease-associated biomarkers is not expected to be large, it is of great importance to take the sparse assumption into account when not all variables contribute to outcome variables. Hence, we assume the linear regression model in equation (1) is sparse and define the index set of the active and inactive predictors by $\mathcal{I}_{1}=\{j:\beta_{j}\neq 0\}$ and $\mathcal{I}_{2}=\{j:\beta_{j}=0\}$, respectively. Our practical goal is to identify which biomarkers in the CSA-AKI datasets are disease-related as well as to estimate the corresponding coefficients. The main idea of sparse PRIME is to replace the estimation equation 6 with the penalized estimation equations as follows:

where $\lambda\gamma|\beta_{j}|^{\gamma-1}sign(\beta_{j})$ is the partial derivative for the penalty function with respect to $\beta_{j}$. The least absolute shrinkage and selection operator (LASSO) estimator is defined to satisfy $\gamma=1$. Optimizing the functions in (7) with $\gamma=1$ is computationally cumbersome because the functions are non-differentiable. Fortunately, the shooting algorithm proposed in Fu [16] can be used to compute the LASSO estimator. Moreover, Fu [16, 17] proved that the unique estimator of (7) is equivalent to the solution of the penalized objective function as follows:

We can solve this penalized regression-form problem and get a sparse estimator $\hat{\bm{\beta}}_{s}$ using the glmnet package in R.

To decrease the processing time, we can use the following sparse random projection with i.i.d. entries:

where Achlioptas [18] used $s_{i}=1$ or $s_{i}=3$ and Li et al. [19] showed that one can use $s_{i}=\sqrt{|A_{i}|}$ or even $s_{i}=\frac{|A_{i}|}{\log|A_{i}|}$ to significantly reduce the computing time with little loss of accuracy.

The data generation model has the linear expression

where $n=100,200$, $p=12$, ${\bm{\beta}}=(1,-0.6,1.5,1,1.2,0.4,-1,-0.7,1.3,0.5,1.1,-1.4,0.9)^{\top}$.

We generate $(X_{i1},\cdots,X_{ip})$ from the multivariate normal distribution $N_{p}(0,\Sigma)$. We set the non-diagonal elements $\rho_{ij}$ of $\Sigma$ equal to $0.5$ and the diagonal elements of $\Sigma$ equal to $1$. For $\varepsilon_{i}$, we use the error distribution $N(0,\sigma^{2})$, where $\sigma^{2}$ changes with $R^{2}=\text{Var}({\bm{X}}_{i}^{\top}{\bm{\beta}})/\{\text{Var}({\bm{X}}_{i}^{\top}{\bm{\beta}})+\sigma^{2}\}$. We consider the cases in which $R^{2}=0.1,0.2,\cdots,0.9$.

Missing data are divided into three generative scenarios or assumptions. Missing completely at random (MCAR) means the missingness is independent of the values of the data. Missing at random (MAR) means the propensity of data to be missing depends on the observed values, whereas missing not at random (MNAR) covers the remaining scenario that the mechanism depends on the unobserved values (the variables that are missing). In our study, we consider situations that differ from the classical MCAR, MAR, and MNAR mechanism.

For each sample, variables $X_{10},X_{11},X_{12}$ are always available. There are twelve “typical” missing patterns considered, the details are shown in Table 1. We divide the 12 patterns into two groups. Specifically, the first group ${\bm{A}}^{(1)}$ consists of $(A_{1},A_{2},\cdots,A_{6})$, and the second group ${\bm{A}}^{(2)}$ consists of the rest missing patterns. We randomly assign the missing patterns in ${\bm{A}}^{(1)}$ to the sample with missing probability $P=a$. Furthermore, we set the missing probability of $i$th unit for the patterns in ${\bm{A}}^{(2)}$ as $P=\{1+\exp(b\varepsilon_{i}+c)\}^{-1}$. Then we randomly assign the patterns in ${\bm{A}}^{(2)}$ to the missing samples. The settings in Table 2 are used for the missing rate (MR).

The MSE results are shown in Figures 3, 6, 9, and 12 for $R^{2}=0.2,0.5,0.8$. The results of optimal rate of MSE are displayed in Figures 4, 7, 10, and 13. To determine the relative performance, we rank the NADs of the five methods at each repetition. The mean NAD ranks are displayed in Figures 2, 5, 8, and 11 for $R^{2}=0.2,0.5,0.8$. The results of the cases not shown here are available in the supplementary materials. In general, they exhibit patterns which is similar to those shown here.

To make the comparison of Full, PRIME, and ILSE easier, the MSE bar plots for CC and ML are manually scaled because these two methods lead to much higher MSEs than the other methods. The main conclusions are as follows:

1. 1.

   When $n$ and $R^{2}$ increase, the MSEs of Full, ILSE, and PRIME generally decrease, as expected. However, the relative performance of these three methods does not change.

2. 2.

   Generally, PRIME outperformes ILSE and ML in terms of NAD, and all three significantly outperforme CC. Surprisingly, PRIME was close or even superior to the Full method in estimating the coefficients of $X_{11}$ and $X_{12}$ when $\text{MR}=90\%,n=100,$ and $R^{2}=0.2$ and in estimating the coefficients of $X_{10}$ and $X_{12}$ when $\text{MR}=90\%,n=200,$ and $R^{2}=0.2$. These results confirm the superiority of the PRIME method.

3. 3.

   The comparison of results shows that the proposed PRIME method performed better than the other three methods (ILSE, CC, and ML). The bias and variance decomposition figures show that ILSE produces more biased estimates than PRIME. The CC method?s estimation has extremely high variance in almost all ranges of $R^{2}$, and CC?s approximation error was larger when $R^{2}$ was not very high. Furthermore, the CC method produces biased estimates, as expected, because of the missing-data mechanism. The ML performance is not stable: ML’s performance is close to our proposed method?s when $n=200,\text{MR}=90\%,$ and $R^{2}=0.8$, but ML has the poorest performance when $n=100,\text{MR}=90\%,$ and $R^{2}=0.8$ because ML estimators can be highly biased when the MAR assumption does not hold.

4. 4.

   The optimal rates of MSE show a proportion over $40\%$ for PRIME, and this indicates that PRIME yields the smallest MSE of all the competitors in more than 40% of the trials in Scenario 1. When PRIME yields the optimal rate, ILSE or ML most often yields the second-smallest MSE. When PRIME does not yield the lowest MSE, ILSE and ML most often do. Not surprisingly, CC rarely produces the smallest MSE except when $\text{MR}=60\%$ and $R^{2}=0.9$.

To compare the methods with different correlations, we consider the correlation between $X_{i}$ and $X_{j}$ in four situations: $p_{1}$ for $\rho_{ij}=0.2$, $p_{2}$ for $\rho_{ij}=0.5$, $p_{3}$ for $\rho_{ij}=0.8$ and $p_{4}$ for $\rho_{ij}=0.8^{|i-j|}$. Here, we set $\sigma^{2}$ with $R^{2}=0.7$. All other aspects remain the same as in Scenario 1. For the missing rate, two settings are considered, where

• •

  $(a,b,c)=(0.1,-2,-1)$, so that the missing rate is approximately 60%,

• •

  $(a,b,c)=(0.7,-3.5,-4)$, so that the missing rate is approximately 90%.

The NAD and MSE results are shown in Figures 14, 15, 17, 18, 20, 21, 23, and 24. The optimal MSE rates are shown in Figures 16, 19, 22, and 25. The main conclusions are as follows:

1. 1.

   Scenario 2 results yield conclusions similar to those of Scenario 1. PRIME produces the smallest NADs and MSEs in almost all cases. As shown in Figures 21 and 24, although CC has the smallest estimation error, its estimation variance is extremely high, which causes problems in the MSE.

2. 2.

   The optimal rate results show that PRIME has obvious advantages over other methods because it produces the smallest MSE in almost all situations except when $\text{MR}=90\%,n=200,$ and $\rho_{ij}=0.8$. However, when $\rho_{ij}$ increases, the gap between PRIME and other methods increases. CC still has the worst MSE performance among the four methods.

In this scenario, we illustrate the proposed SPRIME by studying the data from simulations. We consider penalized Full (denoted as SFull), complete-case analysis with a penalty (denoted as SCC), ILSE, and ML as the alternatives. ILSE in Lin et al. [9] and ML in Jiang et al. [4] were proposed without considering the sparse assumption; hence, we use them directly instead of using the penalized estimation form. We acknowledge that there are other approaches such as those in Xue and Qu [21] that can be used to address a high-dimensional missing-data problem. However, the missing-data patterns in these methods are different from the individual-specific case.

The model used to generate data has the linear expression

where $n=200$, $p=30$, $(\beta_{1},\cdots,\beta_{12})^{\top}=(1,-0.6,1.5,1,1.2,0.4,-1,-0.7,1.3,0.5,1.1,-1.4,0.9)$, and $\beta_{j}=0\ (j=13,\cdots,30)$. We generate $(X_{i1},\cdots,X_{ip})$ from the multivariate normal distribution $N_{p}(\bf{0},\Sigma)$. We set the non-diagonal element $\rho_{ij}$ of $\Sigma$ equal to $0.5$. For $\varepsilon_{i}$, we use the error distribution $N(0,\sigma^{2})$, where $\sigma^{2}$ changes with $R^{2}=\text{Var}({\bm{X}}_{i}^{\top}{\bm{\beta}})/\{\text{Var}({\bm{X}}_{i}^{\top}{\bm{\beta}})+\sigma^{2}\}$. We consider only the case in which $R^{2}=0.7$.

When taking sparse structure into consideration, the coefficients $\beta_{j}\ (j=13,\cdots,30)$ are equal to 0. Thus, the criterion, NAD used in Scenario 1 and 2, is no longer meaningful. So we only use MSE and the optimal rate of MSE to assess the performance of different methods. Several conclusions can be drawn from the figures:

1. 1.

   Scenario 3 results yield conclusions similar to those of Scenarios 1 and 2. From Figure 26 we can see that SPRIME produces the smallest MSEs in almost all cases.

2. 2.

   When MR=60%, the optimal rates of MSE are 0.95, 0.01, 0.00, 0.04 for SPRIME, ILSE, CC and ML, respectively. When MR=90%, the optimal rates of MSE are 0.94, 0.05, 0.01, 0.00 for SPRIME, ILSE, CC and ML, respectively. Not surprisingly, SPRIME that considers the sparse structure performs better than ILSE and ML without the penalty, which reconfirms the superiority of PRIME-type approach.