There are many fine how-to articles describing how to implement TMLE but they avoid the details of the underlying theory. I'm currently working my way through Targeted Learning: Causal Inference for Observational and Experimental Data by Mark J. van der Laan and Sherri Rose. The math isn't terribly complicated but the notation and terminology is a bit confusing.

I understand TMLE's aim of finding an unbiased estimate of the Average Treatment Effect by using machine learning, and am familiar with the theory behind causal inference, the Super Learner algorithm, and doubly robust models, but I hit a brick wall when it comes to calculating the efficient influence curve, the "clever covariate", and guaranteeing the final ATE estimate's unbiasedness with the Central Limit Theorem.

My understanding is that TMLE uses the delta method (1st order Taylor series) to approximate the ATE and then converges to an estimate of the ATE via gradient descent(?). Am I too far off?

  • 1
    $\begingroup$ See also stats.stackexchange.com/questions/407444/…, stats.stackexchange.com/questions/134572/… $\endgroup$ Commented Dec 29, 2019 at 6:01
  • $\begingroup$ @kjetilbhalvorsen I did read those questions but unfortunately they weren't helpful. Possibly targeted learning is still such a young field that aside from Mark J. van der Laan and Sherri Rose there aren't a lot of statisticians who do completely understand TMLE theory? $\endgroup$
    – RobertF
    Commented Dec 29, 2019 at 14:55
  • $\begingroup$ I should also add that I haven't been able to find a copy of Frank R. Hampel's 1974 paper on influence curves, "The Influence Curve and its Role in Robust Estimation" that isn't stuck behind a paywall. $\endgroup$
    – RobertF
    Commented Dec 29, 2019 at 15:06
  • 2
    $\begingroup$ Found an article that does a good job explaining the theory behind using influence functions to estimate ATEs: arxiv.org/pdf/1810.03260.pdf $\endgroup$
    – RobertF
    Commented Jan 1, 2020 at 22:58
  • 1
    $\begingroup$ David Benkeser and Antoine Chambaz wrote a useful explanation of TMLE theory here: achambaz.github.io/tlride $\endgroup$
    – RobertF
    Commented Feb 5, 2021 at 20:52

1 Answer 1


Here is a very tentative (incomplete) answer. I will add more (e.g. the TMLE for ATE) later. It is fairly technical.

One important thing to keep in mind is the separation between the data structure, statistical model & the parameter of interest, and the estimation of the target estimand (e.g. the actual TMLE procedure).

Definition of statistical model

Consider the point-treatment data structure $O=(W,A,Y)$ where $W$ are baseline covariates, $A$ is a binary treatment and $Y$ is an outcome variable. We assume $O$ is drawn from the probability distribution $P_0$. Let $\mathcal{O}$ be the set of possible values of $O$ (e.g. $\mathbb{R}^{d+2}$ if $W$ is d-dimensional).

The statistical model $\mathcal{M}$ is a collection of probability distributions that we assume contains $P_0$. In the context of TMLE, we usually take $\mathcal{M}$ to be nonparametric so that it contains all possible probability distributions. For our data structure, we are only interested in certain nuisance components of $P_0$ rather than the entire probability distribution. In particular, we are interested in the distribution of $W$, $P_W$, the propensity score $g(W):=P_0(A=1|W)$, and the conditional mean of the outcome $Q(A,W):=E_0[Y|A,W]$.

Definition of parameter

Once, we have defined the data structure and statistical model, we can begin talking about our statistical quantity of interest, otherwise known as the target parameter (in this case ATE). Rigorously, the target parameter is a mapping from the statistical model to the real line (or more generally $\mathbb{R}^d$). That is a target parameter $\Psi: \mathcal{M} \mapsto \mathbb{R}$ assigns to each distribution $P \in \mathcal{M}$ the value $\Psi(P)$. The statistical estimand $\psi_0 := \Psi(P_0)$ is defined as the value of the parameter at the true data-generating distribution $P_0$. The ATE parameter is given by $\Psi_{ATE}(P) = E_{P,W}E_P[Y|A=1,W] - E_{P,W}E_P[Y|A=0,W]$.

Theory of efficient influence function

A useful reference is the following: https://arxiv.org/abs/1903.01706 which has a similar but more thorough treatment as below. Another great reference is the following: https://arxiv.org/pdf/2107.00681.pdf

The efficient influence function (EIF) is determined by the statistical model $\mathcal{M}$ and the target parameter $\Psi$. It happens to be a very useful object for efficient estimation and inference, but computing it has nothing to do with estimation. Actually, for each $P$ in the statistical model $\mathcal{M}$, there is an efficient influence function. For this reason, we usually write the EIF as $D_P:O\mapsto D_P(O)$, where it is stressed that the EIF is indexed by $\mathcal{M}$. The efficient influence function $D_P$ can be viewed as the gradient of the parameter $\Psi$ at $P$ in the statistical model $\mathcal{M}$ (analogous to the notion of a gradient of a multivariate function from multivariate calculus). To rigorously define the EIF requires quite a bit of mathematics. In particular, we need to generalize the notions of (directional) derivatives and gradients to functionals (like $\Psi$) defined on infinite dimensional spaces (like $\mathcal{M}$). If you are familiar with aspects of differential geometry or manifold theory, specifically the concept of tangent spaces and derivative operators defined on them, and functional analysis, specifically Hilbert space theory and projections, then that will be very helpful. One key property of the EIF $D_P$ is that it satisfies for $P,P_1 \in \mathcal{M}$ that $$\Psi(P_1) = \Psi(P) + \int D_P(o)d(P_1-P)(o) + R(P_1,P),$$ where $R({P_1,P})$ is a remainder that goes to zero when $P_1$ tends towards $P$ in an appropriate sense. Thus, $D_P$ encodes local (first order derivative) information of $\Psi$ at $P$ and allows for a first order taylor expansion.

I will take a shot at semi-rigorously defining the general notion of an efficient influence function. I will assume that each distribution in $\mathcal{M}$ has a density (with respect to the Lebesgue measure). We can thus view $\mathcal{M}$ as a collection of probability densities and $\Psi$ as a functional of densities. Let $P \in \mathcal{M}$ arbitrary and let $p$ be its density.

First, we need to define the notion of a directional derivative of $\Psi$. This first requires defining paths through a given density $p$. To this end, for each $P \in \mathcal{M}$, define the Tangent Space: $$T_P\mathcal{M} := (O \mapsto h(O): E_{P}h(O)^2 < \infty,\, E_{P}h(O) = 0),$$ which is just the space of mean-zero and finite variance functions of $O$ (i.e. $L^2_0(P_0) \subset L^2(P_0)$). It is well-known that $T_P\mathcal{M}$ is a closed linear space that can be equipped with the inner product $$\langle h, g \rangle_P := E_P[h(O)g(O)]= \int h(o) g(o) p(o)do,$$ in which case it becomes a Hilbert space. For some $h \in T_P\mathcal{M}$, consider the path through $p$, $$p_{\varepsilon,h} = (1+ \varepsilon h)p,$$ where $\varepsilon \in \mathbb{R}$. Note that $$\int p_{\varepsilon,h}(o)do = \int p(o)d0 + \varepsilon E_P[h(O)] = 1 + 0 =1, $$ and for small enough $\varepsilon$, we have $p_{\varepsilon,h} \geq 0$. Therefore $p_{\varepsilon,h}$ is indeed a density corresponding with some distribution $P_{\varepsilon,h} \in \mathcal{M}$ for all $\varepsilon$ small enough. It is crucial that $P_{\varepsilon,h} \in \mathcal{M}$, which is not necessarily true for statistical models that are not fully nonparametric). We can define the derivative/direction/score of the path $P_{\varepsilon,h}$ as $$\frac{d}{d\varepsilon} \log p_{\varepsilon,h} \Big |_{\varepsilon = 0} = \frac{1}{p} \frac{d}{d\varepsilon} (1+ \varepsilon h)p |_{\varepsilon = 0} = \frac{ph}{p} = h.$$ Thus, the score of the path $P_{\varepsilon,h}$ is exactly equal to tangent space element $h$. (We use $\frac{d}{d\varepsilon,h} \log p_{\varepsilon,h} $ instead of $\frac{d}{d\varepsilon} p_{\varepsilon,h} $ due to nice properties of the former). $T_P\mathcal{M}$ is called the tangent space at $P$ because it encodes all possible directions one can move at $P$ while still staying in the statistical model. Given $h \in T_P\mathcal{M}$ and its corresponding path $P_{\varepsilon,h}$, we define the directional derivative of $\Psi$ in the direction $h$ as $$\frac{d}{d\varepsilon} \Psi(P_{\varepsilon,h}) := \lim_{\varepsilon \rightarrow 0} \frac{\Psi(P_{\varepsilon,h}) - \Psi(P_{0,h})}{\varepsilon}.$$ Note that $\varepsilon \mapsto \Psi(P_{\varepsilon,h})$ is just a univarate real-valued function so that this derivative can be computed using standard calculus. We now define the total derivative $d\Psi_P:T_P\mathcal{M} \mapsto \mathbb{R}$ of $\Psi$ pointwise as $$d\Psi_P(h) = \frac{d}{d\varepsilon} \Psi(P_{\varepsilon,h}).$$ Under some assumptions, $d\Psi_P$ is a linear and bounded/continuous map defined on the Hilbert space $T_P\mathcal{M}$. It then follows by Reisz representation theorem that there exists a unique gradient in $D_P \in T_P \mathcal{M}$ such that $$d\Psi_P(h) = \langle h, D_P \rangle_P = E_P h(O)D_P(O).$$ This gradient is exactly the efficient influence function of $\Psi$ at $P$. Under some conditions, one has that $$\Psi(P_{\varepsilon,h}) = \Psi(P) + \langle h, D_P \rangle_P + o(\varepsilon),$$ so that $D_P$ encodes a first order taylor expansion of $\Psi$ at $P$. Since $\mathcal{M}$ is convex at the density level, for any $P_1,P \in \mathcal{M}$, we can choose the path $p_{\varepsilon} = p + \varepsilon(p_1-p)$ which has score $h = \frac{p_1-p}{p}$ to find that $$\Psi(P_1) = \Psi(P) + \int D_P(o)d(P_1-P)(o) + R(P_1,P),$$ where $R(P_1,P)$ should go to zero as $P_1$ gets "close to" $P$.

Key properties of the efficient influence function

The efficient influence function $D_P$ has a number of key statistical and non-statistical properties. Firstly, as shown in the previous section, we have for smooth paths $P_{\varepsilon}$ with $P_{\varepsilon=0} = P$ and score $\frac{d}{d\varepsilon} \log p_{\varepsilon} \Big |_{\varepsilon=0} = h$ that $$d\Psi_P(h) = \frac{d}{d\varepsilon} \Psi(P_{\varepsilon}) \Big |_{\varepsilon=0} = E_P[D_P(O)h(O)].$$ Moreover, for any $P_1,P \in \mathcal{M}$ $$\Psi(P_1) = \Psi(P) + \int D_P(o)d(P_1-P)(o) + R(P_1,P),$$ where $R(P_1,P)$ should go to zero as $P_1$ gets "close to" $P$. Also, the EIF has mean-zero so the above simplifies to $\Psi(P_1) = \Psi(P) + \int D_P(o)dP_1(o) + R(P_1,P).$

Another property of $D_P$ is as follows. Consider the path $P_{\varepsilon}$ induced by the density $p_{\varepsilon} = (1 + \varepsilon \frac{D_P}{\sqrt{E_P[D_P(O)^2]}})p$. We say that this path moves in the direction of $D_P$. Since $D_P$ is mean zero, this is a locally valid path through probability distributions. I claim that this path is in the direction of steepest change of $\Psi$ at $P$ (just as the gradient is in multivariate calculus). To see why this is true, note for any other path with score $h$ (e.g. $\varepsilon \mapsto p(1+\varepsilon h)$ ) normalized to have $E_P[h(O)^2]=1$, we have the derivative of $\Psi$ along this path is $$d\Psi_P(h) = E_P[D_P(O)h(O)].$$ By the Cauchy-Schwartz inequality, we have $$d\Psi_P(h) = E_P[D_P(O)h(O)] \leq \sqrt{E_P[D_P(O)^2]} \sqrt{E_P[h(O)^2]} = \sqrt{E_P[D_P(O)^2]}.$$ Thus, the directional derivative can never exceed $\sqrt{E_P[D_P(O)^2]}$. However, we can compute $$d\Psi_P({D_P}/{\sqrt{E_P[D_P(O)^2]}}) = \frac{1}{\sqrt{E_P[D_P(O)^2]}}E_P[D_P(O)D_P(O)] = \sqrt{E_P[D_P(O)^2]} .$$ Thus, the path $P_{\varepsilon}$ in direction ${D_P}/{\sqrt{E_P[D_P(O)^2]}}$ has the maximal directional derivative along all paths through $P$. More generally, we can show that the normalized local maximal change in $\Psi$ is given by $$\sup_{h \in T_P\mathcal{M}}\frac{d\Psi_P(h)}{\sqrt{E_P[h(O)^2]}} = \sqrt{E_P[D_P(O)^2]}$$ which is achieved when $h = D_P$. (The key idea behind TMLE is to perform parametric maximum likelihood estimation along a path in the direction of $D_P$ as this will give us the largest change in $\Psi$ with the least amount of fitting.)

The key statistical property is that $D_P$ encodes a generalization of the Cramer-Rao lower to bound to the infinite-dimensional statistical model $\mathcal{M}$. In particular, the best possible standardized variance for asymptotically unbiased estimators of $\Psi(P_0)$ (with some regularity conditions) is given by $E_{P_0} D_{P_0}(O)^2$ (The best possible distribution is Gaussian). (This is the backbone of semiparametric efficiency theory. See for example Asymptotic Statistics, van der Vaart, 2000 for an introduction to the relevant theory). Thus, if one is able to construct an estimator $\Psi_n$ such that $\sqrt{n}(\Psi_n - \Psi(P_0))$ has an asymptotic distribution being a mean-zero Gaussian random variable with variance $E_{P_0} D_{P_0}(O)^2$ then one knows that the estimator is efficient. Usually, we seek an estimator that satisfies $$\Psi_n - \Psi(P_0) = \frac{1}{n}\sum_{i=1}^n D_{P_0}(O_i) + o_P(n^{-1/2}).$$ Such an estimator is (when centered) asymptotically equivalent to the sample mean of $D_{P_0}(O)$ and therefore trivially has $\sqrt{n}(\Psi_n - \Psi(P_0))$ goes to a mean-zero Gaussian random variable with variance $E_{P_0} D_{P_0}(O)^2$ as desired. Somewhat confusingly $D_{P_0}$ is then called the "influence function" of the estimator $\Psi_n$. This is also why some sources might be confusing (e.g. robust statistics concerns itself with the influence functions of estimators). This is also why the efficient influence function is sometimes called the "canonical gradient".

TMLE as well as double-robust estimators more generally (e.g. one-step efficient estimator) attempts to construct an estimator that is asymptotically equivalent to the sample mean of $D_{P_0}$.

One Step efficient estimators

The first nonparametric efficient estimator is the so-called one-step estimator. It is motivated by the first order taylor expansion of the previous section. Let $P_{n,0} \in \mathcal{M}$ be an initial estimator of $P_0$. Let $P_n$ denote the empirical measure. Almost always $\Psi(P_{n,0})$ is overly biased and will not be $\sqrt{n}$-consistent and therefore surely not efficient. The first order taylor expansions motivates performing a first-order bias correction (essentially one newton gradient descent step update in the parameter space). We can write (doing the taylor expansion at $P_{n,0}$ in the direction of the score $d(P_0-P_{n,0})/dP_{n,0}$), $$\Psi(P_{n,0}) - \Psi(P_0) = -d\Psi_{P_{n,0}}(d(P_0-P_{n,0})/dP_{n,0}) + R_{n,2}(P_{n,0},P_0) $$ $$= -\int D_{P_{n,0}}(o) dP_0(o) + R_{n,2}(P_{n,0},P_0)$$ where $R_{n,2}(P_{n,0},P_0)$ is a remainder term. Under some conditions, one can show that $R_{n,2}(P_{n,0},P_0) = o_P(n^{-1/2})$ so that it is asymptotically negligible. If we knew the true distribution $P_0$, a natural estimator would be $\Psi_{n,oracle} := \Psi(P_{n,0}) + \int D_{P_{n,0}}(o) dP_0(o) = \Psi(P_{n,0}) + E_{P_0}[D_{P_{n,0}}(O)]$ which would under conditions satisfy $\Psi_{n,oracle} - \Psi(P_0) = o_P(n^{-1/2})$ and be faster than $\sqrt{n}$-consistent. This is just a first-order gradient-descent/newton update for $\Psi$. Unfortunately, this estimator is impossible as we don't know $P_0$. Instead, we will look at an empirical version of this estimator given by $\Psi(P_{n,0}) + E_{P_n}[D_{P_{n,0}}(O)]$, which is tractable.

Now, I am going to add $\int D_{P_{n,0}}(o) dP_n(o)$ to both sides to get $$\Psi(P_{n,0}) - \Psi(P_0) + \int D_{P_{n,0}}(o) dP_n(o) = \int D_{P_{n,0}}(o) d(P_n-P_0)(o) + R_{n,2}(P_{n,0},P_0).$$ This can be further written as $$\Psi(P_{n,0}) - \Psi(P_0) + \int D_{P_{n,0}}(o) dP_n(o) = \int D_{P_{0}}(o) d(P_n-P_0)(o) + \int (D_{P_{n,0}}(o)-D_{P_{0}}(o)) d(P_n-P_0)(o) + R_{n,2}(P_{n,0},P_0).$$ When $P_{n,0}$ is consistent for $P_0$, then $(D_{P_{n,0}}(o)-D_{P_{0}}(o))$ converges to zero in probability. In particular, $\int (D_{P_{n,0}}(o)-D_{P_{0}}(o)) d(P_n-P_0)(o) $ is the sample mean of something converging to zero and is therefore under some conditions $o_P(n^{-1/2})$. Putting it all together, we get $$\Psi(P_{n,0}) - \Psi(P_0) + \int D_{P_{n,0}}(o) dP_n(o) = \int D_{P_{0}}(o) d(P_n-P_0)(o) + o_P(n^{-1/2}) = \frac{1}{n}\sum_{i=1}^n D_{P_0}(O_i) + o_P(n^{-1/2})$$ as desired (last equality uses that $D_{P_0}$ is necessarily mean zero). Thus, $$\Psi(P_{n,0}) + \int D_{P_{n,0}}(o) dP_n(o) - \Psi(P_0) $$ $$= \Psi(P_{n,0}) + \frac{1}{n}\sum_{i=1}^n D_{P_{n,0}}(O_i) - \Psi(P_0) $$ is asymptotically equivalent to the sample mean of $D_{P_0}$. Thus, $\Psi(P_{n,0}) + \int D_{P_{n,0}}(o) dP_n(o)$ is an efficient estimator, under some conditions.

The one-step estimator is given by $$\Psi_{n, onestep} = \Psi(P_{n,0}) + \frac{1}{n}\sum_{i=1}^n D_{P_{n,0}}(O_i).$$ It requires an initial estimator $\Psi(P_{n,0})$ of $\Psi(P_0)$ and an initial estimator $ D_{P_{n,0}}$ of the efficient influence function $ D_{P_{0}}$. It then bias corrects the initial estimator by performing a single gradient-descent-type update. It is important to note that the gradient descent occurs in the parameter space (the space of values of $\Psi(P_0)$, e.g. $\mathbb{R}$) and not in the model $\mathcal{M}$. This has some undesirable properties. First, $\Psi_{n, onestep}$ is usually not a substitution estimator. That is, there usually does not exist a probability distribution $P_{n}^* \in \mathcal{M}$ that satisfies $\Psi(P_n^*) = \Psi_{n, onestep}$. Because of this, $\Psi_{n, onestep}$ might give erratic and impossible values since it does not respect the constraints of the statistical model. For instance if $P \mapsto \Psi(P)$ takes values in $[0,1]$ then the substitution estimator $\Psi(P_{n,0})$ will also, but $ \Psi_{n, onestep}$ does not necessarily satisfy these constraints. This can lead to poor finite sample performance of the one-step estimator.

This brings us to TMLE which allows one to construct an efficient substitution estimator $\Psi(P_n^*)$ of $\Psi(P_0)$. TMLE accomplishes this by performing gradient descent (specifically maximum likelihood estimation) in the statistical model $\mathcal{M}$.


Targeted maximum likelihood estimation (TMLE) constructs an estimator that is both efficient and a substitution estimator. One way of thinking of TMLE is as a very special kind of one-step estimator. Suppose we had an estimator $P_n^*$ of $P_0$. The one-step estimator is given by $$\Psi(P_n^*) + \frac{1}{n}\sum_{i=1}^n D_{P_n^*}(O_i).$$ Now, suppose that $P_n^*$ satisfies $$\frac{1}{n}\sum_{i=1}^n D_{P_n^*}(O_i) = 0.$$ Then, we would have $$\Psi(P_n^*) = \Psi(P_n^*) +0 =\Psi(P_n^*) + \frac{1}{n}\sum_{i=1}^n D_{P_n^*}(O_i),$$ so that the substitution estimator $ \Psi(P_n^*) $ is equal to the one-step estimator and therefore is also efficient. TMLE provides a general template for updating/calibrating/targeting an arbitrary initial estimator $P_{n,0} \in \mathcal{M}$ into an updated estimator $P_n^*$ that satisfies $\frac{1}{n}\sum_{i=1}^n D_{P_n^*}(O_i)=0$. The key idea is the following. Construct a submodel/path $P_{n,\varepsilon}$ with $P_{n, \varepsilon = 0} = P_{n,0}$ that satisfies $\frac{d}{d\varepsilon} \log p_{n,\varepsilon} = D_{P_{n,\varepsilon}}$. Recalling the previous sections, $P_{n,\varepsilon}$ is a path that for each $\varepsilon$ locally moves in the direction of greatest change of $\Psi$. TMLE performs maximum likelihood estimation along this path (where the initial estimator ensures that this parametric model is close enough to $P_0$ to be useful). This can be viewed as functional gradient descent under the constraint that the empirical likelihood always increases. We update an initial estimator by performing gradient descent in the direction of maximum change of $\Psi$ (whether we go forwards to backwards depends on which direction increases the likelihood). We keep performing gradient descent until the empirical likelihood no longer increases. At this point, we have maximized an empirical likelihood and solved a certain score/derivative equation. This equation we solve is exactly $\frac{1}{n}\sum_{i=1}^n D_{P_n^*}(O_i)=0$.

Consider the empirical log likelihood risk function: $$R_{n,Loglik}(P) = \frac{1}{n}\sum_{i=1}^n \log p(O_i).$$ In view of the above, we have $$\frac{d}{d\varepsilon} R_{n,Loglik}(P_{n,\varepsilon}) = \frac{1}{n}\sum_{i=1}^n \frac{d}{d\varepsilon} \log p_{n,\varepsilon}(O_i) \approx \frac{1}{n}\sum_{i=1}^n D_{P_{n,\varepsilon}}(O_i). $$ Let $\varepsilon_n^* = \text{argmin}_{\varepsilon}R_{n,Loglik}(P_{n,\varepsilon}) $ then we must have (since it is a minimizer) $$\frac{d}{d\varepsilon} R_{n,Loglik}(P_{n,\varepsilon}) \Big|_{\varepsilon = \varepsilon_n^*} = \frac{1}{n}\sum_{i=1}^n D_{P_{n,\varepsilon_n^*}}(O_i) = 0.$$ Thus, $P_n^* := P_{n,\varepsilon_n^*}$ has the desired property that $\frac{1}{n}\sum_{i=1}^n D_{P_{n,\varepsilon_n^*}}(O_i) = 0$. We call $P_n^*$ the TMLE.

Such a path $P_{n,\varepsilon_n}$ can be approximately defined as follows. Fix a step size $\delta > 0$. Initialize $k=0$. Define $p_{n,\varepsilon=0} = p_{n,0}$. For $\varepsilon \in [0,\delta] = [k\delta,(k+1)\delta]$, define $p_{n,\varepsilon} = (1+\varepsilon D_{P_{n,0}})p_{n,0}$. Define $p_{n,k} = p_{n,\varepsilon k}$. Now, arguing recursively for $k=k+1$, for $\varepsilon \in [k\delta,(k+1)\delta]$, define $p_{n,\varepsilon} = (1+\varepsilon D_{P_{n,k}})p_{n,k}$. This path satisfies $$\frac{d}{d\varepsilon} \log p_{n,\varepsilon} \Big |_{\varepsilon = k \delta} = D_{P_{n,k\delta}},$$ so the violation of $\frac{d}{d\varepsilon} \log p_{n,\varepsilon} = D_{P_{n,\varepsilon}}$ is usually negligible as long as the step size $\delta$ is chosen small enough.

Note in practice, we don't estimate the entire probability distribution $P_0$. We estimate key nuisance features like conditional means. A version of TMLE (targeted minimum loss estimation) instead performs the above in a lower dimensional feature space (e.g. the space of conditional means). In this case, we perform functional gradient descent in a lower dimensional space (but still infinite dimensional) under the constraint that the risk of the nuisance estimator does not increase (where the risk is given by some custom risk function). This is also where the "clever covariates" come into play (as a computational trick).

A simple TMLE for the CDF of a univariate CDF at a point.

Suppose we observe $n$ iid observations of a real-valued random variable $X \sim P_0$ with CDF $x \mapsto F_0(x)$. The statistical model $\mathcal{M}$ is given by all probability distributions on $\mathbb{R}$. The parameter of interest is the CDF at $t \in \mathbb{R}$, $$\Psi(P) = P(X \leq t)$$ and the estimand of interest is $\Psi(P_0) = F_0(t) = P_0(X \leq t)$.

We will first compute the efficient influence function of $\Psi$. Let $p$ be some density associated with a $P \in \mathcal{M}$. Consider the path $p_{\varepsilon}: \varepsilon \mapsto (1+ \varepsilon h)p$ where $h$ satisfies $\int h(x) p(x)dx = 0$ and has finite variance. Let $P_{\varepsilon}$ be the path in $\mathcal{M}$ corresponding with density path $p_{\varepsilon}$. Note the directional/pathwise derivative is given by $$d\Psi(h) = \frac{d}{d\varepsilon} \Psi(P_{\varepsilon}) \Big |_{\varepsilon = 0} = \frac{d}{d\varepsilon} E_{P_{\varepsilon}}[1(X \leq t)] \Big |_{\varepsilon = 0}$$ $$ = \frac{d}{d\varepsilon} \int 1(x \leq t) (1+ \varepsilon h(x))p(x)dx \Big |_{\varepsilon = 0}$$ $$ = \int 1(x \leq t) h(x)p(x)dx = E_P[1(X\leq t) h(X)] .$$ Note that $ E_P[1(X\leq t) h(X)]$ corresponds with the $L^2(P)$-inner product of $1(X \leq t)$ and $h(X)$, where $h$ is actually an arbitrary element in the tangent space $T_P\mathcal{M}$. This might incorrectly suggest that $1(X\leq t)$ is the efficient influence function. This is incorrect because the efficient influence function must be contained in $T_P\mathcal{M}$ itself. $X \mapsto 1(X\leq t)$ has finite variance but it is not mean-zero, so it is not in the tangent space. In other words, $x \mapsto p(x)(1+\varepsilon 1(x\leq t))$ is not a valid probability density for any $\varepsilon$. Fortunately, there is an easy fix. Note the identity $$ E_P[1(X\leq t) h(X)] = E_P[[1(X\leq t) - P(X\leq t) ]h(X)],$$ where $1(X\leq t) - P(X\leq t) $ being mean-zero is contained in the tangent space. Thus, we have shown $$d\Psi(h) = \langle D_P, h \rangle_{L^2_0(P)} $$ for $D_P(X) := 1(X\leq t) - P(X\leq t)$ and all $h \in T_P\mathcal{M}$. It follows (from uniqueness) that $D_P$ is the efficient influence function of $\Psi$ at $P$.

We can now construct an efficient one-step estimator. Let $p_n$ be an arbitrary initial estimator of the density $p_0$ corresponding with $P_0$. We have (slightly abusing notation) $$\Psi(p_n) = \int 1(x \leq t)p_n(x)dx,$$ $$D_{p_n}(x) = 1(x \leq t) - \int 1(x \leq t)p_n(x)dx.$$ The one-step estimator is given by $$\Psi_{n,onestep} = \Psi(p_n) + \frac{1}{n}\sum_{i=1}^n D_{p_n}(X_i) $$ $$= \int 1(x \leq t)p_n(x)dx + \frac{1}{n}\sum_{i=1}^n [1(X_i \leq t) - \int 1(x \leq t)p_n(x)dx] $$ $$ = + \frac{1}{n}\sum_{i=1}^n 1(X_i \leq t) + 0 = \frac{1}{n}\sum_{i=1}^n 1(X_i \leq t),$$ which is nothing else than the well-known empirical CDF estimator. Thus, the one-step estimator and the empirical CDF estimator are the same regardless of the initial estimator of the density (In this case, we remarkably don't even need $p_n$ to be consistent).

We can also construct a TMLE for $\Psi(P_0)$. Define the empirical log-likelihood risk function: $$R_n(p) = -\frac{1}{n} \sum_{i=1}^n \log p(X_i).$$ Given an initial estimator $p_n$, we define the path $$p_{n,\varepsilon} = (1 + \varepsilon D_{p_n})p_n.$$

The TMLE $p_n^*$ is given by performing maximum likelihood along this path. In this case, $\varepsilon_n^* := \text{argmin}_{\varepsilon} R_n(p_{n,\varepsilon})$ has an analytical solution. The solution solves the following equation $$\frac{1}{n} \sum_{i=1}^n \frac{1}{(1 + \varepsilon_n^* D_{p_n}(X_i))p_n(X_i)} D_{p_n}(X_i)p_n(X_i) = 0$$ which is equivalent to $$\frac{1}{n} \sum_{i=1}^n \frac{ D_{p_n}(X_i)}{(1 + \varepsilon_n^* D_{p_n}(X_i))} = 0.$$ Note that $D_{p_n}$ is special because it is constant on $(-\infty,t]$ and $(t,\infty]$. So we can write $$ = \frac{1}{n} \sum_{i=1}^n 1(X_i \leq t) \frac{ D_{p_n}(X_i)}{(1 + \varepsilon_n^* D_{p_n}(X_i))} + \frac{1}{n} \sum_{i=1}^n 1(X_i >t) \frac{ D_{p_n}(X_i)}{(1 + \varepsilon_n^* D_{p_n}(X_i))} $$ $$ = \frac{ D_{p_n}(t_-)}{(1 + \varepsilon_n^* D_{p_n}(t_-))} \frac{1}{n} \sum_{i=1}^n 1(X_i \leq t) + \frac{ D_{p_n}(t_+)}{(1 + \varepsilon_n^* D_{p_n}(t_+))} \frac{1}{n} \sum_{i=1}^n 1(X_i >t) = 0$$ where $t_- < t < t_+$. This is equivalent to $$ D_{p_n}(t_-)(1 + \varepsilon_n^* D_{p_n}(t_+)) \frac{1}{n} \sum_{i=1}^n 1(X_i \leq t) + D_{p_n}(t_+)(1 + \varepsilon_n^* D_{p_n}(t_-)) \frac{1}{n} \sum_{i=1}^n 1(X_i >t) = 0$$ which is equivalent to $$ \varepsilon_n^*[D_{p_n}(t_-)D_{p_n}(t_+)]\frac{1}{n} \sum_{i=1}^n 1(X_i \leq t) + \varepsilon_n^*[D_{p_n}(t_-)D_{p_n}(t_+)]\frac{1}{n} \sum_{i=1}^n 1(X_i >t) = - D_{p_n}(t_-) \frac{1}{n} \sum_{i=1}^n 1(X_i \leq t) - D_{p_n}(t_+) \frac{1}{n} \sum_{i=1}^n 1(X_i >t)$$

which is equivalent to $$\varepsilon_n^* = -1 \cdot \frac{ D_{p_n}(t_-) \frac{1}{n} \sum_{i=1}^n 1(X_i \leq t) + D_{p_n}(t_+) \frac{1}{n} \sum_{i=1}^n 1(X_i >t)}{ D_{p_n}(t_-)D_{p_n}(t_+)}.$$

Now note that $D_{p_n}(t_-) = 1 - \Psi(p_n)$ and $D_{p_n}(t_+) = - \Psi(p_n)$ to find $$\varepsilon_n^* = -1 \cdot \frac{ ( 1 - \Psi(p_n)) \frac{1}{n} \sum_{i=1}^n 1(X_i \leq t) - \Psi(p_n)\frac{1}{n} \sum_{i=1}^n 1(X_i >t)}{ D_{p_n}(t_-)D_{p_n}(t_+)}$$ which is equivalent to $$\varepsilon_n^* = \frac{ \frac{1}{n} \sum_{i=1}^n 1(X_i \leq t) - \Psi(p_n)}{(1 - \Psi(p_n))(\Psi(p_n))},$$ which is a very elegant looking solution.

One can show that the TMLE $p_n^*:= p_{n, \varepsilon_n^*}$ satisfies $\frac{1}{n} \sum_{i=1}^n D_{p_n^*}(X_i) = 0$ as desired so that the substitution TMLE $\Psi(p_n^*) = \int 1(x\leq t)p_n^*(x)dx$ is efficient, and in fact $\Psi(p_n^*) = \frac{1}{n} \sum_{i=1}^n 1(X_i \leq t)$ must necessarily hold since $\frac{1}{n} \sum_{i=1}^n D_{p_n^*}(X_i) = 0$. Thus, the TMLE is, like the one-step estimator, equal to the empirical CDF estimator.

Here is R code that implements the TMLE and plots how the initial CDF and density estimates are changed by targeting.

n <- 200
X <- rnorm(n , mean = 0, sd = 1)

#fit_density <- density(X)
#ps <- fit_density$y
xs <- seq(min(X), max(X), length = 1000)
xs <- sort(c(X,xs))

ps <- rep(1/(max(xs) - min(xs)), length(xs))

get_density_preds <- function(grid, p_at_grid, x) {
get_CDF <- function(x, p, t) {
  #x: grid for values of X
  #p: density values at x
  #t: value to get CDF at
  sum(diff(x)*p[-1]*(x<=t)[-1]) # Sloppy Reimann sum integration 
get_CDF <- Vectorize(get_CDF, vectorize.args = "t")

num_steps <- 5
max_eps <- 1
t <- 1
p_i <- ps / get_CDF(xs, ps, max(xs))

# Path in direction of EIF
path <- function(p, EIF, eps) {
  (1 + eps*EIF) * p
 # Plot initial
plot(xs, p_i, main = "Initial density estimate")
plot(xs, get_CDF(xs, ps, xs) , main = "Initial CDF estimate")
for(i in 1:3) {
  Psi <- get_CDF(xs, p_i, t)
  EIF <- (xs<=t) -  Psi #EIF evaluated at grid
  EIF_empirical <- (X<=t) - Psi  #EIF evaluated at observations
  # Check for convergence
  if(i>1 & abs(mean(EIF_empirical)) <= sd(EIF_empirical)/sqrt(n)/log(n)) {
  # loglikelihood risk of density at epsilon along path
  risk_function <- function(eps) {
     p_eps <- path(p_i, EIF, eps)
     p_eps <- p_eps / get_CDF(xs, p_i, max(xs))
     p_obs <- sapply(X, get_density_preds, grid = xs, p_at_grid = p_eps)
    mean(  -log(p_obs))
  # Find risk minimizer along path
  best_eps <- optim(0, risk_function, method = "Brent", lower = -max_eps, upper = max_eps)$par
  analytical_eps <- mean(EIF_empirical)/(1-Psi)/Psi
  print("chosen numerical epsilon and analytical epsilon")
  print(c(best_eps, analytical_eps))
   # Update initial estimator
   p_i <- path(p_i, EIF, analytical_eps)
   # Plot updates
   plot(xs, p_i, main = "Updated density estimates")
  plot(xs, get_CDF(xs, p_i, xs), main = "Updated CDF estimates")

Answer to comments:

TMLE aims to construct an asymptotically unbiased estimate of the ATE by first constructing a path through an initial density estimator that is contained in the statistical model (which you should visually think of as a geometric surface but each point is actually a density function/probability distribution). This path has the property that it points in the direction of greatest change of the ATE (i.e. the direction of steepest ascent/descent). So If I am limited to taking only one small step away from the initial density estimator, a step along this path will give me a new density estimator with the most different ATE estimate. Given this path of steepest change of my target parameter, it makes sense to find the point along this path that maximizes the empirical log likelihood, which gives me the best density fit. Because this path is in the direction of steepest change of my parameter, this additional fitting done by the MLE is targeted towards my parameter of interest. It makes sure that all additional fitting that I do is maximally beneficial for giving me a better estimate for my parameter of interest.

  • 1
    $\begingroup$ By maximizing the empirical log likelihood, you can show that the targeted density estimator will be closer to the true density function than the initial density estimator. As you mentioned, this does not necessarily mean that our ATE estimate is “closer” to the true ATE estimate. $\endgroup$
    – user327671
    Commented Jul 16, 2021 at 22:02
  • 1
    $\begingroup$ For example, it could be that I am estimating the conditional density $P_0(Y=y|A=1,W)$. For the ATE parameter, only its average $E_{P_0}[Y|A=1,W]$ is relevant. If for instance, I chose a path that gave me a better fit of the density $P_0(Y=y|A=1,W)$ but this path does not change the estimate of $E_{P_0}[Y|A=1,W]$, then I would have done a lot of fitting with zero change in my ATE estimate. Because of this, there are many paths through an initial estimator that do not change/improve the ATE estimate at all, even though my density estimator is getting better if I perform MLE along it. $\endgroup$
    – user327671
    Commented Jul 16, 2021 at 22:03
  • 1
    $\begingroup$ This is where the motivation for TMLE kicks in. If a path through the initial density estimator leads to maximal change of the ATE estimate then it must be that this path is maximally changing the ATE-relevant parts of the initial density estimator (e.g. the conditional mean). That is, this path is targeted for the ATE parameter. It then turns out that you can rigorously prove that TMLE actually does debias the initial ATE estimator by utilizing the functional first order expansion of the ATE parameter. $\endgroup$
    – user327671
    Commented Jul 16, 2021 at 22:03
  • 1
    $\begingroup$ Also, when I say true ATE. I mean the parameter that identifies it: $E_WE[Y|A=1,W]-E_WE[Y|A=0,W]$. Causal assumptions only identify a parameter of interest from the observed data. TMLE estimates the estimand associated with that parameter, and does not care about causal assumptions. So TMLE will get you closer to $E_WE[Y|A=1,W]-E_WE[Y|A=0,W]$. If the causal assumptions are not true, then TMLE will not necessarily get you any closer to the causal ATE estimand. $\endgroup$
    – user327671
    Commented Jul 16, 2021 at 22:07
  • 2
    $\begingroup$ So TMLE gives a machine-learning substitution estimator of the ATE that 1. converges quickly (rate n^(-1/2)); 2. we understand how it behaves asymptotically so we can get inference, even though the estimation procedure is blackbox; 3. The limiting distribution is optimal so our confidence intervals have optimal widths; 4. We also have double-robustness so even if you estimate the mean-outcome incorrectly, we can still get consistent estimators if we estimate the propensity score correctly. $\endgroup$
    – user327671
    Commented Jul 17, 2021 at 22:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.