Developing and Testing New Methods for Estimating Treatment Effectiveness in Observational Studies Using High-Dimension Data

Background:

Appropriate causal inference methods are required for comparative effectiveness research to produce valid and relevant findings from observational data. Two prominent classes of such methods are based on unconfoundedness or instrumental variable (IV) assumptions. Although extensive research has been done, it remains highly challenging to estimate propensity scores (PSs) and regression functions and to perform subsequent inference about average treatment effects (ATEs). The conventional approach employs an iterative process of model building and fitting, depending on ad hoc modeling choices, where statistical uncertainty is difficult to quantify. Recently, various methods have been proposed that apply off-the-shelf machine learning algorithms but either ignore statistical inference or invoke strong smoothness assumptions to justify consistent estimation of regression functions and subsequent statistical inference about treatment effects.

Objectives:

The objective of our research is to develop and evaluate a new set of statistically rigorous, numerically tractable, and pragmatic methods for drawing inferences about ATEs using PSs or IVs, while fitting PS and regression models with a large number of regressors.

Methods:

We propose regularized calibrated estimation and model-assisted inferences about ATEs under the assumption of no unmeasured confounding or local ATEs under IV assumptions in high-dimensional settings. We derived numerical algorithms to implement the methods, conducted simulation studies to evaluate the methods, and investigated empirical applications of the proposed methods, compared with existing methods.

Results:

The proposed methods are shown to yield valid statistical inference (ie, CIs and hypothesis testing) about treatment parameters under weak technical conditions in high-dimensional settings. Simulation studies and empirical applications demonstrate the advantages of the proposed methods compared with related methods.

Conclusions:

We developed new statistical methods and theory using PSs or IVs for causal inference. Using the proposed methods, PS and regression models can be fitted with a possibly large number of regressors, including main effects and interactions of the covariates, and CIs, and hypothesis tests can be obtained about treatment effects in a numerically tractable and statistically principled manner. Our methods are implemented in the publicly released R package RCAL.

Limitations:

Currently used methods handle cross-sectional studies and assume that the treatment and instrument are binary. It is desirable to extend our methods to handle multivalued treatments and instruments and to analyze longitudinal and survival data.

Causal Inference

PCORI was established with the mission to help patients, clinicians, and other health care stakeholders make better-informed health care decisions and improve health care delivery and outcomes. One of the most frequently sought types of information is whether a specific treatment will cause deterioration or improvement in the outcomes of interest compared with other treatments. To provide such information, it is important to conduct comparative effectiveness research (CER). However, a fundamental challenge in nonrandomized, observational studies is that it can be enormously difficult to disentangle the effects of the treatment from other risk factors, known as confounding variables (or covariates), which might influence the health outcomes and, at the same time, vary between patients who had the intervention treatment and those who had the comparator treatment. The subject of how to address confounding and draw inferences about treatment effects is broadly called causal inference, spanning multiple disciplines, including epidemiology, econometrics, and statistics.

As discussed in the PCORI Methodology Report (PCORI Methodology Committee, 2021), Section III (8), appropriate causal inference methods are required for CER to produce valid and relevant findings from observational data, where treatment choices are self-selected rather than randomly assigned. There is extensive literature on various analytical methods for causal inference (Morgan & Winship 2014; van der Laan & Robins 2003). The PCORI Methodology Report explicitly mentions 2 methods, propensity scores (PSs) (Rosenbaum & Rubin 1983) and instrumental variables (IVs) (Angrist et al. 1996; Hernan & Robins 2006), as being “relatively well-developed and increasingly used in CER.” The Methodology Report provides the following standards for the use of these 2 methods:

• CI-5. Report the assumptions underlying the construction of propensity scores and the comparability of the resulting groups in terms of the balance of covariates and overlap.
• CI-6. Assess the validity of the instrumental variable (ie, how the assumptions are met) and report the balance of covariates in the groups created by the instrumental variable.

The PS in CI-5 is defined as the conditional probability of receiving the treatment given covariates (Rosenbaum & Rubin 1983). For IV analysis (CI-6), the conditional probabilities of the instrument given covariates (called IV PSs) can also be used to achieve covariate balance between instrument groups (Tan 2006b, 2010b).

Underlying the preceding standards, however, are at least 2 nontrivial statistical tasks: first, estimating PSs (or IV PSs), and then using the estimated PSs to conduct inferences about treatment effects. In fact, the exact PSs in an observational study are unknown and need to be estimated from empirical data, depending on a PS model, often specified in the form of logistic regression, where the response variable is treatment status and the regressors may include main effects, nonlinear terms, or interactions of the covariates. How such PS models are built and fitted can have a substantial impact on subsequent estimation of treatment effects. In particular, inverse probability weighted (IPW) estimation using PSs may perform poorly even when the PS model appears to be nearly correct (Kang & Schafer 2007). Moreover, statistical uncertainty from model building and fitting for estimating PSs must be properly taken into account when conducting inferences (including CIs and hypothesis tests) about treatment effects.

To mitigate possible misspecification of PS models, researchers can use doubly robust estimation by combining both a PS model and an outcome regression (OR) model, where the response variable is the outcome of interest, and the regressors may include main effects, nonlinear terms, or interactions of the covariates, similar to the PS model. For example, various doubly robust estimators for the average treatment effect (ATE) have been derived in the form of an augmented IPW estimator (Robins et al. 1994; Tan 2006a, 2010a). In theory, a doubly robust estimator remains consistent (ie, unbiased in large samples) if either the PS model or the OR model is correctly specified. In practice, all models used are only approximations to some extent. Hence, the performance of doubly robust methods still depends on how PS and OR models are built and fitted.

Gaps in Existing Methods

Per the preceding discussion, applications of PS and IV methods for estimating treatment effects involve statistical modeling and estimation of PSs and OR functions from empirical data, in addition to the fact that certain structural assumptions are required for the identification of treatment effects at the population level. Moreover, statistical uncertainty from such modeling and estimation must be properly incorporated to obtain valid CIs and hypothesis tests about treatment effects. Below, we outline various limitations in existing methods.

From the perspective of model building for PSs and OR, existing methods suffer broadly from 3 types of limitations:

1. Some methods resort to simple models, for example, on OR PSs with main effects of the covariates only, and hence are susceptible to model misspecification.
2. Some methods employ an iterative process of model building and fitting (McCullagh & Nelder 1989; Rosenbaum & Rubin 1984), depending on ad hoc modeling choices about nonlinear and interaction terms of the covariates, where statistical uncertainty is difficult to quantify in subsequent inferences about treatment effects.
3. Some methods apply off-the-shelf machine learning algorithms, such as random forests or boosted machines (Hastie et al. 2009; McCaffrey et al. 2004), but either ignore statistical inference or invoke strong smoothness assumptions to justify consistent estimation and statistical inference (Chernozhukov et al. 2018).

These limitations are aggravated in high-dimensional settings, with either a large number of covariates (even with main effects only) or a large number of possible nonlinear and interaction terms for even a moderate number of covariates.

To perform model fitting for PSs and outcome regression, existing methods mainly use maximum likelihood or variants as in conventional predictive modeling, regardless of the roles of the fitted functions in the stage of estimating treatment effects. In particular, PS models are often fitted as logistic regression by maximum likelihood (Rosenbaum & Rubin 1984). However, as increasingly realized among practitioners, the purpose of estimating PSs is not primarily to predict treatment status, but to balance risk factors for the outcome across treatment groups to control for measured confounders (Weitzen et al. 2004; Westreich et al. 2011; Wyss et al. 2014). Hence, there is a disconnection between model fitting (ie, parameter estimation) and model evaluation: the parameters of logistic regression are estimated by the predictive criterion of maximum likelihood, but the adequacy of the fitted model is evaluated in terms of covariate balance achieved.

The preceding viewpoint about model fitting is in consonance with recent methodological developments on PSs, including a logistic calibration weighting method (Folsom 1991) and related methods (Kim & Haziza 2014; Vermeulen & Vansteelandt 2015), a calibrated likelihood method (Tan 2010a), and covariate-balancing PSs (Imai & Ratkovic 2014). All of these methods can be seen to address the issue of model fitting in classical, low-dimensional settings, in providing a different way from maximum likelihood to estimate PSs or, equivalently, the IPWs. However, the fundamental issue of model building (ie, how to build a PS model) remains open and challenging, especially in cases with a limited number of observations but a large number of regressors (including nonlinear and interaction terms of the covariates).

Specific Aims

The objective of our research is to develop and evaluate a new set of statistically rigorous, numerically tractable, and practically useful methods using PSs and IVs for estimating treatment effects. We build on a strong track record of working on theory, methods, and applications of PSs and IVs (eg, Tan 2006a b, 2010a b c; Winterstein et al. 2012; Huybrechts et al. 2013; Gerhard et al. 2014), and exploit related ideas in high-dimensional statistical theory and methods, including least absolutely shrinkage and selection operator (Lasso) regularized estimation (Tibshirani 1996) and debiased inferences (Zhang & Zhang 2014; van de Geer et al. 2014; Javanmard & Montanari 2014).

The specific aims of our research are as follows:

• Aim 1. Develop new statistical theory and methods for estimating PSs and for drawing inference about treatment effects from observational data, with possibly a large number of covariates or regressors.
• Aim 2. Develop new statistical theory and methods using IVs for drawing inferences about treatment effects from observational data, with possibly a large number of covariates or regressors.
• Aim 3. Develop and disseminate user-friendly software, including accessible and transparent documentation, for implementation of the new methods.

Significance

The methodological significance of our research is to develop new methods that will substantially improve on existing methods in removing or alleviating the methodological limitations mentioned in the “Gaps in Existing Methods” section above.

• Our approach allows flexible models, prespecified with a possibly large number of regressors, such as main effects and interactions of the covariates, and employs sparsity-inducing regularized estimation to facilitate model selection in a numerically tractable and statistically principled manner.
• Our approach employs calibrated estimation, where the loss functions used in fitting PS and OR models are carefully chosen such that valid inferences can be obtained about treatment parameters while accommodating possible model misspecification. For fitting PS models, calibrated estimation directly induces covariate balancing.
• Our approach is potentially applicable to a variety of causal inference tasks, including inferences about treatment effects using PSs under unconfoundedness or using IVs under the corresponding assumptions.

Regarding the applied significance, our research is expected to have substantial and wide-ranging impact on CER by facilitating practical implementation of PS and IV methods for causal inference, as mentioned in the PCORI Methodology Report. The methods and software developed can be widely applied in various CER studies to increase the transparency, validity, and accuracy of the estimates of treatment effects, in reducing the bias and variation associated with ad hoc model building or over-optimistic machine learning. As such, the proposed methods are critical in providing patients, clinicians, and other health care stakeholders with the information needed to make better-informed health care decisions and improve health care delivery and outcomes.

Changes in the Research Strategy

There were 2 main changes in the research strategy. First, the statistical approach described in the original research proposal involves combining boosting or additive learning (Schapire & Freund 2012; Buhlmann & Hothorn 2007; Friedman et al. 2009) and calibrated estimation (Folsom 1991; Tan 2010a). This approach deals with mainly estimation of PSs and IV PSs, but the issue of conducting valid inferences after PS estimation is not addressed.

During the research, we decided to pursue the approach of combining Lasso regularized estimation and calibrated estimation, which not only facilitates model selection in a numerically tractable manner when building and fitting PS and regression models but also makes it possible to derive model-assisted inferences (including CIs and hypothesis tests) about treatment effects in high-dimensional settings. In other words, we believe that the current methods developed in our research are statistically and numerically more satisfactory than those originally suggested in the application for funding.

Second, the original research proposal included empirical applications, using PSs or IVs, to the then-ongoing PCORI project, “Comparative Effectiveness of Adaptive Pharmacotherapy Strategies for Schizophrenia.” That study, which is now completed (Stroup et al. 2019), compared the effectiveness of 4 alternative medication strategies beyond antipsychotic monotherapy (addition of a second antipsychotic, an antidepressant, a mood stabilizer, or a benzodiazepine) for people diagnosed with schizophrenia. Health outcomes included times to psychiatric hospitalization, cardiovascular events, and death. However, the current methods developed in our research assume that the treatment is dichotomous and the outcomes of interest are fully observed (hence excluding possibly censored survival data). Although effort has been made to investigate survival analysis in our new approach (Tan 2019) as part of the modified contract, a fully developed method is beyond the scope of the current project. Hence, we conducted empirical applications using 2 existing observational studies, with Connors et al (1996) using PSs and Card (1995) using IVs; these are among the example studies often used in the related literature to evaluate new methods for causal inference.

Setup

Suppose that a simple random sample of n patients is available from a population under study. The observed data consist of independent and identically distributed observations {(T_i, Y_i, X_i):I = 1,…,n}, where T_i is a dichotomous treatment (T_i = 1 if treated or T_i = 0 otherwise), Y_i is an outcome variable, and X_i = (X_1i, …,X_pi) is a collection of measured covariates. In the potential outcomes framework for causal inference (Neyman 1923; Rubin 1974), 2 potential outcomes ($Y_i^0$, $Y_i^1$) are defined on each patient to indicate what the response would be under treatment 0 or 1, respectively. By consistency, the observed outcome Y_i is assumed to be either $Y_i^0$ or $Y_i^1$, depending on whether T_i = 0 or T_i = 1. Two causal parameters commonly studied in CER are the ATE, defined as E(Y¹ − Y⁰) = μ¹ − μ⁰, with μ^t = E(Y^t), and the ATE on the treated (ATT), defined as E(Y¹ − Y⁰ |T = 1). For concreteness, we mainly discuss estimation of ATE.

For causal inference, a fundamental difficulty in estimating ATEs is that for each patient, only 1 potential outcome, Y⁰ or Y¹, is revealed as the observed outcome, and the other is missing (Holland 1986). Nevertheless, the ATE can be identified from the observed data under the 2 assumptions:

1. Unconfoundedness: T and (Y⁰, Y¹) are conditionally independent given X (Rubin 1976)
2. Overlap (or positivity): 0 < π^*(X) < 1, where π^*(X) = P(T = 1|X) is called the PS (Rosenbaum & Rubin 1983)

Under these assumptions, there are broadly 2 modeling approaches for estimating ATE.

One approach is to build a regression model for the OR function m^t* (X) = E (Y|T = t, X):

where ψ(∙) is an inverse link function; g⁰ (X) and g¹ (X) are vectors of known functions of X, allowed to include nonlinear and interaction terms; and (α⁰, α¹) are vectors of unknown parameters. Let ($\left({\widehat{\alpha }}_{ML}^0,{\widehat{\alpha }}_{ML}^1\right)$) be the least-squares or maximum quasi-likelihood estimators of (α⁰, α¹), and let ${\widehat{m}}_{ML}^t\left(X\right)=m^t(X;{\widehat{\alpha }}_{ML}^t)$. If OR model (1) is correctly specified, then μ^t = E(Y^t) can be consistently estimated by a substitution estimator, ${\widehat{\mu }}_{OR}^t=n^{-1}{\displaystyle {\sum }_{i=1}^n{\widehat{m}}_{ML}^t\left(X_i\right)\hspace{0.17em}}$ for t = 0, 1.

An alternative approach is to build a regression model for the PS, π^*(X) = P(T = 1|X) (Rosenbaum & Rubin 1983):

where Π(∙) is an inverse link function; f(X) is a vector of known functions of X, allowed to include nonlinear and interaction terms; and γ is a vector of unknown parameters. Let ${\widehat{\gamma }}_{ML}$ be the maximum likelihood estimator of γ. Various methods have been proposed to estimate ATEs by matching, stratification, or weighting on the fitted PS ${\widehat{\pi }}_{ML}\left(X\right)=\pi \left(X;{\widehat{\gamma }}_{ML}\right)$. We focus on IPW, which is central to rigorous theory of statistical estimation for causal inference and missing-data problems (eg, Tsiatis 2006). Two standard IPW estimators for μ¹ are

Similarly, 2 standard IPW estimators for μ⁰ are

If PS model (2) is correctly specified, then these IPW estimators are asymptotically valid as the sample size n increases to ∞. However, if PS model (2) is misspecified, the IPW estimators can perform poorly. Even when the PS model is correct or mildly misspecified, inverse weighting can be instable due to estimated PSs close to 0 or 1 (eg, Kang & Schafer 2007).

To attain consistency for the mean μ^t, the OR estimator ${\widehat{\mu }}_{OR}^t$ or the IPW estimator ${\widehat{\mu }}_{rIPW}^t({\widehat{\pi }}_{ML})$ relies on correct specification of OR model (1) or PS model (2), respectively. In contrast, there are doubly robust estimators depending on both OR and PS models in the augmented IPW form, which for μ¹ reads as follows:

where ${\widehat{m}}^1(X)$ and $\widehat{\hspace{0.17em}\pi }(X)$ are fitted values of m^1*(X) and π^*(X), respectively, and

For example, the augmented IPW estimator used by Robins et al (1994) is ${\widehat{\mu }}^1({\widehat{m}}_{ML}^1,\hspace{0.17em}{\widehat{\pi }}_{ML})$, using the fitted values ${\widehat{m}}_{ML}^1\left(X\right)$ and ${\widehat{\pi }}_{ML}(X)$ from maximum (quasi-)likelihood estimation. See Kang & Schafer (2007) and Tan (2007, 2010a) for reviews in low-dimensional settings.

In high-dimensional settings, Belloni et al (2014) and Farrell (2015) studied the estimator ${\widehat{\mu }}^1({\widehat{m}}_{RML}^1,\hspace{0.17em}{\widehat{\pi }}_{RML})$ using the fitted values ${\widehat{m}}_{RML}^1(X)$ and ${\widehat{\pi }}_{RML}(X)$ from Lasso regularized maximum likelihood (RML) estimation or similar methods. Their results are mainly of 2 types, each under suitable conditions. The first type shows double robustness: ${\widehat{\mu }}^1({\widehat{m}}_{RML}^1,\hspace{0.17em}{\widehat{\pi }}_{RML})$ remains consistent if either OR model (1) or PS model (2) is correctly specified. The second type establishes valid CIs: ${\widehat{\mu }}^1({\widehat{m}}_{RML}^1,\hspace{0.17em}{\widehat{\mu }}_{RML})$ admits the asymptotic expansion

if both OR model (1) and PS model (2) are correctly specified. To address the gap between these results, it is desirable to develop new methods that lead to valid CIs without requiring both OR and PS models to be correctly specified.

There has been a large, growing literature on causal inference related to our work, in addition to work from Belloni et al (2014) and Farrell (2015). Examples include Robins et al (2017), Chernozhukov et al (2018), Smucler et al (2019), van der Laan et al (2019), Ning et al (2020), Dukes and Vansteelandt (2021), Hirschberg and Wager (2021), and Bradic et al (2021). See Tan (2020a, 2020b) and Ghosh and Tan (2020) for further discussion.

PSs and IPW Estimation

Note: The material in this section is adapted from Tan Z. Regularized calibrated estimation of propensity scores with model misspecification and high-dimensional data. Biometrika. 2020;107;137-158.

We developed regularized calibrated estimation for fitting PS models in high-dimensional settings. For concreteness, assume that PS model (2) is a logistic regression model:

where f(x) = {1, f₁ (x), …, f_p (x)}^T is a vector of known functions including the constant, and γ = (γ₀, γ₁, …, γ_p)^T is a vector of unknown coefficients. The maximum likelihood estimator ${\widehat{\gamma }}_{ML}$ is defined as a minimizer of the likelihood loss

or, equivalently, as a solution to the score equation

Alternatively, calibrated estimation uses a system of estimating equations such that the weighted averages of the covariates in the treated subsample are equal to the simple averages in the overall sample. The calibrated estimator ${\widehat{\gamma }}_{CAL}^1$ is defined as a solution to

Interestingly, the estimator ${\widehat{\gamma }}_{CAL}^1$ can be equivalently defined as a minimizer of the following loss function, called the calibration loss:

Moreover, ℓ_CAL (γ) is convex in γ and is strictly convex and bounded from below under a certain nonseparation condition (Tan 2020a). The preceding idea of calibrated estimation and similar methods have been studied, and sometimes independently (re)derived, in various contexts of causal inference, missing-data problems, and survey sampling (eg, Folsom 1991; Tan 2010a; Graham et al. 2012; Hainmueller 2012; Imai & Ratkovic 2014; Kim & Haziza 2014; Vermeulen & Vansteelandt 2015; Chan et al. 2016).

To compare calibrated and likelihood estimation, we establish a new relationship between the loss functions ℓ_CAL (γ) and ℓ_ML (γ). To allow for misspecification of PS model (8), we write ℓ_ML (γ) = κ_ML (γ^Tf) and ℓ_CAL (γ) = κ_CAL (γ^Tf), where for a function g(x),

Then, κ_ML (g^*) and κ_CAL (g^*) are well defined for the true log odds ratio g^*(x) = log[π^*(x)/{1-π* (x)}], even when PS model (8) is misspecified, that is, when g^*(x) is not of the form γ^T f(x). For 2 probabilities, ρ ∈ (0, 1) and ρ′∈ (0,1), the Kullback-Leibler divergence is

In addition, let K(ρ,ρ′) = ρ′/ρ – 1 − log(ρ′/ρ), which is strictly convex in ρ′/ρ, with a minimum of 0 at ρ′/ρ = 1. Then Tan (2020a, Proposition 1) shows that

Hence, minimizing the expected calibration loss in γ involves reducing both the expected likelihood loss E [L{π(X;γ), π^*(X)}] and an additional term E [K{π(X; γ), π^*(X)}], which can be related to the mean squared relative error (MRSE) between π(X; γ) and π^*(X):

This measure of relative errors in PSs governs the mean squared error (MSE) of the IPW estimator, $MSE\left(\gamma \right)=E\hspace{0.17em}\left[{\{{\widehat{\mu }}_{IPW}^1\left(\gamma \right)-{\mu }^1\}}^2\right]$, where ${\widehat{\mu }}_{IPW}^1\left(\gamma \right)=n^{-1}{\displaystyle {\sum }_{i=1}^n{T}_i{Y}_i/\hspace{0.17em}\pi (X_{i;}\gamma )}$. The following result was obtained by Tan (2020a, Proposition 2):

Suppose that E[(Y¹)²] ≤ c and π^*(X) ≥ δ for some constants c > 0 and δ ∈ (0,1).

• If π(X; γ) ≥ a π^*(X) for some constant a ∈ (0,1/2), then
  $$\text{MRSE}(\gamma )\hspace{0.17em}\le \frac{5}{3a}E\left[K\left\{\pi \left(X;\gamma \right),{\pi }^{\text{*}}\left(X\right)\right\}\right]\le \frac{5}{3a}E\left\{{\ell }_{CAL}\left(\gamma \right)-{\kappa }_{CAL}\left(g^{\text{*}}\right)\right\}.$$

  (11)
  The factor 5/(3a) in general cannot be improved up to a constant, independent of a.
• If π(X; γ) ≥ b for some constant b ∈ (0,1), then
  $$\text{MRSE}(\gamma )\hspace{0.17em}\le \frac{1}{2b^2}E\left[L\left\{\pi \left(X;\gamma \right),{\pi }^{\text{*}}\left(X\right)\right\}\right]\le \frac{1}{2b^2}E\left\{{\ell }_{ML}\left(\gamma \right)-{\kappa }_{ML}\left(g^{\text{*}}\right)\right\}.$$

  (12)
  The factor 1/(2b²) in general cannot be improved up to a divisor of order log(b⁻¹).

The expected calibration loss is inflated in (11) by a leading factor 5/(3a), which remains bounded from above as long as a is bounded away from 0, even when some PSs π(X; γ) are close to 0. In contrast, the expected likelihood loss is inflated in (12) by a leading factor 1/(2b²), which may be arbitrarily large when some PSs π(X; γ) are close to 0. This result demonstrates that the minimization of the expected calibration loss controls the MRSE of PSs more strongly than does the expected likelihood loss. See Tan (2020a), Section 3.2, for further discussion.

The preceding discussion compares calibrated and likelihood estimation through their loss functions. Next, we propose a regularized calibrated estimator of γ in PS model (8). There are 2 motivations. First, regularization is needed in the situation where ${\widehat{\gamma }}_{CAL}^1$ does not exist because the calibration loss (10) is unbounded from below. In our numerical study (Tan 2020a), nonconvergence was found, for example, in 100% of repeated simulations with n = 200 and p = 50. Second, regularization is also needed in high-dimensional settings where the dimension of the regressor vector f (X) is close to or greater than the sample size.

The regularized calibrated estimator, denoted by ${\widehat{\gamma }}_{RCAL}^1$, is defined by minimizing the calibration loss ℓ_CAL (γ) with a Lasso penalty (Tibshirani 1996):

where γ_1:p = (γ₁, …, γ_p)^T excluding the intercept γ₀, ∥⋅∥₁ denotes the L₁ norm such that $\parallel {\gamma }_{1:p}{\parallel }_1={\displaystyle {\sum }_{j=1}^p\left|{\gamma }_j\right|}$, and λ ≥ 0 is a tuning parameter. By the Karush-Kuhn-Tucker (KKT) condition for minimization of (13), the fitted PSs, ${\widehat{\pi }}_{RCAL}^1\left(X\right)=\pi \left(X;{\widehat{\gamma }}_{RCAL}^1\right)$, satisfies

The equality holds in (15) for any j such that the jth estimate ${({\widehat{\gamma }}_{RCAL}^1)}_j$ is nonzero. The weighted average of each covariate f_j(X_i) over the treated group may differ from the overall average of f_j(X_i) by no more than λ. In other words, introducing the Lasso penalty to calibrated estimation yields a relaxation of the equalities (9) to box constraints (15). Moreover, by (14), the IPWs, $1/\hspace{0.17em}{\widehat{\pi }}_{RCAL}^1\left(X_i\right)$ with T_i = 1, sum to the sample size n. Then, the 2 resulting IPW estimators, ${\widehat{\mu }}_{IPW}^1({\widehat{\pi }}_{RCAL}^1)$ and ${\widehat{\mu }}_{rIPW}^1({\widehat{\pi }}_{RCAL}^1)$, obtained from (3) with ${\widehat{\pi }}_{ML}(X)$ replaced by ${\widehat{\pi }}_{RCAL}^1(X)$, are identical to each other.

We present a novel algorithm for computing the proposed estimator ${\widehat{\gamma }}_{RCAL}^1$, that is, minimizing ℓ_RCAL (γ) in (13) for any fixed λ. The basic idea of the algorithm is to iteratively form a quadratic approximation to the calibration loss ℓ_CAL (γ) in (10) and solve a Lasso-penalized weighted least-squares problem, similar to existing algorithms for Lasso-penalized maximum likelihood–based logistic regression (eg, Friedman et al. 2010). However, we construct a suitable quadratic approximation after an additional step of replacing certain sample quantities with model expectations. This idea is known as Fisher scoring and was previously used to derive the iterative reweighted least-squares method for fitting generalized linear models with noncanonical links, such as probit regression (McCullagh & Nelder 1989). Moreover, to reduce computational cost, we exploit the majorization-minimization technique (Wu & Lange 2010; Bohning & Lindsay 1988) in solving the Lasso-penalized weighted least-squares problem. See Tan (2020a), Section 4.2, for a detailed discussion.

We provide a theoretical analysis of the regularized calibrated estimator ${\widehat{\gamma }}_{RCAL}^1$ and the resulting IPW estimator of μ¹, allowing for misspecification of PS model (8), in high-dimensional settings where the number of regressors in f(X) is close to or greater than the sample size. Our analysis requires a sparsity assumption that only a small but unknown subset (relative to the sample size) of “significant” regressors is associated with nonzero coefficients. Such assumptions are commonly made in high-dimensional regression, where Lasso-type-penalized estimation can theoretically be shown to achieve nearly as small errors of estimation as in oracle estimation, where those significant regressors were known and only those were included in the regression model (eg, Buhlmann & van de Geer 2011). Heuristically, sparsity-based methods and theory, including ours, can be motivated according to the “bet on sparsity” principle: “Use a procedure that does well in sparse problems, since no procedure does well in dense problems” (Hastie et al. 2015).

Our analysis of Lasso-penalized M-estimators deals with the convergence of the estimators, ${\widehat{\gamma }}_{RCAL}^1$, to their target values ${\overline{\gamma }}_{CAL}^1$ under model misspecification, which is related to but distinct from previous results on excess prediction errors (eg, Buhlmann & van de Geer 2011). Moreover, our analysis of ${\widehat{\mu }}_{IPW}^1\left({\widehat{\pi }}_{RCAL}^1\right)$, the IPW estimator of μ¹ based on ${\widehat{\gamma }}_{RCAL}^1$, carefully exploits the results described earlier on the calibration loss ℓ_CAL (γ) to obtain convergence under weaker conditions than previously realized. Our results show that the squared difference $|{\widehat{\mu }}_{IPW}^1\left({\widehat{\pi }}_{RCAL}^1\right)-{\widehat{\mu }}_{IPW}^1\left({\overline{\pi }}_{CAL}^1\right){|}^2$ is of the order

provided that $s\sqrt{\frac{\log p}{n}}$ is sufficiently small, for example, tends to 0 as n and p increase, where s denotes the sparsity size of ${\overline{\gamma }}_{CAL}^1$, that is, the number of nonzero elements in ${\overline{\gamma }}_{CAL}^1$. See Tan (2020a, Section 4.3) for a detailed discussion.

So far, our theory and methods have focused mainly on the estimation of μ¹, but they can be directly extended to an estimation of μ⁰, and hence the ATE, μ¹ − μ⁰. For estimation of μ⁰ with PS model (8), the calibrated estimator of γ, denoted by ${\widehat{\gamma }}_{CAL}^0$, is defined as a solution to

By exchanging T with 1 − T and exchanging γ with − γ in (10), the corresponding loss function minimized by ${\widehat{\gamma }}_{CAL}^0$ is

For fixed λ ≥ 0, the regularized calibrated estimator ${\widehat{\gamma }}_{RCAL}^0$ is defined as a minimizer of

The fitted PS, ${\widehat{\pi }}_{RCAL}^0\left(X\right)=\pi \left(X;{\widehat{\gamma }}_{RCAL}^0\right)$, then satisfies (14) and (15) with T_i replaced by 1 − T_i and ${\widehat{\pi }}_{RCAL}^1\left(X_i\right)$ replaced by $1-{\widehat{\pi }}_{RCAL}^0\left(X_i\right)$, and similar interpretations apply as discussed previously. The resulting IPW estimator of μ⁰ is ${\widehat{\mu }}_{IPW}^0({\widehat{\pi }}_{RCAL}^0)$, obtained from (4) with ${\widehat{\pi }}_{ML}(X)$ replaced by ${\widehat{\pi }}_{RCAL}^0(X)$, and that of μ¹ − μ⁰ is ${\widehat{\mu }}_{IPW}^1\left({\widehat{\pi }}_{RCAL}^1\right)-{\widehat{\mu }}_{IPW}^0({\widehat{\pi }}_{RCAL}^0)$.

An interesting feature of our approach is that 2 different estimators of the PS are used when estimating μ¹ and μ⁰. The estimators ${\widehat{\gamma }}_{RCAL}^1$ and ${\widehat{\gamma }}_{RCAL}^0$ may in general have different asymptotic limits when PS model (8) is misspecified, even though their asymptotic limits coincide when model (8) is correctly specified. Such possible differences should not be of concern: the 2 estimators ${\widehat{\mu }}_{IPW}^1(\widehat{\pi })$ and ${\widehat{\mu }}_{IPW}^0(\widehat{\pi })$ are decoupled, involving 2 disjoint subsets of fitted PSs on the treated {i: T_i = 1} and the untreated {i: T_i = 0}, respectively. In fact, separate estimation of PSs and inverse probability weights for the treated and untreated samples can lead to more flexible approximations, and hence potentially less bias in the presence of model misspecification. Furthermore, the discovery of whether substantial differences exist between these separately fitted PSs can be used for diagnosis of the validity of model (8). See Chan et al (2016, Section 2.3) for a related discussion.

The proposed method for estimating PSs is implemented as part of the publicly released R package RCAL (Tan & Sun 2020). The implementation is based on the Fisher scoring descent algorithm in Tan (2020a). At the low level, the least-squares Lasso problem is solved using a variation of the active set algorithm, which enjoys a finite termination property (Osborne et al. 2000; Yang & Tan 2018).

Model-Assisted and Doubly Robust Inferences

Note: The material in this section is adapted from Tan Z. Model-assisted inference for treatment effects using regularized calibrated estimation with high-dimensional data. Ann Statist. 2020;48:811-837.

We develop new methods and theory in high-dimensional settings to obtain not only doubly robust point estimators for ATEs, which remain consistent if either a PS model or an OR model is correctly specified, but also model-assisted CIs, which are valid when the PS model is correctly specified but the OR model may be misspecified.

The development in the “PSs and IPW Estimation” section only involves estimating PSs and applying the IPW estimators. To mitigate possible misspecification of PS models and facilitate constructing CIs, we use a doubly robust estimator of μ¹ in the augmented IPW form (5), depending on both a PS model and an OR model. There are 2 motivations, in theory. First, the double-robustness property ensures consistency of point estimation even when 1 of the OR and PS models is misspecified. Second, more importantly in high-dimensional settings, augmented IPW estimation makes it possible to establish a simple asymptotic expansion of the resulting estimator such that valid and numerically tractable CIs can be obtained for μ¹. In contrast, obtaining valid CIs is difficult based on IPW estimation only, due to penalized estimation in high-dimensional settings, even when the PS model is correctly specified.

To focus on main ideas, consider a logistic PS model (8) as in the “PSs and IPW Estimation” section and a linear OR model,

where ${\alpha }^1={({\alpha }_0^1,\hspace{0.17em}{\alpha }_1^1,\hspace{0.17em}\dots ,\hspace{0.17em}{\alpha }_p^1)}^T$, of the same dimension as γ. That is, OR model (1) is specified with the identity link and the regressor vector g¹ (X) is the same as f (X) in PS model (2). This condition can be satisfied possibly after enlarging model (1) or (2) to reach the same dimension. Our point estimator of μ¹ is

where φ () is defined in (6), ${\widehat{\pi }}_{RCAL}^1\left(X\right)=\pi \left(X;{\widehat{\gamma }}_{RCAL}^1\right)$ with the regularized calibrated estimator ${\widehat{\gamma }}_{RCAL}^1$ as in the “PSs and IPW Estimation” section, and ${\widehat{m}}_{RWL}^1\left(X\right)=m^1(X;{\widehat{\alpha }}_{RWL}^1)$ with ${\widehat{\alpha }}_{RWL}^1$ defined as follows. See Tan (2020b), Section 3.2, for a discussion on how the construction of these estimators is linked to the desired asymptotic expansion (22) for ${\widehat{\mu }}^1\left({\widehat{m}}_{RWL}^1,{\widehat{\pi }}_{RCAL}^1\right)$.

The estimator ${\widehat{\alpha }}_{RWL}^1$ is a regularized weighted least-squares estimator of α¹, defined as a minimizer to the penalized loss function

where ${\alpha }_{1:p}^1={({\alpha }_1^1,\hspace{0.17em}\dots ,\hspace{0.17em}{\alpha }_p^1)}^T$ excluding the intercept ${\alpha }_0^1$, λ ≥ 0 is a tuning parameter, and ${\ell }_{WL}\left({\alpha }^1;\hspace{0.17em}{\widehat{\gamma }}_{RCAL}^1\right)$ is the weighted least-squares loss,

That is, the observations in the treated group are weighted by $\left\{1-{\widehat{\pi }}_{RCAL}^1\left(X_i\right)\right\}/{\widehat{\pi }}_{RCAL}^1(X_i)$, which differs slightly from the commonly used inverse probability weight $1/{\widehat{\pi }}_{RCAL}^1(X_i)$. Larger weights are associated with observations with lower-fitted PSs, thereby requiring model (16) to be better fitted in covariate regions with more missing data.

Similarly as discussed regarding ${\widehat{\gamma }}_{RCAL}^1$, there are simple and interesting implications of the estimator ${\widehat{\alpha }}_{RWL}^1$. By the KKT condition for minimizing (18), the fitted OR function ${\widehat{m}}_{RWL}^1\left(X\right)$ satisfies

In (20), the inequality reduces to equality for any j such that the jth estimate ${({\widehat{\alpha }}_{RWL}^1)}_j$ is nonzero. Equation (19) implies that the estimator ${\widehat{\mu }}^1\left({\widehat{m}}_{RWL}^1,\hspace{0.17em}{\widehat{\pi }}_{RCAL}^1\right)$ can be recast as

which takes the form of linear prediction estimators known in the survey literature (eg, Sarndal et al. 1992): $n^{-1}{\displaystyle {\sum }_{i=1}^n\hspace{0.17em}{T}_i{Y}_i+\left(1-T_i\right)\hspace{0.17em}{\widehat{m}}^1\left(X_i\right)}$, for some fitted OR function ${\widehat{m}}^1\left(X_i\right)$. As a consequence, ${\widehat{\mu }}^1\left({\widehat{m}}_{RWL}^1,{\widehat{\pi }}_{RCAL}^1\right)$ always falls within the range of the observed outcomes {Y_i: T_i = 1, i = 1,…,n} and the predicted values $\{{\widehat{m}}_{RWL}^1\left(X_i\right):\hspace{0.17em}{T}_i=1,\hspace{0.17em}i=1,\hspace{0.17em}\dots ,\hspace{0.17em}n\}$. This boundedness property is not satisfied by the estimator ${\widehat{\mu }}^1\left({\widehat{m}}_{RML}^1,\hspace{0.17em}{\widehat{\pi }}_{RML}\right)$, where ${\widehat{m}}_{RML}^1(X)$ and ${\widehat{\pi }}_{RML}(X)$ are the fitted values similar to ${\widehat{m}}_{ML}^1(X)$ and ${\widehat{\pi }}_{ML}(X)$ in the “Setup” section but based on regularized likelihood estimation in OR and PS models.

We provide a high-dimensional analysis of the proposed estimator ${\widehat{\mu }}^1\left({\widehat{m}}_{RWL}^1,\hspace{0.17em}{\widehat{\pi }}_{RCAL}^1\right)$, allowing for possible model misspecification. See Tan (2020b, Section 3.3) for a detailed discussion. In contrast, with asymptotic expansion (7) for ${\widehat{\mu }}^1\left({\widehat{m}}_{RML}^1,\hspace{0.17em}{\widehat{\pi }}_{RML}\right)$ with correctly specified OR and PS models, our main result shows that under suitable conditions, the estimator ${\widehat{\mu }}^1\left({\widehat{m}}_{RWL}^1,\hspace{0.17em}{\widehat{\pi }}_{RCAL}^1\right)$ admits the asymptotic expansion

where ${\overline{\pi }}_{CAL}^1\left(X\right)=\pi \left(X;{\overline{\gamma }}_{CAL}^1\right)$, ${\overline{m}}_{WL}^1\left(X\right)=m^1(X;{\overline{\alpha }}_{WL}^1)$, and ${\overline{\gamma }}_{CAL}^1$ and ${\overline{\alpha }}_{WL}^1$ are limit or target values defined as follows. With possible model misspecification, the target value ${\overline{\gamma }}_{CAL}^1$ is defined as a minimizer of the expected calibration loss

If PS model (8) is correctly specified, then ${\overline{\pi }}_{CAL}^1\left(X\right)={\pi }^{*}(X)$. Otherwise, ${\overline{\pi }}_{CAL}^1\left(X\right)$ may differ from π^*(X). The target value ${\overline{\alpha }}_{WL}^1$ is defined as a minimizer of the expected loss

If OR model (16) is correctly specified, then ${\overline{m}}_{WL}^1\left(X\right)=m^{1*}(X)$. But ${\overline{m}}_{WL}^1\left(X\right)$ may in general differ from m^1*(X). Then, Tan (2020b), Proposition 1, shows that if either logistic PS model (8) or linear OR model (16) is correctly specified, the following results hold under suitable conditions related to the sparsity of the target values ${\overline{\gamma }}_{CAL}^1$ and ${\overline{\alpha }}_{WL}^1$.

• ${\widehat{\mu }}^1\left({\widehat{m}}_{RWL}^1,\hspace{0.17em}{\widehat{\pi }}_{RCAL}^1\right)$ is consistent for μ¹ and asymptotically normally distributed:
  $$\sqrt{n}\left\{{\widehat{\mu }}^1\left({\widehat{m}}_{RWL}^1,{\widehat{\pi }}_{RCAL}^1\right)-{\mu }^1\right\}{\to }_{D}N\left(0,V\right),$$

  (23)
  where $V=\hspace{0.17em}var\hspace{0.17em}\left\{\phi \left(Y,T,X;\hspace{0.17em}{\overline{m}}_{WL}^1,\hspace{0.17em}{\overline{\pi }}_{CAL}^1\right)\right\}$.
• A consistent estimator of V is
  $$\widehat{V}=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left\{\phi \left(Y_i,T_i,X_i;{\widehat{m}}_{RWL}^1,{\widehat{\pi }}_{RCAL}^1\right)-{\widehat{\mu }}^1\left({\widehat{m}}_{RWL}^1,{\widehat{\pi }}_{RCAL}^1\right)\right\}}^2.$$

  (24)
• An asymptotic (1 − c) CI for μ¹ is
  $${\widehat{\mu }}^1\left({\widehat{m}}_{RWL}^1,{\widehat{\pi }}_{RCAL}^1\right)\pm z_{\frac{c}{2}}\sqrt{\frac{\widehat{V}}{n}},$$

  (25)
  where z_c/2 is the (1 − c/2) quantile of N(0,1). Hence, a doubly robust CI for μ¹ is obtained.

Next, we present model-assisted CIs for μ¹, in the setting where a generalized linear OR model is used together with a logistic PS model. Consider a generalized linear OR model with a canonical link,

that is, model (1) with the vector of covariate functions g¹ (X) taken to be the same as f (X) in model (8). Our point estimator of μ¹ is ${\widehat{\mu }}^1\left({\widehat{m}}_{RWL}^1,\hspace{0.17em}{\widehat{\pi }}_{RCAL}^1\right)$ as defined in (17), where ${\widehat{\pi }}_{RCAL}^1\left(X\right)=\hspace{0.17em}\pi \left(X;{\widehat{\gamma }}_{RCAL}^1\right)$ and ${\widehat{m}}_{RWL}^1\left(X\right)=m^1\left(X;\hspace{0.17em}{\widehat{\alpha }}_{RWL}^1\right)$. The estimator ${\widehat{\gamma }}_{RCAL}^1$ is the regularized calibrated estimator of γ as before. However, ${\widehat{\alpha }}_{RWL}^1$ is a regularized weighted likelihood estimator of α¹, defined as a minimizer of the penalized loss function

where ${\alpha }_{1:p}^1={({\alpha }_1^1,\hspace{0.17em}\dots ,\hspace{0.17em}{\alpha }_p^1)}^T$, excluding the intercept ${\alpha }_0^1$, λ ≥ 0 is a tuning parameter, and ${\ell }_{WL}\left({\alpha }^1;\hspace{0.17em}{\widehat{\gamma }}_{RCAL}^1\right)$ is the weighted likelihood loss,

with $\text{Ψ}\left(u\right)={\displaystyle {\int }_0^u\psi \left(\text{t}\right)\hspace{0.17em}dt}$. The regularized weighted least-squares estimator ${\widehat{\alpha }}_{RWL}^1$ for a linear OR model is recovered in the special case of the identity link, ψ(u) = u and Ψ(u) = u²/2. In addition, the KKT condition for minimizing the penalized loss function (27) remains the same as in (19) and (20); hence, the estimator ${\widehat{\mu }}^1\left({\widehat{m}}_{RWL}^1,\hspace{0.17em}{\widehat{\pi }}_{RCAL}^1\right)$ can be put in the prediction form (21), which ensures that the boundedness property that ${\widehat{\mu }}^1\left({\widehat{m}}_{RWL}^1,\hspace{0.17em}{\widehat{\pi }}_{RCAL}^1\right)$ always falls within the range of the observed outcomes {Y_i: T_i = 1, i = 1,…, n} and the predicted values $\{{\widehat{m}}_{RWL}^1\left(X_i\right):\hspace{0.17em}{T}_i=1,\hspace{0.17em}i=1,\hspace{0.17em}\dots ,\hspace{0.17em}n\}$.

Tan (2020b, Proposition 3) shows that if logistic PS model (8) is correctly specified but OR model (26) may be misspecified, then ${\widehat{\mu }}^1\left({\widehat{m}}_{RWL}^1,\hspace{0.17em}{\widehat{\pi }}_{RCAL}^1\right)$ admits asymptotic expansion (22) with ${\overline{\pi }}_{CAL}^1\left(X\right)={\pi }^{*}(X)$, and results (23) to (25) hold under suitable conditions, where ${\overline{m}}_{WL}^1\left(X\right)=m^1(X;\hspace{0.17em}{\overline{\alpha }}_{WL}^1)$ and the target value ${\overline{\alpha }}_{WL}^1$ is defined as a minimizer of the expected loss $E\left\{{\ell }_{WL}\left({\alpha }^1;\hspace{0.17em}{\overline{\gamma }}_{CAL}^1\right)\right\}$, similar to before. The CIs obtained for μ¹ are said to be PS based, OR assisted. We make the following remarks:

• As mentioned in Tan (2020b, Section 3.4), it is also possible to derive OR-based, PS-assisted CIs for μ¹ with a nonlinear OR model. In that case, an estimator of γ would be defined depending on an estimator of α¹, similar to how ${\widehat{\alpha }}_{RWL}^1$ is defined depending on ${\widehat{\gamma }}_{RCAL}^1$.
• The preceding PS-based, OR-assisted CIs are convenient when ATEs associated with several outcomes need to be estimated using the same set of fitted PSs. Moreover, regularized calibration estimation of PSs enjoys attractive properties, as discussed in the “PSs and IPW Estimation” section.
• With additional complexity, a 2-step procedure was developed by Ghosh and Tan (2020) to obtain doubly robust CIs for μ¹.

Our theory and methods are presented mainly on estimation of μ¹, but they can be directly extended for estimating μ⁰ and hence ATE, that is, μ¹ − μ⁰. Consider a logistic PS model (8) and a generalized linear OR model,

where f(X) is the same vector of covariate functions as in PS model (8) and α⁰ is a vector of unknown parameters. Our point estimator of ATE is ${\widehat{\mu }}^1\left({\widehat{m}}_{RWL}^1,\hspace{0.17em}{\widehat{\pi }}_{RCAL}^1\right)-{\widehat{\mu }}^0\left({\widehat{m}}_{RWL}^0,\hspace{0.17em}{\widehat{\pi }}_{RCAL}^0\right)$, and that of μ⁰ is

where φ() is defined in (6), ${\widehat{\pi }}_{RCAL}^0\left(X\right)=\pi \left(X;{\widehat{\gamma }}_{RCAL}^0\right)$, and ${\widehat{m}}_{RWL}^0\left(X\right)=m^0(X;{\widehat{\alpha }}_{RWL}^0)$. The estimator ${\widehat{\gamma }}_{RCAL}^0$ is the same regularized calibrated estimator as in the “PSs and IPW Estimation” section. The estimator ${\widehat{\alpha }}_{RWL}^0$ is defined similarly as ${\widehat{\alpha }}_{RWL}^1$, but with the loss function ${\ell }_{WL}\left({\alpha }^1;{\widehat{\gamma }}_{RCAL}^1\right)$ in (27) replaced by

Under similar conditions as in the estimation of μ¹, the estimator ${\widehat{\mu }}^0\left({\widehat{m}}_{RWL}^0,\hspace{0.17em}{\widehat{\pi }}_{RCAL}^0\right)$ admits the asymptotic expansion

where ${\overline{\pi }}_{CAL}^0\left(X\right)=\hspace{0.17em}\pi \left(X;{\overline{\gamma }}_{CAL}^0\right)$, ${\overline{m}}_{WL}^0\left(X\right)=m^0(X;\hspace{0.17em}{\overline{\alpha }}_{WL}^0)$, and ${\overline{\gamma }}_{CAL}^0$ and ${\overline{\alpha }}_{WL}^0$ are the target values defined similarly as ${\overline{\gamma }}_{CAL}^1$ and ${\overline{\alpha }}_{WL}^1$. Then, Wald CIs for μ⁰ and ATE can be derived similarly as in the estimation of μ⁰ and shown to be either doubly robust in the case of linear OR models, or valid if PS model (8) is correctly specified but OR models (26) and (28) may be misspecified for nonlinear outcome models.

The proposed method for model-assisted inference about ATEs is implemented as part of the publicly released R package RCAL (Tan & Sun 2020). The package provides both functions for regularized calibrated estimation of PSs and OR functions and for subsequent estimation of ATEs. In addition to a full reference manual, a vignette is also included in the package to give a direct and accessible introduction to the method with simple examples.

To facilitate applications, we focus on the parametric setting as described in the “Setup” section, where OR and PS models are generalized linear regressions in the usual form but prespecified with possibly a large number of regressors. Each regressor is a known function of the covariates, for example, the main effect or a spline term of a single covariate or an interaction of 2 covariates. Our method employs Lasso-type-penalized estimation to perform variable selection (or, more precisely, regressor selection). For technical justification, our theory is developed in Tan (2020a b) with suitable assumptions regarding the exact sparsity (ie, the number of nonzero coefficients) on the limit values of the estimators in OR and PS models, which may be misspecified. Nevertheless, various extensions can be investigated. For example, approximate or weak sparsity can be accommodated, along similar lines as in Negahban et al (2012), Smucler et al (2019), and Bradic et al (2021). Moreover, regularized calibrated estimation can be combined with related methods, such as nonlinear regression models (Hastie et al. 2009), targeted maximum likelihood estimation (van der Laan & Rubin 2006; van der Laan & Rose 2017), and ensemble learning (van der Laan et al. 2007).

Simulation Studies

We conduct extensive simulation studies to evaluate the performances of the proposed methods and compare them with existing methods in both low-dimensional and high-dimensional settings. These studies are reported, in full detail, in Tan (2020a) for PS estimation and in Tan (2020b) for model-assisted inference about ATEs. Below, we describe the simulation study in Tan (2020b, Section 4).

Let X = (X₁, …, X_p) be multivariate normal with means 0 and covariances cov (X_j, X_k) = 2^−|j−k| for 1 ≤ j, k ≤ p. In addition, let $X_j^{†}=\hspace{0.17em}{X}_j+{\left\{max\left(0,\hspace{0.17em}{X}_j+1\right)\right\}}^2$ for j = 1,…, 4. Consider the following data-generating configurations:

C1.

Generate T given X from a Bernoulli distribution with

$$P\left(T=1\text{|}X\right)={\{1+\exp (-1-X_1-0.5X_2-0.25X_3-0.125X_4)\}}^{-1},$$

and, independently, generate Y¹ given X from a normal distribution with variance 1 and mean

$$\text{E}(Y^1\text{|X)}=X_1+0.5X_2+0.25X_3+0.125X_4.$$

C2.

Generate T given X as in (C1), but, independently, generate Y¹ given X from a normal distribution with variance 1 and mean

$$\text{E}(Y^1\text{|X)}=X_1^{†}+0.5X_2^{†}+0.25X_3^{†}+0.125X_4^{†}.$$

C3.

Generate Y¹ given X as in (C1), but, independently, generate T given X from a Bernoulli distribution with

$$P\left(T=1\text{|}X\right)={\{1+\exp (-1-X_1^{†}-0.5X_2^{†}-0.25X_3^{†}-0.125X_4^{†})\}}^{-1}.$$

To study the estimation of μ¹, note that the observed data consist of independent and identically distributed observations {(T_iY_i, Y_i, X_i): i = 1,…, n}. Consider logistic PS model (8) and linear OR model (16), both with f_j(X) = X_j for j = 1,…, p. Then, the 2 models can be classified as follows, depending on the data configuration above: