Two strategies dominate the estimation of treatment effects from observational data. The first models the outcome: fit a regression of the outcome on treatment and covariates, then predict the potential outcomes under treatment and control and average their difference (the g-computation approach). The second models treatment assignment: estimate each unit’s propensity to receive treatment and use inverse-probability weights to construct a weighted population in which treatment assignment is balanced with respect to measured covariates (the inverse-probability-weighting [IPW] approach). Under the usual causal assumptions, each approach can consistently estimate the treatment effect when its corresponding model is correctly specified. G-computation relies on correctly modeling the conditional outcome distribution, whereas IPW relies on correctly modeling the treatment assignment mechanism. Using either approach alone therefore makes the analysis sensitive to misspecification of the model on which the estimator depends (Robins et al. 1994; Bang and Robins 2005).
Doubly robust estimators address this dilemma by combining an outcome model and a propensity model in a single estimator. Under the usual causal identification assumptions, they remain consistent if either the outcome model or the propensity model is correctly specified, even when the other is misspecified (Bang and Robins 2005; Funk et al. 2011). This is the “two chances to be right” property: rather than relying entirely on one modeling strategy, the estimator can remain valid when one of the two nuisance models is wrong.
Among doubly robust estimators, targeted maximum likelihood estimation (TMLE) is distinctive in how it combines these two models. Rather than simply plugging model-based predictions into a formula, TMLE begins with an initial estimate of the outcome regression and then updates it through a targeting step. This update uses information from the estimated treatment mechanism to modify the initial outcome estimate in a direction that is specifically targeted toward the causal parameter of interest (Laan and Rubin 2006; Laan and Rose 2011).
This targeting step gives TMLE two attractive features. First, it is a plug-in estimator: the final estimate is obtained by substituting the updated outcome predictions into the g-computation formula. Because the estimate is constructed from the updated predicted outcomes, it naturally respects the parameter’s bounds—for example, a probability remains in \([0,1]\). Second, TMLE is doubly robust and, when the nuisance functions are estimated consistently under the required regularity conditions, can attain the semiparametric efficiency bound when both models are correctly specified. This combination makes TMLE well suited to incorporating flexible machine-learning methods for estimating the nuisance functions (Laan and Rose 2011; Schuler and Rose 2017).
In this post, I build TMLE step by step for the average treatment effect (ATE) with a binary outcome: the initial outcome model, the propensity model, the targeting update, and the final plug-in estimate. I then demonstrate double robustness directly by deliberately misspecifying one model at a time and showing that TMLE remains close to the truth as long as the other model is correctly specified.
The Estimand and the Ingredients
I target the ATE on a binary outcome, \(\psi = E[Y(1) - Y(0)]\), the difference in the population average outcome if everyone were treated versus if everyone were untreated. TMLE assembles this from two nuisance models (Laan and Rose 2011; Luque-Fernandez et al. 2018).
The first is the outcome model\(\bar{Q}(A, W) = E[Y \mid A, W]\), the expected outcome given treatment \(A\) and covariates \(W\). From an initial fit \(\bar{Q}^0\) we obtain predictions \(\bar{Q}^0(1, W)\) and \(\bar{Q}^0(0, W)\) for each unit under both treatment values. The second is the propensity model\(g(W) = P(A = 1 \mid W)\), the probability of treatment given covariates. These are exactly the two models from the g-computation and IPW approaches, respectively (TMLE uses both).
The bridge between them is the clever covariate, defined for each unit as
\[
H(A, W) = \frac{A}{g(W)} - \frac{1 - A}{1 - g(W)}.
\] This covariate encodes the inverse-probability weights, and it is the vehicle through which propensity information enters the outcome model during targeting. It is the component of the efficient influence function for the ATE that drives the targeting update (Laan and Rose 2011; Schuler and Rose 2017).
The Targeting Step
Working on the logit scale, we fit a logistic regression of the observed outcome \(Y\) on the clever covariate \(H(A, W)\), using the initial prediction \(\bar{Q}^0(A, W)\) as an offset and estimating a single fluctuation parameter \(\epsilon\)(Laan and Rubin 2006; Gruber and Laan 2009):
\[
\text{logit}\,\bar{Q}^1(A, W) = \text{logit}\,\bar{Q}^0(A, W) + \epsilon \, H(A, W).
\] The estimated \(\hat{\epsilon}\) tilts the initial predictions in the direction that removes residual confounding not captured by the outcome model alone. Because \(H\) carries the propensity information, this single-parameter update is what makes the final estimator doubly robust: if the outcome model was already correct, \(\hat{\epsilon}\) is near zero and the update does little; if the outcome model was wrong but the propensity model is right, the update corrects the bias (Laan and Rose 2011). The updated predictions \(\bar{Q}^1(1, W)\) and \(\bar{Q}^1(0, W)\) are then formed by applying the fluctuation with the clever covariate evaluated at \(A = 1\) and \(A = 0\).
The final estimate is the plug-in average of the updated predictions:
The steps are implemented directly. First, data are simulated with a binary outcome, a binary treatment confounded by covariates, and a known treatment effect.
library(dplyr)set.seed(42)n <-2000# CovariatesW1 <-rnorm(n)W2 <-rnorm(n)# Treatment depends on covariates (confounding)g_true <-plogis(-0.4+0.7* W1 -0.5* W2)A <-rbinom(n, 1, g_true)# Outcome depends on treatment and covariatesqbar_true <-plogis(-1+0.8* A +0.6* W1 +0.5* W2 +0.3* W1 * W2)Y <-rbinom(n, 1, qbar_true)dat <-data.frame(Y, A, W1, W2)# True ATE by g-computation on the known outcome modelq1 <-plogis(-1+0.8*1+0.6* W1 +0.5* W2 +0.3* W1 * W2)q0 <-plogis(-1+0.8*0+0.6* W1 +0.5* W2 +0.3* W1 * W2)true_ate <-mean(q1 - q0)round(true_ate, 3)
[1] 0.159
Next, a TMLE function is defined with separate formulas for the outcome model and propensity model. This allows the function to be evaluated under correctly specified and misspecified models to assess double robustness.
tmle_ate <-function(data, q_formula, g_formula) {# Step 1: initial outcome model Qbar^0 q_fit <-glm(q_formula, data = data, family =binomial())# predictions under A = 1 and A = 0 d1 <-transform(data, A =1) d0 <-transform(data, A =0) Q_A <-predict(q_fit, newdata = data, type ="response") # at observed A Q_1 <-predict(q_fit, newdata = d1, type ="response") Q_0 <-predict(q_fit, newdata = d0, type ="response")# Step 2: propensity model g g_fit <-glm(g_formula, data = data, family =binomial()) g_hat <-predict(g_fit, type ="response") g_hat <-pmin(pmax(g_hat, 0.01), 0.99) # bound away from 0/1# Step 3: clever covariate H_A <- data$A / g_hat - (1- data$A) / (1- g_hat) H_1 <-1/ g_hat H_0 <--1/ (1- g_hat)# Step 4: targeting update (fluctuation parameter epsilon)# logistic regression of Y on H with offset = logit(Q_A) eps_fit <-glm( data$Y ~-1+ H_A,offset =qlogis(Q_A),family =binomial() ) eps <-coef(eps_fit)[1]# updated predictions Q1_star <-plogis(qlogis(Q_1) + eps * H_1) Q0_star <-plogis(qlogis(Q_0) + eps * H_0)# Step 5: plug-in ATEmean(Q1_star - Q0_star)}
With both models correctly specified, TMLE should land close to the true ATE.
est_correct <-tmle_ate( dat,q_formula = Y ~ A + W1 + W2 + W1:W2,g_formula = A ~ W1 + W2)data.frame(Estimator ="TMLE (both models correct)",Estimate =round(est_correct, 3),Truth =round(true_ate, 3))
Estimator Estimate Truth
1 TMLE (both models correct) 0.178 0.159
Demonstrating Double Robustness
This is tested directly by deliberately misspecifying each model in turn, dropping the interaction and a covariate to break the model, while keeping the other correct. For contrast, the naive g-computation estimate (outcome model only) is also computed under the same misspecification, which has no second chance to be right.
# Naive g-computation (outcome model only), for comparisongcomp_ate <-function(data, q_formula) { fit <-glm(q_formula, data = data, family =binomial()) Q1 <-predict(fit, newdata =transform(data, A =1), type ="response") Q0 <-predict(fit, newdata =transform(data, A =0), type ="response")mean(Q1 - Q0)}# Correct and misspecified formulasq_correct <- Y ~ A + W1 + W2 + W1:W2q_wrong <- Y ~ A # omits all covariatesg_correct <- A ~ W1 + W2g_wrong <- A ~1# intercept only (ignores confounding)results <-data.frame(Scenario =c("Both correct","Outcome WRONG, propensity correct","Outcome correct, propensity WRONG","Both wrong" ),TMLE =c(tmle_ate(dat, q_correct, g_correct),tmle_ate(dat, q_wrong, g_correct),tmle_ate(dat, q_correct, g_wrong),tmle_ate(dat, q_wrong, g_wrong) ),`g-computation only`=c(gcomp_ate(dat, q_correct),gcomp_ate(dat, q_wrong),gcomp_ate(dat, q_correct),gcomp_ate(dat, q_wrong) ),check.names =FALSE)results$TMLE <-round(results$TMLE, 3)results$`g-computation only`<-round(results$`g-computation only`, 3)results$Truth <-round(true_ate, 3)knitr::kable(results, digits =3)
Scenario
TMLE
g-computation only
Truth
Both correct
0.178
0.181
0.159
Outcome WRONG, propensity correct
0.175
0.207
0.159
Outcome correct, propensity WRONG
0.181
0.181
0.159
Both wrong
0.207
0.207
0.159
The pattern in this table is the whole point. When the outcome model is wrong but the propensity model is correct, TMLE remains close to the truth while g-computation alone is biased—the targeting step, powered by the correct propensity model, rescues the estimate. When the propensity model is wrong but the outcome model is correct, TMLE is again close, because a correct outcome model needs no correction. Only when both models are wrong does TMLE fail, which is the price no estimator can avoid. The single-model g-computation, by contrast, is biased whenever its one model is misspecified, with no second line of defense.
Using the tmle Package
Implementing TMLE by hand is the best way to understand it, but for real analyses the tmle package handles the targeting step, bounds, influence-curve-based standard errors, and lets you plug in flexible machine-learning estimators for the nuisance models through SuperLearner(Gruber and Laan 2012; Laan et al. 2007). Running it on the same data is also a useful check on our from-scratch estimate: the two should agree closely when given the same model specifications.
The tmle() function takes the outcome Y, treatment A, and covariate matrix W directly. To match the hand-rolled version, just supply simple parametric specifications through the Qform and gform arguments.
library(tmle)W <- dat[, c("W1", "W2")]fit_pkg <-tmle(Y = dat$Y,A = dat$A,W = W,Qform = Y ~ A + W1 + W2 + W1:W2, # outcome modelgform = A ~ W1 + W2, # propensity modelfamily ="binomial")# ATE estimate and inference from the packageknitr::kable(data.frame(Estimator =c("TMLE (from scratch)", "TMLE (tmle package)"),Estimate =round(c(est_correct, fit_pkg$estimates$ATE$psi), 3),Truth =round(true_ate, 3)),digits =3)
Estimator
Estimate
Truth
TMLE (from scratch)
0.178
0.159
TMLE (tmle package)
0.177
0.159
Beyond the point estimate, tmle also returns a variance estimate and confidence interval derived from the efficient influence function, which the minimal version omitted.
# Influence-function-based inference from the packageci <- fit_pkg$estimates$ATE$CIdata.frame(Estimate =round(fit_pkg$estimates$ATE$psi, 3),CI_lower =round(ci[1], 3),CI_upper =round(ci[2], 3),p_value =signif(fit_pkg$estimates$ATE$pvalue, 3))
The real strength of the package emerges when the nuisance models are estimated with SuperLearner, an ensemble that combines many candidate learners (splines, random forests, gradient boosting, and more) and lets the data choose their weighting (Laan et al. 2007). This pairs naturally with TMLE: the double-robustness and efficiency theory is what justifies plugging flexible, potentially black-box estimators into the nuisance models while still obtaining valid inference for the target parameter.
Instead of the parametric Qform/gform specification used above, each nuisance model is now supplied with a library of learners. SuperLearner fits the candidate learners, evaluates their performance using cross-validation, and combines them into a weighted ensemble. TMLE then targets the resulting nuisance estimates and produces the corresponding inference.
The SuperLearner-based estimate can be compared directly with the parametric package estimate and the from-scratch implementation. All three target the same ATE; the SuperLearner variant differs only in its more flexible estimation of the nuisance functions.
The three estimates agree closely, since the data-generating process here is well approximated by the parametric models, so the flexible learners have little extra structure to find. The value of SuperLearner is insurance for the realistic case where the nuisance relationships are not known to be simple: the ensemble adapts to nonlinearities and interactions the analyst has not anticipated, and the TMLE targeting step preserves valid inference for the ATE regardless of which learners end up carrying the weight. That combination, flexible nuisance estimation with a targeting step that restores \(\sqrt{n}\) inference for the parameter of interest, is precisely why TMLE with SuperLearner has become a standard workhorse for modern causal effect estimation (Laan and Rose 2011; Schuler and Rose 2017).
In practice, the tmle and tmle3 packages implement these ideas with cross-fitting and machine-learning nuisance estimators, making them a natural next step beyond the hand-rolled version presented here.
Thank you for reading, and I welcome your feedback.
References
Bang, Heejung, and James M. Robins. 2005. “Doubly Robust Estimation in Missing Data and Causal Inference Models.”Biometrics 61 (4): 962–73.
Funk, Michele Jonsson, Daniel Westreich, Chris Wiesen, Til Stürmer, M. Alan Brookhart, and Marie Davidian. 2011. “Doubly Robust Estimation of Causal Effects.”American Journal of Epidemiology 173 (7): 761–67.
Gruber, Susan, and Mark J. van der Laan. 2009. “Targeted Maximum Likelihood Estimation: A Gentle Introduction.”U.C. Berkeley Division of Biostatistics Working Paper Series.
Gruber, Susan, and Mark J. van der Laan. 2012. “Tmle: An r Package for Targeted Maximum Likelihood Estimation.”Journal of Statistical Software 51 (13): 1–35.
Laan, Mark J. van der, Eric C. Polley, and Alan E. Hubbard. 2007. “Super Learner.”Statistical Applications in Genetics and Molecular Biology 6 (1).
Laan, Mark J. van der, and Sherri Rose. 2011. Targeted Learning: Causal Inference for Observational and Experimental Data. Springer.
Laan, Mark J. van der, and Daniel Rubin. 2006. “Targeted Maximum Likelihood Learning.”The International Journal of Biostatistics 2 (1).
Luque-Fernandez, Miguel Angel, Michael Schomaker, Bernard Rachet, and Mireille E. Schnitzer. 2018. “A Tutorial on Targeted Maximum Likelihood Estimation for Causal Inference.”Statistics in Medicine 37 (16): 2530–46.
Robins, James M., Andrea Rotnitzky, and Lue Ping Zhao. 1994. “Estimation of Regression Coefficients When Some Regressors Are Not Always Observed.”Journal of the American Statistical Association 89 (427): 846–66.
Schuler, Megan S., and Sherri Rose. 2017. “Targeted Maximum Likelihood Estimation for Causal Inference in Observational Studies.”American Journal of Epidemiology 185 (1): 65–73.