Estimating Average Treatment Effects with a Double-Index Propensity Score

TL;DR

This paper introduces the Double-index Propensity Score (DiPS), a novel method for estimating average treatment effects in observational studies that combines covariate dimension reduction, smoothing, and regularization for improved efficiency and robustness.

Contribution

The paper proposes the DiPS estimator, integrating adaptive LASSO-based dimension reduction with a smoothing approach to enhance ATE estimation under model misspecification and high-dimensional covariates.

Findings

DiPS maintains double-robustness and local semiparametric efficiency.

Smoothing improves efficiency and robustness under model misspecification.

Simulation studies demonstrate finite sample benefits.

Abstract

We consider estimating average treatment effects (ATE) of a binary treatment in observational data when data-driven variable selection is needed to select relevant covariates from a moderately large number of available covariates $\mathbf{X}$. To leverage covariates among $\mathbf{X}$ predictive of the outcome for efficiency gain while using regularization to fit a parametric propensity score (PS) model, we consider a dimension reduction of $\mathbf{X}$ based on fitting both working PS and outcome models using adaptive LASSO. A novel PS estimator, the Double-index Propensity Score (DiPS), is proposed, in which the treatment status is smoothed over the linear predictors for $\mathbf{X}$ from both the initial working models. The ATE is estimated by using the DiPS in a normalized inverse probability weighting (IPW) estimator, which is found to maintain double-robustness and also local semiparametric efficiency with a fixed number of covariates $p$. Under misspecification of working models, the smoothing step leads to gains in efficiency and robustness over traditional doubly-robust estimators. These results are extended to the case where $p$ diverges with sample size and working models are sparse. Simulations show the benefits of the approach in finite samples. We illustrate the method by estimating the ATE of statins on colorectal cancer risk in an electronic medical record (EMR) study and the effect of smoking on C-reactive protein (CRP) in the Framingham Offspring Study.

Tables

Table 1. Table 1: Bias and RMSE of estimators by n 𝑛 n and model specification scenario for p = 15 𝑝 15 p=15 .

		Both Correct		Misspecified ${\mu }_k\left(𝐱\right)$		Misspecified ${\pi }_1\left(𝐱\right)$
Size	Estimator	Bias	RMSE	Bias	RMSE	Bias	RMSE
n=500	IPW-ALAS	0.029	0.350	0.074	1.754	0.023	0.294
	DR-ALAS	0.002	0.330	0.029	1.684	-0.001	0.285
	DR-SIM	-0.021	0.315	0.127	1.495	0.013	0.287
	OAL	0.008	0.321	0.074	1.484	0.001	0.284
	GLiDeR	0.001	0.299	0.087	1.238	0.006	0.282
	MADR	0.022	0.300	0.172	1.247	0.008	0.282
	DiPS	-0.017	0.319	0.101	1.193	0.013	0.293
n=5,000	IPW-ALAS	0.001	0.111	-0.002	0.588	0.033	0.108
	DR-ALAS	-0.003	0.106	-0.014	0.564	-0.008	0.089
	DR-SIM	-0.012	0.103	0.029	0.516	-0.004	0.089
	OAL	-0.002	0.105	0.000	0.527	-0.007	0.089
	GLiDeR	-0.001	0.098	0.034	0.413	-0.006	0.088
	MADR	0.000	0.099	0.124	0.418	-0.008	0.089
	DiPS	-0.016	0.106	0.041	0.349	-0.003	0.091

Table 2. Table 2: Perturbation performance under correctly specified models. Emp SE: empirical standard error over simulations, ASE: average of standard error estimates based on MAD over perturbations, Cover: Coverage of 95% percentile intervals.

p	n	Emp SE	ASE	Cover
15	500	0.350	0.362	0.966
15	2500	0.151	0.167	0.970
15	5000	0.108	0.119	0.965
30	500	0.348	0.356	0.961
30	2500	0.150	0.167	0.975
30	5000	0.103	0.119	0.973

Table 3. Table 3: Data example on the effect of statins on CRC risk in EMR data and the effect of smoking on logCRP in FOS data. Est: Point estimate, SE: estimated SE, 95% CI: confidence interval, p-val: p-value from Wald test of no effect.

IBD EMR Study					FOS
	Est	SE	95% CI	p-val	Est	SE	95% CI	p-val
None	-0.008	0.004	(-0.017, 0)	0.047	0.180	0.058	(0.065, 0.298)	0.002
IPW-ALAS	-0.022	0.004	(-0.031, -0.015)	0.001	0.182	0.063	(0.053, 0.307)	0.004
DR-ALAS	-0.020	0.005	(-0.029, -0.012)	0.001	0.140	0.063	(0.031, 0.277)	0.026
DR-SIM	-0.023	0.003	(-0.029, -0.018)	0.001	0.143	0.057	(0.044, 0.257)	0.013
OAL	-0.008	0.004	(-0.017, 0)	0.048	0.175	0.061	(0.062, 0.301)	0.004
GLiDeR	-0.031	0.005	(-0.04, -0.022)	0.001	0.147	0.058	(0.045, 0.258)	0.012
MADR	-0.030	0.005	(-0.04, -0.021)	0.001	0.149	0.056	(0.037, 0.258)	0.008
DiPS	-0.024	0.003	(-0.029, -0.017)	0.001	0.141	0.058	(0.039, 0.276)	0.015

Equations

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

$$Y=T{Y}^{(1)}+(1-T)Y^{(0)}\text{ with probability }1$$

$${\pi }_1(\mathbf{x})\in [{ϵ}_{\pi },1-{ϵ}_{\pi }]\text{ for some }{ϵ}_{\pi }>0,\text{ when }\mathbf{x}\in \mathcal{X}$$

$$Y^{(1)}\perp \perp T\mid \mathbf{X}\text{ and }{Y}^{(0)}\perp \perp T\mid \mathbf{X},$$

$${\Delta }^{\ast }=\mathbb{E}\{{\mu }_1(\mathbf{X})-{\mu }_0(\mathbf{X})\}=\mathbb{E}\left\{\frac{I(T=1)Y}{{\pi }_1(\mathbf{X})}-\frac{I(T=0)Y}{{\pi }_0(\mathbf{X})}\right\},$$

$${\pi }_1(\mathbf{x})=g_{\pi }({\alpha }_0+{\mathbf{\alpha }}^{\mathsf{T}}\mathbf{x}),$$

$${\mu }_k(\mathbf{x})=g_{\mu }({\beta }_0+{\beta }_{1}k+{\mathbf{\beta }}_k^{\mathsf{T}}\mathbf{x}),\text{ for }k=0,1,$$

$$({\hat{\alpha }}_0,{\hat{\mathbf{\alpha }}}^{\mathsf{T}}{)}^{\mathsf{T}}=\underset{\vec{\mathbf{\alpha }}}{\mathrm{argmax}}\left\{n^{-1}\sum _{i=1}^n{\ell }_{\pi }(\vec{\mathbf{\alpha }};T_i,{\mathbf{X}}_i)-{\lambda }_{\pi ,n}\sum _{j=1}^p\left|{\alpha }_j\right|/{\left|{\tilde{\alpha }}_j\right|}^{\gamma }\right\}$$

$$({\hat{\beta }}_0,{\hat{\beta }}_1,{\hat{\mathbf{\beta }}}_0^{\mathsf{T}},{\hat{\mathbf{\beta }}}_1^{\mathsf{T}}{)}^{\mathsf{T}}=\underset{\vec{\mathbf{\beta }}}{\mathrm{argmax}}\left\{n^{-1}\sum _{i=1}^n{\ell }_{\mu }(\vec{\mathbf{\beta }};{\mathbf{Z}}_i)-{\lambda }_{\mu ,n}\left({\left|\beta \right|}_1/{\left|{\tilde{\beta }}_1\right|}^{\gamma }+\sum _{k=0}^1\sum _{j=2}^p\left|{\beta }_{k,j}\right|/{\left|{\tilde{\beta }}_{k,j}\right|}^{\gamma }\right)\right\},$$

$$\begin{array}{rl}& ({\bar{\alpha }}_0,{\bar{\mathbf{\alpha }}}^{\mathsf{T}}{)}^{\mathsf{T}}\text{ uniquely maximize }\mathbb{E}\left\{{\ell }_{\pi }(\vec{\mathbf{\alpha }};T_i,{\mathbf{X}}_i)\right\}\\ & ({\bar{\beta }}_0,{\bar{\beta }}_1,{\bar{\mathbf{\beta }}}_0^{\mathsf{T}},{\bar{\mathbf{\beta }}}_1^{\mathsf{T}}{)}^{\mathsf{T}}\text{ uniquely maximize }\mathbb{E}\left\{{\ell }_{\mu }(\vec{\mathbf{\beta }};{\mathbf{Z}}_i)\right\}.\end{array}$$

$$s_{\mathbf{\alpha }},s_{{\mathbf{\beta }}_0}\text{ and }{s}_{{\mathbf{\beta }}_1}\text{ are fixed as }n\to \infty .$$

$${\hat{\pi }}_k(\mathbf{x};{\hat{\mathbf{\theta }}}_k)=\frac{n^{-1}{\sum }_{j=1}^n{K}_h\{(\hat{\mathbf{\alpha }},{\hat{\mathbf{\beta }}}_k{)}^{\mathsf{T}}({\mathbf{X}}_j-\mathbf{x})\}I(T_j=k)}{n^{-1}{\sum }_{j=1}^n{K}_h\{(\hat{\mathbf{\alpha }},{\hat{\mathbf{\beta }}}_k{)}^{\mathsf{T}}({\mathbf{X}}_j-\mathbf{x})\}},\text{ for }k=0,1,$$

$${\hat{\mu }}_k={\left\{\sum _{i=1}^n\frac{I(T_i=k)}{{\hat{\pi }}_k({\mathbf{X}}_i;{\hat{\mathbf{\theta }}}_k)}\right\}}^{-1}{\left\{\sum _{i=1}^n\frac{I(T_i=k)Y_i}{{\hat{\pi }}_k({\mathbf{X}}_i;{\hat{\mathbf{\theta }}}_k)}\right\}}^{-1},\text{ for }k=0,1.$$

$${\bar{\mu }}_k=\mathbb{E}\left\{\frac{I(T_i=k)Y_i}{{\pi }_k({\mathbf{X}}_i;{\bar{\mathbf{\theta }}}_k)}\right\},\text{ for }k=0,1,$$

$${\bar{\mu }}_k=\mathbb{E}\left\{\frac{I(T_i=k)Y_i}{{\pi }_k({\mathbf{X}}_i)}\right\}=\mathbb{E}\left\{Y_i^{(k)}\right\},\text{ for }k=0,1.$$

$${\bar{\mu }}_k$$

$${\hat{\mathcal{W}}}_k$$

$$+n^{-1/2}\sum _{i=1}^n{\mathbf{u}}_{k,{\mathcal{A}}_{\mathbf{\alpha }}}^{\mathsf{T}}{\mathbf{\Psi }}_{i,{\mathcal{A}}_{\mathbf{\alpha }}}+{\mathbf{v}}_{k,{\mathcal{A}}_{{\mathbf{\beta }}_k}}^{\mathsf{T}}{\mathbf{{\rm Y}}}_{i,k,{\mathcal{A}}_{{\mathbf{\beta }}_k}}+O_p(n^{1/2}{h}^q+n^{-1/2}{h}^{-2}),$$

$${\hat{\mathcal{W}}}_k$$

$$+n^{-1/2}\sum _{i=1}^n\left\{\frac{I(T_i=k)}{{\hat{\pi }}_k({\mathbf{X}}_i;{\hat{\mathbf{\theta }}}_k)}-\frac{I(T_i=k)}{{\pi }_k({\mathbf{X}}_i;{\bar{\mathbf{\theta }}}_k)}\right\}(Y_i-{\bar{\mu }}_k)+o_p(1).$$

$${\phi }_{i,k}=\frac{I(T_i=k)Y_i}{{\pi }_k({\mathbf{X}}_i;{\bar{\mathbf{\theta }}}_k)}-\left\{\frac{I(T_i=k)}{{\pi }_k({\mathbf{X}}_i;{\bar{\mathbf{\theta }}}_k)}-1\right\}\mathbb{E}(Y_i\mid {\bar{\mathbf{\alpha }}}^{\mathsf{T}}{\mathbf{X}}_i,{\bar{\mathbf{\beta }}}_k^{\mathsf{T}}{\mathbf{X}}_i,T_i=k)-{\bar{\mu }}_k.$$

$${\varphi }_{i,k}=\frac{I(T_i=k)Y_i}{{\pi }_k({\mathbf{X}}_i)}-\left\{\frac{I(T_i=k)}{{\pi }_k({\mathbf{X}}_i)}-1\right\}{g}_{\mu }({\bar{\beta }}_0+{\bar{\beta }}_{1}k+{\bar{\mathbf{\beta }}}_k^{\mathsf{T}}{\mathbf{X}}_i)-{\bar{\mu }}_k.$$

$$({\hat{\alpha }}_0^{\ast },{\hat{\mathbf{\alpha }}}^{\ast \mathsf{T}}{)}^{\mathsf{T}}=\underset{\vec{\mathbf{\alpha }}}{\mathrm{argmax}}\left\{n^{-1}\sum _{i=1}^n{\ell }_{\pi }(\vec{\mathbf{\alpha }};T_i,{\mathbf{X}}_i)G_i-{\lambda }_{\pi ,n}\sum _{j=1}^p\left|{\alpha }_j\right|/{\left|{\tilde{\alpha }}_j^{\ast }\right|}^{\gamma }\right\}$$

$$({\hat{\beta }}_0^{\ast },{\hat{\beta }}_1^{\ast },{\hat{\mathbf{\beta }}}_0^{\ast \mathsf{T}},{\hat{\mathbf{\beta }}}_1^{\ast \mathsf{T}}{)}^{\mathsf{T}}=\underset{\vec{\mathbf{\beta }}}{\mathrm{argmax}}\left\{n^{-1}\sum _{i=1}^n{\ell }_{\mu }(\vec{\mathbf{\beta }};{\mathbf{Z}}_i)G_i-{\lambda }_{\mu ,n}\sum _{j=1}^p\left|{\beta }_j\right|/{\left|{\tilde{\beta }}_j^{\ast }\right|}^{\gamma }\right\},$$

$${\hat{\pi }}_k^{\ast }(\mathbf{x};{\hat{\mathbf{\theta }}}_k^{\ast })=\frac{{\sum }_{j=1}^n{K}_h\{({\hat{\mathbf{\alpha }}}^{\ast },{\hat{\mathbf{\beta }}}_k^{\ast }{)}^{\mathsf{T}}({\mathbf{X}}_j-\mathbf{x})\}I(T_j=k)G_j}{{\sum }_{j=1}^n{K}_h\{({\hat{\mathbf{\alpha }}}^{\ast },{\hat{\mathbf{\beta }}}_k^{\ast }{)}^{\mathsf{T}}({\mathbf{X}}_j-\mathbf{x})\}{G}_j},\text{ for }k=0,1.$$

$${\hat{\mu }}_k^{\ast }={\left\{\sum _{i=1}^n\frac{I(T_i=k)}{{\hat{\pi }}_k^{\ast }({\mathbf{X}}_i;{\hat{\mathbf{\theta }}}_k^{\ast })}{G}_i\right\}}^{-1}{\left\{\sum _{i=1}^n\frac{I(T_i=k)Y_i}{{\hat{\pi }}_k^{\ast }({\mathbf{X}}_i;{\hat{\mathbf{\theta }}}_k^{\ast })}{G}_i\right\}}^{-1},\text{ for }k=0,1.$$

$$\text{Both correct: }{\pi }_1(\mathbf{x})=g_{\pi }(.2+{\mathbf{\alpha }}^{\mathsf{T}}\mathbf{x}),{\mu }_k(\mathbf{x})=k+{\mathbf{\beta }}^{\mathsf{T}}\mathbf{x}$$

$$\text{Misspecified μk(x): }{\pi }_1(\mathbf{x})=g_{\pi }(.2+{\mathbf{\alpha }}^{\mathsf{T}}\mathbf{x}),{\mu }_k(\mathbf{x})=k+{\mathbf{\beta }}_{[1]}^{\mathsf{T}}\mathbf{x}(1+{\mathbf{\beta }}_{[2]}^{\mathsf{T}}\mathbf{x})+k{\mathbf{\zeta }}^{\mathsf{T}}\mathbf{x}$$

$$\text{Misspecified πk(x): }{\pi }_1(\mathbf{x})=g_{\pi }\left\{.2+{\mathbf{\alpha }}_{[1]}^{\mathsf{T}}\mathbf{x}(1+{\mathbf{\alpha }}_{[2]}^{\mathsf{T}}\mathbf{x})\right\},{\mu }_k(\mathbf{x})=k+{\mathbf{\beta }}^{\mathsf{T}}\mathbf{x},$$

$${\hat{\pi }}_k(\mathbf{x};\mathbf{\theta })={\hat{\pi }}_k(\mathbf{x};\mathbf{\alpha },\mathbf{\beta })=\frac{{\hat{l}}_k(\mathbf{x};\mathbf{\theta })}{\hat{f}(\mathbf{x};\mathbf{\theta })}=\frac{{\sum }_{j=1}{K}_h\{(\mathbf{\alpha },\mathbf{\beta }{)}^{\mathsf{T}}({\mathbf{X}}_j-\mathbf{x})\}I(T_i=k)}{{\sum }_{j=1}{K}_h\{(\mathbf{\alpha },\mathbf{\beta }{)}^{\mathsf{T}}({\mathbf{X}}_j-\mathbf{x})\}}.$$

$$n^{1/2}\left|{\bar{X}}_n-{\mu }_n\right|\le O_p(c_n),$$

$$\mathbb{P}\left(n^{1/2}\left|{\bar{X}}_n-{\mu }_n\right|/c_n\ge k\right)\le \frac{{\sigma }_n^2}{c_n^2{k}^2}.$$

$$\underset{\mathbf{x}}{\sup }\Vert {\hat{\pi }}_k(\mathbf{x};\bar{\mathbf{\theta }})-{\pi }_k(\mathbf{x};\bar{\mathbf{\theta }})\Vert =O_p(a_n),$$

$$a_n=h^q+\{log(n)/(n{h}^2){\}}^{1/2}.$$

$$n^{-1}\sum _{i=1}^n\frac{\partial }{\partial {\mathbf{\alpha }}^{\mathsf{T}}}\frac{g({\mathbf{Z}}_i)}{{\hat{\pi }}_k({\mathbf{X}}_i;\bar{\mathbf{\theta }})}=\mathbb{E}\left[{\dot{K}}_h({\bar{\mathbf{S}}}_{ji}{)}^{\mathsf{T}}\left\{{\pi }_k({\bar{\mathbf{S}}}_i)-I(T_j=k)\right\}\frac{g({\mathbf{Z}}_i)}{{\pi }_k({\bar{\mathbf{S}}}_i)l_k({\bar{\mathbf{S}}}_i)}{\mathbf{X}}_{ji}^{†\mathsf{T}}\right]+O_p(b_n)$$

Full text

1 Introduction

There is growing interest in evaluating medical treatments and policies in large-scale observational data such as electronic medical records (EMR). As with any observational data, in the absence of randomization, adjustment for a sufficient set of pre-treatment covariates $\mathbf{X}$ that satisfy “no unmeasured confounding” is needed when estimating average treatment effects (ATE) to avoid confounding bias. This is routinely done using propensity score (PS), outcome regression, and doubly-robust (DR) methods (Lunceford and Davidian, 2004). These methods were initially developed in settings where $p$, the dimension of $\mathbf{X}$, was small relative to the sample size $n$. But large-scale observational data are increasingly collecting rich measurements in large sets of covariates, and data-driven variable selection approaches are needed due to the lack of sufficient prior knowledge to guide manual variable selection.

Effective variable selection for causal effect estimation involves consideration of dependencies between $\mathbf{X}$ with the treatment status $T\in\{0,1\}$ and outcome $Y$. Let $\mathcal{A}_{\pi}\subseteq\{1,2,\ldots,p\}$ index the subset of $\mathbf{X}$ upon which the PS $\pi_{1}(\mathbf{x})=\mathbb{P}(T=1\mid\mathbf{X}=\mathbf{x})$ depends, and let $\mathcal{A}_{\mu}$ be an analogous index set for $\mathbf{X}$ upon which either $\mu_{1}(\mathbf{x})$ or $\mu_{0}(\mathbf{x})$ depends, where $\mu_{k}(\mathbf{x})=\mathbb{E}(Y\mid\mathbf{X}=\mathbf{x},T=k)$. For any index set $\mathcal{S}\subseteq\{1,2,\ldots,p\}$, let $\mathcal{S}^{c}$ denote its complement in $\{1,2,\ldots,p\}$. When $\mathbf{X}$ is sufficient for no unmeasured confounding, the covariates indexed in $\mathcal{A}_{\pi}$ is a reduced set of covariates that is also sufficient for no unmeasured confounding (De Luna et al., 2011). However, additionally adjusting for purely prognostic covariates in $\mathcal{A}_{\pi}^{c}\cap\mathcal{A}_{\mu}$ can improve the efficiency of PS, outcome regression, and DR estimators (Lunceford and Davidian, 2004; Hahn, 2004; Brookhart et al., 2006).

To exploit this phenomenon, we consider an inverse probability weighting (IPW) estimator where the PS is initially estimated by regularized regression. Since variable selection procedures for the PS model would select out covariates in $\mathcal{A}_{\pi}^{c}\cap\mathcal{A}_{\mu}$, we also estimate a regularized regression model for $\mu_{k}(\mathbf{x})$, for $k=0,1$, to recover variation from covariates in $\mathcal{A}_{\pi}^{c}\cap\mathcal{A}_{\mu}$ to inform estimation of a calibrated PS. The calibration is implemented through smoothing $T$ over the linear predictors for $\mathbf{X}$ from both the initial PS and outcome models, which can be viewed as smoothing over working propensity and prognostic scores (Hansen, 2008a). The resulting IPW estimator maintains double-robustness and achieves the semiparametric efficiency bound when $p$ is fixed, under correctly specified PS and outcome working models. To the best of our knowledge, this is the first proposal in the literature that demonstrates these properties can be achieved through weighting only, without explicit augmentation. We show that the estimator is asymptotically linear and use this to characterize large-sample robustness and efficiency properties. The smoothing results in a refinement of the influence function under misspecification of the outcome model that can potentially result in substantial gains in efficiency relative to traditional DR estimators, which is confirmed in simulations. These properties hold in settings where $p$ is either fixed or allowed to diverge slowly with $n$ assuming fixed sparsity indices.

Data-driven variable selection for causal effect estimation has been considered in screening methods based on marginal associations between $\mathbf{X}$ with $T$ and $Y$ (Schneeweiss et al., 2009), but the results can be misleading because marginal associations need not agree with conditional associations. De Luna et al. (2011) carefully characterized and proposed algorithms to identify minimal subsets of covariates that are sufficient for no unmeasured confounding. Recent works have considered using regularized regression to select variables and post-selection methods that estimate treatment effects through partially linear models (Belloni et al., 2013) and DR estimators (Farrell, 2015; Belloni et al., 2017). These methods focus on delivering uniformly valid inference under high-dimensional regimes assuming approximately sparse models. Others have proposed modifying the regularization penalty itself in a way to select the relevant covariates and estimate treatment effects through IPW (Shortreed and Ertefaie, 2017) and DR estimators (Koch et al., 2018). However, these papers generally do not fully work out the full asymptotic distribution of the final estimator, making efficiency comparisons with established methods difficult. Some of the methods are also only singly-robust. Bayesian model averaging (Cefalu et al. (2017) and references therein) offers a principled alternative for variable selection but encounters burdensome computations that are possibly infeasible for large $p$.

Our proposed double-index PS (DiPS) can be viewed as a simple and intuitive approach to dimension reduction of $\mathbf{X}$ for estimating the PS. The approach for DiPS closely resembles a method proposed for estimating mean outcomes in the presence of data missing at random (Hu et al., 2012), except we use the double-score to estimate a PS instead of an outcome model. In contrast to their results, we show that a higher-order kernel is required due to the two-dimensional smoothing, find explicit efficiency gains under misspecification of the outcome model, and consider $p$ diverging with $n$. There is also some similar intuition shared with collaborative DR methods (van der Laan and Gruber, 2010) in that associations with both treatment and outcome are taken into account when estimating a PS. However, DiPS takes a much different approach to estimating the PS. In the following, we introduce the proposed method and consider its asymptotic properties in Sections 2 and 3. A perturbation-resampling method is proposed for inference in Section 4. Simulations and applications to estimating treatment effects in an EMR study and cohort study are presented in Section 5. We conclude with some additional remarks in Section 6.

2 Method

2.1 Notations and Problem Setup

Let $\mathbf{Z}_{i}=(Y_{i},T_{i},\mathbf{X}_{i}^{\sf\tiny T})^{\sf\tiny T}$ be the observed data for the $i$th subject, where $Y_{i}$ is an outcome that could be modeled by a generalized linear model (GLM), $T_{i}\in\{0,1\}$ a binary treatment, and $\mathbf{X}_{i}$ is a $p$-dimensional vector of covariates with support $\mathcal{X}\subseteq\mathbb{R}^{p}$. Here $p$ is allowed to diverge slowly with $n$ such that $log(p)/log(n)\to\nu$, for $\nu\in[0,1)$, which includes the case where $p$ is fixed by taking $\nu=0$. For a given $n$, the observed data consists of independent and identically distributed (iid) observations $\mathscr{D}=\{\mathbf{Z}_{i}:i=1,\ldots,n\}$ drawn from a distribution $\mathbb{P}_{n}$, which potentially may vary with $n$. We suppress the dependence in the notations, implicitly assuming statements involving $\mathbb{P}$ and associated statistical functionals hold for each $n$. Let $Y_{i}^{(1)}$ and $Y_{i}^{(0)}$ denote the counterfactual outcomes had a subject received treatment or control. Based on $\mathscr{D}$, we want to make inferences about the average treatment effect (ATE):

[TABLE]

For identifiability, we require the following standard causal inference assumptions:

[TABLE]

where $\pi_{k}(\mathbf{x})=\mathbb{P}(T=k\mid\mathbf{X}=\mathbf{x})$, for $k=0,1$. The third condition assumes that $\mathbf{X}$ is a sufficient set of covariates such that no unmeasured confounding holds given the entire $\mathbf{X}$. Under these assumptions, $\Delta$ can be identified from the observed data distribution $\mathbb{P}$ through:

[TABLE]

where $\mu_{k}(\mathbf{x})=\mathbb{E}(Y\mid\mathbf{X}=\mathbf{x},T=k)$, for $k=0,1$. We will consider an estimator based on the IPW form that will nevertheless be doubly-robust so that it is consistent under models where either $\pi_{k}(\mathbf{x})$ or $\mu_{k}(\mathbf{x})$ is correctly specified.

2.2 Parametric Models for Nuisance Functions

We consider parametric modeling as a means to reduce the dimensions of $\mathbf{X}$ when estimating the PS. For reference, let $\mathcal{M}_{np}$ be the nonparametric model for the distribution of $\mathbf{Z}$, $\mathbb{P}$, that has no restrictions on $\mathbb{P}$ except requiring the second moment of $\mathbf{Z}$ to be finite. Let $\mathcal{M}_{\pi}\subseteq\mathcal{M}_{np}$ and $\mathcal{M}_{\mu}\subseteq\mathcal{M}_{np}$ respectively denote parametric working models under which:

[TABLE]

where $g_{\pi}(\cdot)$ and $g_{\mu}(\cdot)$ are known link functions, and $\vec{\boldsymbol{\alpha}}=(\alpha_{0},\boldsymbol{\alpha}^{\sf\tiny T})^{\sf\tiny T}\in\Theta_{\boldsymbol{\alpha}}\subseteq\mathbb{R}^{p+1}$ and $\vec{\boldsymbol{\beta}}=(\beta_{0},\beta_{1},\boldsymbol{\beta}_{0}^{\sf\tiny T},\boldsymbol{\beta}_{1}^{\sf\tiny T})^{\sf\tiny T}\in\Theta_{\boldsymbol{\beta}}\subseteq\mathbb{R}^{2p+2}$ are unknown parameters. In (6) slopes are allowed to differ by treatment arms to allow for heterogeneous effects of $T$ for subjects with different $\mathbf{X}$ even with a linear link. When it is reasonable to assume heterogeneity is weak or nonexistent, it may be beneficial for efficiency to restrict $\boldsymbol{\beta}_{0}=\boldsymbol{\beta}_{1}$.

Regardless of the validity of either working model (i.e. whether $\mathbb{P}\in\mathcal{M}_{\pi}\cup\mathcal{M}_{\mu}$), we first obtain estimates of $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}_{k}$’s through adaptive LASSO (Zou, 2006):

[TABLE]

where $\ell_{\pi}(\vec{\boldsymbol{\alpha}};T_{i},\mathbf{X}_{i})$ denotes the log-likelihood for $\vec{\boldsymbol{\alpha}}$ under $\mathcal{M}_{\pi}$ given $T_{i}$ and $\mathbf{X}_{i}$, $\ell_{\mu}(\vec{\boldsymbol{\beta}};\mathbf{Z}_{i})$ is a log-likelihood for $\vec{\boldsymbol{\beta}}$ from a GLM suitable for the outcome type of $Y$ under $\mathcal{M}_{\mu}$ given $\mathbf{Z}_{i}$, $\widetilde{\alpha}_{j}$, $\widetilde{\beta}_{1}$, and $\widetilde{\beta}_{k,j}$ are initial root-$n$ consistent estimates, $\lambda_{\pi,n}$ is a tuning parmaeter such that $n^{1/2}\lambda_{\pi,n}\to 0$ and $n^{(1-\nu)(1+\gamma)/2}\lambda_{\pi,n}\to\infty$, with $\gamma>2\nu/(1-\nu)$, and similarly for $\lambda_{\mu,n}$ (Zou and Zhang, 2009). We specify adaptive LASSO here to estimate the nuisance parameters for concreteness, but use of other penalized likelihood methods can also be justified, so long as they have an oracle property, as in Theorem 2 of Zou (2006) and described below.

Under model (5) and (6), we assume that $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}_{k}$, for $k=0,1$, are sparse. More generally, regardless of whether working models are correct or misspecified, we assume that there exist least false parameters $(\bar{\alpha}_{0},\bar{\boldsymbol{\alpha}}^{\sf\tiny T})^{\sf\tiny T}$ and $(\bar{\beta}_{0},\bar{\beta}_{1},\bar{\boldsymbol{\beta}}_{0}^{\sf\tiny T},\bar{\boldsymbol{\beta}}_{1}^{\sf\tiny T})^{\sf\tiny T}$ (Lu et al., 2012) such that:

[TABLE]

Let $\mathcal{A}_{\boldsymbol{\alpha}}$ and $\mathcal{A}_{\boldsymbol{\beta}_{k}}$ be respective supports for $\bar{\boldsymbol{\alpha}}$ and $\bar{\boldsymbol{\beta}}_{k}$ and let $s_{\boldsymbol{\alpha}}=\left|\mathcal{A}_{\boldsymbol{\alpha}}\right|$ and $s_{\boldsymbol{\beta}_{k}}=\left|\mathcal{A}_{\boldsymbol{\beta}_{k}}\right|$ be the sparsity indices. We further assume $\bar{\boldsymbol{\alpha}}$ and $\bar{\boldsymbol{\beta}}_{k}$ have fixed sparsity such that:

[TABLE]

For any vector $\mathbf{v}$ of length $p$ and any index set $\mathcal{S}\subseteq\{1,2,\ldots,p\}$, let $\mathbf{v}_{\mathcal{S}}$ denote the subvector of $\mathbf{v}$ restricted to elements indexed in $\mathcal{S}$. Assumption (9) is a high-level assumption that would be required for $\widehat{\boldsymbol{\alpha}}$ and $\widehat{\boldsymbol{\beta}}_{k}$ to maintain an oracle property with respect to the least false parameters $\bar{\boldsymbol{\alpha}}$ and $\bar{\boldsymbol{\beta}}_{k}$ under possibly misspecified working models. Under this assumption using arguments similar to those in Lu et al. (2012) and Zou and Zhang (2009) it can be shown that $\mathbb{P}(\widehat{\boldsymbol{\alpha}}_{\mathcal{A}_{\boldsymbol{\alpha}}^{c}}=\boldsymbol{0})\to 1$ and admits an expansion of the form $n^{1/2}(\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}})_{\mathcal{A}_{\boldsymbol{\alpha}}}=n^{-1/2}\sum_{i=1}^{n}\boldsymbol{\Psi}_{i,\mathcal{A}_{\boldsymbol{\alpha}}}+o_{p}(1)$, which would yield the asymptotic normality results of the oracle property, and similarly for $\widehat{\boldsymbol{\beta}}_{k}$. We rely on these results along with (10) to show that the DiPS IPW is asymptotically linear in Theorem 3.1. In regimes where $\nu>0$, (10) models a setting in which a small number of covariates exhibit non-negligible associations with $T$ and $Y$ and a majority of covariates are noise. Assumption 10 may not be required for asymptotic linearity and can potentially be relaxed allowing $s_{\boldsymbol{\alpha}}$ and $s_{\boldsymbol{\beta}_{k}}$ to diverge slowly, for example, if they are $o(n^{1/3})$. We invoke this assumption to avoid complications of a growing support, which may need triangular array asymptotics to accommodate dependence of the support on $n$.

2.3 Double-Index Propensity Score and IPW Estimator

To mitigate the effects of misspecification of (5), one could perform nonparametric smoothing of $T$ over $\widehat{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{X}$ to calibrate the initial PS estimator $g_{\pi}(\widehat{\alpha}_{0}+\mathbf{X}^{\sf\tiny T}\widehat{\boldsymbol{\alpha}})$. We consider smoothing over not only $\widehat{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{X}$ but also $\widehat{\boldsymbol{\beta}}_{k}^{\sf\tiny T}\mathbf{X}$ as well to allow variation in prognostic covariates indexed in $\mathcal{A}_{\boldsymbol{\beta}_{k}}$ to inform this calibration. Such covariates are reduced into $\widehat{\boldsymbol{\beta}}_{k}^{\sf\tiny T}\mathbf{X}$ to allow for nonparametric kernel smoothing in low (two) dimensions. The DiPS estimator for each treatment is:

[TABLE]

where $\widehat{\boldsymbol{\theta}}_{k}=(\widehat{\boldsymbol{\alpha}}^{\sf\tiny T},\widehat{\boldsymbol{\beta}}_{k}^{\sf\tiny T})^{\sf\tiny T}$, $K_{h}(\mathbf{u})=h^{-2}K(\mathbf{u}/h)$, and $K(\mathbf{u})$ is a bivariate $q$-th order kernel function with $q>2$. A higher-order kernel is required here for the asymptotics to be well-behaved, which is the price for estimating the nuisance functions $\pi_{k}(\mathbf{x})$ using two-dimensional smoothing. This allows for the possibility of negative values for $\widehat{\pi}_{k}(\mathbf{x};\widehat{\boldsymbol{\theta}}_{k})$. Nevertheless, $\widehat{\pi}_{k}(\mathbf{x};\widehat{\boldsymbol{\theta}}_{k})$ are nuisance estimates not of direct interest, and we find that such negative PS estimates typically occur infrequently, occurring on average in simulations in $0.01\%$ to $2.10\%$ of observations depending the size of $n$ and $p$ across scenarios where working models are correct or incorrectly specified (Web Appendix D). As they are infrequent and do not appear to compromise the performance of the final estimator, they can potentially be left as is when encountered in practice. Alternatively, methods that discard or trim PS estimates to handle near-violations of positivity, as in Assumption (3), can be considered (Crump et al., 2009). A monotone transformation of the input scores for each treatment $\widehat{\mathbf{S}}_{k}=(\widehat{\boldsymbol{\alpha}},\widehat{\boldsymbol{\beta}}_{k})^{\sf\tiny T}\mathbf{X}$ can be applied prior to smoothing to improve finite sample performance (Wand et al., 1991). In numerical studies, for instance, we applied a probability integral transform based on the normal cumulative distribution function to the standardized scores to obtain approximately uniformly distributed inputs. The components of $\widehat{\mathbf{S}}_{k}$ can also be scaled such that a common bandwidth $h$ can be used for both components of the score.

With $\pi_{k}(\mathbf{x})$ estimated by $\widehat{\pi}_{k}(\mathbf{x};\widehat{\boldsymbol{\theta}}_{k})$, the estimator for $\Delta$ is given by $\widehat{\Delta}=\widehat{\mu}_{1}-\widehat{\mu}_{0}$, where:

[TABLE]

This is the usual normalized IPW estimator, where the PS is estimated by the DiPS. The intuition for double-robustness of the estimator is as follows. Regardless of the validity of either working model, provided the asymptotics are well-behaved, $\widehat{\mu}_{k}$ is consistent for:

[TABLE]

where $\bar{\boldsymbol{\theta}}_{k}=(\bar{\boldsymbol{\alpha}}^{\sf\tiny T},\bar{\boldsymbol{\beta}}_{k}^{\sf\tiny T})^{\sf\tiny T}$, and $\pi_{k}(\mathbf{x};\bar{\boldsymbol{\theta}}_{k})=\mathbb{P}(T_{i}=k\mid\bar{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{X}_{i}=\bar{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{x},\bar{\boldsymbol{\beta}}_{k}^{\sf\tiny T}\mathbf{X}_{i}=\bar{\boldsymbol{\beta}}_{k}^{\sf\tiny T}\mathbf{x})$. Under $\mathcal{M}_{\pi}$, $\pi_{k}(\mathbf{x};\bar{\boldsymbol{\theta}}_{k})=\pi_{k}(\mathbf{x})$ so that the estimand, under the causal assumptions (2)-(4), reduces to:

[TABLE]

On the other hand, under $\mathcal{M}_{\mu}$, $\mathbb{E}(Y_{i}\mid\bar{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{X}_{i}=\bar{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{x},\bar{\boldsymbol{\beta}}_{k}^{\sf\tiny T}\mathbf{X}_{i}=\bar{\boldsymbol{\beta}}_{k}^{\sf\tiny T}\mathbf{x},T_{i}=k)=\mu_{k}(\mathbf{x})$ so that:

[TABLE]

In the following, we show that $\widehat{\mu}_{k}$ (and thus $\widehat{\Delta}$) are asymptotically linear. We then subsequently examine robustness and efficiency properties using the expansion.

3 Asymptotic Robustness and Efficiency Properties

We directly show in Web Appendix B that $\widehat{\mu}_{k}$ is asymptotically linear for $k=0,1$ in general without assuming either of the working models are correct. Let $\bar{\Delta}=\bar{\mu}_{1}-\bar{\mu}_{0}$ and $\widehat{\mathcal{W}}_{k}=n^{1/2}(\widehat{\mu}_{k}-\bar{\mu}_{k})$ for $k=0,1$ so that $n^{1/2}(\widehat{\Delta}-\bar{\Delta})=\widehat{\mathcal{W}}_{1}-\widehat{\mathcal{W}}_{0}$.

Theorem 3.1

Suppose that causal assumptions (2)-(4), the least false parameter and sparsity assumptions (9)-(10) and regularity conditions in Web Appendix A hold. If $log(p)/log(n)\to\nu$ for $\nu\in[0,1)$, then $\widehat{\mu}_{k}$ is asymptotically linear in that it admits the expansion:

[TABLE]

for $k=0,1$, where $\mathbf{u}_{k,\mathcal{A}_{\boldsymbol{\alpha}}}$ and $\mathbf{u}_{k,\mathcal{A}_{\boldsymbol{\beta}_{k}}}$ are deterministic vectors, $\boldsymbol{\Psi}_{i,\mathcal{A}_{\boldsymbol{\alpha}}}$ and $\boldsymbol{\Upsilon}_{i,k,\mathcal{A}_{\boldsymbol{\beta}_{k}}}$ are influence functions from asymptotic expansions of $\widehat{\boldsymbol{\alpha}}_{\mathcal{A}_{\boldsymbol{\alpha}}}$ and $\widehat{\boldsymbol{\beta}}_{k,\mathcal{A}_{\boldsymbol{\beta}_{k}}}$. Under model $\mathcal{M}_{\pi}$ $\mathbf{v}_{k,\mathcal{A}_{\boldsymbol{\beta}_{k}}}=\boldsymbol{0}$, for $k=0,1$. Under $\mathcal{M}_{\pi}\cap\mathcal{M}_{\mu}$, we additionally have that $\mathbf{u}_{k,\mathcal{A}_{\boldsymbol{\alpha}}}=\boldsymbol{0}$, for $k=0,1$.

Proof sketch: $\widehat{\mathcal{W}}_{k}$ can be decomposed as:

[TABLE]

The first term directly contributes to the expansion. The second term is the contribution from re-estimating the PS through kernel smoothing given $\bar{\boldsymbol{\theta}}_{k}$. We apply a V-statistic projection lemma (Newey and McFadden, 1994) to obtain an asymptotically linear representation. The third term can be expanded by Taylor expansion into terms of the form $\mathbf{u}_{k}^{\sf\tiny T}n^{1/2}(\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}})$ and $\mathbf{v}_{k}^{\sf\tiny T}n^{1/2}(\widehat{\boldsymbol{\beta}}-\bar{\boldsymbol{\beta}})$. Applying the selection consistency that $\mathbb{P}(\widehat{\boldsymbol{\alpha}}_{\mathcal{A}_{\boldsymbol{\alpha}}^{c}}=\boldsymbol{0})\to 1$, $\mathbf{u}_{k}^{\sf\tiny T}n^{1/2}(\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}})=\mathbf{u}_{k,\mathcal{A}_{\boldsymbol{\alpha}}}^{\sf\tiny T}n^{1/2}(\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}})_{\mathcal{A}_{\boldsymbol{\alpha}}}+o_{p}(1)$. Lastly, we use that $n^{1/2}(\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}})_{\mathcal{A}_{\boldsymbol{\alpha}}}=n^{-1/2}\sum_{i=1}^{n}\boldsymbol{\Psi}_{i,\mathcal{A}_{\boldsymbol{\alpha}}}+o_{p}(1)$ and work out the forms of the loading vector $\mathbf{u}_{k,\mathcal{A}_{\boldsymbol{\alpha}}}$ and repeat for $\widehat{\boldsymbol{\beta}}_{k}$ to complete the expansion.

Let $\widehat{\Delta}_{dr}=\widehat{\mu}_{1,dr}-\widehat{\mu}_{0,dr}$ denote the usual doubly-robust estimator, as in Equation (9) of Lunceford and Davidian (2004), with the PS $\pi_{k}(\mathbf{x})$ and mean outcome $\mu_{k}(\mathbf{x})$ estimated in the same way as through (7) and (8). The influence function expansion for $\widehat{\Delta}$ in Theorem 3.1 is nearly identical to that of $\widehat{\Delta}_{dr}$. The terms in (13) would be the same except $\pi_{k}(\mathbf{X}_{i};\bar{\boldsymbol{\theta}}_{k})$ and $\mathbb{E}(Y_{i}\mid\bar{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{X}_{i},\bar{\boldsymbol{\beta}}_{k}^{\sf\tiny T}\mathbf{X}_{i},T_{i}=k)$ replaces asymptotic estimates under parametric models. Terms in (14) analogously represent the additional contributions from estimating the nuisance parameters. No contribution from smoothing is incurred provided the bandwidths are suitably chosen. This similarity in the influence functions yields similar robustness and efficiency properties, which are improved upon under model misspecification due to the smoothing.

3.1 Robustness

As a consequence of Theorem 3.1, $\widehat{\Delta}$ is root-$n$ consistent for $\bar{\Delta}$ so that $\widehat{\Delta}-\bar{\Delta}=O_{p}(n^{-1/2})$ provided that $h=O(n^{-\alpha})$ for $\alpha\in(\frac{1}{2q},\frac{1}{4})$. As discussed in Section 2.3, under $\mathcal{M}_{\pi}\cup\mathcal{M}_{\mu}$, $\bar{\Delta}=\Delta$. Hence $\widehat{\Delta}$ is doubly-robust for $\Delta$ in that $\widehat{\Delta}$ is root-$n$ consistent for $\Delta$ under $\mathcal{M}_{\pi}\cup\mathcal{M}_{\mu}$. Beyond this usual form of double-robustness, if the PS model specification is incorrect, we expect the calibration step to at least partially correct for the misspecfication in large samples since $\pi_{k}(\mathbf{x};\bar{\boldsymbol{\theta}}_{k})$ is closer to the true $\pi_{k}(\mathbf{x})$ than the misspecified parametric model $g_{\pi}(\bar{\alpha}_{0}+\bar{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{x})$. Let $\widetilde{\mathcal{M}}_{\pi}$ denote a model under which $\pi_{1}(\mathbf{x})=\widetilde{g}_{\pi}(\boldsymbol{\alpha}^{\sf\tiny T}\mathbf{x})$ for some unknown link function $\widetilde{g}_{\pi}(\cdot)$ and unknown $\boldsymbol{\alpha}\in\mathbb{R}^{p}$, and $\mathbf{X}$ are known to be elliptically distributed such that $\mathbb{E}(\boldsymbol{a}^{\sf\tiny T}\mathbf{X}\mid\boldsymbol{\alpha}_{*}^{\sf\tiny T}\mathbf{X})$ exists and is linear in $\boldsymbol{\alpha}_{*}^{\sf\tiny T}\mathbf{X}$, where $\boldsymbol{\alpha}_{*}$ denotes the true $\boldsymbol{\alpha}$ (e.g. if $\mathbf{X}$ is multivariate normal). By the results of Li and Duan (1989), it can be shown that $\bar{\boldsymbol{\alpha}}=c\boldsymbol{\alpha}_{*}$ for some scalar $c$ under $\widetilde{\mathcal{M}}_{\pi}$. But since $\widehat{\pi}_{k}(\mathbf{x};\widehat{\boldsymbol{\theta}}_{k})$ is consistent for $\pi_{k}(\mathbf{x};\bar{\boldsymbol{\theta}}_{k})=\mathbb{P}(T=k\mid\bar{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{X}=\bar{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{x},\bar{\boldsymbol{\beta}}_{k}^{\sf\tiny T}\mathbf{X}=\bar{\boldsymbol{\beta}}_{k}^{\sf\tiny T}\mathbf{x})$, it recovers $\pi_{k}(\mathbf{x})$ under $\widetilde{\mathcal{M}}_{\pi}$. Consequently, $\widehat{\Delta}$ also has some mild benefits in robustness in that $\widehat{\Delta}-\Delta=O_{p}(n^{-1/2})$ under the slightly larger model $\mathcal{M}_{\pi}\cup\widetilde{\mathcal{M}}_{\pi}\cup\mathcal{M}_{\mu}$. The same phenomenon also occurs when estimating $\boldsymbol{\beta}_{k}$ under misspecification of the link in (6), if we do not assume $\boldsymbol{\beta}_{0}=\boldsymbol{\beta}_{1}$. In this case, if $\widetilde{\mathcal{M}}_{\mu}$ is an analogous model under which $\mu_{1}(\mathbf{x})=\widetilde{g}_{\mu,1}(\boldsymbol{\beta}_{1}^{\sf\tiny T}\mathbf{x})$ and $\mu_{0}(\mathbf{x})=\widetilde{g}_{\mu,0}(\boldsymbol{\beta}_{0}^{\sf\tiny T}\mathbf{x})$ for some unknown link functions $\widetilde{g}_{\mu,0}(\cdot)$ and $\widetilde{g}_{\mu,1}(\cdot)$ and $\mathbf{X}$ are elliptically distributed, then $\widehat{\Delta}-\Delta=O_{p}(n^{-1/2})$ under the slightly larger model $\mathcal{M}_{\pi}\cup\widetilde{\mathcal{M}}_{\pi}\cup\mathcal{M}_{\mu}\cup\widetilde{\mathcal{M}}_{\mu}$. This does not hold when $\boldsymbol{\beta}_{0}=\boldsymbol{\beta}_{1}$, as $T$ is binary so $(T,\mathbf{X}^{\sf\tiny T})^{\sf\tiny T}$ is not exactly elliptically distributed. But the result may still be expected to hold approximately.

3.2 Efficiency

Let the terms contributed to the influence function for $\widehat{\Delta}$ when $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}_{k}$ are known be:

[TABLE]

Under $\mathcal{M}_{\pi}\cap\mathcal{M}_{\mu}$, $\varphi_{i,k}$ is the full influence function for $\widehat{\Delta}$. This is the efficient influence function for $\Delta^{*}$ under $\mathcal{M}_{np}$ at distributions for $\mathbb{P}$ belonging to $\mathcal{M}_{\pi}\cap\mathcal{M}_{\mu}$ when $p$ is fixed (Robins et al., 1994; Tsiatis, 2007), since $\mathbb{E}(Y_{i}\mid\bar{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{X}_{i}=\bar{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{x},\bar{\boldsymbol{\beta}}_{k}^{\sf\tiny T}\mathbf{X}_{i}=\bar{\boldsymbol{\beta}}_{k}^{\sf\tiny T}\mathbf{x},T_{i}=k)=\mu_{k}(\mathbf{x})$ and $\pi_{k}(\mathbf{x};\bar{\boldsymbol{\theta}}_{k})=\pi_{k}(\mathbf{x})$. When $\nu>0$ so that $p$ diverges with $n$, there are no well-established semiparametric efficiency bounds. However with fixed sparsity indices (10), the asymptotic variance still reaches the same bound had $p$ been fixed.

Beyond this characterization of efficiency that parallels that of $\widehat{\Delta}_{dr}$, there are additional benefits of $\widehat{\Delta}$ under $\mathcal{M}_{\pi}\cap\mathcal{M}_{\mu}^{c}$. In this case, akin to $\widehat{\Delta}_{dr}$, estimating $\boldsymbol{\beta}_{k}$ does not contribute to the asymptotic variance since $\mathbf{v}_{k,\mathcal{A}_{\boldsymbol{\beta}_{k}}}=\boldsymbol{0}$, and a similar $n^{1/2}\mathbf{u}_{k,\mathcal{A}_{\boldsymbol{\alpha}}}^{\sf\tiny T}(\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}})_{\mathcal{A}_{\boldsymbol{\alpha}}}$ term is contributed from estimating $\boldsymbol{\alpha}$. The analogous term in the expansion for $\widehat{\Delta}_{dr}$ contributes the negative of a projection of the preceding terms onto the linear span of the score function for $\boldsymbol{\alpha}$, restricted to components in $\mathcal{A}_{\boldsymbol{\alpha}}$, to its influence function (Section 9.1 of Tsiatis (2007)). The same interpretation of the influence function can be adopted for $\widehat{\Delta}$.

Theorem 3.2

Let $\boldsymbol{U}_{\boldsymbol{\alpha}}$ be the score for $\boldsymbol{\alpha}$ under $\mathcal{M}_{\pi}$ and let $[\boldsymbol{U}_{\boldsymbol{\alpha},\mathcal{A}_{\boldsymbol{\alpha}}}]$ denote the linear span of its components indexed in $\mathcal{A}_{\boldsymbol{\alpha}}$. In the Hilbert space of random variables with mean [math] and finite variance $\mathcal{L}_{2}^{0}$ with inner product given by the covariance, let $\Pi\{V\mid\mathcal{S}\}$ denote the projection of some $V\in\mathcal{L}_{2}^{0}$ into a subspace $\mathcal{S}\subseteq\mathcal{L}_{2}^{0}$. If the assumptions required for Theorem 3.1 hold , under $\mathcal{M}_{\pi}$, $\mathbf{u}_{k,\mathcal{A}_{\boldsymbol{\alpha}}}^{\sf\tiny T}n^{1/2}(\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}})_{\mathcal{A}_{\boldsymbol{\alpha}}}=-n^{-1/2}\sum_{i=1}^{n}\Pi\{\varphi_{i,k}\mid[\boldsymbol{U}_{\boldsymbol{\alpha},\mathcal{A}_{\boldsymbol{\alpha}}}]\}+o_{p}(1)$.

The proof is based on simplifying $\mathbf{u}_{k,\mathcal{A}_{\boldsymbol{\alpha}}}$ and is given in Web Appendix B. This result can be used to show that the asymptotic variance of $\widehat{\Delta}$ is lower than that of $\widehat{\Delta}_{dr}$ under $\mathcal{M}_{\pi}\cap\mathcal{M}_{\mu}^{c}$. Based on this result, under $\mathcal{M}_{\pi}\cap\mathcal{M}_{\mu}^{c}$ the influence function for $\widehat{\mu}_{k}$ is $\varphi_{i,k}-\Pi\left\{\varphi_{i,k}\mid[\boldsymbol{U}_{\boldsymbol{\alpha},\mathcal{A}_{\boldsymbol{\alpha}}}]\right\}$, and for the usual DR estimator $\widehat{\mu}_{k,dr}$ is $\phi_{i,k}-\Pi\left\{\phi_{i,k}\mid[\boldsymbol{U}_{\boldsymbol{\alpha},\mathcal{A}_{\boldsymbol{\alpha}}}]\right\}$, where:

[TABLE]

But since $\mathbb{E}(Y_{i}\mid\bar{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{X}_{i}=\bar{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{x},\bar{\boldsymbol{\beta}}_{k}^{\sf\tiny T}\mathbf{X}_{i}=\bar{\boldsymbol{\beta}}_{k}^{\sf\tiny T}\mathbf{x},T_{i}=k)$ better approximates $\mu_{k}(\mathbf{x})$ than the asymptotic estimate under the misspecified parametric model $g_{\mu}(\bar{\beta}_{0}+\bar{\beta}_{1}k+\bar{\boldsymbol{\beta}}_{k}^{\sf\tiny T}\mathbf{x})$, it can then be shown that $\mathbb{E}(\phi_{i,k}^{2})>\mathbb{E}(\varphi_{i,k}^{2})$ for $k=0,1$. Since the influence functions involve projections onto the same space $[\boldsymbol{U}_{\boldsymbol{\alpha},\mathcal{A}_{\boldsymbol{\alpha}}}]$, it can be seen through geometric argument that $\mathbb{E}\left[\varphi_{i,k}-\Pi\left\{\varphi_{i,k}\mid[\boldsymbol{U}_{\boldsymbol{\alpha},\mathcal{A}_{\boldsymbol{\alpha}}}]\right\}\right]^{2}<\mathbb{E}\left[\phi_{i,k}-\Pi\left\{\phi_{i,k}\mid[\boldsymbol{U}_{\boldsymbol{\alpha},\mathcal{A}_{\boldsymbol{\alpha}}}]\right\}\right]^{2}$, so that $\widehat{\Delta}$ is more efficient than $\widehat{\Delta}_{dr}$ under $\mathcal{M}_{\pi}\cap\mathcal{M}_{\mu}^{c}$. We show in the simulation studies that this improvement can lead to substantial efficiency gains under $\mathcal{M}_{\pi}\cap\mathcal{M}_{\mu}^{c}$ in finite samples. These unique robustness and efficiency properties distinguish $\widehat{\Delta}$ from $\widehat{\Delta}_{dr}$ and its variants. We next consider a perturbation scheme to estimate standard errors (SE) and confidence intervals (CI) for $\widehat{\Delta}$.

4 Perturbation Resampling

Although the asymptotic variance of $\widehat{\Delta}$ can be determined through its influence function specified in Theorem (3.1), a direct empirical estimate based on the influence function is infeasible because it involves functionals of $\mathbb{P}$ that are difficult to estimate. Instead we propose a simple perturbation-resampling procedure. Let $\mathcal{G}=\{G_{i}:i=1,\ldots,n\}$ be a set of non-negative iid random variables with unit mean and variance independent of $\mathscr{D}$. The procedure perturbs each “layer” of the estimation of $\widehat{\Delta}$. Let the perturbed estimates of $\vec{\boldsymbol{\alpha}}$ and $\vec{\boldsymbol{\beta}}$ be:

[TABLE]

where $\widetilde{\alpha}_{j}^{*}$ and $\widetilde{\beta}_{j}^{*}$ are perturbed initial estimates obtained from analogously perturbing its estimating equations. The perturbed DiPS estimates are calculated by:

[TABLE]

Lastly the perturbed estimator is given by $\widehat{\Delta}^{*}=\widehat{\mu}_{1}^{*}-\widehat{\mu}_{0}^{*}$ where:

[TABLE]

It can be shown based on arguments in Jin et al. (2001) that the asymptotic distribution of $n^{1/2}(\widehat{\Delta}-\bar{\Delta})$ coincides with that of $n^{1/2}(\widehat{\Delta}^{*}-\widehat{\Delta})\mid\mathscr{D}$. We can thus approximate the SE of $\widehat{\Delta}$ based on the empirical standard deviation or, as a robust alternative, the mean absolute deviations (MAD) of resamples $\widehat{\Delta}^{*}$ and construct CI’s using percentiles of resamples.

5 Numerical Studies

5.1 Simulation Study

We performed extensive simulations to assess the finite sample bias and relative efficiency (RE) of $\widehat{\Delta}$ (DiPS) compared to alternative estimators. We also assessed the performance of the perturbation procedure. Throughout in implementing the adaptive LASSO, we used ridge regression for the initial estimators $\widetilde{\alpha}_{j}$ and $\widetilde{\beta}_{j}$ where the ridge tuning parameter chosen by minimizing the Akaike information criterion (AIC). The adaptive LASSO tuning parameter was chosen by an extended regularized information criterion (Hui et al., 2015), which exhibited relatively good performance for variable selection. We refitted models with selected covariates to reduce bias, as suggested in Hui et al. (2015). The power parameter $\gamma$ was set as $\lceil\frac{2\nu}{1-\nu}\rceil+1$, where $\nu=log(p)/log(n)$. A Gaussian product kernel of order $q=4$ with a plug-in bandwidth at the optimal order (see Discussion) was used for smoothing. For comparison, we considered alternative standard estimators with nuisances estimated by regularization and recently developed methods for estimating ATE that incorporate variable selection: (1) IPW with $\pi_{1}(\mathbf{x})$ estimated by adaptive LASSO (ALAS), (2) $\widehat{\Delta}_{dr}$ with nuisances estimated by adaptive LASSO (DR-ALAS), (3) Modification of $\widehat{\Delta}_{dr}$ in which $\pi_{1}(\mathbf{x})$ and $\mu_{k}(\mathbf{x})$ are estimated by separate one-dimensional kernel smoothing of $T\sim\widehat{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{X}$ and $Y\sim\widehat{\boldsymbol{\beta}}_{k}^{\sf\tiny T}\mathbf{X}$ among those assigned to $T=k$, for $k=0,1$ (DR-SIM), to allow for estimation of single index models (SIM) for $\pi_{1}(\mathbf{x})$ and $\mu_{k}(\mathbf{x})$, (4) Outcome-adaptive LASSO (OAL) (Shortreed and Ertefaie, 2017), (5) Group Lasso and Doubly Robust Estimation (GliDeR) (Koch et al., 2018), (6) Model averaged doubly-robust estimator (MADR) (Cefalu et al., 2017). OAL and GLiDeR were implemented with default settings from code provided in the Supplementary Materials of the respective papers. MADR was implemented using the madr package with $M=500$ Markov chain Monte Carlo (MCMC) iterations to reduce the computations. Throughout the numerical studies, we specified $g_{\pi}(u)=1/(1+e^{-u})$ for $\mathcal{M}_{\pi}$ and $g_{\mu}(u)=u$ with $\boldsymbol{\beta}_{0}=\boldsymbol{\beta}_{1}$ for $\mathcal{M}_{\mu}$ as the working models.

The covariates were generated to approximate the distribution of the covariates from the statins EMR data from Section 5.2. This was done to allow for non-elliptically distributed covariates that mimic the distribution of a real dataset. Initially we generated $\widetilde{\mathbf{X}}\sim N(\widetilde{\boldsymbol{\mu}},\widetilde{\mathbf{\Sigma}})$ where $\widetilde{\boldsymbol{\mu}}$ and $\widetilde{\mathbf{\Sigma}}$ were the empirical mean and covariance matrix of the 15 covariates, which included $9$ binary, $3$ continuous, and $3$ log-transformed count variables. For binary variables we thresholded the corresponding components of $\widetilde{\mathbf{X}}$ so that its mean matched those in $\widetilde{\boldsymbol{\mu}}$, as in $I\{\widetilde{\sigma}_{j}^{-1}(\widetilde{X}_{j}-\widetilde{\mu}_{j})>\Phi^{-1}(1-\widetilde{\mu}_{j})\}$, where $\widetilde{\sigma}_{j}^{2}$ and $\widetilde{\mu}_{j}$ are the empirical variance and mean of the $j$-th covariate and $\Phi(\cdot)$ is the standard normal cumulative distribution function (CDF). Lastly, we centered and standardized to obtain the final covariates $\mathbf{X}=diag(\widetilde{\mathbf{\Sigma}}^{-1/2})(\widetilde{\mathbf{X}}-\widetilde{\boldsymbol{\mu}})$. The pairwise correlations of $\mathbf{X}$ were generally low, mostly ranging between $-.2$ and $.2$ (full correlation matrix reported in Web Appendix C). For settings with $p>15$, we generated independent groups of the $15$ covariates that maintained the correlation structure within each group.

We subsequently focused on a continuous outcome, generating the data according to $T\mid\mathbf{X}\sim Ber\{\pi_{1}(\mathbf{X})\}$ and $Y\mid\mathbf{X},T\sim N\{\mu_{T}(\mathbf{X}),10^{2}\}$. The simulations varied over scenarios where working models were correct or misspecified in which the true $\pi_{1}(\mathbf{x})$ and $\mu_{k}(\mathbf{x})$ are:

[TABLE]

where the coefficients are $\boldsymbol{\alpha}=.01\cdot(1,2,3,4,5,6,\boldsymbol{0}_{3},3,7,0,7,-5,0,\boldsymbol{0}_{p-15})^{\sf\tiny T}$, $\boldsymbol{\alpha}_{[1]}=\boldsymbol{\alpha},\boldsymbol{\alpha}_{[2]}=(.02,.06,.02,.02,-.1,.02,\boldsymbol{0}_{3},-.14,.1,0,-.1,.14,0,\boldsymbol{0}_{p-15})^{\sf\tiny T}$, $\boldsymbol{\zeta}=(\boldsymbol{0}_{6},1,\boldsymbol{0}_{3},1,\boldsymbol{0}_{2},1,0,\boldsymbol{0}_{p-15})^{\sf\tiny T}$, $\boldsymbol{\beta}=(\mathbf{0}_{3},1,.5,.25,.125,.0625,.03125,0,1,.5,0,.25,.125,\boldsymbol{0}_{p-15})^{\sf\tiny T}$, $\boldsymbol{\beta}_{[1]}=(\boldsymbol{0}_{3},.5,0,.5,\boldsymbol{1}_{3},0,1,$ $2,0,1,2,\boldsymbol{0}_{p-15})^{\sf\tiny T}$, $\boldsymbol{\beta}_{[2]}=(\boldsymbol{0}_{3},-1.5,.75,-1.5,\boldsymbol{0}_{3},0,-1.5,-.75,0,1.5,.75,\boldsymbol{0}_{p-15})^{\sf\tiny T}$, and $\mathbf{a}_{m}$ denotes a $1\times m$ vector that has all its elements as $a$. For the misspecified scenarios, either $\mu_{k}(\mathbf{x})$ or $\pi_{1}(\mathbf{x})$ is a double-index model that includes both linear terms in $\mathbf{x}$ and quadratic and two-way interaction terms among $\mathbf{x}$ that are omitted by linear working models. In the misspecified $\mu_{k}(\mathbf{x})$ case, the second index $\boldsymbol{\beta}_{[2]}^{\sf\tiny T}\mathbf{x}$ has some correlation with the PS index $\boldsymbol{\alpha}^{\sf\tiny T}\mathbf{x}$, modeling a situation in which there exist common latent factors not fully captured by a linear outcome model. The outcome model also includes an interaction term between $\mathbf{x}$ and treatment to allow for treatment effect heterogeneity. The parameters are set such that there are $5$ covariates belonging to each of $\mathcal{A}_{\pi}\cap\mathcal{A}_{\mu}$ (i.e. confounders), $\mathcal{A}_{\pi}\cap\mathcal{A}_{\mu}^{c}$ (instruments), and $\mathcal{A}_{\pi}^{c}\cap\mathcal{A}_{\mu}$ (pure prognostic) when $p=15$. The simulations were run for $R=1,000$ repetitions.

Table 1 presents the bias and root mean square error (RMSE) for $n=500,5,000$ when $p=15$. Among the three scenarios considered, the bias for DiPS is small relative to the RMSE and generally diminishes towards zero as $n$ increases, verifying its double-robustness. There remains some minor bias that persists when $n=5,000$ for DiPS that is likely a result of bias from the smoothing, as DR-SIM also incurs similar residual bias. IPW-ALAS and OAL are singly-robust and the bias does not necessary diminish under the misspecified $\pi_{1}(\mathbf{x})$ scenario, although their bias is also minor in the setting considered. MADR exhibited substantial bias under misspecified $\mu_{k}(\mathbf{x})$ scenario that persisted in large samples, possibly due to selecting out confounders with weak outcome associations in its emphasis on selection of prognostic covariates. The results for bias for $p=50,100$ exhibited similar patterns.

Figure 1 presents the RE under the different scenarios for $n=500,5,000$ and $p=15,50,100$. RE was defined as the ratio of the mean square error (MSE) for DR-ALAS relative to that of each estimator, with RE $>1$ indicating greater efficiency compared to DR-ALAS. Under the “both correct” scenario many of the estimators generally exhibit similar efficiency, which can be expected since many are variants of the usual DR estimator and reach the semiparametric efficiency bound. When $n=500$ and $p=60$, there are some slightly greater differences, with GliDeR and MADR leading in efficiency gains, possibly due to differences in the variable selection performance. These differences in efficiency appear to temper when sample size is increased for $n=5,000$ and $p=60$. The results are similar in the “misspecified $\pi_{1}(\mathbf{x})$” scenario, where most estimators exhibited similar efficiency.

In the “misspecified $\mu_{k}(\mathbf{x})$” scenario, DiPS achieves over 70% efficiency gain compared to GliDeR and MADR and over 140% compared to DR-SIM in the large sample setting when $n=5,000$ and $p=15$. This suggests that expected efficiency gains under misspecified outcome models due to the results of Section 3.2 can be substantial. Even if $\pi_{1}(\mathbf{x})$ and $\mu_{k}(\mathbf{x})$ are estimated under a SIM, there are still gains from DiPS when the PS direction $\bar{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{X}$ is informative of the mean outcome beyond $\bar{\boldsymbol{\beta}}^{\sf\tiny T}_{k}\mathbf{X}$. These gains diminish when $p$ is larger relative to $n$, possibly due to imperfect variable selection. Again GLiDeR and MADR achieve the highest efficiency when $n=500$ and $p=60$, notwithstanding the substantial bias of MADR. Thus the performance of DiPS using adaptive LASSO can be somewhat compromised when $p$ is very large relative to $n$ and the variable selection performance is sub-optimal.

Table 2 presents the performance of perturbation for DiPS when $p=15,30$ under correct working models. SEs for DiPS were estimated using the MAD. The empirical SEs (Emp SE), calculated from the sample standard deviations of $\widehat{\Delta}$ over the simulation repetitions, were generally similar to the average of the SE estimates over the repetitions (ASE), despite some overestimation up to 2-15% of the Emp SE. The coverage of the percentile CI’s (Cover) were close to nominal 95% levels but tended to be somewhat conservative.

5.2 Data Example: Effect of Statins on Colorectal Cancer Risk in EMRs

We applied DiPS to assess the effect of statins, a medication for lowering cholesterol levels, on the risk of colorectal cancer (CRC) among patients with inflammatory bowel disease (IBD) identified using data from EMRs of Partners Healthcare. Previous studies have suggested that statins have a protective effect on CRC, but few studies have considered the effect specifically among IBD patients. The EMR cohort consisted of $n=10,817$ IBD patients, including 1,375 statin users. CRC status and statin use were ascertained by the presence of ICD9 diagnosis and prescription codes. We adjusted for $p=15$ covariates as potential confounders, including age, gender, race, smoking status, indication of elevated inflammatory markers, examination with colonoscopy, use of biologics and immunomodulators, subtypes of IBD, disease duration, and presence of primary sclerosing cholangitis (PSC).

For the working model $\mathcal{M}_{\mu}$, we specified $g_{\mu}(u)=1/(1+e^{-u})$ to accomodate the binary outcome. SEs for other estimators were obtained from the MAD over bootstrap resamples. CIs were calculated from percentile intervals. We also calculated a two-sided p-value from a Wald test for the null that statins have no effect, using the point and SE estimates for each estimator. The unadjusted estimate (None) based on difference in means by statins use was also calculated as a reference. The left side of Table 3 shows that, without adjustment, the naive risk difference is estimated to be -0.8% with a SE of 0.4%. The other methods estimated that statins had a protective effect ranging from around -1% to -3% after adjustment for covariates. DiPS and DR-SIM were the most efficient estimators, with DiPS achieving estimated variance that ranged 34% to 61% lower than that of other estimators.

5.3 Data Example: Framingham Offspring Study

The Framingham Offspring Study (FOS) is a cohort study initiated in 1971 that enrolled 5,124 adult children and spouses of the original Framingham Heart Study. The study collected data over time on participants’ medical history, physician examination, and laboratory tests to examine epidemiological and genetic risk factors of cardiovascular disease (CVD). A subset of the FOS participants also have their genotype from the Affymetrix 500K single-nucleotide polymorphism (SNP) array available through the Framingham SNP Health Association Resource (SHARe) on dbGaP. We assessed the effect of smoking on C-reactive protein (CRP) levels, an inflammation marker highly predictive of CVD risk, while adjusting for potential confounders including gender, age, diabetes status, use of hypertensive medication, systolic and diastolic blood pressure measurements, and HDL and total cholesterol measurements, as well as a large number of SNPs in gene regions previously reported to be associated with inflammation or obesity. While the inflmmation-related SNPs are not likely to impact smoking, we include them as prognostic covariates for efficiency. The analysis includes $n=1,892$ individuals with available information on the CRP and the $p=121$ covariates, of which 113 were SNPs.

Since CRP is heavily skewed, we applied a log transformation so that the linear regression model in $\mathcal{M}_{\mu}$ better fits the data. SEs, CIs, and p-values were calculated in the same way as above. The right side of Table 3 shows that different methods agree that smoking significantly increases logCRP. In general, point estimates tended to attenuate after adjusting for covariates since smokers are likely to have other characteristics that increase inflammation. DiPS, DR-SIM, and MADR were among the most efficient, though efficiency gains are tempered in this setting with larger $p$ relative to $n$.

6 Discussion

In this paper we developed a novel IPW estimator for the ATE that accommodates data-driven variable selection through regularized regression. The estimator retains double-robustness and is locally semiparametric efficient when $\nu=0$. By calibrating the initial PS through smoothing, additional gains in efficiency can potentially be achieved in large samples under misspecification of the working outcome model.

In numerical studies, we used the extended regularized information criterion (Hui et al., 2015) to tune adaptive LASSO, which maintains selection consistency when $log(p)/log(n)\to\nu$, for $\nu\in[0,1)$. Other criteria such as cross-validation can also be used and may exhibit better performance in some cases. To obtain a suitable bandwidth $h$, the bandwidth must be selected such that the dominating errors in the influence function, which are of order $O_{p}(n^{1/2}h^{q}+n^{-1/2}h^{-2})$, converges to [math]. This is satisfied for $h=O(n^{-\alpha})$ for $\alpha\in(\frac{1}{2q},\frac{1}{4})$. The optimal bandwidth $h^{*}$ is one that balances these bias and variance terms and is of order $h^{*}=O(n^{-1/(q+2)})$. In practice we use a plug-in estimator $\widehat{h}^{*}=\widehat{\sigma}n^{-1/(q+2)}$, where $\widehat{\sigma}$ is the sample standard deviation of either $\widehat{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{X}_{i}$ or $\widehat{\boldsymbol{\beta}}_{k}^{\sf\tiny T}\mathbf{X}_{i}$, possibly after applying a monotonic transformation. Cross-validation can also be used to select the the smoothing bandwidth.

The adaptive LASSO estimators $\widehat{\boldsymbol{\alpha}}$ and $\widehat{\boldsymbol{\beta}}_{k}$ are not uniformly root-$n$ consistent when the penalty is tuned to achieve consistent model selection (Pötscher and Schneider, 2009), and its oracle properties derived under fixed parameter asymptotics may fail to capture essential features of finite-sample distributions. For example, they are not root-$n$ consistent when the true parameters are of order $O(n^{-1/2})$, if the true signals are relatively weak. The importance of uniform inference also been recently highlighted for treatment effect estimation in high-dimensional settings (Belloni et al., 2013; Farrell, 2015). It would be of interest to consider alternative variable selection approaches beyond those grounded in oracle properties to achieve uniform inference. Another limitation of relying on adaptive LASSO is that when $p$ is large so that $\nu$ is large, a large power parameter $\gamma$ would be required to maintain the oracle properties, leading to an unstable penalty and poor finite sample performance. It would be of interest to consider modifications of the proposed procedure to accommodate high-dimensional settings with $p\gg n$ and more general sparsity assumptions in future work.

Acknowledgements

The authors would like to thank the editor, associate editor, and two referees for their insightful feedback and suggestions. Most of this work was done when the first author was a graduate student at Harvard University. This work was supported by National Institutes of Health grants T32CA009337 and R01HL089778. The views expressed in this article are those of the authors and do not necessarily reflect the views of the Department of Veterans Affairs.

Supporting Information

Web Appendices referenced in Sections 2, 3, and 5 as well as the R code implementing the procedure are available with this paper at the Biometrics website on Wiley Online Library.

\processdelayedfloats

Supporting Information for “Estimating Average Treatment Effects with a Double-Index Propensity Score” by David Cheng, Abhishek Chkrabortty, Ashwin N. Ananthakrishnan, and Tianxi Cai.

These supplementary materials describe the requisite regularity conditions (Web Appendix A) and provides derivations of the two theorems in the main text (Web Appendix B). Web Appendix C reports the correlation matrix used for the covariates in the simulations. Web Appendix D reports simulation results on the proportion of observations with negative values of DiPS across different scenarios.

The following notations will facilitate the derivations. Throughout this Web Appendix, we suppress the $k$ in $\boldsymbol{\beta}_{k}$, $\widehat{\boldsymbol{\beta}}_{k}$, $\bar{\boldsymbol{\beta}}_{k}$, $\boldsymbol{\theta}_{k}$, $\widehat{\boldsymbol{\theta}}_{k}$, and $\bar{\boldsymbol{\theta}}_{k}$ for ease of notation but implicitly understand these quantities to be defined with respect to treatment $k=0,1$ in general. Let $\bar{\mathbf{S}}=(\bar{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{X},\bar{\boldsymbol{\beta}}^{\sf\tiny T}\mathbf{X})^{\sf\tiny T}$ be $\mathbf{X}$ in the directions of $\bar{\boldsymbol{\alpha}}$ and $\bar{\boldsymbol{\beta}}$, regardless of the adequacy of the working models. Let the true density of $\bar{\mathbf{S}}$ at $\mathbf{s}$ be $f(\mathbf{s})$, the propensity score given $\bar{\mathbf{S}}=\mathbf{s}$ for $k=0,1$ be $\pi_{k}(\mathbf{s})=P(T=k\mid\bar{\mathbf{S}}=\mathbf{s})$, and $l_{k}(\mathbf{s})=\pi_{k}(\mathbf{s})f(\mathbf{s})$. Given a $\mathbf{x}\in\mathbb{R}^{p}$, $\boldsymbol{\alpha},\boldsymbol{\beta}\in\mathbb{R}^{p}$, for $\boldsymbol{\theta}=(\boldsymbol{\alpha}^{\sf\tiny T},\boldsymbol{\beta}^{\sf\tiny T})^{\sf\tiny T}$, let:

[TABLE]

For a $p$ length random vector $\mathbf{V}$, let $\mathbf{V}^{\dagger}=(\mathbf{V},\boldsymbol{0}_{p})$ be the $p\times 2$ matrix of the vector augmented by column of zeros on the right and $\mathbf{V}^{\ddagger}=(\boldsymbol{0}_{p},\mathbf{V})$ similarly by a column of zeros on the left. For any two vectors $\mathbf{V}_{i}$ and $\mathbf{V}_{j}$, let $\mathbf{V}_{ji}=\mathbf{V}_{j}-\mathbf{V}_{i}$. Let $K(\mathbf{u})$ be a bivariate symmetric kernel function of order $q>2$, with a finite $q$-th moment. Let $\dot{K}(\mathbf{u})=\partial K(\mathbf{u})/\partial\mathbf{u}$ and $\dot{K}_{h}(\mathbf{v})=h^{-3}\dot{K}(\mathbf{v}/h)$. For any vector $\mathbf{V}$ of length $p$ and $\mathcal{A}\subseteq\{1,2,\ldots,p\}$, with $\left|\mathcal{A}\right|=p_{0}$, let $\mathbf{V}_{\mathcal{A}}$ denote a $p_{0}$-length vector that is $\mathbf{V}$ restricted to coordinates indexed in $\mathcal{A}$. Similarly, let $\mathbf{V}^{\sf\tiny T}_{\mathcal{A}}$ denote $\mathbf{V}^{\sf\tiny T}$ restricted to coordinates indexed in $\mathcal{A}$.

Web Appendix A: Regularity Conditions

(i) $K(\mathbf{u})$ is a bivariate kernel function of order $q>2$, with a finite $q$-th moment. (ii) $K(\mathbf{u})$ is bounded and continuously differentiable with a compact support. (iii) $\dot{K}(\mathbf{u})$ is bounded, integrable, and Lipshitz continuous. (iv) $\mathcal{X}$ is compact. (v) $f(\mathbf{s})$ is bounded and bounded away from [math] over its support. (vi) $f(\mathbf{s})$, $\pi_{k}(\mathbf{s})$, and $E(Y|\bar{\mathbf{S}}=\mathbf{s},T=k)$ for $k=0,1$ are $q$-times continuously differentiable. (vii) $E(\mathbf{X}|\bar{\mathbf{S}}=\mathbf{s})$, $E(\mathbf{X}|\bar{\mathbf{S}}=\mathbf{s},T=k)$, and $E(\mathbf{X}Y|\bar{\mathbf{S}}=\mathbf{s},T=k)$ are continuously differentiable for $k=0,1$. (viii) There exists $0<k_{1}<k_{2}<\infty$ such that the minimum and maximum eigenvalues of $\frac{1}{n}\sum_{i=1}^{n}\mathbf{X}_{i}\mathbf{X}_{i}^{\sf\tiny T}$ around bounded below by $k_{1}$ and above by $k_{2}$. (ix) $\Theta_{\boldsymbol{\alpha}}$ and $\Theta_{\boldsymbol{\beta}}$ are compact. (x) For all $u\in\mathbb{R}$, $1/M\leq g_{\mu}^{\prime}(u)\leq M$ and $\left|g_{\mu}^{\prime\prime}(u)\right|\leq M$ and for some $0<M<\infty$.

Web Appendix B: Derivations of Theorems 1 and 2

Supporting Lemmas

Lemma .1 identifies the stochastic order of a standardized mean when the variance of the observations is of a known order. It will be useful for controlling certain terms that will emerge in the expansion. Lemma .3 shows the uniform convergence rate for kernel smoothing when $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ are fixed, which is a fundamental result used in our approach. Lemma .5 simplifies the average of the gradients of the average of terms that are inversely weighted by the calibrated PS evaluated at the least false parameters. These terms appear repeatedly in subsequent derivations.

Lemma .1

Let $\{X_{i,n}\}$ be a triangular array such that $X_{1,n},\ldots,X_{n,n}$ are iid for each $n\in\mathbb{N}$. Suppose that $\sigma^{2}_{n}=Var(X_{i,n})=O(c_{n}^{2})$, where $c_{n}$ is some positive sequence. Then:

[TABLE]

where $\bar{X}_{n}=n^{-1}\sum_{i=1}^{n}X_{i,n}$ and $\mu_{n}=E(X_{i,n})$.

Proof .2

By Chebyshev’s inequality, for any $k>0$:

[TABLE]

Let $M=\sup_{n\in\mathbb{N}}\sigma^{2}_{n}/c_{n}^{2}$. For any $\epsilon>0$, the desired result is obtained by taking $k=(M/\epsilon)^{1/2}$.

Lemma .3

The uniform convergence rate for two-dimensional smoothing over $\mathbf{X}$ in the directions $\bar{\boldsymbol{\alpha}}$ and $\bar{\boldsymbol{\beta}}$ is given by:

[TABLE]

where $\bar{\boldsymbol{\theta}}=(\bar{\boldsymbol{\alpha}}^{\sf\tiny T},\bar{\boldsymbol{\beta}}^{\sf\tiny T})^{\sf\tiny T}$, $\pi_{k}(\mathbf{x};\bar{\boldsymbol{\theta}})=P(T=k\mid\bar{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{X}=\bar{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{x},\bar{\boldsymbol{\beta}}^{\sf\tiny T}\mathbf{X}=\bar{\boldsymbol{\beta}}^{\sf\tiny T}\mathbf{x})$, and:

[TABLE]

Proof .4

Smoothing over $\mathbf{X}$ in the directions of $\bar{\boldsymbol{\alpha}}$ and $\bar{\boldsymbol{\beta}}$ is the same as a two-dimensional kernel smoothing since $\bar{\boldsymbol{\alpha}}$ and $\bar{\boldsymbol{\beta}}$ are fixed. See, for example, Hansen (2008b) for the derivation of uniform convergence rates for $d$-dimensional smoothing.

Lemma .5

Let $g(\mathbf{Z})$ denote a real-valued square-integrable transformation of the data $\mathbf{Z}=(\mathbf{X},T,Y)^{\sf\tiny T}$. Under the above regularity conditions and that $E\{g(\mathbf{Z})|\bar{\mathbf{S}}=\mathbf{s}\}$ and $E\{\mathbf{X}g(\mathbf{Z})|\bar{\mathbf{S}}=\mathbf{s}\}$ are continuous in $\mathbf{s}$:

[TABLE]

where $b_{n}=n^{-1/2}h^{-1}+n^{-1}h^{-3}$, for $k=0,1$.

Proof .6

We will show the first equality for the gradient with respect to $\boldsymbol{\alpha}$, with the second equality being analogous. First note each of the gradients can be written:

[TABLE]

Consequently, the average of the gradients can be written:

[TABLE]

where we make repeated use of the uniform convergence of $\widehat{l}_{k}(\mathbf{X}_{i};\bar{\boldsymbol{\theta}})$ and $\widehat{f}(\mathbf{X}_{i};\bar{\boldsymbol{\theta}})$ to $l_{k}(\bar{\mathbf{S}}_{i})$ and $f(\bar{\mathbf{S}}_{i})$, $\mathcal{E}_{n}$ is a term of the same order as the main term so that $O_{p}(a_{n}\mathcal{E}_{n})$ will be a negligible lower-order term, and use that $l_{k}(\mathbf{s})$ is bounded over its support in the last equality. To facilitate application of the V-statistic projection lemma, define:

[TABLE]

We now further evaluate each term. The first term can be simplified through a change-of-variables:

[TABLE]

where the last step follows from bounding the integrand. Similarly for the second term:

[TABLE]

where again the last step follows from bounding the integrand. Now, $\boldsymbol{\varepsilon}_{2,k}=O_{p}(n^{-1}h^{-3})$ from bounding the terms in the expectation, except for $g(\mathbf{Z}_{i})$. The projection lemma thus yields:

[TABLE]

for $k=0,1$, where the last line follows from application of Lemma .1. Re-arrangement of terms and collecting the dominant errors yield the desired result.

Expansion of Normalization Constant

We will first show the normalization constant is $1$ up to some lower order terms, which will allow us to account for the normalization in the expansion. The approach for the analysis parallels that of the main expansion. First note that:

[TABLE]

where:

[TABLE]

The second term is of order:

[TABLE]

where the last step follows from uniform convergence of $\widehat{\pi}_{k}(\mathbf{X}_{i};\bar{\boldsymbol{\beta}})$ to $\pi_{k}(\mathbf{X}_{i};\bar{\boldsymbol{\theta}})$ and noting the remaining sum is $O_{p}(1)$ plus some lower-order term. The third term can be written:

[TABLE]

where the last equality uses that that $\dot{K}(\mathbf{u})$ is Lipshitz continuous and that $\mathcal{E}_{\boldsymbol{\alpha},n}$, $\mathcal{E}_{\boldsymbol{\beta},n}$, and $\mathcal{E}_{\boldsymbol{\alpha}\boldsymbol{\beta},n}$ are terms of the same order as $n^{-1}\sum_{i=1}^{n}\frac{\partial}{\partial\boldsymbol{\alpha}^{\sf\tiny T}}\widehat{\pi}_{k}(\mathbf{X}_{i};\bar{\boldsymbol{\alpha}};\bar{\boldsymbol{\beta}})^{-1}$ so that the error terms will be negligible lower-order terms. Applying Lemma .5, we can simplify:

[TABLE]

Further simplifying the expectation we have:

[TABLE]

where the last step follows from bounding terms in the integrand. Similarly:

[TABLE]

Collecting all the results:

[TABLE]

Main Results

The approach for the expansion will be to decompose $\widehat{\mathcal{W}}_{k}$ into terms representing the variability contributed from smoothing, with known $\boldsymbol{\theta}$, and from estimating $\boldsymbol{\theta}$. The term contributed from smoothing is written in terms of a V-statistic and analyzed using a V-statistic projection lemma (Lemma 8.4 of Newey and McFadden (1994)). The term contributed from estimating $\boldsymbol{\theta}$ is analyzed applying arguments for the oracle properties of adaptive LASSO estimators from Zou and Zhang (2009),Hui et al. (2015), and Lu et al. (2012). First note that:

[TABLE]

where the second step follows from the result in Web Appendix C. Define:

[TABLE]

where:

[TABLE]

We now proceed to further expand the second and third terms. For the second term:

[TABLE]

where the last equality follows from repeated use of uniform convergence of of $\widehat{l}(\mathbf{X}_{i};\bar{\boldsymbol{\theta}})$ and $\widehat{f}(\mathbf{X}_{i};\bar{\boldsymbol{\theta}})$ to $l(\mathbf{X}_{i};\bar{\boldsymbol{\theta}})$ and $f(\mathbf{X}_{i};\bar{\boldsymbol{\theta}})$ and that $n^{-1/2}\sum_{i=1}^{n}\frac{I(T_{i}=k)}{\pi_{k}(\mathbf{X}_{i};\bar{\boldsymbol{\theta}})}\frac{Y_{i}-\bar{\mu}}{l_{k}(\mathbf{X}_{i};\bar{\boldsymbol{\theta}})}=O_{p}(n^{1/2})$. Thus:

[TABLE]

with a centered and a non-centered V-statistic:

[TABLE]

To facilitate application of the projection lemma, let:

[TABLE]

We now evaluate each term. The first term can be simplified through change-of-variables:

[TABLE]

where $\left\lVert\bar{\mathbf{S}}_{j}^{*}-\bar{\mathbf{S}}_{j}\right\rVert\leq h\left\lVert\boldsymbol{\psi}_{j}\right\rVert$ and:

[TABLE]

For the second term, due to the centering:

[TABLE]

For the remaining terms:

[TABLE]

where the order of the second error can be obtained from bounding terms inside the expectation. Now, we apply the projection lemma to find that:

[TABLE]

We used that $\pi(\mathbf{s})$ and $E(Y|\bar{\mathbf{S}}=\mathbf{s},T=k)$ are $q$-times continuously differentiable to bound the remainder error term from $m_{1,ct,2,k}(\mathbf{Z}_{j})$.

We now repeat a similar analysis for $\widetilde{\mathcal{W}}_{nc,2,k}$. Let:

[TABLE]

The first term is:

[TABLE]

where $\bar{\mathbf{S}}_{j}^{*}$ is such that $\left\lVert\bar{\mathbf{S}}_{j}^{*}-\bar{\mathbf{S}}_{j}\right\rVert\leq h\left\lVert\boldsymbol{\psi}_{j}\right\rVert$ and the last equality can be obtained through bounding $\frac{\partial}{\partial\boldsymbol{s}^{\otimes q}}\pi_{k}(\bar{\mathbf{S}}_{j}^{*})$ and $\xi_{k}(h\boldsymbol{\psi}_{j}+\bar{\mathbf{S}}_{j})$. Similarly, for the second term:

[TABLE]

where $\bar{\mathbf{S}}_{i}^{*}$ is such that $\left\lVert\bar{\mathbf{S}}_{i}^{*}-\bar{\mathbf{S}}_{i}\right\rVert\leq h\left\lVert\boldsymbol{\psi}_{i}\right\rVert$ and the last equality could be obtained through bounding $\frac{\partial}{\partial\boldsymbol{s}^{\otimes q}}\pi_{k}(\bar{\mathbf{S}}_{i}^{*})$ and $f(h\boldsymbol{\psi}_{i}+\bar{\mathbf{S}}_{i})$. The errors are:

[TABLE]

where the order of the second error can be obtained from bounding terms inside the expectation. Application of the projection lemma now yields:

[TABLE]

where we use that $Var\{m_{1,nc,2,k}(\mathbf{Z}_{j})\}=O(h^{2q})$ and $Var\{m_{2,nc,2,k}(\mathbf{Z}_{i})\}=O(h^{2q})$ and apply Lemma .1. We now evaluate $m_{nc,2,k}$:

[TABLE]

where $\mathbf{s}^{*}$ is such that $\left\lVert\mathbf{s}^{*}-\mathbf{s}_{1}\right\rVert\leq h\left\lVert\boldsymbol{\psi}_{1}\right\rVert$, and the last equality follows from bounding $\frac{\partial}{\partial\boldsymbol{s}^{\otimes q}}\pi_{k}(\mathbf{s}_{1}^{*})$ and $f(h\boldsymbol{\psi}_{1}+\mathbf{s}_{1})$. We have now have that:

[TABLE]

We now proceed to expand $\widetilde{\mathcal{W}}_{3,k}$. We then first analyze the gradients in general, under model $\mathcal{M}_{\pi}$, and under model $\mathcal{M}_{\pi}\cap\mathcal{M}_{\mu}$, using Lemma .5. First note that:

[TABLE]

using that $\frac{\partial}{\partial\boldsymbol{\alpha}^{\sf\tiny T}}\widehat{\pi}_{k}(\mathbf{X}_{i};\boldsymbol{\theta})^{-1}$ and $\frac{\partial}{\partial\boldsymbol{\beta}^{\sf\tiny T}}\widehat{\pi}_{k}(\mathbf{X}_{i};\boldsymbol{\theta})^{-1}$ are Lipshitz continuous in $\boldsymbol{\theta}$. Now it can be shown that $\mathbb{P}\{\widehat{\boldsymbol{\alpha}}_{\mathcal{A}_{\boldsymbol{\alpha}}^{c}}=\boldsymbol{0}\}\to 1$ and $\mathbb{P}\{\widehat{\boldsymbol{\beta}}_{\mathcal{A}_{\boldsymbol{\beta}}^{c}}=\boldsymbol{0}\}\to 1$, using arguments from Hui et al. (2015) and Zou and Zhang (2009) when working models are correctly specified. It can also be shown that this still holds under misspecified models, provided that the least false parameters $\bar{\boldsymbol{\alpha}}$ and $\bar{\boldsymbol{\beta}}$ exist and are sparse, using arguments similar to those from Lu et al. (2012) and Zou and Zhang (2009). Let:

[TABLE]

Using that $\mathbb{P}\{\widehat{\boldsymbol{\alpha}}_{\mathcal{A}_{\boldsymbol{\alpha}}^{c}}=\boldsymbol{0}\}\to 1$ and $\mathbb{P}\{\widehat{\boldsymbol{\beta}}_{\mathcal{A}_{\boldsymbol{\beta}}^{c}}=\boldsymbol{0}\}\to 1$, we have that $\mathbf{u}_{k,n,\mathcal{A}_{\boldsymbol{\alpha}}^{c}}^{\sf\tiny T}n^{1/2}(\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}})_{\mathcal{A}_{\boldsymbol{\alpha}}^{c}}=o_{p}(1)$, $\mathbf{v}_{k,n,\mathcal{A}_{\boldsymbol{\beta}}^{c}}^{\sf\tiny T}n^{1/2}(\widehat{\boldsymbol{\beta}}-\bar{\boldsymbol{\beta}})_{\mathcal{A}_{\boldsymbol{\beta}}^{c}}=o_{p}(1)$, $n^{1/2}\left\lVert(\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}})_{\mathcal{A}_{\boldsymbol{\alpha}}^{c}}\right\rVert^{2}=o_{p}(1)$, and $n^{1/2}\left\lVert(\widehat{\boldsymbol{\beta}}-\bar{\boldsymbol{\beta}})_{\mathcal{A}_{\boldsymbol{\beta}}^{c}}\right\rVert^{2}=o_{p}(1)$, so that:

[TABLE]

Applying Lemma .5 to the gradient restricted to the respective active sets:

[TABLE]

We now verify that $\mathbf{u}_{k,\mathcal{A}_{\boldsymbol{\alpha}}}$ and $\mathbf{v}_{k,\mathcal{A}_{\boldsymbol{\beta}}}$ are $O_{p}(1)$ in general. First note that:

[TABLE]

where:

[TABLE]

Each of the four terms in $\mathbf{u}_{k,\mathcal{A}_{\boldsymbol{\alpha}}}^{\sf\tiny T}$ can be simplified through change-of-variables:

[TABLE]

For some vector $\mathbf{u}$, let $\dot{K}(\mathbf{u})=(\dot{K}(\mathbf{u})_{\{1\}},\dot{K}(\mathbf{u})_{\{2\}})^{\sf\tiny T}$ be the partial derivatives of $K(\mathbf{u})$ with respect to the first and second components of $\mathbf{u}$, evaluated at $\mathbf{u}$. Similarly, for some $\mathbf{s}_{1}$ and $\mathbf{s}_{2}$, let $\left\{\boldsymbol{\eta}(\mathbf{s}_{1})_{1,u,k}\boldsymbol{\eta}(\mathbf{s}_{2})_{1,u,k}\right\}_{\{(i,j)\}}$ denote the $(i,j)$-th element of $\boldsymbol{\eta}(\mathbf{s}_{1})_{1,u,k}\boldsymbol{\eta}(\mathbf{s}_{2})_{1,u,k}$ evaluated at $\mathbf{s}_{1}$ and $\mathbf{s}_{2}$, for $i=1,2$ and $j=1,\ldots,p+1$. Applying integration by parts, the $j$-th element of the above expectation is:

[TABLE]

where the second to last and last equalities can be shown by bounding terms using that $\mathbb{E}(Y_{i}\mid\bar{\mathbf{S}}_{i}=\mathbf{s},T_{i}=k)$, $\pi_{k}(\mathbf{s})$, $f(\mathbf{s})$, $\mathbb{E}(Y_{i}\mathbf{X}_{i}\mid\bar{\mathbf{S}}_{i}=\mathbf{s},T_{i}=k)$, $\mathbb{E}(\mathbf{X}_{i}\mid\bar{\mathbf{S}}_{i}=\mathbf{s},T_{i}=k)$ are differentiable in $\mathbf{s}$ for $k=0,1$, $\mathbb{E}(\mathbf{X}\mid\bar{\mathbf{S}}=\mathbf{s})$ is continuous in $\mathbf{s}$, $\mathcal{X}$ is compact, and $K(\mathbf{u})$ is a kernel function. Consequently:

[TABLE]

Applying the same argument it can be shown that $\mathbf{v}_{k,\mathcal{A}_{\boldsymbol{\beta}}}^{\sf\tiny T}=O(1)$ for $k=0,1$ as well.

We now consider simplifying $\mathbf{v}_{k,\mathcal{A}_{\beta}}^{\sf\tiny T}$ under $\mathcal{M}_{\pi}$. First we note that under $\mathcal{M}_{\pi}$, $T\perp\!\!\!\perp\mathbf{X}\mid\bar{\boldsymbol{\alpha}}^{\sf\tiny T}\mathbf{X}$. This implies that $T\perp\!\!\!\perp\mathbf{X}\mid\bar{\mathbf{S}}$. Applying this and similar calculations used above for $\mathbf{u}_{k,\mathcal{A}_{\boldsymbol{\alpha}}}^{\sf\tiny T}$:

[TABLE]

We now evaluate:

[TABLE]

where $\mathbf{s}^{*}$ is such that $\left\lVert\mathbf{s}^{*}-\mathbf{s}_{1}\right\rVert\leq h\left\lVert\boldsymbol{\psi}_{1}\right\rVert$ and we use that $f(s)$ and $\mathbb{E}(X\mid\bar{\mathbf{S}}=\mathbf{s})$ are continuously differentiable and that $\pi_{k}(\mathbf{s})$ is twice continuously differentiable, $\mathcal{X}$ is compact, $\dot{K}(\mathbf{u})$ is bounded and integrable, to bound terms in the remainder. After some re-arrangement, this can be further simplified:

[TABLE]

where the second equality follows from that $\int\boldsymbol{\psi}_{1}\dot{K}(\boldsymbol{\psi}_{1})^{\sf\tiny T}d\boldsymbol{\psi}_{1}=-\boldsymbol{I}_{2\times 2}$ by integration by parts. Let the partial derivatives of $\pi_{k}(\mathbf{s})$ with respect to $\mathbf{s}$, evaluated at $\mathbf{s}$, be denoted by $\partial\pi_{k}(\mathbf{s})/\partial\mathbf{s}^{\sf\tiny T}=(\partial\pi_{k}(\mathbf{s})/\partial s_{1},\partial\pi_{k}(\mathbf{s})/\partial s_{2})$. Under $\mathcal{M}_{\pi}$ when the PS model is correct, $\partial\pi_{k}(\mathbf{s})/\partial s_{2}=0$ since $\pi_{k}(\mathbf{s})$ would depend only on the first argument. The last equality follows from noting this and that the first row of $\mathbf{X}^{\ddagger\sf\tiny T}_{i}$ is $\boldsymbol{0}^{\sf\tiny T}$.

Finally, we consider the case under $\mathcal{M}_{\pi}\cap\mathcal{M}_{\mu}$. In this case we have not only that $T\perp\!\!\!\perp\mathbf{X}\mid\bar{\mathbf{S}}$ but also $\mathbb{E}(Y\mid\bar{\mathbf{S}},T=k,\mathbf{X})=g_{\mu}(\bar{\beta}_{0}+\bar{\beta}_{1}k+\bar{\boldsymbol{\beta}}^{\sf\tiny T}\mathbf{X})=\mathbb{E}(Y\mid\bar{\mathbf{S}},T=k)$. Thus in this case:

[TABLE]

Consequently, continuing from an analogous expression for $\mathbf{u}_{k,\mathcal{A}_{\boldsymbol{\alpha}}}$ from (140):

[TABLE]

Evaluating the expression, we obtain that:

[TABLE]

where $\mathbf{s}_{1}^{*}$ and $\mathbf{s}_{1}^{**}$ are values such that $\left\lVert\mathbf{s}_{1}^{*}-\mathbf{s}_{1}\right\rVert\leq h\left\lVert\boldsymbol{\psi}_{1}\right\rVert$ and $\left\lVert\mathbf{s}_{1}^{**}-\mathbf{s}_{1}\right\rVert\leq h\left\lVert\boldsymbol{\psi}_{1}\right\rVert$. The last equality can be shown by bounding terms inside the integral by using that $\pi_{k}(\mathbf{s})$ is continuously differentiable and bounded away from [math], $\mathbb{E}(Y-\bar{\mu}_{k}\mid\bar{\mathbf{S}}=\mathbf{s},T=k)$ is continuous, $\mathbb{E}(\mathbf{X}\mid\bar{\mathbf{S}}=\mathbf{s})$ is continuously differentaible, $f(\mathbf{s})$ is continuous, and $\mathcal{X}$ is compact. The same argument can be applied to show that $\mathbf{v}_{k,\mathcal{A}_{\boldsymbol{\beta}}}^{\sf\tiny T}=O(h)$ for $k=0,1$, under $\mathcal{M}_{\pi}\cap\mathcal{M}_{\mu}$.

We now collect all the results in the main expansion:

[TABLE]

where $\mathbf{u}_{k,\mathcal{A}_{\boldsymbol{\alpha}}}^{\sf\tiny T}$ and $\mathbf{v}_{k,\mathcal{A}_{\boldsymbol{\beta}}}$ are deterministic vectors such that, for $k=0,1$, $\mathbf{v}_{k,\mathcal{A}_{\boldsymbol{\beta}}}=\boldsymbol{0}$ under $\mathcal{M}_{\pi}$ and $\mathbf{u}_{k,\mathcal{A}_{\boldsymbol{\alpha}}}=\mathbf{v}_{k,\mathcal{A}_{\boldsymbol{\beta}}}=\boldsymbol{0}$ under $\mathcal{M}_{\pi}\cap\mathcal{M}_{\mu}$. The final form of the expansion by using that $\widehat{\boldsymbol{\alpha}}_{\mathcal{A}_{\boldsymbol{\alpha}}}$ and $\widehat{\boldsymbol{\beta}}_{\mathcal{A}_{\boldsymbol{\beta}}}$ admit an asymptotically linear expansion, using arguments similar to those from Zou and Zhang (2009) and Lu et al. (2012) so that:

[TABLE]

where $\boldsymbol{\Psi}_{i,\mathcal{A}_{\boldsymbol{\alpha}}}=\mathbb{E}(\boldsymbol{U}_{\boldsymbol{\alpha},\mathcal{A}_{\boldsymbol{\alpha}}}\boldsymbol{U}_{\boldsymbol{\alpha},\mathcal{A}_{\boldsymbol{\alpha}}}^{\sf\tiny T})^{-1}\boldsymbol{U}_{\boldsymbol{\alpha},i,\mathcal{A}_{\boldsymbol{\alpha}}}$ and $\boldsymbol{\Upsilon}_{i,\mathcal{A}_{\boldsymbol{\beta}}}=\mathbb{E}(\boldsymbol{U}_{\boldsymbol{\beta},\mathcal{A}_{\boldsymbol{\beta}}}\boldsymbol{U}_{\boldsymbol{\beta},\mathcal{A}_{\boldsymbol{\beta}}}^{\sf\tiny T})^{-1}\boldsymbol{U}_{\boldsymbol{\beta},i,\mathcal{A}_{\boldsymbol{\beta}}}$, $\boldsymbol{U}_{\boldsymbol{\alpha},i}=\mathbf{X}_{i}\left\{T_{i}-\pi_{1}(\mathbf{X}_{i};\bar{\alpha}_{0},\bar{\boldsymbol{\alpha}})\right\}$, and $\boldsymbol{U}_{\boldsymbol{\beta},i}=\mathbf{X}_{i}\left\{Y_{i}-\mu_{T_{i}}(\mathbf{X}_{i};\bar{\beta}_{0},\bar{\beta}_{1},\bar{\boldsymbol{\beta}})\right\}$, $\pi_{1}(\mathbf{x};\alpha_{0},\boldsymbol{\alpha})=g_{\pi}(\alpha_{0}+\boldsymbol{\alpha}^{\sf\tiny T}\mathbf{x})$, and $\mu_{k}(\mathbf{x};\beta_{0},\beta_{1},\boldsymbol{\beta})=g_{\mu}(\beta_{0}+\beta_{1}k+\boldsymbol{\beta}^{\sf\tiny T}\mathbf{x})$.

We next verify Theorem 2 to characterize the influence function contribution from estimating the PS. Since $[\boldsymbol{U}_{\boldsymbol{\alpha},\mathcal{A}_{\boldsymbol{\alpha}}}]$ is a finite dimensional subspace of $\mathcal{L}_{2}^{0}$ spanned by the components of $\boldsymbol{U}_{\boldsymbol{\alpha},\mathcal{A}_{\boldsymbol{\alpha}}}$, the projection of $\varphi_{i,k}$ onto it is given by population least squares:

[TABLE]

As discussed above, $n^{1/2}(\widehat{\boldsymbol{\alpha}}-\bar{\boldsymbol{\alpha}})_{\mathcal{A}_{\boldsymbol{\alpha}}}=\mathbb{E}(\boldsymbol{U}_{\boldsymbol{\alpha},\mathcal{A}_{\boldsymbol{\alpha}}}\boldsymbol{U}_{\boldsymbol{\alpha},\mathcal{A}_{\boldsymbol{\alpha}}}^{\sf\tiny T})^{-1}n^{-1/2}\sum_{i=1}^{n}\boldsymbol{U}_{\boldsymbol{\alpha},i,\mathcal{A}_{\boldsymbol{\alpha}}}+o_{p}(1)$. It thus suffices to show that $\mathbf{u}_{k,\mathcal{A}_{\boldsymbol{\alpha}}}^{\sf\tiny T}=-\mathbb{E}(\varphi_{i,k}\boldsymbol{U}_{\boldsymbol{\alpha},\mathcal{A}_{\boldsymbol{\alpha}}}^{\sf\tiny T})+o(1)$.

We proceed by simplifying the covariance term:

[TABLE]

First consider the $k=1$ case. Using that $\pi_{1}(\mathbf{X};\bar{\alpha}_{0},\bar{\boldsymbol{\alpha}})=\pi_{1}(\bar{\mathbf{S}})$ and $T\perp\!\!\!\perp\mathbf{X}\mid\bar{\mathbf{S}}$ under $\mathcal{M}_{\pi}$:

[TABLE]

where $\mathbf{u}_{k,\mathcal{A}_{\boldsymbol{\alpha}}}^{\sf\tiny T}$ has same form as derived in (151) for $\mathbf{v}_{k,\mathcal{A}_{\boldsymbol{\beta}}}$, except that $\mathbf{X}^{\ddagger\sf\tiny T}_{i}$ is replaced by $\mathbf{X}^{\dagger\sf\tiny T}_{i}$.

In the $k=0$ case, again using $\pi_{1}(\mathbf{X};\bar{\alpha}_{0},\bar{\boldsymbol{\alpha}})=\pi_{1}(\bar{\mathbf{S}})$ and $T\perp\!\!\!\perp\mathbf{X}\mid\bar{\mathbf{S}}$ under $\mathcal{M}_{\pi}$:

[TABLE]

Web Appendix C: Covariate Correlation Matrix in Simulations

As described in the Simulation Study Section 4.1, the covariates were generated as $\mathbf{X}=diag(\widetilde{\mathbf{\Sigma}}^{-1/2})(\widetilde{\mathbf{X}}-\widetilde{\mu})$, and here we report the values of its covariance matrix $diag(\widetilde{\mathbf{\Sigma}}^{-1/2})\widetilde{\mathbf{\Sigma}}diag(\widetilde{\mathbf{\Sigma}}^{-1/2})$ for each group of $15$ covariates. The covariates are ordered as ulcerative colitis disease subtype, female gender, use of anti-TNF therapy, use of immunomodulator, primary sclerosing cholangitis (PSC), elevated C-reactive protein, race1, race2, counts of ever smoking from NLP, counts of current smoking from NLP, counts of never smoking from NLP, utilization score, disease duration, and age. For simulations when $p>15$, we used block diagonal matrix where we repeated this correlation structure for each group of $15$ covariates.

[TABLE]

Web Appendix D: Proportion of Observations with Negative Values of DiPS

As discussed in Section 2.1, the following reports the proportion of total number of observations ($n$) in simulations that have negative values of DiPS, by varying size of $p$ and model specification scenario. The results are averaged over $R=1,000$ repetitions.

[TABLE]