\section{Introduction}

The standard statistical approach of building a model, extracting one or more coefficients and reporting their estimates and associated measures of uncertainty (e.g. confidence intervals) is increasingly being criticised (see e.g. \cite{van2015statistics}).
This standard practice encourages the use of overly simplistic, but misspecified models in order to maintain a simple interpretation of the end result \citep{breiman2001statistical}.
It moreover makes the meaning and definition of the reported coefficients dependent upon the selected models. Inference for such `data-dependent' parameters is not straightforward; ignoring their data-dependent nature, as is commonly done, induces bias, excess variability that is not acknowledged by default standard error estimators and, as a result, overly simplistic inferences.

Building on important results on nonparametric estimation of statistical functionals \citep{pfanzagl1985contributions,pfanzagl1990estimation,bickel1993efficient}, \cite{van_der_laan_targeted_2006}, \cite{van_der_laan_targeted_2011}, \cite{Robins2008} and more recently \cite{chernozhukov_double/debiased_2018}, showed how the aforementioned concerns can be accommodated by centering a statistical analysis around a predefined nonparametric estimand. This is a model-free functional of the observed data distribution which characterises the quantity one wishes to infer from data \citep{Berk2021}. It follows from the existing literature that root-$n$ estimators with well understood asymptotic behaviour can often be derived (under feasible conditions) by making use of the estimand's so-called efficient influence function or canonical gradient under the nonparametric model. The resulting strategies are known as `targeted learning' or 'debiased' machine learning, because they effectively enable the use of data-adaptive estimation strategies to model the data-generating distribution, such as variable selection procedures and machine learning algorithms, whilst permitting valid inference of the estimand of interest.

These developments are quite revolutionary in that they are changing the way in which - we believe - data will be analysed in the future. In particular, they shift the focus from model building and validation to choosing estimands that are well connected to scientific questions of interest \citep{Petersen2014}. This shift enables the analysis to be specified before data is obtained, rather than deciding which statistical quantities to report once a model has been validated, as is usually the case e.g. following model/variable selection. Furthermore, model based analyses usually assume the final model was known a priori, whereas estimand inference based on efficient influence functions tend to be `honest' in the sense of expressing also the uncertainty around selecting the data-generating model, see e.g. \cite{Robins2006} for a precise definition of confidence set `honesty'.

The derivation of the efficient influence function is often regarded as somewhat of a `dark art'. One reason is that it is not given much attention in textbooks on the topic and neither is it given much focus in statistics education. Textbooks that refer to such derivations often rely on a fluency in concepts from functional analysis (e.g. Hilbert Spaces). A further reason is that the majority of research articles that derive the efficient influence function of a statistical estimand, rely on manipulating a derivative expression into a canonical form, as the integral of a product of an efficient influence function and a score function. These derivations are often complicated, with some steps appearing as if from nowhere to achieve the desired form.

In this tutorial paper, we instead advocate an equivalent approach based on Gateaux derivatives, formalised by \cite{ichimura2015influence}, which is much simpler in our opinion. We will explain this approach and show how to make use of it, while also providing intuitive insight into what an efficient influence function is. We will moreover explain how root-$n$ converging statistical/machine-learning-based estimators can be constructed, using the efficient influence function, and what conditions are needed for these to work well. This tutorial obeys the principles of van der Laan's `roadmap' \citep{van_der_laan_targeted_2011}. It is aimed to be broadly accessible to students and researchers who would like to derive efficient influence functions for all sorts of nonparametric estimands, using simple differentiation methods, such as the chain rule. We use diverse examples first to show the steps in calculating the efficient influence function (Section \ref{CalcIF}), and also to convey the very broad applicability of the theory (Section \ref{derivations}).

\section{Step 1: Defining the estimand of interest}

The starting point of most statistical analyses is a (semi)parametric model, which is then often interpreted as representing how nature has generated the data.
For certain applications, such as in the physical sciences, this model can be the result of a deep theoretical understanding of the data-generating mechanism. However, oftentimes, especially in the spheres of medicine, psychology and economics, the model is chosen for its simplicity and convenience. Many ubiquitous models, such as the generalized linear and Cox proportional hazards models, are commonly used without reference to a mechanistic understanding, rather because the parameters indexing those models provide useful summaries of associations that are of interest to the analysis. This is problematic for various reasons.
First, nature is rarely as simple as we would like it to be. This leaves many data analysts torn between reporting a simple model, which is likely misspecified, versus reporting a complex model, which is difficult to interpret \citep{breiman2001statistical}.
It demands choosing between an analysis result that is likely biased (as a result of model misspecification) versus one that is likely useless (in view of its complexity).
Second, standard statistical theory for (semi)parametric models was developed for settings where the model is a priori justified by some biological, economic, ... theory (so that one can assume it to be correct)
and where moreover the data analyst commits to using that model. The truth is that a given model is rarely known to be correct, and that data analysts therefore do not commit to a single model, by adopting model selection strategies. This invalidates standard statistical theory.
Third, even the common attempt to infer the model from data (for instance, by relying on variable selection strategies) is overly ambitious as many competing models often fit the data nearly equally well \citep{breiman2001statistical}. While this is generally well realised, it is also then systematically `forgotten' in how we report and interpret statistical analysis results.

To accommodate these concerns, we will instead aim to infer so-called nonparametric estimands. These are functionals of the true observed data distribution $\mathcal{P}$, which are well defined without reference to a (semi)parametric model, and target the scientific question of interest.
With interest in the mean outcome $Y$, such estimand is unambiguously defined as
\[\Psi_1(\mathcal{P})=E_{\mathcal{P}}(Y),\]
where the subscript $\mathcal{P}$ explicates that the expectation $E_{\mathcal{P}}$ is calculated w.r.t. the true distribution $\mathcal{P}$ of $Y$.
We will equivalently write this as
\[\Psi_1(\mathcal{P})=\mathcal{P}(Y) = \int y d\mathcal{P}(y).
\]
where $d\mathcal{P}(y)$ denotes integration w.r.t. to the probability measure $\mathcal{P}$ for the random variable $Y$. When $Y$ is continuous, $d\mathcal{P}(y)$ in this expression can be replaced with $f(y)dy$ to recover the Riemann integral over the probability density function of $Y$. For many of the examples in this paper we work with Riemann integrals, which are likely to be familiar to most readers.

As a second example, suppose we are interested in the effect of a dichotomous exposure $X$ (coded 0 or 1) on an outcome $Y$ in the presence of data on a possibly high-dimensional vector of covariates $Z$ that is sufficient to adjust for confounding. Then a relevant (statistical) estimand could be defined as
\[\Psi_2(\mathcal{P})=E_{\mathcal{P}}\left\{E_{\mathcal{P}}(Y|X=1,Z)-E_{\mathcal{P}}(Y|X=0,Z)\right\},\]
where, with a slight abuse of notation, the subscript $\mathcal{P}$ now explicates that the expectation $E_{\mathcal{P}}$ is calculated w.r.t. the true distribution $\mathcal{P}$ of $(Z,X,Y)$. This is known in the causal inference literature as the \emph{average causal effect} or \emph{average treatment effect}.

Alternatively, regardless of whether the exposure is dichotomous or not, its effect on $Y$ can also be expressed using the estimand
\[\Psi_3(\mathcal{P})=\frac{E_{\mathcal{P}}\left[\left\{X-E_{\mathcal{P}}(X|Z)\right\}Y\right]}{E_{\mathcal{P}}\left[\left\{X-E_{\mathcal{P}}(X|Z)\right\}^2\right]},\
which equals the expected conditional covariance between $X$ and $Y$, given $Z$, divided by the expected conditional variance of $X$, given $Z$.
Where we are happy to assume a partially linear model for $E_{\mathcal{P}}(Y|X,Z)$, such as
\[E_{\mathcal{P}}(Y|X,Z)=\beta X+\omega(Z),\]
for some function $\omega(.)$, $\Psi_3(\mathcal{P})$ reduces to $\beta$, but remains well defined outside this model \citep{Robins2008,Vansteelandt2020}.

\section{Step 2: Calculate the estimand's efficient influence function}\label{Sec:IF}

\subsection{Preliminaries}

Throughout, we will assume that we have access to i.i.d.~ observed data on $O_i\equiv (Z_i,X_i,Y_i)$ for subjects $i=1,...,n$.
An estimator of the above estimands is then readily obtained by substituting $\mathcal{P}$ by an estimator $\hat{\mathcal{P}}_n$, where the sub-index $n$ denotes the sample size.
For instance, choosing $\hat{\mathcal{P}}_n$ to equal the empirical distribution, $\mathcal{P}_n$, of the observations $Y_1,...,Y_n$ gives rise to the empirical `plug-in' estimator
\[\Psi_1(\hat{\mathcal{P}}_n) = \hat{\mathcal{P}}_n(Y)=\frac{1}{n}\sum_{i=1}^nY_i.\]

For $\Psi_2(\mathcal{P})$, let $\hat{\mathcal{P}}_n$ be any distribution of $(Z,X,Y)$ such that the marginal distribution of $Z$ is given by its empirical distribution, and that the conditional distribution of $Y$ given $X=x$ for $x=0,1$ and $Z=Z_i$ for $i=1,...,n$ has conditional mean equal to a given estimator $\hat{E}(Y|X=x,Z=Z_i)$, such as the prediction from some machine learning algorithm. Then
\[\Psi_2(\hat{\mathcal{P}}_n)=\frac{1}{n}\sum_{i=1}^n \hat{E}(Y|X=1,Z=Z_i)-\hat{E}(Y|X=0,Z=Z_i).\]
Finally, for $\Psi_3(\mathcal{P})$, let $\hat{\mathcal{P}}_n$ be any distribution of $(Z,X,Y)$ such that the conditional distribution of $X$ given $Z=Z_i$ for $i=1,...,n$ has conditional mean equal to a given estimator $\hat{E}(X|Z=Z_i)$,
and that the marginal distribution of $X-\hat{E}(X|Z)$ and $\left\{X-\hat{E}(X|Z)\right\}Y$ is given by its empirical distribution. Then
\[\Psi_3(\hat{\mathcal{P}}_n)=\frac{\sum_{i=1}^n \left\{X_i-\hat{E}(X|Z=Z_i)\right\}Y_i}{\sum_{i=1}^n \left\{X_i-\hat{E}(X|Z=Z_i)\right\}^2} .\]
The key question now is whether $\Psi(\hat{\mathcal{P}}_n)$ is a good proxy for $\Psi(\mathcal{P})$. To understand this, we will scale their difference by $\sqrt{n}$.
When this scaled difference converges in distribution (to a non-degenerate law), then we can roughly say that $\Psi(\hat{\mathcal{P}}_n)$ differs from $\Psi(\mathcal{P})$ up to a term of the order 1 over root-$n$. We then say that
$\Psi(\hat{\mathcal{P}}_n)$ converges to $\Psi(\mathcal{P})$ at parametric rate, or root-$n$ rate, which is usually the best that we can hope to achieve.

For the sample mean $\Psi_1$, we have that $\hat{\mathcal{P}}_n=\mathcal{P}_n$ so that this scaled difference equals
\begin{align*}
\sqrt{n}\left\{\Psi_1(\hat{\mathcal{P}}_n) - \Psi_1(\mathcal{P})\right\} &= \sqrt{n}(\hat{\mathcal{P}}_n - \mathcal{P})Y \\
&= \sqrt{n}\hat{\mathcal{P}}_n(Y-\Psi_1) = \sqrt{n}\mathcal{P}_n(Y-\Psi_1) \\
&= \frac{1}{\sqrt{n}} \sum_{i=1}^n (Y_i - \mu) \overset{d}{\to} \n{0}{\sigma^2},
\end{align*}
by the classical central limit theorem, where $\mu$ and $\sigma^2$ are the mean and variance of $Y$, respectively.

We are lucky here that the difference $\sqrt{n}\left\{\Psi_1(\hat{\mathcal{P}}_n) - \Psi_1(\mathcal{P})\right\}$ can be written in terms of the operator $\sqrt{n}(\mathcal{P}_n - \mathcal{P})$ applied to some $Y$, but this is not generally the case, for the following reasons.
First, the difference $\sqrt{n}\left\{\Psi(\hat{\mathcal{P}}_n) - \Psi(\mathcal{P})\right\}$ will in general depend on how much $\hat{\mathcal{P}}_n$ differs from $\mathcal{P}$. This was easy in the above example, where $\hat{\mathcal{P}}_n$ refers to the empirical distribution $\mathcal{P}_n$ of $Y$, whose behaviour is easy to understand. It is much harder in more general cases where $\hat{\mathcal{P}}_n$ may involve data-adaptive estimators, such as predictions $\hat{E}(Y|X=x,Z=Z_i)$ or $\hat{E}(X|Z=Z_i)$ obtained via machine learning or via parametric model building procedures. For such predictions, we may at best have access to some overall, marginal measure of prediction error, but will often have a poor understanding of the bias and imprecision in these predictions at specific covariate levels $Z_i$.  Second, the difference $\sqrt{n}\left\{\Psi(\hat{\mathcal{P}}_n) - \Psi(\mathcal{P})\right\}$ will in general also depend on how sensitive the estimand $\Psi(.)$ is to changes in the data-generating distribution. This is also generally poorly understood given that $\mathcal{P}$ indexing $\Psi(\mathcal{P})$ is an infinite-dimensional parameter (apart from exceptional cases where the observed data is discrete).

The situation thus looks a bit hopeless at this stage, and indeed, we will not succeed to understand the difference $\sqrt{n}\left\{\Psi(\hat{\mathcal{P}}_n) - \Psi(\mathcal{P})\right\}$ for arbitrary estimators $\hat{\mathcal{P}}_n$ and arbitrary estimands $\Psi(\mathcal{P})$. However, we will see that progress can be made for specific estimators $\hat{\mathcal{P}}_n$, and for estimands $\Psi(\mathcal{P})$ that are sufficiently smooth in the data-generating law $\mathcal{P}$. Before proceeding, we will first formalise the right level of smoothness that is needed.

\subsection{Parametric submodels}\label{sec:param_sub}

To understand how sensitive $\Psi(.)$ is to changes in the data-generating distribution, we will first take a step back.
Rather than examining how $\Psi(.)$ changes as we slightly perturb $\mathcal{P}$ towards $\hat{\mathcal{P}}_n$, we will study the effect of such perturbation in the direction of a fixed, deterministic distribution, $\widetilde{\mathcal{P}}$, which, for the purpose of this discussion, we shall assume is absolutely continuous with respect to $\mathcal{P}$ (i.e. the support of $\widetilde{\mathcal{P}}$ is contained in the support of $\mathcal{P}$).
There are many ways in which we may change $\Psi(.)$ to $\widetilde{\mathcal{P}}$.
Here, we will focus on perturbations in the direction parameterised via the one-dimensional mixture model
\begin{equation}\label{submodel}
\mathcal{P}_t = t\widetilde{\mathcal{P}}+(1-t)\mathcal{P},
\end{equation}
indexed by $t\in[0,1]$, which is called a \emph{parametric submodel}. This is not a parametric model in the usual sense (given that the true data-generating law $\mathcal{P}$ is unknown), but is used here as a convenient tool to formalise small perturbations away from $\mathcal{P}$ in the direction of
$\widetilde{\mathcal{P}}$. In particular, note that $\mathcal{P}_0 = \mathcal{P}$ and $\mathcal{P}_1 = \widetilde{\mathcal{P}}$.

The sensitivity of $\Psi(.)$ to changes in the data-generating distribution in the direction of $\widetilde{\mathcal{P}}$ can now be formalised in terms of the pathwise or directional derivative,
$$\lim _{t\downarrow 0}\left(\frac{\Psi\left(\mathcal{P}_{t}\right)-\Psi(\mathcal{P})}{t}\right) = \frac{d  \Psi(\mathcal{P}_t)}{d  t}\Big|_{t=0},$$
evaluated at $t=0$,
and in the direction of $\widetilde{\mathcal{P}}$.
When this limit exists (i.e., is finite), it is called a G\^ateaux derivative. This generalises the concept of a directional derivative to functional analysis, describing how to take the derivative of a function with respect to a function.
Informally, when this derivative `exists' for all regular parametric submodels, then we will say that the estimand is pathwise differentiable.
Here, a `regular' parametric submodel is such that its score $\widetilde{\mathcal{P}}(O)/\mathcal{P}(O)-1$ has finite variance, a mild restriction that will be needed to ensure that the derivative $d\mathcal{P}_t/d  t|_{t=0}$ is (or more precisely, inner products with this score are) well-defined.
In the next paragraph, we will formalise this definition of pathwise differentiability. This formalisation will be practically useful, as it will provide insight what the so-called efficient influence function (also referred to as \textit{canonical gradient}, or \textit{influence curve}) is, how it can be calculated, and why it is useful.

As in \cite{Fisher2020}, we develop some intuition by first considering the special case of discrete data $O$ with support $\{o_1,...,o_k\}$. Then
\[\frac{d \Psi(\mathcal{P}_t) }{d t}\Big|_{t=0}=\sum_{j=1}^k
\frac{d \Psi(\mathcal{P}_t)}{d \mathcal{P}_t(o_j)}\Big|_{t=0}\frac{d \mathcal{P}_t(o_j)}{d t}\Big|_{t=0}.\]
Here, $d\Psi(\mathcal{P}_t)/d\mathcal{P}_t(o_j)$ expresses the estimand's sensitivity to small changes in the observed data law.
The second component expresses how the observed data law changes along the considered path. It is easily verified to equal
\[\frac{d \mathcal{P}_t(o_j)}{d t}\Big|_{t=0}=\widetilde{\mathcal{P}}(o_j)-\mathcal{P}(o_j).\]
The resulting identity
\begin{equation}\label{discrete}
\frac{d\Psi(\mathcal{P}_t)}{d t}\Big|_{t=0}=\sum_{j=1}^k
\frac{d \Psi(\mathcal{P}_t)}{d \mathcal{P}_t(o_j)}\Big|_{t=0}\left\{\widetilde{\mathcal{P}}(o_j)-\mathcal{P}(o_j)\right\},\end{equation}
is limiting (by being focussed on discrete data) and ignores that the probabilities $\mathcal{P}_t(o_1),...,\mathcal{P}_t(o_k)$ are not variation-independent (i.e., they sum to 1 and thus cannot be changed in arbitrary ways) \citep{Fisher2020}.

We therefore appeal to Riesz's representation theorem,
according to which this derivative, when it exists, can be obtained via integration of a unique `representer' $\phi(O,\mathcal{P})$ with finite variance under $\mathcal{P}$, w.r.t. some measure:

\begin{eqnarray}
\frac{d  \Psi(\mathcal{P}_t)}{d  t}\Big|_{t=0}
&=& \int \phi(o,\mathcal{P})\left\{d\widetilde{\mathcal{P}}(o)-d\mathcal{P}(o)\right\}=(\widetilde{\mathcal{P}}-\mathcal{P}) \{\phi(O,\mathcal{P})\}. \label{IF_Rieszl}
\end{eqnarray}
Contrasting identity \eqref{IF_Rieszl} with \eqref{discrete}, we learn that the representer $\phi(O,\mathcal{P})$ is a functional derivative which characterises how sensitive the estimand $\Psi(\mathcal{P})$ is to changes in the data-generating distribution $\mathcal{P}$.
It is referred to as the estimand's canonical gradient, efficient influence curve or efficient influence function (under the nonparametric model).
The existence of a representer with finite variance such that \eqref{IF_Rieszl} holds, essentially expresses that the estimand is sufficiently smooth as a functional of the data-generating law (so that the notion of a `derivative' is well-defined); here, the finite-variance condition expresses that the `derivative' of the estimand w.r.t. the data-generating distribution is finite. Since identity (\ref{IF_Rieszl}) is insensitive to constant, additive shifts in $\phi(O,\mathcal{P})$, we will henceforth limit ourselves to mean zero functions (under $\mathcal{P}$) without loss of generality.

We can now more formally define the estimand to be pathwise differentiable when there exists a mean-zero, finite-variance function $\phi(O,\mathcal{P})$ which satisfies \eqref{IF_Rieszl}  for all (regular) parametric submodels. Since the efficient influence function has mean zero, $\mathcal{P} \{\phi(O,\mathcal{P})\}=0$, the derivative of $\Psi(\mathcal{P}_t)$ w.r.t. $t$ can equivalently be represented as the average of the efficient influence function over the distribution $\widetilde{\mathcal{P}}$:
\begin{eqnarray}
\frac{d  \Psi(\mathcal{P}_t)}{d  t}\Big|_{t=0} &=& \widetilde{\mathcal{P}} \left\{ \phi(O,\mathcal{P})\right\} = E_{\widetilde{\mathcal{P}}}\left\{\phi(O,\mathcal{P})\right\}. \label{Riesz_non_canonical}
\end{eqnarray}
This result forms the basis of how we will calculate
the efficient influence function of an estimand.

\subsection{How to calculate the efficient influence function of an estimand}\label{CalcIF}

There are several ways to derive efficient influence functions.
We here advocate the ``point mass contamination'' strategy that we find simplest. In particular, we will perturb the estimand it in the direction $\widetilde{\mathcal{P}}$ of a point mass at single observation $\tilde{o}$. Identity \eqref{Riesz_non_canonical} then gives the efficient influence function at observation $o$ directly as
\[\phi(o,\mathcal{P})=\frac{d  \Psi(\mathcal{P}_t)}{d  t}\Big|_{t=0},\]
where the right-hand side is a so-called G\^ateaux derivative. This has the same properties as ordinary derivatives, familiar from calculus, such as the chain rule. This will facilitate calculations. The following examples illustrate this.

Throughout, for convenience, we will implicitly assume that we work with continuous variables, but the results continue to hold for discrete variables, or a mix of discrete and continuous variables, upon swapping sums with integrals, indicators with Dirac delta functions, and probability mass functions with probability density functions, where needed. For our purposes, $\mathbbm{1}_{\tilde{o}}(o)$ denotes the Dirac delta function w.r.t. $\tilde{o}$; i.e., the density of an idealized point mass at $\tilde{o}$, which equals zero everywhere except at $\tilde{o}$ and which integrates to 1.

\paragraph{Example 1 (population mean).} As a first, simple example, consider the  mean of $Y$:
\begin{align*}
\Psi(\mathcal{P}) &= \mathcal{P}(Y)=E_{\mathcal{P}}(Y)= \int yf(y) dy,
\end{align*}
where $f(y)$ denotes the density function of $Y$ under $\mathcal{P}$ (which we assume to be absolutely continuous w.r.t. the Lebesgue measure, though results hold more generally).
Perturbing in the direction of a single observation $\tilde{y}$,
\begin{align*}
f_t(y) &
= t\mathbbm{1}_{\tilde{y}}(y) + (1-t)f(y),
\end{align*}
one obtains,
\begin{align*}
\Psi(\mathcal{P}_t) &= t\int y\mathbbm{1}_{\tilde{y}}(y) dy + (1-t)E_{\mathcal{P}}(Y)= t\tilde{y} + (1-t)\Psi(\mathcal{P}).
\end{align*}
By the chain rule, taking a derivative with respect to $t$ at $t=0$, gives
\begin{align*}
\frac{d  \Psi(\mathcal{P}_t)}{d  t}\Big|_{t=0} &= \tilde{y} - \Psi(\mathcal{P}).
\end{align*}
Because this has finite variance, we conclude that $E(Y)$ is pathwise differentiable with efficient influence function $Y - \Psi(\mathcal{P})$. $\Box$

\paragraph{Example 2 (density at a point $y$).} Consider next the density at a given value $y$, $\Psi(\mathcal{P}) =f(y)$. Under the parametric submodel of Example 1, we readily find that
\begin{align*}
\Psi(\mathcal{P}_t) &= t\mathbbm{1}_{\tilde{y}}(y) + (1-t)f(y).
\end{align*}
By the chain rule, taking a derivative with respect to $t$ at $t=0$, gives
\begin{align*}
\frac{d  \Psi(\mathcal{P}_t)}{d  t}\Big|_{t=0} &= \mathbbm{1}_{\tilde{y}}(y)- \Psi(\mathcal{P}).
\end{align*}
Because the Dirac delta function is unbounded when $Y$ is absolutely continuous w.r.t. Lebesgue measure, and therefore has infinite variance, we conclude that $f(y)$ is not pathwise differentiable.
This lack of smoothness is the result of insufficient information in the data on the density $f(y)$ at the single point $y$. It generally implies that no root-$n$ converging estimators can be constructed.
 $\Box$

\paragraph{Example 3 (average density).} Consider next the average density of $Y$:
\begin{align*}
\Psi(\mathcal{P}) &= \mathcal{P}\left\{f(Y)\right\}=E_{\mathcal{P}}\left\{f(Y)\right\} = \int f^2(y) dy.
\end{align*}
Under the parametric submodel of Example 1,
\begin{align*}
\Psi(\mathcal{P}_t) &= \int f_t^2(y) dy.
\end{align*}
By the chain rule, taking a derivative with respect to $t$ at $t=0$, gives
\begin{align*}
\frac{d  \Psi(\mathcal{P}_t)}{d  t}\Big|_{t=0} &= \int 2f(y)\frac{d}{d  t}f_t(y)\Big|_{t=0}  dy \\
&= 2 \int f(y)\left\{\mathbbm{1}_{\tilde{y}}(y)-f(y)\right\} dy \\
&= 2\left\{f(\tilde{y}) - \Psi(\mathcal{P})\right\}.
\end{align*}
Here, we use that the Dirac delta function $\mathbbm{1}_{y}(Y)$ has average
\begin{align*}
\mathcal{P} \mathbbm{1}_{y}(Y)&= \int \mathbbm{1}_{y}(\tilde{y}) f(\tilde{y})d\tilde{y} = f(y),
\end{align*}
equal to the density at $y$ under the law $\mathcal{P}$.
Since $2\left\{f(Y) - \Psi(\mathcal{P})\right\}$ has finite variance, we conclude that $E_{\mathcal{P}}\left\{f(Y)\right\}$ is pathwise differentiable with efficient influence function $2\left\{f(Y) - \Psi(\mathcal{P})\right\}$. $\Box$

When perturbing the density $f(o)$ of a vector of observations $O$ in the direction of a point mass at $\tilde{o}$, we have the identity
\[\frac{d f_t(o) }{d  t}\Big|_{t=0} = \mathbbm{1}_{\tilde{o}}(o)-f(o),\]
which was also used in Example 2.
This implies a simple formula for the efficient influence function at $\tilde{o}$ of the estimand $E_{\mathcal{P}}\left\{g(O,\mathcal{P})\right\}$ for some function $g(O,\mathcal{P})$ of $O$ and the true distribution:
\begin{eqnarray}
\frac{d}{d  t}E_{\mathcal{P}_t}\left\{g(O,\mathcal{P}_t)\right\}\Big|_{t=0}
&=&\frac{d}{d  t}\left\{\int g(o,\mathcal{P}_t) f_t(o) do \right\}\Big|_{t=0}\nonumber\\
&=&\left\{\int \frac{d}{d  t}g(o,\mathcal{P}_t) f_t(o) do+\int g(o,\mathcal{P}_t) \frac{d}{d  t}f_t(o) do \right\}\Big|_{t=0}\nonumber\\
&=&E_{\mathcal{P}}\left\{\frac{d}{d  t}g(o,\mathcal{P}_t)\right\}\Big|_{t=0}+g(\tilde{o},\mathcal{P})- E_{\mathcal{P}}\left\{g(O,\mathcal{P})\right\}.\label{expectation}
\end{eqnarray}
We apply this general identity in the following example.

\paragraph{Example 4 (covariance).}
The covariance
\[\Psi(\mathcal{P})= E_{\mathcal{P}}\left[\left\{Y-E_{\mathcal{P}}(Y)\right\}\left\{X-E_{\mathcal{P}}(X)\right\}\right],\]
can be written as $E_{\mathcal{P}}\left\{g(O,\mathcal{P})\right\}$ for $O\equiv (X,Y)$ and $g(o,\mathcal{P}) = \left\{y-E_{\mathcal{P}}(Y)\right\}\left\{x-E_{\mathcal{P}}(X)\right\}$.
Using \eqref{expectation}, we thus find that
\[\frac{d  \Psi(\mathcal{P}_t)}{d  t}\Big|_{t=0}  =
 E_{\mathcal{P}}\left[\frac{d}{d  t}
 \left\{Y-E_{\mathcal{P}_t}(Y)\right\}\left\{X-E_{\mathcal{P}_t}(X)\right\}\Big|_{t=0}\right]+\left\{\tilde{y}-E_{\mathcal{P}}(Y)\right\}\left\{\tilde{x}-E_{\mathcal{P}}(X)\right\}-\Psi(\mathcal{P}).\]
By the chain rule, we further have that
\[\frac{d}{d  t}\left\{Y-E_{\mathcal{P}_t}(Y)\right\}\left\{X-E_{\mathcal{P}_t}(X)\right\}\Big|_{t=0}
= - \frac{d}{d  t}E_{\mathcal{P}_t}(Y)\Big|_{t=0}\left\{X-E_{\mathcal{P}}(X)\right\}- \frac{d}{d  t}E_{\mathcal{P}_t}(X)\Big|_{t=0}\left\{Y-E_{\mathcal{P}}(Y)\right\}.\]
Further applying \eqref{expectation} to $E_{\mathcal{P}_t}(Y)$ and $E_{\mathcal{P}_t}(X)$, we find that
\[\frac{d}{d  t}\left\{Y-E_{\mathcal{P}_t}(Y)\right\}\left\{X-E_{\mathcal{P}_t}(X)\right\}\Big|_{t=0}
= - \left\{\tilde{y}-E_{\mathcal{P}}(Y)\right\}\left\{X-E_{\mathcal{P}}(X)\right\}- \left\{\tilde{x}-E_{\mathcal{P}}(X)\right\}\left\{Y-E_{\mathcal{P}}(Y)\right\},\]
which has mean zero. Since $\left\{Y-E_{\mathcal{P}}(Y)\right\}\left\{X-E_{\mathcal{P}}(X)\right\}$ has finite variance, we conclude that the covariance $\Psi(\mathcal{P})$ is pathwise differentiable with efficient influence function
\[\left\{Y-E_{\mathcal{P}}(Y)\right\}\left\{X-E_{\mathcal{P}}(X)\right\}-\Psi(\mathcal{P}).\] $\Box$

\paragraph{Example 5 (potential outcome mean).}

Let $Y^x$ denote the potential outcome under exposure level $x$, which expresses what value the outcome of a given individual would have taken had his/her exposure been set to $x$ by some intervention. Under the usual identifying assumptions, (positivity, consistency, non interference and conditional exchangeability given $Z$) \citep{hernan2006estimating},
\begin{align*}
\Psi(\mathcal{P})&=E_{\mathcal{P}}\left\{E_{\mathcal{P}}(Y|X=1,Z)\right\}
\end{align*}
is a statistical estimand of the population mean of $Y^1$.

Perturbing $\mathcal{P}$ in the direction of a point mass at $(\tilde{z},\tilde{x},\tilde{y})$, we find that
\begin{align*}
\Psi(\mathcal{P}_t)&= \int y f_t(y|1,z)f_t(z) dydz \\
&= \int y \frac{f_t(y,1,z)f_t(z)}{f_t(1,z)} dydz,
\end{align*}
where $f_t(y|x,z)$ is the conditional density function of $Y$, given $X=x,Z=z$, under the parametric submodel, and  $f_t(y,x,z), f_t(x,z),$ and $f_t(z)$ are the joint density functions of $(Y,X,Z), (X,Z),$ and $Z$, respectively under the parametric submodel. By the chain rule, we thus have that
\begin{align*}
\frac{d  \Psi(\mathcal{P}_t)}{d  t}\Big|_{t=0} =&  \int y \Bigg\{ \frac{f(z)}{f(1,z)}\frac{d}{d  t} f_t(y,1,z)\Big|_{t=0}- \frac{f(y,1,z)f(z)}{f(1,z)^2} \frac{d}{d  t} f_t(1,z)\Big|_{t=0}+  \frac{f(y,1,z)}{f(1,z)} \frac{d}{d  t} f_t(z)\Big|_{t=0} \Bigg\} dydz\\
 =& \int     y \frac{f(y,1,z)f(z)}{f(1,z)} \Bigg( \frac{\mathbbm{1}_{\tilde{y},\tilde{x},\tilde{z}} (y,1,z)}{f(y,1,z)} - \frac{\mathbbm{1}_{\tilde{x},\tilde{z}} (1,z)}{f(1,z)} + \frac{\mathbbm{1}_{\tilde{z}} (z)}{f(z)} - 1   \Bigg) dydz.
\end{align*}
Evaluating the integral gives the canonical gradient of $\Psi(\mathcal{P})$ at $(\tilde{z},\tilde{x},\tilde{y})$:
\begin{align*}
\frac{d  \Psi(\mathcal{P}_t)}{d  t}\Big|_{t=0} =\frac{\mathbbm{1}_{\tilde{x}}(1)}{\pi(\tilde{z},\mathcal{P})} \left\{\tilde{y} - m_1(\tilde{z},\mathcal{P})\right\} + m_1(\tilde{z},\mathcal{P}) - \Psi(\mathcal{P}),
\end{align*}
where $m_1(z,\mathcal{P})=E_{\mathcal{P}}(Y|X=1,Z=z)$  and $\pi(z,\mathcal{P}) = f(1|z) = E_{\mathcal{P}}(X|Z=z)$ is the propensity score.
We conclude that $\Psi(\mathcal{P})$ is pathwise differentiable with the above efficient influence function.

From this, it readily follows that
$\Psi_2(\mathcal{P})$ (the average treatment effect) is pathwise differentiable with the efficient influence function given by
\begin{align*}
\varphi_{1}(O,\mathcal{P}) - \varphi_{0}(O,\mathcal{P}) - \Psi_2(\mathcal{P})
\end{align*}
where $\varphi_{x}(O,\mathcal{P})$ is the `uncentered' efficient influence curve
\begin{align}
\varphi_{x}(O,\mathcal{P}) = \frac{\mathbbm{1}_{X}(x)}{f(x|Z)} \left\{Y - m(x,Z)\right\} + m(x,Z). \label{AIPW_infl}
\end{align} $\Box$

\paragraph{Example 6 (conditional outcome mean).}
Before moving on to more elaborate examples, we finally consider
\begin{align*}
\Psi(\mathcal{P})&=E_{\mathcal{P}}(Y|X=x),
\end{align*}
for a given value $x$, where $X$ may be (absolutely) continuous (w.r.t. Lebesgue measure).
Perturbing $\mathcal{P}$ in the direction of a point mass at $(\tilde{x},\tilde{y})$, we find that
\begin{align*}
\Psi(\mathcal{P}_t)&= \int y \frac{f_t(y,x)}{f_t(x)} dy,
\end{align*}
where $f_t(y,x)$ and $f_t(x)$ are the joint density functions of $(Y,X)$ and $X$, respectively under the parametric submodel. By
the chain rule, we thus have that the canonical gradient is
\begin{align}
\phi(\tilde{o},\mathcal{P}) = \frac{d  \Psi(\mathcal{P}_t)}{d  t}\Big|_{t=0} =&  \int y \Bigg\{ \frac{1}{f(x)}\frac{d}{d  t} f_t(y,x)\Big|_{t=0}- \frac{f(y,x)}{f(x)^2} \frac{d}{d  t} f_t(x)\Big|_{t=0}\Bigg\} dy \nonumber \\
 =& \int    \left[ \frac{y}{f(x)} \left\{ \mathbbm{1}_{\tilde{y},\tilde{x}} (y,x)-f(y,x)\right\} - \frac{yf(y,x)}{f(x)^2} \left\{\mathbbm{1}_{\tilde{x}} (x)-f(x)\right\}\right]dy \nonumber \\
 =& \frac{\mathbbm{1}_{\tilde{x}} (x)}{f(x)} \left\{\tilde{y}-E_{\mathcal{P}}(Y|X=x)\right\}. \label{cond_result}
\end{align}
An issue, however, emerges when one considers the variance of the influence function.
\begin{align*}
\mathrm{var}\left\{\phi(O,\mathcal{P})\right\} &= \int \left(\frac{\mathbbm{1}_{\tilde{x}} (x)}{f(x)}\right)^2 \left\{\tilde{y}-E_{\mathcal{P}}(Y|X=x)\right\}^2 f(\tilde{y}|\tilde{x}) f(\tilde{x}) d\tilde{y} d\tilde{x} \\
&= \frac{\mathbbm{1}_{x} (x)}{f(x)} \int \left\{\tilde{y}-E_{\mathcal{P}}(Y|X=x)\right\}^2 f(\tilde{y}|x) d\tilde{y} \\
&= \frac{\mathbbm{1}_{x} (x)}{f(x)} \mathrm{var}(Y|X=x)
\end{align*}
Since the Dirac delta function $\mathbbm{1}_{x}(x)$ takes an infinitely large value when $X$ is continuous (i.e. when its probability distribution is absolutely continuous w.r.t. Lebesgue measure), we conclude that the conditional mean is not pathwise differentiable in that case.

When $X$ is discrete (as in Example 5), however, then we have that the indicator function $\mathbbm{1}_{x} (x)=1$, so that the variance of the efficient influence function is finite (so long as $\mathrm{var}(Y|X=x) < \infty$ and $f(x)>0$). $\Box$

The approach that we have adopted in the above examples follows the calculation in \cite{hampel1974influence}.
A second, perhaps more common approach, instead uses the following, canonical form
\begin{eqnarray}
\frac{d  \Psi(\mathcal{P}_t)}{d  t}\Big|_{t=0} &=& \int \phi(o,\mathcal{P})\left\{d\widetilde{\mathcal{P}}(o)-d\mathcal{P}(o)\right\}\\
&=& \int \phi(o,\mathcal{P})\left(\frac{d\widetilde{\mathcal{P}}(o)}{d\mathcal{P}(o)}-1\right)d\mathcal{P}(o)  \nonumber\\
&=& \int \phi(o,\mathcal{P})S(o)d\mathcal{P}(o)  \nonumber\\
&=&E_{\mathcal{P}} \left\{ \phi(O,\mathcal{P}) S(O)\right\}=\mathcal{P} \left\{ \phi(O,\mathcal{P}) S(O)\right\}. \label{IF_Riesz}
\end{eqnarray}
The efficient influence function is then calculated as the unique mean zero function $\phi(O,\mathcal{P})$ whose inner product (i.e., covariance) with the score $S(O)$ under a parametric submodel $\mathcal{P}_t$ equals the pathwise derivative $d  \Psi(\mathcal{P}_t)/d  t |_{t=0}$, for all parametric submodels; see \cite{Levy2019} for a tutorial. This can be quite laborious, however, since one must manipulate score functions and integral expressions, and moreover solve a functional equation like (\ref{IF_Riesz}) \citep{ichimura2015influence}.

The latter approach nonetheless appears more commonly used because it lends itself easier to semiparametric modelling, where the scores $S(O)$ can now be confined to the scores of those parametric submodels that obey the semiparametric model restrictions.
A further reason for the greater popularity of this approach may be the apparent limitation of the approach advocated in this tutorial, that certain estimands (e.g., example 3) cannot be evaluated at $\mathcal{P}_t$ because of the use of Dirac delta functions.
\cite{ichimura2015influence} note that this does not invalidate the approach, as it can be resolved by substituting the  Dirac delta function in $\mathcal{P}_t$ by a probability measure, indexed by a bandwidth $h$, that approaches a point mass when the bandwidth converges to 0. This modification justifies the approach that we adopt, but for simplicity it will be left implicit in the remainder of the work.

\section{Step 3: construct an estimator based on the estimand's efficient influence function}

\subsection{Plug-in bias and how to remove it}

The previous results help us to develop insight into the scaled difference
\begin{equation}\label{scalediff}\sqrt{n}\left\{\Psi(\widetilde{\mathcal{P}}) - \Psi(\mathcal{P})\right\}.\end{equation}
In particular, the canonical gradient gave us a way to express the notion of a functional derivative of the estimand w.r.t. directional changes in the data-generating law. This in turn forms the basis of a functional analog to the Taylor expansion, the so-called von Mises expansion, which is essentially derived from the Taylor series expansion of $\Psi(\mathcal{P}_t)$ about the point $t=1$ in the one-dimensional parametric submodel.
\begin{eqnarray*}
\Psi(\mathcal{P})  &=& \Psi(\widetilde{\mathcal{P}}) + \frac{d  \Psi(\mathcal{P}_t)}{d  t}\Big|_{t=1} (0-1) + R(\mathcal{P},\widetilde{\mathcal{P}}),
\end{eqnarray*}
where $R(\mathcal{P},\widetilde{\mathcal{P}})$ is a remainder term of the expansion. This expansion contains the pathwise derivative evaluated at $t=1$, which may be evaluated using an anologue of the Riesz-representation theorem result in \eqref{Riesz_non_canonical},
\begin{eqnarray}
\frac{d  \Psi(\mathcal{P}_t)}{d  t}\Big|_{t=1} &=& - \mathcal{P} \left\{ \phi(O,\widetilde{\mathcal{P}})\right\} = - E_{\mathcal{P}}\left\{\phi(O,\widetilde{\mathcal{P}})\right\}, \label{Riesz_non_canonical_t1}
\end{eqnarray}
details of which are given in the appendix. It follows that the scaled difference of interest, equation \eqref{scalediff}, can be written as
\begin{eqnarray}\label{taylor}
\sqrt{n}\left\{ \Psi(\widetilde{\mathcal{P}}) - \Psi(\mathcal{P}) \right\}&=& -\sqrt{n}\mathcal{P} \left\{ \phi(O,\widetilde{\mathcal{P}})\right\} -\sqrt{n}R(\mathcal{P},\widetilde{\mathcal{P}})
\end{eqnarray}
where we note that this identity is guaranteed to hold, so long as we impose no restrictions on the remainder term, which we will consider later. Now letting $\widetilde{\mathcal{P}}$ equal $\hat{\mathcal{P}}_n$, we thus see that
\begin{align*}
\sqrt{n}\left\{ \Psi(\hat{\mathcal{P}}_n) - \Psi(\mathcal{P}) \right\}&= -\sqrt{n}\mathcal{P} \left\{ \phi(O,\hat{\mathcal{P}}_n)\right\} -\sqrt{n}R(\mathcal{P},\hat{\mathcal{P}}_n)\\
&\approx -\frac{1}{\sqrt{n}}\sum_{i=1}^n\phi(O_i,\hat{\mathcal{P}}_n) -\sqrt{n}R(\mathcal{P},\hat{\mathcal{P}}_n),
\end{align*}
Here, the term
\begin{equation}
-\frac{1}{\sqrt{n}}\sum_{i=1}^n\phi(O_i,\hat{\mathcal{P}}_n) \label{drift}
\end{equation}
does not converge to zero (indeed, it would not even converge to zero if $\mathcal{P}$ were used in lieu of $\hat{\mathcal{P}}_n$) and may sometimes even diverge.
This tends not to cause asymptotic bias in $\Psi(\hat{\mathcal{P}}_n)$ (because the calculation of bias requires further scaling by $1/\sqrt{n}$ and, moreover, $\phi(O,\mathcal{P})$ has mean zero and $\mathcal{P}_n$ is assumed to converge to $\mathcal{P}$).
However, it biases the scaled difference $\sqrt{n}\left\{ \Psi(\hat{\mathcal{P}}_n) - \Psi(\mathcal{P}) \right\}$, thereby invalidating na\"{\i}ve confidence intervals and tests.
Understanding the behaviour of \eqref{drift} is difficult as a result of the non-standard behaviour of statistical/machine learning-based estimators affecting the large sample behaviour of $\hat{\mathcal{P}}_n$, which in turn propagates into the behaviour of $\Psi(\hat{\mathcal{P}}_n)$.
Let us reconsider Example 5 (the potential outcome mean), for instance, where
\begin{align*}
\Psi(\mathcal{P})&=E_{\mathcal{P}}\left\{E_{\mathcal{P}}(Y|X=1,Z)\right\}.
\end{align*}
A plug-in estimator is readily obtained as
\begin{align*}
\Psi(\hat{\mathcal{P}}_n)&=\frac{1}{n}\sum_{i=1}^n m_1(Z_i,\hat{\mathcal{P}}_n).
\end{align*}
Here, $m_1(z,\hat{\mathcal{P}}_n)$ denotes a data-adaptive estimator of $m_1(z,\mathcal{P})=E_{\mathcal{P}}(Y|X=1,Z=z)$, e.g. obtained using parametric regression models with variable selection, or via machine learning algorithms.
The plug-in bias term (\ref{drift}) then equals\footnote{Note that we have evaluated the plug-in bias term at the true propensity score because the considered plug-in estimator does not rely on an estimated propensity score. One may alternatively evaluate the plug-in bias term at an estimated propensity score, which will then only affect the remainder term.}
\begin{align*}
&-\frac{1}{\sqrt{n}}\sum_{i=1}^n \frac{\mathbbm{1}_{1}(X_i)}{\pi(Z_i,\mathcal{P})} \left\{Y_i - m_1(Z_i,\hat{\mathcal{P}}_n)\right\} + m_1(Z_i,\hat{\mathcal{P}}_n) - \Psi(\hat{\mathcal{P}}_n) =-\frac{1}{\sqrt{n}}\sum_{i=1}^n \frac{\mathbbm{1}_{1}(X_i)}{\pi(Z_i,\mathcal{P})} \left\{Y_i - m_1(Z_i,\hat{\mathcal{P}}_n)\right\}\\
&=-\frac{1}{\sqrt{n}}\sum_{i=1}^n \frac{\mathbbm{1}_{1}(X_i)}{\pi(Z_i,\mathcal{P})} \left\{Y_i - m_1(Z_i,\mathcal{P})\right\}+\frac{1}{\sqrt{n}}\sum_{i=1}^n \frac{\mathbbm{1}_{1}(X_i)}{\pi(Z_i,\mathcal{P})} \left\{m_1(Z_i,\mathcal{P})-m_1(Z_i,\hat{\mathcal{P}}_n)\right\},
\end{align*}
where the second term will often follow a non-standard distribution. For instance, when $m_1(z,\hat{\mathcal{P}}_n)$ is obtained using parametric regression models with variable selection, it will often follow a mixture distribution for each $z$ as a result of variation in the selected model across repeated samples.

The extent to which the plug-in bias term (\ref{drift}) causes bias is thus generally poorly understood as it inherits the behaviour of $\hat{\mathcal{P}}_n$, which is complex when data-adaptive methods are used.
Rather than attempting to understand its asymptotic behaviour, a much simpler remedy is therefore to adjust the plug-in estimator in such a way that the this bias is zero. One easy way to do this is by defining a new estimator
\[\Psi(\hat{\mathcal{P}}_n)+\frac{1}{n}\sum_{i=1}^n\phi(O_i,\hat{\mathcal{P}}_n)\]
obtained by subtracting an estimate of the plug-in bias from the plug-in estimator. Then, the scaled difference between this so-called one-step estimator and $\Psi(\mathcal{P})$ is governed by $-\sqrt{n}R(\mathcal{P},\hat{\mathcal{P}}_n)$, which will generally be much smaller.

We will see later that there are other ways of modifying the plug-in estimator so that the resulting estimator has zero plug-in bias.

\subsection{The von Mises expansion}

In the previous section, we have built some intuition into plug-in bias and how it can be removed. In order to understand the behaviour of the scaled difference between the one-step estimator and $\Psi(\mathcal{P})$, a more careful derivation is needed. In particular, because $\mathcal{P}$ is unknown, we substituted it by the empirical distribution function $\mathcal{P}_n$ of the observed data, but did not express the error this is adding to the results. Let us therefore take a step back to identity (\ref{taylor}).
By adding and subtracting $\sqrt{n}(\mathcal{P}_n - \mathcal{P}) \left\{\phi(O,\mathcal{P})\right\}$ and $\sqrt{n}\mathcal{P}_n \left\{\phi(O,\widetilde{\mathcal{P}})\right\}$ to the righthand side, we obtain
\begin{eqnarray*}
\sqrt{n}\left\{\Psi(\widetilde{\mathcal{P}}) - \Psi(\mathcal{P})\right\}
&=&\sqrt{n} (\mathcal{P}_n-\mathcal{P})\left\{\phi(O,\mathcal{P})\right\}-\sqrt{n}\mathcal{P}_n\left\{\phi(O,\widetilde{\mathcal{P}})\right\}\\&&+\sqrt{n} (\mathcal{P}_n-\mathcal{P})\left\{\phi(O,\widetilde{\mathcal{P}})-\phi(O,\mathcal{P})\right\}-\sqrt{n}R(\mathcal{P},\widetilde{\mathcal{P}}).
\end{eqnarray*}
Setting $\widetilde{\mathcal{P}}$ to $\hat{\mathcal{P}}_n$ one can rewrite the plug-in bias in form of the so-called von Mises expansion:
\begin{eqnarray}
\nonumber \sqrt{n}\left\{\Psi(\hat{\mathcal{P}}_n) - \Psi(\mathcal{P})\right\} &=& -\sqrt{n}\mathcal{P} \left\{ \phi(O,\hat{\mathcal{P}}_n)\right\} -\sqrt{n}R(\mathcal{P},\hat{\mathcal{P}}_n)  \\
&=&\frac{1}{\sqrt{n}}\sum_{i=1}^n \phi(O_i,\mathcal{P}) -\frac{1}{\sqrt{n}}\sum_{i=1}^n\phi(O_i,\hat{\mathcal{P}}_n) \nonumber \\&&+\sqrt{n} (\mathcal{P}_n-\mathcal{P})\left\{\phi(O,\hat{\mathcal{P}}_n)-\phi(O,\mathcal{P})\right\}-\sqrt{n}R(\mathcal{P},\hat{\mathcal{P}}_n) \label{vonMises}.
\end{eqnarray}
Here, the first term converges to a normal, mean zero variate by the central limit theorem and the unbiasedness of the canonical gradient.
The empirical process term (i.e., the third term in \eqref{vonMises}) and the remainder term $\sqrt{n}R(\mathcal{P},\hat{\mathcal{P}}_n)$ can often be shown to converge to zero under conditions that we will come back to.

Since the asymptotic behaviour of $\hat{\mathcal{P}}_n$, and therefore also of the second term, is often poorly understood, popular approaches are designed to remove this drift term from the expansion. This can be done in multiple possible ways.

\paragraph{One-step estimator.}
The first is to rewrite the above expansion as
\begin{eqnarray*}
\sqrt{n}\left\{\Psi(\hat{\mathcal{P}}_n)+\frac{1}{n}\sum_{i=1}^n\phi(O_i,\hat{\mathcal{P}}_n) - \Psi(\mathcal{P})\right\}
&=&\frac{1}{\sqrt{n}}\sum_{i=1}^n \phi(O_i,\mathcal{P}) \\&&+\sqrt{n} (\mathcal{P}_n-\mathcal{P})\left\{\phi(O,\hat{\mathcal{P}}_n)-\phi(O,\mathcal{P})\right\}-\sqrt{n}R(\mathcal{P},\hat{\mathcal{P}}_n),
\end{eqnarray*}
and thus to calculate the estimator of $\Psi(\mathcal{P})$ as the \emph{one-step estimator}
\[\Psi(\hat{\mathcal{P}}_n)+\frac{1}{n}\sum_{i=1}^n\phi(O_i,\hat{\mathcal{P}}_n).\]
In Example 3, this delivers
\[\int f^2(y,\hat{\mathcal{P}}_n)dy+\frac{2}{n}\sum_{i=1}^n\left\{f(y_i,\hat{\mathcal{P}}_n)-\int f^2(y,\hat{\mathcal{P}}_n)dy\right\}=\left\{\frac{2}{n}\sum_{i=1}^n f(y_i,\hat{\mathcal{P}}_n)\right\}-\int f^2(y,\hat{\mathcal{P}}_n)dy.\]
where $f(y,\hat{\mathcal{P}}_n)$ is a density estimator.
For Example 5 we consider two different cases. When the propensity score $\pi(Z_i,\mathcal{P})$ is known, for instance in randomized experiments, then one obtains the estimator,
\[\Psi(\hat{\mathcal{P}}_n)+\frac{1}{n}\sum_{i=1}^n \frac{\mathbbm{1}_{1}(X_i)}{\pi(Z_i,\mathcal{P})} \left\{Y_i - m_1(Z_i,\hat{\mathcal{P}}_n)\right\} + m_1(Z_i,\hat{\mathcal{P}}_n) - \Psi(\hat{\mathcal{P}}_n),\]
When the propensity score is unknown, as is the case for observational data, it must also be estimated (e.g. using a data-adaptive estimator $\pi(Z_i,\hat{\mathcal{P}}_n)$) and the one-step estimator recovers the augmented IPW estimator,
\[\frac{1}{n}\sum_{i=1}^n \frac{\mathbbm{1}_{1}(X_i)}{\pi(Z_i,\hat{\mathcal{P}}_n)} \left\{Y_i - m_1(Z_i,\hat{\mathcal{P}}_n)\right\} + m_1(Z_i,\hat{\mathcal{P}}_n),\]
This propensity score estimation has consequences for the remainder term; see below.

\paragraph{Estimating equation estimators.}
The second is to force the drift term to be zero by using it as an estimating equation; that is, to calculate an estimator for $\Psi(\mathcal{P})$ as the solution to an estimating equation given by this drift term:
\begin{align}
0&=\frac{1}{n}\sum_{i=1}^n\phi(O_i,\hat{\mathcal{P}}_n). \label{eq_drift_term}
\end{align}
This is easy in the above examples, where the efficient influence function is linear in $\Psi(\mathcal{P})$. In Example 3, solving the identity
\[0=\frac{2}{n}\sum_{i=1}^n\left\{f(y_i,\hat{\mathcal{P}}_n)-\Psi(\hat{\mathcal{P}}_n)\right\}\]
delivers a different estimator than the one-step estimator, namely
\[\Psi(\hat{\mathcal{P}}_n)=\frac{1}{n}\sum_{i=1}^nf(y_i,\hat{\mathcal{P}}_n),\]
with the advantage that it is guaranteed to be non-negative.
In Example 5, solving the identity
\[0=\frac{1}{n}\sum_{i=1}^n \frac{\mathbbm{1}_{1}(X_i)}{\pi(Z_i,\hat{\mathcal{P}}_n)} \left\{Y_i - m_1(Z_i,\hat{\mathcal{P}}_n)\right\} + m_1(Z_i,\hat{\mathcal{P}}_n) - \Psi(\hat{\mathcal{P}}_n)\]
for $\Psi(\hat{\mathcal{P}}_n)$ delivers the same estimator as the one-step estimator.

\paragraph{Targeted learning.}
The third works instead by tuning the initial estimator $\hat{\mathcal{P}}_n$ such that it forces \eqref{eq_drift_term} to hold, which is the focus of targeted learning approaches \citep{van_der_laan_targeted_2006,van_der_laan_targeted_2011}. For instance, tuning the estimator $\hat{\mathcal{P}}_n$ in Example 5 to a
retargeted estimator $\hat{\mathcal{P}}_n^*$ that satisfies
\[0=\frac{1}{n}\sum_{i=1}^n \frac{\mathbbm{1}_{1}(X_i)}{\pi(Z_i,\hat{\mathcal{P}}^*_n)} \left\{Y_i - m_1(Z_i,\hat{\mathcal{P}}^*_n)\right\},\]
ensures that the one-step estimator reduces to the simple plug-in estimator
\begin{align*}
&\frac{1}{n}\sum_{i=1}^n m_1(Z_i,\hat{\mathcal{P}}^*_n),
\end{align*}
which then has standard asymptotic behaviour. This tuning can be achieved in many ways; for instance, one may leave the propensity score model unchanged by defining $\pi(Z_i,\hat{\mathcal{P}}^*_n)=\pi(Z_i,\hat{\mathcal{P}}_n)$ and tune the outcome model by defining,
\[m_1(Z_i,\hat{\mathcal{P}}^*_n)=m_1(Z_i,\hat{\mathcal{P}}_n)+\hat{\epsilon} \frac{1}{\pi(Z_i,\hat{\mathcal{P}}_n)},\]
where $\hat{\epsilon} $ is chosen to set the plug-in bias to zero, i.e., it is the solution to
\[0=\frac{1}{n}\sum_{i=1}^n \frac{\mathbbm{1}_{1}(X_i)}{\pi(Z_i,\hat{\mathcal{P}}^*_n)} \left\{Y_i - m_1(Z_i,\hat{\mathcal{P}}_n)-\hat{\epsilon} \frac{1}{\pi(Z_i,\hat{\mathcal{P}}_n)}\right\}.\]
Retargeting an initial density estimator in Example 3 is less straightforward because of the difficulty of ensuring that the retargeted density continues to be a proper density.

Under sufficient conditions that ensure the empirical process and remainder terms to converge to zero, it follows from the above expansion that all 3 above approaches deliver an estimator $\Psi(\hat{\mathcal{P}}^*_n)$ whose asymptotic distribution obeys
\begin{equation}\label{adist}
\sqrt{n}\left\{\Psi(\hat{\mathcal{P}}^*_n) - \Psi(\mathcal{P})\right\} \overset{d}{\to} \n{0}{\mathcal{P} \left\{\phi(Y,\mathcal{P})^2\right\}}.
\end{equation}
This is a powerful result, which means that the asymptotic efficiency bound for a nonparametric estimand can be derived as the expected square of the efficient influence function. Heuristically, this bound is a nonparametric analogue of the Cramer-Rao lower bound, and estimators of the type in \eqref{adist} are said to be asymptotically efficient, in the sense that they are asymptotically equivalent to the estimator obtained by solving an estimating equation with known rather than estimated influence function:
\[0=\frac{1}{n}\sum_{i=1}^n\phi(O_i,\mathcal{P}).\]
It thus tells us that the influence function behaves like the score function in parametric estimation \citep{Wasserman2006nonparam}. It motivates why the variance of $\Psi(\hat{\mathcal{P}}_n)$ can be estimated as 1 over $n$ times the sample variance of the efficient influence function (evaluated at $\hat{\mathcal{P}}_n$), without needing to account for the uncertainty in $\hat{\mathcal{P}}_n$.
Identity (\ref{adist}) also motivates why the definition of pathwise differentiability includes the requirement of an efficient influence function with finite variance. Pathwise differentiability of an estimand is therefore tantamount to the existence of (regular) root-$n$ consistent estimators of that estimand.

\subsection{Controlling the empirical process term}

The asymptotic behaviour of the empirical process term
\[\sqrt{n} (\mathcal{P}_n-\mathcal{P})\left\{\phi(O,\hat{\mathcal{P}}_n)-\phi(O,\mathcal{P})\right\}\]
is generally difficult to understand when data-adaptive statistical methods are used.
However, it becomes much simpler to understand when the estimator $\hat{\mathcal{P}}_n$ is trained on an independent dataset, as one can then reason conditional on that estimator.
Reasoning as such, a direct application of Chebyshev's inequality shows that the empirical process term converges to zero in probability when the conditional variance of $\phi(O,\hat{\mathcal{P}}_n)-\phi(O,\mathcal{P})$, i.e.,
\[\mathcal{P}\left[\left\{\phi(O,\hat{\mathcal{P}}_n)-\phi(O,\mathcal{P})\right\}^2\right]\]
given $\hat{\mathcal{P}}_n$, converges to zero in probability. The latter can often be shown to hold when the estimator $\hat{\mathcal{P}}_n$ converges to $\mathcal{P}$ in probability (or even weaker conditions that certain functionals of $\hat{\mathcal{P}}_n$ converge to the corresponding functionals of $\mathcal{P}$ in probability) and certain positivity conditions hold (see for instance \cite{chernozhukov2017double,Vansteelandt2020} for detailed examples). The use of an independent sample in this way is important for shrinking the empirical process term, but contrary to what popular wisdom sometimes seems to suggest, does not eliminate the leading plug-in bias terms on which we have focused.

Because one rarely has independent data available to train $\hat{\mathcal{P}}_n$, \cite{van_der_laan_cross-validated_2011} and \cite{chernozhukov_double/debiased_2018} recommend a cross-fitting procedure, whereby the data is split into $K$ folds.
For each individual $i$ from fold $k=1,...,K$, the efficient influence function for that individual is then evaluated in an estimator $\hat{\mathcal{P}}_n$ trained on the data for all individuals, except those in the $k$-th fold. This usually results in a better asymptotic approximation, as reflected by more accurate standard error estimators obtained as 1 over root-$n$ times the sample standard deviation of those influence functions. However, it may induce some finite-sample bias in the estimator as a result of the data-adaptive estimator $\hat{\mathcal{P}}_n$ being trained on a smaller sample of data.

\subsection{Controlling the remainder term}

To understand the remainder term $\sqrt{n}R(\mathcal{P},\hat{\mathcal{P}}_n)$, we return to the von Mises expansion (\ref{vonMises}), from which it is seen to equal
\begin{eqnarray*}
\sqrt{n}R(\mathcal{P},\hat{\mathcal{P}}_n)
&=&-\sqrt{n} \mathcal{P}\left\{\phi(O,\hat{\mathcal{P}}_n)\right\}-\sqrt{n}\left\{\Psi(\hat{\mathcal{P}}_n) - \Psi(\mathcal{P})\right\}.
\end{eqnarray*}
In Example 5, this is
\begin{eqnarray*}
\sqrt{n}R(\mathcal{P},\hat{\mathcal{P}}_n)
&=&-\sqrt{n}E_{\mathcal{P}}\left[\frac{\mathbbm{1}_{1}(X)}{\pi(Z,\mathcal{P})} \left\{Y - m_1(Z,\hat{\mathcal{P}}_n)\right\} + m_1(Z,\hat{\mathcal{P}}_n) - \Psi(\hat{\mathcal{P}}_n)\right] -\sqrt{n}\left\{\Psi(\hat{\mathcal{P}}_n) - \Psi(\mathcal{P})\right\}\\
&=&-\sqrt{n}E_{\mathcal{P}}\left[\frac{\mathbbm{1}_{1}(X)}{\pi(Z,\mathcal{P})} \left\{Y - m_1(Z,\hat{\mathcal{P}}_n)\right\} + m_1(Z,\hat{\mathcal{P}}_n) - \Psi(\mathcal{P})\right] \\
&=&-\sqrt{n}E_{\mathcal{P}}\left[\left\{\frac{\pi(Z,\mathcal{P})}{\pi(Z,\mathcal{P})} -1\right\}\left\{m_1(Z,\mathcal{P}) - m_1(Z,\hat{\mathcal{P}}_n)\right\}\right]=0.
\end{eqnarray*}
Hence when the propensity score $\pi(Z,\mathcal{P})$ is known, the remainder term is zero. Other estimands are also known to have a zero remainder, such as the average density (see Example 3).

When substituting $\pi(Z,\hat{\mathcal{P}}_n)$ for $\pi(Z,\mathcal{P})$,
by the Cauchy-Schwarz inequality, the remainder can be upper bounded by
\[\sqrt{n}E_{\mathcal{P}}\left[\left\{\frac{\pi(Z,\mathcal{P})}{\pi(Z,\hat{\mathcal{P}}_n)} -1\right\}^2\right]^{1/2} E_{\mathcal{P}}\left[\left\{m_1(Z,\mathcal{P}) - m_1(Z,\hat{\mathcal{P}}_n)\right\}^2\right]^{1/2}.\]
This converges to zero in probability when $\pi(Z,\hat{\mathcal{P}}_n)$ and $m_1(Z,\hat{\mathcal{P}}_n)$ converge to $\pi(Z,\mathcal{P})$ and $m_1(Z,\mathcal{P})$, respectively, at faster than $n$ to the quarter rate (and $\pi(Z,\hat{\mathcal{P}}_n)$ is bounded away from zero), which is a typical requirement in the non/semiparametric literature. In this specific example, the remainder also shrinks to zero under more general conditions; $m_1(Z,\hat{\mathcal{P}}_n)$ can be allowed to converge at a slow rate, so long as $\pi(Z,\hat{\mathcal{P}}_n)$ is fast converging. This additional flexibility is sometimes known as `rate double-robustness', and does not apply to remainder terms in general, although it does apply for many common estimands in causal inference/missing data problems \citep{rotnitzky2021characterization}. To obtain fast rates of convergence with flexible methods, we typically rely on strong smoothness/sparsity assumptions (e.g. when $Z$ is high dimensional, $\pi(Z,\mathcal{P})$ and/or $m_1(Z,\mathcal{P})$ should depend on a small number of the covariates), in addition to well-chosen tuning parameters for the learners.

We refer the reader to \cite{Fisher2020} for a rigorous treatment of the remainder terms of the von Mises expansion, which are usually analysed on a case-by-case basis (see for instance \cite{chernozhukov2017double,Vansteelandt2020} for detailed examples).

\section{Examples}\label{derivations}

In this section, we illustrate the calculation of the canonical gradient for the expected conditional covariance and the average derivative effect, deriving the one-step estimators in both cases. Further examples are provided in Appendix B, which also contains results that readers may find helpful for reference.

\subsection{General results}
\label{standard}

For notational convenience we define an operator, $\partial_t$, applied to an arbitrary function, $g(t)$, as
\begin{align*}
\partial_t g(t) = \frac{d  g(t)}{d  t}\Big|_{t=0}.
\end{align*}
For instance, let $f_t(y,x)$ denote a parametric submodel which disturbs the density $f(y,x)$ of $(Y,X)$ at $(y,x)$ in the direction of a point mass at $(\tilde{y},\tilde{x})$. Then from
\begin{align*}
f_t(y|x) &= \frac{f_t(y,x)}{f_t(x)}
\end{align*}
and using the chain rule and the quotient rule for derivatives, we obtain
\begin{align*}
\partial_t f_t(y|x) &= \partial_t \left\{\frac{f_t(y,x)}{f_t(x)}\right\} \\
&=  \frac{\partial_t f_t(y,x) f(x) - f(y,x)\partial_t f_t(x)}{f^2(x)}\\
&= \frac{1}{f(x)} \left[ \mathbbm{1}_{\tilde{y},\tilde{x}}(y,x)-f(y,x) - \frac{f(y,x)}{f(x)} \left\{\mathbbm{1}_{\tilde{x}}(x)-f(x)\right\}  \right] \\
&= \frac{\mathbbm{1}_{\tilde{x}}(x)}{f(x)} \left\{ \mathbbm{1}_{\tilde{y}}(y)-f(y|x) \right\}.
\end{align*}
Similarly to \eqref{expectation}, this expression may be used to derive the following identity for the conditional expectation of an arbitrary function $g(o,\mathcal{P})$, where $o=(y,x')'$:
\begin{align}
\partial_t E_{\mathcal{P}_t}\left\{g(O,\mathcal{P}_t)|X=x\right\} &= \partial_t \int g(o,\mathcal{P}_t) f_t(y|x) dy \nonumber \\
&= \frac{\mathbbm{1}_{\tilde{x}}(x)}{f(x)} \left[g(\tilde{o},\mathcal{P}) - E_{\mathcal{P}}\left\{g(O,\mathcal{P})|X=x\right\} \right] + E_{\mathcal{P}}\left\{\partial_t g(O,\mathcal{P}_t)|X=x\right\}. \label{conditional_expectation}
\end{align}
Such generic expressions are helpful to relate to, as they can be used to speed up derivations. For instance, for the potential outcome mean,
defining $m_1(Z,\mathcal{P})=E_{\mathcal{P}}(Y|X=1,Z)$,
it readily follows from \eqref{expectation} that
\begin{align*}
\partial_t E_{\mathcal{P}_t}\left\{m_1(Z,\mathcal{P}_t)\right\} &= m_1(\tilde{z},\mathcal{P}) - E_{\mathcal{P}}\left\{m_1(Z,\mathcal{P})\right\} + E_{\mathcal{P}}\left\{\partial_t m_1(Z,\mathcal{P}_t)\right\},
\end{align*}
and by \eqref{conditional_expectation}, that
\begin{align*}
\partial_t m_1(Z,\mathcal{P}_t)\ &=  \frac{\mathbbm{1}_{\tilde{x},\tilde{z}}(1,z)}{f(1,z)} \left\{\tilde{y}-m_1(z,\mathcal{P})\right\}+0.
\end{align*}
Averaging over the distribution of $Z$ then delivers
\begin{align*}
E_{\mathcal{P}}\left\{\partial_t m_1(Z,\mathcal{P}_t)\right\} &=  \frac{\mathbbm{1}_{\tilde{x}}(1)}{f(1|z)} \left\{\tilde{y}-m_1(z,\mathcal{P})\right\}.
\end{align*}
Hence, we recover the same result as before.

\paragraph{Example 7 (expected conditional covariance).}
Consider the expected conditional covariance,
\begin{align*}
\Psi(\mathcal{P}) = E_{\mathcal{P}}\left[\left\{Y-E_{\mathcal{P}}(Y|Z)\right\}\left\{X-E_{\mathcal{P}}(X|Z)\right\}\right]
\end{align*}
which appears in hypothesis testing \citep{Shah2018} and in parameter estimation in generalized linear models \citep{Vansteelandt2020}.
Define
\begin{align*}
\mathrm{cov}_t(Y,X|Z) = E_{\mathcal{P}_t}\left[\left\{Y-E_{\mathcal{P}_t}(Y|Z)\right\}\left\{X-E_{\mathcal{P}_t}(X|Z)\right\}|Z\right].
\end{align*}
Upon noting that $\Psi(\mathcal{P}) = E_{\mathcal{P}}\left\{\mathrm{cov}(Y,X|Z)\right\} $ is of the form in \eqref{expectation}, we find that
\begin{align*}
\partial_t\Psi(\mathcal{P}_t) = \mathrm{cov}(Y,X|\tilde{z}) - \Psi(\mathcal{P}) + E_{\mathcal{P}}\left\{\partial_t \mathrm{cov}_t(Y,X|Z) \right\}.
\end{align*}
The complication is clearly in the final term, which is of the form in \eqref{conditional_expectation}, hence,
\begin{align*}
\partial_t \mathrm{cov}_t(Y,X|z) = \frac{\mathbbm{1}_{\tilde{z}}(z)}{f(z)} \left[\left\{\tilde{y}-E_{\mathcal{P}}(Y|\tilde{z})\right\}\left\{\tilde{x}-E_{\mathcal{P}}(X|\tilde{z})\right\} - \mathrm{cov}(Y,X|z)\right] \\
+ E[\partial_t\left\{Y-E_{\mathcal{P}_t}(Y|Z)\right\}\left\{X-E_{\mathcal{P}_t}(X|Z)\right\}|Z=z].
\end{align*}
Similarly to the covariance example previously, the final term above turns out to be zero. It follows, therefore, that the canonical gradient is
\begin{align*}
\phi(O,\mathcal{P}) = \partial_t\Psi(\mathcal{P}_t) = \left\{Y-E_{\mathcal{P}}(Y|Z)\right\}\left\{X-E_{\mathcal{P}}(X|Z)\right\}- \Psi(\mathcal{P}),
\end{align*}
and since this has finite variance,  the expected conditional covariance is pathwise differentiable.

Constructing a one-step estimator or estimating equations estimator based on the canonical gradient of the expected conditional covariance is relatively straightforward and in fact both methods will provide the same result in this example. The one-step estimator  takes an original plug-in estimator $\Psi(\hat{\mathcal{P}}_n)$ and adds a correction term
\begin{align*}
\Psi(\hat{\mathcal{P}}_n) + \frac{1}{n} \sum_{i=1}^n \phi(O_i,\hat{\mathcal{P}}_n) &= \frac{1}{n}\sum_{i=1}^n \left\{Y_i-\hat{m}(Z_i)\right\}\left\{X_i-\hat{\pi}(Z_i)\right\}.
\end{align*}
where $\hat{m}(z) = E_{\hat{\mathcal{P}}_n}(Y|z)$ and $\hat{\pi}(z) = E_{\hat{\mathcal{P}}_n}(X|z)$. Estimation by this strategy therefore requires additional modelling to obtain the functions $\hat{m}(z)$ and $\hat{\pi}(z)$. $\Box$

\paragraph{Example 8 (average derivative effect).}

This example concerns the average derivative effect estimand of \cite{Hardle1989} with canonical gradient given by \cite{Newey1993}.
We let $m(x,z,\mathcal{P}) = E_{\mathcal{P}}(Y|X=x,Z=z)$ be a conditional response surface which is assumed to be differentiable w.r.t. $x$, with derivative $m^\prime(x,z,\mathcal{P})$, and we also introduce a known weight function, $w(x,z)$. The average derivative effect estimand is written
\begin{equation*}
\Psi(\mathcal{P}) = E_{\mathcal{P}}\left\{w(X,Z)m^\prime(X,Z,\mathcal{P})\right\}.
\end{equation*}
\cite{Powell1989} showed that, for a differentiable function $g(x,z)$, with derivative w.r.t. $x$, $g^\prime(x,z)$,
\begin{equation*}
E_{\mathcal{P}}\left\{w(X,Z)g^\prime(X,Z)\right\} = E_{\mathcal{P}}\left\{l(X,Z,\mathcal{P})g(X,Z)\right\}
\end{equation*}
under regularity conditions, which require that $X$ is a continuous random variable and that $w(x,z)f(x,z)$ is differentiable w.r.t. $x$ and is zero on the boundary of the support of $X$, where $f(x,z)$ used to denote the joint distribution of $(X,Z)$ under $\mathcal{P}$. In the above expression, \[l(x,z,\mathcal{P}) \equiv -w^\prime(x,z) - w(x,z) f^\prime(x,z)/f(x,z),\] and, as before, superscript prime denotes the derivative with respect to $x$. Using \eqref{expectation},
\begin{align*}
\partial_t \Psi(\mathcal{P}_t) = w(\tilde{x},\tilde{z})m^\prime(\tilde{x},\tilde{z}) - \Psi(\mathcal{P}) + E_{\mathcal{P}}\left\{w(X,Z)\partial_t m^\prime(X,Z,\mathcal{P}_t)\right\}.
\end{align*}
For the final term, we rely on Powell's identity:
\begin{align*}
E_{\mathcal{P}}\left\{w(X,Z)\partial_t m^\prime(X,Z,\mathcal{P}_t)\right\} &= \partial_t E_{\mathcal{P}}\left\{w(X,Z) m^\prime(X,Z,\mathcal{P}_t)\right\}\\
&= \partial_t E_{\mathcal{P}}\left\{l(X,Z,\mathcal{P}) m(X,Z,\mathcal{P}_t)\right\}\\
&=  E_{\mathcal{P}}\left[l(X,Z,\mathcal{P}) \frac{\mathbbm{1}_{\tilde{x},\tilde{z}}(X,Z)}{f(X,Z)}\left\{\tilde{y} - m(X,Z,\mathcal{P})\right\}\right] \\
&= l(\tilde{x},\tilde{z},\mathcal{P})\left\{\tilde{y} - m(\tilde{x},\tilde{z},\mathcal{P})\right\}.
\end{align*}
Since this has finite variance, the average derivative effect is pathwise differentiable with canonical gradient
\begin{align*}
\phi(O,\mathcal{P}) = \partial_t \Psi(\mathcal{P}_t) = l(X,Z,\mathcal{P})\left\{Y - m(X,Z,\mathcal{P})\right\} + w(X,Z)m^\prime(X,Z,\mathcal{P}) - \Psi(\mathcal{P}).
\end{align*}
Using this efficient influence function, an efficient estimator may be easily derived following the one-step or estimating equation strategy. In this case both will result in the same estimator. Setting the sample average of $\phi(O_i,\hat{\mathcal{P}}_n)$ to zero results in the estimator
\begin{align*}
\Psi(\hat{\mathcal{P}}_n) &= \frac{1}{n}\sum_{i=1}^n l(X_i,Z_i,\hat{\mathcal{P}}_n)\left\{Y_i - m(X_i,Z_i,\hat{\mathcal{P}}_n)\right\} + w(X_i,Z_i)m^\prime(X_i,Z_i,\hat{\mathcal{P}}_n)
\end{align*}
This estimator therefore requires modelling the functions $m(x,z,\mathcal{P}),m^\prime(x,z,\mathcal{P}) $ and $l(x,z,\mathcal{P})$.

We include some extra examples in the Appendix.
\section{Implementation}

We begin by summarising the steps that need to be followed to go from scientific question to (data-adaptive) estimation described in the previous sections.

\paragraph{Step 1:} Defining the estimand of interest.\\
The estimand $\Psi(\mathcal{P})$ is a nonparametrically defined statistical functional which is chosen with reference to the scientific question of interest. The estimand might be motivated for a variety of reasons, such as with reference to causal inference (e.g. example 5), independence testing (e.g. example 7), variable importance (e.g., \cite{Williamson2021}), etc.

\paragraph{Step 2:} Calculating its efficient influence function (under the nonparametric model). There are several ways  to do this.

\begin{enumerate}
    \item  \textit{Point-mass contamination.}
We compute the G\^{a}teaux/ pathwise derivative of $\Psi(\mathcal{P})$ at $\mathcal{P}$ in direction of a probability point mass $\widetilde{\mathcal{P}}$. We consider a parametric submodel $\mathcal{P}_t = (1-t) \mathcal{P}+t \widetilde{\mathcal{P}}$ for $t\in[0,1]$, which we use to evaluate the efficient influence function,
$$
\phi(o,\mathcal{P})=\frac{d\Psi\left( \mathcal{P}_t \right)}{d t} \Big|_{t=0}.
$$

\item
The most general method is to work from the definition of pathwise differentiability. Define a rich class of submodels $\mathcal{P}_t$ for $t \in (-\epsilon, \epsilon)$ such that $ \mathcal{P}_{0} = \mathcal{P}$ and $S(O)=\frac{d}{d t} \log f_{t}(o)\vert_{t=0}$ is the score function of $t$, where $f_t(o)$ denotes the density/probability point mass of $O$ under $\mathcal{P}_t$. Next one writes the derivative of $\Psi(\mathcal{P}_t)$, as the integral,
$$\frac{d\Psi\left( \mathcal{P}_t \right)}{d t} \Big|_{t=0}=\int \phi(o, \mathcal{P}) S(o) d \mathcal{P}(o).$$
By the Riesz representation theorem, $\phi(o,\mathcal{P})$ is the efficient influence function.

\item For many new estimands, we can  manipulate expressions so they can be expressed as simpler components (by applying the chain and product rules) and use known influence functions as building blocks.
\end{enumerate}

\paragraph{Step 3:} Obtaining an estimator based on the efficient influence function (such as one-step or TMLE) which admits a first-order representation, provided we use sample splitting or cross-fitting.

We recommend cross-fitting in the following way: First the data is partitioned into $K$ \textit{folds}, i.e. $K$ smaller data sets of (roughly) equal size. Next, for each fold $k$, estimate the nuisance functionals (e.g. via data-adaptive methods) using the rest of all the data, excluding that in fold $k$. Use these nuisance functionals to evaluate the efficient influence function for each observations in fold $k$. After repeating for each fold, one is left with an estimate of the efficient influence function for all observations in the dataset. Finally, use these efficient influence function estimates to evaluate the estimator, and the error in the estimator. See \cite{van_der_laan_cross-validated_2011} and \cite{chernozhukov_double/debiased_2018} for more on cross-fitting.

\subsection{Software}
These steps have already been implemented in a variety of R software packages for several estimands from the causal inference literature, e.g., the   average treatment effect, local average treatment effects (using instrumental variables) and the effect curve for a continuous treatment.

For example, the R packages \texttt{npcausal} and \texttt{DoubleML} \citep{DoubleML2021} both have functions to estimate the ATE, the ATT and LATE (based on instrumental variables), and allow the use of sample splitting and a range of machine learners. These packages also implement other estimands that may be of interest to the biostatistics and econometrics literature.

The targeted learning literature is rich in software package implementations of a variety of causal estimands. A new package implementing previous targeted learning estimators, is  \texttt{tlverse}.
The TMLE estimators implemented as part of \texttt{tlverse} are all fitted using cross-validation by default (for example ‘tmle3’).
TMLE cross-validated estimators of counterfactual means and causal effects are also implemented in the R package \texttt{drtmle}.
The R package \texttt{AIPW} \citep{Zhong2021} implements estimation of the ATE by AIPW (corresponding to the estimating equations or one-step estimator as we have seen) and also a TMLE estimator based on machine learning  algorithms.

Other more advanced causal estimands, such as causal effects for a continuous treatment and optimal treatment rules are implemented in both \texttt{npcausal} or \texttt{tmle3}.

\section{Discussion}

Statistical education still focuses primarily on parametric statistical models, which are assumed to reflect how the data is generated. The inferential theory that is taught does not reflect how data is usually analysed, where models are chosen data-adaptively and different models may fit equally well, especially nowadays, given the increased popularity of machine learning. We therefore believe that many courses would be better focused on translating a scientific question into an nonparametric estimand, and basing inference on its efficient influence function under the nonparametric model.

Courses and textbook treatments on the calculus of influence functions often focus on (semi)parametric models \citep{Tsiatis2006}. The resulting derivations can be challenging, as they require one to respect the restrictions that the model imposes on the observed data distribution. Moreover, they show how one can use these restrictions in order to make efficiency gains. Extracting information from modelling assumptions nevertheless comes at the risk of invalid inference when assumptions are violated. By contrast, our focus is on inference under a nonparametric model. This not only makes the resulting inferences more honest, but can dramatically simplify calculations. Additional efficiency gains are then reserved for special cases when restrictions are known to hold by the study design \citep{zhang2008improving}, or reflect strong pre-existing scientific knowledge \citep{liu2021efficient}.

It is difficult, however, to proceed entirely nonparametrically and avoid regularity conditions all together. Indeed, without assumptions on distribution tails, inference of the mean, Example 1 in the current paper, is impossible \citep{Bahadur1956,Bickel1975}. Likewise, many of our examples rely on working models for statistical functionals, necessitating certain regularity \citep{Robins1997}. For instance, the one-step estimator for Example 5 requires estimating $m_1(Z_i,\hat{\mathcal{P}}_n)$. Whilst flexible data-adaptive/ machine learning estimators can be used, these are better thought of as very highly parametric rather than nonparametric, and make assumptions on the true functional $m_1(z,\mathcal{P})$, e.g. that it is smooth in $z$. The crucial difference, however, is that compared with the parametric modelling approach, estimators based on the nonparametric model do not `extract efficiency' from highly parametric modelling assumptions.

Because of the crucial role that efficient influence functions play, we focused on their derivation. Whilst the formal justification of the von Mises expansion relies on concepts from advanced mathematics, calculating the efficient influence function can often be done using techniques covered in a basic calculus course. We have illustrated this for several causal and non-causal statistical functionals (estimands); the method of derivation described can lead to simpler proofs than those in the original research papers.

Influence functions have applications beyond using them to define estimators with zero plug-in bias.
Influence functions capture the stability of estimators to outliers (in fact this is one of their original purposes), which makes them additionally useful to diagnose outliers (as measurements with large influence function values). Recently,  influence functions have started to be used in the machine learning literature too. For example, \cite{Koh2017} used influence functions for interpretability of black-box models, by characterising the impact a data point has on the black-box’s predictions. \cite{curth2020semiparametric} and \cite{kennedy2020optimal} use influence functions as the outcome in machine learning procedures of conditional (e.g. subgroup-specific) estimands.

We hope that our contribution helps demystify the calculation of influence functions and thus encourages their wider adoption.