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 the difference (the g-computation approach). The second models the treatment: estimate each unit’s propensity to be treated and reweight the sample to mimic a randomized experiment (the inverse-probability-weighting [IPW] approach). Each is valid, but each carries a single point of failure. G-computation is consistent only if the outcome model is correctly specified; IPW is consistent only if the propensity model is correct. Choose one, and you have staked the entire analysis on getting that one model right (Robins, Rotnitzky, and Zhao 1994; Bang and Robins 2005).
Doubly robust estimators dissolve this dilemma. They combine an outcome model and a propensity model in a single estimator, remaining consistent if either model is correctly specified, even if the other is wrong (Bang and Robins 2005; Funk et al. 2011). This is the “two chances to be right” property: the analyst no longer has to bet everything on one modeling choice. Among doubly robust estimators, targeted maximum likelihood estimation (TMLE) is distinctive. Rather than plugging model predictions into a formula, TMLE starts from an initial outcome estimate and then updates it, a targeting step that uses information from the propensity model to correct the initial fit 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, so the estimate respects the bounds of the parameter (a probability stays in [0,1]). Second, TMLE is doubly robust and, when both models are estimated well, achieves the efficiency bound, making it a natural home for flexible machine-learning nuisance estimates (Laan and Rose 2011; Schuler and Rose 2017).
In this post, I build TMLE step by step for the average treatment effect (ATE) on 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 stays close to the truth as long as the other model is correct.
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
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
The heart of TMLE is a one-dimensional update to the initial outcome model. 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):
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
An important property of TMLE is that it stays consistent when one of the two models is wrong. 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)results
Scenario TMLE g-computation only Truth
1 Both correct 0.178 0.181 0.159
2 Outcome WRONG, propensity correct 0.175 0.207 0.159
3 Outcome correct, propensity WRONG 0.181 0.181 0.159
4 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 is the practical tool: it 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, Polley, and Hubbard 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 packagedata.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))
The package estimate should be close to the hand-built one, confirming the implementation. 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, Polley, and Hubbard 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.
Estimator Estimate CI_lower CI_upper Truth
1 TMLE (from scratch) 0.178 NA NA 0.159
2 TMLE (package, parametric) 0.177 0.133 0.221 0.159
3 TMLE (package, SuperLearner) 0.177 0.135 0.219 0.159
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.
———. 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. New York: 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.