Controlled Variable Selection from Summary Statistics Only?
A Solution via GhostKnockoffs and Penalized Regression

Modern large-scale studies frequently involve a multitude of explanatory variables potentially associated with an outcome we would like to better understand. Oftentimes, the goal is to select those explanatory variables that are meaningfully associated with the response variable. For instance, with recent advances in genome sequencing technologies and genotype imputation techniques, one can now gather tens of millions of variants from hundreds of thousands of samples in large-scale genetic studies, with the aim of pinpointing which genetic variants are biologically associated with specific diseases. This information could provide mechanistic insights and potentially aid the development of targeted drugs. In statistics, this challenge is typically framed as a multiple testing problem. Further, due to the sheer number of hypotheses considered and the cost of following false leads, it is generally required to control some form of error rate on the false positives.

In this paper, we focus on controlling the false discovery rate (FDR), which is the expected proportion of false selections among all selected variables. Compared to the more stringent familywise error rate (FWER) control, keeping the FDR under a nominal level allows for more discoveries while maintaining a reasonable statistical guarantee on the rate of false positives. Several methods for FDR control have been proposed in the literature, with the Benjamini-Hochberg procedure being particularly popular (Benjamini and Hochberg, 1995). However, these approaches often assume a parametric model or the existence of valid $p$-values, which remains difficult, and even problematic, in high-dimensional settings.

Candès et al. (2018) proposed the model-X knockoffs, a broad and flexible framework which allows the statistician to select variables that retain dependence with the response conditional on all other covariates while maintaining FDR control. Model-X knockoffs differs from previous approaches in that (1) it makes no modeling assumptions on the distribution of the response $Y$ we wish to study conditional on the family of covariates $X$, and (2) it does not require the construction of valid $p$-values. Instead, the crucial assumption is that the distribution of $X$ is known. The main idea in Candès et al. (2018) is to generate fake variables $\widetilde{X}$, knockoffs, which we can view as negative controls and can be used to tease apart variables that do influence the response from those who do not. Model-X knockoffs has proved effective in a number of real-world applications, particularly in GWAS; see Bates et al. (2020), Sesia et al. (2021) and He et al. (2022) for examples.

To deploy model-X knockoffs, researchers must have in hand the covariates and responses from all samples. However, in certain situations, individual-level data that may reveal sensitive personal information is not readily accessible. For example, due to privacy concerns, many GWAS studies only publish summary statistics of the original data (Pasaniuc and Price, 2017). Yet in such cases, we would still like to develop controlled variable selection methods that rely solely on summary statistics. In genetic studies, this would enable us to utilize available summary data from different data centers to conduct meta-analysis, enhancing the effective sample size and improving variable selection power. On this front, He et al. (2022) proposed the framework of GhostKnockoffs, which implements the knockoffs procedure with the marginal correlation difference feature importance statistic directly from summary statistics. As we shall review next, the main idea is to generate knockoff $Z-$scores directly without creating knockoff variables; all that is needed are marginal correlations between the response and the features under study. In details, with $n$ being the sample size and $p$ the number of variables being assayed, the method operates with only $\mathbf{X}^{\top}\mathbf{Y}$ and $\lVert\mathbf{Y}\rVert_{2}^{2}$, where $\mathbf{X}$ is the $n\times p$ matrix of covariates, and $\mathbf{Y}$ is the $n\times 1$ response vector.

In this paper, we extend the family of GhostKnockoffs methods to incorporate feature importance statistics obtained from penalized regression. We first consider in Section 3 the situation in which the empirical covariance of the covariate-response pair $(X,Y)$ is available; with the above notation, this means that the summary statistics $\mathbf{X}^{\top}\mathbf{X}$, $\mathbf{X}^{\top}\mathbf{Y}$, $\lVert\mathbf{Y}\rVert_{2}^{2}$ are available along with the sample size $n$. Unsurprisingly, we observe substantial power improvement over the method of He et al. (2022) because we can now employ far more effective test statistics. Next, in Section 4, we consider the case where the empirical covariance $\mathbf{X}^{\top}\mathbf{X}$ of the features is not available. There, we propose new imputation methods that consistently outperform He et al. (2022) in comprehensive synthetic and semi-synthetic simulations and rigorously control the FDR under suitable conditions. Finally, in Section 5 we apply our methods to a meta-analysis of nine large-scale array-based genome-wide association and whole-exome/-genome sequencing studies of Alzheimer’s disease, in which our methods yield more discoveries than He et al. (2022). We note that existing work in the genetics literature has implemented variable selection methods based on penalized regression with summary statistics, e.g., Mak et al. (2017) and Zou et al. (2022). However, none of these provide any guarantee of FDR control. In fact, as we note in the main text, these methods can be leveraged in our approach to create knockoffs versions that do control the FDR.

The software and example code that reproduce the results presented in this paper can be found at https://github.com/biona001/ghostknockoff-gwas-reproducibility/tree/main/chen_et_al. Simulation results in Section 3.5, Section 4.4.2 and Section 4.4.3 can be exactly reproduced. Due to data accessibility issue, we only provide code without real data for Section 4.4.1 and Section 5.

Given covariates $X\in\mathbb{R}^{p}$ and a response $Y\in\mathbb{R}$, we are interested in understanding which variables influence $Y$. We formulate this selection problem as testing the conditional independence hypotheses $\mathcal{H}_{0}^{j}:X_{j}\perp\!\!\!\perp Y\mid X_{-j}$ for $1\leq j\leq p$, where $X_{-j}$ is a shorthand for all the variables except the $j$th; that is $X_{-j}=\{X_{1},...,X_{j-1},X_{j+1},...,X_{n}\}$. In words, we should reject $\mathcal{H}_{0}^{j}$ if we believe that $X_{j}$ can help better predict the outcome than if we only had available the values of all the other variables. Put differently, $X_{j}$ has information about $Y$ which cannot be subsumed by the information contained in all the other variables. By conditioning on $X_{-j}$, these hypothesis tests aim to weed out variables whose relationship to $Y$ is driven by residual correlations with other covariates.

Let $\mathcal{H}_{0}\subset[p]$ be the set of indices for which the null conditional independence hypothesis $\mathcal{H}_{0}^{j}$ is true, and let ${\mathcal{S}}\subset[p]$ be the set of indices of the hypotheses rejected by a selection procedure. The false discovery rate (FDR) is the expected fraction of false positives among the selected, defined as

with the convention that $0/0=0$. Our goal is to make as many rejections as possible while controlling the FDR below a user-specified level $q$.

In this paper, we consider the setting in which, instead of observing i.i.d. samples from the distribution of ($X,Y$), we only have some summary statistics of the i.i.d. samples. In particular, we will show how one can, quite remarkably, perform tests of conditional independence when we do not directly observe the i.i.d. samples. Throughout this paper, we assume that $X\sim\mathcal{N}(\mathbf{0},\mathbf{\Sigma})$ where $\mathbf{\Sigma}$ is known (or, in practice, can be estimated).

The original model-X knockoffs procedure relies on having access to the covariates and responses from all data points, i.e., the matrix of covariates $\mathbf{X}$ and the response vector $\mathbf{Y}$. Henceforth, we call these individual-level data. In many application scenarios, however, individual-level data are not available due to privacy concerns. Instead, we only have access to some summary statistics of $\mathbf{X}$ and $\mathbf{Y}$, e.g., the empirical covariance matrix of the covariaties and the empirical covariance between each covariate and the response.

He et al. (2022) proposed GhostKnockoffs, which implements the knockoffs procedure with marginal correlation difference statistic when only $\mathbf{X}^{\top}\mathbf{Y}$ and $||\mathbf{Y}||_{2}^{2}$ are available. The key idea of He et al. (2022) is to sample the knockoff $Z$-score $\widetilde{\mathbf{Z}}_{s}$ from $\mathbf{X}^{\top}\mathbf{Y}$ and $||\mathbf{Y}||_{2}^{2}$ directly, in a way such that

where $\widetilde{\mathbf{X}}=\mathcal{G}(\mathbf{X},\mathbf{\Sigma})$ is the knockoff matrix generated by the Gaussian knockoff sampler (Algorithm 1). If we use $\mathbf{W}=\mathbf{Z}_{s}-\widetilde{\mathbf{Z}}_{s}$ (where $\mathbf{Z}_{s}=\mathbf{X}^{\top}\mathbf{Y}$) as the feature importance statistic and run the knockoff filter, the resulting rejection set will have the same distribution as that of the knockoffs procedure with marginal correlation difference statistic. Therefore, the two procedures are statistically identical. In particular, they both control the FDR.

Specifically, He et al. (2022) showed that for $\mathbf{P}$ and $\mathbf{V}$ computed in step 3 of Algorithm 1,

satisfies (5) as detailed in Appendix B. All this is summarized in Algorithm 2. In the following sections, we refer to Algorithm 2 as GhostKnockoffs with marginal correlation difference statistic (GK-marginal).

As we have just seen, GhostKnockoffs-marginal gives a way to test conditional hypotheses while maintaining FDR control when only the summary statistics $\mathbf{\mathbf{X}}^{\top}\mathbf{Y}$ and $\lVert\mathbf{Y}\rVert_{2}^{2}$ are available to the analyst. Now, we consider the setting in which we have knowledge of the empirical covariance matrix $\mathbf{X}^{\top}\mathbf{X}$ and the sample size $n$, in addition to $\mathbf{\mathbf{X}}^{\top}\mathbf{Y}$ and $\lVert\mathbf{Y}\rVert_{2}^{2}$. These quantities only reveal sample averages of relevant quantities, as opposed to all the individual-level information.

In this section, we propose a variable selection method that utilizes only $\mathbf{X}^{\top}\mathbf{X}$, $\mathbf{\mathbf{X}}^{\top}\mathbf{Y}$, $\lVert\mathbf{Y}\rVert_{2}^{2}$, and $n$. Our method achieves FDR control and power comparable to the knockoffs procedure with the cross-validated Lasso coefficient difference statistic defined in Section 2. This is interesting because the latter usually outperforms GhostKnockoffs with the marginal correlation difference statistic by a significant margin. Notably, for a fixed tuning parameter $\lambda$, we show that our procedure is equivalent to a corresponding knockoffs method using the Lasso coefficient difference statistic with the same penalty level $\lambda$.

Recall that in the knockoffs procedure with the Lasso coefficient difference statistic, we solve the optimization problem

where $\widetilde{\mathbf{X}}=\mathcal{G}(\mathbf{X},\mathbf{\Sigma}).$ We then define the Lasso coefficient difference feature importance statistics by $W_{j}=\lvert\hat{\beta}_{j}(\lambda)\rvert-\lvert\hat{\beta}_{j+p}(\lambda)\rvert$ for $1\leq j\leq p$. If we have access to individual-level data, $\lambda$ is usually chosen by cross-validation (Candès et al. (2018) and Weinstein et al. (2020)).^***In the case that $Y$ is binary, one may think that utilizing (penalized) logistic regression would give much better power than Lasso. In Appendix P, we show that this intuition may not be correct through simulations, even when $Y$ is generated according to a logistic regression model.

As a first step, we would like to run a statistically equivalent procedure using $\mathbf{X}^{\top}\mathbf{X},\mathbf{X}^{\top}\mathbf{Y}$, $\lVert\mathbf{Y}\rVert_{2}^{2}$, and $n$ for a fixed $\lambda$. Note that, with $\lambda$ fixed, (7) depends on the data only through

and

Define the Gram matrix of $[\mathbf{X},\widetilde{\mathbf{X}},\mathbf{Y}]$

The Gram matrix can of course be equivalently reconstructed from $(\lVert\mathbf{Y}\rVert_{2}^{2},\mathbf{X}^{\top}\mathbf{Y},\widetilde{\mathbf{X}}^{\top}\mathbf{Y},\mathbf{X}^{\top}\mathbf{X},\widetilde{\mathbf{X}}^{\top}\mathbf{X},\widetilde{\mathbf{X}}^{\top}\widetilde{\mathbf{X}})$. The main idea is to sample from the joint distribution of $\mathcal{T}(\mathbf{X},\widetilde{\mathbf{X}},\mathbf{Y})$ using the Gram matrix of $[\mathbf{X},\mathbf{Y}]$ only. Based on this, we can then generate the solution to the Lasso problem (7) (in distribution) for a fixed $\lambda$.^***Careful readers may realize that the solution of the Lasso problem does not depend on $\lVert\mathbf{Y}\rVert_{2}^{2}$. Here we include $\lVert\mathbf{Y}\rVert_{2}^{2}$ as an input of to be able to make a more general statement later that goes beyond the Lasso. This is achieved via the following Proposition 1, which says in words that if we generate ‘fake’ data matrices $\widecheck{\mathbf{X}}$ and $\widecheck{\mathbf{Y}}$ that lead to the same Gram matrix as that of $\mathbf{X}$ and $\mathbf{Y}$, then the distribution of $\mathcal{T}$ remains unchanged if we replace the original data matrices by the fake data matrices.

Proof of Proposition 1 is provided in Appendix C. Specifically, Proposition 1 suggests that summary statistics ($\mathbf{X}^{\top}\mathbf{X},\mathbf{X}^{\top}\mathbf{Y},||\mathbf{Y}||_{2}^{2}$, $\mathbf{\Sigma}$) are sufficient for sampling the Gram matrix $\mathcal{T}(\mathbf{X},\widetilde{\mathbf{X}},\mathbf{Y})$.

We are now able to write down a procedure, namely, Algorithm 3, which is statistically equivalent to the corresponding individual-level knockoffs procedure using the Lasso coefficient difference statistic (or any statistic defined in Sections 3.3 and 3.4). In step 2, $\widecheck{\mathbf{X}}$ and $\widecheck{\mathbf{Y}}$ can be obtained by performing the eigen-decomposition or Cholesky decomposition of $[\mathbf{X}\ \mathbf{Y}]^{\top}[\mathbf{X}\ \mathbf{Y}]$. Brief procedures to construct $\widecheck{\mathbf{X}}$ and $\widecheck{\mathbf{Y}}$ via eigen-decomposition are provided in Appendix D. All we need to do is to run the knockoffs procedure with $\widecheck{\mathbf{X}}$ and $\widetilde{\widecheck{\mathbf{X}}}$ in lieu of ${\mathbf{X}}$ and $\widetilde{{\mathbf{X}}}$. We say that the procedure is equivalent since the rejection sets have the same distribution. In particular, this proves that Algorithm 3 controls the FDR.

We can easily adapt the method above to accommodate other types of regularization, such as Ridge regression and Elastic Net.

In Section 3.2, we assumed that the tuning parameter $\lambda$ in (7) is fixed. In practice, one may choose the penalty level using information from the Gram matrix of $[\mathbf{X},\mathbf{Y}]$, and the sample size $n$. Since individual-level data is not available, we are unable to use data-splitting approaches such as cross-validation.

An alternative way to define feature importance is to use the square-root Lasso (Belloni et al., 2011), for which the choice of a reasonable tuning parameter is convenient. The square-root Lasso applied to the knockoffs setting solves

and a good choice of $\lambda$ is given by

where $\bm{\epsilon}\sim\mathcal{N}(\mathbf{0},\mathbf{I}_{n})$ and $\kappa$ is a unitless hyperparameter (Tian et al., 2018). This value is a scalar multiple of the expected value of the minimal penalty level required such that all the coefficients are shrunk to zero under the global null model. The square-root Lasso has the benefit that the value of the hyperparameter does not depend on the details of the distribution of $Y$ conditional on $X$. We also found that the performance of our procedure does not depend very sensitively on the choice of $\kappa$. In our data examples, we take $\kappa=0.3$.

In the setting where we only know about values of the summary statistics, we simply replace ($\mathbf{X},\widetilde{\mathbf{X}},\mathbf{Y}$) by ($\widecheck{\mathbf{X}},\widetilde{\widecheck{\mathbf{X}}},\widecheck{\mathbf{Y}})$ in (8). Further, we note that for any orthogonal matrix $\mathbf{Q}$,

where the second equality follows from $\mathbf{Q}^{\top}\bm{\epsilon}\stackrel{{\scriptstyle d}}{{=}}\bm{\epsilon}$. Therefore, the value of the hyperparameter in (9) remains unchanged if we multiply $[\mathbf{X}\;\widetilde{\mathbf{X}}]$ by $\mathbf{Q}$ on the left. This implies that (9) is a deterministic function of $[\mathbf{X}\;\widetilde{\mathbf{X}}]^{\top}[\mathbf{X}\;\widetilde{\mathbf{X}}]$. Hence, the feature importance statistic is a function of $\mathcal{T}(\mathbf{X},\widetilde{\mathbf{X}},\mathbf{Y})$. Following Corollary 1, we can apply the knockoffs procedure with the square-root Lasso and matrices $(\widecheck{\mathbf{X}},\widetilde{\widecheck{\mathbf{X}}})$ in lieu of $({\mathbf{X}},\widetilde{{\mathbf{X}}})$. Upon choosing

we get a procedure, which is statistically indistinguishable from that we would get if we were performing all the same steps with $\mathbf{X}$ and $\widetilde{\mathbf{X}}$. (In practice, we compute the value in (10) via Monte Carlo simulation.) In the sequel, we call the resulting procedure summary statistics GhostKnockoffs with square-root Lasso importance statistic (GK-sqrtlasso). Note that GK-sqrtlasso controls the FDR as the flip-sign property of the feature importance statistic holds. This is because swapping a variable with its knockoff does not change the value of the hyperparameter. Therefore, by Corollary 1, applying the knockoff filter to the square-root Lasso feature importance statistics yields FDR control.

In the standard fixed-X knockoffs setting, cross-validation is also not feasible, since doing so would violate the sufficiency condition required for the feature importance statistics. As one possible alternative, Barber and Candès (2015) considered using as the feature importance statistic the value of $\lambda$ on the Lasso path at which feature $X_{j}$ first enters the model. Formally, they define the feature importance statistic

where $\hat{\bm{\beta}}(\lambda)$ is as in (7). We call this statistic the Lasso-max statistic. Intuitively, a larger penalty level is required to shrink an important feature to zero, so we should expect $W_{j}$ to be large and positive for non-nulls.

By Corollary 1, with the Lasso-max statistic Algorithm 3 produces a rejection set that has the same distribution as the rejection set obtained from the corresponding individual-data-based knockoffs procedure. We call this summary-statistic-based procedure GhostKnockoffs with Lasso-max statistic (GK-lassomax).

We remark that choices of other tuning parameters and feature importance statistics are also possible. For instance, we may choose $\lambda$ to minimize the Stein’s unbiased risk estimate (SURE) associated with (7). We shall however focus on the two approaches we have described.

We consider a variety of simulation settings in which we compare the performance of the proposed GhostKnockoffs with square-root Lasso and Lasso-max statistics (GK-sqrtlasso and GK-lassomax, defined in Sections 3.3 and 3.4), GhostKnockoffs with marginal correlation difference statistic (GK-marginal, defined in Section 2), and the knockoffs procedure with (cross-validated) Lasso coefficient difference statistic with individual-level data (KF-lassocv). Note that the first three are statistically equivalent to the corresponding knockoffs procedures with individual-level data.