{'fraction_non_alphanumeric': 0.03636882555603123, 'fraction_numerical': 0.012598617821539967, 'mean_word_length': 4.612585812356979, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Robust variance estimation and inference for causal effect estimation

We consider a longitudinal data structure consisting of baseline covariates, time-varying treatment variables, intermediate time-dependent covariates, and a possibly time dependent outcome. Previous studies have shown that estimating the variance of asymptotically linear estimators using empirical influence functions in this setting result in anti-conservative estimates with increasing magnitudes of positivity violations, leading to poor coverage and uncontrolled Type I errors. In this paper, we present two alternative approaches of estimating the variance of these estimators: (i) a robust approach which directly targets the variance of the influence function as a counterfactual mean outcome, and (ii) a non-parametric bootstrap based approach that is theoretically valid and lowers the computational cost, thereby increasing the feasibility in non-parametric settings using complex machine learning algorithms. The performance of these approaches are compared to that of the empirical influence function in simulations across different levels of positivity violations and treatment effect sizes.


Introduction
A number of estimators are available for the treatment specific mean outcome parameter (and the corresponding causal contrasts) based on longitudinal data structures, such as inverse probability weighting (IPW) [Horvitz andThompson, 1952, Robins, 1998], double robust augmented IPW (AIPW) [Robins and Rotnitzky, 1992, Robins et al., 1994, Robins and Rotnitzky, 2001, Rotnitzky and Robins, 2005, and targeted minimum loss-based estimation (TMLE) [van der Laan and Gruber, 2011]. Variance estimation for each of these are conventionally achieved by using their corresponding influence functions (IF) on the empirical distribution or by resampling methods such as the non-parametric bootstrap. However, a number of shortcomings exists with these variance estimation approaches. In particular, no theory for exists for the non-parametric bootstrap when using data adaptive methods for estimation of nuisance parameters, and both IF-based and bootstrap based confidence intervals can become anti-conservative with increasing levels of practical positivity violations. For example, van der Laan et al. [van der Laan and Gruber, 2011] found IF-based variance estimates for the intervention specific mean outcome that were anti-conservative when compared with the Monte-Carlo variance of the TMLE, leading to invalid confidence intervals. Petersen el al. [Petersen et al., 2012] found poor coverage for influence function-based confidence intervals, owing to both a result of practical positivity violations and relatively rare outcomes. This behaviour is especially true under sparsity in finite samples, even when the assumptions for asymptotic validity of these estimators hold [Petersen et al., 2012]. As a consequence, statistical inference based on these estimators of variance becomes unreliable when the treatment mechanism practically or theoretically violates the underlying positivity assumption.
Additionally, under sparsity issues, the estimated variance may also fail to raise a red flag for unreliable statistical inference [Petersen et al., 2012]. For example, these estimates of the asymptotic variance are not sensitive to theoretical violations of the positivity assumptions under which the asymptotic variance would be infinity, i.e. when positivity fails. Consequently, it is less likely that the analyst will be able to determine if the data at hand provides insufficient information to estimate the desired causal parameter with any reasonable degree of accuracy.
Previous work [Petersen et al., 2012[Petersen et al., , 2014 proposed estimating the asymptotic variance of the estimator with a parametric bootstrap-based on a fit of the density of the data generating distribution, involving estimation of the treatment mechanism and the G-computation factor of the likelihood. This proposal corresponds with evaluation of the variance of a given estimator using the data at hand as a given data generating experiment. The consistency of this estimator relies on correct specification of both the treatment mechanism and the G-computation factor. This parametric bootstrap integrates over sparse events and therefore will explode the variance. An extremely large number of samples is therefore needed to get the true variance under this Monte Carlo scheme. As a consequence, this parametric bootstrap-based variance estimate was only proposed as a measure to raise a red flag for unreliable statistical inference. In addition, in the context of sparsity, one needs to sample a large number of bootstrap samples and refit the likelihood in each iteration in order to obtain a valid evaluation of the estimator variance, in order to capture the rare observations that nonetheless heavily contribute to this variance. Thus, this semi-parametric bootstrap method is extremely computer intensive, making this Monte Carlo scheme an intractable method for complex estimators and complex data generating distributions.
In this article we use analytic expressions to compute the variance of the efficient influence function (EIF) [Hampel, 1974, Robins andRotnitzky, 1992] which provide the asymptotic variance of estimators solving the estimating equation corresponding to this function. These analytic expressions naturally integrate over the rare observations, and thereby avoid the finite sample bias in variance estimation using standard influence curve or non-parametric bootstrap based methods due to rare observations mentioned above. With this, we construct plug-in type estimators of these asymptotic variances that are consistent if both the treatment mechanism and treatment specific means of specified outcomes are consistently estimated. These estimators require estimation of the treatment mechanism and several treatment specific means of specified outcomes (defined as a function of the observed data structure, indexed by the estimator of the treatment mechanism), which can be estimated with either an estimating equation type IPW estimator or an efficient double robust method such as a targeted minimum loss-based estimator (TMLE). The resulting variance estimator, unlike current alternatives based on taking the variance of the empirical influence function, or using a non-parametric bootstrap, will become very large whenever the estimated treatment mechanism reflects practical or theoretical violations of the positivity assumption.
While this newly presented approach performs well in estimating the asymptotic variance of estimators solving the estimating equation corresponding to the EIF, a lower finite sample variance should be expected for substitution based estimators such as TMLE [van der Laan and Gruber, 2011], due to the guaranteed parameter boundaries provided by the estimator. We therefore additionally present a second bootstrap based approach of estimating the finite sample variance. The approach does not require re-estimation of the treatment mechanisms and the q-factor of the likelihood and therefore reduces the computational burden. It is asymptotically consistent under reasonable assumptions; namely, the same essential assumptions needed for the estimator of the target parameter itself to be asymptotically linear. The resulting reduction in the computational load (compared to a fully non-parametric bootstrap approach which refits the likelihood for each iteration) allows for a more tractable approach at estimating the variance.

Organization of this paper
In Section 2, we formally define the observable data, likelihood, and statistical model for its distribution. Our target parameter of the treatment specific mean is defined along with its EIF. We briefly review the causal model and assumptions under which this statistical quantity corresponds with the desired causal parameter of the counterfactual distribution, along with the currently common approach of influence function (IF) based estimator variance estimation.
Section 3 presents an approach for robust estimation of the variance of the EIF under sparsity. The expression for the variance of the efficient influence function is presented along with both an IPW and TMLE based approach at estimating this parameter. To help illustrate, an example is given for a point treatment setting under a static treatment regime. Advantages of this new approach are covered. The Appendix generalizes the approach to working marginal structural working models and provides proofs.
Section 4 discusses the second approach of estimating the estimator variance using the bootstrap, using a modified TMLE. This bootstrap approach does not require re-estimation of the treatment mechanisms and the q-factor of the likelihood, therefore reducing the computational intensity required.
Section 5 illustrates the performance of the variance estimators presented in Sections 3 and 4 by applying them in simulations to both a single time-point and longitudinal setting. Results show that the robust approach at estimating variance is conservative for substitution based estimators of the mean outcome, while the bootstrap approach results in estimates close to the observed Monte-Carlo variance. The resulting confidence intervals are demonstrated to be valid under the newly proposed variance estimation approaches, while the bootstrap approach is shown to retain higher statistical power.
We conclude with a discussion in Section 6, which reviews the results, benefits of this new approach, potential limitations, and future directions.

Definition of data and statistical estimation
Consider a longitudinal study in which subjects are seen at each time point t from t = 0, 1, . . . , K + 1. The observable data structure on a randomly sampled subject is where L(0) includes all baseline covariates, A(t) denotes an intervention node at time t, and L(t) denotes all timevarying covariates at time point t, measured between the intervention nodes A(t − ) and A(t), where for notational convenience we define t − ≡ t − 1. Our outcome of interest Y = L(K + 1) is an outcome measured after the final treatment A(K). We observe n independent and identically distributed (iid) copies copies O i : i = 1, . . . , n, of O.
The likelihood L(O) for the observable data is the product of conditional probabilities such that the likelihood for subject i is whereX(t) ≡ (X(1), X(2), . . . , X(t)), A(−1) = L(−1) = ∅, and p 0 (o) denotes p 0 (O = o) under the true distribution P 0 where we assume O is discrete for sake of presentation.
The statistical model M for the data involves assumptions, if any, only on the conditional distributions of A(t), given P a(A(t)) = (L(t),Ā(t − )), t = 0, . . . , K. Let denote the G-computation formula for the post-intervention distribution of an intervention that setsĀ(K) = d(l(K)) [Robins, 1986]. We use the notation P L(t) for a conditional distribution of L(t), given P a(L(t)) = (L(t − ),Ā(t − )).
i.e. the mean of the outcome at time K + 1 under this distribution. We note that Ψ : M → IR represents a target parameter mapping on the statistical model to the real line. Defining t + ≡ t + 1, the EIF of Ψ at P is given by [Robins, 2000, Bang and where It should be noted that g 0:t − (Ā(t − ),L(t − )) represents the cumulative probability of treatment up to time t − 1 and is defined by recursive regression, starting at t = K + 1 and moving backwards in time. For notational convenience, we let H 0 = 1 so that

Causal model
Under additional assumptions about our data generating where our target statistical estimand is equal to the mean of the counterfactual outcome Y d under intervention to set the vector of treatment nodes to value d(l(K)). Specifically, for interventions of interest d ∈ D, we assume sequential randomization [Robins, 1986] Y d ⊥ ⊥ A(t)|L(t),Ā(t − ) : t = 0, 1, .., K and positivity [Robins, 1999] P (Ā(t) = d(l(t))|L(t),Ā(t − ) = d(l(t))) > 0 a.e. : t = 0, 1, . . . , K Regarding the assumption of positivity, we note that as P ( ] → ∞ and the previously mentioned practical and theoretical positivity violations.

Review of Influence Function based variance
Recall that an estimatorΨ(P n ) is considered to be asymptotically linear if and only if for some mean 0 finite variance influence function D(P 0 )(O) [Hampel, 1974]. If an estimator is asymptotically linear, then it will be asymptotically normal with variance equal to the variance of the influence function over n. The asymptotic variance of the estimator can therefore be consistently estimated with the variance of the empirical influence function D(P n )(O), i.e.v ar[Ψ(P n )] = var[D(P n )(O)]/n, which implies an asymptotically valid confidence interval.

Targeted minimum loss based estimation (TMLE)
One such estimator that solves the the estimating equation corresponding to the efficient influence function for intervention specific mean outcomes is Targeted Minimum Loss-based Estimation [van der Laan and Gruber, 2011]. This estimator solves the estimating equation by forming an intial fit of the q 0 portion of the likelihood and subsequently perturbing it such that the estimating equation is solved. We assume use of this estimator for our estimation problem, such that our attention is focused on estimation of the estimator's variance. Note that our proposed variance estimators also apply to estimating equation approaches, such as the double robust augmented IPW (AIPW) [Robins and Rotnitzky, 1992, Robins et al., 1994, Robins and Rotnitzky, 2001, Rotnitzky and Robins, 2005.

Semi-targeted estimation of the EIF variance
We can directly target the variance of D(P 0 )(O) as an expectation, allowing us to estimate the variance as a mean. The following describes how to obtain a TMLE of the variance of each component of the EIF σ 2 t in the setting of a scalar parameter. We provide a proof for the more general working MSM setting in the Appendix for the interested reader.

Expression for variance of the EIF for EY d
Under regimens d(l(K)), we have Using the expression for H t (g) from Equation (3), and first taking the conditional expectation w.r.t.Ā(t − ) given X = (L d : d), it follows that this can be written as: where we define g 0 (d(l(−1)),L d (−1)) = 1 so that the term at . This is simply a sum of expectations over t ∈ {0, 1, . . . , K + 1}. For notational convenience, we re-write Equation (7) as for the specified function , the mean of a real valued function (indexed by d(l) itself) ofL d (j), which needs to be estimated based on the longitudinal data structure . . , n, so that we can represent the observed data structure as , and we wish to estimate the statistical target parameter where again we assume l(t) is discrete for sake of presentation.

Estimation of variance of the EIF
With the expression for the variance of the efficient influence function in hand (Equation 8), we can now form estimators which target this parameter.Q 0 and g 0 are not known in practice, though estimatesQ * n and g n will be readily available if estimating EY d using a double robust estimator such as TMLE, thus providing us with the observed outcome S d t (Q * n , g n )(L(t)). Treating this variable as our new time point specific outcome, our goal is to estimate the mean of this variable over the post-intervention distribution ofL ) represent the observed data structure.
One possible approach to estimating each of the components (Equation (9)) is to use a simple IPW estimator [Horvitz and Thompson, 1952] . However, such an estimator would still be subject to underestimation of the variance by ignoring the contribution of observations that selected a likely treatmentĀ i , even though their probability of following d(l) is very small. In other words, subjects i with small probabilities of following d(l) would be unlikely to be observed withĀ i = d(l) resulting in an indicator value of 0 for the numerator and, consequently, a contribution of 0 to the IPW estimator. Therefore, we stress that it is important to use a plug-in estimator such as TMLE [van der Laan and Gruber, 2011] to estimate this parameter. A plug-in estimator will integrate over alll(t) in the support of P d t,n and thus contribute many large values of S d t,n (Q * n , g n ) when there are practical or theoretical positivity assumption violations. In addition, the TMLE is a double robust estimator so that it will yield a consistent estimator of this variance if g n is consistent for the true g 0 . GivenQ 0 , g 0 , we will now provide a succinct summary of the TMLE of σ 2,d 0,t = E P d 0 Z d (t) that is based on iterative sequential regression. Note that this iterative sequential regression approach is similar to the one presented by van der Laan et al. [van der Laan and Gruber, 2011] for the intervention specific mean outcome parameter. Denote the counterfactual of Z d (t) under treatment d with Z d,d (t), and let P d 0 be the G-computation formula [Robins, 1986] , which can be represented as a series of iterated conditional expectations The EIF for this target parameter σ 2,d t is given by In this case, the variance of the efficient influence curve is represented as If using a double robust estimator for the estimation of EY a such as TMLE, we are provided with estimators g n and respectively. The second term in the final expression of Equation (10) is easily estimated with the empirical distribution. Given g 0 andQ 0 , the first term can be represented as the mean of a counterfactual S a (L a (0) (0), A] that respects this range, representing it as a logistic regression fit bounded by (a, b), and updating the latter by fitting a univariate logistic regression with clever covariate I(A = a) and observational weight 1/g 0 (a | L(0)), using the initial fit as an off-set. Regarding the initial fit E n [S a | A = a, L(0)], recall from above that S a is a function of L(0) which results in the initial fit being exactly (Y −Q 0 (a, L(0))) 2 /g 0 (a | L(0)) such that regression is unneeded. Following the update step, the TMLE of 3.2 Advantages of this plug-in estimator of the asymptotic variance of the EIF Since σ 2 0 /n equals the asymptotic variance of an asymptotically efficient estimator, it provides a good measure of the amount of information in the data for the target parameter of interest. Therefore, it is sensible to view σ 2 0 /n as a measure of sparsity for the target parameter of interest. If g n is a good estimator of g 0 , then our proposed plug-in estimatorσ 2 n is much less subject to under-estimation due to sparsity than currently available estimators such as the sample variance of the estimated influence function, and the bootstrap-based estimate of the variance of an efficient estimator. Indeed, the non-parametric bootstrap generally is not valid, except when using a parametric model to estimate g 0 andQ 0 which will never capture a true model in practice. This plug-in estimateσ 2 n represents a variance of the estimate of the EIF which involves the integration of rare combinations of treatment and covariates that are unlikely to occur in the actual sample.
In particular, if there are theoretical violations of the positivity assumption, then this true variance σ 2 equals infinity, and, if g n approximates g 0 well, then also the estimateσ 2 n will generate very large values, demonstrating the lack of identifability and thereby raising a red flag for finite sample sparsity bias in the estimators (beyond the large confidence intervals generated byσ 2 n ). If there are serious practical violations of the positivity assumption, then the estimate of this variance should be imprecise, since it is itself a highly variable estimator of a weakly identifiable parameter.

Variance estimation for substitution based estimators
The plug-in estimator of the asymptotic variance of the EIF presented above is superior to the more common approach of taking the empirical EIF variance over the sample (i.e., var[D * (P n )(O)]), in that there is a much stronger contribution of combinations of treatment and covariates that are unlikely to occur in the actual sample. In finite samples, however, the use of substitution based estimators such as TMLE (which are guaranteed to solve the EIF within a bounded range) are often observed to have smaller variance than their asymptotic variance. This is due to the mere fact that they are guaranteed to respect the global constraints of the statistical model and target parameter mapping. That is, as opposed to estimating equations that tend to result in estimates outside parameter boundaries as the EIF variance increases, the use of substitution estimators in finite samples will retain an estimator variance that is smaller than the EIF variance divided by the sample size, n. Thus, using the newly presented robust EIF variance method can result in over-estimation of the estimator variance for these types of estimators. This therefore motivates us to develop an estimator of variance that is less conservative, i.e. more aligned to the true variance of substitution based estimators (such as TMLE).
One alternative approach at estimating the variance for substitution based estimators is to conduct a non-parametric bootstrap. The n observations are sampled with replacement and used to form an estimate of the parameter over B iterations. However, as stated above, the non-parametric bootstrap is generally invalid and not theoretically supported. Additionally, this is a very computer intensive method that usually requires estimating the full likelihood (i.e., P d 0 ) of the longitudinal data structure within each sampled iteration and is therefore normally infeasible in practice unless conducted within an a priori selected smaller parametric statistical model such as logistic regression.
In this section we present an alternative bootstrap based approach that, unlike the standard non-parametric bootstrap, is both computationally feasible and theoretically valid. That is, this bootstrap approach allows us to estimate the variance of the estimator while avoiding re-estimation of g 0 andQ 0 . To facilitate this, we propose a modification of the usual TMLE such that the targeting step is separated from the initial estimation ofQ 0 . Recall that the typical TMLE, as implemented, pivots between the targeting step and the initial estimator for the next regression (preventing us from separating the intitial fit from the targeting step). We propose a minor modification of the TMLE that separates these steps, first estimating all of the initial regressions and subsequently targeting the fits in a separate step. This modified TMLE can then be boostrapped via only the targeting step. Note that, because the modified TMLE has the same asymptotic behavior as the original TMLE, the bootstrap is theoretically supported and will lead to valid inference.

Modified TMLE for EY d
To reduce the computational burden that bootstrapping requires, we first present the modified TMLE approach at estimating the parameter EY d . This parameter can be estimated by the following steps: 1. Estimate g 0:t − (Ā,L) : t = 1, 2, . . . , K + 1 and denote the fits g 0:t − ,n .
2. Determine the range (a, b) for EY d . Recursively for t = K + 1, K, . . . , 1, estimate the conditional expectation ] respecting this range. Denote the fitsQ d t,n . We stress that this step is cruicially different than the typical TMLE, in that all of the initial regression fits are done simultaneously.
3. For time t = K + 1, target the initial fitQ d K+1,n by using a parametric submodel respecting the range (a, b) by adding the covariate I(Ā(K) = d(l(K))) and observational weight 1/g 0:K,n (on the logistic scale), using the initial fits as off-set, and setting Y as the dependent variable. Denote this updated fit asQ d, * K+1,n . 4. GivenQ d, * K+1,n , we can recursively for t = K, K, K − 1, . . . , 1 target the initial fitsQ d t,n by using parametric submodels respecting the range (a, b), adding the covariates I(Ā(t − ) = d(l(t − ))) and observational weight 1/g 0:t − ,n (on the logistic scale), using the initial fits as off-set, and settingQ d, * t + ,n as the dependent variable. Denote the updated fits asQ d, * t,n . 5. At t = 1, we have the estimateQ d, * 1,n , which now is a function of only L(0). Taking the average ofQ d, * 1,n w.r.t. the empirical distribution of L i (0) gives us the desired TMLE estimate of EY d .
This estimator also solves the EIF and is therefore also asymptotically linear and efficient. We note that the analysis of this TMLE is identical to the typical TMLE presented by van der Laan et al. [van der Laan and Gruber, 2011], with the only difference being the initial estimator fits. Here the initial estimators are the original ones, whereas the previous TMLE is implemented with initial estimators using the targeted fits for the outcome.
We emphasize that this estimator is proposed for the sake of the bootstrap method for variance estimation. It is recursive, in that each fitQ d t,n is dependent upon the fit at t + . As opposed to the TMLE, the recursive nature of this TMLE is self contained within each step. In other words, each estimation step in this TMLE can be performed independently of the other steps. This allows the analyst to form all of the initial fits P n prior to performing any of the targeted updates.

Bootstrapping the modified TMLE
The new TMLE approach presented above can be bootstrapped in a fully non-parametric manner, such that observations are drawn with replacement prior to fitting the full likelihood P d 0 and used to form an estimates of the parameter, leading to an estimate of estimator variance. Our recommendation is to only bootstrap the targeting step. More specifically, once the fits g 0:t − ,n andQ d t,n are formed for t = 1, 2, . . . , K + 1, steps 3-5 above are carried out in the bootstrap such that for b = 1, 2, . . . , B we have for a user selected submodel P ( ). The estimator variance is then estimated by taking the variance over the bootstrapped estimates, i.e., var(Ψ(Q n )) = var[Ψ(Q * n,b )]. We emphasize that this TMLE is provided such that we do not need to re-estimateQ n , g n , i.e. it is not a function of the data. If g n → g 0 andQ n → Q, then this TMLE is asymptotically linear with influence function D * (Q, g 0 ). This is conservative relative to the variance of the actual TMLE that is estimated with g n fitted on the data, when g n is consistent.

Simulations
Simulation studies presented in this section illustrate the performance of the two estimators of variance for the estimation of the effect of treatment in both a point treatment setting, and in a longitudinal observational study setting with three time points (i.e. K + 1 = 3) with time-dependent confounding. To analyze the performance, we first compare the variance estimation approaches covered above in estimating the estimator variance. Both the AIPW and TMLE estimators are considered in order to demonstrate the difference in estimating equations and substitution based estimators, respectively. The mean of the variance estimates are compared to the Monte-Carlo variance of each estimator. Additionally, we present the empirical coverage, Type I, and Type II errors resulting from each variance estimation approach. The Monte-Carlo variance of each variance approach is also reported. All analyses were conducted on R version 3.1.1 [Team, 2014].

Point treatment setting
Consider a point treatment setting, such as patient enrollment into a care program, in which the treatment A(0) is only assigned at a single time point. We are interested in determining whether the treatment of interest has a significant effect on the outcome on an additive scale. Our target parameter is therefore the difference of the mean outcomes under treatment and control, i.e., ψ 0,1 ≡ E Y 1 − E Y 0 . Under this setting, the simulated data were generated as follows: L 2 (0) ∼ N (−0.55 + 0.5W 1 + 0.75W 2 , 0.5 2 ) g 0,0 (P a(A(0))) = logit −1 (β p − (β p + 2.5)W 1 + 1.75W 2 +(β p + 3.2)L 1 (0) − 1.8L 2 (0) + 0.8L 1 (0)L 2 (0))) with a positivity associated parameter β p ranging from −2 (minor positivity violations) to 1 (strong practical positivity violations) and the treatment effect associated parameter β ψ0 ranging from 0 (no treatment effect) to 1 (strong treatment effect). Here, L 1 (0) and L 2 (0) are not time-dependent confounders and are therefore considered baseline covariates along with (W 1 , W 2 ), which affect both the treatment and the outcome.

Longitudinal treatment setting
For the longitudinal setting, we considered a treatment A(t) which was allowed to vary over time as a counting process. That is, if A(t) = 1 then we have that A(t) = 1 where X(t) = (X(t), X(t + 1), . . . , X(K)). Similar to the point treatment setting, we are interested in whether the treatment of interest has a significant effect on the outcome at the final time point t * = 3 on an additive scale. Thus, our target parameter is the difference of the mean outcomes under treatment and control at this final time point, i.e., ψ 0, . Under this setting, data for the first time point was generated in the same manner as the point treatment setting in Section 5.1.1 above. For the remaining two time points, the data were generated conditional on survival (i.e. L 3 (t − ) = 0) as follows: Similar to the point treatment setting, the treatment effect associated parameter β ψ0 also ranged from 0 to 1. We note, however, that the positivity issues faced in this scenario will be even more severe due to the higher number of combinations of treatment over time, which result in smaller probabilities. We therefore considered only β p values from −2 to 0 and imposed a truncation level of 0.001 to the estimates of g 0:t . Figure 1 shows the proportion of observations with truncated g 0:2 as a function of β p at a null effect, i.e. β ψ0 = 0. βp Proportion q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Under these settings, the true parameter values ψ 0 were achieved by generating 8 × 10 7 observations under the counterfactual distribution for each β ψ0 considered. Simulation results were obtained for 500 simulations of size n = 500. Within each simulation, the bootstrap estimates of variance were formed from B = 1000 iterations.

Submodels used
Any submodel and loss function for which its loss-function specific score spans D * (P 0 ) can be chosen in TMLE for both estimation of the mean outcome EY a and the variance of the EIF σ 2 . As these submodels solve the equation corresponding to the EIF, they will all be asymptotically equivalent and thus, all be asymptotically efficient. That is, no difference will be seen between the use of various submodels as the sample size grows to infinity. In the TMLE presented by van der Laan et al. [van der Laan and Gruber, 2011], these submodels are used in the targeting step for eachQ t using a loss L(Q t ) that is indexed byQ t+1 . Specifically, for the targeting step we need a loss and submodel with clever covariate such that the score given solves a desired component of the efficient influence function D * t (P 0 ). The use of various submodels in finite samples can have varying performance. For example, under increasing levels of positivity violations the use of linear submodels which use H t (g) as a covariate can have higher variance due to observations with low probabilities of treatment acting as outliers which result in highly influential points for the estimation of the submodel parameter .
Recall that the catalyst for this work was the anti-conservative estimates of estimator variance resulting from the use of the empirical EIF variance. We therefore wish to establish a robust estimator of the variance of estimators which solve the EIF, particularly under violations or near violations of positivity. In other words, we desire a variance estimator which will asymptotically converge to the true variance of the estimator, but also simultaneously act on the conservative side in finite samples. Keeping this in mind, we used two submodel and loss function combinations for our simulations. For the estimation of the target parameter and the robust estimator of the EIF variance, we used submodels which define H t (g) and H d,t m (g) to be observational weights such that logitQ( ) = logitQ + , acknowledging our slight abuse of notation. Alternatively, in our bootstrap approach at estimating the TMLE variance, we define a clever covariate using H t (g) such that logitQ( ) = logitQ + H t (g).
Both submodels use, as loss function, the negative log-likelihood loss. As stated previously, both of the submodels presented solve the equation corresponding to the EIF and are therefore asymptotically equivalent. Figure 2 shows the Monte-Carlo variance under no treatment effect (β ψ0 = 0) for both the AIPW and TMLE estimators, along with the mean of the variance estimates from each estimation approach. To keep the differences in perspective, we plotted results only for the positivity associated parameter β p ≤ 0. At the lower end of β p where positivity violations are minor, the observed estimator variance is low for both the AIPW and TMLE estimators, with the TMLE approach showing lower variance between the two despite solving the same estimating equation corresponding to the EIF. For example, at β p = −2 the Monte-Carlo variance was 0.0024 and 0.0022 for the AIPW and TMLE estimators, respectively. As β p increased, introducing higher levels of positivity violations, the estimator variance increased for both estimators. Additionally, this occurred at a much higher rate for the AIPW estimator than for TMLE, resulting in an increase in the magnitude of difference between the two estimators. For example, at β p = 0 the simulations resulted in a Monte-Carlo variance of 0.0207 and 0.0085 for the AIPW and TMLE estimators, respectively. For the AIPW estimator, the empirical EIF based approach of estimating variance performed resulted in estimates similar to the Monte Carlo estimates. For example, at β p = 0 the mean of the EIF approach was 0.0208. A slight but consistent underestimation of the variance was observed at higher levels of practical positivity violations. The robust approach of estimating variance appeared to result in conservative estimates of variance.

Point treatment results
In the TMLE estimator, all three approaches to variance estimation performed similarly at low values of β p . For example, at β p = −2, the mean of the estimates was 0.0029, 0.0029, and 0.0032 for the empirical EIF, robust, and bootstrapped based approaches, respectively, compared to the estimator's Monte-Carlo variance of 0.0022. As β p increased, the empirical EIF approach tended to result in anti-conservative estimates of variance, while the bootstrap approach resulted in slightly conservative estimates. The robust EIF approach tended to overestimate the TMLE estimator variance. Figure 3 shows the Monte-Carlo variance for each approach taken at estimating the variance. Lower values in this figure can be interpreted as coming from a variance estimator with more precision. In the AIPW estimator, the empirical EIF approach has noticeably higher variance than the robust approach, with a variance of 2.93 at β p = 0 compared with 1.50 for the robust approach. This implies that the empirical EIF approach to estimating the AIPW estimator variance is less reliable than the robust EIF approach. For the TMLE estimator (Figure 3b), the empirical EIF approach to estimating variance showed much lower Monte-Carlo variance. The bootstrap approach also resulted in very low variance, implying a high reliability of this approach at estimating the variance. We evaluated 95% confidence interval coverage for the TMLE estimator of EY d under the three approaches to variance estimation. Due to the lower variance seen in Figures 2 and 3, we focused only on the TMLE estimator here. Figure  4 shows a heat map overlaid with a contour plot of the resulting coverage estimates (i.e. the observed proportion of times the true parameters were captured by the confidence intervals) over the different combinations of β ψ0 and β p . Additionally, we estimated the power to reject the null hypothesis (at a level of 0.05) corresponding to each variance estimation approach under the range of treatment effect sizes and degrees of positivity violation considered above. Figure 5 shows a heat map overlaid with a contour plot of the resulting power estimates. Results at β ψ0 = 0 can be interpreted as Type I errors, as they inform us of the times that the null hypothesis of no treatment effect is incorrectly rejected.
At low instances of positivity issues, coverage appears valid for all three variance estimation approaches with the proportion of time the true parameter was captured consistently at 0.95 or larger (Figure 4). Where positivity issues were low (β ψ0 < −0.5), the empirical EIF approach maintained nominal to conservative coverage. Where severe positivity violations were present, coverage dropped substantially below 0.95. For example, at β p = 1 coverage for this approach varied from 0.41 to 0.85. In contrast, the robust EIF approach consistently resulted in coverage at around 0.95 − 0.96 at low values of β p and increased with β p , consistent with prior results showing overestimation of the variance under increasing positivity by this approach. For example, at β p = −.7, coverage remained at 0.98 at all values of β ψ0 . At β p ≥ −0.1, the observed coverage was almost always greater than or equal to 0.99 at all values of β ψ0 . The bootstrap based coverage shown in Figure 4c varied the least, with coverage consistently between 0.95 − 0.97 irrespective of the treatment effect (β ψ0 ) and positivity severity (β p ) considered. Regarding the observed power ( Figure 5), the empirical-EIF based variance approach resulted in the highest power among all three variance estimation approaches when an effect was present. For example, at β ψ0 = 1 and β p = −1, the observed power was 0.71, 0.51, and 0.51 for the empirical-EIF, robust-EIF, and bootstrap approaches respectively. While this result implies a more efficient approach, it expectedly came at a cost of higher Type I error which became uncontrolled with an increase in β p . For example, at β p = −2 an observed 4.2% of the simulations incorrectly rejected the null hypothesis. This proportion increased to as high as 15% at β p = 1. Alternatively, the robust EIF estimation approach resulted in low Type I errors (i.e. between 0 − 5.8%) with none of the simulations incorrectly rejecting the null beyond β p = −0.1. The bootstrap approach resulted in an intermediate performance, with higher power than the robust EIF approach when an effect was present while simultaneously retaining appropriate control of the Type I error at all levels of β p when no effect was present. For example, at β p = 1 only 4.8% of the simulations incorrectly rejected the null hypothesis.

Longitudinal treatment results
Results for the longitudinal setting were less stable, though still similar to the point treatment setting. Figure 6 shows the mean of the variance estimates under each approach, overlaid with the Monte-Carlo variance of the intervention specific mean outcome estimators. The same trend over the different levels of positivity was seen as in Figure 3, with the variance increasing with the magnitude of positivity issues. The empirical EIF approach also performed well here at low levels of β p for both the AIPW and TMLE estimators. At high values of β p , the approach more noticeably underestimate the variance of both intervention specific mean outcome estimators. Consistent with the point treatment setting, the robust EIF approach consistently over estimated the variance for both estimators. The bootstrap approach resulted in slightly conservative variance, though were still very similar to the Monte-Carlo variance estimates.  Figure 7 shows the coverage corresponding to each variance estimation approach for the TMLE estimator of the intervention specific mean outcome. Coverage for the empirical EIF approach dropped considerably with an increase in positivity issues. For example, at a null effect (i.e. β ψ0 ) the observed coverage was 0.93 at β p = −2 and 0.78 at β p = 0. For the robust EIF approach, coverage increased with positivity. This became as high as 1.00 (i.e. all simulated confidence intervals captured the true parameter value) at higher levels of positivity issues. For the bootstrap approach, a higher level of coverage was also seen. For example, under a null effect, a coverage of 0.95 was observed at β p = −2 and 0.98 at β p = 0. Results for the Type I error and power were also similar to the point treatment setting. When there was an effect, the empirical EIF approach resulted in the highest power. At β ψ0 = 1 and β p = −2, we observed a power of 0.99. However, the Type I error was also uncontrolled here, becoming as high as 0.22 at β p = 0. While the robust EIF approach maintained valid Type I error rates, the power for this approach when an effect was present was the lowest. For example, for an treatment effect size of β ψ0 = 1 we observed a power ranging from 0.996 at β p = −2 to 0.14 at β p = 0. The bootstrap approach also resulted in controlled Type I error rates, with observed values below 0.05 over all values of β p considered. Power was higher than the robust EIF approach across all values of β ψ0 and β p . For a treatment effect size of β ψ0 = 1, we observed a power ranging from 0.988 at β p = −2 to 0.40 at β p = 0 for the bootstrap approach. Compared with the robust EIF approach, this is almost a 3-fold increase in power.

Discussion
The goal of the current study was to establish a consistent and robust approach of estimating the variance of asymptotically efficient estimators such as TMLE, estimating equations, and one step estimators which, in contrast to the common approach based on the empirical variance of the estimated EIF, do not act anti-conservatively when confronted with positivity violations. We have presented two such approaches at estimating this variance: 1) a robust approach that directly targets the asymptotic variance of the EIF, and 2) a bootstrap approach based on fitting the initial density of the data once, followed by a non-parametric bootstrap of the targeting step. In simulations, the variance of AIPW increases noticeably with the variance of the EIF as positivity increases. The variance of the TMLE was constrained in the face of increasing positivity violations, and as a result, while the empirical EIF approach underestimated variance, the robust EIF approach increasingly over-estimated the variance as the degree of positivity violations increased. In contrast, the bootstrap based approach provided less conservative variance estimation, while maintaining valid Type I error control in the face of extreme positivity violations, both in the point treatment and longitudinal setting.
We emphasize that, as the robust EIF approach directly targets the asymptotic variance of the EIF, extremely large values of estimates from this can be used to raise a red flag for unreliable statistical inference due to sparsity, thereby declaring that the target parameter is practically not identifiable from the data, and that the reported variance estimates (though large) will themselves be imprecise. As such, we recommend that this approach be used if there is concern regarding identifiability of the data for the target parameter of interest.
While the EIF can raise a red flag for lack of identifiability, for substitution estimators such as TMLE we suggest that it is overly conservative for constructing valid confidence intervals and tests in finite sample in the face of substantial positivity violations. In previous work [Petersen et al., 2012], we suggested the parametric bootstrap as a diagnostic tool for sparsity-bias in the point treatment setting. The approach can become cumbersome, as the analyst would need to refit the whole likelihood for each iteration of the bootstrap. Our robust EIF approach is able to avoid estimating the whole likelihood by targeting the required means under the post intervention distribution defined by the longitudinal g-computation formula directly. Even if we use an actual TMLE of P d 0 , our analytic estimate of the variance is still much less computer intensive than the parametric bootstrap method, in particular, in view that one would need to run many replicate samples of the data set in order to pick up the observations that correspond with a rare treatment and thus contribute large numbers to the variance expression. Therefore, we believe that the proposed analytic method will be (at least, practically) superior to the earlier proposed parametric bootstrap method. Our presented bootstrap approach, while more computationally intensive than the robust EIF approach, is also superior to the earlier proposed approach, in that we do not have to refit the entire likelihood under each iteration. This also significantly reduces the computational resources required to obtain estimates of the target parameter, particularly if computationally intensive non-parametric machine learning algorithms are used to estimate these densities.
Further refinements can be applied in an attempt to obtain variance estimates with an even smaller finite sample bias. One such approach is a convex combination of the variance estimators considered above. For example, we noticed in supplementary analyses that takingα nσ 2 eEIF,n + (1 −α n )σ 2 rEIF,n had good performance, whereσ 2 eEIF,n is the variance estimate using the empirical EIF approach,σ 2 rEIF,n is the variance estimate using the robust EIF approach, andα n = |σ 2 rEIF,n −σ 2 eEIF,n |/(σ 2 rEIF,n +σ 2 eEIF,n ). We note, however, that such an approach is somewhat ad-hoc and may lead to varying results in other simulations or distributions. We therefore chose not to present the results here.
A potential limitation of the robust approach at estimating the variance involves the conditions for asymptotic linearity to be met. Note, however, that we always require that our estimator (of the parameter) be asymptotically linear. Thus, this is actually a limitation of our parameter estimator. Given that we have an asymptotically linear estimator, we want a good estimator of its variance. Furthermore, it is also required thatQ d,σ 2 t 0 be estimated both consistently and at a fast enough rate. We limited the computational complexity in our simulations by using simpler parametric models to estimateQ d,σ 2 t 0 , though a more non-parametric approach such as Super Learning could have been applied. This approach can become computationally expensive if there are many time points. In this regard, the bootstrap approach is superior as it does not require the additional estimation ofQ d,σ 2 t 0 .
It would be of interest to further evaluate not only the practical performance of these variance estimation approaches in future studies, but also the application of the approaches to other parameters. The appendix derives the general approach for working marginal structural models. Further research into the practical performance of this approach is needed for this setting. These variance estimation approaches can also be used for more complex parameters, such as mean outcomes under dynamic regimes, stochastic interventions, or marginal structural working models. Future research could also develop a collaborative TMLE [van der Laan and Gruber, 2010] or cross-validated [Zheng and van der Laan, 2010] based approach at robustly estimating the EIF variance.
It is important to note that in expressing our variance this way, we integrate out the indicator of treatment overĀ, i.e. I(Ā = d(L)). By getting rid of this indicator, we no longer rely as heavily on observations from subjects following treatment in estimating the variance of D * K+1 . This particularly helps us when there is a lack of positivity, since subjects with low probabilities of desired treatment simply are not observed.
We can therefore also estimate this variance component σ 2 K+1 by computing the TMLE of the intercept β 0 in the working MSM for our newly defined outcome Z(K + 1) using weights h 1 (d, K + 1), and then multiplying it against d∈D h 1 (d, K + 1).
Let's now define the counterfactual outcomeZ and the corresponding observed outcomeZ We could apply the TMLE to estimate EZ d (d) based on the longitudinal data structure (L(0), A(0), . . . , L(K), A(K),Z(d, K + 1)), for each d ∈ D, and use that In applying TMLE here, we should be using that Before we go to the next conditional expectation we need to target with a parametric submodel, such as In this way, we will only have to run one TMLE for each rule d, which still utilizes that the outcome is a sum of outcomes that are known for histories including that outcome.