#  Some Doubly and Multiply Robust Estimators of Controlled Direct
Effects††thanks: Direct all correspondence to Xiang Zhou, Department of
Sociology, Harvard University, 33 Kirkland Street, Cambridge MA 02138; email:
<EMAIL_ADDRESS>The author thanks Aleksei Opacic for helpful
comments.

Xiang Zhou

(November 18, 2020)

###### Abstract

This letter introduces several doubly, triply, and quadruply robust estimators
of the controlled direct effect. Among them, the triply and quadruply robust
estimators are locally semiparametric efficient, and well suited to the use of
data-adaptive methods for estimating their nuisance functions.

##  1 Introduction

Over the past decade, causal mediation analysis has grown popular in social
and biomedical sciences. A common approach to assessing causal mediation
involves decomposing the total effect of a treatment on an outcome into the
so-called natural direct and indirect effects (NDE and NIE; Robins and
Greenland 1992; Pearl 2001). The NDE and NIE, however, are not
nonparametrically identified in the presence of posttreatment confounders,
i.e., when confounders of the mediator-outcome relationship may be affected by
the treatment itself (Avin et al. 2005; VanderWeele and Vansteelandt 2009). In
such cases, researchers have often focused on estimating the controlled direct
effect (CDE), a quantity that measures the effect of treatment when a mediator
is fixed at a given value for all units (Pearl 2001). Thus a nonzero CDE
implies that the effect of treatment on the outcome does not operate
exclusively through the mediator of interest. Unlike the NDE and NIE, the CDE
is identified provided that all confounders for the treatment-outcome
relationship and for the mediator-outcome relationship are observed, even if
some of the mediator-outcome confounders are affected by the treatment itself.

Estimators of the CDE typically rely on correct specification of (at least)
two nuisance functions about the conditional means/densities of the treatment,
mediator, outcome, or posttreatment confounders. For example, the weighting
estimator proposed by VanderWeele (2009) involves fitting two propensity score
models, one for the treatment and one for the mediator, and the sequential
g-estimator proposed by Vansteelandt (2009) involves fitting two outcome
models, one for the observed outcome given all its antecedent variables and
one for a “demediated” outcome given pretreatment confounders and the
treatment (see also Joffe and Greene 2009). To alleviate bias due to model
misspecification, Goetgeluk et al. (2009) proposed a doubly robust estimator
of the CDE that depends on correct specification of (a) a model for treatment
assignment, and (b) either an outcome model or a propensity score model for
the mediator.

The causal structure underlying identification and estimation of the CDE is
akin to that of estimating treatment effects with longitudinal data in the
presence of time-varying confounders (Robins 1999). For the latter problem,
Bang and Robins (2005) have proposed a doubly robust estimator for the mean of
a potential outcome that depends on correct specification of either (a)
propensity score models for treatment status at all time points or (b) models
for an iteratively imputed outcome at all time points. In a recent paper,
Rotnitzky et al. (2017) point out that the Bang-Robins estimator is actually
“multiply robust” because it is consistent whenever the first $k$ propensity
score models and the last $K-k$ “outcome models” are correctly specified,
where $0\leq k\leq K$, and $K$ is the number of time points (see also Molina
et al. 2017). Moreover, these authors show that the Bang-Robins procedure can
be further improved with a $2^{K}$-robust estimator that requires correct
specification of either the propensity score model or the outcome model at
each time point.111In a separate strand of literature, the term “multiple
robustness” has been used to characterize a class of estimators for the mean
of incomplete data in a cross-sectional setting that are consistent if one of
several models for the propensity score or one of several models for the
outcome is correctly specified (e.g., Han and Wang 2013). Following Molina et
al. (2017) and Rotnitzky et al. (2017), we use “multiple robustness” to
characterize estimators that require modeling multiple parts of the observed
data likelihood and are consistent if one of several sets of the corresponding
models are correctly specified.

Capitalizing on the above work, this letter introduces a set of doubly robust,
triply robust, and quadruply robust estimators of the CDE, which, to the best
of the author’s knowledge, are new to the causal mediation literature. While
some of these estimators ($\hat{\psi}_{am}^{\textup{tr}_{1}}$ and
$\hat{\psi}_{am}^{\textup{qr}}$; see Section 5) are closely related to those
proposed in Bang and Robins (2005) and Rotnitzky et al. (2017) for estimating
time-varying treatment effects, the rest of them have not been discussed
elsewhere. These estimators all involve estimating more than two nuisance
functions; yet, under suitable regularity conditions, they are consistent and
asymptotically normal (CAN) when only two of these nuisance functions are
correctly specified and their estimates are $\sqrt{n}$-consistent. The triply
and quadruply robust estimators are locally efficient, i.e., when all of the
nuisance functions are correctly specified, they attain the semiparametric
efficiency bound in the nonparametric model over observed data. Moreover,
their estimating equations are Neyman orthogonal, encouraging the use of
machine learning methods and cross-fitting to estimate the nuisance functions
(Zheng and van der Laan 2011; Chernozhukov et al. 2018), in which case
estimates of the CDE are semiparametric efficient when estimates of the
nuisance functions, for example, all converge at faster-than-$n^{-1/4}$ rates.

##  2 Notation, Assumptions, and Identification

Let $A$ denote treatment, $M$ the mediator, $Y$ the observed outcome, and
$Y(a,m)$ the potential outcome under treatment status $a$ and mediator value
$m$. We focus on the simple setting where the treatment $A$ and the mediator
$M$ are both discrete with finite support. In addition, we denote by $X$ a
vector of pretreatment variables that may confound the causal effect of
$(A,M)$ on $Y$, and denote by $Z$ a vector of posttreatment variables that may
confound the causal effect of $M$ on $Y$. Note that the posttreatment
confounders $Z$ may themselves be affected by the treatment.

The controlled direct effect (CDE) is defined as the average effect of
switching treatment status from $a^{\prime}$ to $a$ while fixing the mediator
at a given level $m$:

$\textup{CDE}(a,a^{\prime},m)=\mathbb{E}[Y(a,m)-Y(a^{\prime},m)]$

By definition, the CDE is identified when the expected potential outcome
$\mathbb{E}[Y(a,m)]$ is identified for any $a$ and $m$. Thus, we focus on the
latter estimand throughout the paper and denote it as $\psi_{am}$. Since it is
the expected potential outcome when both the treatment and the mediator are
“controlled” at given values, we may refer to it as the controlled response
function (CRF). The CRF can also be used to construct other estimands such as
the controlled mediator effect
$\textup{CME}(a,m,m^{\prime})=\psi_{am}-\psi_{am^{\prime}}$ (Zheng and Zhou
2015) and the treatment-mediator interaction effect
$(\psi_{am}-\psi_{am^{\prime}})-(\psi_{a^{\prime}m}-\psi_{a^{\prime}m^{\prime}})$.

The CRF is identified under the assumptions of consistency, sequential
ignorability, and positivity:

  1. 1.

consistency: for any unit, if $A=a$ and $M=m$, then $Y=Y(a,m)$;

  2. 2.

sequential ignorability: $Y(a,m)\perp\\!\\!\\!\perp A|X,\forall a,m,$ and
$Y(a,m)\perp\\!\\!\\!\perp M|X,A,Z$, $\forall a,m$.

  3. 3.

positivity: $p_{A|X}(a|x)>\epsilon>0$ and $p_{M|X,A,Z}(m|x,a,z)>\epsilon>0$
$\forall a,m$, $x\in\textup{supp}(X)$, and $z\in\textup{supp}(Z|X=x,A=a)$,

where $p(\cdot)$ denotes a probability mass/density function. Under
assumptions 1-3, the CRF (and hence the CDE) can be identified via Robins’s
(1986) g-computation formula:

$\psi_{am}=\iiint ydP(y|x,a,z,m)dP(z|x,a)dP(x),$ (1)

where $P(u|v)$ denotes the cumulative distribution function of $U$ given $V$.

##  3 G-Computation, Imputation, and Weighting

Using the law of iterated expectations, equation (1) can be written in several
different forms, each of which points to a different way of estimating the
CRF:

$\displaystyle\psi_{am}$
$\displaystyle=\iint\mathbb{E}[Y|x,a,z,m]dP(z|x,a)dP(x)\quad\textup{(g-computation)}$
(2)
$\displaystyle=\mathbb{E}_{X}\mathbb{E}_{Z|X,A=a}\mathbb{E}[Y|X,A,Z,M=m]\quad\textup{(pure
imputation)}$ (3)
$\displaystyle=\mathbb{E}\big{[}\frac{\mathbb{I}(A=a)\mathbb{E}[Y|X,A,Z,M=m]}{\Pr[A=a|X]}\big{]}\quad\textup{(imputation-
then-weighting)}$ (4)
$\displaystyle=\text{$\mathbb{E}\big{[}\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)Y}{\Pr[A=a|X]\Pr[M=m|X,A,Z]}\big{]}$\quad{(pure
weighting)}}$ (5)
$\displaystyle=\mathbb{E}_{X}\mathbb{E}\big{[}\frac{\mathbb{I}(M=m)Y}{\Pr[M=m|X,A,Z]}|X,A=a\big{]}\quad\textup{(weighting-
then-imputation)}$ (6)

Equation (2) suggests a procedure akin to Robins’s g-computation algorithm:
(1) fit a parametric model for the conditional distribution of $Z$ given $X$
and $A$; (2) fit a parametric or semiparametric model for the conditional mean
of $Y$ given $X,A,Z$, and $M$; and (3) evaluate the inner integral via Monte
Carlo simulation and the outer integral via the empirical distribution of $X$.
In the particular case where the models for $\mathbb{E}[Z|x,a]$ and
$\mathbb{E}[Y|x,a,z,m]$ are both linear, equation (2) can be evaluated using a
simple “regression-with-residuals” procedure (Zhou and Wodtke 2019). Equation
(3) suggests a “pure imputation” procedure: (1) fit a model for the
conditional mean of $Y$ given $X,A,Z$, and $M$ and obtain predicted values for
all units at $M=m$, $\hat{\mathbb{E}}[Y|X,A,Z,M=m]$; (2) fit a model for the
conditional mean of $\hat{\mathbb{E}}[Y|X,A,Z,M=m]$ given $X$ and $A$ and
obtain its predicted values for all units at $A=a$; (3) average these
predicted values over all units. This procedure is similar to the sequential
g-estimator proposed in Vansteelandt (2009) and Joffe and Greene (2009).
Equation (4) suggests an imputation-then-weighting procedure: (1) fit a model
for the conditional mean of $Y$ given $X,A,Z$, and $M$ and obtain predicted
values for all units at $M=m$, $\hat{\mathbb{E}}[Y|X,A,Z,M=m]$; (2) fit a
propensity score model for treatment status and obtain fitted values
$\widehat{\Pr}[A=a|X]$; (3) compute a weighted average of the predicted
outcomes $\hat{\mathbb{E}}[Y|X,A,Z,M=m]$ with inverse-probability weights
$\mathbb{I}(A=a)/\widehat{\Pr}[A=a|X]$.

Equation (5) suggests a “pure weighting” estimator (VanderWeele 2009): (1) fit
a propensity score model for treatment status and obtain fitted values
$\widehat{\Pr}[A=a|X]$; (2) fit a propensity score model for the mediator and
obtain fitted values $\widehat{\Pr}[M=m|X,A,Z]$; and (3) compute a weighted
average of observed outcomes with inverse-probability weights
$\mathbb{I}(A=a)\mathbb{I}(M=m)/(\widehat{\Pr}[A=a|X]\cdot\widehat{\Pr}[M=m|X,A,Z]\big{)}$.
Finally, equation (6) suggests a weighting-then-imputation procedure: (1) fit
a propensity score model for the mediator and obtain fitted values
$\widehat{\Pr}[M=m|X,A,Z]$; (2) fit a model for the conditional mean of the
inverse-probability-weighted outcome
$\mathbb{I}(M=m)Y/\widehat{\Pr}[M=m|X,A,Z]$ given $X$ and $A$ and obtain
predicted values for all units at $A=a$; (3) average these predicted values
over all units.

All of the above estimators involve estimating two nuisance functions about
the conditional means/distributions of the treatment, mediator, outcome, or
posttreatment confounders. Specifically, the g-computation procedure requires
correctly specified models for $\mathbb{E}[Y|x,a,z,m]$ and $P(z|x,a)$; the
pure imputation estimator requires correctly specified models for
$\mathbb{E}[Y|x,a,z,m]$ and $\mathbb{E}_{Z|x,a}\mathbb{E}[Y|X,A,Z,M=m]$; the
imputation-then-weighting estimator requires correctly specified models for
$\mathbb{E}[Y|x,a,z,m]$ and $\Pr[A=a|x]$; the pure weighting estimator
requires correctly specified models for $\Pr[A=a|x]$ and $\Pr[M=m|x,a,z]$; and
the weighting-then-imputation estimator requires correctly specified models
for $\Pr[M=m|x,a,z]$ and
$\mathbb{E}\big{[}\frac{\mathbb{I}(M=m)Y}{\Pr[M=m|X,A,Z]}|x,a\big{]}$. When
either of the two requisite models is misspecified, the resulting estimator
will be inconsistent. Thus, in empirical applications where the confounders
$X$ and $Z$ have many components, these estimators can be highly prone to
model misspecification bias. In the following section, we introduce four
“doubly robust” estimators, each of which requires correct specification of
one particular nuisance function and either of two other nuisance functions.

##  4 Doubly Robust Estimators

Before proceeding, we introduce the following functions (treating $a$ and $m$
as fixed):

$\displaystyle\mu_{y}(x,z)$ $\displaystyle:=\mathbb{E}[Y|x,a,z,m]$
$\displaystyle\nu_{y}(x)$ $\displaystyle:=\mathbb{E}_{Z|x,a}\mu_{y}(X,Z)$
$\displaystyle\pi_{a}(x)$ $\displaystyle:=\Pr[A=a|x]$
$\displaystyle\pi_{m}(x,z)$ $\displaystyle:=\Pr[M=m|x,a,z],$

Under assumptions 1-3, $\mu_{y}(x,z)=\mathbb{E}[Y(a,m)|x,a,z]$ and
$\nu_{y}(x)=\mathbb{E}[Y(a,m)|x]$. Thus $\mu_{y}(x,z)$ reflects how the
potential outcome $Y(a,m)$ depends on pretreatment confounders $X$ and
posttreatment confounders $Z$ among units with treatment status $a$, and
$\nu_{y}(x)$ reflects how the potential outcome $Y(a,m)$ depends on
pretreatment confounders $X$. Let $\mu_{y}^{\textup{w}}(x,z)$,
$\nu_{y}^{\textup{w}}(x)$, $\pi_{a}^{\textup{w}}(x)$, and
$\pi_{m}^{\textup{w}}(x,z$) denote a set of working models for these nuisance
functions, and let $\hat{\mu}_{y}^{\textup{w}}(x,z)$,
$\hat{\nu}_{y}^{\textup{w}}(x)$, $\hat{\pi}_{a}^{\textup{w}}(x)$, and
$\hat{\pi}_{m}^{\textup{w}}(x,z$) denote their estimates. In particular,
consider three different two-step estimators of $\nu_{y}(x)$:

$\displaystyle\hat{\nu}_{y}^{\textup{w}}(x;\hat{\mu}_{y}^{\textup{w}})$
$\displaystyle=\hat{\mathbb{E}}\big{[}\hat{\mu}_{y}^{\textup{w}}(X,Z)|x,a\big{]}$
(7) $\displaystyle\hat{\nu}_{y}^{\textup{w}}(x;\hat{\pi}_{m}^{\textup{w}})$
$\displaystyle=\hat{\mathbb{E}}\big{[}\frac{\mathbb{I}(M=m)Y}{\hat{\pi}_{m}^{\textup{w}}(X,Z)}|x,a\big{]}$
(8)
$\displaystyle\hat{\nu}_{y}^{\textup{w}}(x;\hat{\mu}_{y}^{\textup{w}},\hat{\pi}_{m}^{\textup{w}})$
$\displaystyle=\hat{\mathbb{E}}\big{[}\hat{\mu}_{y}^{\textup{w}}(X,Z)+\frac{\mathbb{I}(M=m)}{\hat{\pi}_{m}^{\textup{w}}(X,Z)}(Y-\hat{\mu}_{y}^{\textup{w}}(X,Z))|x,a\big{]},$
(9)

where $\hat{\mathbb{E}}[U|x,a]$ denotes estimates of the conditional mean of a
random variable $U$ given $X=x$ and $A=a$. In the above equations, the
notation $\hat{\nu}_{y}^{\textup{w}}(x;\hat{\mu}_{y}^{\textup{w}})$ indicates
that this quantity depends on previous estimates of
$\mu_{y}^{\textup{w}}(x,z)$, and the same applies to
$\hat{\nu}_{y}^{\textup{w}}(x;\hat{\pi}_{m}^{\textup{w}})$ and
$\hat{\nu}_{y}^{\textup{w}}(x;\hat{\mu}_{y}^{\textup{w}},\hat{\pi}_{m}^{\textup{w}})$.
The last expression
$\hat{\nu}_{y}^{\textup{w}}(x;\hat{\mu}_{y}^{\textup{w}},\hat{\pi}_{m}^{\textup{w}})$
can be seen as a doubly robust estimator of $\nu_{y}(x)$: when $\nu_{y}(x)$ is
correctly specified and either $\mu_{y}^{\textup{w}}(x,z)$ or
$\pi_{m}^{\textup{w}}(x,z)$ is correctly specified,
$\hat{\nu}_{y}^{\textup{w}}(x;\hat{\mu}_{y}^{\textup{w}},\hat{\pi}_{m}^{\textup{w}})$
will be consistent for $\nu_{y}(x)$.

Now consider the following estimators of $\psi_{am}$:

$\displaystyle\hat{\psi}_{am}^{\textup{dr}_{1}}$
$\displaystyle=\mathbb{P}_{n}\big{[}\hat{\nu}_{y}^{\textup{w}}(X;\hat{\mu}_{y}^{\textup{w}})+\frac{\mathbb{I}(A=a)}{\hat{\pi}_{a}^{\textup{w}}(X)}\big{(}\hat{\mu}_{y}^{\textup{w}}(X,Z)-\hat{\nu}_{y}^{\textup{w}}(X;\hat{\mu}_{y}^{\textup{w}})\big{)}\big{]}$
$\displaystyle\hat{\psi}_{am}^{\textup{dr}_{2}}$
$\displaystyle=\mathbb{P}_{n}\big{[}\hat{\nu}_{y}^{\textup{w}}(X;\hat{\pi}_{m}^{\textup{w}})+\frac{\mathbb{I}(A=a)}{\hat{\pi}_{a}^{\textup{w}}(X)}\big{(}\frac{\mathbb{I}(M=m)Y}{\hat{\pi}_{m}^{\textup{w}}(X,Z)}-\hat{\nu}_{y}^{\textup{w}}(X;\hat{\pi}_{m}^{\textup{w}})\big{)}\big{]}$
$\displaystyle\hat{\psi}_{am}^{\textup{dr}_{3}}$
$\displaystyle=\mathbb{P}_{n}\Big{[}\frac{\mathbb{I}(A=a)}{\hat{\pi}_{a}^{\textup{w}}(X)}\Big{(}\hat{\mu}_{y}^{\textup{w}}(X,Z)+\frac{\mathbb{I}(M=m)}{\hat{\pi}_{m}^{\textup{w}}(X,Z)}\big{(}Y-\hat{\mu}_{y}^{\textup{w}}(X,Z)\big{)}\Big{)}\Big{]}$
$\displaystyle\hat{\psi}_{am}^{\textup{dr}_{4}}$
$\displaystyle=\mathbb{P}_{n}\Big{[}\hat{\nu}_{y}^{\textup{w}}(X;\hat{\mu}_{y}^{\textup{w}},\hat{\pi}_{m}^{\textup{w}})\Big{]},$

where $\mathbb{P}_{n}[\cdot]=n^{-1}\sum_{i}[\cdot]_{i}.$
$\hat{\psi}_{am}^{\textup{dr}_{1}}$ can be seen as a combination of the pure
imputation estimator and the imputation-then-weighting estimator;
$\hat{\psi}_{am}^{\textup{dr}_{2}}$ a combination of the pure weighting
estimator and the weighting-then-imputation estimator;
$\hat{\psi}_{am}^{\textup{dr}_{3}}$ a combination of the pure weighting
estimator and the imputation-then-weighting estimator; and
$\hat{\psi}_{am}^{\textup{dr}_{4}}$ a combination of the pure imputation
estimator and the weighting-then-imputation estimator. Their double robustness
is given in Proposition 1.

######  Proposition 1.

Under assumptions 1-3 and suitable regularity conditions,

  1. 1.

$\hat{\psi}_{am}^{\textup{dr}_{1}}$ is CAN if $\mu_{y}^{\textup{w}}(x,z)$ is
correctly specified and either $\nu_{y}^{\textup{w}}(x)$ or
$\pi_{a}^{\textup{w}}(x)$ is correctly specified.

  2. 2.

$\hat{\psi}_{am}^{\textup{dr}_{2}}$ is CAN if $\pi_{m}^{\textup{w}}(x,z)$ is
correctly specified and either $\nu_{y}^{\textup{w}}(x)$ or
$\pi_{a}^{\textup{w}}(x)$ is correctly specified.

  3. 3.

$\hat{\psi}_{am}^{\textup{dr}_{3}}$ is CAN if $\pi_{a}^{\textup{w}}(x)$ is
correctly specified and either $\pi_{m}^{\textup{w}}(x,z)$ or
$\mu_{y}^{\textup{w}}(x,z)$ is correctly specified.

  4. 4.

$\hat{\psi}_{am}^{\textup{dr}_{4}}$ is CAN if $\nu_{y}^{\textup{w}}(x)$ is
correctly specified and either $\pi_{m}^{\textup{w}}(x,z)$ or
$\mu_{y}^{\textup{w}}(x,z)$ is correctly specified.

The double robustness of these estimators is due to a similar logic to that of
standard doubly robust estimators for the mean of incomplete data (Scharfstein
et al. 1999; Robins et al. 2007). For example, for
$\hat{\psi}_{am}^{\textup{dr}_{1}}$, when $\mu_{y}^{\textup{w}}(x,z)$ and
$\nu_{y}^{\textup{w}}(x)$ are correctly specified, the second term inside
$\mathbb{P}_{n}[\cdot]$ will have a zero mean (asymptotically), leaving only
$\mathbb{P}_{n}[\hat{\nu}_{y}^{\textup{w}}(X;\hat{\mu}_{y}^{\textup{w}})]$,
i.e., the pure imputation estimator; and when $\mu_{y}^{\textup{w}}(x,z)$ and
$\pi_{a}^{\textup{w}}(x)$ are correctly specified, the terms involving
$\hat{\nu}_{y}^{\textup{w}}(X;\hat{\mu}_{y}^{\textup{w}})$ will have a zero
mean, leaving only
$\mathbb{P}_{n}[\big{(}\mathbb{I}(A=a)/\hat{\pi}_{a}^{\textup{w}}(X)\big{)}\hat{\mu}_{y}^{\textup{w}}(X,Z)]$,
i.e., the imputation-then-weighting estimator.

Among these doubly robust estimators, $\hat{\psi}_{am}^{\textup{dr}_{3}}$ can
be particularly useful in randomized trials where the treatment assignment
mechanism is known. In this case, $\hat{\psi}_{am}^{\textup{dr}_{3}}$ is
consistent as long as either $\pi_{m}^{\textup{w}}(X,Z)$ or
$\mu_{y}^{\textup{w}}(X,Z)$ is correctly specified. In observational studies,
however, none of these nuisance functions is known a priori, and the relative
utility of these estimators will depend on the subject matter knowledge the
investigator might have about the data generating process. For example, if the
investigator has a better understanding of the mediator model than of the
outcome models, $\hat{\psi}_{am}^{\textup{dr}_{2}}$ may be preferred. Yet, in
many applications, little information is available about any part of the data
generating process. In those cases, the multiply robust estimators presented
below will be more useful as they do not hinge on correct specification of any
particular nuisance function. Moreover, as we will see, they are more amenable
to the use of flexible machine learning methods for estimating the nuisance
functions.

##  5 Multiply Robust and Semiparametric Efficient Estimators

Henceforth, let $O=(X,A,Z,M,Y)$ denote the observed data, and
$\mathcal{P}_{\textup{np}}$ a nonparametric model over $O$ wherein all laws
satisfy the positivity assumption described in Section 2. Define the following
of submodels of $\mathcal{P}_{\textup{np}}$:

  * •

$\mathcal{P}_{1}=\\{P\in\mathcal{P}_{\textup{np}}$:
$\mu_{y}^{\textup{w}}(x,z)$ and $\pi_{a}^{\textup{w}}(x)$ are correctly
specified$\\}$

  * •

$\mathcal{P}_{2}=\\{P\in\mathcal{P}_{\textup{np}}$:
$\mu_{y}^{\textup{w}}(x,z)$ and $\nu_{y}^{\textup{w}}(x)$ are correctly
specified$\\}$

  * •

$\mathcal{P}_{3}=\\{P\in\mathcal{P}_{\textup{np}}$:
$\pi_{m}^{\textup{w}}(x,z)$ and $\pi_{a}^{\textup{w}}(x)$ are correctly
specified$\\}$

  * •

$\mathcal{P}_{4}=\\{P\in\mathcal{P}_{\textup{np}}$:
$\pi_{m}^{\textup{w}}(x,z)$ and $\nu_{y}^{\textup{w}}(x)$ are correctly
specified$\\}$

Consider the following estimators of $\psi_{am}$:

$\displaystyle\hat{\psi}_{am}^{\textup{tr}_{1}}$
$\displaystyle=\mathbb{P}_{n}\Big{[}\hat{\nu}_{y}^{\textup{w}}(X;\hat{\mu}_{y}^{\textup{w}})+\frac{\mathbb{I}(A=a)}{\hat{\pi}_{a}^{\textup{w}}(X)}\big{(}\hat{\mu}_{y}^{\textup{w}}(X,Z)-\hat{\nu}_{y}^{\textup{w}}(X;\hat{\mu}_{y}^{\textup{w}})\big{)}+\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)}{\hat{\pi}_{a}^{\textup{w}}(X)\hat{\pi}_{m}^{\textup{w}}(X,Z)}\big{(}Y-\hat{\mu}_{y}^{\textup{w}}(X,Z)\big{)}\Big{]}$
$\displaystyle\hat{\psi}_{am}^{\textup{tr}_{2}}$
$\displaystyle=\mathbb{P}_{n}\Big{[}\hat{\nu}_{y}^{\textup{w}}(X;\hat{\pi}_{m}^{\textup{w}})+\frac{\mathbb{I}(A=a)}{\hat{\pi}_{a}^{\textup{w}}(X)}\big{(}\hat{\mu}_{y}^{\textup{w}}(X,Z)-\hat{\nu}_{y}^{\textup{w}}(X;\hat{\pi}_{m}^{\textup{w}})\big{)}+\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)}{\hat{\pi}_{a}^{\textup{w}}(X)\hat{\pi}_{m}^{\textup{w}}(X,Z)}\big{(}Y-\hat{\mu}_{y}^{\textup{w}}(X,Z)\big{)}\Big{]}$
$\displaystyle\hat{\psi}_{am}^{\textup{qr}}$
$\displaystyle=\mathbb{P}_{n}\Big{[}\hat{\nu}_{y}^{\textup{w}}(X;\hat{\mu}_{y}^{\textup{w}},\hat{\pi}_{m}^{\textup{w}})+\frac{\mathbb{I}(A=a)}{\hat{\pi}_{a}^{\textup{w}}(X)}\big{(}\hat{\mu}_{y}^{\textup{w}}(X,Z)-\hat{\nu}_{y}^{\textup{w}}(X;\hat{\mu}_{y}^{\textup{w}},\hat{\pi}_{m}^{\textup{w}})\big{)}+\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)}{\hat{\pi}_{a}^{\textup{w}}(X)\hat{\pi}_{m}^{\textup{w}}(X,Z)}\big{(}Y-\hat{\mu}_{y}^{\textup{w}}(X,Z)\big{)}\Big{]}$

The triple robustness of $\hat{\psi}_{am}^{\textup{tr}_{1}}$ and
$\hat{\psi}_{am}^{\textup{tr}_{2}}$ and the quadruple robustness of
$\hat{\psi}_{am}^{\textup{qr}}$ are given below.

######  Proposition 2.

Under assumptions 1-3 and suitable regularity conditions,
$\hat{\psi}_{am}^{\textup{tr}_{1}}$ is CAN in
$\mathcal{P}_{1}\cup\mathcal{P}_{2}\cup\mathcal{P}_{3}$,
$\hat{\psi}_{am}^{\textup{tr}_{2}}$ is CAN in
$\mathcal{P}_{1}\cup\mathcal{P}_{3}\cup\mathcal{P}_{4}$, and
$\hat{\psi}_{am}^{\textup{qr}}$ is CAN in
$\mathcal{P}_{1}\cup\mathcal{P}_{2}\cup\mathcal{P}_{3}\cup\mathcal{P}_{4}$. In
addition, $\hat{\psi}_{am}^{\textup{tr}_{1}}$,
$\hat{\psi}_{am}^{\textup{tr}_{2}}$, and $\hat{\psi}_{am}^{\textup{qr}}$ are
all locally efficient in the sense that they attain the semiparametric
efficiency bound of $\mathcal{P}_{\textup{np}}$ at
$\mathcal{P}_{1}\cap\mathcal{P}_{3}$, i.e., when all of the four nuisance
functions are correctly specified.

The multiple robustness of these estimators is due to a similar logic to that
of the doubly robust estimators given previously. For example, for
$\hat{\psi}_{am}^{\textup{qr}}$, when $\mu_{y}^{\textup{w}}(x,z)$ and
$\pi_{a}^{\textup{w}}(x)$ are correctly specified ($\mathcal{P}_{1}$), the
terms involving
$\hat{\nu}_{y}^{\textup{w}}(X;\hat{\mu}_{y}^{\textup{w}},\hat{\pi}_{m}^{\textup{w}})$
and the third term inside $\mathbb{P}_{n}[\cdot]$ will both have a zero mean
(asymptotically), leaving only
$\mathbb{P}_{n}[\big{(}\mathbb{I}(A=a)/\hat{\pi}_{a}^{\textup{w}}(X)\big{)}\hat{\mu}_{y}^{\textup{w}}(X,Z)]$,
the imputation-then-weighting estimator; when $\mu_{y}^{\textup{w}}(x,z)$ and
$\nu_{y}^{\textup{w}}(x)$ are correctly specified ($\mathcal{P}_{2}$), both
the second and third terms inside $\mathbb{P}_{n}[\cdot]$ will have a zero
mean, leaving only
$\mathbb{P}_{n}[\hat{\nu}_{y}^{\textup{w}}(X;\hat{\mu}_{y}^{\textup{w}},\hat{\pi}_{m}^{\textup{w}})]$,
i.e., the doubly robust estimator $\hat{\psi}_{am}^{\textup{dr}_{4}}$ ; when
$\pi_{m}^{\textup{w}}(x,z)$ and $\pi_{a}^{\textup{w}}(x)$ are correctly
specified ($\mathcal{P}_{3}$), the terms involving
$\hat{\mu}_{y}^{\textup{w}}(X,Z)$ and
$\hat{\nu}_{y}^{\textup{w}}(X;\hat{\mu}_{y}^{\textup{w}},\hat{\pi}_{m}^{\textup{w}})$
will both have a zero mean, leaving only
$\mathbb{P}_{n}\big{[}\big{(}\mathbb{I}(A=a)\mathbb{I}(M=m)\big{)}Y/\big{(}\hat{\pi}_{a}^{\textup{w}}(X)\hat{\pi}_{m}^{\textup{w}}(X,Z)\big{)}\big{]}$,
i.e., the pure weighting estimator; and when $\pi_{m}^{\textup{w}}(x,z)$ and
$\nu_{y}^{\textup{w}}(x)$ are correctly specified ($\mathcal{P}_{4}$), the
terms involving $\hat{\mu}_{y}^{\textup{w}}(X,Z)$ and
$\hat{\pi}_{a}^{\textup{w}}(X)$ will both have a zero mean, leaving only
$\mathbb{P}_{n}[\hat{\nu}_{y}^{\textup{w}}(X;\hat{\mu}_{y}^{\textup{w}},\hat{\pi}_{m}^{\textup{w}})]$,
i.e., the doubly robust estimator $\hat{\psi}_{am}^{\textup{dr}_{4}}$.

The asymptotic efficiency of these estimators is due to the fact that they all
solve the estimating equation formed by the efficient influence function of
$\psi_{am}$, which is

$\varphi_{am}(O)=\nu_{y}(X)+\frac{\mathbb{I}(A=a)}{\pi_{a}(X)}\big{(}\mu_{y}(X,Z)-\nu_{y}(X)\big{)}+\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)}{\pi_{a}(X)\pi_{m}(X,Z)}\big{(}Y-\mu_{y}(X,Z)\big{)}-\psi_{am},$
(10)

and the fact that $\mathbb{E}[\varphi_{am}(O;\eta)]$ has a zero derivative
with respect to the nuisance functions
$\eta=(\mu_{y}^{\textup{w}}(x,z),\nu_{y}^{\textup{w}}(x),\pi_{a}^{\textup{w}}(x),\pi_{m}^{\textup{w}}(x,z))$
at the truth (for a derivation of this influence function in the context of
time-varying treatments, see van der Laan and Gruber 2012). The latter
property implies that first step estimation of the nuisance functions will
have no (first-order) effect on the influence function of
$\hat{\psi}_{am}^{\textup{tr}_{1}}$, $\hat{\psi}_{am}^{\textup{tr}_{2}}$, and
$\hat{\psi}_{am}^{\textup{qr}}$. In practice, the nuisance functions can be
estimated via data-adaptive methods combined with cross-fitting (Zheng and van
der Laan 2011; Chernozhukov et al. 2018), in which case estimates of
$\psi_{am}$ (and hence CDE) are semiparametric efficient when estimates of the
nuisance functions, for example, all converge at faster-than-$n^{-1/4}$
rates.222More precisely, $\hat{\psi}_{am}^{\textup{tr}_{1}}$,
$\hat{\psi}_{am}^{\textup{tr}_{2}}$, and $\hat{\psi}_{am}^{\textup{qr}}$ are
semiparametric efficient if
$R_{n}(\hat{\pi}_{a})R_{n}(\hat{\nu}_{y})+R_{n}(\hat{\pi}_{m})R_{n}(\hat{\mu}_{y})=o(n^{-1/2})$,
where $R_{n}(\cdot)$ maps a nuisance function to its $L_{2}(P)$ convergence
rate with respect to the true distribution $P$. See Supporting Material C or
Rotnitzky et al. (2017).

Among the above estimators, $\hat{\psi}_{am}^{\textup{tr}_{1}}$ is akin to the
estimator proposed by Bang and Robins (2005) for the mean of a potential
outcome with time-varying treatments and time-varying confounders.
Specifically, they suggest that
$\mathbb{I}(A=a)\mathbb{I}(M=m)Y/\big{(}\hat{\pi}_{a}^{\textup{w}}(X)\hat{\pi}_{m}^{\textup{w}}(X,Z)\big{)}$
be included as a covariate in a generalized linear model (with canonical link)
for $\hat{\mu}_{y}^{\textup{w}}(x,z)$, and
$\mathbb{I}(A=a)/\hat{\pi}_{a}^{\textup{w}}(X)$ be included as a covariate in
a generalized linear model (with canonical link) for
$\hat{\nu}_{y}^{\textup{w}}(x;\hat{\mu}_{y}^{\textup{w}})$, in which case the
score equations ensure that both the second and third terms inside
$\mathbb{P}_{n}[\cdot]$ have a zero sample mean, thus leaving only
$\mathbb{P}_{n}[\hat{\nu}_{y}^{\textup{w}}(X;\hat{\mu}_{y}^{\textup{w}})]$,
i.e., the pure imputation estimator. Because this procedure estimates
$\psi_{am}$ as a sample average of
$\hat{\nu}_{y}^{\textup{w}}(X;\hat{\mu}_{y}^{\textup{w}})$, which typically
resides in the parameter space of $\psi_{am}$, it tends to be more stable in
finite samples than the unadjusted estimator
$\hat{\psi}_{am}^{\textup{tr}_{1}}$ (Robins et al. 2007).
$\hat{\psi}_{am}^{\textup{tr}_{2}}$ and $\hat{\psi}_{am}^{\textup{qr}}$ differ
from $\hat{\psi}_{am}^{\textup{tr}_{1}}$ only in the response variable they
use to model $\nu_{y}(X)$: $\hat{\psi}_{am}^{\textup{tr}_{2}}$ uses
$\mathbb{I}(M=m)/\hat{\pi}_{m}^{\textup{w}}(X,Z)$ whereas
$\hat{\psi}_{am}^{\textup{qr}}$ uses
$\hat{\mu}_{y}^{\textup{w}}(X,Z)+\mathbb{I}(M=m)(Y-\hat{\mu}_{y}^{\textup{w}}(X,Z))/\hat{\pi}_{m}^{\textup{w}}(X,Z)$,
which adds another layer of robustness. In fact,
$\hat{\psi}_{am}^{\textup{qr}}$ constitutes a special case of the
$2^{K}$-robust estimator proposed by Rotnitzky et al. (2017) in the context of
time-varying treatments. In practice, the Bang-Robins procedure can also be
applied to $\hat{\psi}_{am}^{\textup{tr}_{2}}$ and
$\hat{\psi}_{am}^{\textup{qr}}$ to improve their finite-sample performance.

When flexible machine learning methods (instead of generalized linear models)
are used to estimate the nuisance functions, the Bang-Robins procedure can no
longer ensure a zero sample mean of the second and third terms inside
$\mathbb{P}_{n}[\cdot]$. In this case, the method of targeted maximum
likelihood estimation (TMLE; van Der Laan and Rubin 2006) can be used to
adjust the first step estimates of $\mu_{y}(x,z)$ and $\nu_{y}(x)$ such that
the second and third terms inside $\mathbb{P}_{n}[\cdot]$ have a zero sample
mean. This approach may yield better finite-sample performance than the
unadjusted estimators and more robustness than the Bang-Robins procedure based
on generalized linear models.

For inference of $\hat{\psi}_{am}^{\textup{tr}_{1}}$,
$\hat{\psi}_{am}^{\textup{tr}_{2}}$, $\hat{\psi}_{am}^{\textup{qr}}$, and the
corresponding estimates of the CDE, the nonparametric bootstrap can be used
when the nuisance functions are estimated using parametric models. When data-
adaptive methods are used to estimate the nuisance functions, it will be
reasonable to use the empirical analog of the efficient influence function to
construct standard errors and Wald-type confidence intervals. For example, the
variance of $\widehat{\textup{CDE}}(a,a^{\prime},m)$ can be estimated by
$\mathbb{P}_{n}[(\hat{\varphi}_{am}(O)-\hat{\varphi}_{a^{\prime}m}(O))^{2}]/n$.
When the CDE is defined on the risk ratio or odds ratio scale, corresponding
variance estimates can be obtained using the delta method.

Figure 1: Sampling distributions of the doubly and multiply robust estimators
for $n=2,000$.

##  6 Monte Carlo Evidence

We now present a simulation study to demonstrate the multiple robustness of
the proposed estimators. The data generating process is similar to that used
in Miles et al. (2020) and is described in greater detail in Supporting
Material D. We generate 1,000 Monte Carlo samples of size 2,000, and, without
loss of generality, focus on the estimand $\psi_{01}=\mathbb{E}[Y(0,1)]$. We
examine the sampling distributions of all of the doubly and multiply robust
estimators described above under conditions associated with submodels
$\mathcal{P}_{1}$, $\mathcal{P}_{2}$, $\mathcal{P}_{3}$, $\mathcal{P}_{4}$.
The results are shown in Figure 1, where each panel corresponds to an
estimator, and the y axis is recentered at the true value of $\psi_{01}$. The
shaded box plots highlight the cases under which a given estimator should be
consistent. We can see that all of the doubly and multiply robust estimators
behave as expected. They center around the true value if and only if the
requisite nuisance functions are all correctly specified.

## References

  * Avin et al. (2005) Avin, C., Shpitser, I., Pearl, J., 2005.  Identifiability of path-specific effects, in: Proceedings of the 19th International Joint Conference on Artificial Intelligence, Morgan Kaufmann Publishers Inc.. pp. 357–363.
  * Bang and Robins (2005) Bang, H., Robins, J.M., 2005\.  Doubly robust estimation in missing data and causal inference models.  Biometrics 61, 962–973.
  * Chernozhukov et al. (2018) Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., Robins, J., 2018.  Double/debiased machine learning for treatment and structural parameters.  The Econometrics Journal 21, C1–C68.
  * van Der Laan and Rubin (2006) van Der Laan, M.J., Rubin, D., 2006\.  Targeted maximum likelihood learning.  The International Journal of Biostatistics 2\.
  * Goetgeluk et al. (2009) Goetgeluk, S., Vansteelandt, S., Goetghebeur, E., 2009.  Estimation of controlled direct effects.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 70, 1049–1066.
  * Han and Wang (2013) Han, P., Wang, L., 2013.  Estimation with missing data: Beyond double robustness.  Biometrika 100, 417–430.
  * Joffe and Greene (2009) Joffe, M.M., Greene, T., 2009\.  Related causal frameworks for surrogate outcomes.  Biometrics 65, 530–538.
  * van der Laan and Gruber (2012) van der Laan, M.J., Gruber, S., 2012\.  Targeted minimum loss based estimation of causal effects of multiple time point interventions.  The International Journal of Biostatistics 8\.
  * Miles et al. (2020) Miles, C.H., Shpitser, I., Kanki, P., Meloni, S., Tchetgen Tchetgen, E.J., 2020.  On semiparametric estimation of a path-specific effect in the presence of mediator-outcome confounding.  Biometrika 107, 159–172.
  * Molina et al. (2017) Molina, J., Rotnitzky, A., Sued, M., Robins, J., 2017\.  Multiple robustness in factorized likelihood models.  Biometrika 104, 561–581.
  * Pearl (2001) Pearl, J., 2001.  Direct and indirect effects, in: Proceedings of the Seventeenth Conference on Uncertainty in Artificial Intelligence, Morgan Kaufmann Publishers Inc.. pp. 411–420.
  * Robins (1986) Robins, J., 1986.  A new approach to causal inference in mortality studies with a sustained exposure period-application to control of the healthy worker survivor effect.  Mathematical Modelling 7, 1393–1512.
  * Robins et al. (2007) Robins, J., Sued, M., Lei-Gomez, Q., Rotnitzky, A., 2007\.  Comment: Performance of double-robust estimators when “inverse probability” weights are highly variable.  Statistical Science 22, 544–559.
  * Robins (1999) Robins, J.M., 1999.  Marginal structural models versus structural nested models as tools for causal inference.  Statistical Models in Epidemiology: The Environment and Clinical Trials .
  * Robins and Greenland (1992) Robins, J.M., Greenland, S., 1992\.  Identifiability and exchangeability for direct and indirect effects.  Epidemiology 3, 143–155.
  * Rotnitzky et al. (2017) Rotnitzky, A., Robins, J., Babino, L., 2017.  On the multiply robust estimation of the mean of the g-functional.  arXiv preprint arXiv:1705.08582 .
  * Scharfstein et al. (1999) Scharfstein, D.O., Rotnitzky, A., Robins, J.M., 1999.  Adjusting for nonignorable drop-out using semiparametric nonresponse models.  Journal of the American Statistical Association 94, 1096–1120.
  * VanderWeele (2009) VanderWeele, T.J., 2009.  Marginal structural models for the estimation of direct and indirect effects.  Epidemiology 20, 18–26.
  * VanderWeele and Vansteelandt (2009) VanderWeele, T.J., Vansteelandt, S., 2009\.  Conceptual issues concerning mediation, interventions and composition.  Statistics and its Interface 2, 457–468.
  * Vansteelandt (2009) Vansteelandt, S., 2009.  Estimating direct effects in cohort and case–control studies.  Epidemiology 20, 851–860.
  * Zheng and Zhou (2015) Zheng, C., Zhou, X.H., 2015\.  Causal mediation analysis in the multilevel intervention and multicomponent mediator case.  Journal of the Royal Statistical Society: Series B: Statistical Methodology 77, 581–615.
  * Zheng and van der Laan (2011) Zheng, W., van der Laan, M.J., 2011\.  Cross-validated targeted minimum-loss-based estimation, in: Targeted Learning. Springer, New York, NY, pp. 459–474.
  * Zhou and Wodtke (2019) Zhou, X., Wodtke, G.T., 2019\.  A regression-with-residuals method for estimating controlled direct effects.  Political Analysis 27, 360–369.

##  Appendix A Proof of Equations (2-5).

Starting from equation (6), we have

$\displaystyle\mathbb{E}_{X}\mathbb{E}\big{[}\frac{\mathbb{I}(M=m)Y}{\Pr[M=m|X,A,Z]}|X,A=a\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{X}\Big{[}\mathbb{E}\big{[}\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)Y}{\Pr[A=a|X]\Pr[M=m|X,A,Z]}|X,A=a\big{]}\cdot\Pr[A=a|X]+0\cdot\Pr[A=a|X]\Big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}[\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)Y}{\Pr[A=a|X]\Pr[M=m|X,A,Z]}]$
(11) $\displaystyle=$
$\displaystyle\mathbb{E}\big{[}\frac{\mathbb{I}(A=a)}{\Pr[A=a|X]}\mathbb{E}[\frac{\mathbb{I}(M=m)Y}{\Pr[M=m|X,A,Z]}|X,A,Z]\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\big{[}\frac{\mathbb{I}(A=a)}{\Pr[A=a|X]}\mathbb{E}[\frac{\mathbb{I}(M=m)Y}{\Pr[M=m|X,A,Z]}|X,A,Z,M=m]\Pr[M=m|X,A=a,Z]\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\big{[}\frac{\mathbb{I}(A=a)\mathbb{E}[Y|X,A,Z,M=m]}{\Pr[A=a|X]}\big{]}$
(12) $\displaystyle=$
$\displaystyle\mathbb{E}_{X,A}\mathbb{E}\big{[}\frac{\mathbb{I}(A=a)\mathbb{E}[Y|X,A,Z,M=m]}{\Pr[A=a|X]}|X,A\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{X,A}\Big{[}\mathbb{E}\big{[}\frac{\mathbb{I}(A=a)\mathbb{E}[Y|X,A,Z,M=m]}{\Pr[A=a|X]}|X,A=a\big{]}\Pr[A=a|X]\Big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{X}\mathbb{E}\big{[}\mathbb{E}[Y|X,A,Z,M=m]|X,A=a\big{]}$
(13) $\displaystyle=$
$\displaystyle\iiint\mathbb{E}[Y|x,a,z,m]dP(z|x,a)dP(x).$ (14)

Equations (11), (12), (13), and (14) correspond to equations (5), (4), (3),
and (2), respectively.

##  Appendix B Proof of Proposition 1

Below we show that $\hat{\psi}_{am}^{\textup{dr}_{1}}$ is CAN when (a)
$\mu_{y}^{\textup{w}}(x,z)$ is correctly specified and (b) either
$\nu_{y}^{\textup{w}}(x)$ or $\pi_{a}^{\textup{w}}(x)$ is correctly specified.
The double robustness of $\hat{\psi}_{am}^{\textup{dr}_{2}}$,
$\hat{\psi}_{am}^{\textup{dr}_{3}}$, and $\hat{\psi}_{am}^{\textup{dr}_{4}}$
can be verified analogously.

A first-order Taylor expansion of $\hat{\psi}_{am}^{\textup{dr}_{1}}$ implies
that

$\hat{\psi}_{am}^{\textup{dr}_{1}}=\mathbb{P}_{n}\big{[}\nu_{y}^{*}(X)+\frac{\mathbb{I}(A=a)}{\pi_{a}^{*}(X)}\big{(}\mu_{y}^{*}(X,Z)-\nu_{y}^{*}(X)\big{)}\big{]}+o_{p}(1),$

where $\nu_{y}^{*}(X)$, $\pi_{a}^{*}(X)$ and $\mu_{y}^{*}(X,Z)$ denote the
probability limits of $\hat{\nu}_{y}^{\textup{w}}(X)$,
$\hat{\pi}_{a}^{\textup{w}}(X)$, and $\hat{\mu}_{y}^{\textup{w}}(X,Z)$. Hence
it suffices to show
$\mathbb{E}\big{[}\nu_{y}^{*}(X)+\frac{\mathbb{I}(A=a)}{\pi_{a}^{*}(X)}\big{(}\mu_{y}^{*}(X,Z)-\nu_{y}^{*}(X)\big{)}\big{]}=\psi_{am}$
if $\mu_{y}^{*}(X,Z)=\mu_{y}(X,Z)$ and either $\pi_{a}^{*}(X)=\pi_{a}(X)$ or
$\nu_{y}^{*}(X)=\nu_{y}(X)$. Consistency follows from the law of large
numbers, and asymptotic normality follows from standard regularity conditions
for M-estimators.

When $\mu_{y}^{*}(X,Z)=\mu_{y}(X,Z)$ and $\pi_{a}^{*}(X)=\pi_{a}(X)$,

$\displaystyle\textup{plim}\,\hat{\psi}_{am}^{\textup{dr}_{1}}=$
$\displaystyle\mathbb{E}\big{[}\nu_{y}^{*}(X)+\frac{\mathbb{I}(A=a)}{\pi_{a}(X)}\big{(}\mu_{y}(X,Z)-\nu_{y}^{*}(X)\big{)}\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{X}\mathbb{E}\big{[}\nu_{y}^{*}(X)+\frac{\mathbb{I}(A=a)}{\pi_{a}(X)}\big{(}\mu_{y}(X,Z)-\nu_{y}^{*}(X)\big{)}|X\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{X}\big{[}\nu_{y}^{*}(X)+\mathbb{E}[\mu_{y}(X,Z)-\nu_{y}^{*}(X)|X,A=a]\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{X}\big{[}\nu_{y}^{*}(X)+\mathbb{E}[\mu_{y}(X,Z)|X,A=a]-\nu_{y}^{*}(X)\big{]}$
$\displaystyle=$ $\displaystyle\mathbb{E}_{X}\mathbb{E}[\mu_{y}(X,Z)|X,A=a]$
$\displaystyle=$ $\displaystyle\psi_{am}.$

When $\mu_{y}^{*}(X,Z)=\mu_{y}(X,Z)$ and $\nu_{y}^{*}(X)=\nu_{y}(X)$,

$\displaystyle\textup{plim}\,\hat{\psi}_{am}^{\textup{dr}_{1}}=$
$\displaystyle\mathbb{E}\big{[}\nu_{y}(X)+\frac{\mathbb{I}(A=a)}{\pi_{a}^{*}(X)}\big{(}\mu_{y}(X,Z)-\nu_{y}(X)\big{)}\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{X}\mathbb{E}\big{[}\nu_{y}(X)+\frac{\mathbb{I}(A=a)}{\pi_{a}^{*}(X)}\big{(}\mu_{y}(X,Z)-\nu_{y}(X)\big{)}|X\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{X}\big{[}\nu_{y}(X)+\mathbb{E}[\frac{\mathbb{I}(A=a)}{\pi_{a}^{*}(X)}\big{(}\mu_{y}(X,Z)-\nu_{y}(X)\big{)}|X,A=a]\pi_{a}(X)\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{X}\big{[}\nu_{y}(X)+\frac{\pi_{a}(X)}{\pi_{a}^{*}(X)}\mathbb{E}[\mu_{y}(X,Z)-\nu_{y}(X)|X,A=a]\big{]}$
$\displaystyle=$ $\displaystyle\mathbb{E}_{X}\big{[}\nu_{y}(X)+0\big{]}$
$\displaystyle=$ $\displaystyle\psi_{am}.$

##  Appendix C Proof of Proposition 2

Below we show that $\hat{\psi}_{am}^{\textup{tr}_{1}}$ is CAN in
$\mathcal{P}_{1}\cup\mathcal{P}_{2}\cup\mathcal{P}_{3}$ and locally efficient
in $\mathcal{P}_{1}\cap\mathcal{P}_{3}$. The multiple robustness and local
efficiency of $\hat{\psi}_{am}^{\textup{tr}_{2}}$ and
$\hat{\psi}_{am}^{\textup{qr}}$ can be verified analogously.

A first-order Taylor expansion of $\hat{\psi}_{am}^{\textup{tr}_{1}}$ implies
that

$\hat{\psi}_{am}^{\textup{tr}_{1}}=\mathbb{P}_{n}\big{[}\nu_{y}^{*}(X)+\frac{\mathbb{I}(A=a)}{\pi_{a}^{*}(X)}\big{(}\mu_{y}^{*}(X,Z)-\nu_{y}^{*}(X)\big{)}+\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)}{\pi_{a}^{*}(X)\pi_{m}^{*}(X,Z)}\big{(}Y-\mu_{y}^{*}(X,Z)\big{)}\big{]}+o_{p}(1),$

where $\mu_{y}^{*}(X,Z)$, $\nu_{y}^{*}(X)$, $\pi_{a}^{*}(X)$, and
$\pi_{m}^{*}(X,Z)$ denote the probability limits of
$\hat{\mu}_{y}^{\textup{w}}(X,Z)$, $\hat{\nu}_{y}^{\textup{w}}(X)$,
$\hat{\pi}_{a}^{\textup{w}}(X)$, and $\hat{\pi}_{m}(X,Z)$. Hence it suffices
to show that the expectation of the quantity inside $\mathbb{P}_{n}[\cdot]$
equals $\psi_{am}$ in $\mathcal{P}_{1}\cup\mathcal{P}_{2}\cup\mathcal{P}_{3}$.
Consistency follows from the law of large numbers, and asymptotic normality
follows from standard regularity conditions for M-estimators. First, consider
submodel $\mathcal{P}_{2}$, under which we have
$\mu_{y}^{*}(X,Z)=\mu_{y}(X,Z)$ and $\nu_{y}^{*}(X)=\nu_{y}(X)$. From the
proof of proposition 1, we know that

$\mathbb{E}\Big{[}\nu_{y}(X)+\frac{\mathbb{I}(A=a)}{\pi_{a}^{*}(X)}\big{(}\mu_{y}(X,Z)-\nu_{y}(X)\big{)}\Big{]}=\psi_{am}.$

Thus

$\displaystyle\textup{plim}\,\hat{\psi}_{am}^{\textup{tr}_{1}}$
$\displaystyle=\psi_{am}+\mathbb{E}\big{[}\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)}{\pi_{a}^{*}(X)\pi_{m}^{*}(X,Z)}\big{(}Y-\mu_{y}(X,Z)\big{)}\big{]}$
$\displaystyle=\psi_{am}+\mathbb{E}_{X}\mathbb{E}\big{[}\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)}{\pi_{a}^{*}(X)\pi_{m}^{*}(X,Z)}\big{(}Y-\mu_{y}(X,Z)\big{)}\big{|}X\big{]}$
$\displaystyle=\psi_{am}+\mathbb{E}_{X}\Big{[}\mathbb{E}\big{[}\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)}{\pi_{a}^{*}(X)\pi_{m}^{*}(X,Z)}\big{(}Y-\mu_{y}(X,Z)\big{)}\big{)}|X,A=a\big{]}\pi_{a}(X)\Big{]}$
$\displaystyle=\psi_{am}+\mathbb{E}_{X}\Big{[}\frac{\pi_{a}(X)}{\pi_{a}^{*}(X)}\mathbb{E}\big{[}\frac{\mathbb{I}(M=m)}{\pi_{m}^{*}(X,Z)}\big{(}Y-\mu_{y}(X,Z)\big{)}|X,A=a\big{]}\Big{]}$
$\displaystyle=\psi_{am}+\mathbb{E}_{X}\Big{[}\frac{\pi_{a}(X)}{\pi_{a}^{*}(X)}\mathbb{E}_{Z|X,A=a}\mathbb{E}\big{[}\frac{\mathbb{I}(M=m)\pi_{m}(X,Z)}{\pi_{m}^{*}(X,Z)}\big{(}Y-\mu_{y}(X,Z)\big{)}|X,A=a,Z,M=m\big{]}\Big{]}$
$\displaystyle=\psi_{am}+\mathbb{E}_{X}\Big{[}\frac{\pi_{a}(X)}{\pi_{a}^{*}(X)}\mathbb{E}_{Z|X,A=a}\frac{\pi_{m}(X,Z)}{\pi_{m}^{*}(X,Z)}\underbrace{\mathbb{E}\big{[}\big{(}Y-\mu_{y}(X,Z)\big{)}|X,A=a,Z,M=m\big{]}}_{=0}\Big{]}$
$\displaystyle=\psi_{am}.$

Then, under $\mathcal{P}_{1}$, we have $\mu_{y}^{*}(X,Z)=\mu_{y}(X,Z)$ and
$\pi_{a}^{*}(X)=\pi_{a}(X)$. From the proof of proposition 1, we know that

$\mathbb{E}\Big{[}\nu_{y}^{*}(X)+\frac{\mathbb{I}(A=a)}{\pi_{a}(X)}\big{(}\mu_{y}(X,Z)-\nu_{y}^{*}(X)\big{)}\big{]}=\psi_{am}.$

Thus
$\textup{plim}\,\hat{\psi}_{am}^{\textup{tr}_{1}}=\psi_{am}+\mathbb{E}\big{[}\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)}{\pi_{a}^{*}(X)\pi_{m}^{*}(X,Z)}\big{(}Y-\mu_{y}(X,Z)\big{)}\big{]}=0$
(directly from the above proof for submodel $\mathcal{P}_{2}$). Finally, under
$\mathcal{P}_{3}$, we have $\pi_{a}^{*}(X)=\pi_{a}(X)$ and
$\pi_{m}^{*}(X,Z)=\pi_{m}(X,Z)$.

$\displaystyle\textup{plim}\,\hat{\psi}_{am}^{\textup{tr}_{1}}$
$\displaystyle=\mathbb{E}\big{[}\nu_{y}^{*}(X)+\frac{\mathbb{I}(A=a)}{\pi_{a}(X)}\big{(}\mu_{y}^{*}(X,Z)-\nu_{y}^{*}(X)\big{)}+\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)}{\pi_{a}(X)\pi_{m}(X,Z)}\big{(}Y-\mu_{y}^{*}(X,Z)\big{)}\big{]}$
$\displaystyle=\mathbb{E}_{X}[\nu_{y}^{*}(X)]+\mathbb{E}_{X}\Big{[}\mathbb{E}\big{[}\frac{\mathbb{I}(A=a)}{\pi_{a}(X)}\big{(}\mu_{y}^{*}(X,Z)-\nu_{y}^{*}(X)\big{)}|X,A=a\big{]}\pi_{a}(X)\Big{]}$
$\displaystyle+\mathbb{E}_{X}\mathbb{E}_{Z|X,A=a}\Big{[}\mathbb{E}\big{[}\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)}{\pi_{a}(X)\pi_{m}(X,Z)}\big{(}Y-\mu_{y}^{*}(X,Z)\big{)}|X,A=a,Z\big{]}\pi_{a}(X)\Big{]}$
$\displaystyle=\mathbb{E}_{X}[\nu_{y}^{*}(X)]+\mathbb{E}_{X}\mathbb{E}\big{[}\big{(}\mu_{y}^{*}(X,Z)-\nu_{y}^{*}(X)\big{)}|X,A=a\big{]}$
$\displaystyle+\mathbb{E}_{X}\mathbb{E}_{Z|X,A=a}\Big{[}\mathbb{E}\big{[}\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)}{\pi_{a}(X)\pi_{m}(X,Z)}\big{(}Y-\mu_{y}^{*}(X,Z)\big{)}|X,A=a,Z,M=m\big{]}\pi_{m}(X,Z)\pi_{a}(X)\Big{]}$
$\displaystyle=\mathbb{E}_{X}\Big{[}\mu_{y}^{*}(X,Z)+\mathbb{E}_{Z|X,A=a}\mathbb{E}\big{[}\big{(}Y-\mu_{y}^{*}(X,Z)\big{)}|X,A=a,Z,M=m\big{]}\Big{]}$
$\displaystyle=\mathbb{E}_{X}\Big{[}\mu_{y}^{*}(X,Z)+\mu_{y}(X,Z)-\mu_{y}^{*}(X,Z)\Big{]}$
$\displaystyle=\psi_{am}.$

To show that $\hat{\psi}_{am}^{\textup{tr}_{1}}$ is locally efficient, we
first verify that equation (10) is the efficient influence function of
$\psi_{am}$ in $\mathcal{P}_{\textup{np}}$, i.e.,

$\frac{\partial\psi_{am}(t)}{\partial
t}\biggl{|}_{t=0}=\mathbb{E}[\varphi_{a,m}^{\textup{eff}}(O)S_{0}(O)],$ (15)

where $S_{0}(O)$ is the score function for any one-dimensional submodel
$P_{t}(O)$ evaluated at $t=0$. We first note that $S_{t}(O)$ can be written as
$S_{t}(O)=S_{t}(X)+S_{t}(A|X)+S_{t}(Z|X,A)+S_{t}(M|X,A,Z)+S_{t}(Y|X,A,Z,M)$,
where $S_{t}(u|v)=\partial\log p_{t}(u|v)/\partial t$ and $p_{t}(u|v)$ is the
conditional probability density/mass function of $U$ given $V$. Using equation
(1) and the product rule, the left hand side of equation (15) can be written
as

$\displaystyle\frac{\partial\psi_{am}(t)}{\partial t}\biggl{|}_{t=0}$
$\displaystyle=\frac{\partial\iiint
ydP_{t}(y|x,a,z,m)dP_{t}(z|x,a)dP_{t}(x)}{\partial t}\biggl{|}_{t=0}$
$\displaystyle=\underbrace{\iiint
yS_{0}(x)dP_{0}(y|x,a,z,m)dP_{0}(z|x,a)dP_{0}(x)}_{=:\textup{$\phi_{1}$}}$
$\displaystyle+\underbrace{\iiint
yS_{0}(z|x,a)dP_{0}(y|x,a,z,m)dP_{0}(z|x,a)dP_{0}(x)}_{=:\phi_{2}}$
$\displaystyle+\underbrace{\iiint
yS_{0}(y|x,a,z,m)dP_{0}(y|x,a,z,m)dP_{0}(z|x,a)dP_{0}(x)}_{=:\phi_{3}}$
$\displaystyle=\phi_{1}+\phi_{2}+\phi_{3}$

where the second equality follows from the fact that $\partial
dP_{t}(u|v)/\partial t=S_{t}(u|v)dP_{t}(u|v).$

Before evaluating the right hand side of equation (15), we introduce the
following shorthands:

$\displaystyle\varphi_{1}(X)$ $\displaystyle=\nu_{y}(X),$
$\displaystyle\varphi_{2}(X,A,Z)$
$\displaystyle=\frac{\mathbb{I}(A=a)}{\pi_{a}(X)}\big{(}\mu_{y}(X,Z)-\nu_{y}(X)\big{)},$
$\displaystyle\varphi_{3}(X,A,Z,M,Y)$
$\displaystyle=\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)}{\pi_{a}(X)\pi_{m}(X,Z)}\big{(}Y-\mu_{y}(X,Z)\big{)}.$

Thus
$\varphi_{a,m}^{\textup{eff}}(O)=\varphi_{1}(X)+\varphi_{2}(X,A,Z)+\varphi_{3}(X,A,Z,M,Y)-\psi_{am}$.
We first observe

$\displaystyle\mathbb{E}[\varphi_{1}(X)S_{0}(O)]$ $\displaystyle=$
$\displaystyle\mathbb{E}[\varphi_{1}(X)\big{(}S_{0}(X)+S_{0}(A|X)+S_{0}(Z|X,A)+S_{0}(M|X,A,Z)+S_{0}(Y|X,A,Z,M)\big{)}]$
$\displaystyle=$
$\displaystyle\mathbb{E}[\varphi_{1}(X)S_{0}(X)]+\mathbb{E}[\varphi_{1}(X)S_{0}(A|X)]\ldots+\mathbb{E}[\varphi_{1}(X)S_{0}(Y|X,A,Z,M)]$
$\displaystyle=$
$\displaystyle\mathbb{E}[\varphi_{1}(X)S_{0}(X)]+\mathbb{E}\big{[}\varphi_{1}(X)\underbrace{\mathbb{E}[S_{0}(A|X)|X]}_{=0}\big{]}\ldots+\mathbb{E}\big{[}\varphi_{1}(X)\underbrace{\mathbb{E}[S_{0}(Y|X,A,Z,M)|X,A,Z,M]}_{=0}\big{]}$
$\displaystyle=$ $\displaystyle\mathbb{E}[\varphi_{1}(X)S_{0}(X)]$
$\displaystyle=$
$\displaystyle\mathbb{E}[\mathbb{E}_{Z|X,A=a}\mathbb{E}[Y|X,A=a,Z,M=m]S_{0}(X)]$
$\displaystyle=$ $\displaystyle\phi_{1}$ (16)

where we used the fact that $\mathbb{E}[S(U|V)|V]=0$ for any score function
$S(U,V)$. Second,

$\displaystyle\mathbb{E}[\varphi_{2}(X,A,Z)S_{0}(O)]$ $\displaystyle=$
$\displaystyle\mathbb{E}[\varphi_{2}(X,A,Z)\big{(}S_{0}(X)+S_{0}(A|X)+S_{0}(Z|X,A)+S_{0}(M|X,A,Z)+S_{0}(Y|X,A,Z,M)\big{)}]$
$\displaystyle=$
$\displaystyle\mathbb{E}[\varphi_{2}(X,A,Z)S_{0}(X)]+\mathbb{E}[\varphi_{2}(X,A,Z)S_{0}(A|X)]+\mathbb{E}[\varphi_{2}(X,A,Z)S_{0}(Z|X,A)]$
$\displaystyle+\mathbb{E}\big{[}\varphi_{2}(X,A,Z)\underbrace{\mathbb{E}[S_{0}(M|X,A,Z)|X,A,Z]}_{=0}\big{]}+\mathbb{E}\big{[}\varphi_{2}(X,A,Z)\underbrace{\mathbb{E}[S_{0}(Y|X,A,Z,M)|X,A,Z,M]}_{=0}\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}[\varphi_{2}(X,A,Z)S_{0}(X)]+\mathbb{E}[\varphi_{2}(X,A,Z)S_{0}(A|X)]+\mathbb{E}[\varphi_{2}(X,A,Z)S_{0}(Z|X,A)]$
$\displaystyle=$
$\displaystyle\mathbb{E}\big{[}S_{0}(X)\underbrace{\mathbb{E}[\varphi_{2}(X,A,Z)|X,A]}_{=0}\big{]}+\mathbb{E}\big{[}S_{0}(A|X)\underbrace{\mathbb{E}[\varphi_{2}(X,A,Z)|X,A]}_{=0}\big{]}+\mathbb{E}[\varphi_{2}(X,A,Z)S_{0}(Z|X,A)]$
$\displaystyle=$ $\displaystyle\mathbb{E}[\varphi_{2}(X,A,Z)S_{0}(Z|X,A)]$
$\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\frac{\mathbb{I}(A=a)}{\Pr[A=a|X]}\big{(}\mathbb{E}[Y|X,A=a,Z,M=m]-\varphi_{1}(X)\big{)}S_{0}(Z|X,A)\Big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{X}\Big{\\{}\mathbb{E}\Big{[}\frac{\mathbb{I}(A=a)}{\Pr[A=a|X]}\big{(}\mathbb{E}[Y|X,A=a,Z,M=m]-\varphi_{1}(X)\big{)}S_{0}(Z|X,A)\Big{|}X,A=a\Big{]}\cdot\Pr[A=a|X]\Big{\\}}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{X}\mathbb{E}\Big{[}\big{(}\mathbb{E}[Y|X,A=a,Z,M=m]-\varphi_{1}(X)\big{)}S_{0}(Z|X,A)\Big{|}X,A=a\Big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{X}\mathbb{E}_{Z|X,A=a}\big{[}\mathbb{E}[Y|X,A=a,Z,M=m]S_{0}(Z|X,A)\big{]}$
$\displaystyle=$ $\displaystyle\phi_{2}$ (17)

where the fifth line follows from the fact that

$\displaystyle\mathbb{E}[\varphi_{2}(X,A,Z)|X,A]$
$\displaystyle=\mathbb{E}[\frac{\mathbb{I}(A=a)}{\Pr[A=a|X]}\big{(}\mathbb{E}[Y|X,A=a,Z,M=m]-\varphi_{1}(X)\big{)}|X,A]$
$\displaystyle=\text{$\frac{\mathbb{I}(A=a)}{\Pr[A=a|X]}$}\big{[}\underbrace{\mathbb{E}_{Z|X,A}\mathbb{E}[Y|X,A=a,Z,M=m]}_{=\varphi_{1}(X)}-\varphi_{1}(X)\big{]}$
$\displaystyle=0.$

Third,

$\displaystyle\mathbb{E}[[\varphi_{3}(X,A,Z,M,Y)S_{0}(O)]$ $\displaystyle=$
$\displaystyle\mathbb{E}[\varphi_{3}(X,A,Z,M,Y)\big{(}S_{0}(X)+S_{0}(A|X)+S_{0}(Z|X,A)+S_{0}(M|X,A,Z)+S_{0}(Y|X,A,Z,M)\big{)}]$
$\displaystyle=$
$\displaystyle\mathbb{E}\big{[}\big{(}S_{0}(X)+S_{0}(A|X)+S_{0}(Z|X,A)+S_{0}(M|X,A,Z)\big{)}\underbrace{\mathbb{E}[\varphi_{3}(X,A,Z,M,Y)|X,A,Z,M]}_{=0}\big{]}$
$\displaystyle+\mathbb{E}[\varphi_{3}(X,A,Z,M,Y)S_{0}(Y|X,A,Z,M)]$ (18)
$\displaystyle=$
$\displaystyle\mathbb{E}[\varphi_{3}(X,A,Z,M,Y)S_{0}(Y|X,A,Z,M)]$
$\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)YS_{0}(Y|X,A,Z,M)}{\Pr[A=a|X]\Pr[M=m|X,A=a,Z]}\Big{]}$
$\displaystyle-\mathbb{E}\Big{[}\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)\mathbb{E}[Y|X,A=a,Z,M=m]}{\Pr[A=a|X]\Pr[M=m|X,A=a,Z]}\underbrace{\mathbb{E}\big{[}S_{0}(Y|X,A,Z,M)|X,A,Z,M\big{]}}_{=0}\Big{]}$
(19) $\displaystyle=$
$\displaystyle\mathbb{E}\Big{\\{}\mathbb{E}\Big{[}\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)YS_{0}(Y|X,A,Z,M)}{\Pr[A=a|X]\Pr[M=m|X,A=a,Z]}\Big{|}X,A,Z,M=m\Big{]}\Pr[M=m|X,A,Z]\Big{\\}}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{X,A}\mathbb{E}\Big{[}\frac{\mathbb{I}(A=a)\Pr[M=m|X,A,Z]\mathbb{E}\big{[}YS_{0}(Y|X,A,Z,M)|X,A,Z,M=m\big{]}}{\Pr[A=a|X]\Pr[M=m|X,A=a,Z]}\Big{|}X,A\Big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{X}\bigg{\\{}\mathbb{E}\Big{[}\frac{\mathbb{I}(A=a)\Pr[M=m|X,A,Z]\mathbb{E}\big{[}YS_{0}(Y|X,A,Z,M)|X,A,Z,M=m\big{]}}{\Pr[A=a|X]\Pr[M=m|X,A=a,Z]}\Big{|}X,A=a\Big{]}\Pr[A=a|X]\bigg{\\}}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{X}\bigg{\\{}\mathbb{E}\Big{[}\mathbb{E}\Big{[}YS_{0}(Y|X,A,Z,M)\Big{|}X,A,Z,M=m\Big{]}\Big{|}X,A=a\Big{]}\bigg{\\}}$
$\displaystyle=$ $\displaystyle\phi_{3}$ (20)

where the first equality follows from the fact that

$\displaystyle\mathbb{E}[\varphi_{3}(X,A,Z,M,Y)|X,A,Z,M]$ $\displaystyle=$
$\displaystyle\mathbb{E}[\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)\big{(}Y-\mathbb{E}[Y|X,A=a,Z,M=m]\big{)}}{\Pr[A=a|X]\Pr[M=m|X,A=a,Z]}|X,A,Z,M]$
$\displaystyle=$
$\displaystyle\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)\big{(}\mathbb{E}[Y|X,A,Z,M]-\mathbb{E}[Y|X,A=a,Z,M=m]\big{)}}{\Pr[A=a|X]\Pr[M=m|X,A=a,Z]}$
$\displaystyle=$ $\displaystyle 0.$

Summing up equations (16-20), we have

$\displaystyle\mathbb{E}[\varphi_{a,m}^{\textup{eff}}(O)S_{0}(O)]$
$\displaystyle=$
$\displaystyle\mathbb{E}[\big{(}\varphi_{1}(X)+\varphi_{2}(X,A,Z)+\varphi_{3}(X,A,Z,M,Y)-\psi_{am}\big{)}S_{0}(O)]$
$\displaystyle=$
$\displaystyle\mathbb{E}[\varphi_{1}(X)S_{0}(O)]+\mathbb{E}[\varphi_{2}(X,A,Z)S_{0}(O)]+\mathbb{E}[\varphi_{3}(X,A,Z,M,Y)S_{0}(O)]-\psi_{am}\cdot\mathbb{E}[S_{0}(O)]$
$\displaystyle=$ $\displaystyle\phi_{1}+\phi_{2}+\phi_{3}-0$ $\displaystyle=$
$\displaystyle\frac{\partial\psi_{am}(t)}{\partial t}\biggl{|}_{t=0}.$

One way to show the local efficiency of $\hat{\psi}_{am}^{\textup{tr}_{1}}$ is
to verify that $\mathbb{E}[\varphi_{am}(O;\eta)]$ has a zero derivative with
respect to the nuisance functions
$\eta=(\mu_{y}^{\textup{w}}(X,Z),\nu_{y}^{\textup{w}}(X),\pi_{a}^{\textup{w}}(X),\pi_{m}^{\textup{w}}(X,Z))$
at the truth. Suppose these nuisance functions are parameterized by different
components of a vector-valued parameter $\beta$, where $\beta_{0}$ denotes the
truth. We then have

$\displaystyle\frac{\partial\varphi_{am}(O;\eta)}{\partial\beta}\Big{|}_{\beta=\beta_{0}}=$
$\displaystyle\frac{\mathbb{I}(A=a)}{\pi_{a}(X)}\big{[}1-\frac{\mathbb{I}(M=m)}{\pi_{m}(X,Z)}\big{]}\dot{\mu}_{y}^{w}(X,Z)$
$\displaystyle+\big{[}1-\frac{\mathbb{I}(A=a)}{\pi_{a}(X)}\big{]}\dot{\nu}_{y}^{w}(X)$
$\displaystyle-\frac{\mathbb{I}(A=a)}{\pi_{a}^{2}(X)}\big{[}\mu_{y}(X,Z)-\nu_{y}(X)+\frac{\mathbb{I}(M=m)}{\pi_{m}(X,Z)}\big{(}Y-\mu_{y}(X,Z)\big{)}\big{]}\dot{\pi}_{a}(X)$
$\displaystyle-\frac{\mathbb{I}(A=a)\mathbb{I}(M=m)}{\pi_{a}(X)\pi_{m}^{2}(X,Z)}\big{[}Y-\mu_{y}(X,Z)\big{]}\dot{\pi}_{m}(X,Z)$

where $\dot{\mu}_{y}^{w}(X,Z)$, $\dot{\nu}_{y}^{w}(X)$, $\dot{\pi}_{a}(X)$,
and $\dot{\pi}_{m}(X,Z)$ denote the derivatives of the corresponding functions
with respect to $\beta$ at $\beta_{0}$. It is easy to verify that these
components all have a zero mean. Thus
$\mathbb{E}[\partial\varphi_{am}(O;\beta_{0})/\partial\beta]=0$, implying that
the influence function of $\hat{\psi}_{am}^{\textup{tr}_{1}}$ at
$\mathcal{P}_{1}\cap\mathcal{P}_{3}$ is $\varphi_{a,m}^{\textup{eff}}(O)$.

Alternatively, we can also analyze the asymptotic expansion of
$\hat{\psi}_{am}^{\textup{tr}_{1}}$ to establish weaker conditions for its
semiparametric efficiency. Denote
$m(O;\eta)=\varphi_{am}(O;\eta)+\text{$\psi_{am}$}$, we have

$\displaystyle\hat{\psi}_{am}^{\textup{tr}_{1}}$
$\displaystyle=\mathbb{P}_{n}[m(O;\hat{\eta})]$
$\displaystyle=\mathbb{P}_{n}[m(O;\eta)]+P[m(O;\hat{\eta})-m(O;\eta)]+(\mathbb{P}_{n}-P)[m(O;\hat{\eta})-m(O;\eta)],$
(21)

where $Pg=\int gdP$ denotes the expectation of function $g(O)$ taken at the
truth. In equation (21), the first term can be analyzed with the standard
central limit theorem and has an asymptotic variance of
$\mathbb{E}[(\varphi_{am}(O;\eta))^{2}]$. The last term is an empirical
process term that will be $o_{p}(1/\sqrt{n})$ if either the nuisance functions
fall in a Donsker class or if cross-fitting is used to induce independence
between $\hat{\eta}$ and $O$ (Chernozhukov et al. 2018). By rearranging terms,
using the law of iterated expectations, and applying the Cauchy-Schwartz
inequality, we can rewrite the second term as

$\displaystyle P[m(O;\hat{\eta})-m(O;\eta)]$
$\displaystyle=P\Big{[}\frac{(\hat{\pi}_{a}(X)-\pi_{a}(X))(\hat{\nu}_{y}(X)-\nu_{y}(X))}{\hat{\pi}_{a}(X)}\Big{]}$
$\displaystyle+P\Big{[}\frac{\mathbb{I}(A=a)(\hat{\pi}_{m}(X,Z)-\pi_{m}(X,Z))(\hat{\mu}_{y}(X,Z)-\mu_{y}(X,Z))}{\hat{\pi}_{a}(X)\hat{\pi}_{m}(X,Z)}\Big{]}.$
$\displaystyle\leq
C_{1}\|\hat{\pi}_{a}(X)-\pi_{a}(X)\|\cdot\|\hat{\nu}_{y}(X)-\nu_{y}(X)\|$
$\displaystyle+C_{2}\|\hat{\pi}_{m}(X,Z)-\pi_{m}(X,Z)\|\cdot\|\hat{\mu}_{y}(X,Z)-\mu_{y}(X,Z)\|,$

where $C_{1}$ and $C_{2}$ are on the order of $O_{p}(1)$, and $\|g\|=(\int
g^{T}gdP)^{1/2}$. The last line is due to the positivity assumption that
$\pi_{a}(X)$ and $\pi_{m}(X,Z)$ are bounded away from zero. Thus the second
term in equation (21) is asymptotically negligible if
$R_{n}(\hat{\pi}_{a})R_{n}(\hat{\nu}_{y})+R_{n}(\hat{\pi}_{m})R_{n}(\hat{\mu}_{y})=o(n^{-1/2})$,
where $R_{n}(\cdot)$ maps a nuisance function to its $L_{2}(P)$ convergence
rate. This result implies that if all nuisance functions are consistently
estimated and converge at faster than $n^{1/4}$ rates, then
$\hat{\psi}_{am}^{\textup{tr}_{1}}$ is semiparametric efficient.

##  Appendix D More Details of the Simulation Study

The variables $X,A,Z,M,Y$ in the simulation study are generated via the
following model:

$\displaystyle(U_{XA},U_{XZ},U_{XM},U_{XY})$ $\displaystyle\sim N(0,I_{4})$
$\displaystyle X$ $\displaystyle\sim
N((U_{XA},U_{XZ},U_{XM},U_{XY})\beta_{X},1)$ $\displaystyle A$
$\displaystyle\sim\textup{Bernoulli}\big{(}\textup{logit}^{-1}[(1,U_{XA},|X|)\beta_{A}]\big{)}$
$\displaystyle Z$ $\displaystyle\sim
N\big{(}(1,U_{XZ},X,X^{2},A)\beta_{Z},1\big{)}$ $\displaystyle M$
$\displaystyle\sim\textup{Bernoulli}\big{(}\textup{logit}^{-1}[(1,U_{XM},X,X^{2},A,Z,XA,XZ)\beta_{M}]\big{)}$
$\displaystyle Y$ $\displaystyle\sim
N\big{(}(1,U_{XY},X,X^{2},A,Z,XZ,M,AM)\beta_{Y},1\big{)}.$

The coefficients $\beta_{X},\beta_{A},\beta_{Z},\beta_{M},\beta_{Y}$ are
generated using a set of uniform distributions with certain constraints
designed to create nontrivial degrees of model misspecification.

It can be shown that under the above model, the nuisance functions
$\mu_{y}(x,z)$, $\nu_{y}(x)$, $\pi_{a}(x)$, $\pi_{m}(x,z)$ can be consistently
estimated via the following GLMs:

$\displaystyle\mathbb{E}[Y|X,A,Z,M]$
$\displaystyle=(1,X,X^{2},A,Z,XZ,M,AM)\theta;\quad\mu_{y}(X,Z)=\mathbb{E}[Y|X,A=a,Z,M=m]$
$\displaystyle\mathbb{E}[U|X,A]$
$\displaystyle=(1,X,X^{2},X^{3},A,XA)\eta;\quad\nu_{y}(X)=\mathbb{E}[U|X,A=a]$
$\displaystyle\pi_{a}(X)$ $\displaystyle=\textup{logit}^{-1}[(1,|X|)\alpha]$
$\displaystyle\pi_{m}(X,Z)$
$\displaystyle=\textup{logit}^{-1}[(1,X,X^{2},A,Z,XA,XZ)\gamma].$

Here $U$ is the outcome variable used to fit the model for $\nu_{y}(X)$, as
shown in equations (7-9). To demonstrate the multiple robustness of the
proposed estimators, we also fit a misspecified model for each of the nuisance
functions:

$\displaystyle\mathbb{E}[Y|X,A,Z,M]$
$\displaystyle=(1,X,A,Z,M)\tilde{\theta};\quad\mu_{y}(X,Z)=\mathbb{E}[Y|X,A=a,Z,M=m]$
$\displaystyle\mathbb{E}[U|X,A]$
$\displaystyle=(1,X,A)\tilde{\eta};\quad\nu_{y}(X)=\mathbb{E}[U|X,A=a]$
$\displaystyle\pi_{a}(X)$
$\displaystyle=\textup{logit}^{-1}[(1,X)\tilde{\alpha}]$
$\displaystyle\pi_{m}(X,Z)$
$\displaystyle=\textup{logit}^{-1}[(1,X)\tilde{\gamma}].$

Each of the four cases described in Figure 1 reflects a combination of
estimated nuisance functions from these correctly and incorrectly specified
models. For example, for submodel $\mathcal{P}_{1}$, we use correctly
specified models for $\mu_{y}(x,z)$ and $\pi_{a}(x)$ and incorrectly specified
models for $\nu_{y}(x)$ and $\pi_{m}(x,z)$ for all estimators.


11institutetext: Amadeus SAS, 821 Avenue Jack Kilby, 06270 Villeneuve-Loubet,
France 22institutetext: Universite de la Cote d’Azur

# Optimizing Revenue Maximization and Demand Learning in Airline Revenue
Management

Giovanni Gatti Pinheiro 11 2 2 Michael Defoin-Platel 11 Jean-Charles Regin 22

###### Abstract

Correctly estimating how demand respond to prices is fundamental for airlines
willing to optimize their pricing policy. Under some conditions, these
policies, while aiming at maximizing short term revenue, can present too
little price variation which may decrease the overall quality of future demand
forecasting. This problem, known as earning while learning problem, is not
exclusive to airlines, and it has been investigated by academia and industry
in recent years. One of the most promising methods presented in literature
combines the revenue maximization and the demand model quality into one single
objective function. This method has shown great success in simulation studies
and real life benchmarks. Nevertheless, this work needs to be adapted to
certain constraints that arise in the airline revenue management (RM), such as
the need to control the prices of several active flights of a leg
simultaneously. In this paper, we adjust this method to airline RM while
assuming unconstrained capacity. Then, we show that our new algorithm
efficiently performs price experimentation in order to generate more revenue
over long horizons than classical methods that seek to maximize revenue only.

###### Keywords:

Revenue Management Demand Learning Price experimentation Earning while
learning Exploration-exploitation trade-off Active Learning

##  1 Introduction

When optimizing the pricing policy, modern revenue management systems consider
only the revenue-maximizing objective, ignoring the long-term effects on the
future learning of the demand behavior. For example, in certain cases, the
optimized pricing policy may have little price variation making it difficult
to correctly estimate the customer’s price sensitivity, impacting the quality
of future forecasts. The proper estimation of the customer’s price sensitivity
is fundamental to airlines, in fact, researchers have found that a bias of
$\pm 20\%$ in the estimation of the price sensitivity can reduce revenue by up
to $4\%$ [7].

To learn the customer’s price sensitivity, the revenue management system (RMS)
needs to charge different prices that may not be optimal with respect to
revenue maximization (i.e., price experimentation), compromising short term
revenue with the hope that the information gained about the demand behavior
will allow the RMS to collect more revenue in the long run. Thus, we have two
competing objectives and a trade-off must be found. This is known as the
earning while learning (EWL) problem.

The EWL problem is not exclusive to the airline industry. It applies to any
seller that needs to price its products in face of an unknown demand behavior.
This is one of the reasons why this problem attracted the attention from
academia and industry in recent years. In one of the most recent and promising
methods, the authors propose to unify revenue maximization and demand learning
into one single objective function [4]. They find that their method
outperforms several others in literature over simulation studies and different
real-world benchmarks. However, their method is limited to the optimization of
a single product (i.e., a flight), which makes it not directly exploitable by
airlines.

We design an adaptation of their method to the airline revenue management (RM)
environment, which must control the pricing policy of several active flights
simultaneously. We focus on specific characteristics of airline RM, such as
the aspect that the historical booking data used for demand model learning is
collected by many parallel flights, and thus, the optimization of the new
objective must consider the effects across different flights. Moreover, for
the sake of simplicity, we assume the price optimization problem under
unconstrained capacity (which simplifies significantly the price optimization
process), and that a simplified demand model is estimated (which facilitates
interpretation). Applying our method to real-world systems is beyond the scope
of this work, and more research is encouraged.

Our work is organized as follows. In Section 2, we present the simplified
single leg problem and how RMS typically estimates the demand behavior and
optimizes the pricing policy. Then, in Section 3, we present some related work
on the investigation of the earning while learning problem. Next, in Section
4, we demonstrate how to adapt the work [4] to the single leg problem,
allowing their method to be applied to airline RM. Following, in Section 5, we
present the experimental results comparing the adapted method to RMS in a
simulation environment. We finish our work by discussing future work and
conclusions in Section 6.

##  2 Background in airline RM

We consider the single-leg problem, which the airline operates one leg from
point A to point B, and a new flight depart at every time step. We assume that
the flight’s capacity is unconstrained and the time is finite with horizon
$H$. The system is allowed to choose fares from a fenceless fare structure
with $n$ price points, where $f\in\mathcal{F}=\\{f_{0},\dots,f_{n-1}\\}$ is a
price point. The RMS goal is to select prices that maximize the expected
revenue $R(f)=f\cdot d(f)$ for each active flight, where $d(f)$ represents the
expected number of bookings for price point $f$. A demand model extensively
used in the literature [6, 8] is the negative exponential demand model
$d(f;\nu,\phi)=\nu e^{-\phi(f/f_{0}-1)}$, where customers arrive according to
a Poisson process with mean $\nu$ and $\phi$ represents the customer price
sensitivity. The constant $f_{0}$ is the base price for which any arriving
customer purchases with probability 1. Since we assume unconstrained capacity
and a discrete fare structure, the optimal fare can be trivially obtained with
$f_{*}=\operatorname*{arg\,max}_{f}R(f)$.

As the RMS does not have access to the true demand behavior parameters
$(\nu_{*},\phi_{*})$, it needs to estimate them from historical bookings. For
simplicity, we assume that RMS knows the true arrival rate $\nu_{*}$, and it
needs to estimate only the price sensitivity $\phi_{*}$. Typically, the
parameter estimation relies on statistical techniques such as ordinary least
squares or maximum likelihood estimation.

One specificity of the single-leg problem is related to how the historical
bookings are collected. The selling horizon for each flight is divided into
$H$ time steps. For each time step, RMS needs to take a pricing decision for
each of the $H$ active flights (each departing in a different date). Because
the historical database has a fixed size, when new data are appended, the
oldest flight data are removed (first-in, first-out). Every time step, the
system calibrates the parameters of the demand model from historical bookings
and it uses the model to optimize the pricing policy. Given the simplicity of
our demand model (i.e., no seasonality, nor day of week, etc.), we assume that
the historical database contains the booking data collected for the last $H$
sell dates. The historical database layout is illustrated in Fig. 1.

Figure 1: The historical database scheme. Each entry in the database consist
of one offer (i.e., the selected fare) and the number of observed bookings for
that offer. We keep entries for each day to departure and for each flight. At
the end of every day, the booking data collected for that day are immediately
added to the historical data and the oldest record is erased.

To compute the price sensitivity parameter $\phi$ and its uncertainty
$\sigma$, we represent the historical booking data with the help of two
functions. The first one, $o(t,f)\in\mathbb{N}_{0}$ returns the number of
times fare $f$ is offered across all flights at the sell date $t$. Similarly,
the second one, $b(t,f)\in\mathbb{N}_{0}$ returns the number of bookings made
for fare $f$ at the sell date $t$ across all flights at the sell date $t$.
Note that the number of customer arrivals is not directly observable (it is
not possible to distinguish customers not willing to purchase at price $f$,
from no customer arrival at all).

According to [14], for sell date $t$, we can write the log-likelihood for the
negative exponential model as

$\mathcal{LL}(\phi)=\sum_{f\in\mathcal{F}}{\sum_{\tau=t-H}^{t-1}\left[b(\tau,f)\ln\left(d\left(f\right)\right)-o\left(\tau,f\right)d\left(f\right)\right]}$
(1)

then, we can compute the estimated price sensitivity with
$\hat{\phi}=\operatorname*{arg\,max}_{\phi}\mathcal{LL}(\phi)$.

To compute the estimated error $\sigma(f)$, we can use the inverse of the
Fisher information $\sigma^{2}\geq 1/\mathcal{I}(\phi)$, given by

$\displaystyle\mathcal{I}(\phi)$
$\displaystyle=-\mathbb{E}\left[\frac{\partial^{2}}{\partial^{2}\phi}\mathcal{LL}(\phi)\bigg{|}\phi\right]$
$\displaystyle=\sum_{f\in\mathcal{F}}{\sum_{\tau=t-H}^{t-1}o(\tau,f)d(f)\left(\frac{f}{f_{0}}-1\right)^{2}}$
(2)

##  3 Earning while Learning

The common practice in industry consists in estimating the parameters of a
demand model, and then assume the unknown parameters are identical to the most
recent estimates when optimizing prices is in the roots of EWL problem. It has
been demonstrated theoretically that this approach can lead to incomplete
learning, that is, the true demand model parameter estimates never converge to
their true values, independently of how much historical data is collected by
the system. To solve incomplete learning, the system must follow a “semi-
myopic” policy, that accumulates (Fisher) information in an adequate rate, and
that do not deviate “too much” from the greedy policy [10].

The investigation of the trade-off between earning and learning has generated
a stream of studies and ideas [2, 5, 9, 11, 13]. Every proposed method shares
the two principles of accumulating information in appropriate rates without
deviating too much from the greedy policy. However, correctly finding the
optimal balance proved to be a tricky task.

One of the most effective and influential methods is called controlled
variance pricing (CVP) [3]. The idea behind the CVP algorithm is to impose a
constraint to the greedy pricing policy that the selected prices are not too
close from the average of previously selected prices, guarantying sufficient
price dispersion. The forbidden range around to the average of past selected
prices is called a taboo interval. As the amount of the historical data
increases, the taboo interval shrinks. In principle, if the amount of
historical data grows to infinity, the interval shrinks to zero and the system
converges to the optimal price.

More recently, a novel method [4] combines the revenue maximization and the
accumulation of information into a single objective function as

$U(f)=R(f)-\eta\left(\frac{\sigma(f)}{\phi}\right)$ (3)

where $\eta$ is an exponentially decreasing trade-off parameter, and
$\sigma(f)$ represents the future uncertainty of the demand price sensitivity
parameter (how much information $f$ will provide). The system selects the fare
that maximizes the objective function $U(f)$, i.e.
$f_{*}=\operatorname*{arg\,max}_{f}U(f)$. Maximizing objective function is
responsible for finding the prices that accumulate more information without
deviating too much from the greedy policy. The trade-off parameter regulates
the importance the optimization algorithm must give to collecting information
of the demand behavior. This method is proven very effective, outperforming
the CVP in simulated and real data benchmarks. However, this, like many other
studies in EWL, focus on settings where only one product (i.e., one flight) is
being optimized. That is, it does not consider the possibility that the data
from multiple products can be aggregated to perform model calibration such is
the case in the airline industry, where the seats from many flights from the
same origin and destination are being sold for different departure dates in
parallel. In the next section, we modify this algorithm expressed in eq. (3)
to allow the optimization of several parallel flights in the single-leg
problem.

##  4 Adapting EWL unified objective to single leg problem

As described in Section 2, modern RMSs assume that each active flight is
independent of the others, optimizing each flight according to
$f_{*}=\operatorname*{arg\,max}_{f}R(f)$. However, when applying eq. (3) to
the single leg problem, it becomes evident that this “divide to conquer”
strategy is no longer effective because the value that $\sigma$ will assume in
the future depends on the aggregated pricing decisions across every active
flight. Therefore, we need to adapt eq. (3) to accommodate the possibility of
optimizing several parallel flights at once.

To do so, first, we define the pricing policy as a multinomial distribution
$\pi=[\pi_{0},\pi_{1},\dots,\pi_{n-1}]$. Each $\pi_{i}$ component of this
distribution represents the probability of selecting the price point $f_{i}$
for every active flight for sell date $t$. Then, we re-write the expected
revenue $R(f)$ as a function of $\pi$ as

$R(\pi)=H\sum_{i=0}^{n-1}f_{i}\pi_{i}d(f_{i})$ (4)

and the Fisher information can be written in terms of $\pi$ as

$\displaystyle\mathcal{I}(\pi)$
$\displaystyle=-\mathbb{E}\left[\frac{\partial^{2}}{\partial^{2}\phi}\mathcal{LL}(\phi)\bigg{|}\phi,\pi\right]$
$\displaystyle=\sum_{i=0}^{n-1}{\underbrace{\mathbb{E}\left[\sum_{\tau=t+H-1}^{t}o(\tau,f_{i})\bigg{|}f_{i}\sim\pi_{i}\right]}_{\text{expected
historical offers at
time}\>t+1}d(f_{i})\left(\frac{f_{i}}{f_{0}}-1\right)^{2}}$
$\displaystyle=\sum_{i=0}^{n-1}{\left(\underbrace{\sum_{\tau=t+H-1}^{t-1}o(\tau,f_{i})}_{\text{already
observed}}+\underbrace{\mathbb{E}\left[o(t,f_{i})\bigg{|}f_{i}\sim\pi_{i}\right]}_{\text{expected
future offers}}\right)d(f_{i})\left(\frac{f_{i}}{f_{0}}-1\right)^{2}}$
$\displaystyle=\sum_{i=0}^{n-1}{\left(\sum_{\tau=t-H+1}^{t-1}o(\tau,f_{i})+H\pi_{i}\right)d(f_{i})\left(\frac{f_{i}}{f_{0}}-1\right)^{2}}$
(5)

In eq. (4), the term representing the last year of historical offers in eq.
(2), i.e. $\sum_{\tau=t-H}^{t-1}o(\tau,f)$, is replaced by the historical
offers minus the oldest entry (that will be eliminated on the next sell date)
plus the expected number of offers that the system is going to make at the
current sell date for each price point, i.e. $H\pi_{i}$.

Assuming the lower bound for the uncertainty,
$\sigma(\pi)=1/\sqrt{\mathcal{I}(\pi)}$, we can re-write the objective
function as

$U(\pi)=R(\pi)-\eta\frac{1}{\phi\sqrt{\mathcal{I}(\pi)}}$ (6)

We seek a vector $\pi$ that maximizes $U(\pi)$. This can be obtained by
solving the non-linear optimization problem

$\displaystyle\operatorname*{arg\,max}_{\pi}\>U(\pi)$ subject to
$\displaystyle\pi_{0},\pi_{1},\dots,\pi_{n-1}\geq 0$
$\displaystyle\sum_{i=0}^{n-1}\pi_{i}=1$

Once the probability distribution $\pi$ is obtained, the system simply selects
the fare for each active flight independently according to $\pi$.

##  5 Experiments

In this section, we investigate how the performance of our method defined by
the optimization metric defined by eq. (6) compares to RMS under the scenario
where the price sensitivity parameter is unknown to the system and it must be
estimated from historical bookings. To solve the non-linear optimization
problem $\operatorname*{arg\,max}_{\pi}\>U(\pi)$, we find that sequential
least squares programming [12] is effective and robust. We assume that each
flight has infinite capacity (representing, for example, the scenario which
the flight’s capacity is much larger than the expect number of customer
arrivals), and that demand behaves according to the negative exponential model
presented in Section 2. We consider a single fare family with 10 price points
$\mathcal{F}=\\{\$50,\$70,\$90,\dots,\$230\\}$ and that the booking horizon
has $H=22$ days. We perform model calibration at every time step. Recall that,
for simplicity, we assume that only the demand price sensitivity parameter is
estimated from historical bookings, and the arrival rate parameter $\nu_{*}$
is known by the system at all time steps.

From now on, we refer to the price sensitivity parameter $\phi$ as the ratio
of the lowest fare at which the purchase probability is 50% (frat5) [1], which
we denoted as $F_{5}$, and it can be mapped according to
$\phi=\ln(2)/(F_{5}-1)$. The frat5 parameter has the convenience of presenting
a linear relationship to the price that maximizes $R(f)$ under the negative
exponential model assumption. We evaluate both RMS and our method within the
interval $F_{5}^{*}\in[2.1,3.8]$, in which the corresponding revenue-
maximizing price $f_{*}$ falls within almost the entire fare structure (except
by the most extreme prices, $50 and $230, that are never optimal for revenue
maximization). From the system’s perspective, the range of possible values of
the true frat5 parameter is unknown, and its true value could in principle
assume any positive value larger than one $F_{5}^{*}>1$. However, to avoid
unrealistic estimations for the frat5 parameter, we limit the estimates of
this parameter to the range $[1.5,4.3]$. The real systems often have this and
other guardrails. Furthermore, we investigate under a low number of arrivals,
$\nu_{*}=4/22=0.18$, where the scarcity of booking data makes it particularly
hard to estimate the true price sensitivity from historical bookings, making
efficient price experimentation fundamental to the success in the task.

Figure 2: Optimizing the trade-off parameter $\eta$ (99% confidence level).
(Left) The average normalized revenue normalized with respect to the revenue
maximizing policy that knows the true demand price sensitivity at every time
step and the random policy that selects fares from a uniformly random
distribution. We plot a separate estimation for RMS represented by the dashed
horizontal line. (Right) The mean squared error of the price sensitivity
estimation $\mathrm{MSE}=\frac{1}{n}\sum_{i=1}^{n}(\hat{\phi}-\phi_{*})^{2}$.

We now analyze the experimental setup and start discussing the results. We
perform three experiments. In the first experiment, we demonstrate how the
trade-off parameter $\eta$ influences the revenue output and the quality of
the demand model parameter estimation. In the second experiment, we select the
best trade-off parameter found in the first experiment, and we analyze how our
method and RMS compare with respect to revenue generation and quality of
demand price sensitivity estimation for different values of the true frat5
parameter within the evaluation interval. In the third and last experiment, we
seek to understand the strategy developed by our method by presenting a
detailed comparison between our method and RMS policies and parameter
estimation for three distinct values of the true frat5 parameter.

First, we focus on the correct calibration of the trade-off parameter $\eta$.
Because we are interested in the system’s overall performance within the
evaluation interval, for this first experiment, we take the performance metric
to be the average collected revenue for the system’s pricing policy $\pi$ when
the true frat5 parameter is sampled uniformly within the evaluation interval.
We sample 160 values for the exploration rate randomly in the interval
$\eta\in[0,8000]$, and for each one of the sampled exploration rates, we
sample 2560 values of $F_{5}^{*}$ randomly in the evaluation interval and we
simulate each episode as described earlier for $20\cdot H=440$ time steps. In
other words, the performance metric is the average of the revenue obtained by
the 2560 episodes of 440 time steps, each having a different true frat5
parameter randomly selected in the evaluation interval. Mathematically, we
seek to approximate $\mathbb{E}_{F_{5}^{*}\sim[2.1,3.8]}[R(\pi;F_{5}^{*})]$.
In Fig. 2 (left), we show how the system’s performance varies with respect to
the trade-off parameter. When $\eta=0$, the system greedily chooses the
revenue-maximizing price only, behaving like traditional RMS. As the trade-off
parameter increases, the mean squared error (MSE) of the price sensitivity
estimation decreases Fig. 2 (right). The better price sensitivity estimation
for allows the system to improve its revenue performance. After a point,
increasing the trade-off parameter translates to a loss of revenue due to
excessive price experimentation. The best revenue performance was obtained at
$\eta=2167$, with a normalized expected revenue of $78.0\%$, which represents
an absolute improvement of $7.0\%$ when compared to RMS.

For the second experiment, we set the true frat5 parameter to a fixed value of
the evaluation interval, and for each point we compute the average revenue
performance and the price sensitivity MSE of 4000 independent runs when
$\eta=2167$. The comparative results between RMS and our method can be found
in Fig. 3 (left). We observe that our method performs better than traditional
RMS in the entirety of the evaluation interval, with the largest advantage for
intermediate values of frat5 $F_{5}^{*}\in[2.2,3.0]$. In practice, this range
of the frat5 represents a significant part of the evaluation interval that we
are most interested in, because it overlaps with the center of the fare
structure, where the RMSs are usually calibrated to operate. In Fig. 3
(right), we plot the MSE over the estimation of the price sensitivity. RMS
presents larger errors for lower values of the frat5, decreasing as the frat5
increases. Such effect is due to the fact that the true arrival rate parameter
$\nu_{*}$ is known by the system. As the frat5 increases, so does the optimal
fare. Further the selected price is from the base fare $f_{0}$, the higher is
the amount of information collected by that price, which can be verified with
eq. (2). Therefore, the higher is the value of the frat5 parameter, the easier
the estimation of the price sensitivity becomes, and naturally, there is less
value in price experimentation. Our method, on the other hand, displays much
less average MSE on the estimation of the customer’s price sensitivity. The
effect is more important for lower values of frat5 parameter, where the
system’s policy approaches the base fare $f_{0}$, reducing the amount of
information obtained by the interaction with the demand. Again, as the frat5
increases, the revenue maximizing price matches the information maximizing
price. Consequently, we observe that both methods produce similar results.

Figure 3: The comparison between the average performance of RMS and our method
(99% confidence level). (Left) The normalized revenue obtained by both
methods. (Right) The mean squared error of the price sensitivity estimation.

Lastly, we perform a detailed analysis of the policy and estimated price
sensitivity according to eq. (1) for both RMS and our method for three
distinct frat5. For the intermediate value of the frat5 parameter
$F_{5}^{*}=2.56$ Fig. 4(a), we see on the left figures a comparison between
the two policies. In the upper left chart, we can see that RMS tends to price
lower fares ($50, $70, $90) rather frequently (38% of the time), which is a
mistake from the perspective of revenue maximization (revenue-maximizing fare
$f_{*}=\$110$), but also a mistake from the perspective of information
maximization, because the lower fares provide less information about the
demand price sensitivity than the higher fares (see eq. (2)). Instead, our
method shifts pricing towards higher fares that provide more information about
the customer’s price sensitivity. For the figures on the right, we plot the
estimated price sensitivity according to each method. Generally, we see a
better estimation of demand price sensitivity for our method (lower right
chart), with a narrower distribution around the true value represented by the
vertical dashed line.

For the lower frat5 value, $F_{5}^{*}=2.1$ Fig. 4(b), we see in the upper left
chart that RMS tends to price the base fare $f_{0}=\$50$ very frequently (40%
of the time), which is neither the revenue-maximizing price (optimal fare
$f_{*}=\$70$), nor appropriate to learn anything about the demand price
sensitivity because the base fare provides no information (according to eq.
(2)). On the right upper chart, we observe that the estimated frat5 parameter
is often at the minimum accepted value $F_{5}=1.5$, which explains RMS
preference for the base fare $f_{0}$. However, the base fare is near optimal
from the revenue-maximizing perspective because of its proximity to the
optimal fare. This illustrates the reasons behind RMS surprisingly good
revenue performance at low values of the true frat5 parameter that we see in
Fig. 3. In other words, RMS is being “helped” by the problem settings. Our
method, on the other hand, shifts the policy towards higher fares (left bottom
chart), which in turn improves the estimation of the frat5 parameter (right
bottom chart). As a result, our method provides more revenue on average than
the RMS policy while yielding a better estimation of the demand price
sensitivity.

For the higher value of frat5 $F_{5}^{*}=3.7$ in Fig. 4(c), we observe a close
behavior between RMS and our method. Because the optimal price $f_{*}=\$190$
provides important information about the customer price sensitivity, there is
little gain in experimenting with other prices. In this case, most revenue is
lost when the system underestimates the frat5 parameter, selecting lower fares
which are both suboptimal for revenue and information maximization.

(a) $F_{5}^{*}=2.59$.

(b) $F_{5}^{*}=2.1$.

(c) $F_{5}^{*}=3.8$.

Figure 4: Comparison of the policy and the price sensitivity estimation
between RMS an our method. For convenience, we represent the revenue-
maximizing policy when the true price sensitivity is known by the system at
all times.

##  6 Conclusions and future work

Correctly estimating the demand price sensitivity is important for RMSs.
However, due to little price variation, the RMSs may suffer from incomplete
learning of the demand behavior, and as a consequence, produce suboptimal
policies.

Inspired by the work developed in [4], we present a novel method for jointly
optimizing the revenue maximization and model learning under the monopolistic
single leg airline RM problem. We show the effectiveness of our method under
an unconstrained capacity assumption when estimating a single model parameter
(price sensitivity). Our method improves over RMS, which maximizes revenue
only, demonstrating that, in some cases, there is value in price
experimentation and in the control of model uncertainty.

The value of this research is that it illustrates that airline RMSs may be
improved by considering the uncertainty of the demand model parameter
estimates during the price optimization. However, because of the unconstrained
capacity assumption, this method cannot be applied directly to real-world
systems, and further investigation must be done.

Applying our method to real-life airline scenario is challenging for several
reasons. First, the unconstrained capacity assumption is not realistic. In
practice, airlines have fixed sized aircraft. To maximize revenue under
capacity constraints, airlines often use dynamic programming [8], which is an
optimization technique based on the Bellman optimality equation and value
functions. Thus, integrating the method here presented with dynamic
programming is essential to its applicability in real world problems.

Furthermore, we demonstrate the effectiveness of our method when controlling
the error of a single model parameter. However, in real world, demand models
may have up to 30 parameters to be estimated and optimized. Obtaining an
estimation of the error for each parameter may be not a trivial task. Thus, it
will be important to investigate methods for estimating the uncertainty of
each model parameter.

Lastly, we calibrated the trade-off parameter $\eta$ by brute forcing values
and then selecting the value that works best according to our problem
definition. In real world systems, this is very unlikely to be possible.
Therefore, we believe that it is important to develop methods for choosing the
trade-off parameter that can work under production constraints.

## Acknowledgments

Thanks to Thomas Fiig, Micheal Wittman, Alexender Papen and Tianshu Yang for
the support and useful comments. We would like to also thanks the financial
support from Amadeus SAS and the french Industrial Agreements for Training
through Research (CIFRE) program.

## References

  * [1] Belobaba, P.P., Hopperstad, C.: Algorithms for revenue management in unrestricted fare markets. In: Meeting of the INFORMS Section on Revenue Management, Massachusetts Institute of Technology, Cambridge, MA (2004)
  * [2] Besbes, O., Zeevi, A.: Dynamic pricing without knowing the demand function: Risk bounds and near-optimal algorithms. Operations Research 57(6), 1407–1420 (2009)
  * [3] den Boer, A.V., Zwart, B.: Simultaneously learning and optimizing using controlled variance pricing. Management science 60(3), 770–783 (2014)
  * [4] Elreedy, D., Atiya, A.F., Shaheen, S.I.: Novel pricing strategies for revenue maximization and demand learning using an exploration–exploitation framework. Soft Computing 25(17), 11711–11733 (2021)
  * [5] Ferreira, K.J., Simchi-Levi, D., Wang, H.: Online network revenue management using thompson sampling. Operations research 66(6), 1586–1602 (2018)
  * [6] Fiig, T., Isler, K., Hopperstad, C., Belobaba, P.: Optimization of mixed fare structures: Theory and applications. Journal of Revenue and Pricing Management 9(1), 152–170 (2010)
  * [7] Fiig, T., Weatherford, L.R., Wittman, M.D.: Can demand forecast accuracy be linked to airline revenue? Journal of Revenue and Pricing Management 18(4), 291–305 (2019)
  * [8] Gallego, G., Van Ryzin, G.: Optimal dynamic pricing of inventories with stochastic demand over finite horizons. Management science 40(8), 999–1020 (1994)
  * [9] Harrison, J.M., Keskin, N.B., Zeevi, A.: Bayesian dynamic pricing policies: Learning and earning under a binary prior distribution. Management Science 58(3), 570–586 (2012)
  * [10] Keskin, N.B., Zeevi, A.: Dynamic pricing with an unknown demand model: Asymptotically optimal semi-myopic policies. Operations research 62(5), 1142–1167 (2014)
  * [11] Keskin, N.B., Zeevi, A.: Chasing demand: Learning and earning in a changing environment. Mathematics of Operations Research 42(2), 277–307 (2017)
  * [12] Kraft, D., et al.: A software package for sequential quadratic programming (1988)
  * [13] Lobo, M.S., Boyd, S.: Pricing and learning with uncertain demand. In: INFORMS Revenue Management Conference. Citeseer (2003)
  * [14] Newman, J.P., Ferguson, M.E., Garrow, L.A., Jacobs, T.L.: Estimation of choice-based models using sales data from a single firm. Manufacturing & Service Operations Management 16(2), 184–197 (2014)