Double Robust Variance Estimation
Bonnie E. Shook-Sa Paul N. Zivich Chanhwa Lee Keyi Xue Rachael K. Ross Jessie K. Edwards Jeffrey S. A. Stringer Stephen R. Cole

 Interpretable Machine Learning Models for Predicting the Next Targets of Activist Funds
    [
    April 29, 2024
=======================================================================================


Doubly robust estimators have gained popularity in the field of causal inference due to their ability to provide consistent point estimates when either an outcome or exposure model is correctly specified. However, the influence function based variance estimator frequently used with doubly robust estimators is only consistent when both the outcome and exposure models are correctly specified. Here, use of M-estimation and the empirical sandwich variance estimator for doubly robust point and variance estimation is demonstrated. Simulation studies illustrate the properties of the influence function based and empirical sandwich variance estimators. Estimators are applied to data from the Improving Pregnancy Outcomes with Progesterone (IPOP) trial to estimate the effect of maternal anemia on birth weight among women with HIV. In the example, birth weights if all women had anemia were estimated to be lower than birth weights if no women had anemia, though estimates were imprecise. Variance estimates were more stable under varying model specifications for the empirical sandwich variance estimator than the influence function based variance estimator. 
	

	
	


§ INTRODUCTION

Two common approaches for estimating the causal effect of a point exposure or treatment with observational studies are methods based on inverse probability of treatment weighting (IPTW) and the g-formula. IPTW estimators <cit.> rely on estimation of the probability of exposure given a set of covariates, i.e., the propensity score. Alternatively, the g-formula <cit.> is a direct standardization approach where the outcome is modeled as a function of the exposure and covariates. When estimated nonparametrically, IPTW and g-formula estimators are equivalent <cit.>. However, in high dimensional settings (e.g., settings with continuous covariates), analysts typically resort to parametric models for the exposure and outcome, and the IPTW and g-formula estimators are no longer equivalent due to differing parametric constraints. Both approaches provide consistent estimators of the average causal effect (ACE) if the corresponding model is correctly specified. Because the functional forms of these “nuisance” models are typically unknown, correct specification remains an obstacle to valid effect estimation. 

Augmented inverse probability weighted (AIPW) estimators and targeted maximum likelihood estimators (TMLE) combine estimates from the exposure and outcome models in such a way that resulting estimators are “doubly robust.” That is, they are consistent when either the exposure or outcome model is correctly specified, but not necessarily both  <cit.>. In other words, doubly robust methods give the analyst two opportunities for correct model specification. Two common AIPW estimators are an estimator that incorporates estimated propensity scores from the exposure model and pseudo-outcomes from the outcome model (i.e., the classic AIPW) <cit.> and a g-formula type estimator where the outcome model is weighted by the IPTW (i.e., the weighted regression AIPW) <cit.>. TMLE is a variation of the weighted regression AIPW estimator that instead combines IPTW and the outcome model through a “targeting” model fit using the outcome predictions and weights <cit.>.  Both the weighted regression AIPW estimator and TMLE respect the parameter space of the estimand, while the classic AIPW can produce estimates outside of the range of the parameter space <cit.>.

In addition to the double robustness property, doubly robust estimators also offer precision advantages over IPTW estimators. When both nuisance models are correctly specified, the three doubly robust estimators previously discussed are semiparametric efficient, having the smallest asymptotic variance of all IPTW estimators <cit.>. This means that even in settings where the analyst is confident in the specification of one of the models (e.g., the exposure mechanism is well known), doubly robust estimators can offer improved efficiency relative to IPTW estimators, provided that the outcome model is not grossly misspecified. Given these potential precision advantages and that the true specification of both models is typically unknown, doubly robust approaches have gained popularity <cit.>.

Despite adoption of doubly robust estimators, there is a major limitation with the corresponding influence function based variance estimator commonly used in applications and provided by software. Namely, this variance estimator is only consistent if both the propensity and outcome models are correctly specified  <cit.>. This limitation has been highlighted previously <cit.>, yet the influence function based variance estimator is still commonly used. A systematic review published in October 2023 found that 67% of papers published using TMLE did not indicate how variances were estimated. Of the 33% that did report the method, the majority (73%) relied on the influence function based variance estimator <cit.>. Widespread application of the influence function based variance estimator is aided by available software, as it is used for confidence interval (CI) construction in the R  <cit.> and  <cit.> packages as well as the  <cit.> Python package. 

The influence function based variance estimator's reliance on correctly specifying both exposure and outcome models is particularly concerning in instances where researchers understand either the exposure or outcome mechanism well, but not both. If the outcome process is well known but the exposure process is not well understood, leading to correct specification of the outcome model but misspecification of the propensity model, the influence function based variance estimator is inconsistent <cit.>. In this setting, consistent point estimators may be paired with CIs that do not achieve the nominal level of coverage. Alternatively, if the exposure process is well understood but the outcome process is not, leading to correct specification of the propensity model but misspecification of the outcome model, the influence function based variance estimator is conservative <cit.>. Precision advantages of doubly robust estimators may be lost with application of a conservative variance estimator. For example, we aim to estimate the effect of maternal anemia on birth weight among women with HIV. While most causes of anemia among women with HIV are well understood, doubly robust estimators may be preferred over IPTW estimators due to efficiency advantages. Modeling the relationship between anemia and birth weight among women with HIV is challenging due to complex dependencies between anemia and other drivers of birth weight. Furthermore, low birth weight may result from two separate pathways - preterm birth or fetal growth restrictions - each of which is a complex syndrome with several causes. For these reasons, the outcome model may be susceptible to misspecification, thus application of the influence function based variance estimator may lead to conservative variance estimates in this study.

Alternative, but less commonly used, variance estimation approaches are themselves doubly robust, providing consistent estimates of the variance when only one of the two nuisance models is correctly specified. Two such approaches are the nonparametric bootstrap and the empirical sandwich variance estimator based on M-estimation theory. The nonparametric bootstrap, while doubly robust <cit.>, is computationally burdensome. For large datasets, this computational burden may preclude its use. M-estimation <cit.> has been proposed for use with doubly robust estimators <cit.>, but we are not aware of prior work highlighting the doubly robust property of the empirical sandwich variance estimator. Furthermore, the recent systematic review <cit.> did not mention application of the empirical sandwich variance estimator in any of the papers evaluated. This suggests that the doubly robust property of the empirical sandwich variance estimator is not well known.

This paper demonstrates how M-estimation can be applied for AIPW estimators and TMLE to obtain doubly robust point and variance estimators. In Section <ref>, three common doubly robust estimators are described along with two approaches for inference: the influence function based variance estimator and M-estimation paired with the empirical sandwich variance estimator. Estimators are applied to data from the Improving Pregnancy Outcomes with Progesterone (IPOP) trial to estimate the average causal effect of anemia on birth weight in Section <ref>. The doubly robust property of the empirical sandwich variance estimator is demonstrated in a simulation study and is contrasted with the performance of the influence function based variance estimator in Section <ref>. The benefits and limitations of each variance estimation approach are discussed in Section <ref>. R code for implementing each of the doubly robust estimators is provided on GitHub.



§ METHODS




 §.§ Preliminaries


Consider an observational study that seeks to estimate the effect of a binary exposure (or treatment) X on an outcome Y. Assume n independent and identically distributed copies of O_i=(X_i,Y_i,Z_i) are observed, where Z are a set of baseline covariates. Let the potential outcomes for individual i under exposure and no exposure be denoted by Y_i^1 and Y_i^0 such that by causal consistency Y_i=X_iY_i^1+(1-X_i)Y_i^0. Assume the set of covariates Z satisfy Y^x  X | Z for x ∈{0,1}, i.e., conditional exchangeability. Also assume positivity holds such that (X=x | Z=z)>0 for all z such that dF_Z(z)>0 and x ∈{0,1}, where F_Z is the cumulative distribution function of Z. The estimand is the average causal effect (ACE), μ^1-μ^0 where μ^x=E(Y^x) for x ∈{0,1}. Unless otherwise specified, vectors are assumed to be column vectors.



 §.§ Estimators



Three doubly robust approaches for estimating the ACE are considered: the classic AIPW estimator, the weighted regression AIPW estimator, and TMLE. All three methods combine estimates from a model for the probability of exposure conditional on covariates (propensity model) and a model for the outcome conditional on exposure and covariates (outcome model). The approach for modeling the propensity score is common across the three methods, while there are slight differences in fitting the outcome model. For all three approaches, the propensity score for each participant, e_i=(X_i=1 | Z_i), is estimated. In practice, often the model (e_i)=G_i^Tα is fit, where G_i=g(Z_i) is a vector of predictors for participant i for some user-specified function g of Z_i and α is the vector of regression coefficients from the propensity model. The predicted propensity score for each participant is calculated as ê_i=e(Z_i,α̂)=^-1(G_i^T α̂) where α̂ is the maximum likelihood estimate (MLE) of α. The IPTW is then estimated for each participant as Ŵ_i=X_iê_i^-1+(1-X_i)(1-ê_i)^-1.



  §.§.§ Classic AIPW estimator


For the classic AIPW estimator, the Y | Z,X model is fit and used to predict pseudo-outcomes for each participant under both exposures. That is, a_x(Z_i,γ̂)=Ê(Y_i | Z_i,X_i=x) for x ∈{0,1}, where γ̂ are the MLEs of the parameters from the assumed outcome model. Then, the classic AIPW estimator for the ACE is

    DR_C= μ̂_C^1 - μ̂_C^0
 where μ̂_C^1=n^-1∑_i=1^n  ê_i^-1{X_iY_i-(X_i-ê_i)a_1(Z_i,γ̂)} and μ̂_C^0=n^-1∑_i=1^n  (1-ê_i)^-1{(1-X_i)Y_i+(X_i-ê_i)a_0(Z_i,γ̂)}. The estimator (<ref>) was originally proposed by <cit.> and was further examined in simulation studies, such as <cit.>, <cit.>, and <cit.>. 



  §.§.§ Weighted regression AIPW estimator


An alternative AIPW estimator incorporates the estimated propensity scores when fitting the outcome regression model. The parameters of the Y | Z,X model are estimated with IPTW-weighted maximum likelihood estimation. Pseudo-outcomes are obtained for each participant under both exposures, i.e., b_x(Z_i,β̂)=Ê(Y_i | Z_i,X_i=x) for x ∈{0,1}, where β̂ are the MLEs from the weighted outcome model. The weighted regression AIPW estimator is 

    DR_WR= μ̂_WR^1 - μ̂_WR^0
 where μ̂_WR^x=n^-1∑_i=1^n b_x(Z_i,β̂) for x ∈{0, 1}. The weighted regression AIPW estimator (<ref>) has been evaluated previously <cit.> and is expected to be more stable than (<ref>) when IPTWs are extreme due to bounding within the parameter space <cit.>.



  §.§.§ TMLE


TMLE is an alternative doubly robust estimator that consists of a two-stage process. In the first stage, Y | Z,X is modeled as in Section <ref>, and a_x(Z_i,γ̂) are calculated for x ∈{0,1}. In the second stage, referred to as the targeting stage <cit.>, corrections to outcome regression predictions from stage one are made that incorporate the estimated propensity scores. Different methods have been proposed for fitting the targeting model; here, the weighted regression approach is considered, as it is thought to handle random positivity violations better <cit.>. Specifically, the models (Y_i)=η_x+{a_x(Z_i,γ̂)} are fit for x ∈{0, 1}. Parameters η_x are estimated using weighted maximum-likelihood, with weights (1-X)Ŵ_i and XŴ_i for the x=0 and x=1 targeting model, respectively. In the targeting models, logit-transformed pseudo-outcomes ({a_x(Z_i,γ̂)}) for each participant are included as offset terms, and MLEs η̂_0 and η̂_1 can be thought of as corrections to these pseudo-outcomes <cit.>. The stage two pseudo-outcomes are computed as c_x(O_i, γ̂, η̂_x)=[{a_x(Z_i,γ̂)}+η̂_x]. The TMLE estimator of ACE is 

    DR_TMLE= μ̂_TMLE^1 - μ̂_TMLE^0
 where μ̂_TMLE^x=n^-1∑_i=1^n c_x(O_i, γ̂, η̂_x) for x ∈{0, 1}. Note that for continuous outcomes, such as birth weight considered in Section <ref>, the outcome is scaled prior to implementing TMLE, i.e., Y_i^*=(Y_i-a)/(b-a) is defined for each participant, where (a,b) are the bounds of Y. Then, Y_i is replaced with Y_i^* in the steps above, and c_x(O_i, γ̂, η̂_x)=[{a_x(Z_i,γ̂)}+η̂_x](b-a)+a re-scales the outcome back to the original bounds of Y. 

TMLE methods were originally proposed by <cit.>  and have been extended and applied in a number of settings <cit.>. When parametric models are used in both stages of estimation, (<ref>) is expected to be similar to (<ref>) <cit.>.



 §.§ Variance Estimation


Variance estimation for doubly robust estimators (<ref>), (<ref>), and (<ref>) is not straightforward, because the variance of each estimator depends on estimating the parameters of the nuisance models. Here, we consider two methods for variance estimation: the commonly used influence function based variance estimator and the empirical sandwich variance estimator. Both variance estimators can be used to construct Wald-typed CIs for the ACE. 



  §.§.§ Influence function based variance estimator

Estimators (<ref>), (<ref>), and (<ref>) are suggested by the efficient influence function <cit.>. The commonly used influence function based variance estimator has the form <cit.>

    V̂(DR)=n^-2∑_i=1^n Î_i^2
 where Î_i=ê_i^-1X_iY_i-(1-ê_i)^-1(1-X_i)Y_i-{ê_i^-1(1-ê_i)^-1(X_i-ê_i)}{(1-ê_i)Ŷ_i^1+ê_iŶ_i^0}-DR. Here, pseudo-outcomes Ŷ_i^x equal a_x(Z_i,γ̂), b_x(Z_i,β̂), and c_x(O_i, γ̂, η̂_x) for the classic AIPW, weighted regression AIPW, and TMLE methods, respectively, and DR equals (<ref>), (<ref>), and (<ref>) for the three corresponding estimators. 

The estimators (<ref>) are consistent for the variance of (<ref>), (<ref>), and (<ref>) when both propensity score and outcome models are correctly specified <cit.>. When one of the two nuisance models is misspecified, the influence function based variance estimator is not consistent. As noted in the introduction, the estimator (<ref>) is conservative when the weight model is correctly specified but the outcome model is misspecified <cit.>, and may be conservative or anti-conservative when the outcome model is correctly specified but the weight model is misspecified <cit.>. In empirical studies, the variance estimator (<ref>) for the classic AIPW estimator (<ref>) was found to be unstable for small sample sizes (n<1000), even under correct model specification <cit.>. Similar instability has been reported for (<ref>) with TMLE (n=100) <cit.>.



  §.§.§ Empirical sandwich variance estimator

An alternative to the influence function based variance estimator (<ref>) is the empirical sandwich variance estimator. The estimators (<ref>), (<ref>), and (<ref>) can each be expressed as the solution to a set of unbiased estimating equations. The M-estimator for each set of estimating equations, θ̂, is the solution (for θ) to ∑_i=1^n ψ_q(O_i,θ)=0, where ψ_q(O_i,θ) are the estimating functions corresponding to each of the three estimators q ∈{1,2,3} with parameter vector θ. Parameter vectors and estimating functions for (<ref>), (<ref>), and (<ref>) can be specified as follows.

Classic AIPW estimator

For the classic AIPW estimator, θ=[α, γ, μ_C^1, μ_C^0, DR_C]. The estimating function ψ_1 is specified as

    ψ_1(O_i; θ) = [    ψ_α;    ψ_γ; ψ_μ_C1; ψ_μ_C0;    ψ_C ] = [                                  {X_i-(G_i^Tα)}G_i;                              {Y_i-ϕ^-1(H_i^Tγ)}H_i;         {X_iY_i-(X_i-e_i)a_1(Z_i,γ)}e_i^-1 - μ_C^1; {(1-X_i)Y_i+(X_i-e_i)a_0(Z_i,γ)}(1-e_i)^-1 - μ_C^0;                                 μ_C^1-μ_C^0 - DR_C ]

where ψ_α is the vector of score functions from the propensity score model, ψ_γ is the vector of score functions for the outcome model (with link function ϕ), and H_i=h(X_i,Z_i) is a vector of predictors for the outcome model for participant i for some user-specified function h of X_i and Z_i. The estimating functions ψ_μ_C1 and ψ_μ_C0 are for the causal means μ_C^1 and μ_C^0, respectively, and ψ_C is a delta method transformation of the causal means for estimation of the ACE using the classic AIPW, i.e., DR_C. <cit.> demonstrate that this set of estimating equations is unbiased when at least one of the models is correctly specified.

Weighted regression AIPW estimator

For the weighted regression AIPW estimator, θ=[α, β, μ_WR^1,μ_WR^0, DR_WR]. The estimating function ψ_2 is specified as

    ψ_2(O_i; θ) = [     ψ_α;     ψ_β; ψ_μ_WR1; ψ_μ_WR0;    ψ_WR ] = [        {X_i-(G_i^Tα)}G_i; W_i{Y_i-ϕ^-1(H_i^Tβ)}H_i;       b_1(Z_i,β)- μ_WR^1;       b_0(Z_i,β)- μ_WR^0;    μ_WR^1-μ_WR^0 - DR_WR ]

where ψ_β is the vector of IPTW-weighted score functions for the outcome model. The estimating functions ψ_μ_WR1 and ψ_μ_WR0 are for the causal means μ_WR^1 and μ_WR^0, respectively, and ψ_WR is a delta method transformation of the causal means for estimation of the ACE using the weighted AIPW estimator DR_WR. <cit.> demonstrate unbiasedness of this set of estimating equations when at least one of the nuisance models is correctly specified.

TMLE

With TMLE, θ=[α, γ, η_1, η_0, μ_TMLE^1, μ_TMLE^0, DR_TMLE]. For continuous outcomes, the estimating function ψ_3 is specified as

    ψ_3(O_i; θ) = [       ψ_α;       ψ_γ;      ψ_t1;      ψ_t0; ψ_μ_TMLE1; ψ_μ_TMLE0;    ψ_TMLE ] = [                           {X_i-(G_i^Tα)}G_i;                     {Y^*_i-ϕ^-1(H_i^Tγ)}H_i;         X_ie_i^-1(Y^*_i-[η_1+{a_1(Z_i,γ)}]); (1-X_i)(1-e_i)^-1(Y^*_i-[η_0+{a_0(Z_i,γ)}]);                 c_1(O_i, γ, η_1) - μ_TMLE^1;                 c_0(O_i, γ, η_0) - μ_TMLE^0;                 μ_TMLE^1-μ_TMLE^0 - DR_TMLE ]

where ψ_γ is the vector of score functions for the scaled outcome model. The estimating functions ψ_t1 and ψ_t0 are for the two targeting models, while ψ_μ_TMLE1 and ψ_μ_TMLE0 are for the causal means μ_TMLE^1 and μ_TMLE^0, respectively. Note c_x(O_i, γ, η_x)=[{a_x(Z_i,γ)}+η_x](b-a)+a are the pseudo-outcomes following the targeting step, as outlined in Section <ref>. Finally, ψ_TMLE is a delta method transformation of the causal means for estimation of the ACE using DR_TMLE. Unbiasedness of ψ_3(O_i; θ) when at least one nuisance model is correctly specified is demonstrated in the Appendix. 

Because (<ref>), (<ref>), and (<ref>) are each solutions to unbiased estimating equation vectors, it follows under suitable regularity conditions <cit.> that for each estimator, √(n)(θ̂-θ)  N (0, V(θ) ). Here, V(θ)=A(θ)^-1B(θ){A(θ)^-1}^T where A(θ)=E{-∂ψ_q(O_i; θ) / ∂θ} and B(θ)=E{ψ_q(O_i; θ)ψ_q(O_i; θ)^T}. The quantities A(θ) and B(θ) can be consistently estimated by replacing each expected value with its empirical counterpart, resulting in the empirical sandwich variance estimator V̂(θ̂). The variance of (<ref>), (<ref>), and (<ref>) can be consistently estimated by the bottom right element of its corresponding V̂(θ̂). Importantly, because the estimating equations are unbiased when either the propensity or the outcome model is correctly specified, the empirical sandwich variance estimator is consistent under such model misspecification. This is in contrast to the influence function based variance estimator (<ref>) which requires correct specification of both nuisance models for consistent variance estimation <cit.>. 



§ EXAMPLE: IMPROVING PREGNANCY OUTCOMES WITH PROGESTERONE



Prior research has shown links between maternal anemia and low birth weight <cit.>. The relationship between anemia and birth weight is particularly important to understand for women living with HIV, as they experience anemia more frequently than women without HIV <cit.> and are more prone to adverse birth outcomes <cit.>. The effect of maternal anemia on birth weight among women living with HIV was estimated with data from the IPOP trial. 

The IPOP trial was a randomized, double-blind, placebo-controlled trial conducted in Lusaka, Zambia. Participants were pregnant women aged 18 or older living with HIV on (or intending to initiate) antiretroviral therapy (ART) who had viable singleton pregnancies at less than 24 weeks of gestation at baseline. The 800 enrolled participants were randomized with equal probability to either a placebo or 17 alpha-hydroxyprogesterone caproate (17P) treatment to evaluate the efficacy of 17P for prevention of preterm birth. Additional details about the IPOP trial have been published previously <cit.>. While IPOP was designed to examine the effect of 17P on preterm birth or stillbirth, covariates collected during the trial allow for examination of other causal questions, assuming that the identification assumptions hold.

The three doubly robust estimators described in Section <ref> were applied to IPOP data for estimation of the ACE of maternal anemia on birth weight. The exposure, maternal anemia, was defined as having a baseline hemoglobin level below 10.5 g/dL. The IPOP data were limited to the 782 participants (98%) with no prior preterm births and with measured birth weights and baseline hemoglobin values. Note 14 participants experienced stillbirth. Birth weights for these participants were included in the analysis; we return to this issue in the discussion. It was assumed that potential outcomes are independent of the exposure conditional on the following baseline covariates: trimester (1 or 2), maternal age group (18-24, 25-29, 30-34, 35-39, 40+), number of prior births (0, 1, 2, 3, 4+), alcohol use during pregnancy (yes or no), ART use (yes or no), and maternal height (cm) (Table <ref>). At baseline, no IPOP participants reported tobacco or drug use during pregnancy, so these covariates were not included. Primary analyses were conducted using R and were independently replicated in Python. 

As discussed in the introduction, application of doubly robust estimators for estimating the ACE of maternal anemia on birth weight may offer precision advantages over IPTW estimators, but the correct model specifications are unknown. To demonstrate differences in the influence function based variance estimator and the empirical sandwich variance estimator presented in Section <ref>, each of the doubly robust estimators were applied with three sets of model specifications. First, the covariates in Table <ref> were included in both propensity and outcome models. Then, a naive outcome model was fit that included only the exposure variable (anemia status), while the propensity model included the full set of covariates. Under this specification, the outcome model was reasonably thought to be misspecified due to its simplicity. Finally, a naive propensity model was fit that included only an intercept term, while the outcome model included the full set of covariates. Due to the simplicity of the propensity model, it was reasonably thought to be misspecified. For all estimators, maternal height was modeled using restricted cubic splines with four knots placed at the 5th, 35th, 65th, and 95th percentiles <cit.>. The empirical sandwich variance estimator was computed using  in R <cit.> and  in Python <cit.>. Confidence intervals for the ACE were constructed based on the influence function based variance estimator and the empirical sandwich variance estimator. Estimates for each of the three methods and model specifications along with corresponding 95% CIs are presented in Table <ref>, with 95% CI half-widths compared in Figure <ref>. 

IPOP participants with and without anemia at baseline had similar self-reported alcohol use during pregnancy and similar age and height distributions. Participants without anemia were more likely to have one or more prior births, be in their first trimester, and be on ART. The three doubly robust methods provided similar estimates of the ACE for all model specifications examined, with average birth weights if all women had anemia were approximately 40g lower than birth weights if no women had anemia, though all 95% CIs included zero. 

Precision estimates based on the empirical sandwich variance estimator and the influence function based variance estimator were similar when the full covariate set was included in both the outcome and propensity models. Under the naive outcome model, standard errors for the influence function based variance estimator were larger than those for the empirical sandwich variance estimator, resulting in wider CIs. When the propensity model was naive, standard errors for the influence function based variance estimator were smaller than those for the empirical sandwich variance estimator, resulting in narrower CIs. Note there was more fluctuation in 95% CI half-widths across model specifications for the influence function based variance estimator than the empirical sandwich variance estimator (Figure <ref>). 













§ SIMULATION STUDY


A simulation study was conducted to compare the empirical properties of the influence function based variance estimator and the empirical sandwich variance estimator in conjunction with each of the three doubly robust methods discussed in Section <ref>. 



 §.§ Simulation Setup


Simulations were designed based on the example in Section <ref>, which aimed to estimate the effect of maternal anemia on birth weight among women with HIV. The primary simulations were conducted with n=800 participants, similar to the sample size in the example. Results for a sample size of n=2000 are included in the Supporting Information. The simulation study was conducted in R with code provided on GitHub. 

Three covariates (Z_1, Z_2, Z_3), the exposure (X), and potential outcomes (Y^0 and Y^1) were simulated. The covariate Z_1 was distributed Normal with mean 155 and standard deviation 7.6. Two binary covariates Z_2 and Z_3 were simulated from Bernoulli distributions with means 0.25 and 0.75, respectively. The exposure X was simulated from a Bernoulli distribution with mean (15-0.1 Z_1 + 2.5 Z_2 -1 Z_3 - 0.02 Z_1 Z_2 + 0.005 Z_1 Z_3). Potential outcomes Y^0 and Y^1 under exposure and no exposure, respectively, were simulated from a normal distribution with mean E(Y^x)=1000+ 11.5 Z_1 + 100 Z_2 -15 Z_1 Z_2 +25x-5.5 x Z_1 -30 x Z_2 +20 x Z_1 Z_2 and standard deviation σ=400 for x ∈{0,1}. Under this data generating mechanism, the ACE was approximately -60. Additional simulations were conducted with σ=200 with results in the Supporting Information. 

The estimators DR_C, DR_WR, and DR_TMLE were applied to 5000 simulated samples to estimate the ACE. Estimators were computed four ways: (1) with both propensity and outcome models correctly specified, (2) with propensity models correctly specified but outcome models misspecied, (3) with outcome models correctly specified but propensity models misspecified, and (4) with both models misspecified. Misspecified propensity models included only an intercept and a linear term for (Z_1-155)^2, and misspecified outcome models included only an intercept and linear terms for X and (Z_1-155)^2. 

Both the influence function based variance estimator and the empirical sandwich variance estimator were applied, as described in Section <ref>. For each of the three doubly robust estimators, 95% CIs for ACE were constructed based on each variance estimation approach. Simulation results were summarized by computing empirical bias, relative empirical bias (||/||), average standard error (ASE), empirical standard error (ESE), standard error ratio (SER=ASE/ESE), and empirical 95% CI coverage. Coverage was defined as the proportion of simulated samples in which the CI included ACE. In addition, the ratio of the variance estimate for each simulation and the ESE, i.e., the variance ratio (VR), was computed and summarized. That is, VR=se_s / ESE, where se_s=√(V̂(DR)) for a given doubly robust estimator, model specification, and variance estimator for simulation s.



 §.§ Simulation Results


Both point and variance estimators performed as expected in simulations. The results of the simulation study for n=800 and σ=400 are presented in Table <ref> and Figure <ref>. When at least one model was correctly specified, the classic AIPW, weighted regression AIPW, and TMLE displayed minimal bias, but all were substantially biased when both nuisance models were misspecified. Note bias was further reduced under correct specification of at least one model for n=2000 (Table 4). Under correct specification of both nuisance models, the empirical sandwich and influence function based variance estimators tracked closely with the ESE, resulting in SERs close to one and CIs attaining the nominal level of coverage. When both models were misspecified, both variance estimators tracked with the ESE, but bias was substantial, resulting in CIs with below nominal coverage.  

Differences between the variance estimators are apparent when examining the simulation study results for scenarios where one nuisance model was correctly specified and the other was misspecified. When either the outcome model or the propensity model was misspecified, CIs based on the empirical sandwich variance estimator attained the nominal level of coverage, and estimated standard errors tracked closely with the ESE. In Figure <ref>, note variance ratios for the empirical sandwich estimator are closely clustered around one for these scenarios. In contrast, the influence function based variance estimator was empirically biased when either nuisance model was misspecified. As expected, it tended to overestimate the variance when the outcome model was misspecified, leading to estimated standard errors that tracked above the ESE, SERs above one, and conservative CIs that exceeded the nominal level of coverage. When the propensity model was misspecified but the outcome model was correctly specified, the influence function based variance estimator underestimated the variance, resulting in estimated standard errors below the ESE and SERs below one. The biases observed for the influence function based variance estimator were also present when n=2000 (Table 4 and Figure 3).

In the σ=200 scenarios, the empirical sandwich variance estimator performance was similar to the σ=400 scenario, but the performance of the influence function based variance estimator varied (Tables 5-6 and Figures 4-5). Specifically, when the propensity model was misspecified but the outcome model was correctly specified, the influence function based variance estimator demonstrated minimal empirical bias, resulting in SERs close to one.

In summary, the simulation study demonstrated the theoretical properties explained in Section <ref>. That is, the empirical sandwich variance estimator was consistent when at least one of the two nuisance models was correctly specified, and performed well in simulations. The influence function based variance estimator is not consistent under misspecification of either nuisance model. In line with expectations, the influence function based variance estimator was conservative when the outcome model was misspecified. For this data generating mechanism, it was anti-conservative when the propensity model was misspecified for σ=400 and exhibited minimal bias when σ=200.








§ DISCUSSION


Doubly robust estimators have increased in popularity in the field of causal inference due to their ability to provide consistent point estimates when either an outcome or propensity model is correctly specified <cit.>. In this paper, the commonly-used influence function based variance estimator is compared with the empirical sandwich variance estimator in conjunction with three doubly robust estimators: the classic AIPW estimator, the weighted regression AIPW estimator, and TMLE. These three estimators can be specified as solutions to sets of unbiased estimating equations, and thus their variances can be consistently estimated by the empirical sandwich variance estimator. Because the estimating equations remain unbiased when either the outcome or propensity model is correctly specified, the empirical sandwich variance estimator is a doubly robust variance estimator. In contrast, the influence function based variance estimator is consistent only when both outcome and propensity models are correctly specified. This discrepancy may be of particular importance in settings where the functional form of one of the nuisance models is known but the other is not. For example, with randomized trials, randomization to treatment means that the exposure mechanism is known. Doubly robust estimators can offer improved efficiency in this setting, but some of that efficiency gain may not be realized with application of the influence function based variance estimator, which is conservative if the outcome model is misspecified. 

While consideration in this paper was limited to finite-dimensional parametric modeling approaches, machine learning approaches can also be used with doubly robust estimators <cit.>. Machine learning approaches allow for more complex functional forms for continuous covariates than fully parametric approaches and include higher order interaction terms that may not be included in investigator-specified parametric models, often leading to estimators that are more robust to model misspecification <cit.>. However, convergence of machine learning algorithms is typically slower and consistent variance estimation is more challenging than with parametric modeling approaches. Machine learning methods are typically not compatible with the estimating equations approach discussed in this paper, though alternative methods have been developed for doubly robust variance estimation in this context for TMLE estimators <cit.>. 

In the example, the effect of anemia on birth weight among women with HIV was estimated by applying the three doubly robust methods to data from the IPOP randomized trial. Average birth weights if all women had anemia were estimated to be approximately 40g lower than birth weights if no women had anemia, though all 95% CIs included zero. The analysis included 14 participants who experienced stillbirth. Because stillbirth is a competing outcome for live birth weight, future work could reexamine this problem using alternative approaches for competing events.

Use of the empirical sandwich variance estimator in conjunction with doubly robust estimators is not new, though its doubly robust property has not been emphasized. <cit.> discuss use of the empirical sandwich variance estimator for estimating the variance of the classic AIPW estimator and note that it tends to be more stable than the influence function based variance estimator. However, based on the systematic review discussed in the introduction <cit.>, the influence function based variance estimator is frequently applied and the empirical sandwich variance estimator does not appear to be widely used with TMLE. While some existing software packages compute empirical sandwich variance estimators for doubly robust estimators (e.g., the “causaltrt” procedure in SAS and the “dr” procedure in Stata), other popular software packages compute the influence function based variance estimator but not the empirical sandwich variance estimator (e.g., the  and  packages in R and the  package in Python). We provide example R and Python code for estimating the variance of each of the three doubly robust estimators with the empirical sandwich variance estimator, which we hope will allow for easier implementation of doubly robust point and variance estimation.



§ SUPPLEMENTARY MATERIAL
 
R and Python code for computing the different estimators along with the corresponding standard error estimators is available at https://github.com/bonnieshook/DoublyRobustVariance.



§ BIBLIOGRAPHY


	apalike
	

	


§ ACKNOWLEDGEMENTS
 
This work was supported by US National Institutes of Health under award numbers R01 AI157758, R01 AI085073, P30 AI050410, and K01 AI177102. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. The authors thank Bradley Saul at the University of North Carolina at Chapel Hill for coding support. 


	
	
		
		


§ APPENDIX




 §.§ Consistency proof of TMLE estimator under correct specification of either the outcome or propensity model



Recall the TMLE estimating function ψ_3 is specified as

    ψ_3(O_i; θ) = [       ψ_α;       ψ_γ;      ψ_t1;      ψ_t0; ψ_μ_TMLE1; ψ_μ_TMLE0;    ψ_TMLE ] = [                           {X_i-(G_i^Tα)}G_i;                     {Y^*_i-ϕ^-1(H_i^Tγ)}H_i;         X_ie_i^-1(Y^*_i-[η_1+{a_1(Z_i,γ)}]); (1-X_i)(1-e_i)^-1(Y^*_i-[η_0+{a_0(Z_i,γ)}]);                 c_1(O_i, γ, η_1) - μ_TMLE^1;                 c_0(O_i, γ, η_0) - μ_TMLE^0;                 μ_TMLE^1-μ_TMLE^0 - DR_TMLE ]
 where θ=[α, γ , η_1 , η_0 , μ_TMLE^1 , μ_TMLE^0 , DR_TMLE] is the vector of parameters.

Let


θ̃
= [α̃ , γ̃ , η̃_1 , η̃_0 , μ̃_TMLE^1 , μ̃_TMLE^0 , DR_TMLE] be the solution of 
E{ψ_3(O_i; θ)} = 0 and DR_TMLE be the proposed TMLE estimator of ACE, obtained by solving the estimating equation n^-1∑_i ψ_3(O_i;θ) = 0. Then, if either the outcome model or the propensity score model is correctly specified, √(n) (DR_TMLE - ACE) d⟶ N(0,σ^2(θ̃)), where σ^2(θ̃) is the bottom right element of V(θ̃) = V(θ) |_θ = θ̃, and
V(θ) = A(θ)^-1B(θ)A(θ)^-T,
with
A(θ) = E{-∂ψ_3(O_i; θ) / ∂θ},
B(θ) = E{ψ_3(O_i; θ) ψ_3(O_i; θ)^⊤},
which is a doubly robust variance estimator of DR_TMLE.

Proof.

Let θ̂ be the solution of the estimating equation n^-1∑_i ψ_3(O_i;θ) = 0.
Under suitable regularity conditions <cit.>, √(n) (θ̂ - θ̃)
	d⟶ N(0,V(θ̃)). Note the final element of θ̂ is DR_TMLE and therefore
√(n) (DR_TMLE - DR_TMLE)	d⟶ N(0,σ^2(θ̃)). It remains to show that if either the outcome model or the propensity score model is correctly specified, DR_TMLE = ACE.

First, let 
θ_0
= [α_0, γ_0, η_10, η_00, μ_TMLE0^1, μ_TMLE0^0, DR_TMLE0]
be the true parameter values from the data generating process.
Then,
P(X_i = 1 | Z_i) = e(Z_i, α_0) = (G_i^T α_0)
and
E(Y_i^* | Z_i, X_i) = ϕ^-1(H_i^T γ_0) 
= X_i a_1(Z_i, γ_0) + (1-X_i) a_0(Z_i, γ_0).

Now assume the outcome model is correctly specified, 
i.e., γ̃ = γ_0 and
E(Y_i^* | Z_i, X_i)
= X_i a_1(Z_i, γ_0) + (1-X_i) a_0(Z_i, γ_0).
From E{ψ_t1(O_i; θ̃)} = 0,

    0 
    	=   
    	E{ψ_t1(O_i; θ̃)}
    
    	=   
    	E[E{ψ_t1(O_i; θ̃) | X_i}]
    	
    
    	=   
    	E[E{
    		X_i 
    		e(Z_i, α̃)^-1 
    		(Y_i^* - [η̃_1 + {a_1(Z_i, γ_0)}]) 
    	| X_i}]
    	
    
    	=   
    	E{
    		e(Z_i, α̃)^-1 
    		(Y_i^* - [η̃_1 + {a_1(Z_i, γ_0)}]) 
    	| X_i = 1} P(X_i = 1)
    	
    
    	=   
    	E[
    		E{
    			e(Z_i, α̃)^-1 
    			(Y_i^* - [η̃_1 + {a_1(Z_i, γ_0)}]) 
    			| Z_i, X_i = 1
    		}
    		| X_i = 1
    	] 
    	P(X_i = 1)
    	
    
    	=   
    	E[
    		e(Z_i, α̃)^-1{
    			E(Y_i^* | Z_i, X_i = 1)
    			-
    			[η̃_1 + {a_1(Z_i, γ_0)}]
    		}
    		| X_i = 1
    	] 
    	P(X_i = 1)
    	
    
    	=   
    	E[
    		e(Z_i, α̃)^-1{
    			a_1(Z_i, γ_0)
    			-
    			[η̃_1 + {a_1(Z_i, γ_0)}]
    		}
    		| X_i = 1
    	] 
    	P(X_i = 1)
    	
    
    	=   
    	E[
    		e(Z_i, α̃)^-1
    		F(η̃_1)
    		| X_i = 1
    	] 
    	P(X_i = 1)

where
F(t) = (d) - (t + d) and d = {a_1(Z_i, γ_0)}.
Since the function 
F(t)
is strictly decreasing in t,
it can be shown that the last line in the above is strictly decreasing in η̃_1 and attains 0 only at η̃_1 = 0.
Therefore, η̃_1 = 0.
Then, from E{ψ_μ_TMLE1(O_i; θ̃)} = 0,

    μ̃_TMLE^1 
    	=    
    	E(Y_i, TMLE^1)
    	
    
    	=   
    	E{[η̃_1 + {a_1(Z_i, γ_0)}](b-a) + a
    	}
    
    	=   
    	E{a_1(Z_i, γ_0)}(b-a) + a
    	
    
    	=   
    	E[
    		E{
    			Y_i^* | Z_i, X_i = 1
    		}
    	]
    	(b-a) + a
    	
    
    	=   
    	E[
    		E{
    			(Y_i-a)/(b-a) | Z_i, X_i = 1
    		}
    	]
    	(b-a) + a
    	
    
    	=   
    	E{E(Y_i| Z_i, X_i = 1)}
    
    	=   μ^1.

Similarly, μ̃_TMLE^0  = μ^0.

Next, assume the propensity score model is correctly specified, i.e., α̃ = α_0
and
P(X_i=1|Z_i) = e(Z_i, α_0).
From E{ψ_t1(O_i; θ̃)} = 0,

    0 
    	=   
    	E{ψ_t1(O_i; θ̃)}
    
    	=   
    	E[E{ψ_t1(O_i; θ̃) | Z_i}]
    	
    
    	=   
    	E[
    		E{
    			X_i 
    			e(Z_i, α_0)^-1{(Y_i^1-a)/(b-a) 
    			- 
    			[η̃_1 + {a_1(Z_i, γ̃)}]) }
    			| Z_i
    		}
    	]
    	
    
    	=   
    	E(
    		E(X_i | Z_i)
    		e(Z_i, α_0)^-1 
    		[
    			E{(Y_i^1-a)/(b-a)|Z_i}
    			- 
    			[η̃_1 + {a_1(Z_i, γ̃)}]
    		]
    	)
    	
    
    	=   
    	E[E{(Y_i^1-a)/(b-a)|Z_i}]
    	-
    	E([η̃_1 + {a_1(Z_i, γ̃)}])
    	
    
    	=   
    	(μ^1-a)/(b-a)
    	-
    	E([η̃_1 + {a_1(Z_i, γ̃)}]).

Therefore,
E([η̃_1 + {a_1(Z_i, γ̃)}]) = (μ^1-a)/(b-a).
Then, from E{ψ_μ_TMLE1(O_i; θ̃)} = 0,

    μ̃_TMLE^1 
    	=    
    	E{c_x(O_i, γ, η_x)}
    
    	=   
    	E{[η̃_1 + {a_1(Z_i, γ̃)}](b-a) + a
    	}
    
    	=   μ^1.

Similarly, μ̃_TMLE^0  = μ^0.

In conclusion,
if either the outcome model or the propensity score model is correctly specified, 
DR_TMLE = μ̃_TMLE^1  - μ̃_TMLE^0  =  μ^1 - μ^0 = ACE
which proves the stated result.



	




















		
	
	
	
