Estimating Effects with Rare Outcomes and High Dimensional Covariates: Knowledge is Power

1 Introduction

When the outcome of interest occurs infrequently, effect estimation can be particularly challenging. For example, a recent study sought to examine the impact of planned place of delivery (obstetric unit or not) on perinatal mortality and neonatal morbidities, occurring in 250 of 63,827 births (0.39 %) (Birth-place in England Collaborative Group 2011). Due to the paucity of individual birth events, however, the researchers estimated the effect on a composite outcome measure. Likewise, tuberculosis is a main cause of mortality among HIV+ people (World Health Organization 2013). Evaluating strategies to reduce its transmission are essential, but difficult due to the disease’s relatively low incidence. Along the same lines, international consortiums have been established to investigate the burden and treatment for uncommon cancers (e. g. RARECARENet 2014). While these outcomes are rare, a better understanding of their occurrence is likely to have important policy and health implications.

For binary outcomes or proportions, parametric logistic regression is often used to estimate the conditional odds ratio, given the exposure and measured covariates. Several researchers have investigated the performance of this approach when the outcome is extremely rare (e. g. Concato et al. 1993; Harrell et al. 1996; Peduzzi et al. 1996; King and Zeng 2001; Harrell 2001; Braitman and Rosenbaum 2002; Cepeda et al. 2003; Vittinghoff and McCulloch 2007). For example, simulations by Peduzzi et al. (1996) illustrated that estimates could be biased and inference unreliable if the number of outcomes per independent variable in the regression model was less than 10. The authors also found problems with estimator convergence, statistical power and the validity of significance tests (i. e. type I error rates and confidence interval coverage). Harrell et al. (1996) cautioned against over-fitting and encouraged the use of cross-validation or bootstrapping for model validation. Moreover, King and Zeng (2001) found that standard logistic models could substantially under-estimate the probability of the outcome and offered a bias correction with accompanying software.

When dealing with rare events, several researchers have recommended estimators based on the propensity score, which is the conditional probability of being exposed, given the covariates (Rosenbaum and Rubin 1983). These methods avoid estimation of the conditional mean outcome and thereby are expected to perform well when there are very few outcome events (e. g. Joffe and Rosenbaum 1999; Braitman and Rosenbaum 2002; Patorno et al. 2014). Simulations by Cepeda et al. (2003) suggested that propensity score methods were less biased and more efficient than logistic regression for the mean outcome, when the number of events per independent variable in the regression model was less than 8. The authors also cautioned that the performance of propensity score methods depended on the strength of the relationship between the covariates and the exposure.

Targeted minimum loss-based estimation (TMLE) is a general methodology for the construction of semiparametric, efficient substitution estimators (van der Laan and Rubin 2006; van der Laan and Rose 2011). A TMLE for a single time point exposure can be implemented as follows. First, the conditional expectation of the outcome, given the exposure and covariates, is estimated with parametric regression or with a more flexible approach, such as SuperLearner (van der Laan et al. 2007). Second, information on the exposure-covariate relation (i. e. the propensity score) is incorporated to improve this initial estimator. The propensity score can also be estimated parametrically or with a more flexible approach. Informally, this “targeting” step helps remove some of the residual bias due to incomplete adjustment for confounding. More formally, this targeting step serves to solve the efficient score equation. Finally, the targeted predictions of the outcome under the exposure and under no exposure are averaged over the sample and contrasted on the relevant scale.

Thereby, TMLE requires estimation of both the conditional mean outcome as well as the propensity score, and achieves a number of desirable asymptotic properties (van der Laan and Rose; 2011). The standardized estimator is asymptotically normal with mean 0 and variance given by the variance of its influence curve. The TMLE is also double robust: if either the conditional mean outcome or the propensity score is consistently estimated, we will have a consistent estimate of the parameter of interest. If both functions are consistently estimated (at a fast enough rate), the TMLE will be efficient and achieve the lowest possible asymptotic variance among a large class of estimators. Finally, TMLE is a substitution estimator, providing robustness in the face of positivity violations (when there is no or little variability in the exposure within certain covariate strata) and rare outcomes (e. g. Stitelman and van der Laan 2010; Gruber and van der Laan 2010; Sekhon et al. 2011; Petersen et al. 2012; Gruber and van der Laan 2013; Lendle et al. 2013). Building on the work of Gruber and van der Laan (2010), this paper proposes a new TMLE for the semiparametric statistical model m, which imposes bounds on the conditional mean of the outcome, given the exposure and measured confounders. We focus our discussion on rare binary outcomes and rare bounded continuous outcomes (e. g. proportions). The estimation problem and the theoretical motivation for the new TMLE are outlined in Section 2. In particular, the causal parameter, the corresponding statistical parameter, the statistical model and the efficient influence curve are discussed. Section 3 presents the rationale and procedure for the rare outcomes TMLE. Simulations and the applied analysis are given in the Section 4. The article concludes with a discussion of the advantages and disadvantages of the proposed method as well as areas for future work.

2 The estimation problem

We are interested in estimating the impact of a binary exposure A on the risk of a rare outcome Y. For example, Y might be an indicator that the subject develops tuberculosis with an incidence rate of 255/100,000 per person-year in Sub-Saharan Africa (World Health Organization 2013). Alternatively, Y might be the one-year cumulative incidence of tuberculosis for a given community. In the latter scenario, the outcome Y is a proportion bounded between [0,1]. Suppose we measure some baseline characteristics W that are predictors of both the exposure and outcome. In other words, W represents the set of measured confounders. Let O=(W,A,Y) denote the observed data random variable with distribution P0. Throughout, subscript 0 denotes the true, but unknown distribution.

To translate our scientific question into a causal parameter, let us define Ya as the counterfactual outcome, if possibly contrary-to-fact, the unit received exposure A=a (Neyman 1923; Rubin 1974; Pearl 2000). We assume these quantities exist for all units both under the exposure (A=1) and under no exposure (A=0). To relate the observed outcomes to the counterfactual outcomes, we need the stable unit treatment value assumption (SUTVA) (Rubin 1978; 1980): (1) the counterfactual outcomes for one unit must not be impacted by the treatment assignment of another unit (i. e. no interference), and (2) there must not be multiple versions of the treatment A=a. Under this assumption, we have

In words, we only observe the counterfactual outcome corresponding to the observed treatment YA=Y. Thereby, the observed data can be considered as a time-ordered missing data structure on the full data XF=(W,Y1,Y0)∼PX, with the exposure A as the censoring variable (Neyman 1923; Rubin 1974). Throughout our goal is to estimate and obtain inference for the population average treatment effect:

This causal parameter is the difference in expected counterfactual outcomes if everyone in the target population were exposed and if everyone in that target population were unexposed. For a binary outcome, ΨF(PX) corresponds to the causal risk difference: P(Y1=1)−P(Y0=1).

To express ΨF(PX) as a function of the observed data distribution P0, we need several assumptions. First, there must be no unmeasured confounders of the effect of the exposure on the outcome (Rosenbaum and Rubin 1983; Robins 1986). Secondly, there must be sufficient variability in the treatment assignment. In other words, the propensity score P0(A=1|W) must be bounded away from 0 and 1. This condition is known as the positivity assumption. Under these assumptions, we can express the causal parameter in terms of the difference in the conditional mean outcomes, averaged (standardized) with respect to the covariate distribution (Robins 1986):

where the summation generalizes to the integral for continuous covariates and Qˉ0(A,W)=E0(Y|A,W) denotes the conditional mean outcome, given the exposure and covariates. For a binary outcome Ψ(P0) is sometimes called the marginal risk difference.

The challenge, addressed in this paper, is estimation of Ψ(P0) in the context of extremely rare outcomes and high dimensional covariates W. This challenge is illuminated by studying the efficient influence curve (function) of the target parameter Ψ at the true probability distribution P0 (Bickel et al. 1993; van der Laan and Rose 2011):

where I(⋅) is the indicator function. D∗(P0)(O) is the critical ingredient in the construction of double-robust, semiparametric, efficient estimators (van der Laan and Robins 2003; van der Laan and Rose 2011). Specifically, D∗(P0)(O) can be represented as an estimating function in Ψ(P0) as in augmented inverse probability of treatment weighting (A-IPTW) or used to fluctuate initial estimates in the substitution estimator TMLE.

The information for learning the target parameter is captured by the sample size n divided by the variance of the efficient influence curve at the true distribution P0:

(van de Vaart 1998). Low information can occur when the propensity score P0(A=1|W) is or approaches 0 or 1 for certain treatment-covariate combinations. The impact of violations of the positivity assumption on estimator performance are demonstrated in Petersen et al. (2012), among others. Low information can also occur when there are very few outcome events (i. e. E0(Y) is small) and the conditional mean Qˉ0(A,W) is or approaches 1 for some treatment-covariate combinations. In either case, the variance of the efficient influence curve can become large, and the sample size needed to detect effects with sufficient power correspondingly large. In this paper, we address estimation in the context of low information due to rare outcomes and high dimensional covariates.

We consider the semiparametric statistical model m, which incorporates knowledge that the conditional mean of the rare outcome Qˉ0(A,W) is bounded from above by some u<1. For rare binary outcomes, this knowledge has been noted by other researchers (Beck, King and Zeng, 2000, King and Zeng, 2001). Consider, for example, an unusual outcome with incidence 0.5 %. With this rare of an event, it seems unlikely that the conditional mean outcome exceeds 10% for any combination of the exposure and measured confounders. Instead, the researcher might have knowledge that this conditional mean does not exceed 5% for the study population. Formally, the statistical model m is the set of possible data generating distributions:

for some bounds 0≤ℓ<u<1. The specification of the bounds can be based on subject matter knowledge or selected with cross-validation, as discussed in Section 3.3. We place no restrictions on the marginal distribution of baseline covariates P(W) or on the propensity score P(A=1|W).

The variance of the efficient influence curve (Eq. [1]) and the information for estimating the target parameter (Eq. [2]) depend on whether the true distribution P0 and thereby the true conditional mean Qˉ0(A,W) satisfy the constraints in the semiparametric model m. Likewise, the asymptotic variance of a TMLE (under sampling from P0) depends on whether the limit of the estimator of Qˉ0(A,W) satisfies these constraints. Specifically, the variance of the efficient influence curve and the variance of the influence curve of the TMLE will be smaller when the true conditional mean and the limit of the estimator are in [ℓ,u] for some u<1 (Appendix A). Furthermore, when the outcome is rare and the dimension of the covariates large, we are often forced to use a misspecified estimator for Qˉ0(A,W). Then the influence curve of the TMLE is defined by the efficient influence curve at the misspecified limit Qˉ(A,W). The constraints in m heavily affect the variance of this influence curve and therefore motivate us to develop an estimator that respects the known bounds (Appendix B).

In the following section, we propose a new TMLE, which guarantees that the initial and targeted estimator of Qˉ0(A,W) are constrained within the model bounds. We compare the proposed modification to the standard TMLE algorithm for binary or bounded continuous outcomes (Gruber and van der Laan 2010). The two TMLEs will be asymptotically equivalent if both use a consistent estimator of the conditional mean outcome Qˉ0(A,W) or more generally if their estimators of Qˉ0(A,W) converge to the same limit. In finite samples, however, enforcing the constraints is expected to have substantial impact on performance. Indeed, the finite sample behavior of a TMLE is largely driven by the sample variance of the influence curve at the targeted estimator Qˉn∗(A,W), and this variance is expected to be much smaller when Qˉn∗(A,W) respects the model constraints. As a result, we expect and show in simulations the proposed method will result in tighter confidence intervals and more power in finite samples.

Before introducing the new TMLE, we first caution that the statistical model m must contain the true observed data distribution P0. In other words, the assumed bounds on the conditional mean outcome Qˉ0(A,W) should be given careful consideration and not chosen based on convenience. For example, there might be a scenario, where the marginal expectation of the outcome Y is small (e. g. E0(Y)<.01), but the conditional mean Qˉ0(A,W) approaches 1 for some combinations of the exposure and measured confounders. Consider, for example, a study of the association between very preterm birth (28 to <32 weeks) and neonatal mortality. If the set of measured confounders included socio-economic factors, race, maternal smoking and maternal anemia, then the conditional risk of the outcome could approach 1 for some infants. As a result, it would be inappropriate to assume a semiparametric statistical model m with upper bound u<1. Enforcing a small upper bound may result in bias and lower power in finite samples. Instead, for this example, an upper bound of u=1 would be appropriate, and the new estimator, proposed in the following section, would reduce to the standard TMLE procedure for binary or bounded continuous outcome. For the remainder of the article, we assume that the semiparametric statistical model m characterizes the set of possible observed data distributions in that the true conditional mean outcome Qˉ0(A,W) is bounded in [ℓ,u] for some 0≤ℓ<u<1.

3 Estimating effects with rare outcomes

Suppose we have n observations O=(W,A,Y) drawn independently from P0. A simple substitution estimator (i. e. the G-Computation estimator of Robins (1986)) for Ψ(P0) can be implemented with the following steps. First, we obtain an estimate of the conditional mean outcome Qˉ0(A,W). Let us denote an initial estimator based on n observations as Qˉn(A,W). Then we obtain the predicted outcomes for each observation under the exposure Qˉn(1,W) and under the control Qˉn(0,W). The sample average difference in these predicted outcomes provides a point estimate for Ψ(P0). The final step corresponds to estimating the marginal covariate distribution with the empirical proportion.

In some cases, we can estimate Qˉ0(A,W) with the sample average in each exposure-covariate strata. This non-parametric estimator can quickly become ill-defined when there are many exposure-covariate combinations and can suffer from over-fitting, especially with rare outcomes. In some cases, we may have the background knowledge to support parametric regression for Qˉ0(A,W). Suppose, for example, the true conditional mean outcome could be described by the following parametric regression model.

where logit(x)=log(x/1−x) and W=(W1,…,Wp) denotes the set of confounders. Asymptotically, the estimated coefficients will yield predictions of the mean outcome, which respect the model knowledge, and estimates of the target parameter, which are consistent and efficient. However, with small sample sizes and many confounders, estimation will rely heavily on extrapolation and may yield estimates outside of the model bounds. In this realistic situation, incorporating knowledge of the upper bound u into the estimation algorithm is expected to improve performance. In particular, we conjecture that even if the parametric regression is misspecified, then using the bound u will result in more precision and less bias. Nonetheless, we rarely have the knowledge to specify the form of Qˉ0(A,W) with a finite number of parameters (i. e. β). As a result, we advocate for a more flexible approach that respects our semiparametric statistical model m.

Data-adaptive or machine learning algorithms can incorporate our limited model knowledge, while smoothing over areas of data with weak support. SuperLearner, for example, allows a set of pre-specified algorithms to compete and selects the best algorithm using cross-validation (i. e. data-splitting) (van der Laan et al. 2007). The selected algorithm can be used to predict the outcomes for all units given the exposure and their measured covariates Qˉn(1,W) as well as given no exposure and their measured covariates Qˉn(0,W). The average difference in these predicted outcomes provides a point estimate of the target parameter. Inference, however, must respect the model building process. (Treating the selected algorithm as if it were a priori-specified can result in misleading inference.) Furthermore, estimating the conditional expectation Qˉ0(A,W) is a more ambitious task than estimating a single number Ψ(P0)∈R. Thereby, an estimator of Qˉ0(A,W) will have the wrong bias-variance trade-off for the parameter of interest Ψ(P0). Specifically, the estimator will have too much bias relative to its variance.

TMLE provides a solution to several of these challenges (van der Laan and Rose 2011). TMLE incorporates data-adaptive estimation with an additional targeting step to reduce bias for Ψ(P0) and to attain valid statistical inference. As described in the introduction, a TMLE for Ψ(P0) can be implemented in three steps. First, we estimate of the conditional mean outcome Qˉ0(A,W) with parametric or more flexible methods. Second, we target this initial estimator using information in the estimated propensity score Pn(A=1|W). Finally, we substitute the targeted estimates into the parameter mapping. The updating step is accomplished with a loss function and a carefully selected submodel. For a binary or bounded continuous outcome, we can use the negative log-likelihood loss function (Wedderburn 1974; McCullagh 1983; Gruber and van der Laan 2010):

and the following logistic regression model to update the initial estimator:

Here, ε is the univariate parameter and Pn(A=1|W) denotes an estimate of the propensity score. This combination of loss function and fluctuation model returns the initial estimator Qˉn(A,W) at ε=0 and has score equal to the relevant component of the efficient influence curve at ε=0 (Gruber and van der Laan 2010). To estimate the fluctuation parameter ε, we run logistic regression of the outcome Y on the covariate Hn(A,W) with the logit of the initial estimates as offset. Plugging in the estimated coefficient εn yields the targeted updates. A point estimate can then be obtained by taking the average difference in the targeted predictions under the exposure Qˉn∗(1,W) and under no exposure Qˉn∗(0,W).

The standard TMLE algorithm for a binary or bounded continuous outcome is not optimal for our semiparametric statistical model m. Specifically, this TMLE does not enforce the global constraints on the conditional mean outcome. First, initial estimates Qˉn(A,W) based on the log-likelihood loss function (Eq. [3]) are guaranteed to be within (0,1), but are possibly outside of the model bounds [ℓ,u]. Second, the logistic regression, used to update the initial estimator, also does not respect the constraints implied by m. Therefore, the updated estimates Qˉn∗(A,W) are guaranteed to be between (0,1), but may be outside of the model bounds [ℓ,u]. As a result, the asymptotic and finite sample performance of a TMLE, using this loss function and parametric working model, is expected to be sub-optimal. The algorithm is expected to be unstable (i. e. more variable and less accurate) when the outcome is very rare and there are many confounders.

4 Simulations & data application

The following simulation studies compare the finite sample performance of the standard TMLE with the proposed rare outcomes TMLE (rTMLE). For comparison, we also include a propensity score matching (PSM) estimator (Rosen-baum and Rubin 1983), inverse probability of treatment weighted (IPTW) estimator (Hernán et al. 2000) and augmented inverse probability of treatment weighted (AIPTW) estimator (Robins 2000; van der Laan and Robins 2003). The first two methods rely solely on estimation of the propensity score and may have superior performance with very rare outcomes. AIPTW requires estimation of both the conditional mean outcome and the propensity score, but is double robust and asymptotically efficient under consistent estimation of both functions. AIPTW is an estimating equation (i. e. not a substitution estimator) and thereby can result in impossible parameter estimates (e. g. probabilities less than 0 or greater than 1) (Lendle et al. 2013).

For the PSM estimator, we used the Matching package (Sekhon 2011) for 1:1 matching based on the estimated propensity score and calculated the point estimate as

with YMi denoting the outcome of the observation matched to unit i based on the estimated propensity scores. A point estimate from IPTW is given by the following weighted mean

A point estimate from AIPTW is attained by directly solving the efficient score equation:

where Qˉn(A,W) denotes a non-targeted estimate of the conditional mean outcome. Inference was based on the Abadie-Imbens standard error estimator for the PSM estimator (Abadie and Imbens 2006; Sekhon 2011) and the estimated influence curve for the others. Wald-type confidence intervals were calculated as ψn±1.96sen, where ψn denotes the point estimate and sen denotes the estimated standard error. Likewise, tests of the null hypothesis of no effect were based on test statistic Tn=ψn/sen.

4.1 Simulation 1: individual-level data

The finite sample performance of the estimators was assessed by drawing 2,000 samples of sizes n=1000 and n=2500 according to the following process. First, we generated three baseline covariates:

W1∼Normal(0,0.252),W2∼Uniform(0,1),W3∼Bernoulli(0.5).

The exposure A was drawn from a Bernoulli distribution with probability

P0(A=1|W)=logit−1[−0.5+W1+W2+W3].

With this exposure mechanism, there were no positivity violations; the propensity score was bounded between 20 % and 92 %. Finally, the binary outcome was drawn from a Bernoulli distribution with probability

P0(Y=1|A,W)=logit−1[−3+2∗A+W1+2∗W2−4∗W3+0.5∗A∗W1]/15.

The resulting marginal probability of the outcome was E0(Y)=1.1%. The true bounds on the conditional mean Qˉ0(A,W) were [0%,6.2%], and the true value of the statistical parameter was ψ0=1.3 %. The propensity score P0(A=1|W) was estimated with the correctly specified main terms logistic model and a misspecified regression model failing to adjust for W3. The conditional mean outcome Qˉ0(A,W) was estimated with the correctly specified logistic regression model as well as a misspecified regression with main terms for only A and W1. For the rare outcomes TMLE, the lower bound was set 0 and the upper bound was selected data-adaptively from the set {2.5 %, 5 %, 7.5 %, 10 %} using cross-validation with the log-likelihood loss function.

The results for this simulation are presented in Figures 1 and 2. As expected, the PSM estimator and IPTW exhibited low bias when the regression model for the propensity score P0(A=1|W) was correctly specified, but were biased otherwise. This bias did not disappear with sample size. Likewise, these estimators had poor confidence interval coverage under misspecification of the propensity score. Even under correct specification of this function, the PSM estimator had notably lower power than IPTW, AIPTW or rTMLE. Indeed, neither this PSM estimator nor this IPTW are efficient estimators. (The reader is referred to Abadie and Imbens (2015) for recent work on consistent estimation of the standard error when estimating the propensity score.)

Figure 1:

Boxplots of the point estimates from various algorithms for sample sizes n=1000 (left) and n=2500 (right) for Simulation 1. CC indicates both the outcome regression and the propensity score are correctly specified. CM indicates the outcome regression is correctly specified, but propensity score misspecified. MC indicates the outcome regression is misspecified, but the propensity score correctly specified. The dashed line indicates the true value ψ0=1.3%.

Figure 2:

Attained confidence interval coverage and power for the estimators in Simulation 1. The unfilled and filled points denote sample sizes of n=1000 and n=2500, respectively. Circles indicate both the outcome regression and the propensity score are correctly specified. Triangles indicate the outcome regression is correctly specified, but propensity score misspecified. Squares indicate the outcome regression is misspecified, and propensity score correctly specified. The dashed lines indicate 95 % confidence interval coverage and 80 % power, respectively.

In theory, the other estimators solve the efficient score equation directly (AIPTW) or during the targeting step (TMLE, rTMLE). As a result, these estimators are double robust. For these simulations, however, the standard TMLE exhibited substantial bias when estimating Qˉ0(A,W) with the correctly specified regression model, which included four main terms plus an interaction (Figure 1). Due to the paucity of outcome events, logistic regression of the untransformed outcome Y was unstable. Thereby, the updating step, which involved fitting an additional coefficient ε, did little to reduce bias. As a result, its confidence interval coverage was poor under correct specification of the outcome regression (Figure 2). On the other hand, when the initial estimator of Qˉ0(A,W) was based on the misspecified regression, the performance of the standard TMLE was comparable to AIPTW and rTMLE. Recall the misspecified regression model for Qˉ0(A,W) only included main terms for the exposure A and confounder W1. The smaller adjustment set led to more stable initial estimates, but also complete reliance on consistent estimation of the propensity score for confounding control and consistency. Thereby, a potential approach when estimating effects with rare outcomes is to use a smaller adjustment set for Qˉ0(A,W) and thereby sacrifice double robustness and efficiency for greater stability and potentially reduced bias. This approach was taken by Gruber and van der Laan (2013).

The performance of AIPTW and rTMLE did not suffer when fitting Qˉ0(A,W) in the larger, correctly specified regression model. As shown in Figure 1, both AIPTW and rTMLE had low bias when either the outcome regression, the propensity score or both were correctly specified. Furthermore, these estimators had good confidence interval coverage at both sample sizes and attained power for the larger sample (Figure 2). Therefore, it seemed that the performance of AIPTW and rTMLE were comparable for this set of simulations. It is worth emphasizing, however, that rTMLE is a substitution estimator and thereby guaranteed to respect the global constraints in the model. In particular, AIPTW yielded negative (impossible) estimates for the marginal risk under no exposure E0[E0(Y|A=0,W)](Table 3 in the Appendix).

Under the null, the PSM estimator and IPTW had good type I error control when the propensity was correctly specified (Results not shown). Under misspecification, their type I error rates exceeded 20 %. The standard TMLE also suffered from high type I error rates, when both regression models were correctly specified. AIPTW and rTMLE maintained nominal type I error rates for both sample sizes and all regression specifications.

5 Discussion

In this paper, we proposed a new TMLE for evaluating causal effects and estimating associations with very rare outcomes and high dimensional data. The rare outcomes TMLE (rTMLE) is based on harnessing knowledge in the semiparametric model m, which bounds the conditional mean outcome Qˉ0(A,W) from below by ℓ≥0 and from above by u<1. These bounds can be based on subject matter knowledge or selected with cross-validation. Estimators, incorporating this model knowledge, are expected to be more robust and precise in finite samples.

In simulations, the proposed rare outcomes TMLE performed as well or outperformed the alternative estimators. The PSM estimator and IPTW were biased under misspecification of the propensity score. The standard TMLE algorithm suffered from bias and poor confidence interval coverage when the adjustment set for the conditional mean outcome was large. Both AIPTW and rTMLE were robust to model misspecification and yielded consistent estimates if either the conditional mean outcome or the propensity score were consistently estimated. AIPTW, however, is not a substitution estimator and yielded negative (impossible) risk estimates. In contrast, the proposed TMLE respected the global knowledge in the statistical model. Our simulations further highlighted the potential for data-adaptive estimation to avoid parametric assumptions and to increase power.

We focused on situations where the conditional mean of a rare outcome was bounded from below by ℓ=0 and from above by u<1. The proposed TMLE is equally applicable to situations where the outcome is very common and thereby the conditional mean outcome is bounded from below by 0<ℓ and from above by u=1. Furthermore, the rTMLE algorithm could be applied when we have knowledge of both the lower ℓ and upper bound u and possibly when there are different bounds for different strata of the exposure A and covariates W. We expect that incorporating this knowledge during estimation will increase stability and power especially when the sample size is small and the dimension of the covariates is high. This is an area of future work.

The proposed TMLE is easily generalizable for estimation of other parameters, including the risk ratios, odds ratios and the impacts of longitudinal exposures. The TMLE is also applicable to other sampling designs. Specifically, case-control studies are commonly employed to increase robustness and efficiency of the analysis of rare events. There are several well established methods that correct for selection on the outcome (Anderson 1972; Prentice and Breslow 1978; King and Zeng 2001; Robins 1999; Mansson et al. 2007). For example, van der Laan (2008) presented a general mapping of loss functions, substitution estimators and estimating equations developed for prospective sampling (i. e. cohort sampling) into loss functions, substitution estimators and estimating equations for biased sampling (e. g. case-control sampling). The estimator’s properties, such as double robustness and asymptotic efficiency, are preserved under the mapping. As noted by van der Laan (2008), however, the sample size needed to detect effects on the order of the outcome prevalence will be very large, unless the conditional probability of the outcome is bounded from above by some small constant u. The estimator, proposed in this article, satisfies this condition by construction and thereby has the potential to achieve higher power than its unconstrained counterpart in finite samples. In other words, the mapping provided in van der Laan (2008) will allow us to weight rTMLE appropriately for case-control sampling. We expect the resulting estimator to offer an improvement in terms of stability and power in finite samples.