Multiply Robust Inference of Average Treatment Effects by High-dimensional Empirical Likelihood

TL;DR

This paper introduces a high-dimensional empirical likelihood method for robustly estimating the average treatment effect, accommodating multiple propensity score models and clustered data, with theoretical guarantees and practical validation.

Contribution

It develops a novel multiply robust inference procedure for ATE in high-dimensional settings, allowing for multiple propensity score models and clustered data, with asymptotic validity.

Findings

Method achieves valid confidence intervals under high-dimensional covariates.

Outperforms existing doubly robust methods in simulations.

Demonstrated effectiveness on real medical data.

Abstract

In this paper, we develop a multiply robust inference procedure of the average treatment effect (ATE) for data with high-dimensional covariates. We consider the case where it is difficult to correctly specify a single parametric model for the propensity scores (PS). For example, the target population is formed from heterogeneous sources with different treatment assignment mechanisms. We propose a novel high-dimensional empirical likelihood weighting method under soft covariate balancing constraints to combine multiple working PS models. An extended set of calibration functions is used, and a regularized augmented outcome regression is developed to correct the bias due to non-exact covariate balancing. Those two approaches provide a new way to construct the Neyman orthogonal score of the ATE. The proposed confidence interval for the ATE achieves asymptotically valid nominal coverage under high-dimensional covariates if any of the PS models, their linear combination, or the outcome regression model is correctly specified. The proposed method is extended to generalized linear models for the outcome variable. Specifically, we consider estimating the ATE for data with unknown clusters, where multiple working PS models can be fitted based on the estimated clusters. Our proposed approach enables robust inference of the ATE for clustered data. We demonstrate the advantages of the proposed approach over the existing doubly robust inference methods under high-dimensional covariates via simulation studies. We analyzed the right heart catheterization dataset, initially collected from five medical centers and two different phases of studies, to demonstrate the effectiveness of the proposed method in practice.

Figures

Tables

Table 1. Table 1 : Bias, standard error (Std), root-mean-squared error (RMSE), coverage probability of 95 % 95\% confidence intervals (coverage), and length of 95 % 95\% confidence intervals (CI width) of the proposed method DELW and five existing methods for the average treatment effect under n = 500 n=500 and the PS1 design for propensity scores.

		DELW	CBPS	RC	RB	rAIPW	DS
$p$	Measures	OR1 design for outcome regression
1000	Bias	0.0054	0.0017	-0.0130	0.0033	-0.0022	-0.0025
	Std	0.0967	0.0981	0.1002	0.1096	0.0964	0.1242
	RMSE	0.0968	0.0982	0.1011	0.1097	0.0965	0.1243
	Coverage	0.9500	0.9450	0.9200	0.9600	0.9450	0.9550
	CI width	0.3792	0.3873	0.3659	0.4318	0.3663	0.5182
2000	Bias	0.0039	0.0035	-0.0140	0.0042	-0.0059	-0.0061
	Std	0.0991	0.1009	0.1009	0.1099	0.0964	0.1402
	RMSE	0.0991	0.1010	0.1020	0.1100	0.0965	0.1404
	Coverage	0.9300	0.9350	0.8900	0.9500	0.9200	0.9700
	CI width	0.3761	0.3871	0.3655	0.4278	0.3650	0.5746
		OR2 design for outcome regression
1000	Bias	0.0094	0.0176	0.0191	0.0512	0.0191	0.0543
	Std	0.2027	0.4258	0.4350	0.5148	0.4511	0.4913
	RMSE	0.2029	0.4262	0.4355	0.5173	0.4514	0.4942
	Coverage	0.9200	0.9350	0.9200	0.9400	0.9050	0.9600
	CI width	0.6830	1.4209	1.3134	1.7057	1.3140	1.6734
2000	Bias	0.0105	-0.0246	-0.0167	-0.0008	-0.0334	-0.1299
	Std	0.2239	0.4193	0.4497	0.6513	0.5167	1.0984
	RMSE	0.2241	0.4201	0.4500	0.6513	0.5178	1.1061
	Coverage	0.9050	0.9150	0.8800	0.9300	0.8850	0.9750
	CI width	0.6714	1.3366	1.3042	1.7164	1.3058	2.2246

Table 2. Table 2 : Bias, standard error (Std), root-mean-squared error (RMSE), coverage probability of 95 % 95\% confidence intervals (coverage), and length of 95 % 95\% confidence intervals (CI width) of the proposed method DELW and five existing methods for the average treatment effect under n = 500 n=500 and the PS2 design for propensity scores.

		DELW	CBPS	RC	RB	rAIPW	DS
$p$	Measures	OR1 design for outcome regression
1000	Bias	-0.0111	-0.0060	-0.0149	-0.0072	-0.0092	-0.0066
	Std	0.0976	0.0961	0.1019	0.1064	0.0980	0.1230
	RMSE	0.0980	0.0961	0.1027	0.1064	0.0981	0.1228
	Coverage	0.9400	0.9600	0.9200	0.9500	0.9350	0.9400
	CI width	0.3760	0.3854	0.3692	0.4253	0.3692	0.4959
2000	Bias	-0.0135	-0.0103	-0.0157	-0.0142	-0.0129	-0.0270
	Std	0.1003	0.1012	0.1016	0.1099	0.1006	0.1300
	RMSE	0.1010	0.1015	0.1026	0.1106	0.1012	0.1324
	Coverage	0.9250	0.9400	0.9250	0.9250	0.9200	0.9600
	CI width	0.3734	0.3825	0.3669	0.4182	0.3673	0.5445
		OR2 design for outcome regression
1000	Bias	0.0133	0.2239	0.2056	0.1925	0.2220	0.1952
	Std	0.1832	0.3637	0.3529	0.3281	0.3265	0.3952
	RMSE	0.1832	0.4263	0.4077	0.3797	0.3942	0.4398
	Coverage	0.9300	0.8200	0.7850	0.8300	0.7850	0.9200
	CI width	0.6448	1.1885	1.0534	1.2984	1.0554	1.4325
2000	Bias	0.0024	0.2197	0.1949	0.1888	0.2033	0.1916
	Std	0.1701	0.3069	0.3084	0.3322	0.2957	0.3797
	RMSE	0.1697	0.3768	0.3642	0.3814	0.3582	0.4245
	Coverage	0.9150	0.7900	0.7600	0.8800	0.7650	0.9450
	CI width	0.6220	1.1491	1.0182	1.2987	1.0190	1.5172

Table 3. Table 3 : Bias, standard error (Std), root-mean-squared error (RMSE), coverage probability of 95 % 95\% confidence intervals (coverage), and length of 95 % 95\% confidence intervals (CI width) of the proposed method DELW and five existing methods for the average treatment effect under the clustered data with sample size n = 1000 n=1000 and dimension p = 1000 p=1000 .

	DELW	CBPS	RC	RB	rAIPW	DS
Measures	OR1 design for outcome model
Bias	0.0049	0.0028	0.0025	0.0013	0.0052	-0.0035
Std	0.0706	0.0721	0.0709	0.0740	0.0711	0.0725
RMSE	0.0708	0.0721	0.0709	0.0740	0.0713	0.0726
Coverage	0.9400	0.9300	0.9350	0.9550	0.9400	0.9550
CI width	0.2705	0.2746	0.2856	0.2873	0.2858	0.2847
	OR2 design for outcome model
Bias	0.0113	0.0204	0.0353	0.0402	0.0357	0.0182
Std	0.1172	0.2950	0.2927	0.3358	0.3050	0.3491
RMSE	0.1178	0.2957	0.2948	0.3382	0.3071	0.3496
Coverage	0.9500	0.9150	0.8900	0.8950	0.9000	0.9200
CI width	0.4505	0.9596	0.9565	1.0913	0.9569	0.9838

Table 4. Table 4 : Estimates ± 2 × \pm\ 2\ \times standard errors of one-month survival probabilities (in treatment and control groups) and the ATE by the proposed DELW, RC, and rAIPW methods.

Method	DELW	RC	rAIPW
${\mu }_1$	0.655$\pm$0.022	0.635$\pm$0.021	0.636$\pm$0.021
${\mu }_0$	0.678$\pm$0.016	0.688$\pm$0.016	0.691$\pm$0.016
ATE $\tau$	-0.022$\pm$0.026	-0.053$\pm$0.025	-0.055$\pm$0.025

Equations

$$\tau :={\mu }_1-{\mu }_0=\mathbb{E}\{Y(1)-Y(0)\}.$$

$$\pi (\mathbf{X})=\sum _{k=1}^q{d}_k{\pi }_k(\mathbf{X},{\mathbf{\gamma }}_k)$$

$$\mathbb{E}\{Y(d)\mid \mathbf{X}\}\subset \text{mbox}span\{b_0(\mathbf{X}),\dots ,b_{p_0}(\mathbf{X})\},\text{mbox}where{b}_0(\mathbf{X})=1,$$

$$\hat{\mathbf{b}}(\mathbf{X})=(b_1(\mathbf{X}),\dots ,b_{p_0}(\mathbf{X}),{\pi }_1^{(1)}(\mathbf{X},{\hat{\mathbf{\gamma }}}_1),\dots ,{\pi }_q^{(1)}(\mathbf{X},{\hat{\mathbf{\gamma }}}_q){)}^{\mathrm{T}}$$

$$\{{\tilde{p}}_i:D_i=1\}$$

$$n{\tilde{p}}_i=\{{\tilde{\lambda }}_0+{\tilde{\mathbf{\lambda }}}_1^{\mathrm{T}}({\hat{\mathbf{\pi }}}_i-n_1/n{1}_q)+{\tilde{\mathbf{\lambda }}}_2^{\mathrm{T}}{\tilde{\mathbf{b}}}_i{\}}^{-1},$$

$$\underset{{\lambda }_0,{\mathbf{\lambda }}_1,{\mathbf{\lambda }}_2}{\mathrm{argmin}}\frac{1}{n}\sum _{i=1}^n-D_i\log \{{\lambda }_0+{\mathbf{\lambda }}_1^{\mathrm{T}}({\hat{\mathbf{\pi }}}_i-n_1/n{1}_q)+{\mathbf{\lambda }}_2^{\mathrm{T}}{\tilde{\mathbf{b}}}_i\}+{\lambda }_0+{\omega }_{\mathrm{ps}}\Vert {\mathbf{\lambda }}_2{\Vert }_1.$$

$$\underset{\beta_{0},\boldsymbol{\beta}_{1},\boldsymbol{\beta}_{2}}{\operatorname{argmin}}\frac{1}{n}\sum_{i=1}^{n}\frac{D_{i}}{\widetilde{\pi}_{i}^{2}}\{Y_{i}-\beta_{0}-\boldsymbol{\beta}_{1}^{{\mathrm{\scriptscriptstyle T}}}\widehat{{\mathbf{b}}}({\mathbf{X}}_{i})-\boldsymbol{\beta}_{2}^{{\mathrm{\scriptscriptstyle T}}}\widehat{\boldsymbol{\pi}}_{i}\}^{2}+\omega_{{\mathrm{\scriptstyle or}}}(\|\boldsymbol{\beta}_{1}\|_{1}+\|\boldsymbol{\beta}_{2}\|_{1})\ \ \$$

$${\hat{\mu }}_1=\frac{1}{n}\sum _{i=1}^n\frac{D_i{Y}_i}{{\tilde{\pi }}_i}-\frac{1}{n}\sum _{i=1}^n\frac{D_i-{\tilde{\pi }}_i}{{\tilde{\pi }}_i}\{{\tilde{\beta }}_0+{\tilde{\mathbf{\beta }}}_1^{\mathrm{T}}\hat{\mathbf{b}}({\mathbf{X}}_i)+{\tilde{\mathbf{\beta }}}_2^{\mathrm{T}}{\hat{\mathbf{\pi }}}_i\}.$$

$$\begin{array}{rl}{\hat{\mu }}_1=& \frac{1}{n}\sum _{i=1}^n\frac{D_i}{{\pi }_i}{Y}_i+\frac{1}{n}\sum _{i=1}^n\frac{{\pi }_i-D_i}{{\pi }_i}\{{\beta }_{0,\ast }+{\mathbf{\beta }}_{1,\ast }^{\mathrm{T}}{\mathbf{b}}_{\ast }({\mathbf{X}}_i)+{\mathbf{\beta }}_{2,\ast }^{\mathrm{T}}{\mathbf{\pi }}_{i,\ast }\}+o_p(n^{-1/2});\end{array}$$

$${\hat{\mu }}_1=\frac{1}{n}\sum _{i=1}^n\mathbb{E}(Y_i|{\mathbf{X}}_i,D_i=1)+\frac{1}{n}\sum _{i=1}^n\frac{D_i}{{\pi }_i^{\ast }}\{Y_i-\mathbb{E}(Y_i|{\mathbf{X}}_i,D_i=1)\}+o_p(n^{-1/2});$$

$${\hat{\mu }}_1=\frac{1}{n}\sum _{i=1}^n\mathbb{E}(Y_i|{\mathbf{X}}_i,D_i=1)+\frac{1}{n}\sum _{i=1}^n\frac{D_i}{{\pi }_i}\{Y_i-\mathbb{E}(Y_i|{\mathbf{X}}_i,D_i=1)\}+o_p(n^{-1/2}).$$

$$\widehat{V}_{1}=\frac{1}{n}\sum_{i=1}^{n}\bigg{[}\frac{D_{i}Y_{i}}{\widetilde{\pi}_{i}}-\frac{D_{i}-\widetilde{\pi}_{i}}{\widetilde{\pi}_{i}}\{\widetilde{\beta}_{0}+\widetilde{\boldsymbol{\beta}}_{1}^{{\mathrm{\scriptscriptstyle T}}}\widehat{{\mathbf{b}}}({\mathbf{X}}_{i})+\widetilde{\boldsymbol{\beta}}_{2}^{{\mathrm{\scriptscriptstyle T}}}\widehat{\boldsymbol{\pi}}_{i}\}-\widehat{\mu}_{1}\bigg{]}^{2},$$

$$f(Y(d)\mid{\mathbf{X}})=h(Y(d),\phi)\exp\big{(}[Y(d)\{\boldsymbol{\beta}(d)^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{b}}_{{\mathrm{\scriptstyle or}}}({\mathbf{X}})\}-g\{\boldsymbol{\beta}(d)^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{b}}_{{\mathrm{\scriptstyle or}}}({\mathbf{X}})\}]/a(\phi)\big{)}$$

$$\widehat{\boldsymbol{\beta}}_{ini}=\underset{\boldsymbol{\beta}\in\mathbb{R}^{p_{0}+1}}{{\rm argmin}}\bigg{(}-\frac{1}{n}\sum_{i=1}^{n}D_{i}[Y\{\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{b}}_{{\mathrm{\scriptstyle or}}}({\mathbf{X}})\}-g\{\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{b}}_{{\mathrm{\scriptstyle or}}}({\mathbf{X}})\}]+w_{{\mathrm{\scriptstyle or}}}^{(0)}\|\boldsymbol{\beta}\|_{1}\bigg{)},$$

$$\begin{split}\{\widetilde{p}_{i}:D_{i}=1\}=&\underset{p_{i}}{\operatorname{argmax}}\sum_{i=1}^{n}D_{i}\log(p_{i}),\text{ \ subject to \ }p_{i}\geq 0,\sum_{i=1}^{n}D_{i}p_{i}=1,\\ &\sum_{i=1}^{n}D_{i}p_{i}\{\pi_{k}({\mathbf{X}}_{i},\widehat{\boldsymbol{\gamma}}_{k})-n_{1}/n\}=0\text{ \ for $k=1,\dots,q$ and }\\ &\bigg{|}\sum_{i=1}^{n}D_{i}p_{i}\widetilde{b}_{ij}^{(g)}\bigg{|}\leq\omega_{{\mathrm{\scriptstyle ps}}}\text{ \ for $j=1,\dots,m$},\end{split}$$

$$\begin{split}\{\widetilde{\beta}_{0},\widetilde{\boldsymbol{\beta}}_{1},\widetilde{\boldsymbol{\beta}}_{2}\}=\underset{\beta_{0},\boldsymbol{\beta}_{1},\boldsymbol{\beta}_{2}}{{\rm argmin}}\bigg{(}&\frac{1}{n}\sum_{i=1}^{n}\frac{D_{i}}{\widetilde{\pi}_{i}^{2}}[g\{\beta_{0}+\boldsymbol{\beta}_{1}^{{\mathrm{\scriptscriptstyle T}}}\widehat{{\mathbf{b}}}({\mathbf{X}}_{i})+\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}_{2}\boldsymbol{\widehat{\pi}}_{i}\}\\ &-Y_{i}\{\beta_{0}+\boldsymbol{\beta}_{1}^{{\mathrm{\scriptscriptstyle T}}}\widehat{{\mathbf{b}}}({\mathbf{X}}_{i})+\boldsymbol{\beta}^{{\mathrm{\scriptscriptstyle T}}}_{2}\boldsymbol{\widehat{\pi}}_{i}\}]+w_{{\mathrm{\scriptstyle or}}}(\|\boldsymbol{\beta}_{1}\|_{1}+\|\boldsymbol{\beta}_{2}\|_{1})\bigg{)},\end{split}$$

$${\hat{\mu }}_1=\frac{1}{n}\sum _{i=1}^n\frac{D_i{Y}_i}{{\tilde{\pi }}_i}-\frac{1}{n}\sum _{i=1}^n\frac{D_i-{\tilde{\pi }}_i}{{\tilde{\pi }}_i}\{g^{\prime }({\tilde{\beta }}_0+{\tilde{\mathbf{\beta }}}_1^{\mathrm{T}}\hat{\mathbf{b}}({\mathbf{X}}_i)+{\tilde{\mathbf{\beta }}}_2^{\mathrm{T}}{\hat{\mathbf{\pi }}}_i)\}.$$

$$\widehat{V}_{1}=\frac{1}{n}\sum_{i=1}^{n}\bigg{[}\frac{D_{i}Y_{i}}{\widetilde{\pi}_{i}}-\frac{D_{i}-\widetilde{\pi}_{i}}{\widetilde{\pi}_{i}}g^{\prime}\{\widetilde{\beta}_{0}+\widetilde{\boldsymbol{\beta}}_{1}^{{\mathrm{\scriptscriptstyle T}}}\widehat{{\mathbf{b}}}({\mathbf{X}}_{i})+\widetilde{\boldsymbol{\beta}}_{2}^{{\mathrm{\scriptscriptstyle T}}}\widehat{\boldsymbol{\pi}}_{i}\}-\widehat{\mu}_{1}\bigg{]}^{2}$$

$$({\mathbf{X}}_1,\dots ,{\mathbf{X}}_n{)}^{\mathrm{T}}=\mathbf{U}\mathbf{M}+\mathbf{\epsilon },$$

$${\pi }_{h_1\dots h_v}({\mathbf{X}}_i,{\hat{\mathbf{\gamma }}}_{h_1\dots h_v})=\sum _{\ell =1}^v{\tilde{\pi }}_{h_{\ell }}({\mathbf{X}}_i,{\hat{\mathbf{\gamma }}}_{h_{\ell },\ell })1({\hat{U}}_{i\ell }=1)$$

$$\displaystyle\mathbb{P}(D=1\mid{\mathbf{X}})=\big{[}1+\exp\{-X_{1}+X_{2}/2-X_{3}/4-(X_{4}+X_{5}-X_{6})/10\}\big{]}^{-1},$$

$$Y(0)=1+0.291\sum _{j=5}^{10}{X}_j+{ϵ}_0\text{mbox}andY(1)=2+0.137\sum _{j=5}^8{X}_j+{ϵ}_1,$$

$${\hat{\mu }}_1\approx \frac{1}{n}\sum _{i=1}^n\frac{D_i{Y}_i}{{\tilde{\pi }}_i}-\frac{1}{n}\sum _{i=1}^n\frac{D_i-{\tilde{\pi }}_i}{{\tilde{\pi }}_i}\{{\beta }_{0,\ast }+{\mathbf{\beta }}_{1,\ast }^{\mathrm{T}}\hat{\mathbf{b}}({\mathbf{X}}_i)+{\mathbf{\beta }}_{2,\ast }^{\mathrm{T}}{\hat{\mathbf{\pi }}}_i\},$$

$$\widehat{\mu}_{1}\approx\frac{1}{n}\sum_{i=1}^{n}\bigg{[}{\mathbb{E}}(Y_{i}|{\mathbf{X}}_{i},D_{i}=1)+\frac{D_{i}}{\pi_{i}^{*}}\{Y_{i}-{\mathbb{E}}(Y_{i}|{\mathbf{X}}_{i},D_{i}=1)\}\bigg{]},$$

$${\hat{\mu }}_1$$

$$\frac{1}{n}\sum _{i=1}^n-\frac{D_i{\tilde{\mathbf{a}}}_i^{\mathrm{T}}}{{\pi }_1({\mathbf{X}}_i,{\hat{\mathbf{\gamma }}}_1{)}^2}(\tilde{\mathbf{\lambda }}-{\mathbf{\lambda }}_{\ast })\{Y_i-{\beta }_{0,\ast }-{\mathbf{\beta }}_{1,\ast }^{\mathrm{T}}\mathbf{b}({\mathbf{X}}_i)-{\mathbf{\beta }}_{2,\ast }^{\mathrm{T}}{\mathbf{\pi }}_{i,\ast }\}$$

$$+$$

Full text

Keywords: Calibration; Clustering; Constrained optimization; High dimensionality; Multiply robust inference; Selection bias.

1 Introduction

The estimation of the average treatment effect stands as one of the fundamental problems in causal inference. The concept of propensity score (Rosenbaum and Rubin,, 1983) is central in addressing and mitigating potential selection biases in non-randomized trials. Estimated propensity scores are used in various methods for estimating ATEs, including PS matching, the inverse probability weighted (IPW) estimator, also known as the Horvitz–Thompson estimator, and the augmented IPW (AIPW) estimator that incorporates the outcome regression (OR) model, thereby attaining the doubly robust and semi-parametrically efficient properties (Bang and Robins,, 2005). Covariate balancing, which uses the IPW formulation with the estimated PS under the covariate balancing constraints, is another way to estimate the ATEs (Imai and Ratkovic,, 2014). The AIPW and covariate balancing estimators are asymptotically equivalent under the fixed-dimensional setting, and their properties have been well studied in this regime. See, for example, Cao et al., (2009); Imai and Ratkovic, (2014).

The well-known doubly robust procedures only utilize a single PS and OR model. The two models are likely misspecified when the data are high-dimensional and from multiple heterogeneous sources with different treatment assignment mechanisms. For example, estimating ATEs based on combined datasets is common in clinical trials and medical studies. The right heart catheterization (RHC) dataset (Connors et al.,, 1996) was collected from five medical centers and two studies, a prospective observational study and an interventional study, to examine the effect of the RHC treatment on the survival rate of severely ill patients. However, the released data do not share the subgroup information. Given the unknown subgroup structure, a single PS model is unlikely to be adequate.

1.1 Related work

In fixed-dimensional settings, Han, (2014) proposed a multiply robust inference method by empirical likelihood (EL) that allows multiple working models for both the propensity score and outcome regression. Valid statistical inference can be obtained if any of those models are correctly specified. Chan and Yam, (2014) explored generalized empirical likelihood calibration estimators to achieve multiple robustness for missing data. Chen and Haziza, (2017) proposed a multiply robust imputation method in survey sampling settings. However, those methods do not consider high-dimensional covariates or data with clusters.

In high-dimensional settings, Belloni et al., (2014) and Farrell, (2015) introduced a double selection approach for estimating the ATE. Their approaches employ regularized estimators for both the PS and OR models in the AIPW formulation. Chernozhukov et al., (2018) developed a versatile and flexible framework for conducting inference on a target parameter under high-dimensional nuisance parameters. Their method utilizes the Neyman orthogonal score and data-splitting technique, which includes the AIPW estimator as a special case for estimating the ATE. These methods are doubly robust for estimation of the ATE. However, these methods do not yield doubly robust inference, as valid inference requires correct specification of both the PS and OR models. In a related effort to derive asymptotically valid confidence intervals for ATEs, Athey et al., (2018) required correct specification of a linear outcome model. In recent work on linear outcome models, Ning et al., (2020); Tan, 2020a proposed doubly robust inference methods for ATEs under high-dimensional data by regularized calibrated estimation. Their proposed confidence intervals of ATEs achieve an asymptotically valid coverage if either the PS or OR model is correctly specified. Particularly, Ning et al., (2020) also considered generalized linear models (GLM) for the outcome variable.

All the aforementioned high-dimensional methods do not utilize multiple PS models. The performance of doubly robust estimators can deteriorate, even with slight misspecification in both the PS and OR models. The phenomenon was highlighted in a comprehensive set of numerical studies by Kang and Schafer, (2007), which demonstrated that, in certain scenarios, using two incorrectly specified models does not yield better results than using a single misspecified model. Therefore, it is important to develop multiply robust inference methods under high-dimensional data, especially when we lack confidence in using a single propensity score model to capture complex treatment assignment mechanisms.

1.2 Our contributions

In this paper, we develop a novel method for estimating ATEs that utilizes multiple PS models with high-dimensional covariates. The proposed method achieves multiply robust inference under high-dimensional data; that is, it provides asymptotically valid statistical inference for ATEs if any working PS model, their linear combination, or the outcome regression model is correctly specified. Moreover, the proposed estimator achieves semiparametric efficiency if both the OR model and one of the working PS models (or their linear combination) are correctly specified. The multiply robust inference property is important for clustered data from heterogeneous sources, as a single PS model may not be adequate to reflect the true treatment assignment mechanism. Combining multiple working PS models based on the estimated clusters provides more guarantees of achieving valid statistical inference of the ATE against misspecified clustering algorithms or PS models under a heterogeneous population with subgroup structures.

Our approach comprises high-dimensional empirical likelihood weighting, regularized augmented outcome regression, and high-dimensional clustering. We develop a novel soft covariate balancing procedure by high-dimensional EL to combine estimated propensity scores from multiple working PS models. Our soft calibration EL is related to the regularized EL estimation in Chang et al., (2018, 2021). However, the main goals are different. The regularized EL approach is constructed to estimate high-dimensional estimating equations, while our method is proposed to estimate the propensity scores that achieve approximate covariate balancing. Note that exact covariate balancing is not possible under high-dimensional covariates. To achieve the multiply robust inference property, we apply soft calibration to derivatives of all working PS models and fit augmented outcome regression with regularization to correct the bias due to non-exact covariate balancing. These two approaches provide a new way to construct the Neyman orthogonal score of the ATE, ensuring that the estimation variability of all nuisance parameters does not affect the influence function of the estimated ATE. In this sense, our proposed method is also related to the existing literature on high-dimensional inference of regression models (Zhang and Zhang,, 2014; van de Geer et al.,, 2014). The proposed procedure is also extended to GLMs for the outcome variable. Simulation studies demonstrate the superiority of the proposed method over the existing methods.

1.3 Organization and notations

This paper is organized as follows. Section 2 introduces the problem formulation of estimating ATEs under high-dimensional data and the setting of multiple robust inferences. Section 3 presents the proposed multiply robust procedure under high-dimensional covariates. Section 4 delves into the theoretical properties. Section 5 extends the proposed approach to encompass GLMs for the outcome regression. Section 6 considers the estimation of ATEs under clustered data. Section 7 conducts simulation studies to verify the proposed method and compare it with the existing doubly robust procedures under high dimensional data. Section 8 analyzes the treatment effect of RHC using a real data set with unknown clusters. Section 9 summarizes the paper. A simplified proof to explain the multiply robust property of the proposed method is presented in the Appendix, which explains the importance of soft calibration using augmented covariates and augmented outcome regression to achieve multiply robust inferences under high-dimensional covariates. The theoretical results under GLM, additional numerical results, and all technical proofs are relegated to the supplementary material (SM).

For a $p$-dimensional vector ${\mathbf{x}}=(x_{1},\dots,x_{p})^{{\mathrm{\scriptscriptstyle T}}}$, let $\|{\mathbf{x}}\|_{q}=(\sum_{j=1}^{p}|x|_{j}^{q})^{1/q}$ and $\|{\mathbf{x}}\|_{0}=\sum_{j=1}^{p}1(x_{j}\neq 0)$ be the $\ell_{q}$ and $\ell_{0}$ norms of ${\mathbf{x}}$, respectively, where $1(\cdot)$ denotes the indicator function. For a $p\times p$ symmetric positive semidefinite matrix ${\mathbf{A}}$, let $\lambda_{\text{max}}({\mathbf{A}})$ and $\lambda_{\text{min}}({\mathbf{A}})$ be the maximal and minimal eigenvalues of ${\mathbf{A}}$, respectively. For two sequences of random variables, $\{a_{n}\}$ and $\{b_{n}\}$, $a_{n}\approx b_{n}$ and $a_{n}=o_{p}(b_{n})$ denote $|a_{n}-b_{n}|$ and $a_{n}/b_{n}$ converging to zero in probability, respectively. For two sequences of positive real numbers, $\{\tilde{a}_{n}\}$ and $\{\tilde{b}_{n}\}$, we use $\tilde{a}_{n}\asymp\tilde{b}_{n}$ to indicate that $\tilde{c}_{1}\leq a_{n}/b_{n}\leq\tilde{c}_{2}$ for two positive constants $\tilde{c}_{1}$ and $\tilde{c}_{2}$.

2 Problem and Settings

Let ${\mathbf{X}}=(X_{1},\ldots,X_{p})^{{\mathrm{\scriptscriptstyle T}}}$ be the $p$-dimensional covariate and $D$ be the indicator of treatment assignment, taking the value $1$ if the observation receives the treatment and [math] otherwise. Let $Y(d)$ be the potential outcome under $D=d$ for $d=0,1$, and $Y=DY(1)+(1-D)Y(0)$ be the observed response variable. Let $\mathcal{G}$ be the joint distribution of $(D,{\mathbf{X}}^{{\mathrm{\scriptscriptstyle T}}},Y(0),Y(1))$. Suppose the observed samples $\{(D_{i},{\mathbf{X}}_{i}^{{\mathrm{\scriptscriptstyle T}}},Y_{i}):i=1,\ldots,n\}$ are independent and identically distributed (IID) realizations from the joint distribution of $(D,{\mathbf{X}}^{{\mathrm{\scriptscriptstyle T}}},Y)$. Let $n_{1}=\sum_{i=1}^{n}D_{i}$ and $n_{0}=n-n_{1}$ be the sample sizes of the treated and control groups, respectively. We consider the case of high-dimensional covariates where $p$ can be much larger than $n$. Let $\mu_{1}:={\mathbb{E}}\{Y(1)\}$ and $\mu_{0}:={\mathbb{E}}\{Y(0)\}$, where the expectation is taken with respect to the population distribution $\mathcal{G}$. We are interested in estimating the average treatment effect:

[TABLE]

To identify both $\mu_{0}$ and $\mu_{1}$, we assume that the potential outcomes and treatment assignment are conditionally independent given the covariates, as formalized in the following unconfoundedness assumption.

Condition 1** (Unconfoundedness).**

$Y(1)$* and $Y(0)$ are conditionally independent of $D$ given ${\mathbf{X}}$ such that $Y(d)\perp D\mid{\mathbf{X}}$ for $d=0,1$.*

Under Condition 1, the parameter $\mu_{d}$ can be identified by either the outcome regression of $\mathbb{E}(Y\mid{\mathbf{X}},D=d)$ or inverse weighting of the propensity score $\pi({\mathbf{X}}):=\mathbb{P}(D=1\mid{\mathbf{X}})$. The doubly robust inference methods (Tan, 2020a, ; Ning et al.,, 2020) provide valid statistical inference of ATE when either the OR or PS model is correctly specified. However, specifying a single correct PS model may be challenging, especially when the target population contains multiple heterogeneous subgroups. We consider the problem of multiply robust inference (Han,, 2014; Li et al.,, 2020) in high-dimensional settings, where the true PS model is assumed to lie within the convex hull of multiple candidate models $\{\pi_{k}({\mathbf{X}},\boldsymbol{\gamma}_{k}):k=1,\ldots,q\}$ such that

[TABLE]

for unknown parameters $\{\boldsymbol{\gamma}_{k}\}_{k=1}^{q}$ and weights $\{d_{k}\}_{k=1}^{q}$ where $d_{k}\geq 0$. A mixture model of conditional density, similar to the setting in (2.2), is considered for distributionally robust prediction (Wang et al.,, 2023). Let ${\mathbf{b}}_{{\mathrm{\scriptstyle or}}}({\mathbf{X}})=(b_{0}({\mathbf{X}}),\ldots,b_{p_{0}}({\mathbf{X}}))^{{\mathrm{\scriptscriptstyle T}}}$ be a set of $p_{0}$ functions of ${\mathbf{X}}$, where $p_{0}$ could be larger than $p$ and much larger than $n$. We consider a working linear OR model for ${\mathbb{E}}\{Y(d)\mid{\mathbf{X}}\}$ such that

[TABLE]

where $\mbox{span}\{b_{0}({\mathbf{X}}),\ldots,b_{p_{0}}({\mathbf{X}})\}$ denotes the linear space spanned by ${\mathbf{b}}_{{\mathrm{\scriptstyle or}}}({\mathbf{X}})$. A multiply robust procedure provides asymptotically valid statistical inference of the ATE if either the OR model in (2.3) is correctly specified or any one of the PS models $\{\pi_{k}({\mathbf{X}},\boldsymbol{\gamma}_{k})\}_{k=1}^{q}$ or their linear combination in (2.2) is correctly specified.

The multiply robust procedure is more advantageous than the doubly robust procedure, which only use a single PS model. In many scientific studies, datasets often combine studies from multiple sources. See the right heart catheterization dataset (Connors et al.,, 1996) analyzed in Section 8. The multiply robust property is beneficial particularly when the target population comes from multiple heterogeneous sources and the origin of each observation is unknown. Section 6 introduces the ATE estimation under clustered data, where combining multiple propensity score estimates based on clustering results enables valid statistical inference. The existing high-dimensional doubly robust inference methods cannot be applied to multiple PS models, and the existing multiply robust inference methods are designed only for fixed-dimensional settings.

3 Methodology

We present the proposed multiply robust method to estimate $\mu_{1}={\mathbb{E}}\{Y(1)\}$ under high-dimensional covariates. The estimator of $\mu_{0}={\mathbb{E}}\{Y(0)\}$ can be similarly constructed. Let $\pi_{k}^{(1)}({\mathbf{X}},\boldsymbol{\gamma}_{k})=\partial\pi_{k}({\mathbf{X}},\boldsymbol{\gamma}_{k})/\partial\boldsymbol{\gamma}$ be the first-order derivative of the working PS model $\pi_{k}({\mathbf{X}},\boldsymbol{\gamma}_{k})$, and $\widehat{\boldsymbol{\gamma}}_{k}$ be an estimate of $\boldsymbol{\gamma}_{k}$ for $k=1,\dots,q$. Let

[TABLE]

be the augmented covariate functions, where we use $\widehat{b}_{j}({\mathbf{X}})$ to denote $j$th component of $\widehat{{\mathbf{b}}}({\mathbf{X}})$. The accent in $\widehat{{\mathbf{b}}}({\mathbf{X}})$ indicates the dependence on the estimates $\{\widehat{\boldsymbol{\gamma}}_{k}\}_{k=1}^{q}$. To keep the notations consistent, we also use $\widehat{b}_{j}({\mathbf{X}})$ to denote the first $p_{0}$ components of $\widehat{{\mathbf{b}}}({\mathbf{X}})$, even though they are free of estimates. Let $\overline{{\mathbf{b}}}=(\bar{b}_{1},\ldots,\bar{b}_{m})^{{\mathrm{\scriptscriptstyle T}}}=n^{-1}\sum_{i=1}^{n}\widehat{{\mathbf{b}}}({\mathbf{X}}_{i})$ be the overall sample mean and $\widetilde{{\mathbf{b}}}_{i}=(\widetilde{b}_{i1},\ldots,\widetilde{b}_{im})^{{\mathrm{\scriptscriptstyle T}}}=\widehat{{\mathbf{b}}}({\mathbf{X}}_{i})-\bar{{\mathbf{b}}}$. The treatment group empirical likelihood solves the constrained optimization problem:

[TABLE]

where (3.4) and (3.4) are called internal bias correction (IBC) constraints and soft calibration (soft covariate balancing) constraints, respectively. Here, soft calibration refers to the property that the weighted mean $\sum_{i}D_{i}\tilde{p}_{i}\widehat{b}_{j}({\mathbf{X}}_{i})$ of $\{\widehat{b}_{j}({\mathbf{X}}_{i})\}_{i=1}^{n}$ over the treatment group is approximately, though not necessarily exactly, equal to the overall sample mean $\bar{b}_{j}$ for $j=1,\dots,m$, and $\omega_{{\mathrm{\scriptstyle ps}}}$ is the tolerance parameter for non-exact covariate balancing. Calibrating the derivatives of the PS models is important for the multiply robust inference of the ATE. See the explanations from (A.1) to (A.5) for the asymptotic expansion of the proposed estimator under the misspecified outcome regression model. In the fixed-dimensional setting, Han, (2014) proposed EL weighting with exact covariate balancing. However, in the high-dimensional setting, exact balancing is generally infeasible, as no solution exists for the weights $\{\tilde{p}_{i}\}_{i=1}^{n}$ that simultaneously satisfy all exact balancing constraints.

Let $\widehat{\boldsymbol{\pi}}_{i}=(\pi_{1}({\mathbf{X}}_{i},\widehat{\boldsymbol{\gamma}}_{1}),\ldots,\pi_{q}({\mathbf{X}}_{i},\widehat{\boldsymbol{\gamma}}_{q}))^{{\mathrm{\scriptscriptstyle T}}}$ be the estimated propensity scores for individual $i$ from $q$ working PS models, and $\mathbf{1}_{q}=(1,\ldots,1)^{{\mathrm{\scriptscriptstyle T}}}$ be a $q$-dimensional vector of ones. The solutions of EL $\{\widetilde{p}_{i}:D_{i}=1\}$ for the treatment group take the form

[TABLE]

where $\widetilde{\lambda}_{0}$, $\widetilde{\boldsymbol{\lambda}}_{1}=(\widetilde{\lambda}_{11},\ldots,\widetilde{\lambda}_{1q})^{{\mathrm{\scriptscriptstyle T}}}$, $\widetilde{\boldsymbol{\lambda}}_{2}=(\widetilde{\lambda}_{21},\ldots,\widetilde{\lambda}_{2m})^{{\mathrm{\scriptscriptstyle T}}}$ are Lagrange multipliers such that the solutions $\{\widetilde{p}_{i}:D_{i}=1\}$ satisfy constraints in (3.4), (3.4) and (3.4). Using dual optimization, $(\widetilde{\lambda}_{0},\widetilde{\boldsymbol{\lambda}}_{1},\widetilde{\boldsymbol{\lambda}}_{2})$ can be obtained by the regularized unconstrained optimization problem:

[TABLE]

By the Karush–Kuhn–Tucker (KKT) conditions for problem (3.6), it can be shown that the EL weights $\{\widetilde{p}_{i}\}$, obtained by plugging in the minimizers of (3.6), satisfy the constraints in (3.4)–(3.4). The empirical likelihood weights $\{\widetilde{p}_{i}:D_{i}=0\}$ of the control group can be obtained in a similar manner.

Under the fixed-dimensional setting with exact covariate balancing, namely $\omega_{{\mathrm{\scriptstyle ps}}}=0$ in (3.4), the weighted mean $\sum_{i=1}^{n}D_{i}\widetilde{p}_{i}Y_{i}$ is an implicit AIPW estimator of $\mu_{1}$ (Han,, 2014). However, the equivalence no longer holds in the high-dimensional setting due to the soft calibration in (3.4). To correct the bias from the non-exact covariate balancing, we propose to estimate an augmented outcome regression (AOR) of $Y$ under the treatment group. Let $\{\widetilde{\beta}_{0},\widetilde{\boldsymbol{\beta}}_{1},\widetilde{\boldsymbol{\beta}}_{2}\}$ be the minimizer of the regularized weighted linear regression:

[TABLE]

with the weights $\{\widetilde{\pi}_{i}^{-2}\}_{i=1}^{n}$ on the treatment group, where $\widetilde{\pi}_{i}=(n\widetilde{p}_{i})^{-1}$ and $\omega_{{\mathrm{\scriptstyle or}}}$ is a penalty parameter. In addition to the covariates ${\mathbf{b}}_{{\mathrm{\scriptstyle or}}}({\mathbf{X}})$ under the working OR model in (2.3), the augmented outcome regression also includes the working PS model $\pi_{k}({\mathbf{X}}_{i},\widehat{\boldsymbol{\gamma}}_{k})$ and its derivatives $\pi_{k}^{(1)}({\mathbf{X}},\widehat{\boldsymbol{\gamma}}_{k})$ as additional covariates. The regression coefficient $\boldsymbol{\beta}_{1}$ is composed as $\boldsymbol{\beta}_{1}=(\boldsymbol{\beta}_{10}^{{\mathrm{\scriptscriptstyle T}}},\boldsymbol{\beta}_{11}^{{\mathrm{\scriptscriptstyle T}}},\ldots,\boldsymbol{\beta}_{1q}^{{\mathrm{\scriptscriptstyle T}}})^{{\mathrm{\scriptscriptstyle T}}}$, where $\boldsymbol{\beta}_{10}$ and $\boldsymbol{\beta}_{1k}$ are the coefficients of covariates $\{b_{j}({\mathbf{X}})\}_{j=1}^{p_{0}}$ and $\pi_{k}^{(1)}({\mathbf{X}},\widehat{\boldsymbol{\gamma}}_{k})$ for $k=1,\dots,q$, respectively. The proposed multiply robust estimator of $\mu_{1}$ using the AIPW formulation is

[TABLE]

The multiply robust estimator $\widehat{\mu}_{0}$ of $\mu_{0}$ can be similarly constructed. The proposed estimate of ATE is $\widehat{\tau}=\widehat{\mu}_{1}-\widehat{\mu}_{0}$.

The existing approaches for high-dimensional causal inference (Belloni et al.,, 2014; Ning et al.,, 2020; Tan, 2020a, ) performed calibration only for the covariates of the PS model and apply regularized regression only to the covariates of the OR model. However, such strategies are not sufficient to guarantee the multiple robustness of the estimator. Incorporating augmented covariates in both soft calibration and outcome regression induces Neyman orthogonality for all nuisance parameters, so that the statistical inference of the ATE is multiply robust against misspecified PS models $\{\pi_{k}({\mathbf{X}},\boldsymbol{\gamma}_{k})\}_{k=1}^{q}$ in (2.2) or misspecified OR model in (2.3). We explain the importance of calibrating derivatives of propensity scores and augmented outcome regression in the Appendix.

4 Theoretical Results

In this section, we provide conditions to show the multiply robust properties of the proposed estimator. Let $\{\boldsymbol{\gamma}_{k,\ast}\}_{k=1}^{q}$, $\lambda_{0,\ast}$, $\boldsymbol{\lambda}_{1,\ast}$, $\boldsymbol{\lambda}_{2,\ast}$, $\beta_{0,\ast}$, $\boldsymbol{\beta}_{1,\ast}$, $\boldsymbol{\beta}_{2,\ast}$ be the probability limits of $\{\widehat{\boldsymbol{\gamma}}_{k}\}_{k=1}^{q}$, $\widetilde{\lambda}_{0}$, $\widetilde{\boldsymbol{\lambda}}_{1}$, $\widetilde{\boldsymbol{\lambda}}_{2}$, $\widetilde{\beta}_{0}$, $\widetilde{\boldsymbol{\beta}}_{1}$, $\widetilde{\boldsymbol{\beta}}_{2}$, as $n,p\to\infty$, respectively, under either correctly specified or misspecified models. The rigorous definitions of these parameters are given in Section S3 of the SM. Let $\boldsymbol{\pi}_{i,\ast}=(\pi_{1}({\mathbf{X}}_{i},\boldsymbol{\gamma}_{1,\ast}),\ldots,\pi_{q}({\mathbf{X}}_{i},\boldsymbol{\gamma}_{q,\ast}))^{{\mathrm{\scriptscriptstyle T}}}$, ${\mathbf{b}}_{\ast}({\mathbf{X}}_{i})=(b_{1}({\mathbf{X}}),\ldots,b_{p_{0}}({\mathbf{X}}),\pi_{1}^{(1)}({\mathbf{X}}_{i},\boldsymbol{\gamma}_{1,*})^{{\mathrm{\scriptscriptstyle T}}},\dots,\pi_{q}^{(1)}({\mathbf{X}}_{i},\boldsymbol{\gamma}_{q,*})^{{\mathrm{\scriptscriptstyle T}}})^{{\mathrm{\scriptscriptstyle T}}}$ and $\boldsymbol{\lambda}_{\ast}=(\lambda_{0,\ast},\boldsymbol{\lambda}_{1,\ast}^{{\mathrm{\scriptscriptstyle T}}},\boldsymbol{\lambda}_{2,\ast}^{{\mathrm{\scriptscriptstyle T}}})^{{\mathrm{\scriptscriptstyle T}}}$. Let $s_{1}=\max_{1\leq k\leq q}\|\boldsymbol{\gamma}_{k,*}\|_{0}$ and $s_{2}=\|(\boldsymbol{\beta}_{1,*}^{{\mathrm{\scriptscriptstyle T}}},\boldsymbol{\beta}_{2,*}^{{\mathrm{\scriptscriptstyle T}}})\|_{0}$ denote the sparsity levels of the propensity score models and the outcome regression model, respectively. Let $\pi_{k,j}^{(1)}({\mathbf{X}}_{i},\boldsymbol{\gamma}_{k})$ denote the $j$th component of $\pi_{k}^{(1)}({\mathbf{X}}_{i},\boldsymbol{\gamma}_{k})$. We make the following regularity conditions.

Condition 2**.**

There exists a positive constant $c_{0}\in(0,1/2)$ such that the true propensity score $\pi_{i}=\pi({\mathbf{X}}_{i})$ satisfies $c_{0}\leq\pi_{i}\leq 1-c_{0}$ for $1\leq i\leq n$.

Condition 3**.**

The random variables ${\mathbf{b}}_{\ast}({\mathbf{X}})$ and $\epsilon_{1}=Y(1)-\beta_{0,*}-\boldsymbol{\beta}_{1,*}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{b}}_{\ast}({\mathbf{X}})-\boldsymbol{\beta}_{2,*}^{{\mathrm{\scriptscriptstyle T}}}\boldsymbol{\pi}_{*}$ are sub-Gaussian distributed. Namely, there is a positive constant $C$, such that for every $t\geq 0$ and ${\mathbf{v}}\in\mathbb{R}^{m}$ satisfying $\|\boldsymbol{v}\|_{2}=1$, we have $\mathbb{P}(|{\mathbf{v}}^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{b}}_{\ast}({\mathbf{X}})|\geq t)\leq 2\exp(-t^{2}/C^{2})$ and $\mathbb{P}(|\epsilon_{1}|\geq t)\leq 2\exp(-t^{2}/C^{2})$. When the outcome model is misspecified, the random variable $\epsilon_{1}$ is bounded.

Condition 4**.**

The minimum and maximum eigenvalues of $\boldsymbol{\Sigma}={\rm Cov}\{({\mathbf{b}}_{\ast}^{{\mathrm{\scriptscriptstyle T}}}({\mathbf{X}}),\boldsymbol{\pi}_{*})^{{\mathrm{\scriptscriptstyle T}}}\}$ satisfy $C\leq\lambda_{\text{min}}(\boldsymbol{\Sigma})\leq\lambda_{\text{max}}(\boldsymbol{\Sigma})\leq 1/C$ for a positive constant $C>0$.

Condition 5**.**

(i) The propensity score models satisfy $\mathbb{E}\{\pi_{k}({\mathbf{X}}_{i},\boldsymbol{\gamma}_{k,*})\}=\mathbb{P}(D=1)$, $c_{0}\leq\pi_{k}({\mathbf{X}}_{i},\boldsymbol{\gamma}_{k,*})\leq 1-c_{0}$ for a positive constant $c_{0}\in(0,1/2)$, and are Lipschitz continuous $|\pi_{k}({\mathbf{X}}_{i},\boldsymbol{\gamma}_{k})-\pi_{k}({\mathbf{X}}_{i},\boldsymbol{\gamma}_{k}^{\prime})|\leq C|{\mathbf{X}}_{i}^{\top}(\boldsymbol{\gamma}_{k,*}-\boldsymbol{\gamma}_{k}^{\prime})|$, $|\pi_{k,j}^{(1)}({\mathbf{X}}_{i},\boldsymbol{\gamma}_{k})-\pi_{k,j}^{(1)}({\mathbf{X}}_{i},\boldsymbol{\gamma}_{k}^{\prime})|\leq C|{\mathbf{X}}_{i}^{\top}(\boldsymbol{\gamma}_{k,*}-\boldsymbol{\gamma}_{k}^{\prime})|$ for all $j$ and $|{\mathbf{X}}_{i}^{\top}(\boldsymbol{\gamma}_{k,*}-\boldsymbol{\gamma}_{k}^{\prime})|=o_{p}(1)$. (ii) Estimates $\{\widehat{\boldsymbol{\gamma}}_{k}\}_{k=1}^{q}$ satisfy $\max_{k}\|\widehat{\boldsymbol{\gamma}}_{k}-\boldsymbol{\gamma}_{k,\ast}\|_{1}\leq Cs_{1}\{\log(p)/n\}^{1/2}$, $\max_{k}\|\widehat{\boldsymbol{\gamma}}_{k}-\boldsymbol{\gamma}_{k,\ast}\|_{2}\leq C\{s_{1}(\log p)/n\}^{1/2}$, $|(\boldsymbol{\gamma}_{k,*}-\widehat{\boldsymbol{\gamma}}_{k})^{\top}\pi_{k}^{(2)}({\mathbf{X}}_{i},\boldsymbol{\Delta})(\boldsymbol{\gamma}_{k,*}-\widehat{\boldsymbol{\gamma}}_{k})|\leq C|{\mathbf{X}}_{i}^{\top}(\boldsymbol{\gamma}_{k,*}-\widehat{\boldsymbol{\gamma}}_{k})|^{2}$ for $\boldsymbol{\Delta}\in[\boldsymbol{\gamma}_{k,*},\widehat{\boldsymbol{\gamma}}_{k}]$ with probability approaching to 1. (iii) If all the propensity score models are misspecified and (2.2) does not hold, then the population value $\boldsymbol{\lambda}_{\ast}$ of the optimization problem (3.6), defined in (S3.8) in the SM, satisfies $\|\boldsymbol{\lambda}_{2,*}\|_{0}=s_{\lambda}$, $s_{1}^{1/2}s_{\lambda}\log(m)=o(n^{1/2})$, and $c_{0}\leq\pi_{i}^{*}\leq 1-c_{0}$ for a positive constant $c_{0}\in(0,1/2)$.

Condition 6**.**

Let $s_{0}=\max\{s_{1},s_{2}\}$. Sparsity parameters satisfy $(s_{0}s_{1}s_{2})^{1/2}\log(m)=o(n^{1/2})$.

Condition 2 is the overlapping condition of the treatment and control groups, which is commonly required for causal inference. Conditions 3 and 4 are conventional assumptions in high-dimensional analysis (van de Geer et al.,, 2014; Farrell,, 2015), which are needed for concentration inequalities and the compatibility condition for high-dimensional regularized regression. The assumption of bounded errors $\epsilon_{1}$ under a misspecified outcome model can be relaxed to a sub-Gaussian error if a stronger sparsity condition, $(s_{0}s_{1}s_{2})^{1/2}\log(m)\log(n)=o(n^{1/2})$ with an additional factor $\log(n)$, is imposed in Condition 6. Condition 5 requires certain convergence rates of both $\ell_{1}$ and $\ell_{2}$ norms of the estimators $\{\widehat{\boldsymbol{\gamma}}_{k}\}$ for the working PS models. These convergence rates are well-established in high-dimensional statistics and are satisfied by the Lasso estimator for generalized linear models (van de Geer,, 2008) and the regularized calibration estimator (Tan, 2020b, ). The convergence of the prediction error of $\pi_{k}^{(1)}({\mathbf{X}},\widehat{\boldsymbol{\gamma}}_{k})$ can be achieved through the convergence of $\widehat{\boldsymbol{\gamma}}_{k}$ combined with a Taylor expansion. This condition can be verified for working PS models under high-dimensional clustered data in Section 6. The last part of Condition 5 becomes necessary when all the PS models are misspecified. In this case, we require that the probability limit $\boldsymbol{\lambda}_{\ast}$ from (3.6) to be sparse. Condition 6 represents the regularity condition on the sparsity levels. In Tan, 2020a and Ning et al., (2020), the condition $\max\{s_{1},s_{2}\}\log(m)/n^{1/2}=o(1)$ was assumed, which is slightly weaker than our Condition 6. The additional term $(\min\{s_{1},s_{2}\})^{1/2}$ in our condition arises from the feasibility of high-dimensional empirical likelihood (Chang et al.,, 2021). Moreover, the term $s_{0}^{1/2}$ is necessitated by the inclusion of the plug-in estimates $\pi_{k}^{(1)}({\mathbf{X}},\widehat{\boldsymbol{\gamma}}_{k})$ and $\pi_{k}({\mathbf{X}},\widehat{\boldsymbol{\gamma}}_{k})$ as covariates in the AOR defined in (3.7), which introduces additional estimation errors from $\widehat{\boldsymbol{\gamma}}_{k}$. Nonetheless, incorporating $\pi_{k}^{(1)}({\mathbf{X}},\widehat{\boldsymbol{\gamma}}_{k})$ and $\pi_{k}({\mathbf{X}},\widehat{\boldsymbol{\gamma}}_{k})$ in the AOR is critical for achieving multiply robust inference under a misspecified OR model as demonstrated in (A)–(A.5) in the Appendix.

In the following, we first present the theoretical properties of the proposed empirical likelihood weights $\{\widetilde{p}_{i}:D_{i}=1\}$ in (3.5).

Proposition 1**.**

Under Conditions 1–5, assuming $w_{{\mathrm{\scriptstyle ps}}}\asymp\{s_{1}\log(m)/n\}^{1/2}$ and $s_{1}\log(m)=o(n)$, the following results hold:

(i) if the PS model in (2.2) is correctly specified, then $(\widetilde{\lambda}_{0},\widetilde{\boldsymbol{\lambda}}_{1}^{{\mathrm{\scriptscriptstyle T}}},\widetilde{\boldsymbol{\lambda}}_{2}^{{\mathrm{\scriptscriptstyle T}}})\stackrel{{\scriptstyle p}}{{\to}}(\lambda_{0,\ast},\boldsymbol{\lambda}_{1,\ast}^{{\mathrm{\scriptscriptstyle T}}},\boldsymbol{\lambda}_{2,\ast}^{{\mathrm{\scriptscriptstyle T}}})$, \|(\widetilde{\lambda}_{0},\widetilde{\boldsymbol{\lambda}}_{1}^{{\mathrm{\scriptscriptstyle T}}})-(\lambda_{0,\ast},\boldsymbol{\lambda}_{1,\ast}^{{\mathrm{\scriptscriptstyle T}}})\|_{2}=O_{p}\big{(}\{s_{1}\log(m)/n\}^{1/2}\big{)} and $\widetilde{\boldsymbol{\lambda}}_{2}=\boldsymbol{0}$, with probability approaching $1$, where $\lambda_{0,\ast}=\mathbb{E}(D)\cdot\sum_{k=1}^{q}d_{k}$, $\boldsymbol{\lambda}_{1,\ast}=(d_{1},\ldots,d_{q})^{{\mathrm{\scriptscriptstyle T}}}$, and $\boldsymbol{\lambda}_{2,\ast}={\mathbf{0}}$;

(ii) if the PS model in (2.2) is misspecified, the high-dimensional EL optimization problem in (3.4)–(3.4) still admits a feasible solution with probability at least $1-\exp\{-c\log(m)\}$ for a constant $c>0$.

In the special case of (2.2), where $d_{k}=1$ and $d_{h}=0$ for $h\neq k$, indicating that the $k$th PS model is correctly specified, the first part of Proposition 1 shows that the proposed EL estimator of $\pi({\mathbf{X}})$ successfully selects the true propensity score model. Given the convergence of $\widehat{\boldsymbol{\gamma}}_{k}$, the estimated EL weights converge in probability to the true propensity scores, i.e., $\widetilde{\pi}_{i}\stackrel{{\scriptstyle p}}{{\to}}\pi_{i}$. The second part of Proposition 1 demonstrates that, even when all working PS models are misspecified, the empirical likelihood method still yields a solution that satisfies both the IBC constraints in (3.4) and the soft covariate balancing constraints in (3.4). Consequently, despite the misspecification of all working PS models, the proposed weights retain covariate balancing properties and enable valid statistical inference.

We present the main theorem of the proposed method under three different settings.

Theorem 1**.**

Under Conditions 1–6, $\omega_{{\mathrm{\scriptstyle ps}}}\asymp\{s_{1}\log(m)/n\}^{1/2}$ and $\omega_{{\mathrm{\scriptstyle or}}}\asymp\{\log(m)/n\}^{1/2}$,

(i) if any of the working PS models $\{\pi_{k}({\mathbf{x}},\boldsymbol{\gamma}_{k})\}_{k=1}^{q}$ or their linear combination in (2.2) is correctly specified, we have

[TABLE]

(ii) if the OR model in (2.3) is correctly specified, we have

[TABLE]

(iii) if both the OR model and at least one of the working PS models or their linear combination are correctly specified, we have

[TABLE]

Theorem 1 shows that the influence functions of the proposed estimator $\widehat{\mu}_{1}$ have the same form under either the correctly specified PS model in (2.2) or the correctly specified OR model in (2.3). The influence functions indicate the multiply robust inference property of the proposed method. When both the PS and OR models are correctly specified, the proposed estimator $\widehat{\mu}_{1}$ is locally efficient and attains the semiparametric efficiency bound (Hahn,, 1998).

Based on Theorem 1, the variance of $\widehat{\mu}_{1}$ can be estimated by

[TABLE]

and \text{CI}_{1,\alpha}=\big{(}\widehat{\mu}_{1}-z_{1-\alpha/2}(\widehat{V}_{1}/n)^{1/2},\widehat{\mu}_{1}+z_{1-\alpha/2}(\widehat{V}_{1}/n)^{1/2}\big{)} is the proposed $(1-\alpha)$-confidence interval of $\mu_{1}$, where $z_{1-\alpha/2}$ the upper $\alpha/2$ quantile of the standard normal distribution. The confidence intervals of $\mu_{0}$ and $\tau$ can be constructed similarly. Proposition 2 shows the asymptotic coverage of the proposed confidence interval $\text{CI}_{1,\alpha}$.

Proposition 2**.**

Suppose the conditions in Theorem 1 hold. Furthermore, assume the outcome regression error $\epsilon_{1}$ satisfies $1/c_{1}\leq{\rm Var}(\epsilon_{1}|{\mathbf{X}})\leq c_{1}$ for $c_{1}>0$. Then, for $0<\alpha<1$, $\lim_{n\to\infty}\mathbb{P}(\mu_{1}\in\text{CI}_{1,\alpha})=1-\alpha$.

In Proposition 2, we require an additional assumption on the conditional variance of the error in the OR model. A similar assumption is made in Ning et al., (2020) and Tan, 2020a . Under this assumption and conditions in Theorem 1, the variance of the influence function of the proposed estimator $\widehat{\mu}_{1}$ is bounded away from [math] and $\infty$. The asymptotic normality of $\widehat{\mu}_{1}$ follows from Theorem 1.

5 ATE under Generalized Linear Models

We extend the proposed multiply robust estimator in (3.8) to the setting where the outcome variable follows a generalized linear model. Let $f(Y(d)\mid{\mathbf{X}})$ denote the conditional density of the potential outcome $Y(d)$ given ${\mathbf{X}}$. We assume the following OR model for $Y(d)$:

[TABLE]

where $h(\cdot,\cdot)$, $a(\cdot)$ and $g(\cdot)$ are known functions, $\phi$ is the dispersion parameter, $\boldsymbol{\beta}(d)$ is the unknown regression coefficient and ${\mathbf{b}}_{{\mathrm{\scriptstyle or}}}({\mathbf{X}})$ denotes a set of basis functions. Let $g^{\prime}(\cdot)$ and $g^{\prime\prime}(\cdot)$ denote the first- and second-order derivatives of $g(\cdot)$, respectively. Under model (5.1), the conditional variance of $Y(d)$ given ${\mathbf{X}}$ is $\mathbb{V}(Y(d)\mid{\mathbf{X}})=g^{\prime\prime}\{\boldsymbol{\beta}(d)^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{b}}_{{\mathrm{\scriptstyle or}}}({\mathbf{X}})\}a(\phi)$. Accordingly, we balance the covariates $g^{\prime\prime}\{\boldsymbol{\beta}(d)^{{\mathrm{\scriptscriptstyle T}}}{\mathbf{b}}_{{\mathrm{\scriptstyle or}}}({\mathbf{X}})\}{\mathbf{b}}_{{\mathrm{\scriptstyle or}}}({\mathbf{X}})$ when estimating $\mu_{d}$. In the following, we introduce the procedure to estimate $\mu_{1}$. The estimator of $\mu_{0}$ can be similarly constructed.

First, obtain an initial estimator of $\boldsymbol{\beta}(1)$ via regularized maximum likelihood estimation

[TABLE]

where $w_{{\mathrm{\scriptstyle or}}}^{(0)}$ is a regularization parameter.

Second, construct the EL weights $\{\widetilde{p}_{i}\}$ by solving the following optimization problem:

[TABLE]

where $\widetilde{{\mathbf{b}}}_{i}^{(g)}=(\widetilde{b}_{i0}^{(g)},\widetilde{b}_{i1}^{(g)},\ldots,\widetilde{b}_{im}^{(g)})^{{\mathrm{\scriptscriptstyle T}}}=g^{\prime\prime}(\widehat{\boldsymbol{\beta}}^{{\mathrm{\scriptscriptstyle T}}}_{ini}{\mathbf{b}}_{{\mathrm{\scriptstyle or}}}({\mathbf{X}}_{i}))(1,\widehat{{\mathbf{b}}}({\mathbf{X}}_{i})^{{\mathrm{\scriptscriptstyle T}}})^{{\mathrm{\scriptscriptstyle T}}}-\overline{{\mathbf{b}}}^{(g)}$, $\widehat{{\mathbf{b}}}({\mathbf{X}})$ is defined in (3.1), and $\overline{{\mathbf{b}}}^{(g)}=n^{-1}\sum_{i=1}^{n}g^{\prime\prime}(\widehat{\boldsymbol{\beta}}^{{\mathrm{\scriptscriptstyle T}}}_{ini}{\mathbf{b}}_{{\mathrm{\scriptstyle or}}}({\mathbf{X}}_{i}))(1,\widehat{{\mathbf{b}}}({\mathbf{X}}_{i})^{{\mathrm{\scriptscriptstyle T}}})^{{\mathrm{\scriptscriptstyle T}}}$. Different from the EL formulation for the linear outcome model in (3.4), the soft calibration under GLMs is applied to the weighted basis functions $g^{\prime\prime}(\widehat{\boldsymbol{\beta}}^{{\mathrm{\scriptscriptstyle T}}}_{ini}{\mathbf{b}}_{{\mathrm{\scriptstyle or}}}({\mathbf{X}}_{i}))(1,\widehat{{\mathbf{b}}}({\mathbf{X}}_{i})^{{\mathrm{\scriptscriptstyle T}}})^{{\mathrm{\scriptscriptstyle T}}}$, where the weights correspond to the conditional variance of $Y(1)$ given ${\mathbf{X}}$ under the working GLM in (5.1). The high-dimensional EL optimization in (5.2) can be solved via a regularized unconstrained optimization similar to (3.6).

Third, fit the AOR by minimizing the regularized objective function

[TABLE]

where $\widetilde{\pi}_{i}=(n\widetilde{p}_{i})^{-1}$ is obtained from (5.2) and $w_{{\mathrm{\scriptstyle or}}}$ is the penalty parameter. Finally, the estimator $\widehat{\mu}_{1}$ is constructed as

[TABLE]

Compared to the procedure in Section 3 under the linear OR model, the proposed method under GLMs requires only an initial estimate of $\boldsymbol{\beta}(d)$ under the working OR model specified in (5.1). Following the same procedure, we can construct the estimate $\widehat{\mu}_{0}$ of $\mu_{0}$, and define the estimated ATE as $\widehat{\tau}=\widehat{\mu}_{1}-\widehat{\mu}_{0}$.

The $(1-\alpha)$-confidence interval of $\mu_{1}$ is constructed as \text{CI}_{1,\alpha}=\big{(}\widehat{\mu}_{1}-z_{1-\alpha/2}(\widehat{V}_{1}/n)^{1/2},\widehat{\mu}_{1}+z_{1-\alpha/2}(\widehat{V}_{1}/n)^{1/2}\big{)}, where

[TABLE]

is the estimated variance of $\widehat{\mu}_{1}$. Similar to Theorem 1, the theoretical results of the proposed estimator and the confidence interval under the GLM setting are provided in the SM.

6 ATE under Clustered Data

In this section, we assume that the population distribution $\mathcal{G}$ is from a mixture of $v$ subgroups. Let ${\mathbf{U}}_{i}=(U_{i1},\ldots,U_{iv})^{{\mathrm{\scriptscriptstyle T}}}$ be the unobserved subgroup indicator for the $i$th observation, where $U_{i\ell}=1$ and $U_{ih}=0$ for $h\neq\ell$ if the $i$th observation belongs to the $\ell$th cluster. Let $\boldsymbol{\mu}_{x,\ell}$ be the mean of ${\mathbf{X}}_{i}$ for the $\ell$th cluster, and ${\mathbf{M}}=(\boldsymbol{\mu}_{x,1},\ldots,\boldsymbol{\mu}_{x,v})^{{\mathrm{\scriptscriptstyle T}}}$. We assume the following latent structure:

[TABLE]

where ${\mathbf{U}}=({\mathbf{U}}_{1},\ldots,{\mathbf{U}}_{n})^{{\mathrm{\scriptscriptstyle T}}}$, $\boldsymbol{\varepsilon}=(\boldsymbol{\varepsilon}_{1},\ldots,\boldsymbol{\varepsilon}_{n})^{{\mathrm{\scriptscriptstyle T}}}$, and ${\rm Cov}(\boldsymbol{\varepsilon}_{i})=\boldsymbol{\Sigma}_{x}$. Let $\pi({\mathbf{X}}_{i},\ell)=\mathbb{P}(D_{i}=1\mid{\mathbf{X}}_{i},U_{i\ell}=1)$ be the propensity score of the $i$th observation, given it is from the $\ell$th subgroup. If ${\mathbf{U}}_{i}\not\perp(D_{i},{\mathbf{X}}_{i})$, then the subgroups exhibit heterogeneity in both the distribution of covariates ${\mathbf{X}}$ and the treatment assignment mechanism, resulting in differing propensity scores across clusters.

In the following, we discuss the estimation of the propensity score and show it satisfies Condition 5 under clustered data. Under the model in (6.1), Jin et al., (2017); Abbe et al., (2022) showed that spectral clustering can consistently recover the true class labels under a signal-to-noise ratio condition, where signal and noise are determined by the Euclidean distance among the cluster means $\{\boldsymbol{\mu}_{x,1},\ldots,\boldsymbol{\mu}_{x,v}\}$ and the covariance matrix $\boldsymbol{\Sigma}_{x}$, respectively. As discussed in Section 2, multiple PS models can be fitted based on the estimated cluster labels, followed by the application of the proposed multiply robust inference procedure for estimating the ATE. Note that different clustering methods and multiple candidate values for the number of clusters could also be implemented.

Given a set of candidate PS models $\{\tilde{\pi}_{h}({\mathbf{X}},\boldsymbol{\gamma}_{h}):h=1,\ldots,q_{c}\}$ for the subgroups and the estimated cluster labels $\widehat{{\mathbf{U}}}_{i}=(\widehat{U}_{i1},\ldots,\widehat{U}_{iv})^{{\mathrm{\scriptscriptstyle T}}}$ based on ${\mathbf{X}}_{1},\ldots,{\mathbf{X}}_{n}$, we proceed as follows. For notational simplicity, we consider a single set of estimated clusters, though alternative clustering methods can also be implemented. Let $\widehat{\boldsymbol{\gamma}}_{h,\ell}$ be the estimate of the parameter $\boldsymbol{\gamma}_{h}$ in the $h$th PS model, based on the observations from the $\ell$th estimated cluster with $\widehat{U}_{i\ell}=1$. The estimated working PS model for the whole population is defined as

[TABLE]

for $i=1,\ldots,n$, and $h_{1},\ldots,h_{v}=1,\ldots,q_{c}$, where $\widehat{\boldsymbol{\gamma}}_{h_{1}\ldots h_{v}}^{{\mathrm{\scriptscriptstyle T}}}=(\widehat{\boldsymbol{\gamma}}_{h_{1},1}^{{\mathrm{\scriptscriptstyle T}}},\ldots,\widehat{\boldsymbol{\gamma}}_{h_{v},v}^{{\mathrm{\scriptscriptstyle T}}})$ and $1(\cdot)$ denotes the indicator function. The total number of working PS models for the whole population is $q_{c}^{v}$. To study asymptotic properties of $\pi_{h_{1}\ldots h_{v}}({\mathbf{X}}_{i},\widehat{\boldsymbol{\gamma}}_{h_{1}\ldots h_{v}})$ in (6.2) and verify Condition 5, we make the following conditions on the clustering method.

Condition 7**.**

There exists a positive constant $c_{\min}$ such that $\mathbb{P}_{\mathcal{G}}(U_{i\ell}=1)>c_{\min}$ for $\ell=1,\ldots,v$, where $\mathbb{P}_{\mathcal{G}}(\cdot)$ denotes the probability under the mixture distribution $\mathcal{G}$. Furthermore, the clustering error rate satisfies $n^{-1}\sum_{i=1}^{n}1({\mathbf{U}}_{i}\neq\widehat{{\mathbf{U}}}_{i})=O\{s_{1}\log(p)/n\}$.

Condition 7 assumes that each cluster has a non-trivial probability under the population distribution $\mathcal{G}$. Since separate PS models are fitted within each cluster, the condition ensures a sufficient number of observations for each model. Additionally, Condition 7 requires that the estimated clusters converge at a specific rate, $O_{p}\{s_{1}\log(p)/n\}$, where $s_{1}$ is the sparsity level of the parameter in the PS models. This condition can be satisfied by spectral clustering if there is sufficient separation among the cluster means $\{\boldsymbol{\mu}_{x,1},\ldots,\boldsymbol{\mu}_{x,v}\}$ (Jin et al.,, 2017; Abbe et al.,, 2022). Particularly, Cai et al., (2019) showed that $O_{p}\{s_{3}\log(p)/n\}$ is the minimax convergence rate for two-group clustering with sparse signals in the means, where $s_{3}$ represents the sparsity level of the clustering signals.

The following theorem establishes the convergence of estimated PS models in (6.2) for clustered data under logistic regression within each subgroup. In particular, conditions of Theorem 1, which are required for multiply robust inference, are satisfied.

Theorem 2**.**

Suppose that Condition 7 holds, and ${\mathbf{X}}$ follows a sub-Gaussian distribution. For $\ell=1,\dots,v$, consider the logistic regression model $\tilde{\pi}_{h_{\ell}}({\mathbf{X}}_{i},\boldsymbol{\gamma}_{h_{\ell},\ell})=\{1+\exp(-{\mathbf{X}}^{{\mathrm{\scriptscriptstyle T}}}\boldsymbol{\gamma}_{h_{\ell},\ell})\}^{-1}$, and the estimator $\widehat{\boldsymbol{\gamma}}_{h_{\ell},\ell}$ is obtained via $\ell_{1}$-regularized logistic regression with a penalty of order $\{\log(p)/n\}^{1/2}$, using the data $\{({\mathbf{X}}_{i},D_{i}):\widehat{U}_{il}=1\}$ from the estimated subgroup. If the true parameter $\boldsymbol{\gamma}^{*}_{h_{\ell},\ell}$, defined in (S5.11) of the SM, satisfies $\|\boldsymbol{\gamma}^{*}_{h_{\ell},\ell}\|_{0}\leq s_{1}$ and $c_{0}\leq\tilde{\pi}_{h_{\ell}}({\mathbf{X}},\boldsymbol{\gamma}_{h_{\ell},\ell}^{*})\leq 1-c_{0}$ for a positive constant $c_{0}\in(0,1/2)$, then the PS model defined in (6.2) satisfies the conditions required by Theorem 1.

By Theorem 2, the conclusions of Theorem 1 apply to the proposed procedure using the estimated PS model $\{\pi_{h_{1}\ldots h_{v}}({\mathbf{X}}_{i},\widehat{\boldsymbol{\gamma}}_{h_{1}\ldots h_{v}})\}$ in (6.2), obtained from a clustering algorithm for high-dimensional clustered data.

7 Simulation

In this section, we evaluate the performance of the proposed multiply robust inference procedure—debiased high-dimensional empirical likelihood weighting (DELW)—through simulation studies, and compare it with five existing doubly robust inference methods. The results under GLM settings are reported in the SM.

7.1 Multiple Robustness under non-Clustered Data

We first evaluate the proposed method DELW using multiple PS models under non-clustered data. The covariates ${\mathbf{X}}$ are generated from the normal distribution with mean $\boldsymbol{0}$ and covariance $\boldsymbol{\Sigma}_{x}=(\sigma_{x,j_{1}j_{2}})$, where $\sigma_{x,j_{1}j_{2}}=0.5^{|j_{1}-j_{2}|}$. The treatment indicator $D$ is generated from the logistic regression under two designs of propensity scores,

[TABLE]

where $(Z_{1},\ldots,Z_{p})^{{\mathrm{\scriptscriptstyle T}}}={\mathbf{Z}}({\mathbf{X}})=(\exp(X_{1}/2),X_{2}/\{1+\exp(X_{1})\}+10,(X_{1}X_{3}/25+0.6)^{3},(X_{2}+X_{4}+20)^{2},X_{6},\exp(X_{6}+X_{7}),X_{9}^{2},X_{7}^{3}-20,X_{9},\dots,X_{p})^{{\mathrm{\scriptscriptstyle T}}}$ is a nonlinear transformations of ${\mathbf{X}}$ with the first 8 covariates in ${\mathbf{X}}$ being replaced. Two designs of data generation, OR1 and OR2, are applied to the potential outcomes, where

[TABLE]

and $\epsilon_{0},\epsilon_{1}$ are generated from $N(0,1)$, independent of ${\mathbf{X}}$. The observed response is obtained via $Y=DY(1)+(1-D)Y(0)$. These designs provide four data generation mechanisms for $({\mathbf{X}},Y,D)$. Under each setting, a sample of size $n$ is generated IID. The parameter of interest is the ATE $\tau=\mathbb{E}\{Y(1)-Y(0)\}$. We set $n=500,1000$ and $p=1000,2000$. Each setting of the simulation studies is repeated $500$ times.

We compare the proposed method (DELW) with the high-dimensional covariate balancing propensity score (CBPS) method in Ning et al., (2020), model-assisted regularized calibration (RC) method in Tan, 2020a , approximate residual balancing (RB) method in Athey et al., (2018), the double selection (DS) method in Belloni et al., (2014), and the regularized AIPW (rAIPW) method, which naively incorporates regularized estimates of PS and OR models into the AIPW formulation. All five existing methods utilize a single linear OR model and a single logistical PS model with covariates ${\mathbf{X}}$. The four designs of data generation result in four combinations of model specifications for the doubly robust methods: PS1 and OR1 (both models correct), PS1 and OR2 (only PS model correct), PS2 and OR1 (only OR model correct), PS2 and OR2 (both models incorrect). The proposed method DELW uses the linear regression model with covariate ${\mathbf{X}}$ for OR and two working logistical regression models for PS, one with covariates ${\mathbf{X}}$ and the other with covariates ${\mathbf{Z}}$. For the proposed method, as one of the working PS models is correctly specified under each of the four designs, the linear combination condition for PS in (2.2) is satisfied. However, the outcome model is still misspecified under the design OR2.

Both OR and PS models in all the methods are estimated using $\ell_{1}$-regularized regression. We use the 5-fold cross-validation to determine the penalty parameters in DELW and rAIPW. We use R packages CBPS, balanceHD, and hdm to implement CBPS and RC, RB, and DS, respectively. We use the default values of the tuning parameters in the R packages for those four methods.

The results under $n=500$ are presented in Tables 1 and 2 for the PS1 and PS2 designs, respectively. The results under $n=1000$ are reported in the SM. From the two tables, we observe that the proposed method DELW has a small bias and an accurate coverage rate close to the nominal level of 95% for the confidence interval under all scenarios. We compare its performance with the existing methods under the four data generation mechanisms. Under the design of PS1 and OR1, all methods demonstrate reasonable coverage rates and relatively small bias. Both PS and OR models are correctly specified for all methods. Meanwhile, the RMSE of the proposed method, which allows for an additional estimated PS model, is comparable to that of rAIPW and smaller than the RMSEs of all other methods. Similarly, their performance is comparable under the design of PS2 and OR1, as the outcome model is correctly specified for all methods. However, the proposed method still has a relatively smaller RMSE than most existing methods. Under the design of PS1 and OR2, all methods have a correctly specified PS model but a misspecified outcome model. The proposed method is much more efficient than all other methods, with the smallest standard error and shortest confidence interval. The efficiency gain is attributed to adding additional variables in the soft covariate balancing and augmented outcome regression. Under the design of PS2 and OR2, both the PS and OR models are misspecified for the existing methods. Those methods display significantly larger bias and higher standard error and RMSE, resulting in much wider confidence intervals and coverage rates below the nominal level. In contrast, the proposed method DELW, by incorporating an additional propensity model, exhibits the smallest RMSE and achieves robustness inference with coverage close to the nominal level. In summary, the proposed estimator shows a smaller RMSE, enhanced robustness to model misspecification, and an accurate coverage rate of the confidence interval. Additional simulation studies under GLM for the response variable are provided in the SM.

7.2 Clustered Data

Then, we evaluate the performance of the proposed DELW method under clustered data by simulation. We consider the data generation mechanisms, where the covariates $\{{\mathbf{X}}_{i}\}_{i=1}^{n}$ are independently drawn from two distinct clusters. The first half of samples ($i=1,\dots,n/2$) are from cluster 1, where ${\mathbf{X}}_{i}$ is generated from the multivariate normal distribution with mean $-\boldsymbol{\mu}_{x}$ and covariance $\boldsymbol{\Sigma}_{x}=(\sigma_{x,j_{1}j_{2}})$ with $\sigma_{x,j_{1}j_{2}}=0.5^{|j_{1}-j_{2}|}$. The second half of the samples ($i=n/2+1,\dots,n$) are from cluster 2, where ${\mathbf{X}}_{i}$ is generated from the multivariate normal distribution with mean $\boldsymbol{\mu}_{x}$ and the same covariance $\boldsymbol{\Sigma}_{x}$. The discriminant vector $\boldsymbol{\mu}_{x}=(\mu_{x,1},\ldots,\mu_{x,p})^{{\mathrm{\scriptscriptstyle T}}}$ is sparse with only the second set of $10$ entries being nonzero, where $\mu_{x,j}=1$ for $j=11,\ldots,20$ and $\mu_{x,j}=0$ for all other $j$. This setup is similar to that considered by Cai et al., (2019) for high-dimensional clustering.

Two different designs of the propensity scores are considered for the two clusters, where $\mathbb{P}(D=1\mid{\mathbf{X}})$ follows PS1 for cluster 1 and PS1 with ${\mathbf{X}}$ replaced by $-{\mathbf{X}}$ for cluster 2, or equivalently, $\mathbb{P}(D=0\mid{\mathbf{X}})$ follows PS1 for cluster 2. For the proposed method DELW, we apply the high-dimensional clustering algorithm based on the covariate ${\mathbf{X}}$ by the R package HDclassif with two clusters and fit $\ell_{1}$ regularized logistic regression on the two estimated clusters separately. The clustering algorithm employs the expectation-maximization (EM) algorithm to estimate the model parameters, which are sensitive to the initial values of the parameters. We use two sets of initial values to generate two propensity score models, both incorporated into the proposed procedure. All OR and PS models are estimated by $\ell_{1}$ regularized regression. Each setting of the simulation studies is repeated $500$ times.

The results are presented in Table 3. From this table, we observe that the proposed method DELW has a small bias and an accurate coverage rate close to the nominal level of 95% for the confidence interval under both designs for the potential outcomes. Under the OR1 design, where the outcome regression model is correctly specified for all methods, the results of those methods are comparable, achieving reasonable coverage with a small bias. However, the proposed method DELW has the smallest RMSE and the shortest confidence interval. Under the OR2 design, where the outcome regression model is misspecified for all methods, the five existing methods display significant biases, large standard errors, and inaccurate coverage rates below the nominal level. In contrast, DELW also shows a much smaller RMSE and a shorter confidence interval with a coverage rate closer to the nominal level, indicating the validity of the proposed method under clustered data.

8 Real Data Analysis

This section presents a real-world data example from heterogeneous sources to demonstrate the need for developing a multiply robust method for estimating the ATE. The right heart catheterization dataset was initially reported in Connors et al., (1996) to study the treatment effect of RHC for seriously ill patients in the intensive care unit (ICU) on their subsequent one-month survival rate. The data set contains $5735$ patients with $72$ features from two studies across five medical centers. In the original work, Connors et al., (1996) combined the data sets from those two potentially heterogeneous studies and used a single logistic regression to fit the propensity scores. The following works also took the similar strategy, for example, Hirano and Imbens, (2001); Vermeulen and Vansteelandt, (2015); Tan, 2020a . However, the mechanism of treatment selection could be different between the two studies and the five medical centers, and the subgroup indicators were not released in the dataset.

For each patient, the treatment indicator $D$ takes the value $1$ if the patient received the RHC treatment within 24 hours of admission to the ICU, and $D=0$ otherwise. The outcome variable $Y$ takes the value $1$ if the patient survived up to $30$ days, and $Y=0$ otherwise. The parameter of interest is the ATE of applying RHC on the survival probability. We apply two clustering algorithms by R packages HDclassif and ClusterR with two clusters based on the covariates. We estimate PS by fitting logistic regression models of $D$ on the estimated clusters separately. The results of the first clustering algorithm are presented in Figure 1, which shows that the estimated clusters are separated along the direction of the first principal component. The right panel of Figure 1 shows that the estimated regression coefficients of the PS models differ a lot between the two clusters, indicating the heterogeneity of the propensity scores. The estimated clusters by the second algorithm are presented in the SM.

We use a logistic regression model to fit the binary outcome $Y$ and construct the covariate vector ${\mathbf{X}}$ with a dimension $p=498$, which includes all main effects and two-way interactions of the $72$ original covariates, excluding those with a fraction of nonzero values lower than $10\%$ of the sample size. Table 4 reports the estimates of the one-month survival probabilities in the treatment and control group, along with the estimated ATE obtained by the proposed DELW, RC, and rAIPW methods, where the results for RC and rAIPW are taken from Table 2 in Tan, 2020a . The estimated standard errors of the three methods are similar. According to the results from RC and rAIPW, the ATE is significantly negative, suggesting that taking the RHC treatment would significantly lower the one-month survival probability. The conclusion raises practical concerns. However, the estimated ATE by the proposed method is smaller in absolute value and not significantly different from zero at the $0.05$ significance level.

From a medical perspective, the effectiveness of the RHC test is controversial. Stevenson, (2005) concluded that adding RHC to careful clinical assessment increased anticipated adverse events but did not affect overall mortality and hospitalization. Garan et al., (2020) showed that using RHC is associated with improved survival from cardiogenic shock. Currently, RHC is still recommended in clinical practice (Fleisher et al.,, 2014). To evaluate certain critical cardiopulmonary conditions, RHC is performed when it is believed to be crucial for adjusting treatment plans and when it is unlikely to exacerbate the existing condition. Our analyses of the RHC dataset indicate that RHC treatment had an insignificant impact on early hospital mortality. The finding reflects a positive view on the use of RHC, aligning with its current widespread application.

9 Discussion

This paper proposes a new method for multiply robust inference of ATEs under high-dimensional covariates. The method is particularly useful for clustered data with heterogeneous treatment assignment mechanisms, where correctly specifying a single parametric PS model could be difficult. The proposed procedure complements the high-dimensional doubly robust procedures by utilizing multiple PS models. Although machine learning and nonparametric estimation of propensity scores and outcome regression can be applied to estimate ATEs, these methods generally do not incorporate multiple estimated models, and their properties under general high-dimensional settings remain unclear. The proposed method provides valid statistical inference of the ATE, if any working PS models, their linear combination, or the OR model is correctly specified. The confidence interval can be computed efficiently without requiring bootstrap procedures. In the presence of unobserved clusters, the ATE can be identified under Condition 1 of ignorability, which is an easier case than causal inference under unobserved confounders. An important future direction is to extend multiply robust inference to settings with unobserved confounders where the ignorability assumption doesn’t hold.

Appendix A Multiple Robustness by Augmented Calibration and Outcome Regression

In the following, we outline the key steps to explain the importance of soft calibration using augmented covariates with PS derivatives, and augmented outcome regression to achieve multiply robust inferences under high-dimensional covariates. Recall that the subscript $\ast$ denotes the probability limit of an estimate under either a correctly specified or misspecified model. Calibrating the PS derivatives $\{\pi_{k}^{(1)}({\mathbf{X}},\widehat{\boldsymbol{\gamma}}_{k})\}_{k=1}^{q}$ leads to the approximation:

[TABLE]

where the regression coefficients from the AOR in (3.7) are replaced by their probability limits.

Suppose the OR model in (2.3) is correctly specified, then $\boldsymbol{\beta}_{1k,\ast}={\mathbf{0}}$ for $k=1,\ldots,q$, $\boldsymbol{\beta}_{2,\ast}={\mathbf{0}}$. In this case, $\widehat{\mu}_{1}$ can be asymptotically expressed as

[TABLE]

where $\pi_{i}^{*}$ is the probability limit of $\widetilde{\pi}_{i}$.

Alternatively, suppose one of the working PS models in (2.2), say $\pi_{1}({\mathbf{X}}_{i},\boldsymbol{\gamma}_{1})$, is correctly specified, but the OR model is subject to misspecification. In this case, we can replace $\widehat{{\mathbf{b}}}({\mathbf{X}}_{i})$ and $\widehat{\boldsymbol{\pi}}_{i}$ in (A.1) by their probability limits since ${\mathbb{E}}\{(D/\pi({\mathbf{X}})-1)u({\mathbf{X}})\}=0$ for any integrable function $u({\mathbf{X}})$. This gives:

[TABLE]

where $\boldsymbol{\pi}_{i,\ast}^{{\mathrm{\scriptscriptstyle T}}}=(\pi_{1}({\mathbf{X}}_{i},\boldsymbol{\gamma}_{1,*}),\dots,\pi_{q}({\mathbf{X}}_{i},\boldsymbol{\gamma}_{q,*}))$ contains the limits of the PS models and ${\mathbf{b}}_{\ast}({\mathbf{X}}_{i})$ is the augmented covariate functions using probability limits of PS models for $i=1,\dots,n$. Using Taylor’s expansion, the first term in (A) can be decomposed as

[TABLE]

where $\widetilde{\boldsymbol{\lambda}}=(\widetilde{\lambda}_{0},\widetilde{\boldsymbol{\lambda}}_{1}^{{\mathrm{\scriptscriptstyle T}}},\widetilde{\boldsymbol{\lambda}}_{2}^{{\mathrm{\scriptscriptstyle T}}})^{{\mathrm{\scriptscriptstyle T}}}$, $\boldsymbol{\lambda}_{\ast}=(\lambda_{0,\ast},\boldsymbol{\lambda}_{1,\ast}^{{\mathrm{\scriptscriptstyle T}}},\boldsymbol{\lambda}_{2,\ast}^{{\mathrm{\scriptscriptstyle T}}})^{{\mathrm{\scriptscriptstyle T}}}$, and $\widetilde{{\mathbf{a}}}_{i}=(1,(\widehat{\boldsymbol{\pi}}_{i}-n_{1}/n\mathbf{1}_{q})^{{\mathrm{\scriptscriptstyle T}}},\widetilde{{\mathbf{b}}}_{i}^{{\mathrm{\scriptscriptstyle T}}})^{{\mathrm{\scriptscriptstyle T}}}$. The two terms in (A.4) and (A.5) are negligible due to the construction of the AOR with weights $\{\widetilde{\pi}_{i}^{-2}\}_{i=1}^{n}$ in (3.7). Notably, the influence functions of $\widehat{\mu}_{1}$ are the same in (A.2) and the second term of (A), which establishes the validity of multiply robust confidence intervals for $\mu_{1}$.

From (A.4) and (A.5), the estimation of the nuisance parameters $\{\boldsymbol{\gamma}_{k,\ast}\}_{k=1}^{q}$, $\lambda_{0}$, $\boldsymbol{\lambda}_{1}$, and $\boldsymbol{\lambda}_{2}$ does not impact the inference of the target parameter $\mu_{1}$, nor does the estimation of the regression coefficients $\boldsymbol{\beta}_{1}$ and $\boldsymbol{\beta}_{2}$ from (A.1). The proposed soft calibration and weighted AOR are constructed so that all high-dimensional nuisance parameters are Neyman orthogonal to the ATE. This is achieved by carefully designing the regularized estimation in (3.6) and (3.7) based on KKT conditions.