Debiased Lasso After Sample Splitting for Estimation and Inference in High Dimensional Generalized Linear Models

Recent technological advances have allowed biomedical researchers to collect an extraordinary amount of genetic data from patients. Studies that seek to characterize associations between these genetic factors and medical outcomes face numerous challenges.

Firstly, the standard maximum likelihood estimator (MLE) is not well-defined for classical statistical models such as linear and generalized linear models (GLMs) when the covariates are high dimensional, in particular, $p\gg n$. Several variable selection and estimation procedures are available for sparse models with high dimensional covariates. The commonly-used lasso estimator (Tibshirani,, 1996) introduces a penalty that shrinks many coefficient estimates to exactly zero, thus performing simultaneous variable selection and coefficient estimation. Similar approaches include the adaptive lasso (Zou,, 2006), the smoothly clipped absolute deviation (SCAD) estimator (Fan and Li,, 2001) and elastic net regularization (Zou and Hastie,, 2005). Another common approach, which only performs variable selection, is sure independence screening (Fan and Lv,, 2008; Fan and Song,, 2010).

Secondly, although penalized regression methods can produce coefficient estimates for GLMs, they typically have substantial biases and do not quantify the uncertainty around these estimates. Thus they cannot be used to directly construct confidence intervals or to do hypothesis testing in the high dimensional setting. More recent progress on statistical inference for high dimensional data includes debiased lasso estimators (Zhang and Zhang,, 2014; van de Geer et al.,, 2014), MLE-based sample splitting approaches (Fei and Li,, 2021; Fei et al.,, 2019), and the decorrelated score test (Ning and Liu,, 2017). Other work on hypothesis testing without coefficient estimation includes the sample splitting approach for p-values by Meinshausen et al., (2009) and the test developed by Zhu and Bradic, (2018) for linear models without imposing sparsity assumptions on the regression parameters.

Lastly, we note that some post-model selection inference procedures, such as the simultaneous inference method of Kuchibhotla et al., (2020) and the sample splitting and bootstrap method of Rinaldo et al., (2019), are applicable to high dimensional covariates, although these in particular only address linear models. The goal of such methods is distinct from ours. They are designed to be valid for arbitrary model selection procedures and under potential model misspecification. This robustness comes at a cost of conservative inference results when a model selection procedure, such as the lasso estimator, can select a superset of the correct model with high probability, which is the setting that we consider below.

Our proposed method is based on sample splitting, a two-stage approach where, first, one part of the data is used to select the covariates to be included in the model, and then the other part is used to estimate the now lower-dimensional model ($p<n$ but $p$ may grow with $n$). Similar methods have been recently developed for linear models (Fei et al.,, 2019; Wang et al.,, 2020) and generalized linear models (Fei and Li,, 2021) based on the MLE, as well as for Cox proportional hazards models (Zhang et al.,, 2022) based on the decorrelated score function. Compared to Fei and Li, (2021), we use the lasso instead of sure independence screening for model selection, which allows for more mild theoretical assumptions, as well as the debiased lasso in place of the MLE, which can substantially improve the bias, variance, and confidence interval coverage of the resulting estimators, as shown in our simulations below. Note that debiased lasso estimation would be equivalent to using the MLE in sample splitting approaches for linear models. Our contributions are as follows.

For the model selection component, following Huang and Zhang, (2012) we show that, under mild sufficient conditions, the lasso screens out a large portion of the noise variables and selects a model containing the true signal variables with high probability in the case of random design. Existing work on the selected model size and estimation error of the lasso method for linear regression includes Bickel et al., (2009), Zhang and Huang, (2008), Meinshausen and Yu, (2009), and Zhao and Yu, (2006), under differing assumptions on the covariates. Similar results were shown for generalized linear models by van de Geer, (2008) and Huang and Zhang, (2012), only the latter of which addressed model selection consistency and the number of selected noise variables. Their results imply that the lasso selected model size is of the same order as the number of nonzero regression coefficients, with high probability. This is essential because if the selected model size is too large then its design matrix may not have full rank, even asymptotically, and thus the regression coefficients will still not be estimable via MLE or even the refined debiased lasso of Xia et al., (2021).

For the lower-dimensional model estimation stage, while a naive approach would be to use the standard maximum likelihood estimator with the selected covariates, which is commonly used in sample splitting literature, we instead apply the recently developed refined debiased lasso approach of Xia et al., (2021). This method is better suited for models that contain a substantial number of noise variables, as is typically expected after variable selection. This is a more desirable approach because the conditions for the GLM lasso to not select any noise variables, for example, are known to be quite stringent (Huang and Zhang,, 2012). We illustrate the potentially large difference in performance through simulations, where the MLE exhibits a strong bias that inflates the size of estimated coefficients, and this bias increases with the true signal strength, in contrast to the approximately unbiased estimates from the debiased lasso. Such tendencies are discussed further by Xia et al., (2021) in the lower-dimensional setting.

For a set of prespecified coefficients, since their estimators based on a single sample split suffer from a loss of efficiency due to only using part of the data for estimation, we further investigate the idea of averaging estimates across multiple sample splits. Fei et al., (2019) proposed a multiple splitting procedure for linear models, where they required model selection consistency for their asymptotic theory. Wang et al., (2020) discussed multiple splitting for a single treatment variable in linear models, under much milder assumptions. Recently, Fei and Li, (2021) proposed a multiple splitting procedure for generalized linear models based on the MLE under a partial orthogonality condition, which requires that the signal variables be independent of the noise variables. For our proposed multiple splitting procedure, we apply the debiased lasso estimator instead of the MLE in the estimation stage, and show through simulations that our procedure results in approximately unbiased estimates and substantially reduced variability compared to the MLE-based multiple splitting that is often biased. For the theoretical development, we adapt the mild conditions of Wang et al., (2020) to GLMs and show the asymptotic normality of our proposed approach. As evidenced by simulations, our multiple splitting estimator can produce confidence intervals with the nominal coverage.

We first provide brief overviews of lasso and debiased lasso methods, then introduce our proposed sample splitting debiased lasso methods for GLMs.

For the single split estimator, we make the following assumptions:

1. 1.

   The covariates $\boldsymbol{x}_{i}$ are sub-Gaussian random vectors and $\left\|\boldsymbol{X}\right\|_{\infty}\leq K$ almost surely for some constant $K>0$.

2. 2.

   $\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0}}$ and $E[\boldsymbol{X}^{T}\boldsymbol{X}/n]$ are positive definite with eigenvalues bounded from above and away from zero.

3. 3.

   The derivatives $\dot{\rho}(y,a)=\partial\rho(y,a)/\partial a$ and $\ddot{\rho}(y,a)=\partial^{2}\rho(y,a)/\partial a^{2}$ exist for all $(y,a)$ and there exists a $\delta>0$ and constant $c_{Lip}>0$ such that $\ddot{\rho}$ is locally Lipschitz

   $$\max_{a_{0}\in\{\boldsymbol{x}_{i}^{T}\boldsymbol{\beta}^{0}\}}\sup_{\max\left(|a-a_{0}|,|\widehat{a}-a_{0}|\right)\leq\delta}\sup_{y\in\mathcal{Y}}\frac{|\ddot{\rho}(y,a)-\ddot{\rho}(y,\widehat{a})|}{|a-\widehat{a}|}\leq c_{Lip}.$$

   Also, there exist constants $K_{1},K_{2}>0$ such that the derivatives are bounded:

   $$\max_{a_{0}\in\{\boldsymbol{x}_{i}^{T}\boldsymbol{\beta}^{0}\}}\sup_{y\in\mathcal{Y}}|\dot{\rho}(y,a_{0})|\leq K_{1},$$

   $$\max_{a_{0}\in\{\boldsymbol{x}_{i}^{T}\boldsymbol{\beta}^{0}\}}\sup_{y\in\mathcal{Y}}|\ddot{\rho}(y,a_{0})|\leq K_{2}.$$

4. 4.

   $\left\|\boldsymbol{X}\boldsymbol{\beta}^{0}\right\|_{\infty}$ is bounded above almost surely.

5. 5.

   The sparsity $s_{0}$ satisfies $s_{0}\log(p)/\sqrt{n}=o(1)$ and the regression parameters are large enough that $\min_{j\in S_{0}}|\boldsymbol{\beta}_{j}^{0}|\geq\mathcal{O}\left(s_{0}\sqrt{\log(p)/n}\right)$ as $n\rightarrow\infty$.

Assumptions 1, 2, and 4 are typical in high dimensional literature. Bounded covariates are common in practice, including dummy variables for categorical covariates, minor allele counts, and physical measurements. Assumption 3 is required in order to apply the results for the refined debiased lasso estimator of Xia et al., (2021). All but the boundedness condition on $\dot{\rho}$ are satisfied when the function $A$ is twice continuously differentiable, as is the case with commonly-used GLMs, due to assumption 4. The derivative $\dot{\rho}$ is also bounded if the response variable has bounded support, such as in logistic regression, and this condition can be relaxed to include other common GLMs by instead assuming sub-exponential tails, as in Lemma 3.4 of Ning and Liu, (2017). Assumption 5 guarantees that, with high probability, model selection based on the lasso does not miss any signal variables, so that the selected model is not misspecified. This is a standard condition for sample splitting procedures that is not required by the desparsified lasso (van de Geer et al.,, 2014).

For multiple splitting, we make the following additional assumptions for the set of target covariates $S$, where $\boldsymbol{a}$ is a $(p+1)$-vector such that $\left\|\boldsymbol{a}_{S}\right\|_{2}=1$ and $\left\|\boldsymbol{a}_{S^{c}}\right\|_{2}=0$. Let $\mathcal{Z}=(y_{i},\boldsymbol{x}_{i})_{i=1}^{n}$ denote the entire data set.

1. 6.

   For independent sampling indicator vectors $\boldsymbol{v}$ and $\boldsymbol{\widetilde{v}}$ with corresponding fitted models $\widehat{S}$ and $\widetilde{S}$, define

   $$h_{i,n}=\left[E\left(\boldsymbol{v}_{i}\boldsymbol{a}_{\widehat{S}}^{T}\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0},\widehat{S}}^{-1}\boldsymbol{\widetilde{I}}_{\widehat{S}}\Big{|}\mathcal{Z}\right)-E\left(\boldsymbol{v}_{i}\boldsymbol{a}_{\widetilde{S}}^{T}\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0},\widetilde{S}}^{-1}\boldsymbol{\widetilde{I}}_{\widetilde{S}}\Big{|}\mathcal{Z}\right)\right]\dot{\rho}_{\boldsymbol{\beta}^{0}}(y_{i},\boldsymbol{x}_{i})$$

   where $\boldsymbol{\widetilde{I}}_{S}$ is the $|S|\times(p+1)$ matrix such that $\boldsymbol{\widetilde{I}}_{S}\boldsymbol{b}=\boldsymbol{b}_{S}$ for any $(p+1)$-vector $\boldsymbol{b}$, and the expectations are with respect to the sampling weights, conditional on the data. We assume that $\sum_{i=1}^{n}h_{i,n}/\sqrt{n}=o(1)$ and that $P_{n}\ddot{\rho}_{\boldsymbol{\widetilde{\beta}}^{\lambda}}$ is invertible when computed from any estimation subsample.

2. 7.

   There exists a random $(p+1)$-vector $\boldsymbol{\eta}_{n}$ independent of the data such that $\left\|\boldsymbol{\eta}_{n}\right\|_{\infty}$ is bounded and

   $$\left\|E\left(\boldsymbol{a}_{\widehat{S}}^{T}\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0},\widehat{S}}^{-1}\boldsymbol{\widetilde{I}}_{\widehat{S}}\Big{|}\mathcal{Z}\right)-\boldsymbol{\eta}_{n}^{T}\right\|_{1}={o}_{p}\left(1/\sqrt{\log(p)}\right).$$

Similar conditions were used and discussed further by Wang et al., (2020) in the context of linear models. Since, for low-dimensional GLMs, the refined debiased lasso estimator is asymptotically equivalent to the MLE $\widetilde{\boldsymbol{\beta}}^{ML}$ in the sense that $\sqrt{n}(\boldsymbol{\widetilde{\beta}}^{DL}_{j}-\widetilde{\boldsymbol{\beta}}_{j}^{ML})={o}_{p}\left(1\right)$ for all $j$, our theoretical arguments also apply to multiple splitting with the MLE, using the lasso for model selection. Fei and Li, (2021) used a stronger partial orthogonality condition, where the signal variables are independent of the noise variables, to derive the asymptotic normality of their MLE-based multiple splitting estimator. We further discuss the motivation behind Assumptions 6 and 7 and their relationship to the assumptions used by other work on sample splitting procedures in Appendix D, but provide a brief overview here.

Under Assumptions 1-5 and the sparsity rates required in Theorem 3 below, it can be shown that

Therefore, in order to prove asymptotic normality we need to control the randomness of the vector $\boldsymbol{a}_{\widehat{S}}^{T}\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0},\widehat{S}}^{-1}\boldsymbol{\widetilde{I}}_{\widehat{S}}\boldsymbol{v}_{i}$ averaged across all sample splits that exclude the $i$th sample from model selection (i.e. $\boldsymbol{v}_{i}=1$). The ideal case, when this average is deterministic, can be achieved under model selection consistency, so that $\widehat{S}$ is fixed at $S\cup S_{0}$ for each sample split, or if the set of covariates that are always included in the fitted model are independent of all remaining covariates, essentially giving $\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0}}^{-1}$ a block diagonal structure. Either of these two conditions imply that our less restrictive assumptions hold, where we only require that the effect of always excluding a single sample from model selection is asymptotically negligible (Assumption 6) and that the average converges in probability to a bounded random vector with moderate rate (Assumption 7).

The variable screening properties and asymptotic normality of our estimators are presented in the following three theorems, with their proofs given in Appendix A through Appendix C, respectively. Note that any variable selection procedure with the same screening properties as the lasso may be used in our sample splitting procedures.

We apply the proposed estimating methods with sample splitting in several simulations settings with high dimensional covariates in order to assess the potential benefits of using the lasso for model selection and the debiased lasso in place of the MLE for making statistical inference, and of averaging estimates across multiple sample splits (MS) as opposed to using a single split (SS) estimator for a fixed set of coefficients. Simulation results are all from logistic regression models, where the benefits of using the debiased lasso are particularly apparent. This may be due to the numerical instability of the Hessian matrix of the negative log-likelihood $P_{n}\ddot{\rho}_{\boldsymbol{\beta}}(y_{i},\boldsymbol{x}_{i})=n^{-1}\sum_{i=1}^{n}\hat{p}(\boldsymbol{x}_{i}^{T}\boldsymbol{\beta})\left[1-\hat{p}(\boldsymbol{x}_{i}^{T}\boldsymbol{\beta})\right]\boldsymbol{x}_{i}\boldsymbol{x}_{i}^{T}$, where $\hat{p}(\cdot)$ is the cdf of the standard logistic distribution. Note $\hat{p}(\boldsymbol{x}_{i}^{T}\boldsymbol{\beta})$ will be close to zero or one for large coefficient values, which can result in a near-singular Hessian matrix. Therefore the debiased lasso approach of performing a single-step maximization of the log-likelihood starting from an initial lasso estimator that is biased towards zero can help alleviate this issue.

In all simulations, $n=500$ samples are generated from a logistic regression model with $p=700$ covariates. The covariates are generated from a $N(\boldsymbol{0},\boldsymbol{\Sigma})$ distribution before being truncated at $\pm 3$. The covariance matrix $\boldsymbol{\Sigma}$ has ones on the diagonal and either an AR(1) correlation structure $\boldsymbol{\Sigma}_{jk}=0.5^{|j-k|}$ or a compound symmetry structure $\boldsymbol{\Sigma}_{jk}=0.5$ $(j\neq k)$. Further simulations that demonstrate the robustness of our procedures under higher covariate correlations are presented in Appendix E. There are $s_{0}=6$ signal covariates, and their index set was randomly chosen. The corresponding coefficient values are $\boldsymbol{\beta}^{0}_{S_{0}}=(-1.5,-1,-0.5,0.5,1,1.5)^{T}$.

We consider all signal covariates together with two randomly chosen noise covariates. For each coefficient $\beta_{j}$, we assess the single and multiple splitting debiased lasso estimators fit on $\widehat{S}\cup\{j\}$, as well as the corresponding MLE-based estimators. The splitting proportion is $q=0.5$. We also provide results from the oracle MLE, fit on $S_{0}\cup\{j\}$, for reference, using either the entire sample or half of the sample. Model selection is performed using lasso estimators with the values for $\lambda$ chosen by 10-fold cross validation. The multiple splitting estimators use $B=1,000$ splits. We use the R package glmnet for lasso estimation and the built-in glm function for MLE and debiased lasso estimation, using the control and start arguments to specify a single iteration starting at an initial lasso estimate for the latter.

The simulation results for the AR(1) correlation structure are summarized in Table 1. The single split MLE exhibits a bias that greatly inflates the size of the estimated coefficients, and the magnitude of this bias increases with the signal strength. In contrast, the single split debiased lasso estimator is approximately unbiased. Averaging across multiple splits does not improve the bias of the MLE, which is also apparent to a lesser extent in the logistic regression simulations of Fei and Li, (2021) summarized in their Table 3. The MLE has substantially higher variability than the debiased lasso, even for the multiple splitting estimators. The 95% confidence interval coverage for noise variables is roughly the same across all considered methods. For signal variables, however, the MLE has poor coverage that actually worsens after multiple splitting. This issue appears to be particularly severe in logistic regression models and is mild in some other GLMs such as Poisson regression. In contrast, the debiased lasso after multiple splitting performs well in achieving the nominal 95% coverage for all considered coefficients. All single split standard errors tend to be underestimated for signal variables, leading to slight undercoverage, while the multiple splitting standard errors are approximately unbiased. Multiple splitting drastically lowers the variability of estimates from either the MLE or debiased lasso compared to a single split. This produces a dramatic improvement in rejection rate for small coefficients. Note that the rejection rates for the MLE estimators are inflated due to their bias, which partially offsets the wider confidence interval length. In summary, for each pre-specified coefficient, the multiple splitting debiased lasso estimator provides the best performance in terms of bias, variability, and confidence interval coverage in this simulation setting, where the correlation between covariates decays rather quickly.

The simulation results for compound symmetry correlation structure are summarized in Table 2. The same trends concerning the bias and large variability of the MLE seen in Table 1 are also present in this setting. Again the debiased lasso estimators are approximately unbiased, and the multiple splitting debiased lasso estimator has approximate 95% coverage for each coefficient. Multiple splitting again greatly reduces the variability of single split estimators, resulting in thinner confidence intervals with more power for detecting small coefficient values.

Lastly, we assess the performance of the post-model selection procedure that does not pre-specify any covariate of interest. This simulation setting is identical to that of Table 1 with AR(1) correlation structure, but the estimates are now all based on a single model fit on only the selected covariates. For covariates that are not selected, their coefficient estimate and standard error are both set to zero. The oracle MLE results we present are estimated on $S_{0}$ instead of $S_{0}\cup\{j\}$.

The simulation results are summarized in Table 3. For small coefficients, there is poor confidence interval coverage across all non-oracle estimators due to the randomness associated with their inclusion in each model. For larger coefficients that are nearly always selected by the lasso, the performance resembles that of Table 1. These results demonstrate the importance of only analyzing coefficients that are included in the fitted model. For larger sample sizes, the lasso selection performance can improve substantially, leading to improved coverage and bias correction of the debiased lasso after model selection. See Appendix E for additional simulations results with a larger sample size. We also provide additional simulation results with higher correlations among covariates in Appendix E.

We apply the proposed method to a dataset of single nucleotide polymorphisms (SNPs) from a sample of African-American participants in the Mid-South Tobacco Case-Control study population (Jiang et al.,, 2019; Xu et al.,, 2020; Han et al.,, 2022) to assess genetic risk factors for nicotine dependence. The dataset is publicly available with GEO accession GSE148375. It originally contained 242,901 SNPs measured across 3399 individuals. We excluded SNPs that were insertions or deletions, had a call rate less than 95%, were not in Hardy-Weinberg equilibrium ($p>10^{-6}$), or had a minor allele frequency of less than 0.01. Subjects missing more than 1% of the remaining SNPs were excluded, and missing values were then imputed using the observed allele frequencies for each SNP. After data cleaning, the covariates consist of 32,557 SNPs as well as gender and age. The response variable is a binary indicator of smoking status (1=smoker), where 1607 of the 3317 participants are smokers.

Prior research has identified several genetic regions with SNPs that are associated with nicotine dependence, including 15q25.1 (CHRNA5/A3/B4), 8p11.21 (CHRNB3/A6), and 19q13.2 (CYP2A6/A7/B6, EGLN2, RAB4B, and NUMBL). See, for example, Yang and Li, (2016) and the references therein. We choose the target covariates to be the 16 measured SNPs that lie in these three regions, as well as the demographic variables.

The results from our multiple splitting debiased lasso estimator are presented in Table 4. After adjusting for multiple testing using the Holm procedure (Holm,, 1979), none of the SNPs have a significant association with smoking status. The SNP rs3733829 has the largest coefficient estimate, with an estimated 37% increase in the odds of being a smoker and unadjusted 95% confidence interval (9%, 73%), controlling for all other covariates. The association of rs3733829 with increased cigarette consumption has been discussed by Kim et al., (2010) and Bloom et al., (2014).

This work was supported in part by NIH Grants R01AG056764 and RF1AG075107, and by NSF Grant DMS-1915711.