# A causal approach to functional mediation analysis with application to a
smoking cessation intervention

Donna L. Coffman∗ PhD John J. Dziak PhD Kaylee Litson PhD Yajnaseni
Chakraborti MS Megan E. Piper PhD Runze Li PhD CONTACT D. L. Coffman. Email:
<EMAIL_ADDRESS>Department of Epidemiology and Biostatistics, Temple
University, Philadelphia, PA, 19122, USA  Edna Bennett Pierce Prevention
Research Center, The Pennsylvania State University, University Park, PA,
16802, USA  Instructional Technology & Learning Sciences Department, Utah
State University, Logan, UT, USA  Department of Epidemiology and
Biostatistics, Temple University, Philadelphia, PA, 19122, USA  Center for
Tobacco Research Intervention, University of Wisconsin, Madison, WI, USA
Department of Statistics, The Pennsylvania State University, University Park,
PA, 16802, USA

###### Abstract

The increase in the use of mobile and wearable devices now allows dense
assessment of mediating processes over time. For example, a pharmacological
intervention may have an effect on smoking cessation via reductions in
momentary withdrawal symptoms. We define and identify the causal direct and
indirect effects in terms of potential outcomes on the mean difference and
odds ratio scales, and present a method for estimating and testing the
indirect effect of a randomized treatment on a distal binary variable as
mediated by the nonparametric trajectory of an intensively measured
longitudinal variable (e.g., from ecological momentary assessment). Coverage
of a bootstrap test for the indirect effect is demonstrated via simulation. An
empirical example is presented based on estimating later smoking abstinence
from patterns of craving during smoking cessation treatment. We provide an R
package, `funmediation`, available on CRAN at
https://cran.r-project.org/web/packages/funmediation/index.html, to
conveniently apply this technique. We conclude by discussing possible
extensions to multiple mediators and directions for future research.

###### keywords:

Ecological momentary assessment, functional regression, indirect effect,
mediation, smoking, time-varying effect modeling

This article has been accepted for publication in Multivariate Behavioral
Research, published by Taylor & Francis. Dr. Coffman is now at the University
of South Carolina.

##  1 Introduction

Increasingly, researchers collect intensive longitudinal data on variables
which are hypothesized to mediate the effect of an intervention on an end-of-
study outcome of interest. For example, in a smoking cessation intervention
study, individuals may be randomly assigned to a smoking cessation treatment
such as the nicotine patch or varenicline. The first several weeks post-quit
day are critical to maintaining smoking abstinence. Withdrawal symptoms (e.g.,
cravings, negative affect) may fluctuate even within a single day during these
early weeks, and these fluctuations may have an effect on maintaining smoking
abstinence. An effective smoking cessation treatment may have an effect on
smoking abstinence that is due in part to changes in these mediating variables
(i.e., due to a reduction in withdrawal symptoms).

Intensive longitudinal data can be collected using ecological momentary
assessment (EMA; see Shiffman ., 2008), in which individuals are contacted
multiple times per day for several weeks, and asked to report their cravings,
affect, and other variables, including use behavior. The data need not be
collected using EMA prompts; for example, it may be collected automatically
with a mobile or wearable device. EMA data are also collected in other
substance use cessation studies, and in weight loss and physical activity
studies. In these studies, individuals may be assigned to a behavioral
intervention, mediators are measured intensively over time, and the ultimate
outcome of interest is whether the behavior change in question has
successfully been maintained. The goal of this manuscript is to present and
estimate a mediation model to evaluate the process by which a randomized
intervention affects an end-of-study outcome via an intensively measured
longitudinal process.

Over the last several decades, mediation analysis has been extended from
single time point measures of a treatment, $X$, mediator, $M$, and outcome,
$Y$, to multiple mediators, and to treatment, mediators, and outcomes that are
measured longitudinally (MacKinnon, 2008). More recently, the assumptions
needed to define and identify causal effects in the context of mediation have
received considerable attention. As mediation is inherently about causal
pathways, in which a treatment affects change in a mediator, which in turn
affects change in the outcome, the potential outcomes framework has been
employed to define and identify causal mediation effects (VanderWeele, 2015).

Only in the last decade has mediation with an intensively measured treatment,
mediator, or outcome within the potential outcomes framework been introduced.
Specifically, Lindquist (2012) proposed an approach for a mediator measured
intensively by functional magnetic resonance imaging (fMRI). Other
advancements have included time-varying treatment and time-varying mediators
(VanderWeele Tchetgen Tchetgen, 2017; Park ., 2018; Huang Yuan, 2017), time-
varying mediators and outcomes (Zeng ., 2021; Loh ., 2020; Huang Yuan, 2017),
time-varying treatment and mediator with a survival outcome (Vansteelandt .,
2019; Zheng van der Laan, 2017), and an exogenous (e.g., environmental) time-
varying treatment with time-varying mediator and outcome (Bind ., 2016).
Additionally, Qian . (2019) study the mediation of a short-term effect of a
time-varying treatment (see also K. Liu ., 2020). However, here we focus on a
single-time distal effect of a single-time treatment, as mediated by a process
occurring between those times.

Specifically, we consider a time-invariant treatment, intensively measured
mediator, and end-of-study outcome (see Figure 1). This structure most closely
resembles that of Lindquist (2012), but he considered only a continuous
normally distributed outcome. We propose an extension to a binary outcome,
such as smoking abstinence at follow-up in our motivating study. A binary
outcome in a causal mediation analysis can be somewhat more complicated, even
for a mediator measured at only a single time point (Rijnhart ., 2017, 2019,
2021). Causal effects for numerical outcomes are generally stated as
differences in expected values. For binary outcomes, the difference in
expected values is interpreted as a difference in probabilities (i.e., risk
difference). However, if a risk ratio, odds ratio, or other quantity is
desired, then causal mediation effects are defined differently (VanderWeele
Vansteelandt, 2010).

Although VanderWeele Tchetgen Tchetgen (2017) consider a single time point
measurement of a binary outcome and a time-varying mediator, their approach
only partially applies in our context because we, like Lindquist (2012) and
Zeng . (2021), view the mediation process as functional. One approach to
fitting a mediation model in a setting with intensive longitudinal
measurements is to replace a longitudinal process, $M(t)$ and/or $Y(t)$, with
summary statistics such as means, sums, or estimated linear slopes. This
parametric summary approach could be made sophisticated by treating intercepts
and slopes as subject-level latent random variables (Magidson ., 2015; de
Stavola ., 2022; Sullivan ., 2021). However, as Lindquist (2012) pointed out,
it is also possible to take a less parametric approach, avoiding imposing
linear shapes on trajectories, using ideas from functional regression (Ramsay
Silverman, 2005). Although Zeng . (2021) propose functional mediation for an
intensively measured mediator and intensively measured outcome, they do not
propose a method for evaluating an intensively measured mediator with an end
of study outcome. The functional mediation model proposed in this paper uses
an outcome measured at a single point in time and thus, some of the
assumptions made by Zeng . (2021) are not necessary in the present context.
Nevertheless, Lindquist (2012), VanderWeele Tchetgen Tchetgen (2017), and Zeng
. (2021) provide foundational ideas to the approach proposed here.

This article is organized as follows. In Section 2, we present our motivating
example. In Section 3, we introduce the notation, definitions, and
identification assumptions needed in order for the effects to be interpreted
as causal. In Section 4, we present estimation of the indirect effect in
detail, which we have implemented in an R package, funmediation. In Section 5,
we present a simulation study to assess the coverage and power in data
simulated to resemble to our real-data example. In Section 6, we analyze the
real data using the proposed model. Finally, in Section 7, we discuss
limitations, future directions, and conclusions. In Appendix 1, we review the
proof by Lindquist (2012) for the identity link function and extend it to a
log link function modeling a probability ratio. The more commonly used
logistic link can be approximated by a log link in the case of rare outcomes
(see Feingold ., 2019; VanderWeele Vansteelandt, 2010; Rijnhart ., 2021).

##  2 Motivating Example

Our motivating example comes from the Wisconsin Smokers’ Health Study 2
(WSHS2; Baker ., 2016, ClinicalTrials.gov identifier: NCT01553084). The WSHS2
is a comparative efficacy trial directly comparing the two most effective
smoking cessation therapies (varenicline and combination nicotine replacement
therapy [cNRT]) with one another and with an active comparator treatment
(nicotine patch) (see Baker ., 2016; Le Foll ., 2022).

Smokers ($n=1086$; 52% women, 67% white, mean age = 48, mean number of
cigarettes smoked/day = 17), recruited from Madison and Milwaukee, were
randomly assigned to one of the three pharmacotherapies. The outcome
considered here is carbon-monoxide (CO)-confirmed self-reported 7-day point-
prevalence abstinence at 4 weeks post-target quit day (TQD). Participants
completed one morning EMA prompt, one random prompt during the day, and one
evening prompt every day for the first 2 weeks post-TQD and then every other
day for the remaining 2 weeks of the EMA period. Thus, there are up to 66 EMA
measurements that are not equally spaced due to fluctuations in the timing of
the prompts and to sporadic missingness. Our proposed approach allows for
these irregularly spaced assessments. The goal of our empirical analysis is to
examine whether the intervention of interest (varenicline) has an effect on
smoking cessation (7-day point-prevalence abstinence at 4 weeks post-TQD) that
is mediated by a reduction in craving, compared to any nicotine replacement
therapy. Craving is an important predictor of relapse in many substance use
disorders (Vafaie Kober, 2022; Le Foll ., 2022), and it was assessed at each
EMA prompt during the 4 week post-TQD period. It is treated as a time-varying
mediator in Section 6.

##  3 Proposed Approach

###  3.1 Notation

Consider a sample of $N$ individuals, in which each individual
$i,i=1,\ldots,N$ is randomly assigned to one of two treatments, $X_{i}=0$, the
control or comparison group, or $X_{i}=1$, the active treatment group. Note
that this notation is easily extended to more than two treatment conditions.
Each individual $i$ is assessed at $T_{i}$ measurement occasions, which may
vary across individuals, between times $t_{\mathrm{min}}$ and
$t_{\mathrm{max}}$. At each measurement occasion, $t_{ij},j=1,\ldots,T_{i}$,
such that $t_{\mathrm{min}}\leq t_{1}\leq\ldots\leq t_{T_{i}}$, the mediator,
$M_{i}(t_{ij})$ is assessed. The end-of-study outcome, denoted $Y_{i}$, is
measured at a distal time point, which is after $t_{\mathrm{max}}$ and common
to all individuals. Let
$\boldsymbol{M}_{i}=\\{\boldsymbol{M}_{i}(t):t\in[t_{\mathrm{min}},t_{\mathrm{max}}]\\}$
denote the complete mediator process for individual $i$ from beginning to end
of the study interval, including times at which it is not observed. In the
context of the potential outcomes framework, we additionally use parentheses
to denote the counterfactual outcomes under each treatment condition.
Specifically, let $\boldsymbol{M}_{i}(x)$ be the potential value that the
complete mediator process would take for individual $i$, if that individual,
perhaps counterfactually, received treatment level $x$. Let
$Y_{i}(x,\boldsymbol{M}_{i}({x^{\ast}}))$ be the potential outcome if
individual $i$, perhaps counterfactually, is given treatment $X=x$ and
mediator process $\boldsymbol{M}_{i}({x^{\ast}})$.

We make the stable unit treatment value assumption (SUTVA; Rubin, 1980) of
non-interference between individuals. That is, an individual’s potential
values on the mediator process are not dependent on any other individual’s
treatment assignment and an individual’s potential outcomes are not dependent
on any other individual’s treatment assignment or potential values of their
mediator process. SUTVA also states that there are not multiple versions of
the treatment (see Hernán Robins, 2020). In addition, we make a consistency
assumption, in which the counterfactual outcomes $Y_{i}(x)$ and
$\boldsymbol{M}_{i}(x)$ are equal to the observed outcomes, $Y_{i}$ and
$\boldsymbol{M}_{i}$, respectively, when individual $i$ is assigned to
treatment level $x$, and the composition assumption, i.e.,
$Y_{i}(x)=Y_{i}(x,\boldsymbol{M}_{i}(x))$ (VanderWeele Vansteelandt, 2009,
2010; Lindquist, 2012). Finally, let $Y_{i}(x,\boldsymbol{m})$ denote the
potential outcome for individual $i$ given treatment level $x$ and given
mediator process, $\boldsymbol{m}$. That is,
$\boldsymbol{M}_{i}(x)=\boldsymbol{m}$, a specific mediator process. Where
possible, to avoid confusion, we drop the $i$ subscript for simplicity.

###  3.2 Definitions

In this paper we focus on the natural direct and indirect effects, similar to
VanderWeele Vansteelandt (2010) and Rijnhart . (2021). In contrast,
VanderWeele Tchetgen Tchetgen (2017) and Pinto Pereira . (2020), focus on the
interventional direct and indirect effects; for an explication of differences
between these definitions, see Moreno-Betancur Carlin (2018) and Nguyen .
(2021).

If an identity link is used, the causal effect is most readily expressed in
terms of a mean difference. In this case, the average causal total effect of
the treatment on the distal end-of-study outcome is defined as:

$\mathrm{TE}\equiv
E(Y(1,\boldsymbol{M}(1)))-E(Y(0,\boldsymbol{M}(0)))=E(Y(1))-E(Y(0))$

which can be decomposed into the natural direct effect:

$\mathrm{NDE}(0)=E(Y(1,\boldsymbol{M}(0)))-E(Y(0,\boldsymbol{M}(0)))$

and the natural indirect effect:

$\mathrm{NIE}(1)=E(Y(1,\boldsymbol{M}(1)))-E(Y(1,\boldsymbol{M}(0)))$

obtained by adding and subtracting $E(Y(1,\boldsymbol{M}(0)))$. The natural
direct effect, $\mathrm{NDE}(0)$, defines the effect of the treatment on the
outcome holding the mediator process to the potential value obtained under the
control condition (i.e., $x=0$). The natural indirect effect,
$\mathrm{NIE}(1)$, defines the effect on the outcome of the potential mediator
process under treatment vs. the potential mediator process under control while
holding treatment level fixed to treated (i.e., $x=1$) (see, e.g. VanderWeele,
2015). Alternatively, one could add and subtract $E(Y(0,\boldsymbol{M}(1)))$
to obtain an equivalent decomposition:
$\mathrm{NDE}(1)=E(Y(1,\boldsymbol{M}(1)))-E(Y(0,\boldsymbol{M}(1)))$ and
$\mathrm{NIE}(0)=E(Y(0,\boldsymbol{M}(1)))-E(Y(0,\boldsymbol{M}(0)))$. We use
the former decomposition, $\mathrm{TE}=\mathrm{NDE}(0)+\mathrm{NIE}(1)$,
although the procedures are applicable to the latter as well.

For modeling binary or count data, a log or logistic link may be used instead
of an identity link function. That is,

$\displaystyle g(\mathrm{TE})=$ $\displaystyle
g\big{(}E(Y(1,\boldsymbol{M}(1)))\big{)}-g\big{(}E(Y(0,\boldsymbol{M}(0)))\big{)}$
$\displaystyle g(\mathrm{NDE}(0))=$ $\displaystyle
g\big{(}E(Y(1,\boldsymbol{M}(0)))\big{)}-g\big{(}E(Y(0,\boldsymbol{M}(0)))\big{)}$
$\displaystyle g(\mathrm{NIE}(1))=$ $\displaystyle
g\big{(}E(Y(1,\boldsymbol{M}(1)))\big{)}-g\big{(}E(Y(1,\boldsymbol{M}(0)))\big{)}$

where $g$ represents the log or logit link. In this case, the causal effect
definitions are on the ratio scale (either rate ratio for the log link or odds
ratio for the logit link). For the log link this would be:

$\displaystyle\mathrm{TE}=$ $\displaystyle
E\big{(}Y(1,\boldsymbol{M}(1))\big{)}/E\big{(}Y(0,\boldsymbol{M}(0))\big{)}$
(1) $\displaystyle\mathrm{NDE}(0)=$ $\displaystyle
E\big{(}Y(1,\boldsymbol{M}(0))\big{)}/E\big{(}Y(0,\boldsymbol{M}(0))\big{)}$
$\displaystyle\mathrm{NIE}(1)=$ $\displaystyle
E\big{(}Y(1,\boldsymbol{M}(1))\big{)}/E\big{(}Y(1,\boldsymbol{M}(0))\big{)}.$

Equivalently, the log of each effect is a difference of log probabilities or
log odds, respectively.

###  3.3 Identification Assumptions

In addition to the aforementioned SUTVA, consistency, and composition
assumptions, identifying the potential outcomes from the observed data
requires the following ignorability assumptions given baseline covariates,
$\boldsymbol{C}$, as shown in Figure 2. In addition, we require the positivity
assumption which states that each individual has a positive probability of
being assigned to either treatment and to a given mediator process given their
covariate values, $P(X=x|\boldsymbol{C}=\boldsymbol{c})>0$ and
$P(\boldsymbol{M}_{i}(x)=\boldsymbol{m}|X=x,\boldsymbol{C}=\boldsymbol{c})>0$.
The assumptions allow the identification of the potential outcomes from the
observed data. The assumptions are adapted from VanderWeele Vansteelandt
(2010) (see also Imai ., 2010), but the first three are analogous to
Assumption 1 of Lindquist (2012), and the last is analogous to Assumption 2 of
Lindquist (2012).

  1. 1.

Conditional on $\boldsymbol{C}$, the effect of the treatment on the outcome is
unconfounded (i.e., there are no unmeasured confounders of the treatment and
outcome): $Y_{i}(x,\boldsymbol{m})\perp X_{i}|\boldsymbol{C}_{i}$.

  2. 2.

Conditional on $X$ and $\boldsymbol{C}$, the effect of the mediator on the
outcome is unconfounded (i.e., there are no unmeasured confounders of the
mediator and outcome):
$Y_{i}(x,\boldsymbol{m})\perp\boldsymbol{M}_{i}|X_{i},\boldsymbol{C}_{i}$.

  3. 3.

Conditional on $\boldsymbol{C}$, the effect of the treatment on the mediator
is unconfounded (i.e., there are no unmeasured confounders of the treatment
and mediator): $\boldsymbol{M}_{i}(x)\perp X_{i}|\boldsymbol{C}_{i}$.

  4. 4.

There are no measured or unmeasured confounders of the mediator and outcome
that have been influenced by the treatment:
$Y_{i}(x,\boldsymbol{m})\perp\boldsymbol{M}_{i}({x^{\ast}})|\boldsymbol{C}_{i}$
for all $x$ and ${x^{\ast}}$.

Assumptions 1-3 are standard sequential ignorability assumptions (Imai .,
2010). Given randomization to the treatment levels, the first and third
assumptions are reasonable. In contrast, Assumption 2 may not hold even when
there is randomization to the treatment levels, because there is not
randomization to the mediating process; however, if a rich set of potential
confounders have been measured, then this assumption may still be reasonable
conditional on the measured confounders.

Assumption 4 (often referred to as cross-world independence) may be
problematic. Although the outcome occurs at a distal time point in our
context, and therefore is not itself a time-varying confounder (e.g., as in
Zeng ., 2021), there may be other potential mediating processes which may
violate this assumption. For the present motivating example, other variables
were measured during the EMA period, such as smoking cessation self-efficacy,
self-reported daily smoking, and smoking cessation fatigue. It may be expected
that the treatment may influence these variables, which in turn may affect
both the mediating process of interest (i.e., craving) and the outcome, (i.e.,
abstinence at one month post-TQD). We further discuss possibilities for
allowing these additional mediators in Appendix 2.

There are possible alternative assumptions to Assumption 4 when it may be
violated (see, e.g. VanderWeele, 2015; Robins, 2003). Rather than assuming
that there are no mediator-outcome confounders affected by treatment (i.e.,
Assumption 4), an alternative identifying assumption is to assume that there
are no interactions between treatment assignment and the mediating process.
However, this assumption is not as innocuous as it may seem (Robins, 2003). It
means that the natural direct and indirect effects are the same in both the
treatment and control groups (i.e., $\mathrm{NDE}(0)=\mathrm{NDE}(1)$ and
$\mathrm{NIE}(0)=\mathrm{NIE}(1)$). In this case the controlled direct effect,
defined as $Y(1,\boldsymbol{m})-Y(0,\boldsymbol{m})$ for some specified fixed
trajectory $\boldsymbol{m}$ for the mediator process, (i.e.,
$\boldsymbol{M}(x)=\boldsymbol{m}$), is equal for any level of
$\boldsymbol{m}$ and is also equal to the $\mathrm{NDE}$. Then the
$\mathrm{NIE}$ can be parametrically identified as $\mathrm{TE}-\mathrm{NDE}$,
as in the traditional linear structural equation modeling approach to
mediation. However, this lack of interaction may also be difficult to justify.
Note that Zeng . (2021), Lindquist (2012), and Zhao . (2018) all impose
Assumption 4 or the no-interaction assumption.

Lindquist (2012) discusses an alternative approach to estimating the causal
effect of the mediator on the outcome, without making Assumption 4. He
proposes the intuition of using the ratio of the total effect of $X$ on $Y$ to
the effect of $X$ on $M$ as a measure of the strength of mediation through
$M$. This could be described as using $X$ as an instrumental variable for the
purposes of the causal effect of $M$ on $Y$. However, this alternative
approach requires an exclusion restriction assumption (no direct effect of $X$
on $Y$, i.e., complete mediation) which will often be unrealistic. In the
current example, it does not seem plausible to assume that reduced craving is
the only mechanism by which treatment affects relapse probability.

In any case, identifying the natural direct and indirect effects from the
observed data requires strong assumptions. If Assumption 4, the no-interaction
assumption, or the exclusion restriction are not realistic, then it may be
possible to define interventional direct and indirect effects (see VanderWeele
Tchetgen Tchetgen, 2017), rather than the natural direct and indirect effects
considered here. This approach is rather new, and applying it to functional
regression requires further methodological research, but it is a promising
future direction.

Given Assumptions 1 - 3 and the no-interaction assumption (or Assumption 4),
the natural direct and indirect effects can be identified from the observed
data. The causal effects can be estimated as functions of the expected values
(e.g., $E[Y|X,\boldsymbol{C},\boldsymbol{M})$) of the observed data. Next, we
describe the models we use to estimate these quantities.

##  4 Estimation

In general, the estimation strategy is a two-step process in which first a
model is fit to predict the mediator from treatment and baseline covariates,
and then another model is fit to predict the outcome from the mediator,
treatment, and baseline covariates.

###  4.1 Model for Mediator

For the identity link, the model for the potential values of the mediator at
time $t$ for individual $i$ if given treatment $x$ is:

$M_{i}(t;x)=\alpha_{0}(t)+\boldsymbol{C}_{i}\boldsymbol{\alpha}_{C}(t)+\alpha_{X}(t)x+\epsilon_{i}(t)$
(2)

where $\epsilon_{i}$ is a random process representing measurement-level errors
for the mediator, which are assumed to have $E(\epsilon_{it})=0$ for all $t$
and may be correlated within person, and where $\alpha_{0}(t)$,
$\boldsymbol{\alpha}_{C}(t)$, and $\alpha_{X}(t)$ are smooth functions of
time. By using the consistency assumption, we can simply write the potential
value, $M_{i}(t;x)$ which is a function of both $t$ and treatment level $x$,
as $M(t)$ when expressing an expected value conditional on $x$.

More generally, we can model $M(t)$ as

$g(E(M(t)|X=x,\boldsymbol{C}))=\alpha_{0}(t)+\boldsymbol{C}\boldsymbol{\alpha}_{C}(t)+\alpha_{X}(t)X$
(3)

where $g$ is a link function such as log, logit, or identity.

Model (2) or (3) is an application of concurrent function-on-scalar regression
(Ramsay Silverman, 2005) or equivalently the time-varying effect model (TVEM;
Hastie Tibshirani, 1993; Tan ., 2012). The link function $g$ is treated as
given and known, but the $\alpha$ functions are not assumed to have any simple
parametric form known to the data analyst.

###  4.2 Model for Outcome

The model for the potential outcome using an identity link is given as:

$Y_{i}(x_{i},\boldsymbol{M}_{i}(x_{i}))=\beta_{0}+\boldsymbol{C}_{i}\boldsymbol{\beta}_{C}+\beta_{X}x_{i}+\int_{t_{\mathrm{min}}}^{t_{\mathrm{max}}}\beta_{M}(t)\boldsymbol{M}_{i}(t;x)\mathrm{dt}+\eta_{i},$
(4)

where the $\eta$ are subject-level errors for the response assumed to have
$E(\eta|X=x,\boldsymbol{M}(X)=\boldsymbol{m})=0$. The subject-level errors
$\eta_{i}$ for the outcome in (4) are treated as independent of the
measurement-level errors $\epsilon_{i}(t)$ for the mediator in (2).

In general, we can model the expected value of the observed outcome, $Y$,
given treatment $x$, covariate vector $\boldsymbol{c}$, and mediator
trajectory $\boldsymbol{m}$, as

$g(E(Y|X=x,\boldsymbol{C}=\boldsymbol{c},\boldsymbol{M}=\boldsymbol{m}))=\beta_{0}+\beta_{X}x+\boldsymbol{c}\boldsymbol{\beta}_{C}+\int_{t_{\mathrm{min}}}^{t_{\mathrm{max}}}\beta_{M}(t)\boldsymbol{m}(t)\mathrm{dt}$
(5)

where $g$ is a link function such as an identity, log, or logit. A model such
as (4) and (5) represents scalar-on-function regression (Ramsay Silverman,
2005; Goldsmith ., 2011; Lindquist McKeague, 2009).

###  4.3 Expression for the Natural Indirect Effect

Using the identity link, and under Assumptions 1-3 and the no treatment-
mediator interaction assumption, the $\mathrm{NIE}(0)$ and the
$\mathrm{NIE}(1)$ are equal. Thus, the natural indirect effect is:

$\displaystyle\mathrm{NIE}$
$\displaystyle=E(Y(0,\boldsymbol{M}(1))-Y(0,\boldsymbol{M}(0)))$
$\displaystyle=E(Y(1,\boldsymbol{M}(1))-Y(1,\boldsymbol{M}(0)))$
$\displaystyle=\int_{t_{\mathrm{min}}}^{t_{\mathrm{max}}}\alpha_{X}(t)\beta_{M}(t)\mathrm{dt}$

Lindquist (2012) proved this for the identity link.

For a binary outcome using a log link,

$\displaystyle\mathrm{NIE}$
$\displaystyle=E(Y(0,\boldsymbol{M}(1)))/E(Y(0,\boldsymbol{M}(0)))$
$\displaystyle=E(Y(1,\boldsymbol{M}(1)))/E(Y(1,\boldsymbol{M}(0)))$
$\displaystyle=\exp\Big{(}\log\big{(}E(Y(0,\boldsymbol{M}(1)))\big{)}-\log\big{(}E(Y(0,\boldsymbol{M}(0)))\Big{)}$
$\displaystyle=\exp\Big{(}\int_{t_{\mathrm{min}}}^{t_{\mathrm{max}}}\alpha_{X}(t)\beta_{M}(t)\mathrm{dt}\Big{)},$

and the proof is given in Appendix 1. For binary data using a logistic link,
it is more difficult to express the decomposition of the total effect into
direct and indirect effects in a straightforward way, but the logit link
function approximates a log link function when outcomes are rare (see
VanderWeele Vansteelandt, 2010; Rijnhart ., 2021). In any case, it remains
reasonable to use a significance test for
$\int\alpha_{X}(t)\beta_{M}(t)\mathrm{dt}$ as a test for significant
mediation, by analogy to the product of coefficients method in traditional
mediation analysis. Another option, not pursued here, would be to modify the
definition of direct and indirect effects in order to make the logistic link
easier to handle theoretically, following either the ideas in Geldhof . (2017)
or Doretti . (2021).

The controlled direct effect of $X$ on $Y$ given a fixed mediator trajectory
$\boldsymbol{m}$ is $E(Y(1,\boldsymbol{m}))-E(Y(0,\boldsymbol{m}))=\beta_{X}$
for the identity link or
$E(Y(1,\boldsymbol{m}))/E(Y(0,\boldsymbol{m}))=\exp(\beta_{X})$ for the log
link. Because of the assumption of no treatment-mediator interaction, the
controlled direct effect of $X$ given $\boldsymbol{M}(x)=\boldsymbol{m}$
equals both the $\mathrm{NDE}(0)$ and $\mathrm{NDE}(1)$.

###  4.4 Model for Total Effect

Finally, the total effect using an identity link is given as follows:

$\displaystyle\tau$ $\displaystyle\equiv
E(Y(1,\boldsymbol{M}(1)))-E(Y(0,\boldsymbol{M}(0)))$
$\displaystyle=\beta_{X}+\int_{t_{\mathrm{min}}}^{t_{\mathrm{max}}}\alpha_{X}(t)\beta_{M}(t)\mathrm{dt}$

For the log link, the total effect is additive on a log scale,

$\exp(\tau)=\exp(\beta_{X})+\exp\Big{(}\int_{t_{\mathrm{min}}}^{t_{\mathrm{max}}}\alpha_{X}(t)\beta_{M}(t)\mathrm{dt}\Big{)},$

but is multiplicative on the ratio scale. Again, a logit link is also
possible, although it can make the relationship between direct, indirect, and
total effects difficult to express unless the outcome is rare.

###  4.5 Test for Indirect Effect

To test for an indirect effect, one can test whether a nonparametric
bootstrapped confidence interval for
$\int\alpha_{X}(t)\beta_{M}(t)\mathrm{dt}$ includes zero (Lindquist, 2012). In
this context, the indirect effect is a scalar, even though its individual
effects, $\alpha_{X}(t)$ and $\beta_{M}(t)$, are functions of $t$.

###  4.6 Time-varying effect models

Models (2) and (3) can be fit as longitudinal varying-coefficients models,
i.e., time-varying effect models (TVEM; see Tan ., 2012). In general, a TVEM
model with a time-varying response variable, $M(t)$, is specified as

$g(E(M(t)))=\alpha_{0}(t)+\sum_{j}\alpha_{j}(t)X_{j}(t)$

for a set of longitudinal predictors $X_{j}(t)$ observed on a sparse and
possibly irregular grid. In this general version, both the covariates and
their “effects” (time-specific regression coefficients) are time-varying.
However, any of the $\alpha_{j}$ and/or $X_{j}$ can be optionally constrained
to be constant over time.

Notice that for parsimony and interpretability, and to avoid the curse of
dimensionality, TVEM assumes a linear or generalized linear structure at any
given time, and that the estimates of $\alpha_{0}$ and the $\alpha_{j}$’s are
smooth but may potentially be nonlinear functions of time. Thus, Models (2)
and (3) are included as a special case. The TVEM model is straightforward to
fit as a marginal model using working independence and robust standard errors,
and a spline approximation to the functional coefficients. Thus, a TVEM
marginal model is essentially a generalized estimating equations (Liang Zeger,
1986) model with some of the parametric terms replaced by splines.

We fit the model by representing the functional coefficients as penalized
B-splines (Eilers Marx, 1996), as implemented in the `tvem` function of the
`tvem` R package (Dziak ., 2021). The model is fitted approximately by first
constructing a spline basis for the time interval of interest, and then
performing a large parametric regression including each function of time in
the basis, as well as its product with each coefficient. We use a penalty
function on the nonlinear functions of time in order to avoid overfitting and
control the level of smoothness versus flexibility of the fitted function. By
default, `tvem` uses cubic B-splines with a maximum of 20 interior knots and a
first-order difference penalty on the knot coefficients (see Eilers Marx,
1996). Penalized truncated power splines (see, e.g., Tan ., 2012; Huang Yuan,
2017) are an alternative option.

TVEMs can be viewed as a special case of the varying-coefficient model (Hastie
Tibshirani, 1993) in which the coefficients are indexed by time. The varying
coefficients allow the effect of treatment on the mediator to change as a
nonparametric function of time. Additional time-varying and time-invariant
predictors can also be added to the model. Predictors can be specified to have
either time-varying or time-invariant effects. For example, demographic
variables or treatment group membership may have time-varying relationships
with the outcome variable even if they are themselves not time-varying. In the
tvem package, within-subject correlation is accounted for indirectly by using
sandwich standard errors, similar to working-independence GEE (Liang Zeger,
1986). This can be seen as an application of penalized pseudo-likelihood,
incorporating the penalty required to regularize the spline functions. The
penalty function is tuned by a restricted maximum likelihood (REML) approach
(Wood, 2001). This approach shrinks the changes at each knot, acting
analogously to an empirical Bayes prior. The prior distribution used in the
tvem package is a first- or second-order difference penalty as described in
Eilers & Marx (1996), and it acts to shrink the fitted curve towards a
straight line (possibly even a flat line with zero slope) unless data suggests
otherwise.

###  4.7 Scalar-on-function regression

Although the underlying process, $\boldsymbol{M}_{i}$, is assumed to be a
continuous function, practically it is observed at only a finite number of
points and potentially has measurement error. For computational purposes, it
would be very helpful to have values for $M_{i}(t)$ on a dense, regular grid
without missing values. However, this is not likely to be the case in
practice. One option for working around this difficulty, if the amount of
missing data is not too great, is by linear interpolation between nonmissing
points. Another option is to use spline regression to fit a smoothed curve
$\widehat{\boldsymbol{M}}_{i}$ to the observed data for each participant $i$,
and then using these smoothed values to estimate the $\beta_{M}(t)$ function.
This smoothing could be heuristically justified by supposing that the mediator
is measured with error, and that the true trajectory of the mediator was
smoother than the trajectory which would have been recorded if all possible
measurement occasions on the grid had been observed for all participants. That
is, we let $\widehat{M}_{i}(t)$ at time $t_{ij}$ be a possibly smoothed
prediction of the true latent curve underlying $M_{i}(t)$, rather than the
observed (and presumably error-contaminated) $M_{i}(t)$ itself (Goldsmith .,
2011) by incorporating a presmoothing step in the model (see Appendix 1).

After this presmoothing, the fitted model is

$g(E(Y_{i}|\boldsymbol{M}_{i}(t))\approx\beta_{0}+\boldsymbol{\beta}_{C}\boldsymbol{C}_{i}+\beta_{X}X_{i}+\intop_{t_{\mathrm{min}}}^{t_{\mathrm{max}}}\beta_{M}(t)\widehat{M}_{i}(t)dt.$
(6)

Additional scalar (i.e., time-invariant) predictor variables with scalar
coefficients may be added to Model (6). The coefficients, both scalar and
functional, are estimated following the general method of Goldsmith . (2011)
using REML. The coefficient functions are represented via the penalized thin
plate regression splines approach described by Wood (2003) and Wood (2017).
Notice that the presence of within-subject variability is required, because if
$\boldsymbol{M}_{i}(t)$ were constant over $t$ for each $i$, then
$\beta_{M}(t)$ would not be identified. Accordingly, with substantial and
meaningful within-subject variability, the functional regression coefficient,
$\beta_{M}(t)$, requires regularization using a penalty function and/or a low-
dimensional basis restriction.

###  4.8 Overall indirect effect computation

Once $\hat{\alpha}_{X}(t)$ and $\hat{\beta}_{M}(t)$ have been estimated, the
indirect effect is operationally defined as the integral
$\int_{t_{\mathrm{min}}}^{t_{\mathrm{max}}}\hat{\alpha}_{X}(t)\hat{\beta}_{M}(t)dt$.
This is approximated by taking the mean of
$\hat{\alpha}_{X}(t)\hat{\beta}_{M}(t)$ for an evenly spaced grid of points
$t$ and then rescaling by multiplying it by the width of the time interval of
interest.

Next, a nonparametric bootstrap significance test for an overall $p$-value for
$\int_{t_{\mathrm{min}}}^{t_{\mathrm{max}}}\alpha_{X}(t)\beta_{M}(t)\mathrm{dt}$
is obtained. A simple approximate confidence interval for this overall
quantity using either a normal approximation or the percentile method can be
obtained using the funmediation package.

###  4.9 Difficulties with time-specific indirect effects

We focus on an overall test of mediation, rather than time-specific mediation
estimates, and do not attempt to interpret $\alpha_{X}(t)\beta_{M}(t)$ at a
specific $t$. This is in contrast to Lindquist (2012) who used a “wild
bootstrap” to construct confidence intervals for the point-specific indirect
effect $\alpha_{X}(t)\beta_{M}(t)$ at each of a grid of points $t$. However,
the wild bootstrap is more difficult to apply to binary outcomes, since it
involves permutation of residuals. Also, the practical meaning of the time-
specific indirect effect $\alpha_{X}(t)\beta_{M}(t)$ with EMA data could be
more difficult to interpret because of temporal dependence (see Dziak ., 2019;
VanderWeele Tchetgen Tchetgen, 2017) and because time-specific snapshots may
not adequately represent time-lagged effects (see Maxwell ., 2011); we address
this further in the Discussion.

##  5 Simulations

We performed a simulation study to determine whether the proposed method has
adequate coverage and power under realistic scenarios. We simulated 5000
datasets for each of three sample size scenarios ($N=250$, $N=500$, and
$N=1000$). In each dataset, the treatment $X$ is dichotomous and generated
with equal probability of 0 (representing control group) and 1 (representing
treated group). We imagine an interval from 0 to 1 representing rescaled time
in a smoking cessation program, $X$ as a randomized intervention to reduce
cravings, $M_{i}(t)$ as recent amount of cravings reported at time $t$, and
$Y$ as relapse status by end of program, with $Y=1$ indicating relapse. These
simulated variables map onto the applied example, but the specific values of
the simulated variables represent a much wider array of possible variable
combinations.

###  5.1 Methods

The mediator is generated on a grid of 100 time points on an interval from 0
to 1 (e.g., 0.00, 0.01, …, 0.99). However, 60% of those observations are
randomly deleted, to create uneven person-specific measurement times.
$\boldsymbol{M}_{i}$ as a whole is assumed to follow a Gaussian process with
mean function $E(M_{i}(t)|X_{i})=\alpha_{0}(t)+\alpha_{X}(t)X_{i}$. We set the
time-varying mean of the mediator for the control group to
$\alpha_{0}(t)=\sqrt{t}$, and the time-varying effect of $X$ on the mediator
to $\alpha_{X}(t)=-\sqrt{t/2}$. This means that both the control and treatment
groups start on average with a mean near 0 on the mediator, but the control
group increases nonlinearly over time to near $1$, while the treatment group
increases nonlinearly to near $0.3$. The error process for the mediator is
generated to have a standard deviation of $2$ and an AR(1) structure with
correlation $0.8$ between adjacent measurements on the full 100-point grid.

The outcome $Y_{i}$ was generated as binary, with a logistic link.
Specifically, success probability was given by

$\mathrm{logit}(E(Y_{i}))=\beta_{0}+\beta_{X}X_{i}+\int_{0}^{1}\beta_{M}(t)M_{i}(t)\mathrm{dt}$

where $\beta_{0}=0$, $\beta_{X}=0.2$, and $\beta_{M}=\frac{1}{2}(\exp(t)-1)$.
Therefore, $X$ tends to reduce $M(t)$, and $M(t)$ is positively related to
$Y$, such that the overall indirect effect is negative (a negative coefficient
times a positive coefficient). Specifically, the overall true value of the
indirect effect is taken to be
$\int_{0}^{1}\alpha_{X}(t)\beta_{M}(t)=-0.2113$. Because the treatment is
simulated to be beneficial and the simulated outcome is imagined as
deleterious, the negative effect is reasonable.

The model was fit to each dataset using the proposed combination of TVEM and
scalar-on-function regression. We tested the null hypothesis that the indirect
effect was zero using 199 bootstrap replications. In a real data analysis, 199
would be inadequate, and at least 999 would be recommended. However, we found
that the larger number was computationally infeasible in a simulation setting,
because 1000 simulated datasets per scenario with 999 bootstrap replications
would mean having to fit and tune the model separately for each of
approximately one million datasets per scenario.

The R software environment (R Core Team, 2020) was used, with author-written
functions available in the package `funmediation` available on CRAN. This
package employs the `tvem` package (Dziak ., 2021) and the `pfr` function in
the `refund` package (Goldsmith ., 2020). Each package uses generalized
additive model functionality from the `gam` or `bam` functions in the `mgcv`
package (Wood, 2017, 2021) for back-end calculations. Presmoothing is handled
via interpolation using the `presmooth` option in the `lf` function in the
`refund` R package. Bootstrapping is performed using the `boot` function in
the `boot` package (Canty Ripley, 2020; Davison Hinkley, 1997).

###  5.2 Results

Performance measures for estimates of the relationship of $X$ to $M(t)$ are
shown in Table 1. Estimates for $\alpha_{0}(t)$ and $\alpha_{X}(t)$ had very
little bias and had increasing precision as sample size increased.

Performance measures for estimates of the relationship of $X$ and $M(t)$ to
$Y$ are shown in Table 2. Results indicate a very small bias of the estimated
effects towards zero, which does not seem to vanish with increasing sample
size. Nonetheless, random error decreases with increasing sample size, and
pointwise coverage for the $95\%$ pointwise confidence intervals is
essentially nominal for all sample sizes, so performance can generally be
considered good. The small but non-vanishing bias present may be a side effect
of the shrinkage penalty functions used to prevent overfitting of $\alpha_{X}$
and $\beta_{M}$.

The confidence intervals for $\beta_{M}(t)$ are constructed only for pointwise
coverage (i.e., a 95% long-term probability of containing the true value at
any fixed time point), not nominal familywise coverage (i.e., a 95% long-term
probability of covering the true value for all the time points jointly), so
they would not perform as well on familywise coverage. Lindquist (2012)
recommended false discovery rate correction over values of $t$. However, in
the case of observational studies with potentially diffuse and weak effects
over time, this might lead to a very high Type II error rate. That is,
researchers might misinterpret points which are indicated as significant at a
given $\alpha$ level as being fundamentally different from time points
indicated as nonsignificant, i.e., as times when the exposure or treatment are
active and inactive, respectively, despite both being highly affected by
sample size and sampling error. It is not yet known how best to handle this
issue. Fortunately, it is not a problem for the overall indirect effect
$\int_{0}^{1}\alpha_{X}\beta_{M}\mathrm{dt}$, because it is a scalar and only
tested once.

Performance measures for estimates of the indirect relationship of $X$ to $Y$
are shown in Table 3. There was a small conservative bias towards zero, which
became even smaller for larger sample sizes. Mean squared error also decreased
for larger sample sizes. Power and coverage are shown for the normal-theory
bootstrap, basic bootstrap, and percentile bootstrap. Of the three methods,
the best coverage and power are obtained using the percentile bootstrap.
Coverage is very good, and power increases with larger sample sizes. However,
adequate power requires more than 500 participants.

##  6 Empirical Analysis

In this section, we present the results from analyzing the data described in
Section 2. The R code for implementing this analysis is available in the
supplementary materials.

###  6.1 Sample Characteristics and Missing Data

Our analysis sample, taken from the 1086 participants in the WSHS2 (Baker .,
2016), consisted of 1041 participants who were randomized, provided data to at
least one EMA prompt post-TQD, and whose abstinence or non-abstinence status,
$Y$, could be ascertained. There was very little missingness for $Y$, which
was not subject to noticeable differential missingness across treatment group
(11 missing $Y$ values from among those randomized to patch group, 13 from the
varenicline group, 11 from the cNRT group). Thus, missingness is not a large
concern for the outcome, $Y$.

However, participants varied widely in the proportion of response to the total
66 EMA prompts possible, so the mediator $M(t)$ at any given time $t$ was
subject to noticeable missingness. This may be a concern for assessing the
relationship between $X$ and $M(t)$ or between $M(t)$ and $Y$.

Among the 1086 participants, all of whom have $X$ (randomization status), 45
did not provide any EMA data, and they were excluded from the analysis. The
1041 participants who were included provided between 1 and 66 observations per
participant (mean=41, median=45). There were three EMA prompts per day:
morning, evening, and a random time in between. For participants who responded
to most prompts, the middle, random, prompt was more likely to be missing than
either the morning or evening prompts. The random prompt may have been missing
because the participant merely happened to be temporarily unavailable to
answer. This is likely to be less of a problem than long-term dropout.

Because the current analysis is intended mainly as an illustration of a new
technique, we do not attempt to solve the possible problem of informative
missingness here. However, the serious potential limitations which it poses,
and possible future approaches to address them, are discussed in Section 7. We
conducted an additional analysis in which participants who responded to two or
fewer EMA prompts post-TQD were excluded as a sensitivity analysis to
partially assess the impact of missingness of functional data on estimation
using this technique. The results are similar to those presented here and are
included in Appendix 3.

###  6.2 Total Effect

For comparison, we fit a logistic regression model of the outcome, CO-
confirmed 7-day point-prevalence abstinence (1=abstinent, 0=relapsed), at
four-weeks post-TQD on treatment group (1 = varenicline, 0 = nicotine
replacement therapy). The total effect estimate from this model was not
statistically significant (estimate = 0.023, std. error = 0.135, $z=0.167$,
$p=0.868$). Because this study does not have a no-treatment control group, and
instead focuses on comparative efficacy between active treatments, this null
finding may be due to the true difference between treatment efficacy being
small (i.e., small effect size for comparing treatment levels). However, the
example is still a useful illustration of the proposed procedures.

###  6.3 Direct and Indirect Effects

Next, we examine the direct and indirect effects of craving, defined as the
desire to smoke in the last 15 minutes, throughout the first two weeks post-
TQD using the proposed method. The results of these analyses are unique,
providing insights that are not discernible from any other analytic approach,
including other mediation approaches. The indirect effect estimate was
statistically significant (estimate = 0.044, percentile method 95% bootstrap
confidence interval = (0.006, 0.111)). The estimate of the direct effect was
not statistically significant (estimate = -0.067, std. error = 0.139). The
indirect and direct effects do not add up exactly to the total; this is likely
due to the effects of smoothing of the mediator and of the estimated effects
(see Appendix 1). It is likely that oversmoothing or undersmoothing of the
mediator trajectories may have an effect like that of measurement error, which
is known to cause bias (see Fritz ., 2016; de Stavola ., 2022).

###  6.4 Effect of Treatment on Mediator

Examination of the individual effects that comprise the indirect effect
revealed a more nuanced interpretation. Figure 3 illustrates the time-varying
effect of treatment on the mediator. The effect is such that varenicline
reduces craving and this effect becomes slightly stronger (i.e., more
negative) during the first week post-TQD before leveling off. Other than the
first day post-TQD, the effect is statistically significant throughout the
two-week post-TQD period, at least in the heuristic sense that the pointwise
confidence intervals generally do not include zero.

###  6.5 Effect of Mediator on Outcome

The effect of the mediator on the outcome, CO-confirmed 7-day point prevalence
abstinence at 4 weeks post-TQD, is shown in Figure 4. Results from functional
regression are more difficult to interpret precisely than those from a time-
varying effect model (Dziak ., 2019). Nevertheless, it is evident from the
confidence intervals in Figure 4 that for the most part, the confidence
intervals include 0, with the exception of the end of the first week post-TQD.
All else being equal, craving in the first week (especially around day 4) has
a negative weight in predicting abstinence, while craving in the second week
(especially around day 10) has a weak and perhaps non-significant positive
weight. A negative weight followed by a positive weight may be an indication
that abstinence is more likely in individuals for whom craving increases over
time. However, this explanation is intuitively unlikely because abstainers are
expected to experience a decrease in craving over time, so it might be more
likely that the result simply means that very high early craving can result in
relapse. The wide confidence intervals also make interpretation more
difficult.

##  7 Discussion

In this paper, we proposed and illustrated a method for estimating and testing
the indirect effect of a randomized treatment on a distal binary outcome
variable, as mediated by the nonparametric trajectory of an intensively
measured longitudinal process. This continued the work of Lindquist (2012), in
combining causal analysis in the counterfactual framework, with regularized
overparameterized semiparametric estimation. These two sets of ideas are among
the most important in modern statistics (Gelman Vehtari, 2021). In extending
Lindquist’s ideas to a binary setting, we used simulations and a real-data
example to illustrate the potential use of functional mediation analysis with
binary outcome data. We provide a new R package, `funmediation`, for fitting
this and similar models.

###  7.1 Limitations

A remaining methodological challenge in functional regression modeling
approaches is the issue of missing data and possible informative dropout. If
some observations per subject are missing sporadically and completely at
random, this would not be a major problem. Many EMA designs give prompts at
random times, or at a random number of times, so methods for analyzing them
are constructed to allow non-shared measurement times and unequal coordinates.
This leads to the simple approach used in this paper, which is to treat
occasions of nonresponse to a prompt as though the prompt had never occurred
at all on those occasions (essentially listwise deletion, although on the
scale of the observation, not the entire participant). This approach assumes
that values of $M(t)$ are missing completely at random, which is not ideal.

The situation is more serious if participants drop out for long periods of
time, and especially if this might occur for reasons that relate to the
variables of interest (e.g., because of being embarrassed or discouraged by a
relapse). If the composition of the sample differs at different points in the
time interval of interest, then it may be difficult to interpret the overall
trajectory. About 14% of the participants in our analysis sample provided less
than 20 observations, so even in the absence of informative missingness, their
interpolated $M(t)$ trajectories may not closely match what they would have
been if more data had been provided.

Multiple imputation can be beneficial for handling missing values in
longitudinal data (see Schafer Graham, 2002; Lee ., 2016). More research is
required regarding how to use multiple imputation for functional regression
(Ciarleglio ., 2020), which is likely to be more challenging than in the more
classic case of a few waves of data at fixed times. If informative dropout is
suspected, another possibility might be to jointly model the longitudinal
processes together with a survival process for dropout, possibly treated as a
competing risk for the main outcome of interest (see X. Liu, 2014). However,
this also is likely to be challenging to implement in a functional mediation
context.

In addition to an overall indirect effect integrated over time, Lindquist
(2012) also estimated the product, $\alpha_{X}(t)\beta_{M}$(t), for a grid of
times $t$, leading to a time-specific indirect effect. In principle, bootstrap
tests could be done to test whether individual time-specific indirect effects
are statistically significant, although Lindquist recommends a false discovery
rate adjustment in this case, but this may lead to severe power limitations.
We did not consider time-specific indirect effects here. Functional regression
coefficients can be somewhat difficult to interpret on a time-specific basis
with longitudinal data, at least without additional model assumptions. The
mediator values at different time points are not independent, so that the most
statistically predictive time period need not be the most causally important
time period (see Dziak ., 2019; Green Popham, 2016; de Stavola ., 2022; Gao
Kowal, 2022). Furthermore, even without intensive longitudinal data,
researchers face the challenge of time-varying confounding (VanderWeele
Tchetgen Tchetgen, 2017; de Stavola ., 2022). That is, if $Y$ had also been
time-varying, then earlier values of $Y$ could affect later values of both $M$
and $Y$, becoming a time-varying confounder which would violate Assumption 4.
We focus on testing the overall effect of the mediator process on a single
later time point $Y$ in order to avoid these issues for now, although they
will be important to address if functional mediation is to become more widely
used for causal inference.

Similar to Lindquist (2012) and Zeng . (2021), we assumed no treatment-
mediator interaction in order to identify direct and indirect effects.
Lindquist (2012) does explore an alternative instrumental variables approach,
which requires alternative assumptions. However, this alternative might be
more complicated to extend to models that do not use the identity link. As an
alternative to natural direct and indirect effects, the use of interventional
(see VanderWeele Tchetgen Tchetgen, 2017) or organic (see Lok Bosch, 2020)
effects might make causal inference in the presence of time-varying
confounders more manageable.

We followed Lindquist (2012) in taking a frequentist approach, but a Bayesian
approach to mediation analysis may also be useful (Miočević ., 2018; K. Liu .,
2020). In particular, Huang Yuan (2017) described a Bayesian analogue to the
time-varying mediation method of Lindquist (2012) in the case where the
outcome was time-varying and normally distributed, and the ideas could
presumably be adjusted to the case of a point and/or binary distal outcome.

Finally, we have not considered the possibility of periodic functions, such as
diurnal or seasonal rhythms, as indexes for the varying coefficients in place
of linear time. There are some settings in which such features may be
theoretically and practically important (Adam ., 2017).

###  7.2 Future Directions

Our simulation and example have considered a single time-varying variable as
the mediator. Multiple time-varying mediators may be of interest, and options
for such models are briefly reviewed in Appendix 2.

We have only considered scalar treatment variables. However, another option
would be to have a time-varying exposure $X(t)$ and/or time-varying response
$Y(t)$ in addition to the time-varying mediator $M(t)$. This was explored by
Zhao . (2018) who provided software to fit either concurrent or historical
versions of the model. The concurrent approach could be done by fitting a TVEM
model to predict the mediator, and another TVEM model to predict the response,
and then multiplying the appropriate coefficients pointwise (without
integration). The historical approach is able to consider lagged effects
either within a sliding window (see also Kim ., 2011) or for the entire past.
The concurrent approach is simpler but could miss important information, by
only looking at essentially cross-sectional relationships and ignoring lagged
ones (see Maxwell ., 2011). An alternative, highly flexible but still somewhat
parametric, approach to time-varying mediation could use differential or
integral equations to represent dynamic models in continuous time (Albert .,
2019).

For point treatment and time-varying mediator and response, Cai . (2022) use a
two-step method beginning with a TVEM-like (function-on-scalar) regression to
predict $M$ from $X$, but follow with another function-on-function regression
instead of a scalar-on-function regression when predicting $Y$ from $M$ and
$X$. Cai . (2022) use local polynomial estimation (Fan Zhang, 2000) instead of
splines to handle the estimation of the semiparametric time-varying effect
coefficients, but either approach seems viable.

Functional mediation analysis is a novel method, relevant for both
methodologists and applied researchers interested in understanding how an
intensively measured longitudinal variable can mediate the relationship
between treatment and study outcomes. In the presented approach, we developed
functional mediation, a method for estimating a smooth, nonparametric
trajectory from treatment to an intensively measured longitudinal mediator
using TVEM (Tan ., 2012), as well as estimating a smooth, nonparametric
trajectory from an intensively measured longitudinal mediator to a binary end-
of-study outcome using functional regression (Ramsay Silverman, 2005). More
importantly, we provided a mathematical definition of the causal direct and
indirect effect that combined TVEM and functional regression in terms of
potential outcomes and provided guidance on how to estimate effects using the
R package `funmediation` as well as how to interpret effects.

Although our approach builds on prior work examining TVEM relationships in
regression (Tan ., 2012) and functional mediation (Lindquist, 2012), the
present work is innovative in that prior work has not considered a binary,
end-of-study outcome measured at only one time point. Binary end-of-study
outcomes are often a primary interest in intervention studies, such as smoking
cessation, weight maintenance, or abstinence from substance use at follow-up.
Despite binary outcomes being a common in intervention studies, methods to
evaluate the combination of intensively measured longitudinal data as well as
an end of study outcome following treatment are only beginning to be
developed. Functional mediation analysis, as presented here is only one of
hopefully many future methods that will be developed to model nonparametric
longitudinal processes as mediators to better understand how treatment impacts
behavior in a causal mediation framework.

## Acknowledgements

Research was supported by National Institutes of Health grants P50 DA039838
from the National Institute of Drug Abuse and 1R01 CA229542-01 from the
National Cancer Institute and the NIH Office of Behavioral and Social Science
Research. Content is solely the responsibility of the authors and does not
necessarily represent the official views of the funding institutions.

## Declaration of interest statement

The authors report there are no competing interests to declare.

## References

  * Adam . (2017) Adam2017Adam, EK., Quinn, ME., Tavernier, R., McQuilland, MT., Dahlkee, KA. Gilbert, KE.  2017\.  Diurnal Cortisol Slopes and Mental and Physical Health Outcomes: A Systematic Review and Meta-analysis Diurnal cortisol slopes and mental and physical health outcomes: A systematic review and meta-analysis.  Psychoneuroendocrinology25–41.
  * Albert . (2019) Albert2018Albert, JM., Li, Y., Sun, J., Woyczynski, WA. Nelson, S.  2019\.  Continuous-time causal mediation analysis Continuous-time causal mediation analysis.  Statistics in Medicine384334–4347.
  * Baker . (2016) Baker2016Baker, TB., Piper, ME., Stein, JH., Smith, SS., Bolt, DM., Fraser, DL. Fiore, MC.  2016\.  Effects of nicotine patch vs. varenicline vs. combination nicotine replacement therapy on smoking cessation at 26 weeks: A randomized clinical trial Effects of nicotine patch vs. varenicline vs. combination nicotine replacement therapy on smoking cessation at 26 weeks: A randomized clinical trial.  JAMA3154371-379.
  * Bind . (2016) BindBind, MA., Vanderweele, TJ., Coull, BA. Schwartz, JD.  2016\.  Causal mediation analysis for longitudinal data with exogenous exposure Causal mediation analysis for longitudinal data with exogenous exposure.  Biostatistics171122-134.
  * Cai . (2022) Xizhen2020ArxivCai, X., Coffman, DL., Piper, ME. Li, R.  2022\.  Estimation and Inference for the Mediation Effect in a Time-varying Mediation Model Estimation and inference for the mediation effect in a time-varying mediation model.  BMC Medical Research Methodology22113.
  * Canty Ripley (2020) boot1Canty, A. Ripley, BD.  2020\.  boot: Bootstrap R (S-Plus) Functions boot: Bootstrap R (S-Plus) functions [].  R package version 1.3-25
  * Chen . (2016) Chen2016Chen, YH., Mukherjee, B., Ferguson, KK., Meeker, JD. VanderWeele, TJ.  2016\.  Mediation formula for a binary outcome and a time-varying exposure and mediator, accounting for possible exposure-mediator interaction Mediation formula for a binary outcome and a time-varying exposure and mediator, accounting for possible exposure-mediator interaction.  American Journal of Epidemiology184157–159.
  * Ciarleglio . (2020) ciarleglioCiarleglio, A., Petkova, E. Harel, O.  2020\.  Multiple imputation in functional regression with applications to EEG data in a depression study. Multiple imputation in functional regression with applications to eeg data in a depression study.  arXiv:2001.08175
  * Daniel . (2015) Daniel2015Daniel, RM., de Stavola, BL., Cousens, SN. Vansteelandt, S.  2015\.  Causal Mediation Analysis with Multiple Mediators Causal mediation analysis with multiple mediators.  Biometrics1-14.
  * Davison Hinkley (1997) boot2Davison, AC. Hinkley, DV.  1997\.  Bootstrap Methods and Their Applications Bootstrap methods and their applications.  CambridgeCambridge University Press.
  * de Stavola . (2022) DeStavola2021de Stavola, B., Herle, M. Pickles, A.  2022\.  Framing Causal Questions in Life Course Epidemiology Framing causal questions in life course epidemiology.  Annual Review of Statistics and Its Application91223–248.
  * Derkach . (2019) Derkach2019Derkach, A., Pfeiffer, RM., Chen, T. Sampson, JN.  2019\.  High dimensional mediation analysis with latent variables High dimensional mediation analysis with latent variables.  Biometrics75745–756.
  * Doretti . (2021) Doretti2021Doretti, M., Raggi, M. Stanghellini, E.  2021\.  Exact parametric causal mediation analysis for a binary outcome with a binary mediator Exact parametric causal mediation analysis for a binary outcome with a binary mediator.  Statistical Methods and Applications.  https://doi.org/10.1007/s10260-021-00562-w
  * Dziak . (2021) package_tvemDziak, JJ., Coffman, DL., Li, R., Litson, K. Chakraborti, Y.  2021\.  tvem R package (Time-Varying Effect Models) tvem R package (time-varying effect models) [].  https://cran.r-project.org/web/packages/tvem/tvem.pdf
  * Dziak . (2019) Dziak2019Dziak, JJ., Coffman, DL., Reimherr, M., Petrovich, J., Li, R., Shiffman, S. Shiyko, MP.  2019\.  Scalar-on-function regression for predicting distal outcomes from intensively gathered longitudinal data: Interpretability for applied scientists Scalar-on-function regression for predicting distal outcomes from intensively gathered longitudinal data: Interpretability for applied scientists.  Statistics Surveys13150–180.
  * Eilers Marx (1996) EilersMarx1996Eilers, PHC. Marx, BD.  1996\.  Flexible smoothing with B-splines and penalties Flexible smoothing with B-splines and penalties.  Statistical Science11289 – 121.
  * Fan Zhang (2000) Fan2000Fan, J. Zhang, JT.  2000\.  Two-step estimation of functional linear models with applications to longitudinal data Two-step estimation of functional linear models with applications to longitudinal data.  Journal of the Royal Statistical Society: Series B (Statistical Methodology)622303-322.
  * Feingold . (2019) Feingold2019Feingold, A., MacKinnon, DP. Capaldi, DM.  2019\.  Mediation analysis with binary outcomes: Direct and indirect effects of proalcohol influences on alcohol use disorders Mediation analysis with binary outcomes: Direct and indirect effects of proalcohol influences on alcohol use disorders.  Addictive Behaviors9426-35.
  * Fritz (2014) Fritz2014Fritz, MS.  2014\.  An Exponential Decay Model for Mediation An exponential decay model for mediation.  Prevention Science15611-622.
  * Fritz . (2016) Fritz2016Fritz, MS., Kenny, DA. MacKinnon, DP.  2016\.  The Combined Effects of Measurement Error and Omitting Confounders in the Single-Mediator Model The combined effects of measurement error and omitting confounders in the single-mediator model.  Multivariate Behavioral Research515681-697.
  * Gao Kowal (2022) Gao2022Gao, Y. Kowal, DR.  2022\.  Bayesian adaptive and interpretable functional regression for exposure profiles. Bayesian adaptive and interpretable functional regression for exposure profiles.  https://arxiv.org/abs/2203.00784 arXiv:1903.04697
  * Geldhof . (2017) Geldhof2017Geldhof, GJ., Anthony, KP., Selig, JP. Mendez-Luck, CA.  2017\.  Accommodating binary and count variables in mediation Accommodating binary and count variables in mediation.  International Journal of Behavioral Development42300–308.
  * Gelman Vehtari (2021) Gelman2021Gelman, A. Vehtari, A.  2021\.  What are the Most Important Statistical Ideas of the Past 50 Years? What are the most important statistical ideas of the past 50 years?  Journal of the American Statistical Association1165362087-2097.
  * Goldsmith . (2011) GoldsmithBobb2011Goldsmith, J., Bobb, J., Crainiceanu, CM., Caffo, B. Reich, D.  2011\.  Penalized functional regression Penalized functional regression.  Journal of Computational and Graphical Statistics20830–851.
  * Goldsmith . (2020) package_refundGoldsmith, J., Scheipl, F., Huang, L., Wrobel, J., Di, C., Gellar, J.Staicu, AM.  2020\.  refund R package (Regression with Functional Data) refund R package (regression with functional data) [].  https://cran.r-project.org/web/packages/refund/refund.pdf
  * Green Popham (2016) Green2016Green, MJ. Popham, F.  201612\.  Life course models: improving interpretation by consideration of total effects Life course models: improving interpretation by consideration of total effects.  International Journal of Epidemiology4631057-1062.
  * Hastie Tibshirani (1993) Hastie1993Hastie, T. Tibshirani, R.  1993\.  Varying-coefficient models Varying-coefficient models.  Journal of the Royal Statistical Society, B55757–796.
  * Hernán Robins (2020) HernanWhatIfHernán, MA. Robins, JM.  2020\.  Causal Inference: What If Causal inference: What if.  Boca Raton: Chapman & Hall/CRC.
  * Huang Yuan (2017) Huang2017Huang, J. Yuan, Y.  2017\.  Bayesian Dynamic Mediation Analysis Bayesian dynamic mediation analysis.  Psychological Methods22667–686.
  * Imai . (2010) Imai2010aImai, K., Keele, L. Yamamoto, T.  2010\.  Identification, Inference and Sensitivity Analysis for Causal Mediation Effects Identification, inference and sensitivity analysis for causal mediation effects.  Statistical Science2551–71.  10.1214/10-STS321
  * Kim . (2011) Kim2011Kim, K., Şentürk, D. Li, R.  2011\.  Recent history of functional linear models for sparse longitudinal data Recent history of functional linear models for sparse longitudinal data.  Journal of Statistical Planning and Inference14141554-1566.
  * Le Foll . (2022) LeFollLe Foll, B., Piper, ME., Fowler, CD., Tonstad, S., Bierut, L., Lu, L.Hall, WD.  2022\.  Tobacco and nicotine use Tobacco and nicotine use.  Nature Reviews Disease8Primers 19.
  * Lee . (2016) Lee2016Lee, KJ., Roberts, G., Doyle, LW., Anderson, PJ. Carlin, JB.  2016\.  Multiple imputation for missing data in a longitudinal cohort study: a tutorial based on a detailed case study involving imputation of missing outcome data Multiple imputation for missing data in a longitudinal cohort study: a tutorial based on a detailed case study involving imputation of missing outcome data.  International Journal of Social Research Methodology19575-591.
  * Liang Zeger (1986) LiangZeger1986Liang, KY. Zeger, SL.  1986\.  Longitudinal data analysis using generalized linear models Longitudinal data analysis using generalized linear models.  Biometrika73113-22.
  * Lindquist (2012) Lindquist2012Lindquist, MA.  2012\.  Functional causal mediation analysis with an application to brain connectivity Functional causal mediation analysis with an application to brain connectivity.  Journal of the American Statistical Association1071297–1309.
  * Lindquist McKeague (2009) Lindquist2009Lindquist, MA. McKeague, IW.  2009\.  Logistic regression with Brownian-like predictors Logistic regression with Brownian-like predictors.  Journal of the American Statistical Association1041575–1585.
  * K. Liu . (2020) KuanLiu2020Liu, K., Saarela, O., Feldman, BM. Pullenayegum, E.  2020\.  Estimation of causal effects with repeatedly measured outcomes in a Bayesian framework Estimation of causal effects with repeatedly measured outcomes in a bayesian framework.  Statistical Methods in Medical Research0962280219900362.
  * X. Liu (2014) Liu2014Liu, X.  2014\.  Joint modeling of longitudinal and survival data: new models, computing techniques and applications Joint modeling of longitudinal and survival data: new models, computing techniques and applications Penn State.  https://etda.libraries.psu.edu/catalog/23528
  * Loh . (2020) Loh2020Loh, WW., Moerkerke, B., Loeys, T., Poppe, L., Crombez, G. Vansteelandt, S.  2020\.  Estimation of Controlled Direct Effects in Longitudinal Mediation Analyses with Latent Variables in Randomized Studies Estimation of controlled direct effects in longitudinal mediation analyses with latent variables in randomized studies.  Multivariate Behavioral Research555763-785.
  * Lok Bosch (2020) Lok2020Lok, JJ. Bosch, RJ.  2020\.  Causal organic indirect and direct effects: closer to Baron and Kenny, with a product method for binary mediators. Causal organic indirect and direct effects: closer to Baron and Kenny, with a product method for binary mediators.  arXiv:1903.04697
  * MacKinnon (2008) MacKinnon2008MacKinnon, DP.  2008\.  Introduction to statistical mediation analysis Introduction to statistical mediation analysis.  Routledge.
  * Magidson . (2015) Magidson2015Magidson, JF., Blashill, AJ., Safren, SA. Wagner, GJ.  2015\.  Depressive symptoms, lifestyle structure, and ART adherence among HIV-infected individuals: a longitudinal mediation analysis Depressive symptoms, lifestyle structure, and ART adherence among HIV-infected individuals: a longitudinal mediation analysis.  AIDS and Behavior1934–40.  10.1007/s10461-014-0802-3
  * Maxwell . (2011) Maxwell2011Maxwell, SE., Cole, DA. Mitchell, MA.  2011\.  Bias in Cross-Sectional Analyses of Longitudinal Mediation: Partial and Complete Mediation Under an Autoregressive Model Bias in cross-sectional analyses of longitudinal mediation: Partial and complete mediation under an autoregressive model.  Multivariate Behavioral Research465816-841.
  * Miočević . (2018) Miocevic2018Miočević, M., Gonzalez, O., Valente, MJ. MacKinnon, DP.  2018\.  A tutorial in Bayesian potential outcomes mediation analysis A tutorial in bayesian potential outcomes mediation analysis.  Structural Equation Modeling25121-136.
  * Moreno-Betancur Carlin (2018) Moreno2018Moreno-Betancur, M. Carlin, JB.  2018\.  Understanding Interventional Effects: A More Natural Approach to Mediation Analysis? Understanding interventional effects: A more natural approach to mediation analysis?  Epidemiology295614–617.
  * Nguyen . (2021) Nguyen2021clarifyingNguyen, TQ., Schmid, I. Stuart, EA.  2021\.  Clarifying causal mediation analysis for the applied researcher: Defining effects based on what we want to learn. Clarifying causal mediation analysis for the applied researcher: Defining effects based on what we want to learn.  Psychological Methods262255–271.
  * Park . (2018) ParkPark, S., Steiner, PM. Kaplan, D.  2018\.  Identification and Sensitivity Analysis for Average Causal Mediation Effects with Time-Varying Treatments and Mediators: Investigating the Underlying Mechanisms of Kindergarten Retention Policy Identification and sensitivity analysis for average causal mediation effects with time-varying treatments and mediators: Investigating the underlying mechanisms of kindergarten retention policy.  Psychometrika832298-320.
  * Pinto Pereira . (2020) PintoPereira2020Pinto Pereira, SM., De Stavola, BL., Rogers, NT., Hardy, R., Cooper, R. Power, C.  2020\.  Adult obesity and mid-life physical functioning in two British birth cohorts: investigating the mediating role of physical inactivity Adult obesity and mid-life physical functioning in two british birth cohorts: investigating the mediating role of physical inactivity.  International Journal of Epidemiology49845–856.
  * Qian . (2019) Qian2019Qian, T., Yoo, H., Klasnja, P., Almirall, D. Murphy, SA.  2019\.  Estimating time-varying causal excursion effect in mobile health with binary outcomes. Estimating time-varying causal excursion effect in mobile health with binary outcomes.  arXiv:1906.00528
  * R Core Team (2020) RManualR Core Team.  2020\.  R: A Language and Environment for Statistical Computing R: A language and environment for statistical computing [].  Vienna, Austria.  https://www.R-project.org/
  * Ramsay Silverman (2005) RamsaySilverman2005Ramsay, JO. Silverman, BW.  2005\.  Functional Data Analysis Functional data analysis (2 ).  New YorkSpringer.
  * Rijnhart . (2017) RijnhartTwisk2017Rijnhart, JJM., Twisk, JW., Chinapaw, MJM., de Boer, MR. Heymans, MW.  2017\.  Comparison of methods for the analysis of relatively simple mediation models Comparison of methods for the analysis of relatively simple mediation models.  Contemporary Clinical Trials Communications7130-135.
  * Rijnhart . (2019) RijnhartTwisk2019Rijnhart, JJM., Twisk, JWR., Eekhout, I. Heymans, MW.  2019\.  Comparison of logistic-regression based methods for simple mediation analysis with a dichotomous outcome variable Comparison of logistic-regression based methods for simple mediation analysis with a dichotomous outcome variable.  BMC Medical Research Methodology19.
  * Rijnhart . (2021) RijnhartValente2021Rijnhart, JJM., Valente, MJ., MacKinnon, DP., Twisk, JWR. Heymans, MW.  2021\.  The use of traditional and causal estimators for mediation models with a binary outcome and exposure-mediator interaction The use of traditional and causal estimators for mediation models with a binary outcome and exposure-mediator interaction.  Structural Equation Modeling: A Multidisciplinary Journal345-355.
  * Robins (2003) Robins2003Robins, JM.  2003\.  Semantics of causal DAG models and the identiﬁcation of direct and indirect effects Semantics of causal DAG models and the identiﬁcation of direct and indirect effects.  Highly Structured Stochastic Systems70–81.
  * Rubin (1980) Rubin1980Rubin, D.  1980\.  Randomization analysis of experimental data: The Fisher randomization test (comment) Randomization analysis of experimental data: The Fisher randomization test (comment).  Journal of the American Statistical Association75591–593.
  * Schafer Graham (2002) Schafer2002Schafer, JL. Graham, JW.  2002\.  Missing data: our view of the state of the art Missing data: our view of the state of the art.  Psychological Methods7147-77.
  * Shiffman . (2008) Shiffman2008Shiffman, S., Stone, AA. Hufford, MR.  2008\.  Ecological momentary assessment Ecological momentary assessment.  Annual Review of Clinical Psychology41–32.
  * Sullivan . (2021) Sullivan2021Sullivan, AJ., Gunzler, DD., Morris, N. VanderWeele, TJ.  2021\.  Longitudinal Mediation Analysis with Latent Growth Curves. Longitudinal mediation analysis with latent growth curves.  arXiv:2103.05765
  * Tan . (2012) Tan2012Tan, X., Shiyko, MP., Li, R., Li, Y. Dierker, L.  2012\.  A time-varying effect model for intensive longitudinal data A time-varying effect model for intensive longitudinal data.  Psychological Methods1761–77.
  * Vafaie Kober (2022) VafaieVafaie, N. Kober, H.  202206\.  Association of Drug Cues and Craving With Drug Use and Relapse: A Systematic Review and Meta-analysis Association of Drug Cues and Craving With Drug Use and Relapse: A Systematic Review and Meta-analysis.  JAMA Psychiatry.  https://doi.org/10.1001/jamapsychiatry.2022.1240 10.1001/jamapsychiatry.2022.1240
  * VanderWeele (2015) VanderWeele2015VanderWeele, TJ.  2015\.  Explanation in causal inference: Methods for mediation and interaction Explanation in causal inference: Methods for mediation and interaction.  New YorkOxford University Press.
  * VanderWeele Tchetgen Tchetgen (2017) VanderWeele2017VanderWeele, TJ. Tchetgen Tchetgen, EJ.  2017\.  Mediation analysis with time varying exposures and mediators Mediation analysis with time varying exposures and mediators.  J. R. Statist. Soc. B79917–938.
  * VanderWeele Vansteelandt (2009) VanderWeeleVansteelandt2009VanderWeele, TJ. Vansteelandt, S.  2009\.  Conceptual issues concerning mediation, interventions, and composition Conceptual issues concerning mediation, interventions, and composition.  Statistics and Its Interface2457-468.
  * VanderWeele Vansteelandt (2010) VanderWeeleVansteelandt2010VanderWeele, TJ. Vansteelandt, S.  2010\.  Odds ratios for mediation analysis for a dichotomous outcome Odds ratios for mediation analysis for a dichotomous outcome.  American Journal of Epidemiology172121339-1348.
  * Vansteelandt Daniel (2017) Vansteelandt2017Vansteelandt, S. Daniel, RM.  2017\.  Interventional effects for mediation analysis with multiple mediators Interventional effects for mediation analysis with multiple mediators.  Epidemiology28258–265.
  * Vansteelandt . (2019) Vansteelandt2019Vansteelandt, S., Linder, M., Vandenberghe, S., Steen, J. Madsen, J.  2019\.  Mediation of time-to-event endpoints accounting for repeatedly measured mediators subject to time-varying confounding Mediation of time-to-event endpoints accounting for repeatedly measured mediators subject to time-varying confounding.  Statistics in Medicine384828-4840.
  * Wang . (2013) Wang2013Wang, W., Nelson, S. Albert, JM.  2013\.  Estimation of causal mediation effects for a dichotomous outcome in multiple-mediator models using the mediation formula Estimation of causal mediation effects for a dichotomous outcome in multiple-mediator models using the mediation formula.  Statistics in Medicine324211–4228.
  * Wood (2003) Wood2003Wood, SN.  2003\.  Thin plate regression splines Thin plate regression splines.  Journal of the Royal Statistical Society, B6595-114.
  * Wood (2017) mgcvBookWood, SN.  2017\.  Generalized Additive Models: An Introduction with R (2nd edition) Generalized additive models: An introduction with R (2nd edition).  New YorkChapman and Hall/CRC.
  * Wood (2021) package_mgcvWood, SN.  2021\.  mgcv R package (Mixed GAM Computation Vehicle with Automatic Smoothness Estimation) mgcv R package (Mixed GAM Computation Vehicle with Automatic Smoothness Estimation) [].  https://cran.r-project.org/web/packages/mgcv/index.html
  * Zeng . (2021) ZengFanLi2021Zeng, S., Rosenbaum, S., Archie, E., Alberts, S. Li, F.  2021\.  Causal Mediation Analysis for Sparse and Irregular Longitudinal Data. Causal mediation analysis for sparse and irregular longitudinal data.  arXiv:2007.01796
  * Zhao . (2020) Zhao2020Zhao, Y., Lindquist, MA. Caffo, BS.  2020\.  Sparse principal component based high-dimensional mediation analysis Sparse principal component based high-dimensional mediation analysis.  Computational Statistics and Data Analysis142106835.
  * Zhao . (2018) zhao2018functionalZhao, Y., Luo, X., Lindquist, M. Caffo, B.  2018\.  Functional Mediation Analysis with an Application to Functional Magnetic Resonance Imaging Data. Functional mediation analysis with an application to functional magnetic resonance imaging data.
  * Zheng van der Laan (2017) ZhengZheng, W. van der Laan, M.  2017\.  Longitudinal Mediation Analysis with Time-varying Mediators and Exposures, with Application to Survival Outcomes Longitudinal mediation analysis with time-varying mediators and exposures, with application to survival outcomes.  J Causal Inference52.

## Appendix 1: Effect Decomposition with Log Link

### Proof for Log Link

Suppose that a log link function is being used, with the effect definitions
given in expression (1) and the assumptions given in section 3.3. The mediator
is continuous as given by model (2). The response is binary or count using
model (5), with the log function for $g$. The proof follows VanderWeele
Vansteelandt (2010), Lindquist (2012), and VanderWeele (2015).

By definition,

$\displaystyle\log(E(Y(x,\boldsymbol{M}({x^{\ast}}))|\boldsymbol{C}=\boldsymbol{c}))$
$\displaystyle=$ $\displaystyle\log\Big{(}\int
E\big{(}Y(x,M({x^{\ast}}))|M({x^{\ast}})=\boldsymbol{m},\boldsymbol{C}=\boldsymbol{c})P(M({x^{\ast}})=\boldsymbol{m}|\boldsymbol{C}=\boldsymbol{c})\mathrm{d{\bf
m}}\Big{)}$

By the consistency and composition assumptions and assumption 4, we can drop
the conditioning of $Y(x,\boldsymbol{M}({x^{\ast}})$ on
$\boldsymbol{M}({x^{\ast}})$ so that

$\displaystyle\log\Big{(}E(Y(x,\boldsymbol{M}({x^{\ast}}))|\boldsymbol{C}=\boldsymbol{c}))$
$\displaystyle=$ $\displaystyle\log\Big{(}\int
E(Y(x,\boldsymbol{M}({x^{\ast}}))|\boldsymbol{C}=\boldsymbol{c},\boldsymbol{M}({x^{\ast}})=\boldsymbol{m}))P(\boldsymbol{M}({x^{\ast}})=\boldsymbol{m}|\boldsymbol{C}=\boldsymbol{c})\mathrm{d{\bf
m}}\Big{)}$

By the first three assumptions, we can equate the potential outcomes with
conditional expectations so that

$\displaystyle\log(E(Y(x,\boldsymbol{M}({x^{\ast}})))$ $\displaystyle=$
$\displaystyle\log\Big{(}\int
E(Y|\boldsymbol{C}=\boldsymbol{c},X=x,\boldsymbol{M}=\boldsymbol{m}))P(\boldsymbol{M}=\boldsymbol{m}|\boldsymbol{C}=\boldsymbol{c},x={x^{\ast}})\mathrm{d{\bf
m}}\Big{)}$

Then using model (5) in which there is no interaction between the treatment
and the mediator, the above becomes

$\displaystyle\log\Bigg{(}\int\exp\Big{(}\beta_{0}+\boldsymbol{c}\boldsymbol{\beta}_{C}+\beta_{X}x+\int\beta(t)\boldsymbol{m}(t)\mathrm{dt}\Big{)}P\Big{(}\boldsymbol{M}=\boldsymbol{m}|\boldsymbol{C}=\boldsymbol{c},X={x^{\ast}}\Big{)}\mathrm{d{\bf
m}}\Bigg{)}$ $\displaystyle=$
$\displaystyle\int\big{(}\beta_{0}+\boldsymbol{c}\boldsymbol{\beta}_{C}+\beta_{X}x\big{)}P(\boldsymbol{M}=\boldsymbol{m}|\boldsymbol{C}=\boldsymbol{c},X={x^{\ast}})\mathrm{d{\bf
m}}$
$\displaystyle+\int\Big{(}\int\beta(t)\boldsymbol{m}(t)\mathrm{dt}\Big{)}P(\boldsymbol{M}=\boldsymbol{m}|\boldsymbol{C}=\boldsymbol{c},X={x^{\ast}})\mathrm{d{\bf
m}}$ (7) $\displaystyle=$
$\displaystyle\beta_{0}+\boldsymbol{c}\boldsymbol{\beta}_{C}+\beta_{X}x+E\Big{(}\int\beta_{M}(t)\boldsymbol{m}(t)\mathrm{dt}\Big{|}{x^{\ast}}\Big{)}.$

Using model (3),

$\displaystyle
E\big{(}\int\beta_{M}(t)\boldsymbol{m}(t)\mathrm{dt}\Big{|}\boldsymbol{c},{x^{\ast}}\big{)}$
$\displaystyle=\int\beta_{M}(t)E\big{(}\boldsymbol{m}(t)|\boldsymbol{c},{x^{\ast}}\big{)}\mathrm{dt}$
$\displaystyle=\int\alpha_{0}+\boldsymbol{c}\boldsymbol{\alpha}_{C}(t)+\alpha_{X}(t){x^{\ast}}\mathrm{dt}$

Therefore,

$\displaystyle\log(E(Y(x,\boldsymbol{M}({x^{\ast}})|\boldsymbol{c}))$
$\displaystyle=\beta_{0}+\boldsymbol{c}\boldsymbol{\beta}_{C}+\beta_{X}X+\int\beta_{t}\big{(}\alpha_{0}+\boldsymbol{c}\boldsymbol{\alpha}_{C}(t)+\alpha_{X}(t){x^{\ast}}\big{)}\mathrm{dt}.$

Now we can compute the total effect, natural direct effect, and natural
indirect effect. First,

$\displaystyle\log\mathrm{TE}=$ $\displaystyle\log
E\Big{(}Y\big{(}1,\boldsymbol{M}(1)\big{)}\Big{)}-\log
E\Big{(}Y\big{(}0,\boldsymbol{M}(0)\big{)}\Big{)}$
$\displaystyle=\Big{(}\beta_{0}+\boldsymbol{c}\boldsymbol{\beta}_{c}+\beta_{X}+\int\beta_{M}(t)\big{(}\alpha_{0}(t)+\alpha_{X}(t)\big{)}\mathrm{dt}\Big{)}-\Big{(}\beta_{0}+\boldsymbol{c}\boldsymbol{\beta}_{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}C}+\int\beta_{M}(t)\alpha_{0}(t)\mathrm{dt}\Big{)}$
$\displaystyle=\beta_{X}+\int\alpha_{X}(t)\beta_{M}(t)\mathrm{dt}.$

Next,

$\displaystyle\log\mathrm{DE}=$ $\displaystyle\log
E\Big{(}Y\big{(}1,\boldsymbol{M}(0)\big{)}\Big{)}-\log
E\Big{(}Y\big{(}0,\boldsymbol{M}(0)\big{)}\Big{)}$ $\displaystyle=$
$\displaystyle\Big{(}\beta_{0}+\boldsymbol{c}\boldsymbol{\beta}_{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}C}+\beta_{X}+\int\beta_{M}(t)\alpha_{X}(t)\mathrm{dt}\Big{)}-\Big{(}\beta_{0}+\boldsymbol{c}\boldsymbol{\beta}_{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}C}+\int\beta_{M}(t)\alpha_{X}(t)\mathrm{dt}\Big{)}=\beta_{X}.$

Last,

$\displaystyle\log\mathrm{IE}$ $\displaystyle=\log
E\Big{(}Y\big{(}0,\boldsymbol{M}(1)\big{)}\Big{)}-\log
E\Big{(}Y\big{(}0,\boldsymbol{M}(0)\big{)}\Big{)}$
$\displaystyle=\Big{(}\beta_{0}+\boldsymbol{c}\boldsymbol{\beta}_{C}+\int\beta_{M}(t)\big{(}\alpha_{0}+\alpha_{X}(t)\big{)}\mathrm{dt}\Big{)}-\Big{(}\beta_{0}+\boldsymbol{c}\boldsymbol{\beta}_{C}+\int\beta_{M}(t)\alpha_{0}(t)\mathrm{dt}\Big{)}$
$\displaystyle=\int\alpha_{X}(t)\beta_{M}(t)\mathrm{dt}.$

Comparing the above results, it is clear that the direct and indirect effects
add (on the log scale) to the total effects, as the definitions of these
effects would imply. The proof above would not work exactly for a logit link,
because the equivalent of (Proof for Log Link) would not hold; it is generally
not the case that $\mathrm{logit}(a)+\mathrm{logit}(b)=\mathrm{logit}(ab)$.

A further limitation is that, following Lindquist (2012), this proof treats
$\boldsymbol{M}(t)$ as fully observed, so that it can be replaced with its
expectation $\alpha_{0}+\alpha_{X}(t)X$ without any bias. In practice, the
sample estimates do not exactly decompose, presumably because of smoothing.

## Appendix 2: Extending the Model to Multiple Treatments or Multiple
Mediators

So far, we did not consider the possibility of more than one mediator process.
One way to handle multiple mediators is to study the joint mediated effect
through all of the mediators (see VanderWeele, 2015; Wang ., 2013), which
could be estimated by subtracting the natural direct effect from the total
effect. In this case, the natural direct effect represents the effect that
does not go through any of the mediators. If the no-interaction assumption is
made, then the controlled direct effect is assumed to equal the natural direct
effect, and thus the joint mediated effect may be estimated by subtracting the
controlled direct effect from the total effect. Of course, this approach does
not reveal which of the mediators, individually, is most important.

Separately estimating the effects of multiple treatment variables and multiple
mediator variables is also possible, although the number of parameters
required increases quickly. Suppose there are $p$ treatment variables, each
either dichotomous or numerical. For example, a set of $p$ binary treatment
variables might represent dummy codes for a categorical treatment variable
with $p+1$ levels, such as $p$ treatment groups being compared to a single
control group.

Suppose there are $q$ mediator variables. Then the model for the mediator
could be specified as

$E(\boldsymbol{M}(t))=\boldsymbol{A}_{0}(t)+\boldsymbol{A}(t)\mathbf{X}$

for the vector of time-varying mediators, where $\boldsymbol{A}$ is a
$q_{2}\times p$ matrix of functions of $t$, and the model for the outcome
could be specified as

$g(E(Y))=\beta_{0}+\sum_{k=1}^{p}\beta_{x_{k}}x_{k}+\sum_{\ell=1}^{q}\int_{t}\beta_{M_{\ell}}(t)M_{\ell}(t)\mathrm{dt}.$
(8)

If mediator, $\ell$, is not time-varying and does not have time-varying
effects, the integral $\int\beta_{M_{\ell}}(t)M_{\ell}(t)\mathrm{dt}$ is not
needed and could simply be replaced by a product, $\beta_{M_{\ell}}M_{\ell}$,
for some scalar, $\beta_{M_{\ell}}$. Likewise, in this case, the $\ell$th
column in $\mathbf{A}(t)$ would equal a constant $A_{\ell}$ across all
observation times.

In total, if there are $p$ treatment variables and $q$ mediator variables,
then there are $pq$ distinct indirect effect pathways. As a caveat, obtaining
identifiable estimates and useful confidence intervals from this model might
be difficult unless $p$ and $q$ are very small or the sample size is very
large. This is because the underlying model for $Y$ is essentially a
regression with many covariates (considerably larger than $pq$, depending on
the effective number of time points after smoothing) but only one observation
per subject for the outcome. One option for making this problem easier would
be to impose parametric assumptions, such as linear or quadratic trajectories,
for the mediator or for the $\alpha$ and $\beta$ coefficients. A decaying
exponential trajectory for the effects may be more interpretable than a
quadratic trajectory (Fritz, 2014). Alternatively, the effects of a mediator
might be constrained to be exactly cumulative, which amounts to its $\beta$
being constant over time (see Chen ., 2016).

The model shown in Equation (8) treats each mediator as a parallel pathway
from treatment to outcome, so that mediators never interact or affect each
other. If the number of mediators and number of measurement times were both
very small, then it might be feasible to go beyond this restriction. One could
then estimate the path-specific effects (Daniel ., 2015). However, this will
generally not be feasible or interpretable for functional data, because the
problems described earlier would be increased by the large number of pathways.
Models with large numbers of mediators, such as in genomics studies, might be
more challenging to combine with EMA data because of the resulting
multiplication of the number of parameters. However, it should be noted that
functional mediation models share some similarities with high-dimensional
mediation models, in that both involve some form of regularization or data
reduction to make the mediator tractable (Derkach ., 2019; Zhao ., 2020).
Therefore, future work in both areas might lead to further adaptations of
ideas from one to the other.

As an alternative to using multiple mediators, it might be beneficial to use
interventional direct and indirect effects (Chen ., 2016; Vansteelandt Daniel,
2017), which do not require assuming that no confounders of the mediator and
outcome have been affected by the treatment. A disadvantage is that
interventional direct and indirect effects do not necessarily sum to the total
effects. The controlled direct effect is a special case of an interventional
direct effect and thus, it is also identified in the presence of a mediator-
outcome confounder that has been affected by the treatment.

## Appendix 3: Analysis excluding participants with two or fewer EMA prompts
post-TQD

We reanalyzed our empirical example by excluding the 21 participants who had
two or fewer EMA prompts post-TQD. The overall conclusions remain unchanged.
The total effect estimate was not statistically significant (estimate =
-0.016, std. error = 0.136, $z=-0.119$, $p=0.906$). The indirect effect
estimate was statistically significant (estimate = 0.043, percentile method
95% bootstrap confidence interval = (0.005, 0.114)). The estimate of the
direct effect was not statistically significant (estimate = -0.067, std. error
= 0.139). The impact of excluding participants with greater than 97% missing
EMA responses is reflected in the total effect estimate but not in the direct
and indirect effect estimates. This could be interpreted as the proposed
technique’s strength in estimating functional mediation processes in spite of
missingness. However, the impact of missingness needs to be thoroughly
investigated before claiming this interpretation as scientific evidence of the
technique’s functionality.

## Supplementary Material: R Code

The code below was used to fit the model to the data set eos1. Further details
on the function arguments can be found in the vignette tutorial in the
funmediation R package.

        eos1 <- read.csv("eos1.csv")
        funresult <- funmediation(data=eos1,
                               treatment=NRT1,
                               mediator=WantToSmokeLst15minC,
                               outcome=abstDay28,
                               id=SubjectID,
                               time=timeseq,
                               binary_outcome = TRUE,
                               interpolate = FALSE,
                               nboot=500)

        print(funresult)
        plot(funresult, what_plot="pfr")
        plot(funresult, what_plot="coefs")
        plot(funresult, what_plot="tvem")

Table 1:  Simulation Results for TVEM Predicting Mediator from Treatment Sample size | 250 | 500 | 1000
---|---|---|---
Intercept function $\alpha_{0}(t)$
Bias | -0.0010 | -0.0002 | -0.0006
Mean squared error | 0.0105 | 0.0054 | 0.0029
Root mean squared error | 0.1022 | 0.0736 | 0.0537
Mean estimated SE | 0.1038 | 0.0739 | 0.0535
Pointwise coverage | 0.9461 | 0.9464 | 0.9460
Familywise coverage | 0.6422 | 0.6270 | 0.6276
Treatment effect function $\alpha_{X}(t)$
Bias | 0.0001 | 0.0008 | 0.0001
Mean squared error | 0.0203 | 0.0106 | 0.0057
Root mean squared error | 0.1423 | 0.1031 | 0.0758
Mean estimated SE | 0.1444 | 0.1032 | 0.0747
Pointwise coverage | 0.9485 | 0.9531 | 0.9510
Familywise coverage | 0.7170 | 0.6956 | 0.6412
Table 2:  Simulation Results for Functional Regression Predicting Outcome from Mediator Sample size | 250 | 500 | 1000
---|---|---|---
Intercept $\beta_{0}$
Bias | 0.0277 | 0.0236 | 0.0301
Mean squared error | 0.0592 | 0.0284 | 0.0145
Root mean squared error | 0.2432 | 0.1685 | 0.1203
Mean estimated SE | 0.2379 | 0.1668 | 0.1174
Coverage | 0.9492 | 0.9504 | 0.9488
Treatment $\beta_{X}$
Bias | -0.0161 | -0.0176 | -0.0241
Mean squared error | 0.0835 | 0.0410 | 0.0203
Root mean squared error | 0.2889 | 0.2024 | 0.1424
Mean estimated SE | 0.2833 | 0.1986 | 0.1399
Coverage | 0.9486 | 0.9434 | 0.9504
Mediator conditional effect function $\beta_{M}(t)$
Bias | -0.0186 | -0.0180 | -0.0219
Mean squared error | 0.1329 | 0.0651 | 0.0339
Root mean squared error | 0.3645 | 0.2552 | 0.1841
Pointwise coverage | 0.9603 | 0.9577 | 0.9552
Familywise coverage | 0.8598 | 0.8492 | 0.8364
Table 3:  Simulation Results for Indirect Effect Sample size | 250 | 500 | 1000
---|---|---|---
Performance of point estimate
Bias | 0.0220 | 0.0165 | 0.0160
Mean squared error | 0.0140 | 0.0070 | 0.0036
Root mean squared error | 0.1182 | 0.0839 | 0.0603
Mean estimated SE | 0.1325 | 0.0869 | 0.0593
Bootstrap coverage by method
Normal | 0.9648 | 0.9500 | 0.9340
Basic | 0.9682 | 0.9516 | 0.9348
Percentile | 0.9576 | 0.9512 | 0.9452
Bootstrap power by method
Normal | 0.2366 | 0.6350 | 0.9284
Basic | 0.1570 | 0.5654 | 0.9104
Percentile | 0.3704 | 0.6906 | 0.9420

$X$$M$$Y$$\alpha$$\beta$$\gamma$ |  | $X$$M(t)$$Y$$\alpha(t)$$\beta(t)$$\gamma$
---|---|---

Figure 1:  Simple (left) versus functional (right) mediation with single-time
treatment and outcome.

$X$$\boldsymbol{C}$$M(t)$$Y$$\alpha(t)$$\beta(t)$

Figure 2:  Mediation model with functional $M$, including baseline covariates
$\boldsymbol{C}$.

Figure 3:  Time-varying effect of varenicline (vs. nicotine patch only) on the
mediator, craving to smoke in the last 15 min. Left and right panes represent
time-varying intercept function and time-varying effect of treatment,
respectively. Upper and lower dotted lines represent pointwise approximate 95%
confidence intervals.

Figure 4:  Functional effect of the mediator process on CO-confirmed 7-day
point-prevalence abstinence 4-weeks post-TQD. Dotted lines represent
confidence intervals.

## Figure Captions

  1. 1.

Simple (left) versus functional (right) mediation with single-time treatment
and outcome.

  2. 2.

Mediation model with functional $M$, including baseline covariates
$\boldsymbol{C}$.

  3. 3.

Time-varying effect of varenicline (vs. nicotine patch only) on the mediator,
craving to smoke in the last 15 min. Left and right panes represent time-
varying intercept function and time-varying effect of treatment, respectively.
Upper and lower lines represent pointwise approximate 95% confidence
intervals.

  4. 4.

Functional effect of the mediator process on CO-confirmed 7-day point-
prevalence abstinence 4-weeks post-TQD. Dotted lines represent confidence
intervals.