# Quadruply robust estimation of marginal structural models in observational
studies subject to covariate-driven observations

Janie Coulombe  Department of Mathematics and Statistics, Université de
Montréal, Montreal, Quebec, H3T 1J4 Canada  Shu Yang  Department of
Statistics, North Carolina State University, 2311 Stinson Drive, 5109 SAS
Hall, Raleigh, United States

###### Abstract

Electronic health records and other sources of observational data are
increasingly used for drawing causal inferences. The estimation of a causal
effect using these data not meant for research purposes is subject to
confounding and irregular covariate-driven observation times affecting the
inference. A doubly-weighted estimator accounting for these features has
previously been proposed that relies on the correct specification of two
nuisance models used for the weights. In this work, we propose a novel
consistent quadruply robust estimator and demonstrate analytically and in
large simulation studies that it is more flexible and more efficient than its
only proposed alternative. It is further applied to data from the Add Health
study in the United States to estimate the causal effect of therapy
counselling on alcohol consumption in American adolescents.

##  1 Introduction

The study of causes and effects is an important scientific topic and an
essential component of learning healthcare systems aimed at improving public
health (see e.g., Krumholz (2014); Dahabreh & Kent (2014)). When prescribing a
treatment, physicians need to know about the causal effect of that treatment
on the outcome to improve, instead of the mere association (or statistical
dependence) between the treatment and outcome. Statistical methods with good
properties must, therefore, be developed to estimate causal effects
consistently with the least bias and variance possible. This manuscript
proposes a novel, consistent and efficient estimator for the marginal causal
effect of a treatment (exposure) on a longitudinal outcome using observational
data.

Randomized controlled trials are the gold standard for drawing causal
inferences. The randomization of participants to treatments ensures a balance
in patient characteristics between treatment groups at study entry, allowing a
fair comparison of clinical outcomes across treatment groups (Rubin, 2008).
Often, randomized controlled trials also have clear protocols for the timing
of patients’ visits (i.e., observation times) at which patient health status
is measured.

It is not always possible to conduct a randomized controlled trial designed
specifically for answering a causal question, such that researchers often turn
to observational data (Black, 1996). In this work, we focus on the particular
features often met in observational data from electronic health records.
Electronic health records data may contain information on patients’
demographic variables, diagnoses, health system usage, drug dispensations, and
outcomes, making them highly useful to infer drug safety and effectiveness.

While electronic health records are increasingly available for analysis, they
are not meant for research purposes. The treatments measured in electronic
health records are not randomized to patients. Their observation rather
depends on real-world prescription mechanisms, which depend on patients’
characteristics. This leads to spurious associations in the data called
confounding (Greenland & Morgenstern, 2001). These data are also measured
irregularly across patients. They are not collected under a controlled scheme.
Instead, each patient follows their pattern in how they access care, also
likely to depend on their characteristics (Lokku et al., 2020). For instance,
sicker patients tend to interact with the healthcare system more often than
healthier patients. Statistically, this creates a long-term dependence
structure between the outcome and the visit processes that can result in
biased estimators of causal or associational parameters (see e.g., Lin & Ying
(2001); Lipsitz et al. (2002); Farewell (2010); Pullenayegum & Lim (2016) and
most recently Coulombe et al. (2021a, b, 2022); Yang (2022) in a context of
causal inference). When aiming for a causal effect, this bias can be due to
confounding by the visit process or, if the visit indicators act as colliders
(i.e. are themselves affected by the treatment prescribed and the study
outcome), to collider-stratification bias described in Greenland (2003).

Under a set of causal assumptions, the causal marginal effect of a binary
treatment on a longitudinal continuous outcome can generally be inferred by
estimating the parameters of a marginal structural model fitted on the data
from a pseudo-population that is free of confounding and other types of
spurious associations, such as collider-stratification bias (Robins et al.,
2000). Previous work has tackled this problem in settings with covariate-
driven observation times and confounding, leading to the Flexible Inverse
Probability of Treatment and Monitoring weighted (FIPTM) estimator (Coulombe
et al., 2021b, 2022). However, that method suffers from two important
problems. First, it relies on the correct specification of the treatment and
outcome observation models as a function of patient characteristics. This
implies both a correct specification of the variables to include in each model
and their functional forms. When one or both models are not correctly
specified, the FIPTM can be biased. Secondly, that estimator can be variable
due to its inverse weights and could be made more efficient by using the
geometry of influence functions and by augmenting the estimating equations
(Robins et al., 1994; Tsiatis, 2006).

To address the two issues raised above, we propose the first quadruply robust
estimator for the causal marginal effect of a binary treatment on a
longitudinal, continuous outcome, that accounts for confounding and irregular
covariate-driven observation times of the outcome simultaneously.

##  2 Methods

###  2.1 Notation

Our interest is in the causal marginal effect of a binary treatment or
exposure (a choice between two possible options, say treatment with aspirin
versus treatment with the absence of aspirin) on a longitudinal continuous
outcome that is measured repeatedly. That effect could consist of a “contrast”
causal effect between two different active drugs, or of the causal effect of
an active drug when compared with a placebo.

We assume working with a random sample of size $n$ from a larger population
and denote by $i$ the patient index and by $t\in[0,\tau]$ the time, with
$\tau$ a maximum follow-up time in the cohort. Let $A_{i}(t)$ represent the
binary treatment taking values in $\left\\{0,1\right\\}$ and $Y_{i}(t)$ be the
continuous, longitudinal study outcome for patient $i$ at time $t$, that is
only observed at certain points in time.

The type of data-generating mechanism we focus on is presented in the left
panel of Fig. 1. The proposed approach is not specific to that data-generating
mechanism and the novel estimator could be tailored according to various data-
generating mechanisms. In Fig. 1, variables $K_{i}(t)$ are confounders that
create (open) backdoor paths from the treatment $A_{i}(t)$ to the outcome
$Y_{i}(t)$. Associations that are not due to the causal effect of the
treatment on the outcome can pass through them (Pearl, 2009). The set
$M_{i}(t)$ denotes a set of mediators for the treatment effect on the outcome.
The set $P_{i}(t)$ contains pure predictors of the outcome that potentially
also affect the observation of a patient outcome $Y_{i}(t)$.

We use a counting process to model observation times. Only the outcome process
is assumed to be measured sporadically (one may think of weight or systolic
blood pressure measured irregularly according to an observation scheme
depending on patient characteristics such as a change in medication or smoking
status). All the other variables necessary to the estimation of the marginal
causal effect of treatment are assumed to be available at all times during
follow-up, for each patient. While this may seem like a strong assumption, the
drugs and comorbidities are typically recorded in electronic health records
data anytime there is a new diagnosis or a new prescription is made or
dispensed. Let $N_{i}(t)$ be the total number of observation times of the
outcome $Y_{i}(t)$ between times 0 and $t$ for individual $i$. The indicator
$\rm{d}N_{i}(t)$ equals to 1 when there is a jump in the process at time $t$,
i.e., the observation of the outcome $Y_{i}(t)$, and 0 otherwise. We denote by
$T_{i,1},...,T_{i,Q_{i}}$ the observation times of the outcome, with $Q_{i}$
the total number of observation times of individual $i$. The set $V_{i}(t)$
includes all the variables causing observation times (i.e., causing
$\rm{d}N_{i}(t)$ in Fig. 1) and also including all the confounders of the
treatment-outcome relationship. We must include $K_{i}(t)$ in the set
$V_{i}(t)$ for the proposed methodology to be consistent. In Fig. 1, the set
$V_{i}(t)$ contains the treatment, the mediator(s) of the treatment effect on
the outcome, the confounder(s), and the pure predictor(s).

$A_{i}(t)$$M_{i}(t)$$Y_{i}(t)$$\rm{d}N_{i}(t)$$K_{i}(t)$$P_{i}(t)$ (a) Assumed
causal diagram at time $t$ in patient $i$, postulated to be common across all
patients

Indicator ofcounselling($t$)Depressivemood$(t)$Alcohol
consumption$(t)$$\rm{d}N(t)$Confounders$(t)$ (b) Causal diagram for the Add
Health study data at time $t$ common across all adolescents

Figure 1: Causal diagram illustrating the assumed data generating mechanism

Patients are allowed to have different follow-up times $C_{i}$, assumed to be
non-informative conditional on the outcome model design matrix, an assumption
denoted by $Y_{i}(t)\perp C_{i}\mid A_{i}(t),V_{i}(t)$. Finally,
$\xi_{i}(t)=\mathbf{1}\\{C_{i}>t\\}$ is an indicator that patient $i$ is still
in the study at time $t$, meaning they are not lost to follow-up.

###  2.2 Causal Estimand

The potential outcome framework (Neyman, 1923; Rubin, 1976) is used to define
our estimand. Denote by $Y^{1}_{i}(t)$ the potential outcome of individual $i$
at time $t$ if they received treatment option $1$, and $Y^{0}_{i}(t)$ for
treatment $0$. The causal marginal effect of a binary treatment on a
continuous outcome is defined as
$\beta_{1}=E\left[Y^{1}_{i}(t)-Y^{0}_{i}(t)\right]$. The treatment and the
outcome are allowed to vary in time, but our interest lies in a cross-
sectional effect $\beta_{1}$ that does not vary in time.

Suppose a certain time discretization for which there can be only one jump in
the counting process $N(t)$. For instance, suppose visits at the doctor’s
office can occur daily such that the time granularity is the day. If we had
access to all potential outcomes under both treatments and at each time $t$
(for the time granularity chosen above) for a random sample of participants of
size $n$, then we could estimate $\beta_{1}$ straightforwardly using sample
means. This is impossible in practice and sometimes referred to as the
fundamental problem of causal inference (Holland, 1986).

On the other hand, by conducting a randomized controlled trial and randomly
allocating patients to one of the two treatment options, and observing
patients at prespecified visit times, there is no reason why patients
allocated to, e.g., treatment $1$, would differ than the others before
receiving the treatment. One could then use the marginal structural model to
estimate $\beta_{1}$:

$E\left[Y^{a}_{i}(t)\right]=E\left[Y^{a}_{i}(t)\mid
A_{i}(t)=a\right]=\beta_{0}+\beta_{1}a.$ (1)

Under this model, $E\left[Y^{1}_{i}(t)\right]=\beta_{0}+\beta_{1}$ and
$E\left[Y^{0}_{i}(t)\right]=\beta_{0}$, and the causal contrast of interest is
the parameter $\beta_{1}$. An outcome consistency assumption (discussed in
$\mathsection$2.3) is required to estimate $\beta_{1}$ with the marginal
structural model above, but it generally does not require any adjustment for
confounding under random treatment allocation.

In observational data from electronic health records, unfortunately, data are
affected by a confounding mechanism, whereby we observe the potential outcomes
$Y^{1}_{i}(t)$ in those who had greater chances of being treated with
$A_{i}(t)=1$, and the potential outcomes $Y^{0}_{i}(t)$ in those who had
greater chances of being treated with $A_{i}(t)=0$ (as a consequence,
$E\left[Y^{a}_{i}(t)\mid A_{i}(t)=a\right]\neq E\left[Y^{a}_{i}(t)\right]$ in
general). In addition to that confounding, the potential outcomes for an
individual $i$ are only observed at times $t$ when $\rm{d}N_{i}(t)=1$, which
may also depend on patient characteristics. Therefore, we do not have access
to all potential outcomes and require a set of causal assumptions to equate
the estimand, which depends on $Y_{i}^{a}(t)$, to equations depending on
$Y_{i}(t)$.

###  2.3 Causal Assumptions

Five causal assumptions are required for consistent estimation of the causal
marginal effect of treatment (Table 1). Other (non-causal) assumptions on the
correct specification of the nuisance models are also required and discussed
in $\mathsection$2.4.

First, it requires outcome consistency, which implies that
$Y_{i}(t)=Y^{1}_{i}(t)$ if patient $i$ received treatment 1 at time $t$ and
$Y_{i}(t)=Y^{0}_{i}(t)$ otherwise and means that the treatment is well defined
such that the observed outcome effectively corresponds to one of the two
potential outcomes. An example of consistency violation is one in which a
potential outcome is defined as the outcome under the treatment physical
exercise 5 times a week and in which the actual treatment received corresponds
to having done exercise 3 times a week. The outcome observed under that
regimen would not correspond to the defined potential outcome.

We assume positivity of treatment, meaning that anyone should have a chance of
receiving any of the two treatment options, and positivity of observation,
such that they had a chance to have their outcome observed at any time given
their characteristics. For instance, scenarios in which some patient
characteristics used as predictors in the treatment model are only represented
in one of the two treatment groups would violate the positivity of treatment
assumption, and similarly for the observed and non-observed groups.

Finally, we assume conditional exchangeability, with includes the assumptions
of 1) no unmeasured confounder, i.e., all confounders $K_{i}(t)$ of the
relationship between the treatment and the outcome are available in the
analysis; and 2) conditional independence of the observation indicators, i.e.,
adjusting for $V_{i}(t)$ makes the observation indicator independent of other
variables in the analysis. In our setting, $K_{i}(t)\subset V_{i}(t)$ and
exchangeability is conditional on $V_{i}(t)$.

These five assumptions allow the use of so-called G-methods (Naimi et al.,
2017; Robins & Hernán, 2008, Ch. 23, p. 553) for consistent estimation of
causal effects. The methods discussed next can be thought of as being part of
a larger G-methods framework.

Table 1: Causal assumptions required for the proposed estimator to be
consistent

Assumption | Definition
---|---
Outcome consistency | $Y_{i}(t)=A_{i}(t)Y_{i}^{1}(t)+\left\\{1-A_{i}(t)\right\\}Y_{i}^{0}(t)$
Positivity of treatment | $0<pr\\{A_{i}(t)\mid K_{i}(t)\\}<1$
Positivity of observation | $0<E[\rm{d}N_{i}(t)\mid V_{i}(t)]<1$
No unmeasured confounder | $\left\\{Y_{i}^{0}(t),Y_{i}^{1}(t)\right\\}\perp A_{i}(t)\mid K_{i}(t)$
Conditional exchangeability |  $\left\\{Y_{i}^{0}(t),Y_{i}^{1}(t)\right\\}\perp A_{i}(t)\mid K_{i}(t)$ and $\\{A_{i}(t),Y^{0}_{i}(t),Y^{1}_{i}(t)\\}\bot\rm{d}N_{i}(t)\mid V_{i}(t)$

###  2.4 Novel Estimator

The conditional exchangeability can be recovered by breaking the spurious
associations due to the treatment and observation mechanisms via inverse
weights (marginal approach), by conditioning on the sets
$A_{i}(t),K_{i}(t),V_{i}(t)$ in a regression model for the outcome and using
methods such as g-computation (Robins, 1986) (standardization approach), or by
using both approaches simultaneously to make it more robust to models
misspecification, which is our proposal.

Using the marginal approach corresponds to using the FIPTM estimator proposed
in Coulombe et al. (2021b). It consists of a doubly-weighted least squares
estimator that incorporates inverse probability of treatment weights (Horvitz
& Thompson, 1952; Rosenbaum & Rubin, 1983; Rosenbaum, 1987) and inverse
intensity of visit weights (Lin et al., 2004). The inverse probability of
treatment weights are functions of the confounders $K_{i}(t)$ and the IIV
weights are functions of the visit predictors $V_{i}(t)$. The estimator is
consistent for $\beta_{1}$ when both weights are correctly specified. A
parametric model can be used to model the treatment, and the inverse
probability of treatment weights be obtained as follows:

$\displaystyle\mathbf{1}\\{A_{i}(t)=a\\}/pr\\{A_{i}(t)=a\mid K_{i}(t);\psi\\}$
(2)

where $pr\\{A_{i}(t)=1\mid K_{i}(t);\psi\\}$, the propensity score, is the
probability of receiving the treatment $1$ as a function of predictors
$K_{i}(t)$ and parameters $\psi$ (Rosenbaum & Rubin, 1983). A logistic
regression can be used to compute an estimated propensity score. The inverse
intensity of visit weights, on the other hand, can be obtained from the
proportional rate model:

$\displaystyle E[\rm{d}N_{i}(t)\mid
V_{i}(t);\gamma]=\xi_{i}(t)\exp\\{\gamma^{{\mathrm{\scriptscriptstyle
T}}}V_{i}(t)\\}\lambda_{0}(t)\rm{d}t.$ (3)

The baseline rate of observation $\lambda_{0}(t)$ in (3) consists of the visit
rate when all variables $V_{i}(t)$ are set to their reference level. With the
FIPTM estimator, the baseline rate can be dropped from the inverse intensity
of visit weights without affecting the causal marginal effect of treatment
estimate since removing it would still make the weights in (4) proportional to
the intensity of being observed as a function of $V_{i}(t)$. The inverse
intensity of visit weights can also be stabilized, in which case the baseline
rate cancels automatically in the weights and need not be estimated (Buzkova &
Lumley, 2009). This leads to the following intensity of visit weights (we take
the inverse in the estimating equations), from which $\gamma$ parameters can
be estimated using the Andersen & Gill (1982) model:

$\displaystyle E[\rm{d}N_{i}(t)\mid
V_{i}(t);\gamma]=\exp\\{\gamma^{{\mathrm{\scriptscriptstyle T}}}V_{i}(t)\\}.$
(4)

Then, the FIPTM estimator solves the following equation:

$\displaystyle
E_{n}\left[\int_{0}^{\tau}\frac{\frac{\mathbf{1}\\{A_{i}(t)=a\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}Y_{i}(t)-\zeta_{i}(t;\beta_{a})}{E\\{\rm{d}N_{i}(t)\mid
V_{i}(t);\hat{\gamma}\\}}\rm{d}N_{i}(t)\right]=0,$ (5)

where $\zeta_{i}(t;\beta_{a})=\beta_{0}+\beta_{1}a$ is the structural model
from equation (1) and $E_{n}$ stands for the empirical mean. A drawback of the
doubly-weighted estimator is that it requires both the treatment and the
observation models to be correctly specified. This is not easy in practice.

We propose the augmented AAIIW estimator (which acronym stands for doubly
augmented and doubly inverse weighted) that is more flexible and allows two
out of four different models (see Table 2) to be misspecified while the
estimator remains consistent. The estimator uses the theory introduced in
Robins et al. (1994), very well laid out in Funk et al. (2011), Tsiatis &
Davidian (2007) and Cao et al. (2009), and further related to model-assisted
estimation from the survey sampling field (see e.g., the discussion in Robins
& Rotnitzky (1998)). Using the model-assisted estimation approach to justify
the construction of the novel estimator may be more intuitive to the reader
than using the semiparametric theory of influence functions (see e.g., Jiang
et al. (2022) for related discussions), so we briefly discuss that framework.
Grosso modo, the estimating equations of the FIPTM estimator are to be
transformed twice following the model-assisted estimation approach. Based on
their work in Robins et al. (1994), Robins & Rotnitzky (1998) introduce the
equation

$\displaystyle
N\hat{T}_{diff}(\mu)=\sum_{i=1}^{N}A_{i}\mu(X_{i})+\sum_{i:A_{i}=a}A_{i}\left\\{Y_{i}-\mu(X_{i})\right\\}/\pi_{i}$

and connect it to model-assisted estimation, where $\hat{T}_{diff}(\mu)$ is a
designed-based standard difference estimator for the parameter $E[A_{i}Y_{i}]$
of interest, and $\mu(X)$ is a function of $X$. Theorem 1 in Robins &
Rotnitzky (1998) implies that the class of such estimators
$\hat{T}_{diff}(\mu)$ contains all the semiparametric estimators and that the
asymptotic variance of the estimator evaluated at $\mu_{eff}(x)=E[Y_{i}\mid
X_{i}=x]$ in their notation leads to the smallest variance possible. We use
that construction twice, starting first with
$A_{i}=\mathbf{1}\\{A_{i}(t)=a\\}$, $\pi_{i}=pr\\{A_{i}(t)=a\mid K_{i}(t)\\}$
and $``\mu(X_{i})"=E\left[Y_{i}(t)\mid A_{i}(t)=a,K_{i}(t)\right]$ using their
notation. Once this projection is obtained, we use the strategy a second time,
with the “$Y_{i}$” now corresponding to the previous expression obtained, and
with $``A_{i}"=\rm{d}N_{i}(t)$, $``\pi_{i}"=pr\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t)\\}$ and a novel $``\mu(X_{i})"$ function corresponding to the
expectation of the previous term before augmentation (in equation (6) below,
this expectation corresponds to $E[\eta_{i}(t)\mid
A_{i}(t)=a,K_{i}(t),V_{i}(t)]$). By construction, the novel proposed estimator
is the most efficient among its class of semiparametric estimators, which also
includes the FIPTM estimator.

The correspondence between the model-based estimation approach and the
geometry of influence functions is in that the novel estimator’s influence
function corresponds to sequential projections of the FIPTM’s influence
function onto spaces orthogonal to the residuals from the weight models and
orthogonal to projections due to the outcome mean models onto the spaces of
$K_{i}(t)$ and $V_{i}(t)$. The novel estimator is obtained by solving the
following estimating equations, which are augmented versions of the equations
in 5:

$\displaystyle
E_{n}\left[\int_{0}^{\tau}\frac{\eta_{i}(t)}{E\\{\rm{d}N_{i}(t)\mid
V_{i}(t);\hat{\gamma}\\}}\rm{d}N_{i}(t)\right]-E_{n}\left[\int_{0}^{\tau}\frac{\rm{d}M_{i}(t)E\\{\eta_{i}(t)\mid
A_{i}(t)=a,K_{i}(t),V_{i}(t)\\}}{E\\{\rm{d}N_{i}(t)\mid
V_{i}(t);\hat{\gamma}\\}}\right]$ $\displaystyle=0,$ (6)

where the nuisance terms are estimated using parametric models, with

$\displaystyle\eta_{i}(t)=$
$\displaystyle\frac{\mathbf{1}\\{A_{i}(t)=a\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}Y_{i}(t){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-\frac{\mathbf{1}\\{A_{i}(t)=a\\}-pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}\mu_{a}\\{K_{i}(t);\hat{\alpha}_{K}\\}}-\zeta_{i}(t;\beta_{a}),$

with
$\rm{d}M_{i}(t)=\rm{d}N_{i}(t)-\xi_{i}(t)\exp\\{\hat{\gamma}^{{\mathrm{\scriptscriptstyle
T}}}V_{i}(t)\\}\hat{\lambda}_{0}(t)\rm{d}t$ the martingale residual for the
observation process. The conditional outcome mean models in the augmentated
terms are $\mu_{a}\left\\{K_{i}(t);\alpha_{K}\right\\}=E[Y_{i}(t)\mid
A_{i}(t)=a,K_{i}(t);\alpha_{K}]$ and
$\mu_{a}\left\\{V_{i}(t);\alpha_{V}\right\\}=E[Y_{i}(t)\mid
A_{i}(t)=a,V_{i}(t);\alpha_{V}]$. The latter model arises when taking the
expectation $E\left[Y_{i}(t)\mid A_{i}(t)=a,V_{i}(t)\right]$ in the term
$E\\{\eta_{i}(t)\mid A_{i}(t)=a,K_{i}(t),V_{i}(t)\\}$ in equation 6 (more
details are found in Supplementary Material A and B discussed later). For the
novel estimator, the baseline rate $\lambda_{0}(t)$ in (3) must be estimated
before calculating the inverse intensity of visit weights. The inverse
intensity of visit weights in the equations for the AAIIW are the inverse of
$E[\rm{d}N_{i}(t)\mid
V_{i}(t);\hat{\gamma}]=\hat{\lambda}_{0}(t)\exp\\{\hat{\gamma}^{{\mathrm{\scriptscriptstyle
T}}}V_{i}(t)\\}$. We use the following Breslow’s estimator:

$\displaystyle\hat{\lambda}_{0}(t)=\frac{\sum_{i=1}^{n}\mathbf{1}\\{\rm{d}N_{i}(t)=1\\}}{\sum_{i=1}^{n}\mathbf{1}\\{\rm{d}N_{i}(t)=1\\}\exp\left\\{\hat{\gamma}^{T}V_{i}(t)\right\\}},$

(Cox, 1972). Table 2 shows the possible combinations of correctly specified
models leading to a consistent AAIIW estimator. For the estimator to be
unbiased, at least one of the two models related to confounders (either the
treatment or the outcome mean model conditional on the confounders) and at
least one of the two models related to the observation predictors (either the
observation or the outcome mean model conditional on the observation
predictors) must be correctly specified. Correct specification for a nuisance
model requires that the corresponding data-generating mechanism can be
modelled parametrically and that there exists a true set of parameters leading
to the actual data-generating mechanism that we can estimate consistently (we
denote the true sets by $\psi_{0}$, $\gamma_{0}$, $\alpha_{K0}$ and
$\alpha_{V0}$ for the treatment, observation, and two conditional mean outcome
models, respectively). For instance, for the treatment model, it means that
the data generating mechanism $pr\\{A_{i}(t)=a\mid
K_{i}(t)\\}=pr\\{A_{i}(t)=a\mid K_{i}(t);\psi_{0}\\}$, that we can model this
data generating mechanism using the correct functional formats for covariates
$K_{i}(t)$ in the model, and that estimators $\hat{\psi}$ converge in
probability to the true parameters $\psi_{0}$.

Naturally, the AAIIW estimator is unbiased under specific combinations of two
correctly specified models (Table 2, see the proof of multiple robustness in
Supplementary Material A) but also when three out of the four models in Table
2 are correctly specified, or when all models are correctly specified. Thus,
using the same nuisance models as those used in the FIPTM estimator, if they
are correctly specified, leads to an unbiased estimator, but the advantage of
the AAIIW estimator is that it also has several other opportunities to be
unbiased.

Table 2: Quadruple robustness of AAIIW: AAIIW is consistent under Scenarios
(a)–(d): $\checked$ means correctly specified and $X$ means no requirement

Scenario | $pr\\{A_{i}(t)=a\mid K_{i}(t);\psi\\}$ | $\mu_{a}\\{K_{i}(t);\alpha_{K}\\}$ | $E\left[\rm{d}N_{i}(t)=1\mid V_{i}(t);\gamma\right]$ | $\mu_{a}\\{V_{i}(t);\alpha_{V}\\}$
---|---|---|---|---
(a) | $\checked$ | $X$ | $\checked$ | $X$
(b) | $X$ | $\checked$ | $X$ | $\checked$
(c) | $X$ | $\checked$ | $\checked$ | $X$
(d) | $\checked$ | $X$ | $X$ | $\checked$

In practice, the AAIIW estimator can be computed by first estimating the
parameters from all nuisance models (using a logistic regression, a
proportional rate model with e.g., coxph in R, and two linear models for the
outcome means conditional on $K_{i}(t)$ or $V_{i}(t)$). The estimates can be
plugged into the estimating equations of the AAIIW. A root solver such as
uniroot in R can be used then to estimate $\beta_{0}$ and that estimate be
further plugged into the second estimating equation, which is solved for
$\beta_{1}$.

###  2.5 Efficiency and Asymptotic Properties

The asymptotic variance of the AAIIW estimator can be derived using theory on
two-step estimators (Newey & McFadden, 1994) or deriving the variance of its
influence function (Tsiatis, 2006). We use the latter approach. Derivations
are shown in Supplementary Material B, in which we also show that the AAIIW
estimator is more efficient than the FIPTM estimator when all nuisance models
are correctly specified for both estimators. We find that under correct
nuisance models specification, the FIPTM asymptotic variance equals to

$\displaystyle\sigma^{2}_{FIPTM}=E\left[\frac{\\{Y^{1}_{i}(t)-\mu_{1}\\}^{2}}{\rho\\{V_{i}(t)\\}e_{1}\\{K_{i}(t)\\}}\right]+E\left[\frac{\\{Y^{0}_{i}(t)-\mu_{0}\\}^{2}}{\rho\\{V_{i}(t)\\}e_{0}\\{K_{i}(t)\\}}\right],$

where $\mu_{a}=E[Y^{a}_{i}(t)]$, $\rho\\{V_{i}(t)\\}=E[\rm{d}N_{i}(t)\mid
V_{i}(t);\gamma_{0}]$, and $e_{a}\\{K_{i}(t)\\}=pr\\{A_{i}(t)=a\mid
K_{i}(t);\psi_{0}\\}$. The augmented AAIIW estimator is more efficient, with
asymptotic variance

$\displaystyle\sigma^{2}_{AAIIW}=\sigma^{2}_{FIPTM}-\mu_{1}^{2}E\left[\frac{1+e_{1}\\{K_{i}(t)\\}}{e_{1}\\{K_{i}(t)\\}}\right]-\mu_{0}^{2}E\left[\frac{1+e_{0}\\{K_{i}(t)\\}}{e_{0}\\{K_{i}(t)\\}}\right].$

Furthermore, the theory of influence function can be used to show the
asymptotic normality of the estimator. The AAIIW estimator is consistent and
converges to the true causal effect of the binary treatment, so its limiting
distribution is normally distributed around $\beta_{1}$.

##  3 Simulation study

###  3.1 Comparators

We compared four different estimators in large simulation studies: an ordinary
least squares estimator that does not account for confounding or informative
observation (OLS); an inverse probability of treatment-weighted estimator that
accounts for the confounding process with a propensity score that is either
correctly specified (IPTc) or not correct (IPTnc); a doubly-weighted estimator
that corresponds to the FIPTM estimator from Coulombe et al. (2021b) with
either both the inverse intensity of visit and inverse probability of
treatment weights correctly specified (DWc), only the inverse probability of
treatment weights correctly specified (DWiptc), only the inverse intensity of
visit weights correctly specified (DWiivc), or both weights misspecified
(DWnc); and the novel AAIIW estimator with either all four nuisance models
correctly specified (AAIIWc), both weight models correctly specified and both
conditional outcome models misspecified (AAIIWs.a, with s.a referring to
scenario a in Table 2), both conditional outcome models correctly specified
and both weights misspecified (AAIIWs.b), the inverse intensity of visit
weights and the outcome model conditional on confounders being the only
correctly specified models (AAIIWs.c), and the inverse probability of
treatment weights and the outcome model conditional on observation predictors
being the only correctly specified models (AAIIWs.d).

The data-generating mechanism was strongly inspired by similar simulation
studies presented in Buzkova & Lumley (2009); Coulombe et al. (2021b, 2022)
and is described in much more detail in Supplementary Material C. The data-
generating mechanism included a set of confounders at baseline repeated
through follow-up, a time-varying binary treatment, a set of observation
predictors that varied in time, and irregular observation of the outcome. The
main results for 1000 simulations using a nonhomogeneous Poisson rate to
simulate the observation times of the outcome and a sample of size 1000 are
presented in the following section. In another set of simulations, we replaced
the nonhomogeneous Poisson rate with a nonhomogeneous Bernoulli probability
for the observation indicator and used a logistic regression instead of the
Andersen and Gill model to fit the probability of observation at each time
point. These results are presented in Supplementary Material D (Suppl. Fig. 2
and 3), along with the results under a sample of size 250 instead of 1000
(Suppl. Fig. 1) and all Monte Carlo biases and mean square errors (Suppl.
Table 1). In both settings using either the Poisson rate of the Bernoulli
probability, we tested four different sets of $\gamma$ parameters in the
observation model, including one set of zeros (which we call “set 1” in the
results) corresponding to uninformative observation.

###  3.2 Results

The distributions of 1000 estimates obtained with each estimator using a
sample of size 1000 patients are presented in Fig. 2.

The results are as expected. First, the ordinary least squares estimator is
strongly biased in all $\gamma$ parameter settings 1) to 4). In scenario 1) in
which we expected no bias due to the visit process, the inverse probability of
treatment-weighted estimator IPTc is empirically unbiased and the inverse
probability of treatment-weighted estimator using a wrong treatment model,
IPTnc, exhibits bias. In scenarios 2) to 4), both inverse probability of
treatment-weighted estimators are biased since they do not account properly
for the outcome observation process.

The doubly-weighted estimator DWc is consistently unbiased as it accounts
properly for both types of bias. When the visit process is uninformative, in
scenario 1) for the observation process, it is also unbiased even when the
inverse intensity of visit weights are not correctly specified, as long as the
inverse probability of treatment weights are correctly specified. In scenarios
2) to 4), the doubly-weighted estimator is only unbiased in settings in which
its two weight models are correctly specified.

The multiply robust AAIIW estimator is empirically unbiased in all scenarios
1) to 4) for the observation process, whenever using one of the four
combinations of correctly specified models shown in Table 2 or when all four
nuisance models are correctly specified. It exhibits particularly small
variance when the two conditional outcome mean models are correctly specified
(scenario b from Table 2) or, as expected when all four models are correctly
specified.

Results for the second set of simulations using the Bernoulli probability to
simulate the observations, and those for a sample of size 250 are presented in
Supplementary Material D. As expected, the estimators are more variable when
using a sample size of 250, although the same patterns in the comparison of
estimators are observed (Suppl. Fig. 1). Similar results are observed when
using the Bernoulli probability instead of the Poisson rate for the simulation
of observation indicators (Suppl. Fig. 2 and 3). The simulations using the
Bernoulli probability did not require the use of Breslow’s estimator for the
baseline rate, which would partly explain the smaller variances observed
overall (e.g., compare Fig. 1 and Suppl. Fig. 2).

Figure 2: Results of the simulation studies with a sample size of 1000 using a
nonhomogeneous Poisson rate to simulate the observation indicators and the
Andersen and Gill model with Breslow estimator to estimate the inverse
intensity of visit weights. Each boxplot represents the distribution of 1000
estimates for the corresponding estimator. The dashed line represents the gold
standard, i.e., the true value for the marginal effect of exposure that equals
to 1. Different strengths of the visit process on covariates are represented
with scenarios 1) $\gamma=(0,0,0,0,0,-5)$ (i.e., no bias due to the visit
process expected); 2) $\gamma=(0.5,0.3,-0.5,-2,0,-3)$; 3)
$\gamma=(0.5,-0.5,-0.2,-1,1,-3)$; and 4) $\gamma=(-1,-0.8,0.1,0.3,-1,-3)$.

##  4 Motivating Example

We applied the proposed AAIIW estimator and different more naive comparators
to data from the Add Health study in the United State (Harris, 2009, 2013;
Harris & Udry, 2022). It consists of a longitudinal study with multiple waves.
The study started in 1994 when a pool of adolescents representative from the
United States was selected. These adolescents grew up to become adults during
the study. At each wave, they were asked to fill out questionnaires.

We have access to data from the first four waves of the Add Health study,
corresponding to the years 1994-1995, 1996, 2001-2002 and 2008-2009
respectively. Although these data do not consist of electronic health records,
and that adolescents in the Add Health study were observed at pre-specified
observation times, they tended not to be observed at each wave for all
variables measured in the study, and one could consider that the four waves
represent four consecutive time points when there could be or no a visit or an
observation time, just like in medical records. Various types of information
such as demographics, health status, nutrition, family dynamics, sexual
activity and substance use were collected for the study. Some questions varied
across the four waves and we focused on a causal research question for which
we have data at all four waves. Our goal was to estimate the marginal causal
effect of counselling (psychotherapy) on alcohol consumption. We think that
the effect of counselling on alcohol consumption is mediated by the depressive
mood of adolescents and that their mood can be affected by counselling and may
in its turn affect alcohol consumption (see the assumed data generating
mechanism in Fig. 1(b)). Two important challenges we wished to consider in the
analysis are the irregular observation of the outcome and, because the study
is observational, the potential confounding of the psychotherapy-alcohol
consumption relationship.

We selected potential confounders for that relationship, which included the
teen’s age, sex, socioeconomic status, weight in pounds, and whether they
smoked at least once in the previous month. The socioeconomic status was
computed by summing two variables that we transformed, namely the parents’
total income in 1994 before taxes and one of the parents’ education, usually
that of the resident mother (Harris, 2009). The parent’s total income was
transformed into quintiles (1 to 5 with 5 being the highest). One of the
parents’ education was categorized in 5 levels corresponding to 1- 8th grade
or less or never went to school, 2-more than 8th grade but did not graduate
high school, 3-went to a business, trade or vocational school instead of high
school, high school graduate, or completed a general educational development
program, 4–went to a business, trade or vocational school after high school,
went to college but did not graduate, or 5-graduated from college or
university, or professional training beyond a 4-year college or university
training. Socioeconomic status was defined as the sum of the two transformed
5-category variables.

The analysis dataset contained several missing values. Unless we had enough
information in the dataset to replace missing values in the variables age and
sex (e.g., if age was measured at a previous wave and it could be used for
extrapolation), we used multiple imputations by chained equations (Rubin,
1988; Azur et al., 2011) to impute missing values in these variables as well
as in variables socioeconomic status, smoking status, weight, depressive mood,
and the exposure to counselling. We used these variables as predictors each
time to impute each variable using fully conditional specification. The
alcohol outcome was kept as missing when it was not measured.

The outcome was defined using the question Think of all the times you had a
drink during the past 12 months. How many drinks did you usually have each
time?. It consisted of a self-assessed number of drinks the adolescent would
consume, on average, each time they consumed alcohol. It ranged from 0 to 90.
In this application, the outcome tended to be assessed at each of the four
waves (i.e., not irregularly and with most data not being missing in the
outcome). To assess the advantage of our approach, we simulated missingness in
the outcome and assessed the different estimators in that setting, knowing the
true underlying missingness mechanism.

Assuming that all potential confounders as well as the mediator (depressive
mood) and the exposure (counselling) affect the chance of observing the
alcohol consumption outcome, the outcome observation (i.e., the opposite of
missingness) was simulated using the following model across the four waves:

$\displaystyle E[\rm{d}N_{i}(t)\mid\text{age, sex, counselling, depressive
mood, socioeconomic status}]=\text{expit}\\{18-$ $\displaystyle 0.3\hskip
2.84544pt\text{age}(t)+0.8\hskip 2.84544pt\text{sex}+1.8\hskip
2.84544pt\text{counselling}(t)-3\hskip 2.84544pt\text{depressive
mood}(t)-1.3\hskip 2.84544pt\text{socioeconomic status}\\}$

for $t\in 1,2,3,4$, where
expit$(\cdot)=\exp(\cdot)/\left\\{1+\exp(\cdot)\right\\}$.

We conducted the analysis and fit a propensity score model as a function of
age, sex, weight, socioeconomic status and smoking. We fit two different
proportional rate models for the observation of the outcome, one correctly
specified (as a function of age, sex, counselling, depressive mood, and
socioeconomic status) and one that was not correctly specified (as a function
of the sinus of age and the depressive mood, therefore not including the right
format for the age variable and missing some important variables in the
model). The estimators compared in the application are a standard ordinary
least squares estimator (not adjusted for confounding nor for the observation
process), an inverse intensity of visit-weighted estimator that does not
account for confounding but does account for the observation process (we
tested the two sets of the inverse intensity of visit weights), a doubly-
weighted estimator corresponding to the FIPTM estimator (incorporating the
inverse probability of treatment weights based on our assumptions on the
potential confounders, and inverse intensity of visit weights - we tested the
two sets of the inverse intensity of visit weights here again), and the AAIIW
estimator in which we incorporated the inverse probability of treatment
weights and the two different sets of the inverse intensity of visit weights,
one at a time.

We found an important difference in the proportion of males and females across
both exposure groups, with females reporting higher rates of counselling,
$14\%$ more smokers in the counselling group than the other, and a greater
depressive mood in those receiving counselling (Supplementary Material E,
Suppl. Table 2). These variables were, therefore, considered potential
confounders in our analysis, except for the mediator depressive mood that
should not be conditioned upon in the confounding set. After inverse
probability of treatment-weighting, the two exposure groups are similar with
respect to all potential confounders (Suppl. Table 2).

In the outcome observation model, we found modest differences in the
depressive mood and counselling rates between those for whom the alcohol
consumption was observed and the others (Supplementary Material E, Suppl.
Table 3). There was also a slight difference in age means and female sex
proportions, as expected (since we simulated the missingness mechanism
ourselves). After inverse intensity of visit weighting, most differences
vanished, with age means, depressive mood means, and female proportions that
are much closer (Suppl. Table 3).

Table 3:  Estimates (95% bootstrap percentiles confidence intervals) of the
marginal effect of counselling on the average number of alcoholic beverages
consumed, Add Health study, United States, 1996-2008

OLS | IPTϕ | IIV† | IIV‡
---|---|---|---
0.62 (0.39, 0.75) | 0.34 (0.15, 0.48) | 0.64 (0.40, 0.77) | 0.72 (0.49, 0.92)
FIPTMϕ,† | FIPTMϕ,‡ | AAIIWϕ,† | AAIIWϕ,‡
0.35 (0.12, 0.50) | 0.46 (0.24, 0.67) | 0.35 (0.12, 0.51) | 0.28 (0.04, 0.54)

Acronyms: IPT, inverse probability of treatment; IIV: inverse intensity of
visit; FIPTM: the flexible inverse probability of treatment and monitoring;
AAIIW: the doubly augmented, doubly inverse weighted.
ϕ. Note we do not know the true data-generating mechanism for the treatment
mechanism in the application.
†. This estimator uses a correctly specified generating mechanism for outcome
missingness. ‡. This estimator uses a wrongly specified generating mechanism
for outcome missingness.

In this application using data from the Add Health study, both the adjustment
for confounding and the one for outcome missingness bring the estimates for
the marginal effect of exposure to counselling towards the null. For instance,
the inverse probability of treatment-weighted estimator is closer to $0$ than
the standard ordinary least squares estimator that does not adjust for
confounding or the observation process. The estimator using the correctly
adjusted inverse intensity of visit weights (IIV†) slightly brings the
estimate towards the null when compared with using the wrong inverse intensity
of visit weights (IIV‡). Combining both adjustments, the estimate for the
causal marginal effect of counselling on alcohol consumption goes from $0.62$
($95\%$ confidence interval $0.39$, $0.75$) with no adjustment at all, to
$0.35$ ($0.12$, $0.51$) when using the correct inverse intensity of visit
weights (and a propensity score based on our assumptions) in the multiply
robust AAIIW estimator.

Of most interest is the comparison of the FIPTM and the AAIIW estimators in
this application. Using a wrong model for the calculation of inverse intensity
of visit weights leads to an estimate of $0.46$ ($95\%$ confidence interval
$0.24$, $0.67$) with the FIPTM and $0.28$ ($95\%$ confidence interval $0.04$,
$0.54$) with the AAIIW estimator. Both estimators lead to similar estimates
when using the correctly specified inverse intensity of visit weights. Those
results indicate that in a setting in which we would not know the true
observation mechanism, the AAIIW estimator might still lead to an estimate of
the causal effect closer to the null (around $0.2$ or $0.3$ here) while we
know that the FIPTM risks being biased when its weights are not correctly
modelled. This indicates that the true effect is probably closer to $0.2$ or
$0.3$ in this application (and indeed, the $95\%$ confidence interval for the
AAIIW estimator using the wrong inverse intensity of visit weights, which may
be protected against bias since its outcome mean model conditional on the
observation predictors has a chance at being correctly specified, is realistic
with a lower bound closer to $0$).

##  5 Discussion

This work proposed the first quadruply robust estimator that is consistent
when only two out of four nuisance models (the treatment, the visit, and two
conditional outcome mean models) are correctly specified, as long as the
combination of correctly specified models is one shown in Table 2. The
proposed estimator is particularly useful for observational settings subject
to confounding and irregular observation times when researchers suspect they
can specify correctly at least one of two models related to confounders
(either the prescription mechanism, or the outcome biological mechanism), and
one of two models related to the outcome observation (either the observation
of the outcome or the outcome biological mechanism). The estimator also allows
mediators and other variables that could be on the causal path from the
treatment to the outcome to affect the observation model and yet be adjusted
for.

In addition to being more robust than its only alternative, the FIPTM which is
a doubly-weighted estimator, the AAIIW estimator is also the most efficient
estimator in its class of semiparametric estimators. In simulation studies,
the AAIIW was demonstrated to be robust and empirically as efficient as the
FIPTM when the two weight models are correctly specified but much more
efficient in other scenarios (such as when all four models used in its
estimating equations are correctly specified or when the two outcome mean
models conditional on confounders or observation predictors are correctly
specified).

In an application to the Add Health study in the United States, we have found
a difference between more naive estimators and the multiply robust AAIIW
estimator in the estimation of the causal marginal effect of therapy
counselling on alcohol consumption (e.g., a causal effect of 0.35 more drinks
with counselling therapy (95% confidence interval 0.12, 0.51) versus 0.62 more
drinks (95% confidence interval 0.39, 0.75) with the most naive estimator). It
is possible that unmeasured confounding remains and that a better adjustment
would have brought these estimates even closer to the null. Sensitivity
analyses could be used to assess the effect of unmeasured confounding or visit
predictors that were not accounted properly in the estimator (see Smith et al.
(2022) for informative observation or see e.g., Schneeweiss (2006);
VanderWeele & Arah (2011) for unmeasured confounding).

The consistency of our proposed estimator relies on specific combinations of
correctly specified nuisance models listed in Table 2 and some classical
causal assumptions mentioned in $\mathsection$2. An analyst using the novel
approach, in collaboration with an expert from the substantive research field,
should identify the confounders of the relationship between the exposure and
outcome and the outcome observation predictors at risk of creating spurious
associations between the exposure and the outcome (conditional exchangeability
assumption). The use of a causal diagram can help in depicting the open paths
by which dependencies that are not due to causal effects arise. It is not
enough to include these variables in the weight models (inverse probability of
treatment or inverse intensity of visit weights) or in the conditional outcome
mean models discussed in this manuscript. A model that is correctly specified
implies that the functional format of the predictors in the model is correctly
specified. Second, the analyst should ensure that the observations included in
the analysis meet the two positivity assumptions for the treatment and the
observation models. They might decide to remove from the analysis a patient
with a set of characteristics that are only represented in one of the two
exposure groups (otherwise, the positivity of treatment assumption may be
violated) or with a set of characteristics that are only represented when the
outcome is not observed (otherwise violating the positivity of observation
assumption). Third, the trivial causal assumption of consistency of the
outcome must be met. In our application to the Add Health study, for instance,
we assume that the outcome observed in those who claimed to have received
counselling is truly equal to their potential outcome under counselling, and
vice-versa. Measurement errors, or patients not having filled in the
information properly or truthfully could alter the estimates.

Interesting future avenues of work include the use of more flexible methods,
such as machine learning methods, to correctly model the outcome mean as a
function of confounders or observation predictors. The framework of targeted
maximum likelihood estimation (Van Der Laan & Rubin, 2006; Schuler & Rose,
2017) could be used for that purpose. Another advantage of the proposed
estimator is that it is based on the general framework of generalized
estimating equations (Zeger & Liang, 1986) and could therefore be extended
straightforwardly to account for other types of outcomes (e.g., binary
outcomes). The extension to categorical or continuous exposure is also
possible.

## Acknowledgement

We thank Professor Marie Davidian at North Carolina State University for the
enriching discussions we had as part of a virtual research internship of
author JC in the Department of Statistics at North Carolina State University
in 2021-2022.

This research uses data from Add Health, funded by grant P01 HD31921 (Harris)
from the Eunice Kennedy Shriver National Institute of Child Health and Human
Development (NICHD), with cooperative funding from 23 other federal agencies
and foundations. Add Health is currently directed by Robert A. Hummer and
funded by the National Institute on Aging cooperative agreements U01 AG071448
(Hummer) and U01AG071450 (Aiello and Hummer) at the University of North
Carolina at Chapel Hill. Add Health was designed by J. Richard Udry, Peter S.
Bearman, and Kathleen Mullan Harris at the University of North Carolina at
Chapel Hill.

## Supplementary material

Supplementary material available at the end of this document includes:

  * –

Supplementary material A: Proof of the multiple robustness of the novel
estimator under the different scenarios presented in Table 2.

  * –

Supplementary material B: Sketch-proof for the efficiency of the AAIIW
estimator.

  * –

Supplementary material C: Details of the simulation studies.

  * –

Supplementary material D: Results of the simulation studies for a sample size
of 250 with the nonhomogeneous Poisson rate for the observation indicators,
for sample sizes 250 or 1000 using the Bernoulli probability for the
observation of the outcome, and Monte Carlo bias and mean square error of each
estimator, in all scenarios tested.

  * –

Supplementary material E: Tables of characteristics in the Add Health study
stratified by two weighting strategies.

## References

  * Andersen & Gill (1982) Andersen, P. K. & Gill, R. D. (1982).  Cox’s regression model for counting processes: A large sample study.  Ann. Stat. 10, 1100–1120.
  * Azur et al. (2011) Azur, M. J., Stuart, E. A., Frangakis, C. & Leaf, P. J. (2011).  Multiple imputation by chained equations: What is it and how does it work?  Int. J. Methods Psychiatr. Res. 20, 40–49.
  * Black (1996) Black, N. (1996).  Why we need observational studies to evaluate the effectiveness of health care.  BMJ 312, 1215–1218.
  * Buzkova & Lumley (2009) Buzkova, P. & Lumley, T. (2009).  Semiparametric modeling of repeated measurements under outcome-dependent follow-up.  Stat. Med. 28, 987–1003.
  * Cao et al. (2009) Cao, W., Tsiatis, A. A. & Davidian, M. (2009).  Improving efficiency and robustness of the doubly robust estimator for a population mean with incomplete data.  Biometrika 96, 723–734.
  * Coulombe et al. (2021a) Coulombe, J., Moodie, E. E. & Platt, R. W. (2021a).  Estimating the marginal effect of a continuous exposure on an ordinal outcome using data subject to covariate-driven treatment and visit processes.  Stat. Med. 40, 5746–5764.
  * Coulombe et al. (2021b) Coulombe, J., Moodie, E. E. & Platt, R. W. (2021b).  Weighted regression analysis to correct for informative monitoring times and confounders in longitudinal studies.  Biometrics 77, 162–174.
  * Coulombe et al. (2022) Coulombe, J., Moodie, E. E., Platt, R. W. & Renoux, C. (2022).  Estimation of the marginal effect of antidepressants on body mass index under confounding and endogenous covariate-driven monitoring times.  Ann. Appl. Stat. 16, 1868–1890.
  * Cox (1972) Cox, D. R. (1972).  Regression models and life-tables.  J. R. Statist. Soc. B 34, 187–202.
  * Dahabreh & Kent (2014) Dahabreh, I. J. & Kent, D. M. (2014).  Can the learning health care system be educated with observational data?  JAMA 312, 129–130.
  * Farewell (2010) Farewell, D. (2010).  Marginal analyses of longitudinal data with an informative pattern of observations.  Biometrika 97, 65–78.
  * Funk et al. (2011) Funk, M. J., Westreich, D., Wiesen, C., Stürmer, T., Brookhart, M. A. & Davidian, M. (2011).  Doubly robust estimation of causal effects.  Am. J. Epidemiol. 173, 761–767.
  * Greenland (2003) Greenland, S. (2003).  Quantifying biases in causal models: Classical confounding vs collider-stratification bias.  Epidemiology 14, 300–306.
  * Greenland & Morgenstern (2001) Greenland, S. & Morgenstern, H. (2001).  Confounding in health research.  Annu. Rev. Public Health 22, 189–212.
  * Harris (2009) Harris, K. M. (2009).  The National Longitudinal Study of Adolescent to Adult Health (Add Health), Waves I & II, 1994–1996; Wave III, 2001–2002; Wave IV, 2007-2009 [machine-readable data file and documentation].  Chapel Hill: Carolina Population Center, University of North Carolina at Chapel Hill .
  * Harris (2013) Harris, K. M. (2013).  The Add Health study: Design and accomplishments.  Chapel Hill: Carolina Population Center, University of North Carolina at Chapel Hill , 1–22.
  * Harris & Udry (2022) Harris, K. M. & Udry, J. R. (2022).  The National Longitudinal Study of Adolescent to Adult Health (Add Health), Waves I & II, 1994–1996; Wave III, 2001–2002; Wave IV, 2007-2009 [machine-readable data file and documentation].  Inter-university Consortium for Political and Social Research [distributor]. https://www.icpsr.umich.edu/web/DSDR/studies/21600/versions/V25 (accessed March 10, 2023) .
  * Holland (1986) Holland, P. W. (1986).  Statistics and causal inference.  J. Am. Statist. Assoc. 81, 945–960.
  * Horvitz & Thompson (1952) Horvitz, D. G. & Thompson, D. J. (1952).  A generalization of sampling without replacement from a finite universe.  J. Am. Statist. Assoc. 47, 663–685.
  * Jiang et al. (2022) Jiang, Z., Yang, S. & Ding, P. (2022).  Multiply robust estimation of causal effects under principal ignorability.  J. R. Statist. Soc. B 84, 1423–1445.
  * Krumholz (2014) Krumholz, H. M. (2014).  Big data and new knowledge in medicine: The thinking, training, and tools needed for a learning health system.  Health Aff. 33, 1163–1170.
  * Lin & Ying (2001) Lin, D. Y. & Ying, Z. (2001).  Semiparametric and nonparametric regression analysis of longitudinal data.  J. Am. Statist. Assoc. 96, 103–126.
  * Lin et al. (2004) Lin, H., Scharfstein, D. O. & Rosenheck, R. A. (2004).  Analysis of longitudinal data with irregular, outcome-dependent follow-up.  J. R. Statist. Soc. B 66, 791–813.
  * Lipsitz et al. (2002) Lipsitz, S. R., Fitzmaurice, G. M., Ibrahim, J. G., Gelber, R. & Lipshultz, S. (2002).  Parameter estimation in longitudinal studies with outcome-dependent follow-up.  Biometrics 58, 621–630.
  * Lokku et al. (2020) Lokku, A., Lim, L. S., Birken, C. S. & Pullenayegum, E. M. (2020).  Summarizing the extent of visit irregularity in longitudinal data.  BMC Medical Res. Methodol. 20, 1–9.
  * Naimi et al. (2017) Naimi, A. I., Cole, S. R. & Kennedy, E. H. (2017).  An introduction to G methods.  Int. J. Epidemiol. 46, 756–762.
  * Newey & McFadden (1994) Newey, W. K. & McFadden, D. (1994).  Large sample estimation and hypothesis testing.  Handbook of econometrics 4, 2111–2245.
  * Neyman (1923) Neyman, J. S. (1923).  On the application of probability theory to agricultural experiments. Essay on principles. Section 9 (translation published in 1990).  Stat. Sci. 5, 472–480.
  * Pearl (2009) Pearl, J. (2009).  Causality.  Cambridge: Cambridge university press.
  * Pullenayegum & Lim (2016) Pullenayegum, E. M. & Lim, L. S. (2016).  Longitudinal data subject to irregular observation: A review of methods with a focus on visit processes, assumptions, and study design.  Stat. Methods Med. Res. 25, 2992–3014.
  * 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.  Math. Model. 7, 1393–1512.
  * Robins & Rotnitzky (1998) Robins, J. & Rotnitzky, A. (1998).  Discussion of “Robust models in probability sampling” by Firth and Bennett.  J. R. Statist. Soc. B , 51–52.
  * Robins & Hernán (2008) Robins, J. M. & Hernán, M. A. (2008).  Estimation of the causal effects of time-varying exposures.  In Longitudinal data analysis, chap. 3. Boca Raton: CRC press, p. 579.
  * Robins et al. (2000) Robins, J. M., Hernán, M. A. & Brumback, B. (2000).  Marginal structural models and causal inference in epidemiology.  Epidemiology 11, 550–560.
  * Robins et al. (1994) Robins, J. M., Rotnitzky, A. & Zhao, L. P. (1994).  Estimation of regression coefficients when some regressors are not always observed.  J. Am. Statist. Assoc. 89, 846–866.
  * Rosenbaum (1987) Rosenbaum, P. R. (1987).  Model-based direct adjustment.  J. Am. Statist. Assoc. 82, 387–394.
  * Rosenbaum & Rubin (1983) Rosenbaum, P. R. & Rubin, D. B. (1983).  The central role of the propensity score in observational studies for causal effects.  Biometrika 70, 41–55.
  * Rubin (1976) Rubin, D. B. (1976).  Inference and missing data.  Biometrika 63, 581–592.
  * Rubin (1988) Rubin, D. B. (1988).  An overview of multiple imputation.  In Proceedings of the survey research methods section of the American Statistical Association, vol. 79. Princeton: Citeseer.
  * Rubin (2008) Rubin, D. B. (2008).  For objective causal inference, design trumps analysis.  Ann. Appl. Stat. 2, 808–840.
  * Schneeweiss (2006) Schneeweiss, S. (2006).  Sensitivity analysis and external adjustment for unmeasured confounders in epidemiologic database studies of therapeutics.  Pharmacoepidemiol. Drug. Saf. 15, 291–303.
  * Schuler & Rose (2017) Schuler, M. S. & Rose, S. (2017).  Targeted maximum likelihood estimation for causal inference in observational studies.  Am. J. Epidemiol. 185, 65–73.
  * Smith et al. (2022) Smith, B., Yang, S., Apter, A. J. & Scharfstein, D. O. (2022).  Trials with irregular and informative assessment times: A sensitivity analysis approach.  arXiv:2204.11979 .
  * Tsiatis (2006) Tsiatis, A. A. (2006).  Semiparametric theory and missing data.  New York: Springer.
  * Tsiatis & Davidian (2007) Tsiatis, A. A. & Davidian, M. (2007).  Comment: Demystifying double robustness: A comparison of alternative strategies for estimating a population mean from incomplete data.  Stat. Sci. 22, 569–573.
  * Van Der Laan & Rubin (2006) Van Der Laan, M. J. & Rubin, D. (2006).  Targeted maximum likelihood learning.  Int. J. Biostat. 2, 1–38.
  * VanderWeele & Arah (2011) VanderWeele, T. J. & Arah, O. A. (2011).  Bias formulas for sensitivity analysis of unmeasured confounding for general outcomes, treatments, and confounders.  Epidemiology 22, 42–52.
  * Yang (2022) Yang, S. (2022).  Semiparametric estimation of structural nested mean models with irregularly spaced longitudinal observations.  Biometrics 78, 937–949.
  * Zeger & Liang (1986) Zeger, S. L. & Liang, K.-Y. (1986).  Longitudinal data analysis for discrete and continuous outcomes.  Biometrics 42, 121–130.

Supplementary Material for “Quadruply robust estimation of marginal structural
models in observational studies subject to covariate-driven observations”

Supplementary Material A

Proof of the multiple robustness of the novel estimator under the different
scenarios presented in Table 2 in the main manuscript

The probability limiting estimating equations of the AAIIW estimator are

$\displaystyle
pr\left[\int_{0}^{\tau}\frac{\frac{\mathbf{1}\\{A_{i}(t)=a\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}Y_{i}(t){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-\frac{\mathbf{1}\\{A_{i}(t)=a\\}-pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}\mu_{a}\\{K_{i}(t);\hat{\alpha}_{K}\\}}-\zeta_{i}(t;\beta_{a})}{E\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\hat{\gamma}\right\\}}\rm{d}N_{i}(t)\right]$ $\displaystyle-
pr\left[\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\int_{0}^{\tau}\frac{\rm{d}M_{i}(t)\left(\frac{\mu_{a}\\{V_{i}(t);\hat{\alpha}_{V}\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}-\frac{\mathbf{1}\\{A_{i}(t)=a\\}-pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}\mu_{a}\\{K_{i}(t);\hat{\alpha}_{K}\\}-\zeta_{i}(t;\beta_{a})\right)}{E\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\hat{\gamma}\right\\}}\right]=0.$

In this proof, denote the correct models using an asterisk onto the
parameters, i.e., the correct model for the treatment is denoted by
$pr\left\\{A_{i}(t)=a\mid K_{i}(t);\psi^{*}\right\\}$.

Under scenario a) (Table 2) we have

$\displaystyle
pr\left[\int_{0}^{\tau}\frac{\frac{\mathbf{1}\\{A_{i}(t)=a\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\psi^{*}\\}}Y_{i}(t){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-\frac{\mathbf{1}\\{A_{i}(t)=a\\}-pr\\{A_{i}(t)=a\mid
K_{i}(t);\psi^{*}\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\psi^{*}\\}}\mu_{a}\\{K_{i}(t);\hat{\alpha}_{K}\\}}-\zeta_{i}(t;\beta_{a})}{pr\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\gamma^{*}\right\\}}\rm{d}N_{i}(t)\right]$ $\displaystyle-
pr\left[\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\int_{0}^{\tau}\frac{(\rm{d}N_{i}(t)-pr\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\gamma^{*}\right\\})\left(\frac{\mu_{a}\\{V_{i}(t);\hat{\alpha}_{V}\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\psi^{*}\\}}-\frac{\mathbf{1}\\{A_{i}(t)=a\\}-pr\\{A_{i}(t)=a\mid
K_{i}(t);\psi^{*}\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\psi^{*}\\}}\mu_{a}\\{K_{i}(t);\hat{\alpha}_{K}\\}\right)}{pr\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\gamma^{*}\right\\}}\right]$
$\displaystyle+pr\left[\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\int_{0}^{\tau}\frac{(\rm{d}N_{i}(t)-pr\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\gamma^{*}\right\\})\zeta_{i}(t;\beta_{a})}{pr\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\gamma^{*}\right\\}}\right]$
$\displaystyle=pr\left[\int_{0}^{\tau}\frac{\frac{\mathbf{1}\\{A_{i}(t)=a\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\psi^{*}\\}}Y_{i}(t)-\zeta_{i}(t;\beta_{a})}{pr\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\gamma^{*}\right\\}}\rm{d}N_{i}(t)\right]$ $\displaystyle=0$

which is clearly unbiased for the parameter of interest since both weights are
correctly specified in scenario a).
In scenario b) (Table 2), we rewrite the original estimating equations as
follows:

$\displaystyle
pr\left[\int_{0}^{\tau}\frac{\frac{\mathbf{1}\\{A_{i}(t)=a\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}(Y_{i}(t){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-\mu_{a}\\{K_{i}(t);\hat{\alpha}_{K}\\}})+\mu_{a}\\{K_{i}(t);\hat{\alpha}_{K}\\}-\zeta_{i}(t;\beta_{a})}{E\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\hat{\gamma}\right\\}}\rm{d}N_{i}(t)\right]$ $\displaystyle-
pr\left[\int_{0}^{\tau}\frac{\frac{\mathbf{1}\\{A_{i}(t)=a\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}(\mu_{a}\\{V_{i}(t);\hat{\alpha}_{V}\\}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-\mu_{a}\\{K_{i}(t);\hat{\alpha}_{K}\\}})+\mu_{a}\\{K_{i}(t);\hat{\alpha}_{K}\\}-\zeta_{i}(t;\beta_{a})}{E\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\hat{\gamma}\right\\}}\rm{d}N_{i}(t)\right]$
$\displaystyle+pr\left[\int_{0}^{\tau}\frac{\frac{\mathbf{1}\\{A_{i}(t)=a\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}(\mu_{a}\\{V_{i}(t);\hat{\alpha}_{V}\\}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-\mu_{a}\\{K_{i}(t);\hat{\alpha}_{K}\\}})+\mu_{a}\\{K_{i}(t);\hat{\alpha}_{K}\\}-\zeta_{i}(t;\beta_{a})}{E\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\hat{\gamma}\right\\}}E\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\hat{\gamma}\right\\}\right].$

Then, using the asterisk notation for the correct models once again, we obtain

$\displaystyle
pr\left[\int_{0}^{\tau}\frac{\frac{\mathbf{1}\\{A_{i}(t)=a\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}(Y_{i}(t){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-\mu_{a}\\{K_{i}(t);\alpha^{*}_{K}\\}})+\mu_{a}\\{K_{i}(t);\alpha^{*}_{K}\\}-\zeta_{i}(t;\beta_{a})}{E\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\hat{\gamma}\right\\}}\rm{d}N_{i}(t)\right]$ $\displaystyle-
pr\left[\int_{0}^{\tau}\frac{\frac{\mathbf{1}\\{A_{i}(t)=a\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}(\mu_{a}\\{V_{i}(t);\alpha^{*}_{V}\\}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-\mu_{a}\\{K_{i}(t);\alpha^{*}_{K}\\}})+\mu_{a}\\{K_{i}(t);\alpha^{*}_{K}\\}-\zeta_{i}(t;\beta_{a})}{E\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\hat{\gamma}\right\\}}\rm{d}N_{i}(t)\right]$
$\displaystyle+\small{pr\left[\int_{0}^{\tau}\frac{\frac{\mathbf{1}\\{A_{i}(t)=a\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}(\mu_{a}\\{V_{i}(t);\alpha^{*}_{V}\\}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-\mu_{a}\\{K_{i}(t);\alpha^{*}_{K}\\}})+\mu_{a}\\{K_{i}(t);\alpha^{*}_{K}\\}-\zeta_{i}(t;\beta_{a})}{E\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\hat{\gamma}\right\\}}E\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\hat{\gamma}\right\\}\right]}$
$\displaystyle=pr\left[\int_{0}^{\tau}\frac{\frac{\mathbf{1}\\{A_{i}(t)=a\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}(Y_{i}(t){-\mu_{a}\\{V_{i}(t);\alpha^{*}_{V}\\}})}{E\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\hat{\gamma}\right\\}}\rm{d}N_{i}(t)\right]\text{(by combining the
first two lines).}$

The last equation is unbiased since the model
$\mu_{a}\\{V_{i}(t);\alpha^{*}_{V}\\}$ is correctly specified .

In scenario c) (Table 2), we have the following equations once plugging-in the
correct specified models:

$\displaystyle
pr\left[\int_{0}^{\tau}\frac{\frac{\mathbf{1}\\{A_{i}(t)=a\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}\left[Y_{i}(t)-\mu_{a}\\{K_{i}(t);\alpha^{*}_{K}\\}\right]+\mu_{a}\\{K_{i}(t);\alpha^{*}_{K}\\}-\zeta_{i}(t;\beta_{a})}{pr\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\gamma^{*}\right\\}}\rm{d}N_{i}(t)\right]$ $\displaystyle-
pr\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\int_{0}^{\tau}\frac{\rm{d}N_{i}(t)-pr\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\gamma^{*}\right\\}}{pr\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\gamma^{*}\right\\}}\left(\frac{\mu_{a}\\{V_{i}(t);\hat{\alpha}_{V}\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}\right)$
$\displaystyle\small{+pr\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\int_{0}^{\tau}\frac{\rm{d}N_{i}(t)-pr\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\gamma^{*}\right\\}}{pr\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\gamma^{*}\right\\}}\left(\frac{\mathbf{1}\\{A_{i}(t)=a\\}-pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}\mu_{a}\\{K_{i}(t);\alpha^{*}_{K}\\}-\zeta_{i}(t;\beta_{a})\right)}$
$\displaystyle=pr\left[\int_{0}^{\tau}\frac{\mu_{a}\\{K_{i}(t);\alpha^{*}_{K}\\}-\zeta_{i}(t;\beta_{a})}{pr\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\gamma^{*}\right\\}}\rm{d}N_{i}(t)\right]=0.$

The two last rows in the first development above cancel out because the
martingale residuals are 0-mean under correctly specified IIV weights. The
last equation above is unbiased for the causal effect since the outcome model
conditional on the confounders and the IIV weights are correctly specified.

Finally, in scenario d), we transform the original equations into:

$\displaystyle
pr\left[\int_{0}^{\tau}\frac{\frac{\mathbf{1}\\{A_{i}(t)=a\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}\left[Y_{i}(t)-\mu_{a}\\{V_{i}(t);\hat{\alpha}_{V}\\}\right]}{E\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\hat{\gamma}\right\\}}\rm{d}N_{i}(t)\right]$
$\displaystyle+pr\left[\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\int_{0}^{\tau}\left(\frac{\mathbf{1}\\{A_{i}(t)=a\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\hat{\psi}\\}}\left[\mu_{a}\\{V_{i}(t);\hat{\alpha}_{V}\\}-\mu_{a}\\{K_{i}(t);\hat{\alpha}_{K}\\}\right]+\mu_{a}\\{K_{i}(t);\hat{\alpha}_{K}\\}-\zeta_{i}(t;\beta_{a})\right)\rm{d}t\right]=0.$

Once replacing with the correctly specified models in scenario d), we obtain

$\displaystyle
pr\left[\int_{0}^{\tau}\frac{\frac{\mathbf{1}\\{A_{i}(t)=a\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\psi^{*}\\}}\left[Y_{i}(t)-\mu_{a}\\{V_{i}(t);\alpha^{*}_{V}\\}\right]}{E\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\hat{\gamma}\right\\}}\rm{d}N_{i}(t)\right]$
$\displaystyle+pr\left[\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\int_{0}^{\tau}\left(\frac{\mathbf{1}\\{A_{i}(t)=a\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\psi^{*}\\}}\left[\mu_{a}\\{V_{i}(t);\alpha^{*}_{V}\\}-\mu_{a}\\{K_{i}(t);\hat{\alpha}_{K}\\}\right]+\mu_{a}\\{K_{i}(t);\hat{\alpha}_{K}\\}-\zeta_{i}(t;\beta_{a})\right)\rm{d}t\right]$
$\displaystyle=pr\left[\int_{0}^{\tau}\frac{\frac{\mathbf{1}\\{A_{i}(t)=a\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\psi^{*}\\}}\left[Y_{i}(t)-\mu_{a}\\{V_{i}(t);\alpha^{*}_{V}\\}\right]}{E\left\\{\rm{d}N_{i}(t)=1\mid
V_{i}(t);\hat{\gamma}\right\\}}\rm{d}N_{i}(t)\right]$
$\displaystyle+pr\left[\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\int_{0}^{\tau}\left(\frac{\mathbf{1}\\{A_{i}(t)=a\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\psi^{*}\\}}\mu_{a}\\{V_{i}(t);\alpha^{*}_{V}\\}\right)dt\right]$
$\displaystyle-
pr\left[\left(\frac{\mathbf{1}\\{A_{i}(t)=a\\}-pr\\{A_{i}(t)=a\mid
K_{i}(t);\psi^{*}\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\psi^{*}\\}}\mu_{a}\\{K_{i}(t);\hat{\alpha}_{K}\\}+\zeta_{i}(t;\beta_{a})\right)\rm{d}t\right]$
$\displaystyle=pr\left[\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\int_{0}^{\tau}\left(\frac{\mathbf{1}\\{A_{i}(t)=a\\}}{pr\\{A_{i}(t)=a\mid
K_{i}(t);\psi^{*}\\}}\mu_{a}\\{V_{i}(t);\alpha^{*}_{V}\\}-\zeta_{i}(t;\beta_{a})\right)dt\right]=0.$

The last equation above is unbiased for the causal effect.

Supplementary Material B
Variance estimation and proof of efficiency of the AAIIW estimator compared to
the FIPTM.

To make the demonstrations lighter, denote

$\displaystyle I_{a}=$ $\displaystyle\mathbf{1}\\{A_{i}(t)=a\\}$
$\displaystyle e_{a}=$ $\displaystyle pr\\{A_{i}(t)=a\mid
K_{i}(t);\psi_{0}\\}$ $\displaystyle\rm{d}N=$
$\displaystyle\mathbf{1}\\{\rm{d}N_{i}(t)=1\\}$ $\displaystyle\rho=$
$\displaystyle E\\{\rm{d}N_{i}(t)=1\mid V_{i}(t);\gamma_{0}\\}$
$\displaystyle\rm{d}M=$ $\displaystyle\rm{d}N_{i}(t)-\rho$
$\displaystyle\mu_{aK}=$ $\displaystyle\mu_{a}\\{K_{i}(t);\alpha_{0K}\\}$
$\displaystyle\mu_{aV}=$ $\displaystyle\mu_{a}\\{V_{i}(t);\alpha_{0V}\\}$
$\displaystyle\mu_{0}=$ $\displaystyle E[Y^{0}_{i}(t)]=\beta_{0}$
$\displaystyle\mu_{1}=$ $\displaystyle E[Y^{1}_{i}(t)]=\beta_{0}+\beta_{1}.$

With the simpler notation, and under all nuisance models correctly specified,
the FIPTM (in this case, we do not consider modelling flexibly the intercept
as in the original paper proposing the FIPTM but rather consider a standard
doubly weighted least squares estimator for the FIPTM) and AAIIW estimators
respectively correspond to the solutions of the following sets of estimating
equations:

$\displaystyle
E_{n}\left[\int_{0}^{\tau}\frac{\rm{d}N}{\rho}\frac{I_{a}}{e_{a}}\left\\{Y_{i}(t)-\zeta_{i}(t;\beta_{a})\right\\}\right]=0$
(7)

and

$\displaystyle
E_{n}\left[\int_{0}^{\tau}\left\\{\frac{I_{a}}{e_{a}}Y_{i}(t){-\left(\frac{I_{a}-e_{a}}{e_{a}}\right)\mu_{aK}}-\zeta_{i}(t;\beta_{a})\right\\}\frac{\rm{d}N}{\rho}\right]-E_{n}\left[\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\int_{0}^{\tau}\frac{\rm{d}M}{\rho}\left\\{\frac{\mu_{aV}}{e_{a}}-\left(\frac{I_{a}-e_{a}}{e_{a}}\right)\mu_{aK}-\zeta_{i}(t;\beta_{a})\right\\}\right]$
$\displaystyle=0.$ (8)

Note that each set of estimating equations contains two equations, one for
$\beta_{0}$ and one for $\beta_{1}$ once $\beta_{0}$ is evaluated. We use the
influence function of each estimator to derive its asymptotic variance. The
variance of each estimator equals to the variance of its influence function.

Denote by $\beta_{a}=[\beta_{0}\hskip 2.84544pt\beta_{1}]$ the vector of true
parameters. We prove the relative efficiency of the AAIIW estimator when
compared with the FIPTM estimator under the assumption that all nuisance
models are correctly specified for both estimators. The influence function of
the $\hat{\beta}_{DW}$ estimator, denoted
$\psi(o_{i})=\psi(o_{i};\beta_{a},\delta_{0})$ for $o_{i}$ the observation
data of the ith individual, and
$\delta_{0}=\left\\{\psi_{0},\gamma_{0}\right\\}$ the true nuisance parameters
(for the treatment and observation models) must satisfy

$\sqrt{n}\left(\hat{\beta}_{DW}-\beta_{a}\right)=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\psi(o_{i})+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(1).$

By replacing terms in equation 7 above, we obtain a vector of influence
functions the same size as the vector of parameters of interest $\beta_{a}$,
that is:

$\psi(o_{i})=\left[\begin{matrix}\int_{t=0}^{\tau}\frac{\rm{d}N}{\rho}\frac{I_{0}}{e_{0}}\left\\{Y_{i}(t)-\mu_{0}\right\\}\\\
\int_{t=0}^{\tau}\frac{\rm{d}N}{\rho}\left\\{\left(\frac{I_{1}}{e_{1}}-\frac{I_{0}}{e_{0}}\right)Y_{i}(t)+\mu_{0}I_{0}-\mu_{1}I_{1}\right\\}\end{matrix}\right].$

For the variance of the DW estimator, our interest is in $\beta_{1}$ and so we
focus on the second equation (the only equation that depends on it), assuming
we already have an estimate for $\beta_{0}$ from the first equation. Using
similar developments as those that follow, we could show that the variance of
the estimate for $\beta_{0}$ is also smaller when using the AAIIW estimator as
compared with the FIPTM estimator, but this is omitted in what follows.

The variance of the FIPTM estimator equals to the variance of its influence
function. We denote that variance by $\sigma_{FIPTM}^{2}$. In further
derivations, we drop the integral that sums the terms over all times when
there is an observation of the outcome (as this integral sum is the same for
the FIPTM and the AAIIW estimator). We, therefore, focus on only one term in
the integral, yielding:

$\displaystyle\sigma_{FIPTM}^{2}=$ $\displaystyle
E[\psi(o_{i})^{2}]+E[\psi(o_{i})]^{2}$ $\displaystyle=$ $\displaystyle
E[\psi(o_{i})^{2}]+0$ $\displaystyle=$ $\displaystyle
E\left(\left[\frac{\rm{d}N}{\rho}\left\\{\left(\frac{I_{1}}{e_{1}}-\frac{I_{0}}{e_{0}}\right)Y_{i}(t)+\mu_{0}I_{0}-\mu_{1}I_{1}\right\\}\right]^{2}\right)$
$\displaystyle=$ $\displaystyle
E\left[\overbrace{\left\\{\frac{\rm{d}N}{\rho}\left(\frac{I_{1}}{e_{1}}\right)(Y_{i}(t)-\mu_{1})\right\\}^{2}}^{A}\right]+E\left[\overbrace{\left\\{\frac{\rm{d}N}{\rho}\left(\frac{I_{0}}{e_{0}}\right)(Y_{i}(t)-\mu_{0})\right\\}^{2}}^{B}\right]-2E\left[E\left\\{A\times
B\mid V_{i}(t)\right\\}\right]$ $\displaystyle=$ $\displaystyle
E\left[\frac{\rm{d}N}{\rho^{2}}\frac{I_{1}}{e_{1}^{2}}(Y_{i}(t)-\mu_{1})^{2}\right]+E\left[\frac{\rm{d}N}{\rho^{2}}\frac{I_{0}}{e_{0}^{2}}(Y_{i}(t)-\mu_{0})^{2}\right]-2E\left\\{E\left[A\mid
V_{i}(t)\right]\times E\left[B\mid V_{i}(t)\right]\right\\}$ with the third
term true since both terms are independent given V (treated and untreated
patients). Now, the third term is 0 because each expectation conditional on
$V_{i}(t)$ is zero. $\displaystyle=E\left[\frac{E[\rm{d}N\mid
V_{i}(t)]}{\rho^{2}}\frac{E[I_{1}\mid
V_{i}(t)]}{e_{1}^{2}}E\left[(Y_{i}(t)-\mu_{1})^{2}\mid
V_{i}(t),A_{i}(t)=1,\rm{d}N_{i}(t)=1\right]\right]$ $\displaystyle\hskip
11.38092pt+E\left[\frac{E[\rm{d}N\mid V_{i}(t)]}{\rho^{2}}\frac{E[I_{0}\mid
V_{i}(t)]}{e_{0}^{2}}E\left[(Y_{i}(t)-\mu_{0})^{2}\mid
V_{i}(t),A_{i}(t)=0,\rm{d}N_{i}(t)=1\right]\right]+0$
$\displaystyle=E\left[\frac{1}{\rho}\frac{1}{e_{1}}E\left[(Y_{i}(t)-\mu_{1})^{2}\mid
V_{i}(t),A_{i}(t)=1,\rm{d}N_{i}(t)=1\right]\right]$ $\displaystyle\hskip
11.38092pt+E\left[\frac{1}{\rho}\frac{1}{e_{0}}E\left[(Y_{i}(t)-\mu_{0})^{2}\mid
V_{i}(t),A_{i}(t)=0,\rm{d}N_{i}(t)=1\right]\right]$
$\displaystyle=E\left[\frac{(Y^{1}_{i}(t)-\mu_{1})^{2}}{\rho
e_{1}}\right]+E\left[\frac{(Y^{0}_{i}(t)-\mu_{0})^{2}}{\rho e_{0}}\right].$
(9)

For the AAIIW estimator, we need to show that the variance is smaller than
(9). First we derive the influence function
$\phi(o_{i})=\phi(o_{i};\beta_{a},\delta_{0})$ for $o_{i}$ the observation
data of the ith individual and $\delta_{0}$ now defined as the vector of true
parameters $\left\\{\psi_{0},\gamma_{0},\alpha_{K0},\alpha_{V0}\right\\}$
(which now includes the parameters of the two conditional outcome mean
models). We have that the influence function satisfies

$\sqrt{n}(\hat{\beta}_{AAIIW}-\beta_{a})=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\phi(o_{i})+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(1)$

or, using the formula for M-estimators provided in Tsiatis (Chapter 3, section
3.2) and noting that the AAIIW estimator is an M-estimator and is the solution
to

$\displaystyle\sum_{i=1}^{n}m(o_{i};\hat{\beta}_{AAIIW})=0$ (10)

with the m-function that can be derived directly from (10), then we can
compute the influence function by first computing the partial derivatives of
the $m$ function w. r. t. to the two parameters in $\beta_{a}$ (use, e.g.,
formula 3.6 in Tsiatis) and find that

$\displaystyle\phi(o_{i})=$
$\displaystyle\psi(o_{i})+\left[\begin{matrix}\int_{t=0}^{\tau}-\frac{\rm{d}N}{\rho}\left(\frac{I_{0}-e_{0}}{e_{0}}\right)\mu_{0K}-\frac{\rm{d}M}{\rho}\left\\{\frac{\mu_{0V}}{e_{0}}-\left(\frac{I_{0}-e_{0}}{e_{0}}\right)\mu_{0K}-\mu_{0}\right\\}\\\
\int_{t=0}^{\tau}-\frac{\rm{d}N}{\rho}\left(\frac{I_{1}-e_{1}}{e_{1}}\right)\mu_{1K}-\frac{\rm{d}M}{\rho}\left\\{\frac{\mu_{1V}}{e_{1}}-\left(\frac{I_{1}-e_{1}}{e_{1}}\right)\mu_{1K}-\mu_{1}\right\\}\\\
\hskip
28.45274pt+\frac{\rm{d}N}{\rho}\left(\frac{I_{0}-e_{0}}{e_{0}}\right)\mu_{0K}+\frac{\rm{d}M}{\rho}\left\\{\frac{\mu_{0V}}{e_{0}}-\left(\frac{I_{0}-e_{0}}{e_{0}}\right)\mu_{0K}-\mu_{0}\right\\}\end{matrix}\right]$

with the second part added to $\psi(o_{i})$ corresponding to the augmented
terms. Note that the second row of $\phi(o_{i})$ contains two terms, one being
the augmented term for the treated and one for the untreated. The change in
sign for both terms is due to the subtraction between the augmented terms due
to the treated and the untreated patients.

To compute the variance of the AAIIW, that we denote by $\sigma_{AAIIW}^{2}$,
we compute the variance of its influence function. We again focus on the
second influence function assuming that $\beta_{0}$ is already estimated, and
again we drop the integral sign in the following derivations. We have:

$\displaystyle\sigma_{AAIIW}^{2}=$ $\displaystyle E[\phi(o_{i})^{2}]$
$\displaystyle=$ $\displaystyle E\left[\psi(o_{i})^{2}\right]$
$\displaystyle\hskip
11.38092pt+E\left[\left(\overbrace{-\left\\{\frac{I_{1}-e_{1}}{e_{1}}\right\\}\mu_{1K}-\frac{\rm{d}M}{\rho}\left\\{\frac{\mu_{1V}}{e_{1}}-\mu_{1}\right\\}}^{C_{1}}+\overbrace{\left\\{\frac{I_{0}-e_{0}}{e_{0}}\right\\}\mu_{0K}+\frac{\rm{d}M}{\rho}\left\\{\frac{\mu_{0V}}{e_{0}}-\mu_{0}\right\\}}^{C_{2}}\right)^{2}\right]$
$\displaystyle\hskip 34.14322pt+2E[\psi(o_{i})\times(C_{1}+C_{2})].$

The first term, $E\left[\psi(o_{i})^{2}\right]$, is the variance of the FIPTM
estimator. Denote the second and third terms by $D=E[(C_{1}+C_{2})^{2}]$ and
$F=2E[\psi(o_{i})\times(C_{1}+C_{2})]$. We must show that $D+F<0$ to show
relative efficiency of the AAIIW estimator. We have that

$\displaystyle D=E[(C_{1}+C_{2})^{2}]=$ $\displaystyle
E[C_{1}^{2}+C_{2}^{2}+2C_{1}C_{2}]=E[C_{1}^{2}]+E[C_{2}^{2}]+E[2C_{1}C_{2}]$

where

$\displaystyle E[C_{1}^{2}]=$ $\displaystyle
E\left[\frac{\rm{d}M^{2}}{\rho^{2}}\left(\frac{\mu_{1V}}{e_{1}}-\mu_{1}\right)^{2}\right]+E\left[\left(\frac{I_{1}-e_{1}}{e_{1}}\right)^{2}\mu^{2}_{1K}\right]$
$\displaystyle-2E\left[E\left[-\frac{\rm{d}M}{\rho}\left\\{\frac{\mu_{1V}}{e_{1}}-\mu_{1}\right\\}\left(\frac{I_{1}-e_{1}}{e_{1}}\right)\mu_{1K}\mid
V_{i}(t)\right]\right]$ $\displaystyle=$ $\displaystyle
E\left[\frac{\rm{d}M^{2}}{\rho^{2}}\left(\frac{\mu_{1V}}{e_{1}}-\mu_{1}\right)^{2}\right]+E\left[\frac{I_{1}+e_{1}^{2}-2I_{1}e_{1}}{e_{1}^{2}}\mu^{2}_{1K}\right]$
$\displaystyle\text{where the third term was removed since}\hskip
5.69046ptE\left[\rm{d}M\mid V_{i}(t)\right]=0$ $\displaystyle=$ $\displaystyle
E\left[\frac{E[(\rm{d}N+\rho^{2}-2\rm{d}N\rho)\mid
V_{i}(t)]}{\rho^{2}}\left(E[Y^{1}_{i}(t)]-\mu_{1}\right)^{2}\
\right]+\mu_{1}^{2}E\left[\frac{E\left[I_{1}+e_{1}^{2}-2I_{1}e_{1}\mid
V_{i}(t)\right]}{e_{1}^{2}}\right]$ $\displaystyle=$ $\displaystyle
0+\mu_{1}^{2}E\left[\frac{1-e_{1}}{e_{1}}\right]\hskip
5.69046pt\text{since}\hskip 5.69046ptE[Y^{1}_{i}(t)]-\mu_{1}=0.$

Similarly, we have

$\displaystyle E[C_{2}^{2}]=$
$\displaystyle\mu_{0}^{2}E\left[\frac{1-e_{0}}{e_{0}}\right].$

The third term, $E[2C_{1}C_{2}]$, cancels out, after using iterated
expectation and conditioning on e.g., $V_{i}(t)$, as the martingale residuals
and the treatment residuals are zero-mean conditional on that set. Thus,

$\displaystyle D=$
$\displaystyle\mu_{1}^{2}E\left[\frac{1-e_{1}}{e_{1}}\right]+\mu_{0}^{2}E\left[\frac{1-e_{0}}{e_{0}}\right].$

We are left with calculting $F$. We have:

$\displaystyle F=$ $\displaystyle 2E[\psi(o_{i})\times(C_{1}+C_{2})]$
$\displaystyle=$ $\displaystyle 2E[\psi(o_{i})\times
C_{1}]+2E[\psi(o_{i})\times C_{2}].$

We start with the left term and have

$\displaystyle E[\psi(o_{i})\times C_{1}]=$ $\displaystyle
E\left[\frac{\rm{d}N}{\rho}\left\\{\left(\frac{I_{1}}{e_{1}}-\frac{I_{0}}{e_{0}}\right)Y_{i}(t)+\mu_{0}I_{0}-\mu_{1}I_{1}\right\\}\times\left(-\left\\{\frac{I_{1}-e_{1}}{e_{1}}\right\\}\mu_{1K}-\frac{\rm{d}M}{\rho}\left\\{\frac{\mu_{1V}}{e_{1}}-\mu_{1}\right\\}\right)\right]$
$\displaystyle=$
$\displaystyle-E\left[\frac{\rm{d}N}{\rho}\left\\{\left(\frac{I_{1}}{e_{1}}-\frac{I_{0}}{e_{0}}\right)Y_{i}(t)+\mu_{0}I_{0}-\mu_{1}I_{1}\right\\}\times\left(\left\\{\frac{I_{1}-e_{1}}{e_{1}}\right\\}\mu_{1K}+\frac{\rm{d}M}{\rho}\left\\{\frac{\mu_{1V}}{e_{1}}-\mu_{1}\right\\}\right)\right]$
$\displaystyle=$
$\displaystyle-E\left[\frac{\rm{d}N}{\rho}\left\\{\frac{I_{1}}{e_{1}}Y_{i}(t)-\mu_{1}\right\\}\times\left\\{\frac{I_{1}-e_{1}}{e_{1}}\right\\}\mu_{1K}\right]+E\left[\frac{\rm{d}N}{\rho}\left\\{\frac{I_{0}}{e_{0}}Y_{i}(t)-\mu_{0}\right\\}\times\left\\{\frac{I_{1}-e_{1}}{e_{1}}\right\\}\mu_{1K}\right]$
$\displaystyle\hskip
7.11317pt-E\left[\frac{\rm{d}N}{\rho}\left\\{\frac{I_{1}}{e_{1}}Y_{i}(t)-\mu_{1}\right\\}\times\frac{\rm{d}M}{\rho}\left\\{\frac{\mu_{1V}}{e_{1}}-\mu_{1}\right\\}\right]+E\left[\frac{\rm{d}N}{\rho}\left\\{\frac{I_{0}}{e_{0}}Y_{i}(t)-\mu_{0}\right\\}\times\frac{\rm{d}M}{\rho}\left\\{\frac{\mu_{1V}}{e_{1}}-\mu_{1}\right\\}\right]$
$\displaystyle=-E\left[\frac{\rm{d}N}{\rho}\left\\{\frac{I_{1}-e_{1}I_{1}}{e_{1}^{2}}Y_{i}(t)-\left\\{\frac{I_{1}-e_{1}}{e_{1}}\right\\}\mu_{1}\right\\}\times\mu_{1K}\right]$
$\displaystyle\hskip
14.22636pt+E\left[\frac{\rm{d}N}{\rho}\left\\{\frac{I_{0}(I_{1}-e_{1})}{e_{0}e_{1}}Y_{i}(t)-\left\\{\frac{I_{1}-e_{1}}{e_{1}}\right\\}\mu_{0}\right\\}\times\mu_{1K}\right]$
$\displaystyle\hskip
28.45274pt-E\left[\frac{\rm{d}N-\rm{d}N\rho}{\rho^{2}}\left\\{\frac{I_{1}}{e_{1}}Y_{i}(t)-\mu_{1}\right\\}\times\left\\{\frac{\mu_{1V}}{e_{1}}-\mu_{1}\right\\}\right]$
$\displaystyle\hskip
42.67912pt+E\left[\frac{\rm{d}N-\rm{d}N\rho}{\rho^{2}}\left\\{\frac{I_{0}}{e_{0}}Y_{i}(t)-\mu_{0}\right\\}\times\left\\{\frac{\mu_{1V}}{e_{1}}-\mu_{1}\right\\}\right].$

Taking iterated expectation and conditioning on $V_{i}(t)$, several residuals
cancel out and we are left with

$\displaystyle=-E\left[\frac{1-e_{1}}{e_{1}}\mu_{1}^{2}\right]-E\left[\mu_{1}^{2}\right]-0+0$
$\displaystyle=-\mu_{1}^{2}E\left[\frac{1}{e_{1}}\right]$

Similarly, we find

$\displaystyle E[\psi(o_{i})\times
C_{2}]=-\mu_{0}^{2}E\left[\frac{1}{e_{0}}\right].$

The variance of the AAIIW is therefore obtained by summing:

$\displaystyle\sigma^{2}_{AAIIW}$ $\displaystyle=\sigma^{2}_{FIPTM}+D+F$
$\displaystyle=\sigma^{2}_{FIPTM}+\mu_{1}^{2}E\left[\frac{1-e_{1}}{e_{1}}\right]+\mu_{0}^{2}E\left[\frac{1-e_{0}}{e_{0}}\right]-2\mu_{1}^{2}E\left[\frac{1}{e_{1}}\right]-2\mu_{0}^{2}E\left[\frac{1}{e_{0}}\right]$
$\displaystyle\hskip 56.9055pt\iff$
$\displaystyle\sigma^{2}_{AAIIW}-\sigma^{2}_{FIPTM}$
$\displaystyle=\mu_{1}^{2}E\left[\frac{1-e_{1}}{e_{1}}-\frac{2}{e_{1}}\right]+\mu_{0}^{2}E\left[\frac{1-e_{0}}{e_{0}}-\frac{2}{e_{0}}\right].$

And, since $\mu_{1}^{2},\mu_{0}^{2}>0$, and

$\frac{1-e_{1}}{e_{1}}-\frac{2}{e_{1}}=\frac{-1-e_{1}}{e_{1}}\leq\frac{-1}{e_{1}}<0\hskip
5.69046pt\text{for all}\hskip 5.69046pte_{1}\hskip 5.69046pt\text{s. t.}\hskip
5.69046pt0\leq e_{1}\leq 1$

and

$\frac{1-e_{0}}{e_{0}}-\frac{2}{e_{0}}=\frac{-1-e_{0}}{e_{0}}\leq\frac{-1}{e_{0}}<0\hskip
5.69046pt\text{for all}\hskip 5.69046pte_{0}\hskip 5.69046pt\text{s. t.}\hskip
5.69046pt0\leq e_{0}\leq 1$

we find that

$\sigma^{2}_{AAIIW}-\sigma^{2}_{FIPTM}<0$
$\iff\sigma^{2}_{AAIIW}<\sigma^{2}_{FIPTM}$

under all nuisance models correctly specified for both estimators.

Supplementary Material C

Simulation studies details

Data were simulated to emulate a setting in which time is continuous and
covariates can be measured or updated at any point in time. For that, we used
a time grid starting at time 0 spanning up to time $\tau=2$ over which the
variables were simulated at each 0.1-width time bin. First, the confounders,
the treatment, the mediator, the pure predictor and the outcome were simulated
over each of these small time bins, after which the observation process was
simulated (and outcomes were removed). The individual index is omitted in what
follows.

For each individual, three baseline confounders were simulated as $K_{1}\sim
N(1,1)$,
$K_{2}\sim\text{Bernoulli}(0.55)$ and $K_{3}\sim N(0,1)$ and were repeated
throughout follow-up time. At each time $t$, a time-varying binary treatment
was simulated as $A(t)\sim\text{Bernoulli}(p_{t})$ with
$p_{t}=\text{expit}(-0.5+0.8\hskip 2.84544ptK_{1}-0.4\hskip
2.84544ptK_{2}-0.4\hskip 2.84544ptK_{3})$ where
expit$(j)=\exp(j)/\left\\{1+\exp(j)\right\\}$ and the confounders are kept
time-fixed (i.e., we emulate a setting in which they are measured at
baseline). At each time $t$, a mediator $M(t)$ of the causal effect of $A(t)$
on the outcome $Y(t)$ was further simulated as $M(t)\mid A(t)=1\sim N(2,1)$
and $M(t)\mid A(t)=0\sim N(4,2)$ where the parameters in the Normal
respectively correspond to the mean and variance of the random variable. A
time-varying pure predictor was simulated as $P(t)\sim N(0.5,0.09)$.

A few details are particularly important in the simulation of the outcome
$Y(t)$. First, since we would like to assess the performance of the different
approaches when a mediator affects the outcome observation process, the
outcome must depend on the mediator $M(t)$ and the mediator must affect the
outcome observation. However, we would like to know the gold standard for the
causal effect in the simulations, i.e., the true causal marginal effect of
treatment once we marginalize the outcome mean model over $M(t)$. For that, we
first modeled $E[M(t)\mid A(t),K_{1},K_{2},K_{3}]$ using a least squares
linear model. Then, we simulated the outcome as

$Y(t)=\kappa+1\hskip 4.26773ptA(t)+0.4\hskip 4.26773ptK_{1}+0.05\hskip
4.26773ptK_{2}-0.6\hskip 4.26773ptK_{3}+$ $\hskip
56.9055pt3\left\\{M(t)-E[M(t)\mid A(t),K_{1},K_{2},K_{3}]\right\\}+0.3\hskip
4.26773ptP(t)+\epsilon(t)$

where $\kappa=0.5$, $\epsilon(t)\sim N(\phi,0.01)$, and $\phi\sim N(0,0.04)$.
We let $\phi$ vary by individual. When marginalizing over $M(t)$, we obtain a
true marginal causal effect of 1 for treatment $A(t)$. This would not be the
case if we had interactions between the treatment $A(t)$ and the confounders
$K$ in the outcome generating mechanism above.

The outcome above was first simulated once for each point on the grid. The
outcome values were then set to missing according to a covariate-dependent
observation process with an observation rate (Poisson model) or probability
(Bernoulli model) depending on the treatment, the mediator, the pure
predictor, and the confounders. For each time point, this was done by
simulating an observation indicator (one minus a missingness indicator). In a
first set of simulation studies, we used a non-homogeneous Poisson process
with rate denoted by

$\lambda(t\mid
A(t),M(t),P(t),K;\gamma^{\dagger})=0.25(t+0.05)\text{exp}\left\\{\gamma_{8}A(t)+\gamma_{9}M(t)+\gamma_{10}K_{1}+\gamma_{11}K_{2}+\gamma_{12}K_{3}+\gamma_{13}P(t)\right\\}.$

We tested the four following combinations for the parameters
$\gamma^{\dagger}$: $1)\gamma^{\dagger}=(0,0,0,0,0,-5)$ (i.e., no bias due to
the visit process expected); 2) $\gamma^{\dagger}=(0.5,0.3,-0.5,-2,0,-3)$; 3)
$\gamma^{\dagger}=(0.5,-0.5,-0.2,-1,1,-3)$; and 4)
$\gamma^{\dagger}=(-1,-0.8,0.1,0.3,-1,-3)$. In a second set of simulation
studies, we used Bernoulli random variables with probabilities denoted by
$\lambda(t\mid A(t),M(t),P(t),K;\gamma)$
$=\text{expit}\left\\{\gamma_{1}+\gamma_{2}A(t)+\gamma_{3}M(t)+\gamma_{4}K_{1}+\gamma_{5}K_{2}+\gamma_{6}K_{3}+\gamma_{7}P(t)\right\\}$
to simulate the indicators. We tested the four following combinations for the
parameters $\gamma$: $1)\gamma=(0.4,0,0,0,0,0,-5)$ (i.e., no bias due to the
visit process expected); 2) $\gamma=(0.4,1,-1,-0.5,-2,0,-3)$;
3) $\gamma=(0.4,0.5,-0.5,-0.2,-1,1,-3)$; and 4)
$\gamma=(0.4,-0.5,0.8,0.1,0.3,-1,-3)$. In both sets of simulations, the
observation indicators could be obtained by simulating Bernoulli random
variables with probability proportional to the Poisson rate or the Bernoulli
probability. When, at a given time point, the simulated observation indicator
was simulated as 0, the corresponding outcome value was set to missing.

The ordinary least squares estimator, the inverse probability of treatment-
weighted estimators and the FIPTM estimators were obtained by fitting standard
linear regressions (e.g., with the function lm in R) using the weight
statement whenever applicable for the inverse probability of treatment or the
inverse intensity of visit weights. The conditional outcome mean models were
obtained by fitting two different linear models, one conditioning on the
confounders and the treatment (corresponding to
$\mu_{a}\\{K,A=a;\hat{\alpha}_{K}\\}$ in the main manuscript) and one
conditional on the visit predictors and the treatment (corresponding to
$\mu_{a}\\{V,A=a;\hat{\alpha}_{V}\\}$).

To assess the performance of the AAIIW estimator when the model
$\mu_{a}\\{K,A=a;\alpha_{K}\\}$ is correctly specified, we had to know the
true model $\mu_{a}\\{K,A=a;\alpha_{K0}\\}$. Based on the outcome generating
mechanism, that true model is

$\displaystyle E[Y(t)\mid A(t),K]$ $\displaystyle=E[\kappa+1\hskip
4.26773ptA(t)+0.4\hskip 4.26773ptK_{1}+0.05\hskip 4.26773ptK_{2}-0.6\hskip
4.26773ptK_{3}$ $\displaystyle\hskip 28.45274pt+3\left\\{M(t)-E[M(t)\mid
A(t),K_{1},K_{2},K_{3}]\right\\}+0.3\hskip 4.26773ptP(t)+\epsilon(t)\mid
A(t),K]$ $\displaystyle=\kappa+A(t)+0.4K_{1}+0.05K_{2}-0.6K_{3}+0.3P(t).$

The variable $P(t)$ is not mandatory to adjust for in our setting, since it
should not affect the contrast in the outcome mean across treated ($A(t)=1$)
and untreated patients ($A(t)=0$) (unless there was an interaction term
between $P(t)$ and the treatment in the outcome generating mechanism). The
correctly specified model $\mu_{a}\\{V,A=a;\alpha_{V}\\}$ merely contains all
the predictors that were used in generating the outcome above (i.e., the
treatment, the re-centered mediator, the pure predictor, and all the
confounders).

Supplementary Material D

Additional results of the simulation studies: sample of size 250 with the use
of a nonhomogeneous Poisson rate, or the use of a Bernoulli probability
(sample sizes 250 and 1000) to simulate the observation indicators and Monte
Carlo empirical bias and mean square error (MSE) in the simulation studies

Suppl. Fig. 1: Results of the simulation studies with a sample size of 250
using a nonhomogeneous Poisson rate to simulate the observation indicators and
the Andersen and Gill model with Breslow estimator to estimate the IIV
weights. Each boxplot represents the distribution of 1000 estimates for the
corresponding estimator. The dashed line represents the gold standard, i.e.,
the true value for the marginal effect of exposure that equals to 1. Different
strengths of the visit process on covariates are represented with scenarios 1)
$\gamma=(0,0,0,0,0,-5)$ (i.e., no bias due to the visit process expected); 2)
$\gamma=(0.5,0.3,-0.5,-2,0,-3)$; 3) $\gamma=(0.5,-0.5,-0.2,-1,1,-3)$; and 4)
$\gamma=(-1,-0.8,0.1,0.3,-1,-3)$. Suppl. Fig. 2: Results of the simulation
studies using Binomial probabilities with a sample size of 1000. Each boxplot
represents the distribution of 1000 estimates for the corresponding estimator.
The dashed line represents the gold standard, i.e., the true value for the
marginal effect of exposure that equals to 1. Different strengths of the visit
process on covariates are represented with scenarios 1)
$\gamma=(0.4,0,0,0,0,0,-5)$ (i.e., no bias due to the visit process expected);
2) $\gamma=(0.4,1,-1,-0.5,-2,0,-3)$; 3) $\gamma=(0.4,0.5,-0.5,-0.2,-1,1,-3)$;
and 4) $\gamma=(0.4,-0.5,0.8,0.1,0.3,-1,-3)$. Suppl. Fig. 3: Results of the
simulation studies using Binomial probabilities with a sample size of 250.
Each boxplot represents the distribution of 1000 estimates for the
corresponding estimator. The dashed line represents the gold standard, i.e.,
the true value for the marginal effect of exposure that equals to 1. Different
strengths of the visit process on covariates are represented with scenarios 1)
$\gamma=(0.4,0,0,0,0,0,-5)$ (i.e., no bias due to the visit process expected);
2) $\gamma=(0.4,1,-1,-0.5,-2,0,-3)$; 3) $\gamma=(0.4,0.5,-0.5,-0.2,-1,1,-3)$;
and 4) $\gamma=(0.4,-0.5,0.8,0.1,0.3,-1,-3)$. Suppl. Table 1: Simulations
results, empirical bias and mean squared error (MSE)

| Average $N(\tau)$ | Gamma set | Poisson rate | Bernoulli probability
---|---|---|---|---
Estimator | $(A_{i}(t)=0,1)$ | $\gamma$ | $n=1000$ | $n=250$ | $n=1000$ | $n=250$
| Poisson $\dagger$ | Poisson $\dagger$ | Bias | MSE | Bias | MSE | Bias | MSE | Bias | MSE
| or Bernoulli $\ddagger$ | or Bernoulli $\ddagger$ |  |  |  |  |  |  |  |
OLS | (12, 12)† | 1† | 0.47 | 0.22 | 0.47 | 0.24 | 0.47 | 0.22 | 0.46 | 0.22
IPTc | (18, 18)‡ | 1‡ |  $<$0.01 | 0.01 | 0.01 | 0.03 |  $<$0.01 |  $<$0.01 | 0.01 | 0.01
IPTnc |  |  | 0.23 | 0.06 | 0.23 | 0.08 | 0.23 | 0.05 | 0.22 | 0.05
DWc |  |  | 0.01 | 0.13 | 0.01 | 0.44 |  $<$0.01 | 0.01 |  $<$0.01 | 0.03
DWiptc |  |  |  $<$0.01 | 0.01 | 0.01 | 0.03 |  $<$0.01 |  $<$0.01 | 0.01 | 0.01
DWiivc |  |  | 0.23 | 0.18 | 0.22 | 0.47 | 0.23 | 0.06 | 0.22 | 0.08
DWnc |  |  | 0.23 | 0.06 | 0.23 | 0.08 | 0.23 | 0.05 | 0.22 | 0.05
AAIIWc |  |  |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$ 0.01 |  $<$ 0.01
AAIIWs.a. |  |  | 0.01 | 0.15 |  $<$0.01 | 0.73 |  $<$0.01 | 0.01 |  $<$0.01 | 0.03
AAIIWs.b. |  |  |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01
AAIIWs.c. |  |  |  $<$0.01 | 0.14 |  $<$0.01 | 0.84 |  $<$0.01 | 0.01 |  $<$0.01 | 0.03
AAIIWs.d. |  |  |  $<$0.01 | 0.10 |  $<$0.01 | 0.46 |  $<$0.01 | 0.01 |  $<$0.01 | 0.02
OLS | (22, 17)† | 2† | 0.22 | 0.05 | 0.23 | 0.07 | 1.30 | 1.68 | 1.29 | 1.67
IPTc | (1, 7)‡ | 2‡ | 0.75 | 0.56 | 0.75 | 0.57 | 0.93 | 0.86 | 0.92 | 0.87
IPTnc |  |  | 0.52 | 0.27 | 0.52 | 0.29 | 1.13 | 1.27 | 1.12 | 1.27
DWc |  |  | 0.05 | 0.08 | 0.04 | 0.29 |  $<$0.01 | 0.02 |  $<$0.01 | 0.06
DWiptc |  |  | 0.82 | 0.68 | 0.82 | 0.71 | 0.94 | 0.89 | 0.93 | 0.88
DWiivc |  |  | 0.27 | 0.15 | 0.26 | 0.37 | 0.22 | 0.06 | 0.22 | 0.10
DWnc |  |  | 0.61 | 0.38 | 0.61 | 0.41 | 1.14 | 1.31 | 1.13 | 1.30
AAIIWc |  |  |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01
AAIIWs.a. |  |  | 0.01 | 0.10 | 0.02 | 0.49 |  $<$0.01 | 0.02 |  $<$0.01 | 0.07
AAIIWs.b. |  |  |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01
AAIIWs.c. |  |  | 0.03 | 0.10 | 0.04 | 0.51 |  $<$0.01 | 0.02 |  $<$0.01 | 0.06
AAIIWs.d. |  |  | 0.05 | 0.08 | 0.05 | 0.42 |  $<$0.01 | 0.01 |  $<$0.01 | 0.05
OLS | (3, 8)† | 3† | 2.02 | 4.09 | 2.00 | 4.08 | 1.85 | 3.43 | 1.85 | 3.43
IPTc | (6, 14)‡ | 3‡ | 1.62 | 2.63 | 1.61 | 2.65 | 1.51 | 2.27 | 1.51 | 2.28
IPTnc |  |  | 1.79 | 3.22 | 1.78 | 3.23 | 1.64 | 2.69 | 1.64 | 2.70
DWc |  |  | 0.01 | 0.67 | 0.18 | 1.63 | 0.01 | 0.09 | 0.04 | 0.29
DWiptc |  |  | 1.56 | 2.59 | 1.55 | 2.62 | 1.55 | 2.41 | 1.55 | 2.44
DWiivc |  |  | 0.21 | 0.57 | 0.34 | 1.55 | 0.23 | 0.12 | 0.24 | 0.28
DWnc |  |  | 1.77 | 3.17 | 1.75 | 3.24 | 1.70 | 2.91 | 1.70 | 2.92
AAIIWc |  |  |  $<$0.01 |  $<$0.01 |  $<$0.01 | 0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01
AAIIWs.a. |  |  | 0.38 | 2.68 | 0.31 | 11.13 |  $<$0.01 | 0.14 |  $<$0.01 | 0.91
AAIIWs.b. |  |  |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01
AAIIWs.c. |  |  | 0.34 | 1.83 | 0.28 | 9.62 |  $<$0.01 | 0.10 | 0.01 | 0.42
AAIIWs.d. |  |  | 0.45 | 1.80 | 0.43 | 8.67 |  $<$0.01 | 0.09 | 0.01 | 0.37
OLS | (2, 5)† | 4† | 2.78 | 7.77 | 2.78 | 7.80 | 0.62 | 0.38 | 0.61 | 0.38
IPTc | (81, 58)‡ | 4‡ | 2.26 | 5.12 | 2.25 | 5.19 | 0.28 | 0.08 | 0.28 | 0.08
IPTnc |  |  | 2.54 | 6.45 | 2.53 | 6.48 | 0.38 | 0.15 | 0.38 | 0.14
DWc |  |  | 0.01 | 0.84 | 0.10 | 1.93 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01
DWiptc |  |  | 2.35 | 5.57 | 2.35 | 5.79 | 0.28 | 0.08 | 0.28 | 0.08
DWiivc |  |  | 0.25 | 0.91 | 0.37 | 2.08 | 0.23 | 0.05 | 0.22 | 0.05
DWnc |  |  | 2.59 | 6.78 | 2.59 | 6.98 | 0.37 | 0.14 | 0.37 | 0.14
AAIIWc |  |  |  $<$0.01 |  $<$0.01 |  $<$0.01 | 0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01
AAIIWs.a. |  |  | 0.34 | 4.77 | 0.21 | 13.54 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01
AAIIWs.c. |  |  | 0.38 | 6.53 | 0.25 | 14.89 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01
AAIIWs.d. |  |  | 0.64 | 6.33 | 0.41 | 12.58 |  $<$0.01 |  $<$0.01 |  $<$0.01 |  $<$0.01

$\dagger$. 1: $\gamma$=(0, 0, 0, 0, 0, -5); 2: $\gamma$=(0.5, 0.3, -0.5, -2,
0, -3); 3: $\gamma$=(0.5, -0.5, -0.2, -1, 1, -3); 4: $\gamma$=(-1, -0.8, 0.1,
0.3, -1, -3).
$\ddagger$. 1: $\gamma$=(0.4, 0, 0, 0, 0, 0, -5); 2: $\gamma$=(0.4, 1, -1,
-0.5, -2, 0, -3); 3: $\gamma$=(0.4, 0.5, -0.5, -0.2, -1, 1, -3); 4:
$\gamma$=(0.4, -0.5, 0.8, 0.1, 0.3, -1, -3).

Supplementary Material E
Tables of characteristics in the Add Health study stratified by weighting
strategy

Suppl. Table 2: Longitudinal characteristics stratified by adolescents
receiving counselling therapy or not in the first dataset imputed with
multiple imputations by chained equations, Add Health study, United States,
1996-2008

| Before inverse probability | After inverse probability
---|---|---
| of treatment-weighting | of treatment-weighting
Variable, N (%) | Counselling | No counselling | Counselling | No counselling
Age, mean (SD) | 20.3 (5.8) | 20.8 (5.7) | 20.9 (6.0) | 20.8 (5.7)
Female sex | 1579 (60.8) | 11,848 (50.6) | 13,528 (51.8) | 13,427 (51.6)
Weight, mean (SD) | 156.2 (42.6) | 160.7 (44.4) | 160.6 (44.4) | 160.3 (44.3)
Socioeconomic status, mean (SD) | 6.5 (2.2) | 6.4 (2.2) | 6.4 (2.2) | 6.5 (2.2)
Smoking | 1090 (42.0) | 6476 (27.7) | 7451 (28.6) | 7563 (29.1)
Depressive mood, mean (SD)† | 1.9 (0.9) | 1.4 (0.7) | 1.8 (0.9) | 1.4 (0.7)

$\dagger$. This is not considered as a confounder but as a mediator in our
analyses. Therefore, depressive mood was not included in the inverse
probability of treatment weights. Acronym: SD, standard deviation.

Suppl. Table 3: Longitudinal characteristics stratified by the alcohol
consumption being observed or not in the first dataset imputed with multiple
imputations by chained equations and after inducing outcome missingness, Add
Health study, United States, 1996-2008

| Before inverse intensity | After inverse intensity
---|---|---
| of visit-weighting | of visit-weighting
Variable, N (%) | Observed | Not observed | Observed | Not observed
Age, mean (SD) | 20.9 (5.7) | 20.7 (5.7) | 20.9 (5.7) | 21.0 (5.7)
Female sex | 5283 (50.0) | 8144 (52.7) | 2777 (50.2) | 2496 (50.0)
Weight, mean (SD) | 160.3 (44.3) | 160.2 (45.0) | 159.9 (44.1) | 161.7 (45.0)
Socioeconomic status, mean (SD) | 6.6 (2.1) | 6.5 (2.1) | 6.6 (2.1) | 6.6 (2.1)
Smoking | 3159 (29.9) | 4912 (31.8) | 1686 (30.5) | 1412 (28.0)
Depressive mood, mean (SD) | 1.1 (0.4) | 1.6 (0.8) | 1.1 (0.3) | 1.2 (0.4)
Counselling | 0.2 (0.4) | 0.1 (0.2) | 0.3 (0.4) | 0.1 (0.3)

Acronym: SD, standard deviation.