Sure independence screening in the presence of missing data

Abstract

Variable selection in ultra-high dimensional data sets is an increasingly prevalent issue with the readily available data arising from, for example, genome-wide associations studies or gene expression data. When the dimension of the feature space is exponentially larger than the sample size, it is desirable to screen out unimportant predictors in order to bring the dimension down to a moderate scale. In this paper we consider the case when observations of the predictors are missing at random. We propose performing screening using the marginal linear correlation coefficient between each predictor and the response variable accounting for the missing data using maximum likelihood estimation. This method is shown to have the sure screening property. Moreover, a novel method of screening that uses additional predictors when estimating the correlation coefficient is proposed. Simulations show that simply performing screening using pairwise complete observations is out-performed by both the proposed methods and is not recommended. Finally, the proposed methods are applied to a gene expression study on prostate cancer.

Appendix

Proof of Theorem 2.1

Proof

Recall that the log-likelihood of ${\varvec{\phi }}_j$ is

$$\begin{aligned}&\ell (\phi _j| \{(X_i, Y_i)\}_{i = 1}^{n_j}, \{Y_k\}_{k = n_j+1}^{n})\\&\quad = -\frac{1}{2\sigma _{j \cdot y}^2}\sum _{i=1}^{n_j}(X_{ij} - \mu _{j \cdot y} - \beta _{jy}Y_i)^2 - \frac{n_j\log (\sigma _{j \cdot y}^2)}{2}\\&\qquad -\frac{1}{2\sigma _{y}^2}\sum _{i=1}^n(Y_i - \mu _{y})^2 - \frac{n\log (\sigma _{y}^2)}{2}. \end{aligned}$$

The inverted hessian of the log-likelihood evaluated at the estimated parameters is

$$\begin{aligned} H^{-1}_{{\varvec{\phi }}_j}|_{{\hat{{\varvec{\phi }}}}_j} = \begin{bmatrix} {\hat{\sigma }}_{y}^2/n&0&0&0&0 \\ 0&2{\hat{\sigma }}_{y}^4/n&0&0&0 \\ 0&0&{\hat{\sigma }}_{j \cdot y}^2(1 + {\bar{y}}^2/s_{y}^2)/n_j&-{\bar{y}}{\hat{\sigma }}_{j \cdot y}^2/(n_js_{y}^2)&0 \\ 0&0&-{\bar{y}}{\hat{\sigma }}_{j \cdot y}^2/(n_js_{y}^2)&{\hat{\sigma }}_{j \cdot y}^2/(n_js_{y}^2)&0\\ 0&0&0&0&2{\hat{\sigma }}_{j \cdot y}^4/n_j \end{bmatrix}, \end{aligned}$$

so that the large sample covariance matrix for ${\varvec{\theta }}_j$ can be written as $D(\rho _j)H^{-1}_{{\varvec{\phi }}_j}|_{{\hat{{\varvec{\phi }}}}_j}D(\rho _j)^T$, where $D(\rho _j) = \left( \frac{\partial \rho _j}{\partial \mu _y}, \frac{\partial \rho _j}{\partial \sigma _y^2}, \frac{\partial \rho _j}{\partial \mu _{j \cdot y}}, \frac{\partial \rho _j}{\partial \beta _{j \cdot y}}, \frac{\partial \rho _j}{\partial \sigma _{j \cdot y}^2}\right)$. It can be shown that

$$\begin{aligned}&D(\rho _j) {=} \left( 0, \frac{\sigma _{j \cdot y}^2\beta _{j \cdot y}}{2\sqrt{\sigma _{y}^2}(\beta _{j \cdot y}^2\sigma _{y}^2 + \sigma _{j \cdot y}^2)^{3/2}}, 0, \frac{\sigma _{j \cdot y}^2\sqrt{\sigma _{y}^2}}{(\beta _{j \cdot y}^2\sigma _{y}^2 + \sigma _{j \cdot y}^2)^{3/2}}, -\frac{\beta _{j \cdot y}\sqrt{\sigma _{y}^2}}{2(\beta _{j \cdot y}^2\sigma _{y}^2 + \sigma _{j \cdot y}^2)^{3/2}}\right) , \end{aligned}$$

and hence one finds

$$\begin{aligned}&D(\rho _j)H^{-1}_{{\varvec{\phi }}_j}|_{{\hat{{\varvec{\phi }}}}_j}D(\rho _j)^T \\&\quad = \frac{{{\hat{\sigma }}}_{j \cdot y}^4{{\hat{\beta }}}_{j \cdot y}^2}{4{{\hat{\sigma }}}_{y}^2({{\hat{\beta }}}_{j \cdot y}^2{{\hat{\sigma }}}_{y}^2 + {{\hat{\sigma }}}_{j \cdot y}^2)^{3}}\frac{2{{\hat{\sigma }}}_y^4}{n} + \frac{{{\hat{\sigma }}}_{j \cdot y}^4{{\hat{\sigma }}}_{y}^2}{({{\hat{\beta }}}_{j \cdot y}^2{{\hat{\sigma }}}_{y}^2 + {{\hat{\sigma }}}_{j \cdot y}^2)^{3}}\frac{{{\hat{\sigma }}}_{j \cdot y}^2}{n_js_{y}^2} \\&\qquad + \frac{{{\hat{\beta }}}_{j \cdot y}^2{{\hat{\sigma }}}_{y}^2}{4({{\hat{\beta }}}_{j \cdot y}^2{{\hat{\sigma }}}_{y}^2 + {{\hat{\sigma }}}_{j \cdot y}^2)^{3}}\frac{2{{\hat{\sigma }}}_{j \cdot y}^4}{n_j}\\&\quad = \Big [2{{\hat{\sigma }}}_y^4s_y^2n_j(s_j^2-s_{jy}^2/s_y^2)^2s_{jy}^2/s_y^4 + 4n{{\hat{\sigma }}}_y^4(s_j^2-s_{jy}^2/s_y^2)^3\\&\qquad + 2ns_y^2{{\hat{\sigma }}}_y^4(s_j^2-s_{jy}^2/s_y^2)^2s_{jy}^2/s_y^4\Big ]/\left[ 4nn_j{{\hat{\sigma }}}_y^2s_y^2({{\hat{\sigma }}}_y^2s_{jy}^2/s_y^4+s_j^2-s_{jy}^2/s_y^2)^3\right] \\&\quad = \frac{(1-{{\tilde{\rho }}})^2{{\hat{\sigma }}}_y^4\left[ 2s_y^2n_j{{\tilde{\rho }}}^2s_j^6/s_y^2 + 4ns_j^6(1-{{\tilde{\rho }}}^2) + 2ns_y^2s_j^6{{\tilde{\rho }}}^2/s_y^2\right] }{4nn_js_y^2{{\hat{\sigma }}}_y^2s_j^6({{\tilde{\rho }}}^2({{\hat{\sigma }}}_y^2/s_y^2)+1)^3}\\&\quad = (1 - {\tilde{\rho }}_j^2)^2\left( \frac{{\hat{\sigma }}_y^2}{s_y^2}\right) \left( \frac{1}{nn_j}\right) \left( \frac{{\tilde{\rho }}_j^2(n_j-n)/2 + n}{\left( {\tilde{\rho }}_j^2(\frac{{\hat{\sigma }}_y^2}{s_y^2} - 1) + 1\right) ^3}\right) . \end{aligned}$$

Since ${\hat{\rho }}_j$ is the maximum likelihood estimator computed from a Normal distribution, it follows that $[D(\rho _j)H^{-1}_{{\varvec{\phi }}}|_{{\hat{{\varvec{\phi }}}}_j}D(\rho _j)^T]^{-1/2}({\hat{\rho }}_j - \rho _j)$ converges to a standard Normal distribution. $\square$

Proof of Theorem 2.2

Proof

First note that the estimated covariance $s_{jy}$ based on the completely observed pairs is, except for a scale of $(n_j-1)/(n_j)$, a U-statistic (Kowalski and Tu 2007)

$$\begin{aligned} s_{jy}= & {} \frac{1}{n_j}\sum _{i=1}^{n_{j}}(Y_{i} - {\bar{Y}})(X_{ij} - {\bar{X}}_j) = \frac{n_j -1}{n_j}{n_j \atopwithdelims ()2}^{-1}\sum _{i\ne k}^{n_j}\frac{1}{2}(Y_i - Y_k)(X_{ij} - X_{kj})\\= & {} \frac{n_j -1}{n_j} \frac{1}{(n_j)(n_j-1)}\sum _{i\ne k}^{n_j} h_j(Y_i, Y_k, X_{ij}, X_{kj}) := \frac{n_j -1}{n_j}s_{jy}^*, \end{aligned}$$

where ${\bar{X}}_j = \sum _{i=1}^nX_{ij}$ and $h_j(Y_i, Y_k, X_{ij}, X_{kj}) = (Y_i - Y_k)(X_{ij} - X_{kj})$ is the kernel of the U-statistic $s_{jy}^*$. Note that $E(s_{jy}^*) = \sigma _{jy} := \sigma _j\sigma _y\rho _j$.

We follow steps similar to those in Li et al. (2012a). First write

$$\begin{aligned}&s_{jy}^* = s_{jy,1}^{*} + s_{jy,2}^{*}\\&\quad := \frac{1}{n_j(n_j-1)}\sum _{i\ne k}^{n_j} h_j(Y_i, Y_k, X_{ij}, X_{kj})I(h_j(Y_i, Y_k, X_{ij}, X_{kj}) \le M) \\&\qquad + \frac{1}{n_j(n_j-1)}\sum _{i\ne k}^{n_j} h_j(Y_i, Y_k, X_{ij}, X_{kj})I(h_j(Y_i, Y_k, X_{ij}, X_{kj}) > M), \end{aligned}$$

and define

$$\begin{aligned} \sigma _{jy,1}:= & {} E(s_{jy,1}^{*}) = E[ h_j(Y_i, Y_k, X_{ij}, X_{kj})I(h_j(Y_i, Y_k, X_{ij}, X_{kj}) \le M)], \\ \sigma _{jy,2}:= & {} E(s_{jy,2}^{*}) = E[h_j(Y_i, Y_k, X_{ij}, X_{kj})I(h_j(Y_i, Y_k, X_{ij}, X_{kj}) > M)]. \end{aligned}$$

Because $s_{jy,1}^{*}$ can be written as an average of averages of i.i.d. random variables (Serfling 1980 - sec. 5.1.6), for any $t > 0$ and $\epsilon > 0$ we have

$$\begin{aligned}&P(s_{jy,1}^{*} - \sigma _{yj,1} \ge \epsilon ) \\&\quad \le \exp (-t\epsilon )\exp (-t\sigma _{jy,1})E(\exp (ts_{jy,1}^{*}))\\&\quad = \exp (-t\epsilon )\exp (-t\sigma _{jy,1}) E\left( \exp \left( t \frac{1}{n_j!}\sum _{n_j!}\frac{1}{m}\sum _{m}h_j^{(m)}I(h_j^{(m)} \le M)\right) \right) \\&\quad \le \exp (-t\epsilon )\exp (-t\sigma _{jy,1})E^m\left( \exp \left( \frac{1}{m}th_j^{(m)}I(h_j^{(m)} \le M)\right) \right) \\&\quad = \exp (-t\epsilon )E^m\left( \exp \left( \frac{1}{m}t\left( h_j^{(m)}I(h_j^{(m)} \le M) - \sigma _{yj,1}\right) \right) \right) , \end{aligned}$$

where $m = [n_j/2]$ and the last inequality follows from Theorem 5.6.1A in Serfling (1980). Choose $t = 4\epsilon m/M^2$ so that $P(s_{jy,1}^{*} - \sigma _{jy,1} \ge \epsilon ) \le \exp (-2\epsilon ^2 m/M^2)$ and by symmetry of the U-statistics

$$\begin{aligned} P(|s_{jy,1}^{*} - \sigma _{jy,1}| \ge \epsilon ) \le 2\exp (-2\epsilon ^2 m/M^2). \end{aligned}$$

(3)

Now we deal with $s_{jy,2}^*$. Note that using Cauchy–Schwarz and Markov inequalities we have

$$\begin{aligned} \sigma _{jy, 2}^2\le & {} E[(Y_i - Y_k)^2(X_{ij} - X_{kj})^2]P[(Y_i - Y_k)(X_{ij} - X_{kj}) \ge M]\\\le & {} E[(Y_i - Y_k)^2(X_{ij} - X_{kj})^2]E[\exp (s(Y_i - Y_k)(X_{ij} - X_{kj}))]\exp (-sM) \end{aligned}$$

for any $s > 0$. Using assumptions C1, if we choose $M = cn_j^\gamma$ for $0< \gamma < 1/2 - k$, then $\sigma _{jy,2} \le \epsilon /2$ when $n_j$ is sufficiently large. Consequently,

$$\begin{aligned}&P(|s_{jy,2}^{*} - \sigma _{jy,2}|> \epsilon )\\&\quad \le P(|s_{jy,2}^{*}|> \epsilon /2) \le P(\cup \{(Y_i - Y_k)(X_{ij} - X_{kj})> M\}\\&\quad \le n_jP((Y_i - Y_k)(X_{ij} - X_{kj})> M)\\&\quad = n_jP[\exp (s(Y_i - Y_k)(X_{ij} - X_{kj})) > \exp (sM)]\\&\quad \le n_j\exp (-sM)E(\exp \{s(Y_i - Y_k)(X_{ij} - X_{kj})\}) = n_jC\exp (-sM), \end{aligned}$$

for any $s > 0$. Hence

$$\begin{aligned}&P(|s_{jy}^{*} - \sigma _{jy}|> 2\epsilon )\nonumber \\&\quad = P(|s_{jy,1}^{*} + s_{jy,2}^{*} - \sigma _{jy,1} - \sigma _{jy,2}| \ge 2\epsilon )\nonumber \\&\quad \le P(|s_{jy,1}^{*} - \sigma _{jy,1}|> \epsilon ) + P(|s_{jy,2}^{*} - \sigma _{jy,2}| > \epsilon )\nonumber \\&\quad \le O(\exp (-c_1\epsilon ^2n_j^{1-2\gamma }) + n_j\exp (-c_2n_j^\gamma )). \end{aligned}$$

(4)

Recall ${{\hat{\rho }}}_j = \frac{s_{jy}}{s_{j}s_y}\frac{{\hat{\sigma }}_y}{s_{y}}\frac{s_j}{{{\hat{\sigma }}}_j}$. Using similar arguments, one can show that the convergence rate of $s_{y}, s_j, {\hat{\sigma }}_y$ and ${\hat{\sigma }}_j$ have the same form of (4) and hence by Lemma S4 in Liu et al. (2014) so does ${{\hat{\rho }}}_j$, so that we have

$$\begin{aligned}&P(|{\hat{\rho }}_j - \rho _{j} | \ge cn_j^{-\kappa })\\&\quad \le P(|{\hat{\rho }}_j - \rho _{j}| \ge c n_j^{-\kappa })\\&\quad = O([\exp (-c_1n_j^{1-2(\gamma +\kappa )}) + n_j\exp (-c_2n_j^\gamma )]). \\&P(|{\hat{\rho }}_j - \rho _{j} | \ge cn_j^{-\kappa }, \text { for all } j)\\&\quad \le \sum _{j = 1}^d P(|{\hat{\rho }}_j - \rho _{j}| \ge c n_j^{-\kappa })\\&\quad = \sum _{j = 1}^dO([\exp (-c_1n_j^{1-2(\gamma +\kappa )}) + n_j\exp (-c_2n_j^\gamma )]). \end{aligned}$$

Letting $\epsilon = cn_j^{-\kappa }$ we have

$$\begin{aligned}&P(\max _{j = 1, \ldots , d}|{\hat{\rho }}_j - \rho _{j} | \ge c n_j^{-\kappa })\\&\quad \le d \max _{j = 1, \ldots , d} P(|{\hat{\rho }}_j - \rho _{j}| \ge cn_j^{-\kappa })\\&\quad = d\max _{j = 1, \ldots , d}O(\exp (-c_1 n_j^{-2\kappa }n_j^{1-2\gamma }) + n_j\exp (-c_2n_j^\gamma ))\\&\quad \le O(d\exp (-c_1\min _jn_j^{1-2(\gamma +\kappa )}) + \max _j\{n_j\exp (-c_2n_j^\gamma )\}). \end{aligned}$$

If ${\mathcal {A}} \not \subseteq \hat{{\mathcal {A}}}$, then there exists a $j \in {\mathcal {A}}$ such that ${\hat{\rho }}_j < cn_j^{-\kappa }$. From condition C2 it follows that $|{\hat{\rho }}_j - \rho _j| > cn_j^{-\kappa }$ for some $j \in {\mathcal {A}}$. This implies that $\{{\mathcal {A}} \not \subseteq \hat{{\mathcal {A}}}\} \subseteq \{|{\hat{\rho }}_j - \rho _j| > cn_j^{-\kappa } \text { for some } j \in {\mathcal {A}}\}$. Then

$$\begin{aligned}&P({\mathcal {A}} \subseteq \hat{{\mathcal {A}}})\\&\quad \ge P(|{\hat{\rho }}_j - \rho _{j}| \le cn_j^{-\kappa }, \text { for all } j \in {\mathcal {A}}) \\&\quad = 1 - P(|{\hat{\rho }}_j - \rho _{j}|> cn_j^{-\kappa }, \text { for some } j \in {\mathcal {A}})\\&\quad \ge 1 - \sum _{j \in {\mathcal {A}}}P(|{\hat{\rho }}_j - \rho _{j}| > cn_j^{-\kappa })\\&\quad = 1 - \sum _{j \in {\mathcal {A}}}O(\exp (-c_1\min _jn_j^{1-2(\gamma +\kappa )}) + \max _j\{n_j\exp (-c_2n_j^\gamma )\}). \end{aligned}$$

$\square$