# Mediation Analysis with Graph Mediator

Yixi Xu1 and Yi Zhao1
1Department of Biostatistics and Health Data Science
Indiana University School of Medicine

###### Abstract

This study introduces a mediation analysis framework when the mediator is a
graph. A Gaussian covariance graph model is assumed for graph representation.
Causal estimands and assumptions are discussed under this representation. With
a covariance matrix as the mediator, parametric mediation models are imposed
based on matrix decomposition. Assuming Gaussian random errors, likelihood-
based estimators are introduced to simultaneously identify the decomposition
and causal parameters. An efficient computational algorithm is proposed and
asymptotic properties of the estimators are investigated. Via simulation
studies, the performance of the proposed approach is evaluated. Applying to a
resting-state fMRI study, a brain network is identified within which
functional connectivity mediates the sex difference in the performance of a
motor task.

Keywords: Common diagonalization; Covariance regression; Decomposition method;
Gaussian covariance graph model

##  1 Introduction

Mediation analysis has been widely used in clinical and biomedical studies to
delineate the intermediate effect of a third variable, called mediator, in the
causal pathway between the exposure/treatment and the outcome. It helps
dissect the underlying causal mechanism by decomposing the effect of the
exposure on the outcome into the part through the mediator and the part not
through the mediator. Since the introduction of the classic Baron and Kenny
framework (Baron and Kenny,, 1986), mediation analysis has been extensively
studied over decades, including under the causal inference framework (see
review articles by Imai et al.,, 2010; VanderWeele,, 2015, 2016, among
others). With the advances of biological technologies, mediation analysis has
been extended to study cases of high-dimensional mediator (such as Huang and
Pan,, 2016; Chén et al.,, 2017; Derkach et al.,, 2019; Zhao and Luo,, 2022,
among others) and complex data output, including functional data (Lindquist,,
2012; Zhao et al.,, 2018; Zeng et al.,, 2021), time series data (Gu et al.,,
2014; Zhao and Luo,, 2019), network mediator (Zhao et al.,, 2022), image
mediator (Jiang and Colditz,, 2023; Chen and Zhou,, 2023), and so on. In this
study, we introduce a framework considering a graph as the mediator.

Our study is motivated by resting-state functional magnetic resonance imaging
(fMRI) experiments. A primary interest in resting-state fMRI is to portray
brain coactivation patterns, the so-called brain functional connectivity or
connectome, when the participant is at rest without any external stimulus.
These coactivation patterns are captured by the covariance matrices (or
correlation matrices after standardizing the data) of the fMRI signals
extracted from brain voxels or regions of interest (ROI), where each ROI is a
group of voxels defined by a chosen brain parcellation atlas (Friston,, 2011).
Considering a cognitive behavior, it is also of great interest to understand
the brain mechanism related to the divergence in the behavior under various
exposure/treatment conditions. Formulating it into a mediation analysis
framework, brain functional connectivity is considered as the mediator between
the exposure and behavior. Straightforward approaches include pairwise
mediation analysis, where mediation analysis is repeated for each pair of
functional connectivity followed by a multiple-testing correction, and high-
dimensional mediation analysis, where the mediator is the vectorization of the
upper triangular portion of the connectivity matrix. The former ignores the
dependence between the correlations and suffers from power deficiency due to
multiplicity and the latter disregards the structure property and the positive
definiteness of a covariance matrix.

In general, resting-state neuronal fluctuations measured by fMRI can be
considered as Gaussian processes (Lindquist,, 2008). The functional
connectivity matrix can thus be modeled via a Gaussian covariance graph model.
As a subclass of graphical models, the Gaussian covariance graph model defines
a correspondence between the graph and the correlation pattern embedded in the
covariance matrix (Richardson and Spirtes,, 2002; Chaudhuri et al.,, 2007). A
missing connection between two nodes in the graph is in concordance with zero
correlation between the two corresponding variables in the covariance matrix,
which also coincides with the marginal independence. In this study, we will
consider such a connectivity graph as the mediator to articulate the
intermediate effect of brain connectome on the pathway between exposure and
cognitive behavior. Figure 1 presents a conceptual diagram, where the graph is
the mediator (denoted as $\mathcal{G}$). Closely related to this concept, Zhao
et al., (2022) recently proposed a Bayesian network mediation analysis. In
their framework, stochastic block models (SBMs) were employed on individual
elements in the adjacency matrix of brain regions and the vectorized model
parameters were treated as the latent mediators to integratively impose
feature selection and effect estimation. As the network mediator was
represented by the adjacency matrix, only symmetry was imposed on the matrix
structure and less stringent conditions were required for modeling. In this
study, instead of assuming underlying models on the adjacency matrix, it is
proposed to directly model the covariance matrix as the mediator. It extends
the linear structural equation modeling (LSEM) framework with a network
mediator and at the same time preserves the positive definite property of a
covariance matrix in parsimonious modeling. Here, we would also like to
distinguish the proposed framework from the case of image mediator considered
by Jiang and Colditz, (2023) and Chen and Zhou, (2023), where an image is
(multidimensional) scalar output with spatial information.

Regression models for covariance matrix outcomes have been long before studied
to capture the heterogeneity in the covariance structure across
individuals/populations. In order to maintain the positive definiteness in
modeling, example parametric approaches include modeling the covariance matrix
as a linear combination of symmetric matrices (Anderson,, 1973), modeling the
logarithmic transformation of matrix elements (Chiu et al.,, 1996), modeling
the covariance matrix via a quadratic link function (Hoff and Niu,, 2012;
Seiler and Holmes,, 2017) or a linear combination of similarity matrices (Zou
et al.,, 2017) of the covariates, and decomposition approaches based on the
common diagonalization assumption via either eigendecomposition (Flury,, 1984;
Boik,, 2002; Hoff,, 2009; Franks and Hoff,, 2019; Zhao et al.,, 2021) or
Cholesky decomposition (Pourahmadi et al.,, 2007). In this study, we adapt the
covariance regression framework introduced by Zhao et al., (2021) to mediation
analysis. In their framework, orthogonal common linear projections were
assumed on the covariance matrices and data heteroscedasticity in the
projection space was revealed via a log-linear model. The advantages are to
preserve the positive definiteness of the covariance matrices, offer high
flexibility in modeling, and enable a network-level interpretation, where each
projection can be considered a subnetwork after proper thresholding or
sparsifying. Integrating this framework with LSEM, common linear projections
are assumed on the covariance matrices and the logarithmic transformation of
data variance in the projection space is considered as the mediator
measurement in the structural equation models. The objective is to
simultaneously identify the linear projections and the causal parameters.
Contributions of this study include  (i) to enrich mediation analysis
literature with graph mediator;  (ii) to extend covariance regression modeling
to SEM framework;  (iii) and to offer a better understanding of brain
mechanisms when studying the exposure effect on cognitive behaviors.

The rest of the manuscript is organized as follows. Section 2 discusses the
causal definitions and assumptions and the proposed parametric mediation model
when the mediator is a graph. Likelihood-based estimators are proposed for
estimating model parameters, an efficient algorithm is introduced for
computation, and a resampling approach is considered for inference. Asymptotic
properties of the estimators are investigated under regularity conditions. Via
simulation studies in Section 3, the performance of the proposed framework is
evaluated and compared with competing approaches. Section 4 applies the
proposed mediation analysis to data acquired from the National Consortium on
Alcohol and NeuroDevelopment in Adolescence (NCANDA) to dissect the mediation
effect of brain functional connectome on sex difference in the performance of
a motor task. Section 5 summarizes the manuscript with a discussion.

Figure 1: A conceptual mediation diagram with graph mediator ($\mathcal{G}$).
$X$: exposure; $Y$: outcome; $W$: confounding factor.

##  2 Model and Methods

###  2.1 Causal definitions and assumptions

Let $\mathcal{G}(x)=(\mathcal{V},\mathcal{E}(x))$ denote the graph object when
the exposure is at level $x$, where $\mathcal{V}$ is the set of vertices/nodes
and $\mathcal{E}(x)$ is the set of edges under the exposure level. Denote
$\mathscr{G}=\\{\mathcal{G}\\}$ as the set of such graph objects. A type of
graph, covariance graph, is considered in this study, which can be represented
by a covariance matrix embedded with association structure between vertices
(Chaudhuri et al.,, 2007). Denote
$\mathbf{M}(x)=(M_{1}(x),\dots,M_{p}(x))^{\top}\in\mathbb{R}^{p}$ as the
potential outcome of the $p$-dimensional random variables in the vertex set
$\mathcal{V}=\\{1,\dots,p\\}$ and
$\boldsymbol{\Sigma}(x)=(\sigma_{jk}(x))\in\mathbb{R}^{p\times p}$ as the
corresponding covariance matrix when the exposure is $x$. As the study focus
is on the covariance matrix, without loss of generality, the variable means
are assumed to be zero. Assuming $\mathbf{M}(x)$ follows a multivariate
Gaussian distribution, a Gaussian covariance graph model is considered, where
a missing edge between two vertices is equivalent to the marginal independence
between the corresponding random variables (Edwards,, 2012): for
$\forall~{}j,k\in\mathcal{V}$ and $j\neq k$,

$(j,k)\notin\mathcal{E}(x)\quad\Leftrightarrow\quad
M_{j}(x)~{}\rotatebox[origin={c}]{90.0}{$\models$}~{}M_{k}(x)\quad\Leftrightarrow\quad\sigma_{jk}(x)=\sigma_{kj}(x)=0.$
(1)

Let $Y(x,\mathcal{G}(x))\in\mathbb{R}$ denote the potential outcome of $Y$
when the exposure is $x$ and the mediator graph is $\mathcal{G}(x)$. Let
$\mathscr{X}$ denote the possible exposure values. Assume a binary exposure,
$\mathscr{X}=\\{0,1\\}$. The following defines the average total effect (ATE)
of the exposure on the outcome.

$\tau_{\text{ATE}}=\mathbb{E}\left\\{Y(1,\mathcal{G}(1))-Y(0,\mathcal{G}(0))\right\\}.$
(2)

With a graph mediator, the above total effect can be decomposed into direct
and indirect effects. Extending the framework in Imai et al., (2010), the
following defines the average indirect effect (AIE) under exposure level $x$
and the average direct effect (ADE) when the mediator graph is at level
$\mathcal{G}(x)$.

$\displaystyle\tau_{\text{AIE}}(x)$ $\displaystyle=$
$\displaystyle\mathbb{E}\left\\{Y(x,\mathcal{G}(1))-Y(x,\mathcal{G}(0))\right\\},\quad
x=0,1;$ (3) $\displaystyle\tau_{\text{ADE}}(x)$ $\displaystyle=$
$\displaystyle\mathbb{E}\left\\{Y(1,\mathcal{G}(x))-Y(0,\mathcal{G}(x))\right\\},\quad
x=0,1.$ (4)

The AIE quantifies the difference between the potential outcome corresponding
to the graph of $\mathcal{G}(1)$ and that under $\mathcal{G}(0)$ fixing the
exposure at level $x$. It is also called the natural indirect effect (Pearl,,
2001) or the pure indirect effect for $\tau_{\text{AIE}}(0)$ and total
indirect effect for $\tau_{\text{AIE}}(1)$ (Robins and Greenland,, 1992). The
ADE quantifies the exposure effect on the outcome not through the mediator. It
is easy to show that the ATE is the sum of AIE and ADE,

$\tau_{\text{ATE}}=\tau_{\text{AIE}}(x)+\tau_{\text{ADE}}(1-x),\quad x=0,1.$
(5)

The following discusses the identification assumptions of the causal estimands
defined above. For mediation analysis, the identification assumptions have
been extensively discussed under various settings in the literature (examples
include Robins and Greenland,, 1992; Pearl,, 2001; Imai et al.,, 2010;
VanderWeele,, 2015, and so on). In this study, we extend the assumptions in
Imai et al., (2010) to the scenario of a graph mediator.

Assumption A1 (SUTVA)

Stable unit treatment value assumption.

Assumption A2 (Positivity)

Let $X$ denote the actual treatment assignment, then $\mathbb{P}(X=x)>0$ for
$\forall~{}x\in\mathscr{X}$.

Assumption A3 (Ignorability)

Let $Y(x,\boldsymbol{\mathcal{G}})$ denote the potential outcome of $Y$ under
the exposure level $x$ and the mediator graph taking the value of
$\boldsymbol{\mathcal{G}}\in\mathscr{G}$ and $\mathbf{W}\in\mathbb{R}^{q}$
denote a vector of $q$-dimensional observed confounding factors.

$\left\\{Y(1,\boldsymbol{\mathcal{G}}),Y(0,\boldsymbol{\mathcal{G}}),\mathbf{M}(1),\mathbf{M}(0)\right\\}~{}\rotatebox[origin={c}]{90.0}{$\models$}~{}X~{}|~{}\mathbf{W},$
(6)

for any $\boldsymbol{\mathcal{G}}\in\mathscr{G}$.

Assumption A4 (Sequential ignorability)

$Y(x,\boldsymbol{\mathcal{G}})~{}\rotatebox[origin={c}]{90.0}{$\models$}~{}\mathbf{M}(x)~{}|~{}X=x,\mathbf{W},$
(7)

for any $x\in\mathscr{X}$ and $\boldsymbol{\mathcal{G}}\in\mathscr{G}$.

Assumption A5 (Parametric LSEM)

The following linear structural equation models (LSEMs) with independent
errors are assumed.

$\displaystyle\log(\boldsymbol{\theta}^{\top}\boldsymbol{\Sigma}(x)\boldsymbol{\theta})$
$\displaystyle=$ $\displaystyle\alpha_{0}+\alpha
x+\mathbf{W}^{\top}\boldsymbol{\phi}_{1}+\eta,$ (8) $\displaystyle
Y(x,\boldsymbol{\mathcal{G}})$ $\displaystyle=$
$\displaystyle\gamma_{0}+\gamma
x+\beta\log(\boldsymbol{\theta}^{\top}\boldsymbol{\Sigma}\boldsymbol{\theta})+\mathbf{W}^{\top}\boldsymbol{\phi}_{2}+\epsilon,$
(9)

where $\boldsymbol{\theta}\in\mathbb{R}^{p}$ is a projection vector with unit
$L_{2}$-norm (i.e. $\|\boldsymbol{\theta}\|_{2}=1$),
$\\{\alpha_{0},\alpha,\boldsymbol{\phi}_{1},\gamma_{0},\gamma,\beta,\boldsymbol{\phi}_{2}\\}$
are model coefficients, and $\eta$ and $\epsilon$ are independent model errors
with mean zero.

The SUTVA assumption (Rubin,, 1980) assumes that the treatment assignment
regime remains the same across units and the outcome of one unit is not
influenced by the treatment assignment of other units. In the application,
this no interference assumption holds as the imaging scan and cognition
evaluation were performed for each individual independently. The positivity
assumption is standard in causal inference literature when possible treatment
assignments are finite. Assumption A3 assumes that there is no unmeasured
mediator-exposure or outcome-exposure confounding given the observed
covariates. In the application study, sex is considered as the preceding
exposure, which is assumed to be naturally randomized among the population of
the same age. Assumption A4 assumes the conditional independence between the
potential outcome of $Y$ and the potential outcome of $\mathbf{M}$ under the
so-called single-world intervention graphs (Richardson and Robins,, 2013),
which is a relaxation of the cross-world independence assumption (Robins and
Richardson,, 2010). Assumption 5 assumes parametric linear structural equation
models, where $\boldsymbol{\mathcal{G}}$ and $\boldsymbol{\Sigma}$ have a one-
to-one correspondence under the Gaussian covariance graph model. With
Assumptions A3–A5, it is sufficient to identify the ADE and AIE (Andrews and
Didelez,, 2020). Considering a Gaussian covariance graph model, the covariance
matrix ($\boldsymbol{\Sigma}$) has a one-to-one correspondence with the graph
object ($\mathcal{G}$). Under Assumption A5, the potential outcome of $Y$
under a multiple-worlds model can be expressed as

$\displaystyle Y(x,\mathcal{G}(x^{\prime}))$ $\displaystyle=$
$\displaystyle\gamma_{0}+\gamma x+\beta(\alpha_{0}+\alpha
x^{\prime}+\mathbf{W}^{\top}\boldsymbol{\phi}_{1}+\eta)+\mathbf{W}^{\top}\boldsymbol{\phi}_{2}+\epsilon$
$\displaystyle=$ $\displaystyle\gamma x+\alpha\beta
x^{\prime}+(\gamma_{0}+\alpha_{0}\beta+\beta\mathbf{W}^{\top}\boldsymbol{\phi}_{1}+\mathbf{W}^{\top}\boldsymbol{\phi}_{2})+\beta\eta+\epsilon.$

The following theorem gives the parametric representation of the causal
estimands.

######  Theorem 1.

Under Assumptions A1–A5,

$\displaystyle\tau_{\mathrm{ATE}}=\mathbb{E}\left\\{Y(1,\mathcal{G}(1))-Y(0,\mathcal{G}(0))\right\\}=\gamma+\alpha\beta;$
$\displaystyle\tau_{\mathrm{AIE}}(x)=\mathbb{E}\left\\{Y(x,\mathcal{G}(1))-Y(x,\mathcal{G}(0))\right\\}=\alpha\beta,\quad\text{for
}x=0,1;$
$\displaystyle\tau_{\mathrm{ADE}}(x)=\mathbb{E}\left\\{Y(1,\mathcal{G}(x))-Y(0,\mathcal{G}(x))\right\\}=\gamma,\quad\text{for
}x=0,1.$

###  2.2 Method

This section introduces an approach to estimate model parameters in (8) and
(9) and draw inference based on resampling techniques. The parameter set
includes not only the model coefficients, but also the projection vector,
$\boldsymbol{\theta}$. In addition, the covariance matrices,
$\boldsymbol{\Sigma}$’s, are not directly observable. Rather, the observed
data is the realization of the vertices. Let
$\mathbf{M}_{it}=(M_{it1},\dots,M_{itp})^{\top}\in\mathbb{R}^{p}$ denote the
$t$th replicate of the $p$ vertices acquired from unit $i$, for
$t=1,\dots,T_{i}$ and $i=1,\dots,n$, where $T_{i}$ is the number of
realizations in unit $i$ and $n$ is the number of units. Let $X_{i}$ denote
the actual treatment assignment, $Y_{i}$ denote the observed outcome, and
$\mathbf{W}_{i}\in\mathbb{R}^{q}$ denote the $q$-dimensional observed
confounding factors of unit $i$. Assume the model errors, $\eta_{i}$ and
$\epsilon_{i}$, are normally distributed with mean zero and variances
$\pi^{2}$ and $\sigma^{2}$, respectively. Models (8) and (9) can be viewed as
a multilevel model. Thus, it is proposed to estimate the parameters by
maximizing the hierarchical likelihood, which has been shown to be
asymptotically equivalent to maximizing the marginal likelihood (Lee and
Nelder,, 1996). To simplify the notations, let

$\mathbf{X}_{i}=\begin{pmatrix}X_{i},\\\
\mathbf{W}_{i}\end{pmatrix}\in\mathbb{R}^{q+1},~{}\boldsymbol{\alpha}=\begin{pmatrix}\alpha\\\
\boldsymbol{\phi}_{1}\end{pmatrix}\in\mathbb{R}^{q+1},~{}\text{and
}\boldsymbol{\gamma}=\begin{pmatrix}\gamma\\\
\boldsymbol{\phi}_{2}\end{pmatrix}\in\mathbb{R}^{q+1}.$

The following gives the negative hierarchical likelihood (with a constant
difference).

$\displaystyle\ell$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\frac{T_{i}}{2}\left\\{(\alpha_{0i}+\mathbf{X}_{i}^{\top}\boldsymbol{\alpha})+(\boldsymbol{\theta}^{\top}\mathbf{S}_{i}\boldsymbol{\theta})\exp(-\alpha_{0i}-\mathbf{X}_{i}^{\top}\boldsymbol{\alpha})\right\\}$
$\displaystyle+\sum_{i=1}^{n}\frac{1}{2}\left\\{\log\sigma^{2}+\frac{1}{\sigma^{2}}(Y_{i}-\gamma_{0}-\mathbf{X}_{i}^{\top}\boldsymbol{\gamma}-\beta\log(\boldsymbol{\theta}^{\top}\boldsymbol{\Sigma}_{i}\boldsymbol{\theta}))^{2}\right\\}+\sum_{i=1}^{n}\frac{1}{2}\left\\{\log\pi^{2}+\frac{1}{\pi^{2}}(\alpha_{0i}-\alpha_{0})^{2}\right\\},$

where $\alpha_{0i}=\alpha_{0}+\eta_{i}$ and
$\mathbf{S}_{i}=T_{i}^{-1}\sum_{t=1}^{T_{i}}\mathbf{M}_{it}\mathbf{M}_{it}^{\top}$.
The first component in the likelihood corresponds to the conditional
likelihood of $\mathbf{M}_{it}$ given $\alpha_{0i}$ and $\mathbf{X}_{i}$ under
model (8); the second component corresponds to the conditional likelihood of
$Y_{i}$ given $\alpha_{0i}$, $\mathbf{X}_{i}$, and the graph,
$\boldsymbol{\Sigma}_{i}$, under model (9); and the third component
corresponds to the likelihood of the random intercept, $\alpha_{0i}$. For the
second component in (2.2), the covariance matrices,
$\boldsymbol{\Sigma}_{i}$’s, are not directly observed. A sample counterpart
needs to be introduced and utilized instead in practice for estimation. For
fixed $p$, denote $\hat{\boldsymbol{\Sigma}}_{i}$ as a consistent estimator of
$\boldsymbol{\Sigma}_{i}$. For example, the sample covariance matrix,
$\mathbf{S}_{i}$, is a consistent estimator of $\boldsymbol{\Sigma}_{i}$.
Denote $\hat{\ell}$ as the sample counterpart of $\ell$ by replacing
$\boldsymbol{\Sigma}_{i}$ with $\hat{\boldsymbol{\Sigma}}_{i}$.

$\displaystyle\hat{\ell}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\frac{T_{i}}{2}\left\\{(\alpha_{0i}+\mathbf{X}_{i}^{\top}\boldsymbol{\alpha})+(\boldsymbol{\theta}^{\top}\mathbf{S}_{i}\boldsymbol{\theta})\exp(-\alpha_{0i}-\mathbf{X}_{i}^{\top}\boldsymbol{\alpha})\right\\}$
$\displaystyle+\sum_{i=1}^{n}\frac{1}{2}\left\\{\log\sigma^{2}+\frac{1}{\sigma^{2}}(Y_{i}-\gamma_{0}-\mathbf{X}_{i}^{\top}\boldsymbol{\gamma}-\beta\log(\boldsymbol{\theta}^{\top}\hat{\boldsymbol{\Sigma}}_{i}\boldsymbol{\theta}))^{2}\right\\}+\sum_{i=1}^{n}\frac{1}{2}\left\\{\log\pi^{2}+\frac{1}{\pi^{2}}(\alpha_{0i}-\alpha_{0})^{2}\right\\}.$

######  Lemma 1.

For fixed $p$ and $n$, assume $\hat{\boldsymbol{\Sigma}}_{i}$ is a consistent
estimator of $\boldsymbol{\Sigma}_{i}$ as $\min_{i}T_{i}\rightarrow\infty$,
for $\forall~{}i\in\\{1,\dots,n\\}$. Then,

$\hat{\ell}\rightarrow\ell,\quad\text{as }\min_{i}T_{i}\rightarrow\infty.$

With Lemma 1, it is proposed to estimate model parameters by minimizing
$\hat{\ell}$. A coordinate descent algorithm is considered to solve for
solution over the parameter set
$\boldsymbol{\Theta}=(\boldsymbol{\theta},\alpha_{0i},\alpha_{0},\boldsymbol{\alpha},\gamma_{0},\boldsymbol{\gamma},\beta,\pi^{2},\sigma^{2})$.
For $\boldsymbol{\theta}$, a constraint is required to identify a unique
solution. The following optimization problem is considered.

minimize $\displaystyle\hat{\ell},$ such that
$\displaystyle\boldsymbol{\theta}^{\top}\mathbf{H}\boldsymbol{\theta}=1,$ (12)

where $\mathbf{H}\in\mathbb{R}^{p\times p}$ is a positive definite matrix. The
choice of $\mathbf{H}$ can be subject dependent. Similar to the choice in PCA
or CCA, one can set $\mathbf{H}$ to the $p$-dimensional identity matrix. As
$\mathbf{M}_{it}$’s are assumed to be normally distributed, one can set
$\mathbf{H}=\bar{\mathbf{S}}=\sum_{i=1}^{n}\sum_{j=1}^{T_{i}}\mathbf{M}_{it}\mathbf{M}_{it}^{\top}/\sum_{i=1}^{n}T_{i}$
to incorporate distributional information of the data. This type of choice was
also considered in the study of common PCA (Krzanowski,, 1984) and a
covariance regression model by Zhao et al., (2021).

Algorithm 1 summarizes the optimization steps of problem (12). For parameters
without an analytic solution, including $\alpha_{0i}$ and
$\boldsymbol{\alpha}$, the Newton-Raphson algorithm is employed to find the
update. For
$\\{\alpha_{0},\gamma_{0},\boldsymbol{\gamma},\beta,\pi^{2},\sigma^{2}\\}$,
explicit solutions can be obtained and utilized for update. For
$\boldsymbol{\theta}$, with the constraint, the following gives the Lagrangian
form.

$\mathcal{L}(\boldsymbol{\theta},\lambda)=\frac{1}{2}\sum_{i=1}^{n}\left\\{\boldsymbol{\theta}^{\top}(T_{i}U_{i}\mathbf{S}_{i})\boldsymbol{\theta}+\frac{1}{\sigma^{2}}(V_{i}-\beta\log(\boldsymbol{\theta}^{\top}\hat{\boldsymbol{\Sigma}}_{i}\boldsymbol{\theta}))^{2}\right\\}-\lambda(\boldsymbol{\theta}^{\top}\mathbf{H}\boldsymbol{\theta}-1),$
(13)

where $U_{i}=\exp(-\alpha_{0i}-\mathbf{X}_{i}^{\top}\boldsymbol{\alpha})$,
$V_{i}=Y_{i}-\gamma_{0}-\mathbf{X}_{i}^{\top}\boldsymbol{\gamma}$, and
$\lambda$ is the Lagrangian parameter. Take the partial derivative over
$\boldsymbol{\theta}$ and set it to zero, it gives

$\frac{\partial\mathcal{L}(\boldsymbol{\theta},\lambda)}{\partial\boldsymbol{\theta}}=\sum_{i=1}^{n}\left[(T_{i}U_{i}\mathbf{S}_{i})\boldsymbol{\theta}-\frac{1}{\sigma^{2}}\left\\{V_{i}-\beta\log(\boldsymbol{\theta}^{\top}\hat{\boldsymbol{\Sigma}}_{i}\boldsymbol{\theta})\right\\}\frac{2\beta\hat{\boldsymbol{\Sigma}}_{i}\boldsymbol{\theta}}{\boldsymbol{\theta}^{\top}\hat{\boldsymbol{\Sigma}}_{i}\boldsymbol{\theta}}\right]-2\lambda\mathbf{H}\boldsymbol{\theta}=\boldsymbol{\mathrm{0}}.$
(14)

The above equation is hard to derive the explicit form to solve for
$\boldsymbol{\theta}$. It is observed that the quantity of
$\boldsymbol{\theta}^{\top}\hat{\boldsymbol{\Sigma}}_{i}\boldsymbol{\theta}$
is in both a logarithmic function and the denominator. To simplify the
computation, it is proposed to plug in the value of $\boldsymbol{\theta}$ from
the previous step $s$ into
$\boldsymbol{\theta}^{\top}\hat{\boldsymbol{\Sigma}}_{i}\boldsymbol{\theta}$
denoted as
$\xi_{i}^{(s)}=\boldsymbol{\theta}^{(s)\top}\hat{\boldsymbol{\Sigma}}_{i}\boldsymbol{\theta}^{(s)}$.
Let

$\mathbf{A}^{(s)}=\sum_{i=1}^{n}\left\\{T_{i}U_{i}^{(s)}\mathbf{S}_{i}-\frac{2\beta^{(s)}(V_{i}^{(s)}-\beta^{(s)}\log\xi_{i}^{(s)})}{\sigma^{2(s)}\xi_{i}^{(s)}}\hat{\boldsymbol{\Sigma}}_{i}\right\\},$

where $U_{i}^{(s)}$, $V_{i}^{(s)}$, $\beta^{(s)}$, and $\sigma^{2(s)}$ are the
solutions from the $s$th step. It is proposed to solve the following instead
of the equation in (14),

$\mathbf{A}^{(s)}\boldsymbol{\theta}-\lambda\mathbf{H}\boldsymbol{\theta}=\boldsymbol{\mathrm{0}},$
(15)

where the solution $(\boldsymbol{\theta},\lambda)$ is the pair of eigenvector
and eigenvalue of $\mathbf{A}^{(s)}$ with respect to $\mathbf{H}$ that
minimizes $\mathcal{L}(\boldsymbol{\theta},\lambda)$. More details of the
optimization algorithm is presented in Section B of the supplementary
materials. After obtaining an estimate of model parameters, causal estimands
are estimated following Theorem 1.

Algorithm 1 The optimization algorithm for problem (12).

1:$\\{X_{i},\\{\mathbf{M}_{i1},\dots,\mathbf{M}_{iT_{i}}),Y_{i},\mathbf{W}_{i}~{}|~{}i=1,\dots,n\\}$

2:initialization:
$(\boldsymbol{\theta}^{(0)},\alpha_{0i}^{(0)},\alpha_{0}^{(0)},\boldsymbol{\alpha}^{(0)},\gamma_{0}^{(0)},\boldsymbol{\gamma}^{(0)},\beta^{(0)},\pi^{2(0)},\sigma^{2(0)})$

3:repeat for iteration $s=0,1,2,\dots$

4: update $\boldsymbol{\alpha}$ and $\\{\alpha_{0i}\\}$ using the Newton-
Raphson algorithm, denoted as $\boldsymbol{\alpha}^{(s+1)}$ and
$\\{\alpha_{0i}^{(s+1)}\\}$,

5: update
$(\alpha_{0},\gamma_{0},\boldsymbol{\gamma},\beta,\pi^{2},\sigma^{2})$ with
$\alpha_{0}^{(s+1)}=\frac{1}{n}\sum_{i=1}^{n}\alpha_{0i}^{(s+1)},\quad\pi^{2(s+1)}=\frac{1}{n}\sum_{i=1}^{n}\left(\alpha_{0i}^{(s+1)}-\alpha_{0}^{(s+1)}\right)^{2},$
$\boldsymbol{\mu}^{(s+1)}=\left(\sum_{i=1}^{n}\mathbf{Z}_{i}^{(s)}\mathbf{Z}_{i}^{(s)\top}\right)^{-1}\left(\sum_{i=1}^{n}Y_{i}\mathbf{Z}_{i}^{(s)}\right),\quad\sigma^{2(s+1)}=\frac{1}{n}\sum_{i=1}^{n}\left(Y_{i}-\mathbf{Z}_{i}^{(s)\top}\boldsymbol{\mu}^{(s+1)}\right)^{2},$
where  $\mathbf{Z}_{i}^{(s)}=\begin{pmatrix}1\\\ \mathbf{X}_{i}\\\
\log(\boldsymbol{\theta}^{(s)\top}\mathbf{S}_{i}\boldsymbol{\theta}^{(s)})\end{pmatrix}\text{
and }\boldsymbol{\mu}^{(s+1)}=\begin{pmatrix}\gamma_{0}^{(s+1)}\\\
\boldsymbol{\gamma}^{(s+1)}\\\ \beta^{(s+1)}\end{pmatrix},$

6: update $\boldsymbol{\theta}$ by solving (15), denoted as
$\boldsymbol{\theta}^{(s+1)}$,

7:until the objective function in (12) converges;

8:consider a random series of initializations, repeat Steps 1–6, and choose
the estimates with the minimum objective value.

9:$(\hat{\boldsymbol{\theta}},\hat{\alpha}_{0i},\hat{\alpha}_{0},\hat{\boldsymbol{\alpha}},\hat{\gamma}_{0},\hat{\boldsymbol{\gamma}},\hat{\beta},\hat{\pi}^{2},\hat{\sigma}^{2})$

Algorithm 1 offers an approach to identify one mediation component based on
the likelihood criterion. To identify higher-order components, it is proposed
to remove the identified components from the data first and then replace the
input in Algorithm 1 with the new data to identify the next component. Let
$\hat{\boldsymbol{\Theta}}^{(k)}=(\hat{\boldsymbol{\theta}}_{1},\dots,\hat{\boldsymbol{\theta}}_{k})\in\mathbb{R}^{p\times
k}$ denote the first $k$ identified components. For $i=1,\dots,n$, set

$\hat{\mathbf{M}}_{i}^{(k+1)}=\mathbf{M}_{i}-\mathbf{M}_{i}\hat{\boldsymbol{\Theta}}^{(k)}\hat{\boldsymbol{\Theta}}^{(k)\top}\text{
and
}\hat{Y}_{i}^{(k+1)}=Y_{i}-\sum_{j=1}^{k}\hat{\beta}_{j}\log(\hat{\boldsymbol{\theta}}_{j}^{\top}\hat{\boldsymbol{\Sigma}}_{i}\hat{\boldsymbol{\theta}}_{j})$
(16)

as the new data, where
$\mathbf{M}_{i}=(\mathbf{M}_{i1},\dots,\mathbf{M}_{iT_{i}})^{\top}\in\mathbb{R}^{T_{i}\times
p}$ is the mediator outcome of unit $i$ and $\hat{\beta}_{j}$ is the estimate
of $\beta$ in the $j$th component, for $j=1,\dots,k$.
$\hat{\mathbf{M}}_{i}^{(k+1)}$ is created following a similar strategy as in
the principal component analysis to identify orthogonal components, which is
also employed in a recently introduced covariance regression model (Zhao et
al.,, 2021). It is also proposed to remove the identified mediation effect
from the outcome in order to identify a new mediation component. In this
sense, the mediation mechanism of the identified components are considered
parallel. Fitting the mediation components marginally is then equivalent to
fitting them jointly (VanderWeele,, 2015). To determine the number of
mediation components, a criterion called the average deviation from
diagonality introduced in Zhao et al., (2021) is utilized, which is defined as

$\mathrm{DfD}(\hat{\boldsymbol{\Theta}}^{(k)})=\prod_{i=1}^{n}\left[\frac{\det\left\\{\mathrm{diag}(\hat{\boldsymbol{\Theta}}^{(k)\top}\hat{\boldsymbol{\Sigma}}_{i}\hat{\boldsymbol{\Theta}}^{(k)})\right\\}}{\det(\hat{\boldsymbol{\Theta}}^{(k)\top}\hat{\boldsymbol{\Sigma}}_{i}\hat{\boldsymbol{\Theta}}^{(k)})}\right]^{T_{i}/\sum_{i}T_{i}},$
(17)

where $\mathrm{diag}(\mathbf{A})$ is a diagonal matrix taking the diagonal
elements in a square matrix $\mathbf{A}$ and $\det(\mathbf{A})$ is the
determinant of $\mathbf{A}$. The above $\mathrm{DfD}$ metric achieves its
minimum value of one when $\hat{\boldsymbol{\Theta}}^{(k)}$ commonly
diagonalizes $\hat{\boldsymbol{\Sigma}}_{i}$ for all $i\in\\{1,\dots,n\\}$. As
$k$ increases, it gets more difficult to diagonalize all estimated covariance
matrices and the value of the $\mathrm{DfD}$ metric may increase dramatically.
The number of components can be then chosen by setting a threshold. For
example, choose $k$ such that
$\mathrm{DfD}(\hat{\boldsymbol{\Theta}}^{(k)})\leq 2$ as recommended in Zhao
et al., (2021).

###  2.3 Inference

To draw inference on the causal estimands, a resampling procedure is
considered. Particularly, inference on the average indirect effect, which is
denoted as the product of $\alpha$ and $\beta$ according to Theorem 1, raises
concerns. In general, even under the normality assumption, the distribution of
the product of two estimates can be far from Gaussian in finite sample. Thus,
the following nonparametric bootstrap procedure is introduced.

Step 0.

Using all sample units, $k$ mediation components are identified and denote the
estimated linear projections as $\hat{\boldsymbol{\theta}}_{j}$, for
$j=1,\dots,k$.

Step 1.

Generate a bootstrap sample of size $n$ by sampling with replacement, denoted
as
$\\{X_{i}^{*},(\mathbf{M}_{i1},\dots,\mathbf{M}_{iT_{i}})^{*},Y_{i}^{*},\mathbf{W}_{i}^{*}~{}|~{}i=1,\dots,n\\}$.

Step 2.

For $j=1,\dots,k$, using the resampled data, estimate model coefficients and
variances following Algorithm 1 with
$\boldsymbol{\theta}=\hat{\boldsymbol{\theta}}_{j}$.

Step 3.

Repeat Steps 1–2 for $B$ times.

Step 4.

Construct bootstrap confidence intervals for the causal estimands under a
prespecified significance level.

As the study focus is on the causal parameters, the resampling is conducted at
the unit level. Once a unit is sampled, all mediator observations remain for
estimation. The above procedure is implemented fixing the $k$ identified
projections. Thus, inference on $\boldsymbol{\theta}$ is unachievable. If one
seeks to perform inference on $\boldsymbol{\theta}$ via a bootstrap procedure,
it requires a matching step on the estimates across bootstrap samples as the
number of components and the order of identifying the components may differ.
In addition, the inference performance can be sensitive to the quality of the
matching step and the metrics used for matching. On the other hand,
considering the computation cost, it will be significantly inflated when
adding the estimation of $\boldsymbol{\theta}$ and matching. Thus, we leave
the inference on $\boldsymbol{\theta}$ to future research.

###  2.4 Asymptotic properties

In this section, we study the asymptotic properties of the proposed estimator.
We first discuss the regularity conditions on the graph mediator. As discussed
in Section 2.1, under the Gaussian covariance graph model, the graph and the
covariance matrix have one-to-one correspondence. Considering the covariance
matrix, $\boldsymbol{\Sigma}_{i}$ (for $i=1,\dots,n$), it is assumed that it
has the eigendecomposition of
$\boldsymbol{\Sigma}_{i}=\boldsymbol{\Pi}_{i}\boldsymbol{\Lambda}_{i}\boldsymbol{\Pi}_{i}^{\top}$,
where
$\boldsymbol{\Pi}_{i}=(\boldsymbol{\pi}_{i1},\dots,\boldsymbol{\pi}_{ip})\in\mathbb{R}^{p\times
p}$ is an orthonormal matrix and
$\boldsymbol{\Lambda}_{i}=\mathrm{diag}\\{\lambda_{i1},\dots,\lambda_{ip}\\}\in\mathbb{R}^{p\times
p}$ is a diagonal matrix of corresponding eigenvalues. Let
$\boldsymbol{\zeta}_{it}=\boldsymbol{\Pi}_{i}^{\top}\mathbf{M}_{it}\in\mathbb{R}^{p}$,
and then $\mathrm{Cov}(\boldsymbol{\zeta}_{it})=\boldsymbol{\Lambda}_{i}$. The
elements in $\boldsymbol{\zeta}_{it}$ are uncorrelated. Under the normality
assumption, they are mutually independent. In addition, the following
assumptions are imposed.

Assumption B1

Let $T=\min_{i}T_{i}$. $p\ll T$ is fixed.

Assumption B2

The eigenvectors of $\boldsymbol{\Sigma}_{i}$ are identical, that is
$\boldsymbol{\Pi}_{i}=\boldsymbol{\Pi}$, for $i=1,\dots,n$.

Assumption B3

For $\forall~{}i=1,\dots,n$, there exists (at least) one column in
$\boldsymbol{\Pi}_{i}$ indexed by $k_{i}$, such that
$\boldsymbol{\theta}=\boldsymbol{\pi}_{ik_{i}}$ and the parametric model
assumption (Assumption A5) holds.

Assumption B1 assumes a low-dimensional scenario for the Gaussian covariance
graph model. Under this assumption, the sample covariance matrices are well-
conditioned and consistent estimators. Replace $\hat{\boldsymbol{\Sigma}}_{i}$
with the sample covariance matrix, Lemma 1 holds. Assumption B2 is a common-
diagonalization assumption assuming all the covariance matrices have the same
set of eigenvectors. However, the order of the eigenvectors may vary if
following the descending order of the corresponding eigenvalues. Assumption B3
assumes the parametric models in Assumption A5 are correctly specified. Under
these assumptions, one can choose the eigenvectors of $\bar{\mathbf{S}}$ as
the initial value of $\boldsymbol{\theta}$ in Algorithm 1, where
$\bar{\mathbf{S}}=\sum_{i=1}^{n}\sum_{j=1}^{T_{i}}\mathbf{M}_{it}\mathbf{M}_{it}^{\top}/\sum_{i=1}^{n}T_{i}$
is the average sample covariance across all units. The following proposition
demonstrates the consistency of the proposed estimator.

######  Proposition 1.

Assume Assumptions B1–B3 hold. As $n,T\rightarrow\infty$, the proposed
estimator of model parameters is asymptotically consistent.

##  3 Simulation Study

This section compares the performance of the proposed mediation framework with
a competing approach via simulation studies. As no existing approach can be
directly implemented for mediation analysis with a graph mediator, or when the
mediator is a covariance matrix, an approach integrating the covariate
assisted principal (CAP) regression (Zhao et al.,, 2021) and regression-based
mediation analysis, named as CAP-Med, is considered as the competing method.
The CAP-Med approach has two steps.  (i) Perform CAP analysis under Model (8)
to identify the mediator components.  (ii) For each identified component in
Step (i), perform the regression-based mediation analysis as in Imai et al.,
(2010).  The proposed graph mediation approach is named as GMed.

Two simulation studies are conducted:  (1) to examine and compare the
performance of GMed to CAP-Med; and  (2) to examine the robustness of GMed to
model misspecification.  In both studies, common eigenstructure is assumed
across covariance matrices and the covariance matrices are generated from the
eigendecomposition,
$\boldsymbol{\Sigma}_{i}=\boldsymbol{\Pi}\boldsymbol{\Lambda}_{i}\boldsymbol{\Pi}^{\top}$,
where
$\boldsymbol{\Pi}=(\boldsymbol{\pi}_{1},\dots,\boldsymbol{\pi}_{p})\in\mathbb{R}^{p\times
p}$ is an orthonormal matrix and
$\boldsymbol{\Lambda}_{i}=\mathrm{diag}\\{\lambda_{i1},\dots,\lambda_{ip}\\}\in\mathbb{R}^{p\times
p}$ is a diagonal matrix of individual eigenvalues, for $i=1,\dots,n$. In
Simulation (1), two components, the second component (D2) and the forth
component (D4), are chosen to satisfy the model assumption, Assumption A5. In
both components, the magnitude of model parameters are set to one and both the
direct effect ($\tau_{\text{ADE}}$) and indirect effect ($\tau_{\text{AIE}}$)
are with the value of one. The exposure ($X$) is generated from a Bernoulli
distribution with probability $0.5$ of being one. The random errors, $\eta$
and $\epsilon$, are generated from normal distributions with mean zero and
standard deviation $0.1$. For the rest dimensions, the eigenvalues are
generated from a normal distribution with mean value exponentially decaying
from $3$ to $-1$ and standard deviation $0.1$. With the generated covariance
matrices, $\mathbf{M}_{it}$ are generated from the multivariate normal
distribution with mean zero and $Y_{i}$ is generated from Model (9). In this
simulation study, no covariate ($W$) is considered. Two scenarios of data
dimension are considered, $p=10$ and $p=50$, and the sample sizes are set to
be $(n,T)=(500,100)$ and $(n,T)=(500,500)$, where the number of observations
within each unit is set to be identical ($T_{i}=T$). The number of units ($n$)
is set to $500$ to be comparable with the sample size in the application study
in Section 4. In Simulation (2), two covariates ($q=2$),
$\mathbf{W}=(W_{1},W_{2})^{\top}$, one continuous and one binary, are
considered and the corresponding model coefficients with the magnitude of
$0.5$. The continuous covariate ($W_{1}$) is generated from a normal
distribution with mean zero and standard deviation $0.5$ and the binary
covariate ($W_{2}$) is generated from a Bernoulli distribution with
probability $0.5$ of being one. The rest model parameters are the same as in
Simulation (1) and data are generated following (8) and (9). In this
simulation, to evaluate the performance to model misspecification, the
proposed GMed approach is applied ignoring the covariates, named as GMed-Mis.
The result under the correctly specified model is denoted as GMed. To evaluate
the performance of identifying mediation components, the absolute value of the
inner product of the estimated projection and the truth, denoted as
$|\langle\hat{\boldsymbol{\theta}},\boldsymbol{\pi}_{j}\rangle|$ ($j=2,4$), is
used as a metric of similarity between two unit-norm vectors. For both CAP and
GMed approaches, the number of components is determined using the
$\mathrm{DfD}\leq 2$ criterion. Simulations are repeated for $200$ times.

Table 1 presents the performance of estimating mediation components and the
corresponding average indirect effect. In Simulation (1), considering Model
(8) only, the CAP approach shall achieve the optimal performance in estimating
$\boldsymbol{\theta}$ and model coefficients as the method is likelihood-based
and assumptions required are satisfied. The proposed approach yields
comparable performance to the CAP-Med approach. As the number of observations
within each unit increases, the performance of GMed in estimating
$\boldsymbol{\theta}$ improves more significantly compared to the CAP-Med
approach. In addition, the GMed approach performs consistently better in
estimating $\boldsymbol{\theta}$ in the second mediation component (D4). This
suggests that incorporating the outcome information in the estimating
procedure helps improve performance. To evaluate the finite sample performance
of the GMed approach, Figure 2 presents the performance with $p=10$ and
various sample size combinations. From the figure, as both the number of units
and the number of observations within each unit increase, the similarity of
$\hat{\boldsymbol{\theta}}$ to the truth converges to one and the bias and
mean squared error (MSE) of the estimated average indirect effect converge to
zero.

Simulation (2) aims to evaluate the robustness of GMed to model
misspecification. Under the current setting, it also examines the robustness
of GMed to the existence of an additive unmeasured mediator-outcome
confounding not induced by the exposure. From Table 1, when the models are
misspecified, the projections are still correctly identified though the
estimate of the average indirect effect is off. This demonstrates the
robustness of the GMed approach in estimating relevant mediation components.
Utilizing this property, Section C.1 of the supplementary materials suggests
an approach for sensitivity analysis.

Table 1: Performance of estimating target mediation components ($\boldsymbol{\theta}$) and corresponding average indirect effect ($\tau_{\text{AIE}}$) over $200$ replicates in the simulation study. SE: standard error; MSE: mean squared error. |  |  |  |  | $\hat{\boldsymbol{\theta}}$ |  | $\hat{\tau}_{\mathrm{AIE}}$
---|---|---|---|---|---|---|---
Case | $p$ | $(n,T)$ | Method |  |  $|\langle\hat{\boldsymbol{\theta}},\boldsymbol{\pi}_{j}\rangle|$ (SE) |  | Bias | MSE
|  |  |  | D2 |  $0.879$ ($0.113$) |  | $-0.677$ | $0.460$
|  |  | GMed | D4 |  $0.929$ ($0.140$) |  | $-0.653$ | $0.428$
|  |  |  | D2 |  $0.886$ ($0.138$) |  | $-0.664$ | $0.443$
|  | $(500,100)$ | CAP-Med | D4 |  $0.905$ ($0.120$) |  | $-0.635$ | $0.406$
|  |  |  | D2 |  $0.962$ ($0.047$) |  | $-0.292$ | $0.088$
|  |  | GMed | D4 |  $0.977$ ($0.096$) |  | $-0.276$ | $0.080$
|  |  |  | D2 |  $0.911$ ($0.134$) |  | $-0.243$ | $0.067$
| $10$ | $(500,500)$ | CAP-Med | D4 |  $0.919$ ($0.121$) |  | $-0.129$ | $0.038$
|  |  |  | D2 |  $0.841$ ($0.101$) |  | $-0.665$ | $0.444$
|  |  | GMed | D4 |  $0.910$ ($0.142$) |  | $-0.659$ | $0.436$
|  |  |  | D2 |  $0.841$ ($0.104$) |  | $-0.660$ | $0.438$
|  | $(500,100)$ | CAP-Med | D4 |  $0.898$ ($0.115$) |  | $-0.632$ | $0.402$
|  |  |  | D2 |  $0.915$ ($0.095$) |  | $-0.289$ | $0.087$
|  |  | GMed | D4 |  $0.985$ ($0.073$) |  | $-0.277$ | $0.080$
|  |  |  | D2 |  $0.905$ ($0.109$) |  | $-0.253$ | $0.071$
(1) | $50$ | $(500,500)$ | CAP-Med | D4 |  $0.904$ ($0.137$) |  | $-0.129$ | $0.041$
|  |  |  | D2 |  $0.993$ ($0.005$) |  | $-0.669$ | $0.449$
|  |  | GMed | D4 |  $0.999$ ($0.001$) |  | $-0.674$ | $0.456$
|  |  |  | D2 |  $0.993$ ($0.005$) |  | $-1.319$ | $1.863$
|  | $(500,100)$ | GMed-Mis | D4 |  $0.999$ ($0.001$) |  | $-0.114$ | $0.014$
|  |  |  | D2 |  $0.999$ ($0.001$) |  | $-0.291$ | $0.087$
|  |  | GMed | D4 |  $1.000$ ($0.000$) |  | $-0.279$ | $0.081$
|  |  |  | D2 |  $0.999$ ($0.001$) |  | $-1.370$ | $2.076$
(2) | $10$ | $(500,500)$ | GMed-Mis | D4 |  $1.000$ ($0.000$) |  | $-0.002$ | $0.001$

((a)) D2: $|\langle\hat{\boldsymbol{\theta}},\boldsymbol{\pi}_{j}\rangle|$

((b)) D2: Bias of $\hat{\tau}_{\mathrm{AIE}}$

((c)) D2: MSE of $\hat{\tau}_{\mathrm{AIE}}$

((d)) D4: $|\langle\hat{\boldsymbol{\theta}},\boldsymbol{\pi}_{j}\rangle|$

((e)) D4: Bias of $\hat{\tau}_{\mathrm{AIE}}$

((f)) D4: MSE of $\hat{\tau}_{\mathrm{AIE}}$

Figure 2: Performance of estimating target mediation components
($\boldsymbol{\theta}$) and corresponding average indirect effect
($\tau_{\text{AIE}}$) over $200$ replicates in the simulation study as the
sample size $(n,T)$ varies with $p=10$. MSE: mean squared error.

##  4 Application

We apply the proposed approach to data acquired from the National Consortium
on Alcohol and NeuroDevelopment in Adolescence (NCANDA). The consortium aims
to identify the deterministic effect of alcohol use on the developing
adolescent brain. A core battery of measurements, including structural and
functional brain scans and cognitive testing, has been developed. In
adolescence, the sex difference in cognitive behaviors has been identified in
different functional domains, including spatial, verbal, math abilities,
social recognition and so on (Lauer et al.,, 2019; Meinhardt-Injac et al.,,
2020; Esnaola et al.,, 2020). This difference has been shown to be partially
mediated by brain functional connectivity captured by resting-state fMRI
(Alarcón et al.,, 2018). In this study, we aim to identify the brain
subnetwork within which functional connectivity mediates the sex difference in
cognition. To exclude the impact of alcohol consumption, $n=621$ subjects
($312$ Males and $309$ Females) aged between $12$ and $22$ without excessive
drinking are analyzed. Among these subjects, a significantly lower median
response time in the motor speed test is observed in males compared to females
after adjusting for age ($\mathrm{ATE}=-0.324$, $p\text{-value}<0.001$), where
the motor task measures sensorimotor ability via having the participant use
the mouse to click on a shrinking box when it moves to a new position on the
screen. We further apply the proposed approach for mediation analysis, where
sex is the binary exposure ($X$, $X=1$ for male), resting-state fMRI data are
the mediator ($\mathbf{M}$), the $z$-score of median response time for correct
responses of the motor speed task is the outcome ($Y$), and age is the
confounding factor ($W$). After preprocessing, fMRI time courses are extracted
from $p=75$ brain regions ($60$ cortical and $15$ subcortical regions)
spanning the whole brain using the Harvard-Oxford Atlas in FSL (Smith et al.,,
2004). These regions are grouped into $10$ functional modules, which will be
used for an ad hoc procedure of sparsifying the loading profile using the
fused lasso (Tibshirani et al.,, 2005) to impose local smoothness and
constancy within each module. To remove the temporal dependence in the time
courses, a subsample is taken with an effective sample size of $T_{i}=T=125$
and denote the subsampled data as $\mathbf{M}_{i}\in\mathbb{R}^{T_{i}\times
p}$, for $i=1,\dots,n$.

Using the $\mathrm{DfD}\leq 2$ criterion in (17), the proposed approach
identifies four components. Table 2 presents the estimated AIE, as well as
$\alpha$ and $\beta$ coefficient. The confidence intervals are obtained from
$500$ bootstrap samples. Among these identified components, the third
component ($M_{3}$) shows a significantly positive AIE with both $\alpha$ and
$\beta$ negative. Figure 3 presents the sparsified loading profile of
$\boldsymbol{\theta}_{3}$ and the corresponding brain map. Section C.2 of the
supplementary materials presents the results from the CAP-Med approach
introduced in Section 3. The identified component with a significant mediation
effect is consistent with $M_{3}$ identified by the proposed approach. In
$M_{3}$, four regions with a non-zero loading are all in the limbic-system
network, including the temporal pole (left and right) and the medial
orbitofrontal cortex (left and right). Compared to females, (weighted)
functional connectivity within this network is lower in males, while this
lower functional connectivity increases the reaction time. The temporal pole
has been found associated with high-level cognitive functions (Herlin et al.,,
2021). Though no direct relation to reaction time has been established, an
indirect influence was hypothesized by contributing to processes involving
decision-making, response selection, and emotion evaluation (Pessoa,, 2010).
One of the primary functions of the medial orbitofrontal cortex is to
integrate emotional reaction with sensory and/or contextual stimuli playing a
role in reward processing and value-based decision making, which allows
individuals to make adaptive responses to stimuli based on emotional
significance (Rudebeck and Murray,, 2014). Thus, indirectly, activation in the
area may prolong the reaction time. Regional sex differences of the temporal
and frontal cortices have been observed in the developing brain using multiple
imaging modalities. Sex difference in the brain was also suggested to be
relevant to the symptomatic sex difference in psychiatric disorders
(Kaczkurkin et al.,, 2019). Via a mediation analysis, the proposed approach
offers a way of articulating the underlying mechanism.

Table 2: Estimated average indirect effect (AIE) and $\alpha$ and $\beta$
coefficient of the identified mediator components in the NCANDA dataset.
Confidence intervals are constructed from $500$ bootstrap samples. Est.:
estimate; SE: standard error; CI: confidence interval.

| AIE |  | $\alpha$ |  | $\beta$
---|---|---|---|---|---
| Est. (SE) |  $p$-value |  $95\%$ CI |  | Est. (SE) |  $p$-value |  $95\%$ CI |  | Est. (SE) |  $p$-value |  $95\%$ CI
$M_{1}$ |  $-0.018$ ($0.031$) | $0.574$ | $(-0.079,0.044)$ |  |  $-0.325$ ($0.069$) | $<0.001$ | $(-0.460,-0.190)$ |  |  $0.053$ ($0.094$) | $0.572$ | $(-0.131,0.237)$
$M_{2}$ |  $-0.014$ ($0.032$) | $0.676$ | $(-0.077,0.050)$ |  |  $-0.286$ ($0.065$) | $<0.001$ | $(-0.414,-0.157)$ |  |  $0.051$ ($0.110$) | $0.646$ | $(-0.166,0.267)$
$M_{3}$ |  $0.066$ ($0.027$) | $0.014$ | $(0.013,0.119)$ |  |  $-0.211$ ($0.064$) | $<0.001$ | $(-0.336,-0.086)$ |  |  $-0.318$ ($0.094$) | $0.001$ | $(-0.503,-0.134)$
$M_{4}$ |  $-0.035$ ($0.033$) | $0.287$ | $(-0.100,0.030)$ |  |  $0.256$ ($0.042$) | $<0.001$ | $(0.173,0.338)$ |  |  $-0.138$ ($0.127$) | $0.279$ | $(-0.388,0.112)$

((a)) Sparsified loading profile of $\boldsymbol{\theta}_{3}$

((b)) Brain map of $\boldsymbol{\theta}_{3}$

Figure 3: The sparsified loading profile and brain map of the component with a
significant AIE ($M_{3}$).

##  5 Discussion

This study introduces a mediation analysis framework when the mediator is a
graph. A Gaussian covariance graph model is assumed for graph representation.
Causal estimands and assumptions are discussed under this representation. With
a covariance matrix as the mediator, parametric mediation models are
considered based on matrix decomposition. Assuming Gaussian random errors,
likelihood-based estimators are introduced to simultaneously identify the
decomposition and causal parameters. An efficient computational algorithm is
proposed and the asymptotic properties of the estimators are investigated. Via
simulation studies, the performance of the proposed approach is evaluated.
Applying to a resting-state fMRI study, a brain network is identified within
which functional connectivity mediates the sex difference in the performance
of a motor task.

In causal mediation analysis, an essential while untestable assumption is the
assumption of no unmeasured mediator-outcome confounding. A sensitivity
analysis is usually conducted to justify the validity of the conclusion to
this assumption. For parametric approaches, one type of commonly used approach
is to parametrize the confounding effect to evaluate the causal effects under
various values, such as the one proposed in Imai et al., (2010). Using
simulation studies, it is demonstrated that the proposed approach is robust to
the existence of unmeasured mediator-outcome confounding in identifying the
mediation component, that is in estimating the projection vector
$\boldsymbol{\theta}$. With given $\boldsymbol{\theta}$, one can employ the
approach in Imai et al., (2010) for sensitivity analysis.

The asymptotic consistency of the proposed estimator requires the common
diagonalization assumption on the covariance matrices. Via simulation studies,
Zhao et al., (2021) pointed out that this assumption can be relaxed to partial
common diagonalization. As this study also introduces a likelihood-based
procedure, it is expected that the proposed approach is robust to this
relaxation to partial common diagonalization.

Considering a graph mediator, under the Gaussian covariance graph model, this
study assumes the number of nodes in the graph is fixed and low dimensional.
The sample covariance matrices are thus well-conditioned and a likelihood-
based approach is introduced to estimate model parameters. In many practical
settings, for example, in voxel-level fMRI analysis, data dimension can be
even higher than the number of fMRI data points. A well-conditioned estimator
of the covariance matrix is required and we leave the introduction of such an
estimator and the study of theoretical results as one of future directions. As
discussed in Section 2.3, inference on projection vectors is not
straightforward and requires rigorous theoretical and numerical
investigations, which we leave to future research.

.

## Appendix

This Appendix collects the technical proof of the theorems in the main text,
additional technical details, and additional data analysis results.

##  Appendix A Theory and Proof

###  A.1 Proof of Theorem 1

###### Proof.

As discussed in Section 2.1, the potential outcome of $Y$ under a multiple-
worlds model is expressed as

$Y(x,\mathcal{G}(x^{\prime}))=\gamma x+\alpha\beta
x^{\prime}+(\gamma_{0}+\alpha_{0}\beta+\beta\mathbf{W}^{\top}\boldsymbol{\phi}_{1}+\mathbf{W}^{\top}\boldsymbol{\phi}_{2})+\beta\eta+\epsilon.$
$\displaystyle\tau_{\mathrm{ATE}}$ $\displaystyle=$
$\displaystyle\mathbb{E}\left\\{Y(1,\mathcal{G}(1))-Y(0,\mathcal{G}(0))\right\\}$
$\displaystyle=$
$\displaystyle\left\\{\gamma+\alpha\beta+(\gamma_{0}+\alpha_{0}\beta+\beta\mathbf{W}^{\top}\boldsymbol{\phi}_{1}+\mathbf{W}^{\top}\boldsymbol{\phi}_{2})\right\\}-\left\\{(\gamma_{0}+\alpha_{0}\beta+\beta\mathbf{W}^{\top}\boldsymbol{\phi}_{1}+\mathbf{W}^{\top}\boldsymbol{\phi}_{2})\right\\}$
$\displaystyle=$ $\displaystyle\gamma+\alpha\beta.$
$\displaystyle\tau_{\mathrm{AIE}}(x)$ $\displaystyle=$
$\displaystyle\mathbb{E}\left\\{Y(x,\mathcal{G}(1))-Y(x,\mathcal{G}(0))\right\\}$
$\displaystyle=$ $\displaystyle\left\\{\gamma
x+\alpha\beta+(\gamma_{0}+\alpha_{0}\beta+\beta\mathbf{W}^{\top}\boldsymbol{\phi}_{1}+\mathbf{W}^{\top}\boldsymbol{\phi}_{2})\right\\}-\left\\{\gamma
x+(\gamma_{0}+\alpha_{0}\beta+\beta\mathbf{W}^{\top}\boldsymbol{\phi}_{1}+\mathbf{W}^{\top}\boldsymbol{\phi}_{2})\right\\}$
$\displaystyle=$ $\displaystyle\alpha\beta.$
$\displaystyle\tau_{\mathrm{ADE}}(x)$ $\displaystyle=$
$\displaystyle\mathbb{E}\left\\{Y(1,\mathcal{G}(x))-Y(0,\mathcal{G}(x))\right\\}$
$\displaystyle=$ $\displaystyle\left\\{\gamma+\alpha\beta
x+(\gamma_{0}+\alpha_{0}\beta+\beta\mathbf{W}^{\top}\boldsymbol{\phi}_{1}+\mathbf{W}^{\top}\boldsymbol{\phi}_{2})\right\\}-\left\\{\alpha\beta
x+(\gamma_{0}+\alpha_{0}\beta+\beta\mathbf{W}^{\top}\boldsymbol{\phi}_{1}+\mathbf{W}^{\top}\boldsymbol{\phi}_{2})\right\\}$
$\displaystyle=$ $\displaystyle\gamma.$

The above proves Theorem 1. ∎

###  A.2 Proof of Lemma 1

From the lemma assumption, $\hat{\boldsymbol{\Sigma}}_{i}$ is a consistent
estimator of $\boldsymbol{\Sigma}_{i}$. With given $\boldsymbol{\theta}$, the
likelihood function is a continuous function of
$\hat{\boldsymbol{\Sigma}}_{i}$ (for $i=1,\dots,n$). Thus, $\hat{\ell}$
converges to $\ell$ as $\min_{i}T_{i}\rightarrow\infty$.

###  A.3 Proof of Proposition 1

Assumption B1 assumes a low-dimensional scenario and the dimension of the
mediator is fixed. Lemma 1 demonstrates the consistency of the approximated
likelihood function used for optimization. Asymptotically, the proposed
estimator is likelihood-based. Thus, under the imposed assumptions, the
consistency of the estimator follows.

##  Appendix B Details of Algorithm 1

$\displaystyle\hat{\ell}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\frac{T_{i}}{2}\left\\{(\alpha_{0i}+\mathbf{X}_{i}^{\top}\boldsymbol{\alpha})+(\boldsymbol{\theta}^{\top}\mathbf{S}_{i}\boldsymbol{\theta})\exp(-\alpha_{0i}-\mathbf{X}_{i}^{\top}\boldsymbol{\alpha})\right\\}$
$\displaystyle+\sum_{i=1}^{n}\frac{1}{2}\left\\{\log\sigma^{2}+\frac{1}{\sigma^{2}}(Y_{i}-\gamma_{0}-\mathbf{X}_{i}^{\top}\boldsymbol{\gamma}-\beta\log(\boldsymbol{\theta}^{\top}\hat{\boldsymbol{\Sigma}}_{i}\boldsymbol{\theta}))^{2}\right\\}+\sum_{i=1}^{n}\frac{1}{2}\left\\{\log\pi^{2}+\frac{1}{\pi^{2}}(\alpha_{0i}-\alpha_{0})^{2}\right\\}.$

First considering $\boldsymbol{\alpha}$ and $\alpha_{0i}$ ($i=1,\dots,n$),
analytical solution is not available, thus the Newton-Raphson algorithm is
employed. For $\boldsymbol{\alpha}$,

$\frac{\partial\hat{\ell}}{\partial\boldsymbol{\alpha}}=\frac{1}{2}\sum_{i=1}^{n}T_{i}\left\\{1-(\boldsymbol{\theta}^{\top}\mathbf{S}_{i}\boldsymbol{\theta})\exp(-\alpha_{0i}-\mathbf{X}_{i}^{\top}\boldsymbol{\alpha})\right\\}\mathbf{X}_{i},$
$\frac{\partial^{2}\hat{\ell}}{\partial\boldsymbol{\alpha}\partial\boldsymbol{\alpha}^{\top}}=\frac{1}{2}\sum_{i=1}^{n}T_{i}(\boldsymbol{\theta}^{\top}\mathbf{S}_{i}\boldsymbol{\theta})\exp(-\alpha_{0i}-\mathbf{X}_{i}^{\top}\boldsymbol{\alpha})\mathbf{X}_{i}\mathbf{X}_{i}^{\top},$
$\boldsymbol{\alpha}^{(s+1)}=\boldsymbol{\alpha}^{(s)}-\left(\left.\frac{\partial^{2}\hat{\ell}}{\partial\boldsymbol{\alpha}\partial\boldsymbol{\alpha}^{\top}}\right|_{\boldsymbol{\alpha}^{(s)}}\right)^{-1}\left(\left.\frac{\partial\hat{\ell}}{\partial\boldsymbol{\alpha}}\right|_{\boldsymbol{\alpha}^{(s)}}\right).$

For $\alpha_{0i}$ ($i=1,\dots,n$),

$\frac{\partial\hat{\ell}}{\partial\alpha_{0i}}=\frac{1}{2}\left[T_{i}\left\\{1-(\boldsymbol{\theta}^{\top}\mathbf{S}_{i}\boldsymbol{\theta})\exp(-\alpha_{0i}-\mathbf{X}_{i}^{\top}\boldsymbol{\alpha})\right\\}+\frac{2}{\tau^{2}}(\alpha_{0i}-\alpha_{0})\right],$
$\frac{\partial^{2}\hat{\ell}}{\partial\alpha_{0i}^{2}}=\frac{1}{2}\left\\{T_{i}(\boldsymbol{\theta}^{\top}\mathbf{S}_{i}\boldsymbol{\theta})\exp(-\alpha_{0i}-\mathbf{X}_{i}^{\top}\boldsymbol{\alpha})+\frac{2}{\tau^{2}}\right\\},$
$\alpha_{0i}^{(s+1)}=\alpha_{0i}^{(s)}-\left(\left.\frac{\partial^{2}\hat{\ell}}{\partial\alpha_{0i}^{2}}\right|_{\alpha_{0i}^{(s)}}\right)^{-1}\left(\left.\frac{\partial\hat{\ell}}{\partial\alpha_{0i}}\right|_{\alpha_{0i}^{(s)}}\right)$

For $\alpha_{0}$ and $\tau$, explicit form for update is available,

$\alpha_{0}^{(s+1)}=\frac{1}{n}\sum_{i=1}^{n}\alpha_{0i}^{(s+1)},\quad\tau^{2(s+1)}=\frac{1}{n}\sum_{i=1}^{n}\left(\alpha_{0i}^{(s+1)}-\alpha_{0}^{(s+1)}\right)^{2}.$

For $(\gamma_{0},\boldsymbol{\gamma},\beta)$, the analytical solution for
update can be found jointly. Let

$\mathbf{Z}_{i}=\begin{pmatrix}1\\\ \mathbf{X}_{i}\\\
\log(\boldsymbol{\theta}^{\top}\hat{\boldsymbol{\Sigma}}_{i}\boldsymbol{\theta})\end{pmatrix},\quad\boldsymbol{\mu}=\begin{pmatrix}\gamma_{0}\\\
\boldsymbol{\gamma}\\\ \beta\end{pmatrix},$
$\frac{\partial\hat{\ell}}{\partial\boldsymbol{\mu}}=\frac{1}{2}\sum_{i=1}^{n}\frac{2}{\sigma^{2}}(Y_{i}-\mathbf{Z}_{i}^{\top}\boldsymbol{\mu})(-\mathbf{Z}_{i})=\boldsymbol{\mathrm{0}},$
$\Rightarrow\quad\boldsymbol{\mu}^{(s+1)}=\left(\sum_{i=1}^{n}\mathbf{Z}_{i}^{(s)}\mathbf{Z}_{i}^{(s)\top}\right)^{-1}\left(\sum_{i=1}^{n}Y_{i}\mathbf{Z}_{i}^{(s)}\right)=(\mathbf{Z}^{(s)\top}\mathbf{Z}^{(s)})^{-1}\mathbf{Z}^{(s)\top}\mathbf{Y},$

where $\mathbf{Z}=(\mathbf{Z}_{1},\dots,\mathbf{Z}_{n})^{\top}$. For $\sigma$,

$\sigma^{2(s+1)}=\frac{1}{n}\sum_{i=1}^{n}(Y_{i}-\mathbf{Z}_{i}^{(s)\top}\boldsymbol{\mu}^{(s+1)})^{2}.$

For $\boldsymbol{\theta}$, to find the solution to (15), it is equivalent to
find the eigenvectors and eigenvalues of $\mathbf{A}$ with respect to
$\mathbf{H}$. We first assume $\boldsymbol{\theta}_{0}$ is a solution
eigenvector with unit norm, $\|\boldsymbol{\theta}_{0}\|_{2}=1$. Since
$\mathbf{H}$ is positive definite, let
$\boldsymbol{\theta}=\mathbf{H}^{-1/2}\boldsymbol{\theta}_{0}$, then

$\boldsymbol{\theta}^{\top}\mathbf{H}\boldsymbol{\theta}=\boldsymbol{\theta}_{0}^{\top}\mathbf{H}^{-1/2}\mathbf{H}\mathbf{H}^{-1/2}\boldsymbol{\theta}_{0}=\boldsymbol{\theta}_{0}^{\top}\boldsymbol{\theta}_{0}=1,$

which satisfies the constraint condition. Replace $\boldsymbol{\theta}$ with
$\boldsymbol{\theta}=\mathbf{H}^{-1/2}\boldsymbol{\theta}_{0}$ in (15),

$\mathbf{A}\mathbf{H}^{-1/2}\boldsymbol{\theta}_{0}-\lambda\mathbf{H}\mathbf{H}^{-1/2}\boldsymbol{\theta}_{0}=\boldsymbol{\mathrm{0}}.$
$\Rightarrow\quad\mathbf{H}^{-1/2}\mathbf{A}\mathbf{H}^{-1/2}\boldsymbol{\theta}_{0}=\lambda\boldsymbol{\theta}_{0}.$

Therefore, $\boldsymbol{\theta}_{0}$ is the eigenvector of matrix
$\mathbf{H}^{-1/2}\mathbf{A}\mathbf{H}^{-1/2}$ and $\lambda$ is the
corresponding eigenvalue. Then,
$\boldsymbol{\theta}=\mathbf{H}^{-1/2}\boldsymbol{\theta}_{0}$ is the update
for $\boldsymbol{\theta}$.

##  Appendix C Additional Results of the NCANDA Study

###  C.1 Sensitivity Analysis

As discussed in Section 3, the proposed approach is robust to the existence of
an additive unmeasured mediator-outcome confounding in identifying the
projection vectors. Based on this property, a sensitivity analysis introduced
in Imai et al., (2010) is implemented where for an additive unmeasured
mediator-outcome confounding, the correlation between the model errors is
considered as the sensitivity parameter. In the NCANDA application study,
$M_{3}$ shows a significant mediation effect. With the estimated
$\boldsymbol{\theta}$ of $M_{3}$, Figure C.1 presents the sensitivity analysis
plot over the range of the sensitivity parameter ($\rho$). In the figure, the
$95\%$ confidence intervals are constructed from $500$ bootstrap samples. From
the figure, the $95\%$ confidence interval covers zero for the $\rho$ value
between $-0.20$ and $-0.05$. When $\rho>-0.05$, the average indirect effect is
positive and significant; when $\rho<-0.20$, the average indirect effect is
negative and significant.

Figure C.1: Sensitivity analysis of $M_{3}$ in the NCANDA study.

###  C.2 Results from the CAP-Med approach

We also apply the CAP-Med approach introduced in Section 3 to the NCANDA data.
The CAP step identifies $8$ components using the criterion of
$\mathrm{DfD}\leq 2$. After running the mediation analysis, the seventh
component (C7) shows a significantly positive average indirect effect with
$\mathrm{AIE}=0.038$ and $95\%$ CI $(0.004,0.078)$. Figure C.2 presents the
sparsified loading profile of this component and the corresponding regions in
the brain map. Compared to the profile of $M_{3}$ identified by the proposed
approach, GMed, a high similarity between the two is observed with a
similarity of
$|\langle\hat{\boldsymbol{\theta}}_{3}^{\text{(GMed)}},\hat{\boldsymbol{\theta}}_{7}^{\text{(CAP-
Med)}}\rangle|=0.864$. For the rest components identified by GMed, a
corresponding component from CAP-Med with high similarity can be found:
$|\langle\hat{\boldsymbol{\theta}}_{1}^{\text{(GMed)}},\hat{\boldsymbol{\theta}}_{8}^{\text{(CAP-
Med)}}\rangle|=0.805$,
$|\langle\hat{\boldsymbol{\theta}}_{2}^{\text{(GMed)}},\hat{\boldsymbol{\theta}}_{6}^{\text{(CAP-
Med)}}\rangle|=0.772$, and
$|\langle\hat{\boldsymbol{\theta}}_{4}^{\text{(GMed)}},\hat{\boldsymbol{\theta}}_{2}^{\text{(CAP-
Med)}}\rangle|=0.862$. Though both approaches identify these consistent
components, the proposed approach offers an integrated way of targeting the
ones demonstrating a mediation effect rather than performing a two-step
approach.

((a)) Sparsified loading profile of C7

((b)) Brain map of C7

Figure C.2: The sparsified loading profile and brain map of the component with
a significant AIE (C7) using the CAP-Med approach.

## References

  * Alarcón et al., (2018) Alarcón, G., Pfeifer, J. H., Fair, D. A., and Nagel, B. J. (2018).  Adolescent gender differences in cognitive control performance and functional connectivity between default mode and fronto-parietal networks within a self-referential context.  Frontiers in Behavioral Neuroscience, 12:73.
  * Anderson, (1973) Anderson, T. (1973).  Asymptotically efficient estimation of covariance matrices with linear structure.  The Annals of Statistics, 1(1):135–141.
  * Andrews and Didelez, (2020) Andrews, R. M. and Didelez, V. (2020).  Insights into the cross-world independence assumption of causal mediation analysis.  Epidemiology, 32(2):209–219.
  * Baron and Kenny, (1986) Baron, R. M. and Kenny, D. A. (1986).  The moderator–mediator variable distinction in social psychological research: Conceptual, strategic, and statistical considerations.  Journal of Personality and Social Psychology, 51(6):1173.
  * Boik, (2002) Boik, R. J. (2002).  Spectral models for covariance matrices.  Biometrika, 89(1):159–182.
  * Chaudhuri et al., (2007) Chaudhuri, S., Drton, M., and Richardson, T. S. (2007).  Estimation of a covariance matrix with zeros.  Biometrika, 94(1):199–216.
  * Chen and Zhou, (2023) Chen, M. and Zhou, Y. (2023).  Causal mediation analysis with a three-dimensional image mediator.  arXiv preprint arXiv:2303.06560.
  * Chén et al., (2017) Chén, O. Y., Crainiceanu, C., Ogburn, E. L., Caffo, B. S., Wager, T. D., and Lindquist, M. A. (2017).  High-dimensional multivariate mediation with application to neuroimaging data.  Biostatistics, 19(2):121–136.
  * Chiu et al., (1996) Chiu, T. Y., Leonard, T., and Tsui, K.-W. (1996).  The matrix-logarithmic covariance model.  Journal of the American Statistical Association, 91(433):198–210.
  * Derkach et al., (2019) Derkach, A., Pfeiffer, R. M., Chen, T.-H., and Sampson, J. N. (2019).  High dimensional mediation analysis with latent variables.  Biometrics, 75(3):745–756.
  * Edwards, (2012) Edwards, D. (2012).  Introduction to graphical modelling.  Springer Science & Business Media.
  * Esnaola et al., (2020) Esnaola, I., Sesé, A., Antonio-Agirre, I., and Azpiazu, L. (2020).  The development of multiple self-concept dimensions during adolescence.  Journal of Research on Adolescence, 30:100–114.
  * Flury, (1984) Flury, B. N. (1984).  Common principal components in $k$ groups.  Journal of the American Statistical Association, 79(388):892–898.
  * Franks and Hoff, (2019) Franks, A. M. and Hoff, P. (2019).  Shared subspace models for multi-group covariance estimation.  Journal of Machine Learning Research, 20(171):1–37.
  * Friston, (2011) Friston, K. J. (2011).  Functional and effective connectivity: a review.  Brain Connectivity, 1(1):13–36.
  * Gu et al., (2014) Gu, F., Preacher, K. J., and Ferrer, E. (2014).  A state space modeling approach to mediation analysis.  Journal of Educational and Behavioral Statistics, 39(2):117–143.
  * Herlin et al., (2021) Herlin, B., Navarro, V., and Dupont, S. (2021).  The temporal pole: From anatomy to function—a literature appraisal.  Journal of Chemical Neuroanatomy, 113:101925.
  * Hoff, (2009) Hoff, P. D. (2009).  A hierarchical eigenmodel for pooled covariance estimation.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71(5):971–992.
  * Hoff and Niu, (2012) Hoff, P. D. and Niu, X. (2012).  A covariance regression model.  Statistica Sinica, 22(2):729–753.
  * Huang and Pan, (2016) Huang, Y.-T. and Pan, W.-C. (2016).  Hypothesis test of mediation effect in causal mediation model with high-dimensional continuous mediators.  Biometrics, 72(2):402–413.
  * Imai et al., (2010) Imai, K., Keele, L., and Yamamoto, T. (2010).  Identification, inference and sensitivity analysis for causal mediation effects.  Statistical Science, 25(1):51–71.
  * Jiang and Colditz, (2023) Jiang, S. and Colditz, G. A. (2023).  Causal mediation analysis using high-dimensional image mediator bounded in irregular domain with an application to breast cancer.  Biometrics.
  * Kaczkurkin et al., (2019) Kaczkurkin, A. N., Raznahan, A., and Satterthwaite, T. D. (2019).  Sex differences in the developing brain: insights from multimodal neuroimaging.  Neuropsychopharmacology, 44(1):71–85.
  * Krzanowski, (1984) Krzanowski, W. (1984).  Principal component analysis in the presence of group structure.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 33(2):164–168.
  * Lauer et al., (2019) Lauer, J. E., Yhang, E., and Lourenco, S. F. (2019).  The development of gender differences in spatial reasoning: A meta-analytic review.  Psychological Bulletin, 145(6):537.
  * Lee and Nelder, (1996) Lee, Y. and Nelder, J. A. (1996).  Hierarchical generalized linear models.  Journal of the Royal Statistical Society. Series B (Methodological), pages 619–678.
  * Lindquist, (2008) Lindquist, M. A. (2008).  The statistical analysis of fMRI data.  Statistical Science, 23(4):439–464.
  * Lindquist, (2012) Lindquist, M. A. (2012).  Functional causal mediation analysis with an application to brain connectivity.  Journal of the American Statistical Association, 107(500):1297–1309.
  * Meinhardt-Injac et al., (2020) Meinhardt-Injac, B., Daum, M. M., and Meinhardt, G. (2020).  Theory of mind development from adolescence to adulthood: Testing the two-component model.  British Journal of Developmental Psychology, 38(2):289–303.
  * Pearl, (2001) Pearl, J. (2001).  Direct and indirect effects.  In Proceedings of the seventeenth conference on uncertainty in artificial intelligence, pages 411–420. Morgan Kaufmann Publishers Inc.
  * Pessoa, (2010) Pessoa, L. (2010).  Emotion and cognition and the amygdala: from “what is it?” to “what’s to be done?”.  Neuropsychologia, 48(12):3416–3429.
  * Pourahmadi et al., (2007) Pourahmadi, M., Daniels, M. J., and Park, T. (2007).  Simultaneous modelling of the Cholesky decomposition of several covariance matrices.  Journal of Multivariate Analysis, 98(3):568–587.
  * Richardson and Spirtes, (2002) Richardson, T. and Spirtes, P. (2002).  Ancestral graph markov models.  Annals of Statistics, pages 962–1030.
  * Richardson and Robins, (2013) Richardson, T. S. and Robins, J. M. (2013).  Single world intervention graphs (swigs): A unification of the counterfactual and graphical approaches to causality.  Center for the Statistics and the Social Sciences, University of Washington Series. Working Paper, 128(30):2013.
  * Robins and Greenland, (1992) Robins, J. M. and Greenland, S. (1992).  Identifiability and exchangeability for direct and indirect effects.  Epidemiology, pages 143–155.
  * Robins and Richardson, (2010) Robins, J. M. and Richardson, T. S. (2010).  Alternative graphical causal models and the identification of direct effects.  Causality and psychopathology: Finding the determinants of disorders and their cures, pages 103–158.
  * Rubin, (1980) Rubin, D. B. (1980).  Randomization analysis of experimental data: The Fisher randomization test comment.  Journal of the American Statistical Association, 75(371):591–593.
  * Rudebeck and Murray, (2014) Rudebeck, P. H. and Murray, E. A. (2014).  The orbitofrontal oracle: cortical mechanisms for the prediction and evaluation of specific behavioral outcomes.  Neuron, 84(6):1143–1156.
  * Seiler and Holmes, (2017) Seiler, C. and Holmes, S. (2017).  Multivariate heteroscedasticity models for functional brain connectivity.  Frontiers in Neuroscience, 11:696.
  * Smith et al., (2004) Smith, S. M., Jenkinson, M., Woolrich, M. W., Beckmann, C. F., Behrens, T. E., Johansen-Berg, H., Bannister, P. R., De Luca, M., Drobnjak, I., Flitney, D. E., et al. (2004).  Advances in functional and structural MR image analysis and implementation as FSL.  NeuroImage, 23:S208–S219.
  * Tibshirani et al., (2005) Tibshirani, R., Saunders, M., Rosset, S., Zhu, J., and Knight, K. (2005).  Sparsity and smoothness via the fused lasso.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67(1):91–108.
  * VanderWeele, (2015) VanderWeele, T. J. (2015).  Explanation in Causal Inference: Methods for Mediation and Interaction.  Oxford University Press.
  * VanderWeele, (2016) VanderWeele, T. J. (2016).  Mediation analysis: a practitioner’s guide.  Annual review of public health, 37:17–32.
  * Zeng et al., (2021) Zeng, S., Rosenbaum, S., Alberts, S. C., Archie, E. A., and Li, F. (2021).  Causal mediation analysis for sparse and irregular longitudinal data.  The Annals of Applied Statistics, 15(2):747–767.
  * Zhao et al., (2022) Zhao, Y., Chen, T., Cai, J., Lichenstein, S., Potenza, M. N., and Yip, S. W. (2022).  Bayesian network mediation analysis with application to the brain functional connectome.  Statistics in Medicine, 41(20):3991–4005.
  * Zhao and Luo, (2019) Zhao, Y. and Luo, X. (2019).  Granger mediation analysis of multiple time series with an application to functional magnetic resonance imaging.  Biometrics, 75(3):788–798.
  * Zhao and Luo, (2022) Zhao, Y. and Luo, X. (2022).  Pathway lasso: Pathway estimation and selection with high dimensional mediators.  Statistics and Its Interface, 15(1):39–50.
  * Zhao et al., (2018) Zhao, Y., Luo, X., Lindquist, M., and Caffo, B. (2018).  Functional mediation analysis with an application to functional magnetic resonance imaging data.  arXiv preprint arXiv:1805.06923.
  * Zhao et al., (2021) Zhao, Y., Wang, B., Mostofsky, S. H., Caffo, B. S., and Luo, X. (2021).  Covariate assisted principal regression for covariance matrix outcomes.  Biostatistics, 22(3):629–645.
  * Zou et al., (2017) Zou, T., Lan, W., Wang, H., and Tsai, C.-L. (2017).  Covariance regression analysis.  Journal of the American Statistical Association, 112(517):266–281.