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

TL;DR

This paper introduces a sample splitting approach using the debiased lasso for estimation and inference in high-dimensional generalized linear models, improving accuracy and efficiency over traditional methods.

Contribution

It proposes a novel multiple splitting procedure with the debiased lasso that achieves asymptotic normality and reduces bias and variance in high-dimensional settings.

Findings

Debiased lasso estimates follow a normal distribution asymptotically.

Multiple splitting improves efficiency over single splitting.

Simulation shows reduced bias and variance compared to maximum likelihood estimates.

Abstract

We consider random sample splitting for estimation and inference in high dimensional generalized linear models, where we first apply the lasso to select a submodel using one subsample and then apply the debiased lasso to fit the selected model using the remaining subsample. We show that, no matter including a prespecified subset of regression coefficients or not, the debiased lasso estimation of the selected submodel after a single splitting follows a normal distribution asymptotically. Furthermore, for a set of prespecified regression coefficients, we show that a multiple splitting procedure based on the debiased lasso can address the loss of efficiency associated with sample splitting and produce asymptotically normal estimates under mild conditions. Our simulation results indicate that using the debiased lasso instead of the standard maximum likelihood estimator in the estimation stage can vastly reduce the bias and variance of the resulting estimates. We illustrate the proposed multiple splitting debiased lasso method with an analysis of the smoking data of the Mid-South Tobacco Case-Control Study.

Tables

Table 1. Table 1: Logistic regression simulation results for n = 500 𝑛 500 n=500 , p = 700 𝑝 700 p=700 , s 0 = 6 subscript 𝑠 0 6 s_{0}=6 , and AR(1) correlation structure with parameter 0.5. Selection results refer to the lasso in the single split estimator, where the average selected model size was 41. The fitted model for estimating each β j subscript 𝛽 𝑗 \beta_{j} was S ^ ∪ { j } ^ 𝑆 𝑗 \widehat{S}\cup\{j\} . Nominal confidence interval coverage probabilities are 0.95.

	Estimator	${\beta }_{489}$	${\beta }_{130}$	${\beta }_{680}$	${\beta }_{488}$	${\beta }_{476}$	${\beta }_{190}$	${\beta }_{510}$	${\beta }_{336}$
${𝜷}_j^0$		-1.50	-1.00	-0.50	0.00	0.00	0.50	1.00	1.50
Selection Rate		1.00	1.00	0.64	0.13	0.07	0.69	1.00	1.00
Bias	Debiased SS	0.01	0.06	0.00	-0.04	-0.01	-0.03	0.00	0.00
	MLE SS	-0.52	-0.30	-0.24	-0.06	-0.01	0.18	0.37	0.54
	Debiased MS	0.01	0.04	0.02	-0.02	0.00	-0.02	-0.01	0.00
	MLE MS	-0.59	-0.37	-0.21	-0.02	-0.01	0.21	0.40	0.60
	Oracle $\left(n_2\right)$	-0.09	-0.05	0.00	-0.01	0.00	0.04	0.07	0.09
	Oracle $\left(n\right)$	-0.05	-0.01	-0.01	-0.02	0.00	0.01	0.04	0.05
Coverage	Debiased SS	0.84	0.89	0.92	0.94	0.95	0.88	0.87	0.89
	MLE SS	0.73	0.86	0.84	0.90	0.92	0.86	0.77	0.69
	Debiased MS	0.94	0.94	0.93	0.95	0.97	0.96	0.96	0.95
	MLE MS	0.39	0.64	0.81	0.94	0.94	0.79	0.58	0.40
	Oracle $\left(n_2\right)$	0.95	0.96	0.94	0.95	0.96	0.94	0.96	0.97
	Oracle $\left(n\right)$	0.96	0.95	0.94	0.92	0.94	0.95	0.96	0.96
Rejection Rate	Debiased SS	1.00	0.99	0.74	0.06	0.05	0.62	1.00	1.00
($H_0:{\beta }_j=0$)	MLE SS	1.00	0.99	0.77	0.10	0.08	0.71	1.00	1.00
	Debiased MS	1.00	1.00	0.92	0.05	0.03	0.94	1.00	1.00
	MLE MS	1.00	1.00	0.94	0.06	0.06	0.95	1.00	1.00
	Oracle $\left(n_2\right)$	1.00	1.00	0.74	0.05	0.04	0.80	1.00	1.00
	Oracle $\left(n\right)$	1.00	1.00	0.97	0.09	0.06	0.98	1.00	1.00
Standard Error	Debiased SS	0.22	0.20	0.19	0.21	0.18	0.19	0.20	0.22
	MLE SS	0.37	0.30	0.28	0.29	0.25	0.27	0.31	0.37
	Debiased MS	0.18	0.15	0.14	0.14	0.12	0.14	0.16	0.18
	MLE MS	0.27	0.22	0.20	0.22	0.19	0.20	0.22	0.27
	Oracle $\left(n_2\right)$	0.26	0.22	0.20	0.22	0.19	0.20	0.22	0.26
	Oracle $\left(n\right)$	0.18	0.15	0.14	0.15	0.13	0.14	0.15	0.18
Empirical SD	Debiased SS	0.30	0.23	0.21	0.22	0.18	0.24	0.26	0.29
	MLE SS	0.58	0.41	0.36	0.38	0.31	0.40	0.48	0.56
	Debiased MS	0.17	0.14	0.14	0.14	0.12	0.14	0.14	0.18
	MLE MS	0.29	0.22	0.22	0.23	0.20	0.21	0.24	0.30
	Oracle $\left(n_2\right)$	0.28	0.22	0.20	0.22	0.19	0.20	0.22	0.26
	Oracle $\left(n\right)$	0.17	0.14	0.14	0.16	0.14	0.14	0.14	0.17

Table 2. Table 2: Logistic regression simulation results for n = 500 𝑛 500 n=500 , p = 700 𝑝 700 p=700 , s 0 = 6 subscript 𝑠 0 6 s_{0}=6 , and compound symmetry correlation structure with parameter 0.5. Selection results refer to the lasso in the single split estimator, where the average selected model size was 37. The fitted model for estimating each β j subscript 𝛽 𝑗 \beta_{j} was S ^ ∪ { j } ^ 𝑆 𝑗 \widehat{S}\cup\{j\} . Nominal confidence interval coverage probabilities are 0.95.

	Estimator	${\beta }_{489}$	${\beta }_{130}$	${\beta }_{680}$	${\beta }_{488}$	${\beta }_{476}$	${\beta }_{190}$	${\beta }_{510}$	${\beta }_{336}$
${𝜷}_j^0$		-1.50	-1.00	-0.50	0.00	0.00	0.50	1.00	1.50
Selection Rate		1.00	0.96	0.51	0.07	0.05	0.44	0.98	1.00
Bias	Debiased SS	0.01	0.04	0.04	-0.02	0.00	0.00	-0.02	-0.05
	MLE SS	-0.42	-0.27	-0.13	-0.02	-0.01	0.19	0.30	0.37
	Debiased MS	0.02	0.03	0.03	-0.01	0.00	-0.02	-0.02	-0.02
	MLE MS	-0.39	-0.26	-0.14	-0.01	0.00	0.15	0.26	0.39
	Oracle $\left(n_2\right)$	-0.05	-0.04	-0.02	0.00	0.00	0.01	0.04	0.07
	Oracle $\left(n\right)$	-0.03	-0.02	-0.02	-0.01	0.00	0.02	0.02	0.03
Coverage	Debiased SS	0.88	0.92	0.92	0.94	0.96	0.92	0.90	0.89
	MLE SS	0.80	0.89	0.91	0.90	0.93	0.86	0.85	0.82
	Debiased MS	0.94	0.96	0.94	0.96	0.96	0.94	0.94	0.96
	MLE MS	0.68	0.84	0.90	0.96	0.94	0.92	0.81	0.76
	Oracle $\left(n_2\right)$	0.94	0.94	0.95	0.91	0.94	0.92	0.96	0.94
	Oracle $\left(n\right)$	0.96	0.98	0.96	0.92	0.92	0.96	0.96	0.98
Rejection Rate	Debiased SS	1.00	0.99	0.49	0.06	0.04	0.57	0.98	0.99
($H_0:{\beta }_j=0$)	MLE SS	1.00	1.00	0.56	0.10	0.07	0.63	0.98	0.99
	Debiased MS	1.00	1.00	0.82	0.04	0.04	0.87	1.00	1.00
	MLE MS	1.00	1.00	0.84	0.04	0.06	0.88	1.00	1.00
	Oracle $\left(n_2\right)$	1.00	1.00	0.66	0.09	0.06	0.66	1.00	1.00
	Oracle $\left(n\right)$	1.00	1.00	0.92	0.07	0.07	0.92	1.00	1.00
Standard Error	Debiased SS	0.25	0.23	0.23	0.23	0.23	0.23	0.23	0.25
	MLE SS	0.36	0.32	0.30	0.28	0.29	0.30	0.32	0.36
	Debiased MS	0.19	0.17	0.16	0.15	0.15	0.16	0.17	0.19
	MLE MS	0.27	0.23	0.21	0.21	0.21	0.21	0.23	0.27
	Oracle $\left(n_2\right)$	0.27	0.24	0.22	0.22	0.22	0.22	0.24	0.27
	Oracle $\left(n\right)$	0.18	0.17	0.15	0.15	0.15	0.15	0.17	0.18
Empirical SD	Debiased SS	0.29	0.25	0.26	0.22	0.24	0.26	0.27	0.29
	MLE SS	0.56	0.43	0.39	0.34	0.38	0.40	0.50	0.53
	Debiased MS	0.18	0.16	0.16	0.15	0.15	0.16	0.17	0.18
	MLE MS	0.28	0.24	0.22	0.22	0.22	0.21	0.24	0.27
	Oracle $\left(n_2\right)$	0.28	0.25	0.23	0.25	0.23	0.23	0.25	0.29
	Oracle $\left(n\right)$	0.18	0.15	0.16	0.16	0.16	0.15	0.16	0.17

Table 3. Table 3: Logistic regression simulation results for n = 500 𝑛 500 n=500 , p = 700 𝑝 700 p=700 , s 0 = 6 subscript 𝑠 0 6 s_{0}=6 , and AR(1) correlation structure with parameter 0.5. Selection results refer to the lasso in the single split estimator, where the average selected model size was 41. The fitted model for estimating each β j subscript 𝛽 𝑗 \beta_{j} was S ^ ^ 𝑆 \widehat{S} . Nominal confidence interval coverage probabilities are 0.95.

	Estimator	${\beta }_{489}$	${\beta }_{130}$	${\beta }_{680}$	${\beta }_{488}$	${\beta }_{476}$	${\beta }_{190}$	${\beta }_{510}$	${\beta }_{336}$
${𝜷}_j^0$		-1.50	-1.00	-0.50	0.00	0.00	0.50	1.00	1.50
Selection Rate		1.00	1.00	0.64	0.13	0.07	0.69	1.00	1.00
Bias	Debiased SS	0.00	0.06	0.18	-0.01	0.00	-0.18	0.00	0.00
	MLE SS	-0.66	-0.35	0.02	-0.02	0.01	-0.03	0.44	0.66
	Oracle $\left(n_2\right)$	-0.09	-0.05	0.00	0.00	0.00	0.04	0.07	0.09
Coverage	Debiased SS	0.83	0.90	0.57	0.99	1.00	0.59	0.87	0.89
	MLE SS	0.70	0.83	0.51	0.98	1.00	0.57	0.75	0.66
	Oracle $\left(n_2\right)$	0.95	0.96	0.94	1.00	1.00	0.94	0.96	0.97
Rejection Rate	Debiased SS	1.00	0.99	0.47	0.01	0.00	0.42	1.00	1.00
($H_0:{\beta }_j=0$)	MLE SS	1.00	0.99	0.49	0.02	0.00	0.45	1.00	1.00
	Oracle $\left(n_2\right)$	1.00	1.00	0.74	0.00	0.00	0.80	1.00	1.00
Standard Error	Debiased SS	0.22	0.20	0.12	0.03	0.01	0.13	0.20	0.22
	MLE SS	0.40	0.32	0.18	0.04	0.02	0.20	0.33	0.40
	Oracle $\left(n_2\right)$	0.26	0.22	0.20	0.00	0.00	0.20	0.22	0.26
Empirical SD	Debiased SS	0.32	0.23	0.30	0.10	0.04	0.30	0.26	0.29
	MLE SS	0.83	0.49	0.47	0.18	0.07	0.50	0.55	0.75
	Oracle $\left(n_2\right)$	0.28	0.22	0.20	0.00	0.00	0.20	0.22	0.26

Table 4. Table 4: Multiple splitting debiased lasso estimates for a logistic regression model of smoking status on standardized age, gender, and 32,557 SNPs, based on a case-control sample of 3317 individuals. There are 16 measured SNPs that lie in the regions of interest after data cleaning. P-values are adjusted for multiple testing using the Holm procedure and truncated at one, while confidence intervals are not adjusted.

Covariate	Gene	$\hat{\beta }$	SE	Holm P-value	Odds Ratio (95% CI)
Intercept		-0.01	0.18	1.00	0.99 (0.70, 1.41)
Age		0.03	0.03	1.00	1.03 (0.96, 1.09)
Male		0.45	0.07	$1.76\times {10}^{-10}$	1.58 (1.38, 1.79)
rs35327613	CHRNB3	-0.06	0.11	1.00	0.94 (0.76, 1.17)
rs61740655	CHRNA5	-0.10	0.13	1.00	0.90 (0.70, 1.17)
rs79109919	CHRNA5	-0.06	0.16	1.00	0.94 (0.68, 1.29)
rs16969968	CHRNA5	-0.09	0.11	1.00	0.91 (0.74, 1.13)
rs938682	CHRNA3	0.10	0.10	1.00	1.11 (0.91, 1.35)
rs8042374	CHRNA3	-0.18	0.12	1.00	0.83 (0.66, 1.05)
rs61737502	CHRNB4	0.04	0.10	1.00	1.04 (0.85, 1.28)
rs56218866	CHRNB4	-0.10	0.14	1.00	0.91 (0.68, 1.21)
rs950776	CHRNB4	-0.02	0.08	1.00	0.98 (0.83, 1.15)
rs12440298	CHRNB4	-0.03	0.06	1.00	0.97 (0.87, 1.09)
rs3865452	COQ8B	-0.06	0.06	1.00	0.94 (0.84, 1.04)
rs3733829	EGLN2	0.32	0.12	0.13	1.37 (1.09, 1.73)
rs75152309	CYP2A7	0.07	0.09	1.00	1.07 (0.90, 1.28)
rs73032311	CYP2A7	-0.09	0.08	1.00	0.91 (0.78, 1.07)
rs28399499	CYP2B6	0.00	0.10	1.00	1.00 (0.82, 1.22)
rs7260329	CYP2B6	-0.16	0.08	0.78	0.86 (0.73, 1.00)

Table 5. Table 5: Logistic regression simulation results for n = 2000 𝑛 2000 n=2000 , p = 2800 𝑝 2800 p=2800 , s 0 = 6 subscript 𝑠 0 6 s_{0}=6 , and AR(1) correlation structure with parameter 0.5. Selection results refer to the lasso in the single split estimator, where the average selected model size was 78. The fitted model for estimating each β j subscript 𝛽 𝑗 \beta_{j} was S ^ ^ 𝑆 \widehat{S} . Nominal confidence interval coverage probabilities are 0.95.

	Estimator	${\beta }_{489}$	${\beta }_{130}$	${\beta }_{680}$	${\beta }_{488}$	${\beta }_{137}$	${\beta }_{190}$	${\beta }_{510}$	${\beta }_{336}$
${𝜷}_j^0$		-1.50	-1.00	-0.50	0.00	0.00	0.50	1.00	1.50
Selection Rate		1.00	1.00	1.00	0.11	0.03	1.00	1.00	1.00
Bias	Debiased SS	-0.01	-0.01	0.01	-0.00	0.00	-0.01	0.02	0.01
	MLE SS	-0.24	-0.16	-0.08	-0.00	0.00	0.07	0.17	0.24
	Oracle $\left(n_2\right)$	-0.03	-0.02	-0.00	0.00	0.00	0.01	0.03	0.03
Coverage	Debiased SS	0.89	0.90	0.98	0.99	1.00	0.91	0.93	0.92
	MLE SS	0.64	0.71	0.90	0.99	1.00	0.86	0.74	0.65
	Oracle $\left(n_2\right)$	0.96	0.93	0.97	1.00	1.00	0.94	0.95	0.95
Rejection Rate	Debiased SS	1.00	1.00	1.00	0.01	0.00	1.00	1.00	1.00
($H_0:{\beta }_j=0$)	MLE SS	1.00	1.00	1.00	0.01	0.00	1.00	1.00	1.00
	Oracle $\left(n_2\right)$	1.00	1.00	1.00	0.00	0.00	1.00	1.00	1.00
Standard Error	Debiased SS	0.11	0.10	0.09	0.01	0.00	0.09	0.10	0.11
	MLE SS	0.15	0.13	0.11	0.01	0.00	0.11	0.13	0.15
	Oracle $\left(n_2\right)$	0.12	0.11	0.09	0.00	0.00	0.09	0.11	0.12
Empirical SD	Debiased SS	0.14	0.11	0.10	0.03	0.01	0.11	0.11	0.13
	MLE SS	0.21	0.17	0.12	0.05	0.02	0.14	0.16	0.22
	Oracle $\left(n_2\right)$	0.12	0.11	0.09	0.00	0.00	0.10	0.10	0.12

Table 6. Table 6: Logistic regression simulation results for n = 500 𝑛 500 n=500 , p = 700 𝑝 700 p=700 , s 0 = 6 subscript 𝑠 0 6 s_{0}=6 , and AR(1) correlation structure with parameter 0.8. Selection results refer to the lasso in the single split estimator, where the average selected model size was 38. The fitted model for estimating each β j subscript 𝛽 𝑗 \beta_{j} was S ^ ∪ { j } ^ 𝑆 𝑗 \widehat{S}\cup\{j\} . Nominal confidence interval coverage probabilities are 0.95.

	Estimator	${\beta }_{489}$	${\beta }_{130}$	${\beta }_{680}$	${\beta }_{488}$	${\beta }_{476}$	${\beta }_{190}$	${\beta }_{510}$	${\beta }_{336}$
${𝜷}_j^0$		-1.50	-1.00	-0.50	0.00	0.00	0.50	1.00	1.50
Selection Rate		1.00	0.97	0.56	0.24	0.05	0.64	0.98	1.00
Bias	Debiased SS	0.04	0.03	0.00	-0.01	-0.03	-0.02	-0.02	-0.02
	MLE SS	-0.54	-0.37	-0.24	0.02	-0.05	0.21	0.39	0.55
	Debiased MS	0.01	0.03	0.01	-0.02	-0.01	-0.01	-0.01	-0.01
	MLE MS	-0.54	-0.35	-0.23	-0.01	-0.01	0.22	0.39	0.55
	Oracle $\left(n_2\right)$	-0.08	-0.03	-0.02	-0.03	-0.01	0.02	0.05	0.08
	Oracle $\left(n\right)$	-0.04	-0.02	-0.02	-0.01	-0.01	0.01	0.03	0.04
Coverage	Debiased SS	0.89	0.93	0.92	0.96	0.96	0.90	0.91	0.90
	MLE SS	0.77	0.82	0.85	0.93	0.93	0.88	0.81	0.74
	Debiased MS	0.93	0.92	0.93	0.96	0.95	0.94	0.90	0.92
	MLE MS	0.58	0.73	0.86	0.97	0.94	0.85	0.70	0.53
	Oracle $\left(n_2\right)$	0.93	0.96	0.95	0.93	0.94	0.96	0.94	0.95
	Oracle $\left(n\right)$	0.94	0.96	0.92	0.96	0.96	0.94	0.93	0.96
Rejection Rate	Debiased SS	0.99	0.92	0.59	0.04	0.04	0.57	0.95	1.00
($H_0:{\beta }_j=0$)	MLE SS	0.99	0.90	0.64	0.07	0.07	0.60	0.94	0.99
	Debiased MS	1.00	0.99	0.90	0.04	0.05	0.88	1.00	1.00
	MLE MS	1.00	0.99	0.90	0.03	0.06	0.89	1.00	1.00
	Oracle $\left(n_2\right)$	1.00	1.00	0.74	0.07	0.06	0.78	1.00	1.00
	Oracle $\left(n\right)$	1.00	1.00	0.98	0.04	0.04	0.98	1.00	1.00
Standard Error	Debiased SS	0.27	0.26	0.24	0.32	0.20	0.23	0.25	0.27
	MLE SS	0.43	0.39	0.33	0.44	0.27	0.33	0.38	0.44
	Debiased MS	0.19	0.17	0.15	0.20	0.12	0.15	0.17	0.19
	MLE MS	0.30	0.26	0.22	0.32	0.19	0.22	0.26	0.30
	Oracle $\left(n_2\right)$	0.26	0.22	0.20	0.32	0.19	0.20	0.22	0.25
	Oracle $\left(n\right)$	0.17	0.15	0.14	0.22	0.13	0.14	0.15	0.17
Empirical SD	Debiased SS	0.33	0.30	0.26	0.30	0.18	0.26	0.29	0.33
	MLE SS	0.72	0.59	0.43	0.55	0.31	0.44	0.52	0.69
	Debiased MS	0.21	0.20	0.17	0.21	0.12	0.17	0.19	0.20
	MLE MS	0.35	0.32	0.26	0.33	0.19	0.26	0.30	0.34
	Oracle $\left(n_2\right)$	0.29	0.22	0.20	0.34	0.19	0.20	0.23	0.26
	Oracle $\left(n\right)$	0.19	0.15	0.14	0.23	0.13	0.15	0.17	0.18

Table 7. Table 7: Logistic regression simulation results for n = 500 𝑛 500 n=500 , p = 700 𝑝 700 p=700 , s 0 = 6 subscript 𝑠 0 6 s_{0}=6 , and compound symmetry correlation structure with parameter 0.8. Selection results refer to the lasso in the single split estimator, where the average selected model size was 23. The fitted model for estimating each β j subscript 𝛽 𝑗 \beta_{j} was S ^ ∪ { j } ^ 𝑆 𝑗 \widehat{S}\cup\{j\} . Nominal confidence interval coverage probabilities are 0.95.

	Estimator	${\beta }_{489}$	${\beta }_{130}$	${\beta }_{680}$	${\beta }_{488}$	${\beta }_{476}$	${\beta }_{190}$	${\beta }_{510}$	${\beta }_{336}$
${𝜷}_j^0$		-1.50	-1.00	-0.50	0.00	0.00	0.50	1.00	1.50
Selection Rate		0.98	0.62	0.24	0.01	0.03	0.22	0.68	0.96
Bias	Debiased SS	0.04	0.04	0.03	-0.02	0.00	-0.05	-0.08	-0.06
	MLE SS	-0.14	-0.09	-0.04	-0.02	0.01	0.01	0.04	0.11
	Debiased MS	0.03	0.08	0.02	-0.02	0.00	-0.02	-0.07	-0.07
	MLE MS	-0.17	-0.07	-0.07	-0.03	-0.01	0.06	0.08	0.13
	Oracle $\left(n_2\right)$	-0.07	0.03	-0.07	-0.03	0.00	0.05	0.03	0.02
	Oracle $\left(n\right)$	-0.03	0.02	-0.03	-0.01	0.00	0.03	0.01	0.01
Coverage	Debiased SS	0.90	0.91	0.91	0.91	0.90	0.94	0.91	0.91
	MLE SS	0.89	0.89	0.91	0.91	0.89	0.93	0.88	0.90
	Debiased MS	0.90	0.94	0.91	0.96	0.96	0.93	0.92	0.90
	MLE MS	0.86	0.90	0.93	0.96	0.95	0.93	0.88	0.88
	Oracle $\left(n_2\right)$	0.94	0.94	0.94	0.96	0.96	0.94	0.96	0.96
	Oracle $\left(n\right)$	0.94	0.94	0.92	0.96	0.93	0.94	0.92	0.94
Rejection Rate	Debiased SS	1.00	0.82	0.33	0.09	0.10	0.33	0.77	0.99
($H_0:{\beta }_j=0$)	MLE SS	1.00	0.85	0.36	0.09	0.11	0.34	0.78	0.98
	Debiased MS	1.00	0.98	0.59	0.04	0.04	0.56	0.98	1.00
	MLE MS	1.00	0.98	0.58	0.04	0.05	0.58	0.98	1.00
	Oracle $\left(n_2\right)$	1.00	0.88	0.47	0.04	0.04	0.44	0.90	1.00
	Oracle $\left(n\right)$	1.00	1.00	0.67	0.04	0.07	0.71	1.00	1.00
Standard Error	Debiased SS	0.32	0.32	0.31	0.31	0.31	0.31	0.32	0.32
	MLE SS	0.38	0.37	0.35	0.34	0.35	0.35	0.36	0.38
	Debiased MS	0.23	0.23	0.22	0.21	0.21	0.22	0.23	0.23
	MLE MS	0.27	0.26	0.25	0.25	0.25	0.26	0.26	0.26
	Oracle $\left(n_2\right)$	0.34	0.31	0.30	0.31	0.31	0.31	0.32	0.33
	Oracle $\left(n\right)$	0.23	0.22	0.21	0.21	0.21	0.21	0.22	0.23
Empirical SD	Debiased SS	0.36	0.37	0.34	0.33	0.35	0.35	0.39	0.38
	MLE SS	0.51	0.46	0.40	0.40	0.41	0.42	0.50	0.49
	Debiased MS	0.26	0.25	0.24	0.20	0.22	0.23	0.25	0.26
	MLE MS	0.35	0.32	0.29	0.24	0.26	0.27	0.31	0.34
	Oracle $\left(n_2\right)$	0.32	0.33	0.33	0.29	0.30	0.34	0.32	0.34
	Oracle $\left(n\right)$	0.24	0.23	0.24	0.20	0.23	0.22	0.24	0.24

Equations

$${\rho }_{\mathbf{\beta }}(y,\mathbf{x})=\rho (y,{\mathbf{x}}^T\mathbf{\beta })=A({\mathbf{x}}^T\mathbf{\beta })-y{\mathbf{x}}^T\mathbf{\beta }$$

$${\tilde{\mathbf{\beta }}}^{\lambda }=\underset{\mathbf{\beta }}{\mathrm{argmin}}{P}_n{\rho }_{\mathbf{\beta }}+\lambda {\Vert \mathbf{\beta }\Vert }_1.$$

$${\tilde{\mathbf{\beta }}}^{DL}={\tilde{\mathbf{\beta }}}^{\lambda }-{\hat{\mathbf{\Theta }}}_{{\tilde{\mathbf{\beta }}}^{\lambda }}{P}_n{\dot{\rho }}_{{\tilde{\mathbf{\beta }}}^{\lambda }}$$

$$\widehat{Var}\left({\hat{\mathbf{\beta }}}_j\right)={\hat{V}}_j-\frac{n(n_2)}{B^2(n-n_2)}\sum _{b=1}^B{\left({\tilde{\mathbf{\beta }}}_{j,b}^{DL}-{\hat{\mathbf{\beta }}}_j\right)}^2,$$

$${\hat{V}}_j=\frac{n(n-1)}{(n-n_2{)}^2}\sum _{i=1}^n{\left\{\frac{1}{B}\sum _{b=1}^B\left({\mathbf{v}}_{b,i}-{\bar{\mathbf{v}}}_i\right)({\tilde{\mathbf{\beta }}}_{j,b}^{DL}-{\hat{\mathbf{\beta }}}_j)\right\}}^2,$$

$$\widehat{Var}\left({\hat{\mathbf{\beta }}}_S\right)=\frac{n(n-1)}{(n-n_2{)}^2}\sum _{i=1}^n{\hat{\mathbf{C}}}_i{\hat{\mathbf{C}}}_i^T-\frac{n(n_2)}{B^2(n-n_2)}\sum _{b=1}^B{\hat{\mathbf{D}}}_b{\hat{\mathbf{D}}}_b^T$$

$${\hat{\mathbf{C}}}_i=\frac{1}{B}\sum _{b=1}^B\left({\mathbf{v}}_{b,i}-{\bar{\mathbf{v}}}_i\right)\left({\tilde{\mathbf{\beta }}}_{S,b}^{DL}-{\hat{\mathbf{\beta }}}_S\right),{\hat{\mathbf{D}}}_b={\tilde{\mathbf{\beta }}}_{S,b}^{DL}-{\hat{\mathbf{\beta }}}_S.$$

$$\underset{a_0\in \{{\mathbf{x}}_i^T{\mathbf{\beta }}^0\}}{\max }\underset{\max \left(|a-a_0|,|\hat{a}-a_0|\right)\le \delta }{\sup }\underset{y\in \mathcal{Y}}{\sup }\frac{|\ddot{\rho }(y,a)-\ddot{\rho }(y,\hat{a})|}{|a-\hat{a}|}\le c_{Lip}.$$

$$\underset{a_0\in \{{\mathbf{x}}_i^T{\mathbf{\beta }}^0\}}{\max }\underset{y\in \mathcal{Y}}{\sup }|\dot{\rho }(y,a_0)|\le K_1,$$

$$\underset{a_0\in \{{\mathbf{x}}_i^T{\mathbf{\beta }}^0\}}{\max }\underset{y\in \mathcal{Y}}{\sup }|\ddot{\rho }(y,a_0)|\le K_2.$$

$$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})$$

$$\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).$$

$$\sqrt{n}{\mathbf{a}}_S^T({\hat{\mathbf{\beta }}}_S-{\mathbf{\beta }}_S^0)=\frac{(1-q{)}^{-1}}{\sqrt{n}}\sum _{i=1}^{n}E\left[{\mathbf{a}}_{\hat{S}}^T{\mathbf{\Sigma }}_{{\mathbf{\beta }}^0,\hat{S}}^{-1}{\tilde{\mathbf{I}}}_{\hat{S}}{\mathbf{v}}_i|\mathcal{Z}\right]{\dot{\rho }}_{{\mathbf{\beta }}^0}(y_i,{\mathbf{x}}_i)+o_p\left(1\right).$$

$$P\left(\hat{S}\supseteq S_0,|\hat{S}|\le k{s}_0\right)\ge 1-c_1/p^{c_2}$$

$$\frac{\sqrt{n_2}{\mathbf{a}}_{\hat{S}}^T({\tilde{\mathbf{\beta }}}_{\hat{S}}^{DL}-{\mathbf{\beta }}_{\hat{S}}^0)}{({\mathbf{a}}_{\hat{S}}^T{\hat{\mathbf{\Theta }}}_{{\tilde{\mathbf{\beta }}}^{\lambda },\hat{S}}{\mathbf{a}}_{\hat{S}}{)}^{1/2}}\to N(0,1)$$

$$\frac{\sqrt{n}{\mathbf{a}}_S^T({\hat{\mathbf{\beta }}}_S-{\mathbf{\beta }}_S^0)}{({\mathbf{\eta }}_n^T{\mathbf{\Sigma }}_{{\mathbf{\beta }}^0}{\mathbf{\eta }}_n{)}^{1/2}}\to N(0,1)$$

$$\underset{i\le n}{\max }\left|\log A^{\prime \prime }({\mathbf{x}}_i^T{\mathbf{\beta }}^0+t)-\log A^{\prime \prime }({\mathbf{x}}_i^T{\mathbf{\beta }}^0)\right|$$

$$=\underset{i\le n}{\max }\left|\log \frac{A^{\prime \prime }({\mathbf{x}}_i^T{\mathbf{\beta }}^0+t)}{A^{\prime \prime }({\mathbf{x}}_i^T{\mathbf{\beta }}^0)}\right|$$

$$\le \underset{i\le n}{\max }\frac{|A^{\prime \prime }({\mathbf{x}}_i^T{\mathbf{\beta }}^0+t)-A^{\prime \prime }({\mathbf{x}}_i^T{\mathbf{\beta }}^0)|}{\min \left(A^{\prime \prime }({\mathbf{x}}_i^T{\mathbf{\beta }}^0),A^{\prime \prime }({\mathbf{x}}_i^T{\mathbf{\beta }}^0+t)\right)}$$

$$\le \frac{c_{Lip}|t|}{{\inf }_{i\le n,|t|\le \delta }{A}^{\prime \prime }({\mathbf{x}}_i^T{\mathbf{\beta }}^0+t)}$$

$${\kappa }_{+}(s)=\underset{|B|=s,B\cap S=\varnothing ,\mathbf{b}=\in \mathcal{C}\left(\xi ,S\right),M_3{\Vert {\mathbf{b}}_S\Vert }_1\le {\eta }^{\ast }}{\sup }{\lambda }_{max}\left({\int }_0^1\frac{1}{n}\sum _{i=1}^n{A}^{\prime \prime }({\mathbf{x}}_i^T{\mathbf{\beta }}^0+t{\mathbf{x}}_i^T\mathbf{b}){\left[{\mathbf{x}}_i{\mathbf{x}}_i^T\right]}_{B}dt\right)$$

$$P\left(\underset{2\le j\le p+1}{\max }\left|\frac{1}{n}\sum _{i=1}^n{x}_{ij}^2-{\mathbf{\Sigma }}_{jj}\right|\ge t\right)$$

$$P(\mathcal{A})$$

$$\sqrt{n_2}{\mathbf{a}}_{\hat{S}}^T\left({\tilde{\mathbf{\beta }}}_{\hat{S}}^{DL}-{\mathbf{\beta }}_{\hat{S}}^0\right)=-\frac{1}{\sqrt{n_2}}{\mathbf{a}}_{\hat{S}}^T{\mathbf{\Sigma }}_{{\mathbf{\beta }}^0,\hat{S}}^{-1}\sum _{i\in {\mathcal{D}}_2}{\dot{\rho }}_{{\mathbf{\beta }}_{\hat{S}}^0}(y_i,({\mathbf{x}}_i{)}_{\hat{S}})+r_n$$

$$P\left(|r_n|\ge \epsilon \right)$$

$$\sqrt{n_2}{\mathbf{a}}_S^T({\tilde{\mathbf{\beta }}}_{S,b}^{DL}-{\mathbf{\beta }}_S^0)$$

$$\sqrt{n}{\mathbf{a}}_S^T({\hat{\mathbf{\beta }}}_S-{\mathbf{\beta }}_S^0)$$

$$(1-q{)}^{-1}E\left[{\mathbf{a}}_{{\hat{S}}_b}^T{\mathbf{\Sigma }}_{{\mathbf{\beta }}^0,{\hat{S}}_b}^{-1}{\tilde{\mathbf{I}}}_{{\hat{S}}_b}{v}_{b,i}|\mathcal{Z}\right]-{\mathbf{\eta }}_n^T$$

$$=(1-q{)}^{-1}E\left[{\mathbf{a}}_{{\tilde{S}}_b}^T{\mathbf{\Sigma }}_{{\mathbf{\beta }}^0,{\tilde{S}}_b}^{-1}{\tilde{\mathbf{I}}}_{{\tilde{S}}_b}{v}_{b,i}|\mathcal{Z}\right]-{\mathbf{\eta }}_n^T$$

$$+(1-q{)}^{-1}\left(E\left[{\mathbf{a}}_{{\hat{S}}_b}^T{\mathbf{\Sigma }}_{{\mathbf{\beta }}^0,{\hat{S}}_b}^{-1}{\tilde{\mathbf{I}}}_{{\hat{S}}_b}{v}_{b,i}|\mathcal{Z}\right]-E\left[{\mathbf{a}}_{{\tilde{S}}_b}^T{\mathbf{\Sigma }}_{{\mathbf{\beta }}^0,{\tilde{S}}_b}^{-1}{\tilde{\mathbf{I}}}_{{\tilde{S}}_b}{v}_{b,i}|\mathcal{Z}\right]\right),$$

$$r_n=\underset{r_{n,1}}{\underbrace{\left(E\left[{\mathbf{a}}_{\tilde{S}}^T{\mathbf{\Sigma }}_{{\mathbf{\beta }}^0,\tilde{S}}^{-1}{\tilde{\mathbf{I}}}_{\tilde{S}}|\mathcal{Z}\right]-{\mathbf{\eta }}_n^T\right)\sqrt{n}{P}_n{\dot{\rho }}_{{\mathbf{\beta }}^0}}}+\underset{r_{n,2}}{\underbrace{\frac{(1-q{)}^{-1}}{\sqrt{n}}\sum _{i=1}^n{h}_{i,n}}}$$

$$P\left(\underset{j}{\max }\left|\sqrt{n}{\left(P_n{\dot{\rho }}_{{\mathbf{\beta }}^0}\right)}_j\right|\ge t\sqrt{\log (p)}\right)$$

$$|r_{n,1}|$$

Full text

Keywords: Asymptotic normality, genetic marker, high-dimensional inference, single nucleotide polymorphism (SNP), sparse regression.

1 Introduction

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.

2 Methods

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

2.1 Notation

For a positive integer $p$, we denote the size of any index set $S\subseteq\{1,\ldots,p\}$ by $|S|$. For a $p$-dimensional vector $\boldsymbol{b}$ and a $p\times p$ dimensional matrix $\boldsymbol{B}$, let $\boldsymbol{b}_{S}$ denote the $|S|$-dimensional subvector with entries indexed by $S$, and similarly $\boldsymbol{B}_{S}$ denote the $|S|\times|S|$ submatrix with rows and columns indexed by $S$. The $\ell_{q}$ norm is denoted as $\left\|\boldsymbol{b}\right\|_{q}$ for $q\geq 1$. For positive sequences $a_{n}$ and $b_{n}$, we write $a_{n}=\mathcal{O}\left(b_{n}\right)$ if there exists a constant $c>0$ and $N>0$ such that $a_{n}/b_{n}<c$ for all $n>N$, we write $a_{n}=o(b_{n})$ if $a_{n}/b_{n}\rightarrow 0$ as $n\rightarrow\infty$, and we write $a_{n}\sim b_{n}$ if $a_{n}=\mathcal{O}\left(b_{n}\right)$ and $b_{n}=\mathcal{O}\left(a_{n}\right)$.

Let $(y_{i},\boldsymbol{\widetilde{x}}_{i})$, $i=1,\ldots,n$, be independent and identically distributed copies of $(y,\boldsymbol{\widetilde{x}})$, where $y$ is a scalar-valued response variable and $\boldsymbol{\widetilde{x}}$ is a $p$-dimensional vector of covariates. We consider the high dimensional setting where $p$ can be much larger than $n$. Let $\boldsymbol{x}_{i}=(1,\boldsymbol{\widetilde{x}}_{i}^{T})^{T}$ and denote the $n\times(p+1)$ design matrix by $\boldsymbol{X}$, with $i$th row $\boldsymbol{x}_{i}^{T}$ and $j$th column $\boldsymbol{X}_{j}$. We assume without loss of generality that $\boldsymbol{\widetilde{x}}$ has mean zero. For any function $f(y,\boldsymbol{x})$, we define $P_{n}f=n^{-1}\sum_{i=1}^{n}f(y_{i},\boldsymbol{x}_{i})$.

We consider generalized linear models (McCullagh and Nelder,, 2019) with canonical link function and known dispersion parameter. Denote

[TABLE]

for a known twice-differentiable function $A(\cdot)$, and its first and second derivatives with respect to $\boldsymbol{\beta}$ as $\dot{\rho}_{\boldsymbol{\beta}}$ and $\ddot{\rho}_{\boldsymbol{\beta}}$, respectively. The negative log-likelihood for $\boldsymbol{\beta}$ is then $P_{n}\rho_{\boldsymbol{\beta}}=n^{-1}\sum_{i=1}^{n}\rho_{\boldsymbol{\beta}}(y_{i},\boldsymbol{x}_{i})$ with score function $P_{n}\dot{\rho}_{\boldsymbol{\beta}}=n^{-1}\sum_{i=1}^{n}\{A^{\prime}(\boldsymbol{x}_{i}^{T}\boldsymbol{\beta})-y_{i}\}\boldsymbol{x}_{i}$ and Hessian matrix $P_{n}\ddot{\rho}_{\boldsymbol{\beta}}(y_{i},\boldsymbol{x}_{i})=n^{-1}\sum_{i=1}^{n}A^{\prime\prime}(\boldsymbol{x}_{i}^{T}\boldsymbol{\beta})\boldsymbol{x}_{i}\boldsymbol{x}_{i}^{T}$. Denote $\boldsymbol{\Sigma}_{\boldsymbol{\beta}}=E\{P_{n}\ddot{\rho}_{\boldsymbol{\beta}}(y_{i},\boldsymbol{x}_{i})\}$. The $(p+1)$-dimensional unknown true coefficient vector $\boldsymbol{\beta}^{0}$ is assumed to be sparse, and we denote the index set of signal variables by $S_{0}=\{j:\boldsymbol{\beta}_{j}^{0}\neq 0\}$, with $s_{0}=|S_{0}|$. The quantities $p$, $\boldsymbol{\beta}^{0}$, $s_{0}$, and $S_{0}$ are allowed to change with $n$. For practical applications, one would generally consider $s_{0}$ and $\boldsymbol{\beta}^{0}_{S_{0}}$ to be fixed so that the regression coefficients have the usual interpretation in terms of the conditional mean $E[y|\boldsymbol{x}]$ as $p$ grows with $n$. Letting $s_{0}$ grow with $n$ does not maintain the log odds ratio interpretation of $\boldsymbol{\beta}^{0}$, for example, in a logistic regression model.

2.2 The lasso estimator for generalized linear models

The lasso estimator $\boldsymbol{\widetilde{\beta}}^{\lambda}$ for the GLM parameters $\boldsymbol{\beta}^{0}$ is, given a tuning parameter $\lambda>0$, a minimizer of the penalized negative log likelihood

[TABLE]

In general there is no closed form solution for $\boldsymbol{\widetilde{\beta}}^{\lambda}$, but the objective function is convex and can be optimized efficiently with widely available software (Friedman et al.,, 2010). The penalization shrinks the estimated coefficients towards zero, with many of their values set exactly to zero, resulting in a sparse estimator $\boldsymbol{\widetilde{\beta}}^{\lambda}$. Model selection may then be performed by keeping only the covariates with nonzero estimated coefficients. In practice, a small subset with a finite number of coefficients may be left unpenalized so that they are never excluded from the selected model. In our simulations and data analysis, we do not penalize the intercept.

2.3 The debiased lasso estimator

Zhang and Zhang, (2014) and van de Geer et al., (2014) proposed debiased, also called desparsified, lasso estimators for linear models and GLMs, respectively. This procedure may be used to obtain coefficient estimates and confidence intervals across the entire parameter vector $\boldsymbol{\beta}^{0}$ in high dimensional settings. It does not rely on sample splitting or variable selection. Instead, using the lasso estimator $\boldsymbol{\widetilde{\beta}}^{\lambda}$ as an initial value, it performs a single Newton iteration for minimizing the negative log likelihood $P_{n}\rho_{\boldsymbol{\beta}}$, using an approximate inverse for the Hessian matrix $P_{n}\ddot{\rho}_{\boldsymbol{\widetilde{\beta}}^{\lambda}}$ such as a nodewise lasso estimator. The resulting desparsified lasso estimator is no longer sparse and is, under further conditions, asymptotically unbiased with a normal distribution.

A key assumption for accurately estimating the inverse of the Hessian matrix in high dimensional settings, however, is the sparsity assumption on $\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0}}^{-1}=[E\{A^{\prime\prime}(\boldsymbol{x}^{T}\boldsymbol{\beta}^{0})\boldsymbol{x}\boldsymbol{x}^{T}\}]^{-1}$. For generalized linear models, each entry of this matrix generally depends on all the signal variables due to the non-identity link function, so the sparsity assumption is difficult to interpret and often may not hold in practice.

Such issues were discussed further by Xia et al., (2021), who proposed a refined debiased lasso method for lower dimensional GLMs with a diverging number of covariates. In that setting and under some standard conditions, $\boldsymbol{\widehat{\Sigma}}_{\boldsymbol{\widetilde{\beta}}^{\lambda}}=P_{n}\ddot{\rho}_{\boldsymbol{\widetilde{\beta}}^{\lambda}}$ is directly invertible with probability going to one, and the resulting estimator

[TABLE]

with estimated covariance matrix $\boldsymbol{\widehat{\Theta}}_{\boldsymbol{\widetilde{\beta}}^{\lambda}}/n=\boldsymbol{\widehat{\Sigma}}_{\boldsymbol{\widetilde{\beta}}^{\lambda}}^{-1}/n$ was shown to be asymptotically normal. Additionally, it outperforms both the original desparsified lasso estimator and the MLE in several simulation settings of lower-dimensional GLMs with sparse regression coefficients. It can also be computed using standard GLM software packages for the lasso and MLE.

2.4 The debiased lasso after single sample splitting

We propose a debiased lasso after sample splitting procedure that uses the lasso for model selection and the refined debiased lasso for estimation and inference. For a fixed splitting proportion $q\in(0,1)$, we randomly split the sample $(y_{i},\boldsymbol{x}_{i})_{i=1}^{n}$ into two subsets: $\mathcal{D}_{1}$ containing $n_{1}=qn$ individuals to be used for variable selection, and $\mathcal{D}_{2}$ containing $n_{2}=n-n_{1}$ individuals to be used for estimation. For a pre-specified fixed set of coefficient indices $S$ that is allowed to be an empty set, do the following:

using the subsample $\mathcal{D}_{1}$, fit a GLM lasso estimator $\boldsymbol{\widetilde{\beta}}^{\lambda}$ to obtain the selected model $\widehat{S}=S\cup\{j:\boldsymbol{\widetilde{\beta}}^{\lambda}_{j}\neq 0\}$; 2. 2.

using the subsample $\mathcal{D}_{2}$, compute the debiased lasso estimator $\boldsymbol{\widetilde{\beta}}^{DL}_{\widehat{S}}$ for $\boldsymbol{\beta}^{0}_{\widehat{S}}$ based on a GLM with covariates $\widehat{S}$, and use the estimated covariance matrix of $\boldsymbol{\widetilde{\beta}}^{DL}_{\widehat{S}}$, $\boldsymbol{\widehat{\Theta}}_{\boldsymbol{\widetilde{\beta}}^{\lambda},\widehat{S}}/n_{2}$, to construct $(1-\alpha)$% confidence intervals for any contrast $\boldsymbol{a}_{\widehat{S}}^{T}\boldsymbol{\beta}^{0}_{\widehat{S}}$ based on a normal distribution.

Because the two subsamples are independent, we can make asymptotically valid statistical inference on any covariate in $\widehat{S}$ as long as this set contains the true signal variables when $n$ is large enough. In practice, covariates with small coefficients may not be selected in the first stage, so making conclusions about the statistical significance of omitted variables from step 1 is not advised. The large bias that can potentially occur when $S$ is a null set is illustrated through simulations in Section 4. Furthermore, although the size of the fitted model $\widehat{S}$ must be low-dimensional, this does not prevent one from obtaining estimates for each individual coefficient. For example, step 2 may be repeated for each individual coefficient, i.e. for $S=\{1\},\{2\},\ldots,\{p\}$, along with an appropriate correction for multiple testing when computing p-values and confidence intervals.

2.5 The debiased lasso after multiple sample splitting

As previously noted, the single sample splitting estimator suffers from a loss of efficiency due to only using a fraction of the total sample to estimate $\boldsymbol{\beta}^{0}$. The following multiple splitting procedure addresses this issue. For a pre-specified fixed set of coefficient indices $S$, we now generate $B$ random splits of the sample $(y_{i},\boldsymbol{x}_{i})_{i=1}^{n}$, denoted $(\mathcal{D}_{1,b},\mathcal{D}_{2,b})$ for $b=1,\ldots,B$, and do the following:

for $b=1,\ldots,B$,

1. (a)

using the subsample $\mathcal{D}_{1,b}$, fit a GLM lasso estimator $\boldsymbol{\widetilde{\beta}}^{\lambda}_{b}$ to obtain the selected model $\widehat{S}_{b}=S\cup\{j:\boldsymbol{\widetilde{\beta}}^{\lambda}_{b,j}\neq 0\}$; 2. (b)

using the subsample $\mathcal{D}_{2,b}$, compute the debiased lasso estimator $\boldsymbol{\widetilde{\beta}}^{DL}_{S,b}$ for $\boldsymbol{\beta}^{0}_{S}$ based on a GLM with covariates $\widehat{S}_{b}$; 2. 2.

finally the multiple splitting estimators are $\boldsymbol{\widehat{\beta}}_{S}=B^{-1}\sum_{b=1}^{B}\boldsymbol{\widetilde{\beta}}^{DL}_{S,b}$.

Note that all target coefficients to estimate must be pre-specified, in contrast to the single split procedure. For exploratory analysis, where step 1b may be repeated for each individual coefficient, the multiple splitting estimator typically has a substantially more strict significance threshold compared to single split inference on $\widehat{S}$ due to correcting for $p$, $p>|\widehat{S}|$, multiple comparisons.

We apply the same variance estimator as Fei and Li, (2021), based on the nonparametric delta method with a bias correction for controlling the Monte Carlo error (Efron,, 2014; Wager et al.,, 2014). Letting $\boldsymbol{v}_{b}$ denote the $n$-vector of sampling indicators such that $\boldsymbol{v}_{b,i}\in\{0,1\}$ is equal to one if $i\in\mathcal{D}_{2,b}$, and is zero otherwise,

[TABLE]

[TABLE]

where $\bar{\boldsymbol{v}}_{i}=B^{-1}\sum_{b=1}^{B}\boldsymbol{v}_{b,i}$, and $(1-\alpha)$% confidence intervals may be constructed based on a normal approximation. Similarly, the estimated covariance matrix is

[TABLE]

[TABLE]

Both single and multiple sample splitting procedures require that the matrices $\boldsymbol{\widehat{\Sigma}}_{\boldsymbol{\widetilde{\beta}}^{\lambda},\widehat{S}}$ computed from each estimation split are invertible. One way to encourage this is to set a hard limit for the selected model size, which would set a lower bound for the parameter $\lambda$ in the lasso. In our simulations and real data analysis we do not restrict $\lambda$, which is chosen by cross validation, and instead allow the glm function in R to follow its default behavior of automatically dropping covariates from the model when the design matrix is ill-conditioned.

As the number of random splits $B$ grows we typically have diminishing gains in efficiency relative to the increased computational cost. Although our theoretical results consider the maximum possible value of $B={n\choose n_{2}}$, our simulation results indicate that $B=1,000$ is sufficient for many practical settings with up to several hundreds of individuals in the sample.

3 Theoretical Results

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

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

[TABLE]

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

[TABLE]

[TABLE] 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.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

Theorem 1**.**

(Variable screening properties for the GLM lasso.) For a choice of tuning parameter $\lambda\sim\sqrt{\log(p)/n}$ and under assumptions 1-5, the lasso estimator $\boldsymbol{\widetilde{\beta}}^{\lambda}$ and its selected model $\widehat{S}=\{j:\widetilde{\beta}_{j}^{\lambda}\neq 0\}$ satisfy

[TABLE]

for some positive constants $k,c_{1},c_{2}$.

Theorem 2**.**

(Asymptotic normality for single sample splitting.) For the single split debiased lasso estimator $\boldsymbol{\widetilde{\beta}}^{DL}_{\widehat{S}}$ with both lasso tuning parameters of the order $\sqrt{\log(p)/n}$, $s_{0}\log(s_{0})\sqrt{s_{0}/n}=o(1)$, $|S|=\mathcal{O}\left(s_{0}\right)$, and under assumptions 1-5, $\boldsymbol{\widehat{\Sigma}}_{\boldsymbol{\widetilde{\beta}}^{\lambda},\widehat{S}}$ is invertible with probability going to one, and for any $(p+1)$-vector $\boldsymbol{a}$ such that $\left\|\boldsymbol{a}\right\|_{2}=\left\|\boldsymbol{a}_{\widehat{S}}\right\|_{2}=1$, we have

[TABLE]

in distribution as $n\rightarrow\infty$.

Theorem 3**.**

(Asymptotic normality for multiple sample splitting.) For the multiple splitting debiased lasso estimator $\boldsymbol{\widehat{\beta}}_{S}$ with all lasso tuning parameters of the order $\sqrt{\log(p)/n}$, $s_{0}\log(s_{0})\sqrt{s_{0}/n}=o(1)$, $|S|=\mathcal{O}\left(s_{0}\right)$, and under assumptions 1-7, for any $(p+1)$-vector $\boldsymbol{a}$ such that $\left\|\boldsymbol{a}\right\|_{2}=\left\|\boldsymbol{a}_{S}\right\|_{2}=1$, we have

[TABLE]

in distribution as $n\rightarrow\infty$.

4 Simulations

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.

5 Real Data Example: The Mid-South Tobacco Case-Control Study

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).

Acknowledgements

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

Appendix Appendix A Proof of Theorem 1: Results for GLM lasso selected model

In order to show the model selection properties of the GLM lasso under Assumptions 1-5, we consider the results on false negatives and false positives (i.e. selected model size) separately in Lemmas 1 and 2 below. These are each based on oracle inequalities from Huang and Zhang, (2012), specifically using the unweighted lasso, target vector $\boldsymbol{\beta}^{0}$, and target set $S=S_{0}\cup\{1\}$ that includes the intercept. We use $c_{1}$ and $c_{2}$ as generic constants in our probability bounds and denote the sample size as $n$.

Lemma 1**.**

Under Assumptions 1-5, if the diagonal entries of $\boldsymbol{X}^{T}\boldsymbol{X}$ are bounded above and away from zero and the compatibility constant (Huang and Zhang,, 2012) is bounded below, then the GLM lasso has no false negatives with probability at least $1-c_{1}\exp(-c_{2}\log(p))$ for constants $c_{1},c_{2}>0$.

This lemma is a direct result of Theorem 9 (iii) of Huang and Zhang, (2012), which gives an $\ell_{1}$ bound for the GLM lasso estimation error. We proceed by verifying that our assumptions imply the required conditions in Theorem 9 of Huang and Zhang, (2012) are satisfied.

We first use the bound $\log x\in[1-1/x,x-1]$, to verify the required Lipschitz condition

[TABLE]

for all $|t|\leq\delta$, where the denominator is bounded below since $\left\|\boldsymbol{X}\boldsymbol{\beta}^{0}\right\|_{\infty}$ is bounded above almost surely and $A^{\prime\prime}$ is a positive function.

Next, we can choose penalty levels $\lambda>\lambda_{1}=\lambda_{0}$ each of order $\sqrt{\log(p)/n}$ so that (31) of Huang and Zhang, (2012) is satisfied for large $n$ when the compatibility constant $\kappa_{*}(\xi,S)$ is bounded below and $|S|\sqrt{\log(p)/n}=o(1)$, as we have assumed.

Lastly, to verify condition (29) of Huang and Zhang, (2012), note that the positive quantities $A^{\prime\prime}(\boldsymbol{x}_{i}^{T}\boldsymbol{\beta}^{0}),i=1,\ldots,n$ are bounded almost surely, so the diagonal elements of $P_{n}\ddot{\rho}_{\boldsymbol{\beta}^{0}}=\frac{1}{n}\sum_{i=1}^{n}A^{\prime\prime}(\boldsymbol{x}_{i}^{T}\boldsymbol{\beta}^{0})\boldsymbol{x}\boldsymbol{x}_{i}^{T}$ are bounded above and away from zero due to the assumed corresponding bounds on $\frac{1}{n}\sum_{i=1}^{n}x_{ij}^{2}$ for all $j=1,\ldots,p+1$.

Now that we have verified the conditions, by Theorem 9 of Huang and Zhang, (2012) using the seminorm $\phi(\boldsymbol{b})=\left\|\boldsymbol{b}_{S}\right\|_{1}/|S|$, with probability at least $1-c_{1}\exp(-c_{2}\log(p))$ for constants $c_{1},c_{2}>0$, we have $\left\|\boldsymbol{\widetilde{\beta}}^{\lambda}_{S}-\boldsymbol{\beta}_{S}^{0}\right\|_{1}\leq k|S|\sqrt{\log(p)/n}$ for some constant $k>0$. Therefore, on this event, there are no false negatives if $\min_{j\in S}|\boldsymbol{\beta}_{j}^{0}|>k(s_{0}+1)\sqrt{\log(p)/n}$, as we have assumed.

Lemma 2**.**

Under the conditions of Lemma 1, the GLM lasso selected model size has an upper bound proportional to the true model size $s_{0}$, with probability at least $1-c_{1}\exp(c_{2}n)$ for some constants $c_{1},c_{2}>0$.

This lemma is a direct result of Theorem 19 of Huang and Zhang, (2012), which gives an upper bound for the number of false positives selected by the GLM lasso in terms of a compatibility constant and restricted upper eigenvalue.

We begin by showing that the restricted upper eigenvalue

[TABLE]

over the cone $\mathcal{C}\left(\xi,S\right)=\{\boldsymbol{b}\in R^{p+1}:\left\|\boldsymbol{b}_{S^{c}}\right\|_{1}\leq\xi\left\|\boldsymbol{b}_{S}\right\|_{1}\neq 0\}$ is bounded above with probability tending to one when $s$ is proportional to $s_{0}$. For large $n$ and small $\eta^{*}$, the term involving $A^{\prime\prime}$ is bounded almost surely by our Lipschitz assumptions since $\left\|\boldsymbol{b}\right\|_{1}=\left\|\boldsymbol{b}_{S}\right\|_{1}+\left\|\boldsymbol{b}_{S^{c}}\right\|_{1}\leq(1+\xi)\left\|\boldsymbol{b}_{S}\right\|_{1}$ so $|\boldsymbol{x}_{i}^{T}\boldsymbol{b}|\leq\left\|\boldsymbol{x}_{i}\right\|_{\infty}\left\|\boldsymbol{b}\right\|_{1}\leq K(1+\xi)\eta^{*}$ in the cone $\mathcal{C}\left(\xi,S\right)$. Since the intercept is included in $S$, and the eigenvalues of $E[\boldsymbol{x}\boldsymbol{x}^{T}]$ are bounded above, we can apply Lemma 15 of Loh and Wainwright, (2012) for zero-mean subgaussian random vectors to show that $\sup_{|B|\leq s,B\cap S=\emptyset}\lambda_{\max}\left(\frac{1}{n}\left[\boldsymbol{X}^{T}\boldsymbol{X}\right]_{B}\right)^{2}$ is bounded above by a constant with probability at least $1-c_{1}\exp(c_{2}n)$ for some constants $c_{1},c_{2}>0$ when $s\propto s_{0}$ and $s_{0}\log(p)/n=o(1)$. Using this upper bound for the upper restricted eigenvalue and the assumed lower bound for the compatibility constant, Theorem 19 of Huang and Zhang, (2012) implies that the number of false positives, and thus the selected model size, is at most of order $s_{0}$.

We now are ready to prove Theorem 1. We begin by providing probability bounds for the events assumed in Lemmas 1 and 2. First, van de Geer and Bühlmann, (2009) showed that the compatibility constant is bounded away from zero with probability at least $1-c_{1}\exp(-c_{2}n)$ for some constants $c_{1},c_{2}>0$.

Next, we seek to bound $\frac{1}{n}\sum_{i=1}^{n}x_{ij}^{2}$ from above and away from zero uniformly over $j=1,\ldots,p+1$. This is trivial for the intercept, and the rest of the covariates are each bounded almost surely with mean zero and a variance of at least $min_{j}\boldsymbol{\Sigma}_{jj}\geq\lambda_{min}(\boldsymbol{\Sigma})>0$. We apply Hoeffding’s inequality and a union bound to get

[TABLE]

Since $\log(p)/n={o}_{p}\left(1\right)$, for large $n$ we can fix a $t\in(0,\lambda_{min}(\boldsymbol{\Sigma}))$ and the probability bound will be of order $c_{1}\exp(-c_{2}n)$ for some positive constants $c_{1}$ and $c_{2}$.

Lastly, let $\mathcal{B}$ denote the union of events that the compatibility constant is bounded, the diagonal entries of $\boldsymbol{X}^{T}\boldsymbol{X}$ are bounded, and the upper restricted eigenvalue (see Lemma 2) is bounded. Since the probability bound on each of these events has the same order, by a union bound and rescaling of the constants involved we have $P(\mathcal{B}^{c})\leq c_{1}\exp(-c_{2}n)$ for some $c_{1},c_{2}>0$. Then for the event $\mathcal{A}=\{\widehat{S}\supseteq S_{0},|\widehat{S}|\leq ks_{0}\}$ with $k>0$ not depending on $n$ and some constants $c_{3},c_{4},c_{5},c_{6}>0$ we have

[TABLE]

for $n$ sufficiently large, since $\log(p)/n=o(1)$.

Appendix Appendix B Proof of Theorem 2: Asymptotic results for debiased lasso after single sample split

Suppose the sample of $n$ individuals is split into a subsample $\mathcal{D}_{1}$ of size $n_{1}=qn$ and a subsample $\mathcal{D}_{2}$ of size $n_{2}=(1-q)n$, where we assume $n_{1}$ is an integer and $q\in(0,1)$. We apply a model selection procedure to $\mathcal{D}_{1}$ such that the selected model $\widehat{S}$ satisfies the event $\mathcal{A}=\{\widehat{S}\supseteq S_{0},|\widehat{S}|\leq o(\sqrt{n})\}$ with probability tending to one as $n\rightarrow\infty$, such as the lasso. For the selected model $\widehat{S}$, which includes the pre-specified set $S$, let $\boldsymbol{\widetilde{\beta}}^{DL}_{\widehat{S}}$ denote the debiased lasso estimator of $\boldsymbol{\beta}^{0}_{\widehat{S}}$ based on the subsample $\mathcal{D}_{2}$, and similarly denote the asymptotic covariance estimator by $\boldsymbol{\widehat{\Theta}}_{\widehat{S}}$. We decompose

[TABLE]

where, conditional on any sequence of selected models $\widehat{S}$ independent of $\mathcal{D}_{2}$ and satisfying the event $\mathcal{A}$, $r_{n}\rightarrow_{p}0$ by the theoretical results of Xia et al., (2021). Hence for any $\varepsilon>0$,

[TABLE]

where $P\left(|r_{n}|\geq\varepsilon|\mathcal{A}\right)\rightarrow 0$ as $n\rightarrow\infty$ by Xia et al., (2021) as mentioned above, so $r_{n}={o}_{p}\left(1\right)$ unconditionally as $n\rightarrow\infty$. Asymptotic normality, conditional on $\mathcal{A}$, also follows from the low-dimensional results of Xia et al., (2021) using the Lindeberg-Feller Central Limit Theorem.

Appendix Appendix C Proof of Theorem 3: Asymptotic results for debiased lasso after multiple sample splits

First we apply the results from the proof of Theorem 2 to the $(p+1)$-vector $\boldsymbol{a}$ with $\left\|\boldsymbol{a}_{S}\right\|_{2}=\left\|\boldsymbol{a}\right\|_{2}=1$, i.e. a contrast of only the pre-specified covariates $S$. The single split debiased lasso estimator $\boldsymbol{\widetilde{\beta}}^{DL}_{b}$ satisfies

[TABLE]

where $\widehat{S}_{b}\supseteq S_{0}\cup S$ is the fitted model, $\boldsymbol{\widetilde{I}}_{\widehat{S}_{b}}$ is the $|\widehat{S}_{b}|\times(p+1)$ matrix such that $\boldsymbol{\widetilde{I}}_{\widehat{S}_{b}}\boldsymbol{b}=\boldsymbol{b}_{\widehat{S}_{b}}$ for any $(p+1)$-vector $\boldsymbol{b}$, and $v_{b,i}$ is the indicator variable equal to one for $i\in D_{2,b}$.

Now consider the multiple splitting estimator $\boldsymbol{\widehat{\beta}}$ obtained by averaging across all possible splits $(D_{1,b},D_{2,b})_{b=1}^{B}$, $B={n\choose n_{2}}$. Again letting $\widehat{S}_{b}$ denote the fitted model from data $\mathcal{Z}=(y_{i},\boldsymbol{x}_{i})_{i=1}^{n}$ split by sampling indicators $\boldsymbol{v}_{b}$, the multiple splitting estimator satisfies

[TABLE]

where the expectations are with respect to the random splits $b=1,\ldots,B$, conditional on the data.

To establish asymptotic normality, we first show that $r_{n}={o}_{p}\left(1\right)$. Let $\boldsymbol{\widetilde{v}}_{b}$ be another sampling indicator vector independent of $\boldsymbol{v}_{b}$, with fitted model $\widetilde{S}_{b}$. The expectations in the following decomposition are taken with respect to $\boldsymbol{v}_{b}$ and $\boldsymbol{\tilde{v}}_{b}$, i.e. over $b=1,\ldots,B$, conditional on the data:

[TABLE]

and, by independence, $E\left[\boldsymbol{a}_{\widetilde{S}_{b}}^{T}\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0},\widetilde{S}_{b}}^{-1}\boldsymbol{\widetilde{I}}_{\widetilde{S}_{b}}v_{b,i}|\mathcal{Z}\right]=E\left[\boldsymbol{a}_{\widetilde{S}_{b}}^{T}\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0},\widetilde{S}_{b}}^{-1}\boldsymbol{\widetilde{I}}_{\widetilde{S}_{b}}|\mathcal{Z}\right]P(v_{b,i}=1)$ with $P(v_{b,i}=1)=n_{2}/n=1-q$. Since this quantity does not depend on $i$, we can write

[TABLE]

with h_{i,n}=\left[E\left(\boldsymbol{a}_{\widehat{S}_{b}}^{T}\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0},\widehat{S}_{b}}^{-1}\boldsymbol{\widetilde{I}}_{\widehat{S}_{b}}v_{b,i}\Big{|}\mathcal{Z}\right)-E\left(\boldsymbol{a}_{\widetilde{S}_{b}}^{T}\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0},\widetilde{S}_{b}}^{-1}\boldsymbol{\widetilde{I}}_{\widetilde{S}_{b}}v_{b,i}\Big{|}\mathcal{Z}\right)\right]\dot{\rho}_{\boldsymbol{\beta}^{0}}(y_{i},\boldsymbol{x}_{i}) so $r_{n,2}={o}_{p}\left(1\right)$ by Assumption (6). Also, since $P_{n}\dot{\rho}_{\boldsymbol{\beta}^{0}}$ has bounded entries with mean zero, we apply Hoeffding’s inequality and a union bound to obtain

[TABLE]

for some constant $c>0$ and $t>0$, so $\left\|\sqrt{n}P_{n}\dot{\rho}_{\boldsymbol{\beta}^{0}}\right\|_{\infty}=\mathcal{O}_{p}\left(\sqrt{\log(p)}\right)$ and, applying Assumption (7), we have

[TABLE]

In summary,

[TABLE]

so we have asymptotic normality by the Lindeberg-Feller Central Limit Theorem. Note that

[TABLE]

since $\boldsymbol{a}$ is a unit vector, $\boldsymbol{\Sigma}^{-1}_{\boldsymbol{\beta}^{0}}$ has bounded eigenvalues, and $|\widehat{S}_{b}|\leq ks_{0}$ with probability tending to one. Thus the Lindeberg condition can be verified by standard arguments for bounded covariates.

Appendix Appendix D Further discussion of the multiple splitting assumptions

The asymptotic representation used in our multiple splitting proof is distinct from the one used by Fei and Li, (2021) who, in their theoretical arguments, applied a different decomposition of $\sqrt{n}\boldsymbol{a}_{S}^{T}(\boldsymbol{\widehat{\beta}}_{S}-\boldsymbol{\beta}_{S}^{0})$ based on comparing the single split estimators with their respective oracle estimators. Under their partial orthogonality condition, which requires that the signal variables are independent of the noise variables and thus that $\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0}}$ is block diagonal (reordering the covariates without loss of generality), they argue that $\boldsymbol{\widehat{\beta}}_{S}$ has the same asymptotic normal distribution as the oracle estimator fit on $S\cup S_{0}$. Under different assumptions than Fei and Li, (2021) we also showed that $\sqrt{n}\boldsymbol{a}_{S}^{T}(\boldsymbol{\widehat{\beta}}_{S}-\boldsymbol{\beta}_{S}^{0})$ is asymptotically normal, but with variance $\boldsymbol{\eta}_{n}^{T}\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0}}\boldsymbol{\eta}_{n}$ that is possibly larger than the oracle variance $\boldsymbol{a}_{S\cup S_{0}}^{T}\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0},S\cup S_{0}}^{-1}\boldsymbol{a}_{S\cup S_{0}}$. Our assumptions instead concern the limiting behavior of the average (across all random splits that exclude $(y_{i},\boldsymbol{x}_{i})$ from model selection, i.e. such that $v_{b,i}=1$) $E\left[\boldsymbol{a}_{\widehat{S}_{b}}^{T}\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0},\widehat{S}_{b}}^{-1}\boldsymbol{\widetilde{I}}_{\widehat{S}_{b}}v_{b,i}|\mathcal{Z}\right]$, for $i=1,\ldots,n$. Note that, for each $i$, the corresponding average is independent of $(y_{i},\boldsymbol{x}_{i})$. Since the samples are i.i.d., the effect of always excluding the $i$th observation specifically from model selection should be asymptotically negligible (Assumption 6), and furthermore each average should converge in probability to a bounded random vector $\boldsymbol{\eta}_{n}^{T}$, which is independent of the data, with moderate rate $\left\|E\left[\boldsymbol{a}_{\widehat{S}_{b}}^{T}\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0},\widehat{S}_{b}}^{-1}\boldsymbol{\widetilde{I}}_{\widehat{S}_{b}}|\mathcal{Z}\right]-\boldsymbol{\eta}_{n}^{T}\right\|_{1}={o}_{p}\left(1/\sqrt{\log(p)}\right)$ (Assumption 7).

We can achieve the oracle-level result of Fei and Li, (2021) under either of two stronger assumptions that also imply our convergence conditions in Assumptions (6) and (7). First, under model selection consistency we have $\widehat{S}_{b}=S\cup S_{0}$ for all random splits, so $\boldsymbol{a}_{\widehat{S}_{b}}^{T}\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0},\widehat{S}_{b}}^{-1}\boldsymbol{\widetilde{I}}_{\widehat{S}_{b}}=\boldsymbol{a}_{S\cup S_{0}}^{T}\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0},S\cup S_{0}}^{-1}\boldsymbol{\widetilde{I}}_{S\cup S_{0}}$, a deterministic vector, for all $b=1,\ldots,B$ and we can directly set $\boldsymbol{\eta}_{n}^{T}=E\left[\boldsymbol{a}_{\widehat{S}_{b}}^{T}\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0},\widehat{S}_{b}}^{-1}\boldsymbol{\widetilde{I}}_{\widehat{S}_{b}}|\mathcal{Z}\right]$. By (Appendix C) we get $\sqrt{n}\boldsymbol{a}_{S}^{T}(\boldsymbol{\widehat{\beta}}_{S}-\boldsymbol{\beta}_{S}^{0})=\sqrt{n}\boldsymbol{\eta}_{n}^{T}P_{n}\dot{\rho}_{\boldsymbol{\beta}^{0}}+{o}_{p}\left(1\right)$. Then, using the fact that $\boldsymbol{\widetilde{I}}_{S\cup S_{0}}\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0}}\boldsymbol{\widetilde{I}}_{S\cup S_{0}}^{T}=\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0},S\cup S_{0}}$, the asymptotic variance becomes $\boldsymbol{\eta}_{n}^{T}\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0}}\boldsymbol{\eta}_{n}=\boldsymbol{a}_{S\cup S_{0}}^{T}\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0},S\cup S_{0}}^{-1}\boldsymbol{a}_{S\cup S_{0}}$. Alternatively, if the covariates $S\cup S_{0}$ are independent of $(S\cup S_{0})^{C}$ then $\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0}}=E[A^{\prime\prime}(\boldsymbol{x}_{i}^{T}\boldsymbol{\beta}^{0})\boldsymbol{x}_{i}\boldsymbol{x}_{i}^{T}]$ has a block diagonal structure that also implies $\boldsymbol{a}_{\widehat{S}_{b}}^{T}\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0},\widehat{S}_{b}}^{-1}\boldsymbol{\widetilde{I}}_{\widehat{S}_{b}}=\boldsymbol{a}_{S\cup S_{0}}^{T}\boldsymbol{\Sigma}_{\boldsymbol{\beta}^{0},S\cup S_{0}}^{-1}\boldsymbol{\widetilde{I}}_{S\cup S_{0}}$ for all $b=1,\ldots,B$, using the fact that $\boldsymbol{a}_{(S\cup S_{0})^{C}}$ is a zero-vector. Therefore this independence assumption also implies oracle-level variance. Note that, in general, this independence assumption is difficult to justify in practice when the true signal variables are unknown. An exception is in linear models, where it is sufficient to assume the pre-specified covariates $S$ are independent of $S^{C}$, e.g. by randomized treatment assignment, since $A^{\prime\prime}(\cdot)$ is a constant. In this sense our Assumptions (6) and (7), adapted from Wang et al., (2020) for linear models, are less restrictive than those used by others to justify multiple splitting procedures such as Fei and Li, (2021) for GLMs and Fei et al., (2019) for linear models.

Appendix Appendix E Additional simulation results

Table 5 summarizes results from the same simulation setting as the post-model selection procedure of Table 3, but with a larger sample size of $n=2000$, and $p=2800$ covariates. We see that the selection rate for small coefficients improves, resulting in lower bias and better confidence interval coverage, and the debiased lasso estimator still has a substantial advantage over the MLE.

We also present additional simulation results that have stronger correlation between covariates. Table 6 summarizes results from an AR(1) correlation structure with parameter 0.8, and Table 7 summarizes results from a compound symmetry correlation structure with parameter 0.8. We see that the multiple splitting debiased lasso estimator is fairly robust to high correlations between covariates, even as the model selection performance deteriorates. The debiased lasso continues to outperform the MLE, having lower bias and variance as well as better confidence interval coverage.