# Multimodal Data Integration via Mediation Analysis with High-Dimensional
Exposures and Mediators

Yi Zhao1, Lexin Li2, for the Alzheimer’s Disease Neuroimaging
Initiative111Data used in preparation of this article were obtained from the
Alzheimer’s Disease Neuroimaging Initiative (ADNI) database
(adni.loni.usc.edu). As such, the investigators within the ADNI contributed to
the design and implementation of ADNI and/or provided data but did not
participate in analysis or writing of this report. A complete list of ADNI
investigators can be found at: http://adni.loni.usc.edu/wp-
content/uploads/how_to_apply/ADNI_Acknowledgement_List.pdf
1Department of Biostatistics and Health Data Science, Indiana University
School of Medicine
2Department of Biostatistics and Epidemiology, University of California,
Berkeley

###### Abstract

Motivated by an imaging proteomics study for Alzheimer’s disease (AD), in this
article, we propose a mediation analysis approach with high-dimensional
exposures and high-dimensional mediators to integrate data collected from
multiple platforms. The proposed method combines principal component analysis
with penalized least squares estimation for a set of linear structural
equation models. The former reduces the dimensionality and produces
uncorrelated linear combinations of the exposure variables, whereas the latter
achieves simultaneous path selection and effect estimation while allowing the
mediators to be correlated. Applying the method to the AD data identifies
numerous interesting protein peptides, brain regions, and protein-structure-
memory paths, which are in accordance with and also supplement existing
findings of AD research. Additional simulations further demonstrate the
effective empirical performance of the method.

Keywords: Alzheimer’s disease; Mediation analysis; Multimodal data
integration; Neuroimaging; Principal components analysis

##  1 Introduction

Alzheimer’s disease (AD) is an irreversible neurodegenerative disorder and is
characterized by progressive impairment of cognitive and bodily functions and
ultimate death. It is currently affecting over 5.8 million American adults
aged 65 years or older. Meanwhile, its prevalence continues to grow and is
projected to reach 13.8 million by 2050 (Alzheimer’s Association,, 2020).
Multimodal technologies have transformed AD research in recent years, by
collecting different types of data from the same group of subjects and
enabling the investigation of complex interrelated mechanisms underlying AD
development. Notable examples include multimodal neuroimaging studies of the
joint impact of brain structure and function on the disorders (Liu et al.,,
2015; Higgins et al.,, 2018), and imaging genetics studies of the impact of
genetic variants on the brain then the disease outcome (Nathoo et al.,, 2019),
among others.

Our motivation is an imaging proteomics study, which is part of the
Alzheimer’s Disease Neuroimaging Initiative (ADNI) that aims to identify
biomarkers for early detection and tracking of AD and to assist the
development of prevention and intervention strategies. Amyloid-$\beta$ is a
microscopic brain protein fragment, denotes peptides of 36 to 43 amino acids,
and is part of a larger protein called amyloid precursor protein. Tau is a
group of microtubule-associated proteins predominantly found in brain cells
and performs the function of stabilizing microtubules. Amyloid-$\beta$ is the
main component of amyloid plaques, while tau is the main component of
neurofibrillary tangles, both of which are commonly found in the brains of AD
patients. Models of AD pathophysiology hypothesize a temporal sequence, in
which accumulations of amyloid-$\beta$ plaques and neurofibrillary tangles
disrupt cell-to-cell communications and destroy brain cells, leading to brain
structural atrophy in regions such as the hippocampus, and ultimately a
clinical decline in cognition (Mormino et al.,, 2009). However, it remains
unclear how these two proteins interact with each other and with other
proteins in the cerebrospinal fluid (CSF), and how those proteins together
subsequently affect brain atrophy and disease progression. In our study, we
aim to investigate simultaneously the interrelations of multiple protein
peptides in the CSF, along with multiple brain regions of the whole brain, and
their impact on memory.

The problem can be formulated as a mediation analysis, where the goal is to
identify and explain the mechanism, or path, that underlies an observed
relationship between an exposure and an outcome variable, through the
inclusion of an intermediate variable known as a mediator. It decomposes the
effect of exposure on the outcome into a direct effect and an indirect effect,
the latter of which indicates whether the mediator is on a path from the
exposure to the outcome. In our multimodal AD study, the measurements of the
amount of multiple protein peptides serve as the exposure variables, the
volumetric measurements of multiple brain regions serve as the potential
mediators, and a composite memory score serves as the outcome. See Section 2
for more details about the study and the data. Our objective is to identify
paths from proteins to brain regional atrophies that lead to memory decline.

Mediation analysis was first proposed with a single exposure and a single
mediator (Baron and Kenny,, 1986). See VanderWeele, (2016) for a review of
mediation analysis and many references therein. In our setting, both the
exposure variables and mediators are multivariate and potentially high-
dimensional. While there have been numerous extensions of mediation analysis
to account for multiple mediators (see, e.g., Zhao and Luo,, 2016; Chén et
al.,, 2017; Song et al.,, 2018, among many others), there have been very few
works studying multivariate exposures, or both multivariate exposures and
mediators. Recently, Zhang, (2019); Aung et al., (2020); Long et al., (2020)
proposed new approaches for mediation analysis of multivariate exposures and
mediators. In particular, Zhang, (2019) developed two regularization
procedures and applied them to a mouse f2 dataset for diabetes, taking SNP
genotypes as the exposures, islet gene expressions as the mediators, and
insulin level as the outcome. However, they required the mediators to be
independent, which hardly holds in our setting, as different brain regions are
generally believed to influence each other. Aung et al., (2020) studied
environmental toxicants on pregnancy outcomes, taking toxicants as the
exposures, endogenous biomarkers such as inflammation and oxidative stress as
the mediators, and gestational age at delivery as the outcome. A key strategy
of their analysis was to reduce the exposure dimension by creating
environmental risk scores for a small number of groups based on the domain
knowledge. They showed that the between-group correlation in the reduced
exposures is negligible. However, such prior domain knowledge may not always
be available. Long et al., (2020) proposed a general mediation framework to
identify proteins that mediate the effect of metabolic gene expressions on
survival for a type of kidney cancer, taking mRNA levels as the exposures,
protein measures as the mediators, and survival time as the outcome.
Nevertheless, they implicitly required the dimensions of the exposures and
mediators cannot be too high, and thus their method is not directly applicable
to our setting, where the number of exposures and mediators can both be
potentially larger than the sample size.

In this article, we propose a mediation analysis approach, with both high-
dimensional exposures and high-dimensional mediators, for multimodal data
analysis. The method integrates principal components analysis (PCA) with
penalized least squares estimation for a set of linear structural equation
models. The former reduces the dimensionality and produces uncorrelated linear
combinations of the exposure variables, whereas the latter achieves path
selection and effect estimation while allowing the multivariate mediators to
be potentially correlated. We apply this approach to the imaging proteomics
study of AD to integrate CSF proteomics, brain volumes, and a memory measure
of mild cognitive impairment (MCI) subjects in ADNI. We identify several
interesting protein peptides, brain regions, and protein-structure-memory
paths that are in accordance with and also supplement the existing knowledge
of AD. Additional simulations further demonstrate the efficacy of the method.
Similar to Zhang, (2019); Aung et al., (2020); Long et al., (2020), our
approach is among the first attempts to conduct mediation analysis where both
the exposures and mediators are high-dimensional. But unlike the existing
solutions, we do not restrict the dimensionality or the correlation structures
and do not require additional domain knowledge of the exposures or mediators.
Moreover, although focusing on a multimodal neuroimaging study in this
article, our proposed method is equally applicable to a wide range of
multimodal data integration problems, e.g., the multi-omics data analysis
(Richardson et al.,, 2016), and the multimodal healthcare study (Cai et al.,,
2019). As such, our proposal makes a useful addition to the general toolbox of
both mediation analysis and multimodal data integration.

The rest of the article is organized as follows. Section 2 introduces the
motivating imaging proteomics data of AD. Section 3 presents the proposed
model and estimation approach. Section 4 analyzes the AD dataset, with a
detailed discussion on the identified protein peptides, brain regions, and
path. Section 5 complements with additional simulation results to demonstrate
the empirical performance of the method.

##  2 Motivating example

While Alzheimer’s disease is becoming a major public health challenge as the
population ages, there is no effective treatment for AD that is capable of
stopping or slowing the associated cognitive and neuronal degradation.
Therefore, understanding the disease pathology, identifying biological
markers, and finding early diagnosis and intervention strategies are of
critical importance (Alzheimer’s Association,, 2020). Among numerous AD-
related proteins in the CSF, amyloid-$\beta$ and tau are two major proteins
that are consistently identified in the brains of AD patients, and their
abnormal abundance generally indicates AD pathology (Jagust,, 2018). Even
though there has been evidence suggesting a pathological connection between
amyloid-$\beta$ deposition, hippocampus atrophy, and memory decline (Mormino
et al.,, 2009), it remains largely unknown how amyloid-$\beta$ and tau
interact with each other, how they interact with other proteins in the CSF,
and how these proteins together affect the downstream brain atrophy and
cognitive outcome. In our study, we aim to delineate the regulatory
relationships among multiple CSF proteins, structural atrophy of the whole
brain, and cognitive behavior, and to identify important biological paths.

The data used in our study are obtained from the Alzheimer’s Disease
Neuroimaging Initiative (ADNI, adni.loni.usc.edu). The CSF proteomics data
were obtained using targeted liquid chromatography multiple reaction
monitoring mass spectrometry, which is a highly specific, sensitive, and
reproducible technique for quantifying targeted proteins. A list of protein
fragments, or peptides, was sent to the detector. The samples then went
through peak integration, outliers detection, normalization, quantification,
and quality control using test/re-test samples. This procedure results in the
intensity measures of 320 peptides that are annotated from 142 proteins. The
brain imaging data were obtained using anatomical magnetic resonance imaging
(MRI). Each image was first preprocessed following the standard pipeline, then
mapped to an atlas consisting of 145 brain regions-of-interest to extract the
volumetric measures (Doshi et al.,, 2016). The atlas used in the study spans
the entire brain and was actually built on multiple atlases. Individual
atlases were first warped to the target image using a nonlinear registration
method, followed by a spatially adaptive weighted voting strategy to fuse into
a final segmentation. Moreover, the volume of each brain region was
standardized by the total intracranial volume to account for variations of
individual brain size. The cognitive outcome is a composite memory score,
ADNI-MEM, that involves a battery of neuropsychological tests. In our study,
we focus on 135 subjects diagnosed as mild cognitive impairment (MCI) patients
at recruitment. MCI is a prodromal stage of AD, with a slight but noticeable
and measurable decline in cognitive abilities. A person with MCI is at an
increased risk of developing AD or other dementia. Understanding the
pathologic mechanism underlying MCI provides important clues of onset of the
disorder as well as a useful guide for early diagnosis and intervention.

##  3 Model and Method

We first present the proposed model, then an estimation method integrating
principal components analysis and penalized estimation.

###  3.1 Model

Suppose there are totally $n$ subjects. Let
$\mathbf{X}_{i}=(X_{i1},\ldots,X_{ir})^{\top}\in\mathbb{R}^{r}$ denote the
$r$-dimensional vector of exposure variables,
$\mathbf{M}_{i}=(M_{i1},\ldots,M_{ip})^{\top}\in\mathbb{R}^{p}$ denote the
$p$-dimensional vector of mediators, and $Y_{i}\in\mathbb{R}$ denote the
univariate outcome variable, for subjects $i=1,\dots,n$. In our imaging
proteomics study, $\mathbf{X}_{i}$ denotes the protein peptide measures with
$r=320$, $\mathbf{M}_{i}$ denotes the brain volumetric measures with $p=145$,
$Y_{i}$ denotes the memory score, and the sample size $n=135$.

Figure 1: The schematic diagram of the proposed model with $q$ exposure
variables $\tilde{X}_{1},\dots,\tilde{X}_{q}$, $p$ mediators
$M_{1},\dots,M_{p}$, and the outcome variable $Y$.

The first step of our method is to perform a principal components analysis on
$\mathbf{X}_{i}$ to produce uncorrelated composite exposures. If
$\mathbf{X}_{i}$ further follows a multivariate normal distribution, then the
produced composite exposures are independent. Let
$\tilde{\mathbf{X}}_{i}=(\tilde{X}_{i1},\ldots,\tilde{X}_{iq})^{\top}\in\mathbb{R}^{q}$
denote the first $q$ principal components. We then continue to model the path
relations among $\tilde{\mathbf{X}}_{i},\mathbf{M}_{i}$ and $Y_{i}$ via the
following set of linear structural equation models,

$\displaystyle\begin{split}\mathbf{M}&=\tilde{\mathbf{X}}\boldsymbol{\alpha}+\boldsymbol{\epsilon},\\\
\mathbf{Y}&=\tilde{\mathbf{X}}\boldsymbol{\gamma}+\mathbf{M}\boldsymbol{\beta}+\boldsymbol{\eta},\end{split}$
(1)

where
$\tilde{\mathbf{X}}=(\tilde{\mathbf{X}}_{1},\dots,\tilde{\mathbf{X}}_{n})^{\top}\in\mathbb{R}^{n\times
q}$,
$\mathbf{M}=(\mathbf{M}_{1},\dots,\mathbf{M}_{n})^{\top}\in\mathbb{R}^{n\times
p}$, $\mathbf{Y}=(Y_{1},\dots,Y_{n})^{\top}\in\mathbb{R}^{n}$ stack the
composite exposures, mediators, and outcome across all subjects, respectively,
$\boldsymbol{\epsilon}=(\boldsymbol{\epsilon}_{1},\dots,\boldsymbol{\epsilon}_{n})^{\top}\in\mathbb{R}^{n\times
p}$, with
$\boldsymbol{\epsilon}_{i}=(\epsilon_{i1},\dots,\epsilon_{ip})^{\top}\in\mathbb{R}^{p}$,
and $\boldsymbol{\eta}=(\eta_{1},\dots,\eta_{n})^{\top}\in\mathbb{R}^{n}$ are
measurement errors. Suppose both error terms follow some zero mean normal
distribution, and $\boldsymbol{\epsilon}$ is independent of
$\tilde{\mathbf{X}}$, $\boldsymbol{\eta}$ is independent of
$\tilde{\mathbf{X}}$ and $\mathbf{M}$, and $\boldsymbol{\epsilon}$ and
$\boldsymbol{\eta}$ are independent of each other. The parameters
$\boldsymbol{\alpha}=(\alpha_{jk})\in\mathbb{R}^{q\times p}$,
$\boldsymbol{\beta}=(\beta_{1},\dots,\beta_{p})^{\top}\in\mathbb{R}^{p}$, and
$\boldsymbol{\gamma}=(\gamma_{1},\dots,\gamma_{q})^{\top}\in\mathbb{R}^{q}$
capture the path effects. Model (1) is similar to that used in Zhao and Luo,
(2016); Zhao et al., (2020), but none of those can handle multivariate
exposure variables. Besides, we introduce some different forms of penalty
functions in our parameter estimation.

Figure 1 shows a schematic description of model (1). Under this model, we
define the direct effect of $\tilde{X}_{j}$ on $Y$ as
$\mathrm{DE}(\tilde{X}_{j})=\gamma_{j}$, the indirect effect of
$\tilde{X}_{j}$ on $Y$ through $M_{k}$ as
$\mathrm{IE}(\tilde{X}_{j},M_{k})=\alpha_{jk}\beta_{k}$, and the total
indirect effect of $\tilde{X}_{j}$ on $Y$ as
$\mathrm{IE}(\tilde{X}_{j})=\sum_{k=1}^{p}\alpha_{jk}\beta_{k}$, for
$j=1,\dots,q$. The total effect of $\tilde{X}_{j}$ satisfies that
$\mathrm{TE}(\tilde{X}_{j})=\mathrm{IE}(\tilde{X}_{j})+\mathrm{DE}(\tilde{X}_{j})=\sum_{k=1}^{p}\alpha_{jk}\beta_{k}+\gamma_{j}$.

A key characteristic of model (1) is that it allows the multivariate mediators
to be conditionally dependent given the exposures. To better illustrate this,
we consider a simple example of model (1), where $q=1,p=3$, as shown in Figure
2. In this example, Figure 2(a) outlines the sequential influences among all
the mediators, while Figure 2(b) is the proposed model (1). We see that, for
the first mediator, $M_{1}$, $\alpha_{11}=a_{11},\beta_{1}=b_{1}$; for the
second mediator, $M_{2}$, $\alpha_{12}=a_{11}d_{12}+a_{12},\beta_{2}=b_{2}$;
and for the third mediator, $M_{3}$,
$\alpha_{13}=a_{11}d_{13}+a_{11}d_{12}d_{23}+a_{12}d_{23}+a_{13},\beta_{3}=b_{3}$.
As such, $\alpha_{1k}$ consolidates the effects through the $k$th mediator
$M_{k}$, and the indirect effect
$\mathrm{IE}(\tilde{X}_{1},M_{k})=\alpha_{1k}\beta_{k}$ can be viewed as the
consolidated indirect effect through $M_{k}$, $k=1,2,3$.

((a))  ((b))

Figure 2: A model example with $q=1$ exposure variable and $p=3$ sequentially
ordered mediators.

###  3.2 Estimation

We propose to estimate the parameters in model (1) through the penalized
ordinary least squares,

$\underset{\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma}}{\text{minimize}}\;\;\frac{1}{2}\mathcal{L}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma})+\lambda_{1}\mathcal{R}_{1}(\boldsymbol{\alpha},\boldsymbol{\beta})+\lambda_{2}\mathcal{R}_{2}(\boldsymbol{\alpha},\boldsymbol{\beta})+\lambda_{3}\mathcal{R}_{3}(\boldsymbol{\gamma}),$
(2)

where the loss function is the usual least squares loss,

$\mathcal{L}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma})=\mathrm{tr}\left\\{(\mathbf{M}-\tilde{\mathbf{X}}\boldsymbol{\alpha})^{\top}(\mathbf{M}-\tilde{\mathbf{X}}\boldsymbol{\alpha})\right\\}+(\mathbf{Y}-\tilde{\mathbf{X}}\boldsymbol{\gamma}-\mathbf{M}\boldsymbol{\beta})^{\top}(\mathbf{Y}-\tilde{\mathbf{X}}\boldsymbol{\gamma}-\mathbf{M}\boldsymbol{\beta}).$

$\mathcal{R}_{1},\mathcal{R}_{2},\mathcal{R}_{3}$ are three penalty functions,
with the tuning parameters $\lambda_{1},\lambda_{2},\lambda_{3}$,
respectively. We next discuss each penalty function in detail.

The first penalty function $\mathcal{R}_{1}$ is of the form,

$\mathcal{R}_{1}(\boldsymbol{\alpha},\boldsymbol{\beta})=\sum_{j=1}^{q}\sum_{k=1}^{p}\left\\{|\alpha_{jk}\beta_{k}|+c_{0}\left(\alpha_{jk}^{2}+\beta_{k}^{2}\right)\right\\}+c_{1}\left(\sum_{j=1}^{q}\sum_{k=1}^{p}|\alpha_{jk}|+\sum_{k=1}^{p}|\beta_{k}|\right),$

for some parameters $c_{0}$ and $c_{1}$. It is a generalization of the pathway
Lasso penalty of Zhao and Luo, (2016) to $q$ exposure variables, and is to
facilitate selection of individual mediators. Specifically, for a given
mediator $M_{k}$, the term $\sum_{j=1}^{q}|\alpha_{jk}\beta_{k}|$ is a product
Lasso penalty, and encourages all the paths going through $M_{k}$ to be shrunk
to zero, which in effect achieves the goal of mediator selection. The term
$c_{0}(\alpha_{jk}^{2}+\beta_{k}^{2})$ is to make the penalty a convex
function, with a proper choice of the parameter $c_{0}$. It is straightforward
to show that, when $c_{0}\geq 1/2$, the sum
$|\alpha_{jk}\beta_{k}|+c_{0}(\alpha_{jk}^{2}+\beta_{k}^{2})$ is convex. In
our implementation, we fix $c_{0}=2$. The last term in $\mathcal{R}_{1}$ is
the sum of usual Lasso penalty that further penalizes individual path effects
$\alpha_{jk},\beta_{k}$, with $c_{1}$ being an additional tuning parameter. It
is found that this additional penalty helps further improves the selection
accuracy (Zhao and Luo,, 2016).

The second penalty function $\mathcal{R}_{2}$ is of the form,

$\mathcal{R}_{2}(\boldsymbol{\alpha},\boldsymbol{\beta})=\sum_{j=1}^{q}\sqrt{p}\sqrt{\sum_{k=1}^{p}(\alpha_{jk}\beta_{k})^{2}}.$

It is a group Lasso penalty and is to facilitate the selection of individual
exposure. Specifically, for a given exposure $\tilde{X}_{j}$, the penalty
$\\{\sum_{k=1}^{p}(\alpha_{jk}\beta_{k})^{2}\\}^{1/2}$ encourages all the
paths originating from $\tilde{X}_{j}$ to be shrunk to zero, which in effect
achieves the goal of exposure selection.

The third penalty function $\mathcal{R}_{3}$ is of the form,

$\mathcal{R}_{3}(\boldsymbol{\gamma})=\sum_{j=1}^{q}|\gamma_{j}|.$

This is simply the usual Lasso penalty and is to facilitate selection of
direct effects between the exposures and the outcome.

We next discuss how to solve the minimization problem (2). We note that (2)
involves the penalties on the product terms $\alpha_{jk}\beta_{k}$, making it
difficult to derive the analytical solutions. As such, we first introduce a
new parameter, $\mu_{jk}=\alpha_{jk}\beta_{k}$, which turns (2) to an
equivalent problem of solving a sparse group lasso that has an explicit form
of solution (Simon et al.,, 2013). That is, letting
$\boldsymbol{\mu}=(\mu_{jk})\in\mathbb{R}^{q\times p}$, we turn to the
equivalent optimization problem,

$\displaystyle\begin{split}\underset{\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{\mu}}{\text{minimize}}\;\;\frac{1}{2}\mathcal{L}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma})+\lambda_{1}\mathcal{R}_{1}(\boldsymbol{\mu},\boldsymbol{\alpha},\boldsymbol{\beta})+\lambda_{2}\mathcal{R}_{2}(\boldsymbol{\mu})+\lambda_{3}\mathcal{R}_{3}(\boldsymbol{\gamma}),\\\
\text{ such that }\;\;\mu_{jk}=\alpha_{jk}\beta_{k},\;\;\text{ for
}j=1,\ldots,q\text{ and }k=1,\ldots,p.\end{split}$ (3)

Let $\boldsymbol{\mu}_{j}=(\mu_{j1},\dots,\mu_{jp})^{\top}\in\mathbb{R}^{p}$,
$\boldsymbol{\alpha}_{j}=(\alpha_{j1},\dots,\alpha_{jp})^{\top}\in\mathbb{R}^{p}$,
and introduce the augmented Lagrangian parameter
$\boldsymbol{\tau}_{j}=(\tau_{j1},\dots,\tau_{jp})^{\top}\in\mathbb{R}^{p}$,
for $j=1,\dots,q$, and $\boldsymbol{\tau}=(\tau_{jk})\in\mathbb{R}^{q\times
p}$. Then, the augmented Lagrangian form of (3) is,

$\displaystyle\begin{split}\underset{\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{\mu},\boldsymbol{\tau}}{\text{minimize}}\;\;\frac{1}{2}\mathcal{L}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma})+\lambda_{1}\mathcal{R}_{1}(\boldsymbol{\mu},\boldsymbol{\alpha},\boldsymbol{\beta})+\lambda_{2}\mathcal{R}_{2}(\boldsymbol{\mu})+\lambda_{3}\mathcal{R}_{3}(\boldsymbol{\gamma})\\\
+\sum_{j=1}^{q}\left(\langle\boldsymbol{\mu}_{j}-\boldsymbol{\alpha}_{j}\circ\boldsymbol{\beta},\boldsymbol{\tau}_{j}\rangle+\frac{\rho}{2}\|\boldsymbol{\mu}_{j}-\boldsymbol{\alpha}_{j}\circ\boldsymbol{\beta}\|_{2}^{2}\right),\end{split}$
(4)

where $\rho>0$ is the augmented Lagrangian constant that we set $\rho=1$ in
our implementation, $\circ$ is the Hadamard product,
$\langle\cdot,\cdot\rangle$ is the inner product, and $\|\cdot\|_{2}$ is the
$L_{2}$-norm. We next solve (4) by updating
$\boldsymbol{\mu},\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma}$
and $\boldsymbol{\tau}$ iteratively.

More specifically, we first fix
$\boldsymbol{\alpha}^{(s)},\boldsymbol{\beta}^{(s)},\boldsymbol{\gamma}^{(s)},\boldsymbol{\tau}^{(s)}$
at iteration $s$, and update $\boldsymbol{\mu}_{j}$ by solving

$\underset{\boldsymbol{\mu}_{j}}{\text{minimize}}~{}\frac{\rho}{2}\|\boldsymbol{\mu}_{j}-\boldsymbol{\alpha}_{j}^{(s)}\circ\boldsymbol{\beta}^{(s)}\|_{2}^{2}+\boldsymbol{\tau}_{j}^{(s)\top}\left(\boldsymbol{\mu}_{j}-\boldsymbol{\alpha}_{j}^{(s)}\circ\boldsymbol{\beta}^{(s)}\right)+\lambda_{1}\|\boldsymbol{\mu}_{j}\|_{1}+\lambda_{2}\sqrt{p}\|\boldsymbol{\mu}_{j}\|_{2},$

for $j=1,\ldots,q$, where $\|\cdot\|_{1}$ is the $L_{1}$-norm. There is a
closed-form solution,

$\begin{split}\mu_{jk}^{(s+1)}=\begin{cases}\left\\{\|\mathcal{S}(\boldsymbol{\nu}_{j},\lambda_{1}/\rho)\|_{2}-\lambda_{2}\sqrt{p}/\rho\right\\}_{+}\frac{\mathcal{S}(\nu_{jk},\lambda_{1}/\rho)}{\|\mathcal{S}(\boldsymbol{\nu}_{j},\lambda_{1}/\rho)\|_{2}},&\textrm{
if }\;\;\|\mathcal{S}(\boldsymbol{\nu}_{j},\lambda_{1}/\rho)\|_{2}\neq 0,\\\
0&\textrm{ otherwise},\end{cases}\end{split}$ (5)

for $j=1,\ldots,q,k=1,\ldots,p$, where
$\boldsymbol{\nu}_{j}=\boldsymbol{\alpha}_{j}^{(s)}\circ\boldsymbol{\beta}^{(s)}-\boldsymbol{\tau}_{j}^{(s)}/\rho$,
$\mathcal{S}(a,\lambda)=\mathrm{sgn}(a)\max\\{|a|-\lambda,0\\}$ is the soft-
thresholding function with $\mathrm{sgn}(a)$ denoting the sign of $a$ and
$a_{+}=\max\\{a,0\\}$, and $\mathcal{S}(\mathbf{a},\lambda)$ denotes the
element-wise soft-thresholding of a vector $\mathbf{a}$.

We next fix
$\boldsymbol{\mu}^{(s+1)},\boldsymbol{\beta}^{(s)},\boldsymbol{\gamma}^{(s)},\boldsymbol{\tau}^{(s)}$,
and update $\boldsymbol{\alpha}_{j}$ by solving

$\underset{\boldsymbol{\alpha}_{j}}{\text{minimize}}~{}\mathbf{V}_{j}\boldsymbol{\alpha}_{j}+\lambda_{1}c_{1}\mathrm{sgn}(\boldsymbol{\alpha}_{j})-\mathbf{w}_{j},$

where
$\mathbf{V}_{j}=\rho\mathbf{D}_{\boldsymbol{\beta}^{(s)}}^{2}+(4\lambda_{1}+\tilde{\mathbf{x}}_{j}^{\top}\tilde{\mathbf{x}}_{j})\boldsymbol{\mathrm{I}}_{p}$,
$\mathbf{w}_{j}=(\mathbf{M}-\sum_{l\neq
j}\tilde{\mathbf{x}}_{l}\boldsymbol{\alpha}_{l}^{(s)\top})^{\top}\tilde{\mathbf{x}}_{j}+\mathbf{D}_{\boldsymbol{\beta}^{(s)}}\boldsymbol{\tau}_{j}^{(s)}+\rho\mathbf{D}_{\boldsymbol{\beta}^{(s)}}\boldsymbol{\mu}_{j}^{(s+1)}$,
$\mathbf{D}_{\boldsymbol{\beta}^{(s)}}$ is a diagonal matrix with
$\boldsymbol{\beta}^{(s)}$ as the diagonal elements,
$\tilde{\mathbf{x}}_{j}\in\mathbb{R}^{n}$ is the $j$th column of
$\tilde{\mathbf{X}}$, and $\boldsymbol{\mathrm{I}}_{p}$ is the $p$-dimensional
identity matrix. The solution is

$\boldsymbol{\alpha}_{j}^{(s+1)}=\mathbf{V}_{j}^{-1}\mathcal{S}(\mathbf{w}_{j},\lambda_{1}c_{1}),\quad
j=1,\ldots,q.$ (6)

We next fix
$\boldsymbol{\mu}^{(s+1)},\boldsymbol{\alpha}^{(s+1)},\boldsymbol{\gamma}^{(s)},\boldsymbol{\tau}^{(s)}$,
and update $\boldsymbol{\beta}$ by solving

$\underset{\boldsymbol{\beta}}{\text{minimize}}~{}\mathbf{V}_{\boldsymbol{\beta}}\boldsymbol{\beta}+\lambda_{1}\mathrm{sgn}(\boldsymbol{\beta})-\mathbf{w}_{\boldsymbol{\beta}},$

where
$\mathbf{V}_{\boldsymbol{\beta}}=\mathbf{M}^{\top}\mathbf{M}+\rho\sum_{j=1}^{q}\mathbf{D}_{\boldsymbol{\alpha}_{j}^{(s+1)}}^{2}+4\lambda_{1}q\boldsymbol{\mathrm{I}}_{p}$,
$\mathbf{w}_{\boldsymbol{\beta}}=\mathbf{M}^{\top}(\mathbf{Y}-\tilde{\mathbf{X}}\boldsymbol{\gamma}^{(s)})+\sum_{j=1}^{q}\mathbf{D}_{\boldsymbol{\alpha}_{j}^{(s+1)}}\boldsymbol{\tau}_{j}^{(s)}$
$+\rho\sum_{j=1}^{q}\mathbf{D}_{\boldsymbol{\alpha}_{j}^{(s+1)}}\boldsymbol{\mu}_{j}^{(s+1)}$,
and $\mathbf{D}_{\boldsymbol{\alpha}_{j}^{(s+1)}}$ is a diagonal matrix with
$\boldsymbol{\alpha}_{j}^{(s+1)}$ as the diagonal elements. The solution is

$\boldsymbol{\beta}^{(s+1)}=\mathbf{V}_{\boldsymbol{\beta}}^{-1}\mathcal{S}(\mathbf{w}_{\boldsymbol{\beta}},\lambda_{1}c_{1}).$
(7)

We then fix
$\boldsymbol{\mu}^{(s+1)},\boldsymbol{\alpha}^{(s+1)},\boldsymbol{\beta}^{(s+1)},\boldsymbol{\tau}^{(s)}$,
and update $\boldsymbol{\gamma}$ by solving

$\underset{\boldsymbol{\gamma}}{\text{minimize}}~{}\mathbf{V}_{\boldsymbol{\gamma}}\boldsymbol{\gamma}+\lambda_{3}\mathrm{sgn}(\boldsymbol{\gamma})-\mathbf{w}_{\boldsymbol{\gamma}},$

where
$\mathbf{V}_{\boldsymbol{\gamma}}=\tilde{\mathbf{X}}^{\top}\tilde{\mathbf{X}}$
and
$\mathbf{w}_{\boldsymbol{\gamma}}=\tilde{\mathbf{X}}^{\top}(\mathbf{Y}-\mathbf{M}\boldsymbol{\beta}^{(s+1)})$.
The solution is

$\boldsymbol{\gamma}^{(s+1)}=\mathbf{V}_{\boldsymbol{\gamma}}^{-1}\mathcal{S}(\mathbf{w}_{\boldsymbol{\gamma}},\lambda_{3}).$
(8)

Finally, we fix
$\boldsymbol{\mu}^{(s+1)},\boldsymbol{\alpha}^{(s+1)},\boldsymbol{\beta}^{(s+1)},\boldsymbol{\gamma}^{(s+1)}$,
and update $\boldsymbol{\tau}$ by

$\boldsymbol{\tau}_{j}^{(s+1)}=\boldsymbol{\tau}_{j}^{(s)}+\rho\left(\boldsymbol{\mu}_{j}^{(s+1)}-\boldsymbol{\alpha}_{j}^{(s+1)}\circ\boldsymbol{\beta}^{(s+1)}\right),\quad
j=1,\ldots q.$ (9)

We stop the iterations until some stopping criterion is met. In our implement,
we stop when the difference of two consecutive objective values is smaller
than $10^{-6}$. We summarize the above optimization procedure in Algorithm 1.

Algorithm 1 The optimization algorithm for (4).

1:$(\tilde{\mathbf{X}},\mathbf{M},\mathbf{Y})$ and the tuning parameters
$\lambda_{1},\lambda_{2},\lambda_{3}$ and $c_{1}$.

2:initialization:
$\left\\{\boldsymbol{\alpha}^{(0)},\boldsymbol{\beta}^{(0)},\boldsymbol{\gamma}^{(0)},\boldsymbol{\mu}^{(0)},\boldsymbol{\tau}^{(0)}\right\\}$.

3:repeat

4: update $\mu_{jk}^{(s+1)}$ given
$\boldsymbol{\alpha}^{(s)},\boldsymbol{\beta}^{(s)},\boldsymbol{\gamma}^{(s)},\boldsymbol{\tau}^{(s)}$
by (5), for $j=1,\ldots,q,k=1,\ldots,p$.

5: update $\boldsymbol{\alpha}_{j}^{(s+1)}$ given
$\boldsymbol{\mu}^{(s+1)},\boldsymbol{\beta}^{(s)},\boldsymbol{\gamma}^{(s)},\boldsymbol{\tau}^{(s)}$
by (6), for $j=1,\ldots,q$.

6: update $\boldsymbol{\beta}^{(s+1)}$ given
$\boldsymbol{\mu}^{(s+1)},\boldsymbol{\alpha}^{(s+1)},\boldsymbol{\gamma}^{(s)},\boldsymbol{\tau}^{(s)}$
by (7).

7: update $\boldsymbol{\gamma}^{(s+1)}$ given
$\boldsymbol{\mu}^{(s+1)},\boldsymbol{\alpha}^{(s+1)},\boldsymbol{\beta}^{(s+1)},\boldsymbol{\tau}^{(s)}$
by (8).

8: update $\boldsymbol{\tau}_{j}^{(s+1)}$ given
$\boldsymbol{\mu}^{(s+1)},\boldsymbol{\alpha}^{(s+1)},\boldsymbol{\beta}^{(s+1)},\boldsymbol{\gamma}^{(s+1)}$
by (9), for $j=1,\ldots,q$.

9:until the stopping criterion is met.

10:$\left\\{\hat{\boldsymbol{\alpha}},\hat{\boldsymbol{\beta}},\hat{\boldsymbol{\gamma}}\right\\}$.

We tune the parameters in (4) using the Bayesian information criterion (BIC),

$\mathrm{BIC}=-2\log{\mathcal{L}\left(\hat{\boldsymbol{\alpha}},\hat{\boldsymbol{\beta}},\hat{\boldsymbol{\gamma}}\right)}+\log(n)|\hat{\mathcal{A}}|,$

where
$\hat{\boldsymbol{\alpha}},\hat{\boldsymbol{\beta}},\hat{\boldsymbol{\gamma}}$
are the estimates under a given set of tuning parameters
$\lambda_{1},\lambda_{2},\lambda_{3}$ and $c_{1}$,
$\mathcal{A}=\left\\{(j,k):\alpha_{jk}\beta_{k}\neq 0\right\\}$ denotes the
active set, and $|\mathcal{A}|$ is the cardinality. In our implementation, we
adopt the tuning strategy of Zou and Hastie, (2005), by tuning the ratios
$\lambda_{2}/\lambda_{1},\lambda_{3}/\lambda_{1}$ along with $c_{1}$ in a grid
search, and choose the best set of parameters that minimizes the BIC.

##  4 AD Imaging Proteomics Study Revisited

We apply the proposed method to the ADNI imaging proteomics data, taking the
CSF peptide measures as the exposures, the brain volumetric measures as the
mediators, and the memory score as the outcome. Moreover, we adjust the
exposures, mediators, and outcome for age, gender, ApoE4, and years of
education to remove potential confounding effects (Rosenbaum,, 2002). We first
summarize the identified paths with nonzero effects, then discuss the relevant
proteins and brain regions in detail. In summary, our findings are consistent
with the existing knowledge of AD. Moreover, our method also suggests a few
potentially interesting protein-structure-memory paths that may deserve
further examination and verification.

###  4.1 Paths with nonzero effects

We first apply principal components analysis to the peptide data. The top 20
principal components (PCs) account for about $85\%$ of total data variation.
We thus focus on those $q=20$ top PCs and feed them as the exposure variables
into the subsequent penalized path analysis. Figure 3 presents all the
identified paths with a nonzero indirect path effect. Table 1 presents the
estimated path effects including the estimated $\alpha$ and $\beta$ of each
path, and Table 2 presents the indirect, direct, and total effect of each
exposure PC.

Figure 3:  The estimated paths for the AD imaging proteomics study. The red
nodes denote the principal components of the peptides as exposures, the green
nodes the brain regions as mediators, and the blue node the memory score as
outcome. The red arrows indicate positive path effects, and the blue arrows
negative path effects. Table 1:  The brain regions with a nonzero indirect
effect ($\mathrm{IE}=\alpha\beta$) for the AD imaging proteomics study.

|  |  | Principal components of peptides as exposures |
---|---|---|---|---
| Brain regions as mediators |  | PC1 | PC2 | PC4 | PC5 | PC7 | PC9 | PC19 | $\beta$ ($\times 10^{-2}$)
|  | $\alpha$ | - | - | - | -0.17 | - | - | - |
R41 | Left cerebellum white matter |  IE ($\times 10^{-3}$) | - | - | - | -1.30 | - | - | - | 0.76
|  | $\alpha$ | - | - | 0.11 | 0.13 | 0.13 | - | - |
R47 | Right hippocampus |  IE ($\times 10^{-3}$) | - | - | 1.12 | 1.60 | 1.52 | - | - | 1.17
|  | $\alpha$ | 0.11 | - | 0.13 | 0.13 | 0.22 | 0.12 | - |
R48 | Left hippocampus |  IE ($\times 10^{-3}$) | 1.34 | - | 1.59 | 1.76 | 3.40 | 1.55 | - | 1.20
|  | $\alpha$ | -0.25 | 0.15 | -0.26 | -0.18 | - | -0.16 | - |
R49 | Temporal horn of right lateral ventricle |  IE ($\times 10^{-3}$) | 2.06 | -1.01 | 2.03 | 1.28 | - | 1.08 | - | -0.66
|  | $\alpha$ | -0.29 | - | -0.23 | -0.25 | - | -0.21 | - |
R50 | Temporal horn of left lateral ventricle |  IE ($\times 10^{-3}$) | 2.55 | - | 1.78 | 2.05 | - | 1.74 | - | -0.71
|  | $\alpha$ | -0.36 | - | - | - | - | - | - |
R51 | Right lateral ventricle |  IE ($\times 10^{-3}$) | 1.06 | - | - | - | - | - | - | -0.30
|  | $\alpha$ | -0.36 | - | - | - | - | - | - |
R52 | Left lateral ventricle |  IE ($\times 10^{-3}$) | 1.15 | - | - | - | - | - | - | -0.27
|  | $\alpha$ | - | 0.14 | - | - | - | - | - |
R73 | Cerebellar vermal lobules VIII-X |  IE ($\times 10^{-3}$) | - | 1.88 | - | - | - | - | - | 1.00
|  | $\alpha$ | - | -0.17 | 0.20 | 0.25 | - | - | - |
R103 | Left anterior insula |  IE ($\times 10^{-3}$) | - | -1.13 | 1.27 | 1.76 | - | - | - | 0.56
|  | $\alpha$ | - | - | - | - | 0.15 | - | -0.18 |
R106 | Right angular gyrus |  IE ($\times 10^{-3}$) | - | - | - | - | 1.03 | - | -1.41 | 0.78
|  | $\alpha$ | - | - | - | 0.17 | - | - | - |
R117 | Left entorhinal areas |  IE ($\times 10^{-3}$) | - | - | - | 1.15 | - | - | - | 0.76
|  | $\alpha$ | - | - | 0.16 | - | - | - | - |
R120 | Right frontal pole |  IE ($\times 10^{-3}$) | - | - | 1.12 | - | - | - | - | 0.75
|  | $\alpha$ | - | - | - | - | - | - | -0.16 |
R121 | Left frontal pole |  IE ($\times 10^{-3}$) | - | - | - | - | - | - | -1.12 | 0.77
|  | $\alpha$ | 0.16 | - | - | - | - | - | - |
R122 | Right fusiform gyrus |  IE ($\times 10^{-3}$) | 1.32 | - | - | - | - | - | - | 1.02
|  | $\alpha$ | 0.19 | - | - | - | - | - | - |
R123 | Left fusiform gyrus |  IE ($\times 10^{-3}$) | 1.12 | - | - | - | - | - | - | 0.66
|  | $\alpha$ | - | - | 0.21 | - | - | 0.19 | - |
R154 | Right middle temporal gyrus |  IE ($\times 10^{-3}$) | - | - | 1.68 | - | - | 1.66 | - | 0.74
|  | $\alpha$ | 0.13 | - | 0.14 | 0.18 | 0.14 | 0.15 | - |
R155 | Left middle temporal gyrus |  IE ($\times 10^{-3}$) | 1.01 | - | 1.09 | 1.63 | 1.05 | 1.35 | - | 0.79
|  | $\alpha$ | - | - | 0.15 | - | - | - | - |
R169 | Left precuneus |  IE ($\times 10^{-3}$) | - | - | 1.10 | - | - | - | - | 0.82
|  | $\alpha$ | 0.13 | -0.12 | 0.22 | 0.11 | 0.11 | - | - |
R172 | Right posterior insula |  IE ($\times 10^{-3}$) | 1.44 | -1.30 | 2.67 | 1.03 | 1.00 | - | - | 1.03
|  | $\alpha$ | - | -0.17 | 0.24 | 0.15 | - | - | - |
R173 | Left posterior insula |  IE ($\times 10^{-3}$) | - | -1.46 | 2.18 | 1.09 | - | - | - | 0.82
|  | $\alpha$ | - | - | - | - | - | - | -0.24 |
R182 | Right precentral gyrus |  IE ($\times 10^{-3}$) | - | - | - | - | - | - | 1.91 | -0.51

Table 2:  The estimated indirect effects (IE), direct effects (DE), and total
effects (TE) of the top principal components. The PCs with zero IE and DE are
not presented in the table.

| PC1 | PC2 | PC4 | PC5 | PC6 | PC7 | PC9 | PC11 | PC14 | PC15 | PC16 | PC19 | Total
---|---|---|---|---|---|---|---|---|---|---|---|---|---
IE | 0.013 | -0.003 | 0.018 | 0.012 | - | 0.008 | 0.007 | - | - | - | - | -0.001 | 0.054
DE | 0.138 | - | - | 0.066 | -0.035 | 0.168 | 0.065 | -0.018 | 0.102 | -0.007 | 0.156 | - | 0.634
TE | 0.151 | -0.003 | 0.018 | 0.078 | -0.035 | 0.176 | 0.072 | -0.018 | 0.102 | -0.007 | 0.156 | -0.001 | 0.688

###  4.2 Proteins

Among the 20 PCs, seven have nonzero indirect effects on memory. Next, we
focus on PC1, PC4, and PC5 as they account for a higher proportion of total
data variation and demonstrate a relatively higher indirect path effect on the
outcome. Table 3 lists the top proteins in PC1, PC4, and PC5, and the
corresponding gene name. We also include the regulation directions found in
the AD literature, where an upregulation compared to cognitive normal controls
indicates a higher protein abundance in MCI/AD patients, as well as the
direction of correlations with the CSF amyloid-$\beta$ and tau, the two well-
established AD protein biomarkers (Wesenhagen et al.,, 2020). We next discuss
the identified proteins by their relevance in the amyloid-$\beta$ and tau
pathology.

Table 3: Proteins with top loading magnitude in PC1, PC4 and PC5. For each
protein, direction of protein level in MCI/AD compared to normal control and
correlation with CSF tau and amyloid reported in the literature are provided.

|  |  |  | Correlation
---|---|---|---|---
Protein | Loading | Gene | Direction | tau | amyloid
PC1 |  |  |  |  |
ProSAAS | 0.075 | PCSK1N | $\Updownarrow$ |  |
Neuronal growth regulator 1 | 0.075 | NEGR1 | $\downarrow$ |  |
Cell adhesion molecule 3 | 0.075 | CADM3 | $\downarrow$ |  |
Neuroblastoma suppressor of tumorigenicity 1 | 0.073 | NBL1 | $\uparrow$ |  |
Spondin-1 | 0.073 | SPON1 | $\uparrow$ | $\uparrow$ | $\downarrow$
Prostagiandin-H2 D-isomerase | 0.073 | PTGDS | $\downarrow$ |  | $\downarrow$
Monocyte differentiation antigen CD14 | 0.071 | CD14 | $\uparrow$ |  |
VPS10 domain-containing receptor SorCS1 | 0.069 | SORCS1 |  | $\uparrow$ | $\downarrow$
PC4 |  |  |  |  |
Neuronal pentraxin-2 | 0.152 | NPTX2 | $\downarrow$ |  | $\uparrow$
Insulin-like growth factor-binding protein 2 | -0.146 | IGFBP2 | $\Updownarrow$ | $\uparrow$ |
Beta-2-microglobulin | -0.125 | B2M | $\Updownarrow$ | $\downarrow$ |
Neurexin-2 | 0.116 | NRXN2 | $\Updownarrow$ |  |
Apolipoprotein D | -0.095 | APOD | $\Updownarrow$ | $\uparrow$ |
Neuronal pentraxin-1 | 0.093 | NPTX1 | $\downarrow$ |  |
Kallikrein-6 | -0.083 | KLK6 | $\uparrow$ | $\uparrow$ | $\uparrow$
Cystatin-C | -0.066 | CST3 | $\Updownarrow$ | $\Updownarrow$ | $\uparrow$
PC5 |  |  |  |  |
Complement C4-A | -0.180 | C4A | $\uparrow$ |  |
Ectonucleotide pyrophosphatase/phosphodiesterase family member 2 | -0.144 | ENPP2 | $\uparrow$ | $\downarrow$ |
Superoxide dismutase [Cu-Zn] | 0.129 | SOD1 | $\downarrow$ | $\uparrow$ | $\downarrow$
Complement factor B | 0.110 | CFB | $\uparrow$ |  |
Glial fibrillary acidic protein | -0.106 | GFAP | $\uparrow$ |  |
Chromogranin-A | 0.105 | CHGA | $\Updownarrow$ | $\uparrow$ | $\uparrow$
Mimecan | -0.094 | OGN | $\Updownarrow$ |  |
Neurosecretory protein VGF | 0.083 | VGF | $\downarrow$ |  | $\downarrow$
Alpha-1B-glycoprotein | 0.075 | A1BG | $\Updownarrow$ |  |

$\uparrow$: consistently upregulated in MCI/AD or positively correlated;
$\downarrow$: consistently downregulated in MCI/AD or negatively correlated;
$\Updownarrow$: inconsistent reports.

Proteins related to amyloid pathology. Among the top-loaded proteins, SPON1,
SORCS1, PTGDS, CST3, NPTX2, VGF, and CHGA have been found to be related to
amyloid-$\beta$ pathology in AD. The accumulation of amyloid-$\beta$ is
generally considered a hallmark of AD, which is derived from the amyloid
precursor protein (APP) through sequential cleavages by beta-site amyloid
precursor protein cleaving enzyme 1 (BACE1) and $\gamma$-secretase (Vassar et
al.,, 1999). Blocking BACE1 can potentially reduce the abundance in
amyloid-$\beta$, however, this may prohibit the other functions of BACE1 in
psychological activities. For SPON1, using an in vivo AD mouse model, it was
found that, by injecting SPON1, the amount of amyloid-$\beta$ was
significantly reduced, and subsequently, the ameliorated cognitive dysfunction
and memory impairment were improved, suggesting SPON1 to be a potential AD
therapy target (Park et al.,, 2020). Interacting with APOE, human SPON1
suppresses amyloid-$\beta$ level through the APP transgene, and has an impact
on working memory performance through the activation of the triangular part of
the right inferior frontal gyrus (Liu et al.,, 2018). For NPTX1 and NPTX2,
both belong to the family of long neuronal pentraxins. Together with NPTXR,
they bind AMPA type glutamate receptors and contribute to multiple forms of
developmental and adult synaptic plasticity. Using an AD mouse model,
reduction in NPTX2 together with amyloidosis was found to induce a synergistic
reduction of inhibitory circuit function. In AD subjects, the level of NTPX2
was found to be related to hippocampal volume, as well as cognitive decline
(Xiao et al.,, 2017). For CST3, cysteine proteases, including cathepsin B
(CatB), is a recently discovered amyloid-$\beta$-degrading enzyme. Using a
mouse model, CST3 was discovered to be a key inhibitor of CatB-induced
amyloid-$\beta$ degradation in vivo. Genetic ablation of CST3 significantly
reduced soluble amyloid-$\beta$ levels, and attenuated associated cognitive
deficits and behavioral abnormalities, and restored synaptic plasticity in
hippocampus (Sun et al.,, 2008). For VGF, through a mouse model, over-
expression of neuropeptides precursor VGF was found to partially rescue
amyloid-$\beta$-mediated memory impairment and neuropathology, suggesting a
possible causal role of VGF in protecting against AD pathogenesis and
progression (Beckmann et al.,, 2020). For SORCS1, through a meta-analysis of
16 SORCS1-single nucleotide polymorphisms (SNPs) in six independent datasets,
it was found that over-expression of SORCS1 can reduce $\gamma$-secretase
activity and amyloid-$\beta$ levels, and the suppression of SORCS1 can
increase $\gamma$-secretase processing of APP and the levels of
amyloid-$\beta$ (Reitz et al.,, 2011). For PTGDS, it is one of the most
abundant proteins in the CSF, which binds and transports small lipophilic
molecules such as amyloid-$\beta$, and thus has been considered as the
endogenous amyloid-$\beta$ chaperone (Kanekiyo et al.,, 2007), and is believed
to play an important role in AD development. For CHGA, compared to the normal
controls, the level of CHGA was significantly higher in the CSF of patients
with MCI, especially with MCI progressing to AD (Duits et al.,, 2018). CHGA is
the major soluble protein in catecholamine storage vesicles, abnormalities of
which may play a central role in memory deficits in AD. Elevation of CHGA was
observed in AD brains, and was believed to play a role in amyloid-$\beta$
pathology (O’Connor et al.,, 1993; Mattsson et al.,, 2013). It has also been
found that CHGA is negatively associated with hippocampal and entorhinal
volume (Khan et al.,, 2015).

Proteins related to tau pathology. For IGFBP2, it is an abundant cerebral
insulin-like growth factor signaling protein associated with the AD
biomarkers. In both AD mouse models and AD patients, IGFBP2 was observed to be
associated with CSF tau levels and brain atrophy in non-hippocampal regions,
suggesting that it is relevant in neurodegeneration through tau pathology
(Bonham et al.,, 2018).

Proteins related to both amyloid and tau pathology. There was evidence showing
that proteins KLK6 and SOD1 were relevant in both amyloid and tau pathology.
For SOD1, using an APP-overexpressing mouse model, SOD1 deficiency was found
to accelerate amyloid-$\beta$ oligomerization, induce tau phosphorylation and
lower levels of synaptophysin, and consequently memory impairment (Murakami et
al.,, 2011). Kallikrein-related peptidases (KLKs) represent the largest family
of secreted serine proteases. Human KLK6 is the most abundant KLKs in the
spinal cord, brain stem, cerebral cortex including the hippocampus and
thalamus. It has been found that KLK6 cleaves APP and mediates cleavage of
laminin and collagen, which has implications for APP processing and
amyloid-$\beta$ mediated neurotoxicity (Small et al.,, 1993; Angelo et al.,,
2006). In AD patients, the level of KLK6 in CSF is significantly elevated and
is associated with levels of CSF tau suggesting a potential marker of tau
pathology (Goldhardt et al.,, 2019).

Other AD related proteins. NRXN2 is another protein marker that was found to
be up-regulated among MCI patients, especially with MCI progression to AD
(Duits et al.,, 2018). APOD was found to be elevated in the prefrontal cortex
associated with cognitive decline (Thomas et al.,, 2003). GFAP
immunohistochemistry is a marker to assess the oxidative stress and glial cell
activation expressed in astrocytes. Focusing on the human entorhinal cortex
and hippocampus, the GFAP expression was observed in the hippocampus of AD
patients (Hol et al.,, 2003). B2M is a component of major histocompatibility
complex class 1 molecules. Increased soluble B2M has been discovered in the
CSF of patients with AD, and was associated with cognitive decline (Carrette
et al.,, 2003). Using mouse models, elevated B2M was observed in the
hippocampus of aged mice. Injecting exogenous B2M locally in the hippocampus,
impaired hippocampal-dependent cognitive function and neurogenesis were
observed in young mice. The findings suggest that the accumulation of B2M
increases the risk of age-related cognitive dysfunction and neurogenesis
impairment (Smith et al.,, 2015).

Proteins related to brain structure/atrophy. NEGR1 is a member of the
immunoglobulin superfamily of cell adhesion molecules, and is involved in
cortical layering. Using a NEGR1-targeted mouse model, brain morphological
analysis revealed NEGR1-related neuroanatomical abnormalities, including
enlargement of ventricles and decrease in the volume of the whole brain,
corpus callosum, globus pallidus, and hippocampus (Singh et al.,, 2019). CST3
was discovered to be related to a higher hippocampal atrophy rate (Paterson et
al.,, 2014), and atrophy in the entorhinal cortex (Mattsson et al.,, 2014).
APOD and NPTX2 were found to be related to medial temporal lobe atrophy
(Mattsson et al.,, 2014; Swanson et al.,, 2016).

###  4.3 Brain regions

While Table 1 lists the brain regions with nonzero path effects induced by
PC1, PC4, and PC5, Figure 4 visualizes those regions on a template brain. The
identified brain regions include the hippocampus, the entorhinal cortex,
cortical regions on the temporal, parietal and frontal lobes, the lateral
ventricles, and the cerebellum. Brain structural atrophy occurs early in the
medial temporal lobe, including the hippocampus and entorhinal cortex, then
extends soon after to the rest of the cortical areas, usually following a
temporal, parietal, frontal trajectory, whereas the motor areas are affected
toward late stages. (Pini et al.,, 2016). We next discuss those identified
brain regions roughly following this trajectory.

The hippocampus and entorhinal cortex. The hippocampus is a major component of
the human brain located in the medial temporal lobe, and is functionally
involved in response inhibition, episodic memory, and spatial cognition.
Hippocampal atrophy is the best established and validated biomarker across the
entire disease spectrum (Jack Jr et al.,, 2011). The entorhinal cortex also
locates in the medial temporal lobe. It connects the neocortex and the
hippocampus that receives information from the neocortex and projects to the
hippocampus through the perforant pathway (Insausti et al.,, 1995). It has
been consistently reported that, compared to the healthy controls, entorhinal
atrophy was observed in the MCI patients, and more severe atrophy in the AD
patients (Pini et al.,, 2016). The hippocampus and entorhinal cortex, as well
as the anatomically related parahippocampal and perirhinal cortices, are parts
of the medial temporal lobe memory system. Impairments of this system are
responsible for the deficit in episodic memory, and are early hallmark of AD
(Nadel and Hardt,, 2011).

The lateral temporal, parietal, and frontal cortex. The gray matter loss in
the lateral temporal cortex, dorsal parietal, parietal angular and frontal
cortex occurs during the progression from incipient to mild AD. During this
period, cognitive deficits have been observed in both memory and non-memory
domains, including language, visuo-spatial and executive function (Frisoni et
al.,, 2009). Moreover, a higher amount of tau deposition has been observed in
the middle temporal cortex, fusiform gyrus, and entorhinal cortex (Schultz et
al.,, 2018). The fusiform gyrus is critical in facial recognition. Alterations
of gene expression specific to the fusiform gyrus were discovered in AD
patients (Ma et al.,, 2020). The left middle temporal gyrus is related to the
recognition of known faces and accessing word meaning while reading (Acheson
and Hagoort,, 2013). The precuneus, a hub of the default mode network, has
been found to be related to episodic memories (Sadigh-Eteghad et al.,, 2014).
Atrophy in the entorhinal cortex, fusiform, middle temporal gyrus, precuneus,
and precentral has been noted in AD (Parker et al.,, 2018). The association
between atrophy in the insular cortex and memory deficits in AD has been
reported too (Lin et al.,, 2017).

The lateral ventricles. The ventricles are one of the interests in brain
atrophy research as the volumetric measurement is robust to automatic
segmentation due to the sharp contrast between the CSF in the ventricles and
surrounding tissue in T1-weighted images. Thus, as a complement metric of
hemispheric atrophy rates, enlargement in the lateral ventricles is an
important marker of AD progression (Kruthika et al.,, 2019).

The cerebellum. The cerebellum is involved in cognition and emotion and
communicates with cerebral cortices in a topographically organized manner.
Based on existing evidence of cerebellar modulation of cognition and emotion,
it was hypothesized that there exists cerebellar contribution to the cognitive
and neuropsychiatric deficits in AD. However, more research is required to
validate the hypothesis and to understand cerebrocerebellar interactions in AD
pathology (Jacobs et al.,, 2018).

((a)) PC1

((b)) PC4

((c)) PC5

Figure 4: Brain regions with a nonzero mediation effect in (a) PC1, (b) PC4,
and (c) PC5.

##  5 Simulation Study

We complement our data analysis with some additional simulation studies to
further examine the empirical performance of the proposed method.

We generate $\mathbf{X}_{i}\in\mathbb{R}^{r}$ ($i=1,\dots,n$) from a
multivariate normal distribution with mean zero and a covariance matrix whose
eigenvalues exponentially decay. After applying PCA, we obtain
$\tilde{\mathbf{X}}\in\mathbb{R}^{n\times q}$, where $q$ is chosen such that
the top $q$ PCs account for over $80\%$ of total data variation. We then
generate $\mathbf{M}$ and $\mathbf{Y}$ following model (1) given
$\tilde{\mathbf{X}}$. We set $5\%$ of the path effects to be nonzero. We
consider two sets of data dimension, $r=100,p=100$, and $r=350,p=150$, the
latter of which has a similar data dimension as in the ADNI dataset. We also
consider three sample sizes, $n=100,500,1000$.

Table 4 presents the estimated total indirect effects and the indirect effects
of the top six PCs, and Table 5 presents the estimated number of PCs and the
sensitivity and specificity of the identified nonzero path effects. Among all
cases, the estimated number of PCs is 6, which agrees with the truth. From the
tables, we observe that the proposed method achieves a competitive
performance, and the performance improves, with a lower estimation error and a
higher selection accuracy, as the sample size increases.

Table 4: The estimation bias and mean squared error (MSE) of estimating the total indirect effect and indirect effect of top PCs. |  |  |  | $n=100$ |  | $n=500$ |  | $n=1000$
---|---|---|---|---|---|---|---|---
$r$ | $p$ |  | Truth | Bias | MSE |  | Bias | MSE |  | Bias | MSE
|  | Total | -20 | 9.030 | 128.080 |  | -1.827 | 16.593 |  | -0.921 | 13.508
|  | PC1 | -12 | 7.172 | 61.226 |  | 0.033 | 1.843 |  | 0.035 | 2.883
|  | PC2 | 0 | -0.883 | 15.244 |  | -1.410 | 9.968 |  | -0.856 | 4.364
|  | PC3 | -8 | 3.205 | 20.100 |  | -0.168 | 2.531 |  | -0.032 | 2.248
|  | PC4 | 0 | -0.477 | 11.943 |  | -0.272 | 4.242 |  | -0.067 | 1.963
|  | PC5 | 0 | -0.049 | 5.168 |  | -0.078 | 1.328 |  | 0.002 | 0.611
$100$ | $100$ | PC6 | 0 | 0.085 | 2.809 |  | 0.054 | 0.794 |  | 0.019 | 0.370
|  | Total | 8 | -6.317 | 80.520 |  | -1.164 | 37.593 |  | -1.284 | 19.271
|  | PC1 | -8 | 6.461 | 54.644 |  | 0.723 | 13.823 |  | -0.511 | 3.616
|  | PC2 | 12 | -9.528 | 102.378 |  | -0.933 | 20.122 |  | 0.816 | 3.045
|  | PC3 | 4 | -3.446 | 27.204 |  | -1.398 | 10.500 |  | -1.740 | 7.489
|  | PC4 | 0 | -0.093 | 8.015 |  | 0.147 | 2.595 |  | -0.033 | 1.403
|  | PC5 | 0 | 0.153 | 5.445 |  | 0.060 | 1.615 |  | 0.066 | 0.781
$350$ | $150$ | PC6 | 0 | 0.137 | 6.489 |  | 0.238 | 1.054 |  | 0.117 | 0.750
Table 5: The estimated number of PC (and the standard error, SE) in the PCA step and sensitivity and specificity of identifying paths with a nonzero path effect. $r$ | $p$ | $n$ | # PC (SE) | Sensitivity | Specificity
---|---|---|---|---|---
|  | $100$ | 6.03 (0.17) | 0.84 | 0.53
|  | $500$ | 5.21 (0.41) | 1.00 | 0.89
$100$ | $100$ | $1000$ | 6.26 (0.44) | 1.00 | 0.91
|  | $100$ | 5.99 (0.10) | 0.80 | 0.57
|  | $500$ | 6.00 (0.00) | 1.00 | 0.89
$350$ | $150$ | $1000$ | 6.00 (0.00) | 1.00 | 0.91

## Acknowledgments

Data collection and sharing for this project was funded by the Alzheimer’s
Disease Neuroimaging Initiative (ADNI) (National Institutes of Health Grant
U01 AG024904) and DOD ADNI (Department of Defense award number
W81XWH-12-2-0012). ADNI is funded by the National Institute on Aging, the
National Institute of Biomedical Imaging and Bioengineering, and through
generous contributions from the following: AbbVie, Alzheimer’s Association;
Alzheimer’s Drug Discovery Foundation; Araclon Biotech; BioClinica, Inc.;
Biogen; Bristol-Myers Squibb Company; CereSpir, Inc.; Cogstate; Eisai Inc.;
Elan Pharmaceuticals, Inc.; Eli Lilly and Company; EuroImmun; F. Hoffmann-La
Roche Ltd and its affiliated company Genentech, Inc.; Fujirebio; GE
Healthcare; IXICO Ltd.; Janssen Alzheimer Immunotherapy Research &
Development, LLC.; Johnson & Johnson Pharmaceutical Research & Development
LLC.; Lumosity; Lundbeck; Merck & Co., Inc.; Meso Scale Diagnostics, LLC.;
NeuroRx Research; Neurotrack Technologies; Novartis Pharmaceuticals
Corporation; Pfizer Inc.; Piramal Imaging; Servier; Takeda Pharmaceutical
Company; and Transition Therapeutics. The Canadian Institutes of Health
Research is providing funds to support ADNI clinical sites in Canada. Private
sector contributions are facilitated by the Foundation for the National
Institutes of Health (www.fnih.org). The grantee organization is the Northern
California Institute for Research and Education, and the study is coordinated
by the Alzheimer’s Therapeutic Research Institute at the University of
Southern California. ADNI data are disseminated by the Laboratory for Neuro
Imaging at the University of Southern California.

## References

  * Acheson and Hagoort, (2013) Acheson, D. J. and Hagoort, P. (2013).  Stimulating the brain’s language network: syntactic ambiguity resolution after TMS to the inferior frontal gyrus and middle temporal gyrus.  Journal of Cognitive Neuroscience, 25(10):1664–1677.
  * Alzheimer’s Association, (2020) Alzheimer’s Association (2020).  2020 Alzheimer’s disease facts and figures.  Alzheimer’s & Dementia, 16(3):391–460.
  * Angelo et al., (2006) Angelo, P. F., Lima, A. R., Alves, F. M., Blaber, S. I., Scarisbrick, I. A., Blaber, M., Juliano, L., and Juliano, M. A. (2006).  Substrate specificity of human Kallikrein 6 salt and glycosaminoglycan activation effects.  Journal of Biological Chemistry, 281(6):3116–3126.
  * Aung et al., (2020) Aung, M. T., Song, Y., Ferguson, K. K., Cantonwine, D. E., Zeng, L., McElrath, T. F., Pennathur, S., Meeker, J. D., and Mukherjee, B. (2020).  Application of an analytical framework for multivariate mediation analysis of environmental data.  Nature Communications, 11(1):1–13.
  * 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–1182.
  * Beckmann et al., (2020) Beckmann, N. D., Lin, W.-J., Wang, M., Cohain, A. T., Charney, A. W., Wang, P., Ma, W., Wang, Y.-C., Jiang, C., Audrain, M., et al. (2020).  Multiscale causal networks identify VGF as a key regulator of Alzheimer’s disease.  Nature Communications, 11(1):1–19.
  * Bonham et al., (2018) Bonham, L. W., Geier, E. G., Steele, N. Z., Holland, D., Miller, B. L., Dale, A. M., Desikan, R. S., Yokoyama, J. S., Initiative, A. D. N., et al. (2018).  Insulin-like growth factor binding protein 2 is associated with biomarkers of Alzheimer’s disease pathology and shows differential expression in transgenic mice.  Frontiers in Neuroscience, 12:476.
  * Cai et al., (2019) Cai, Q., Wang, H., Li, Z., and Liu, X. (2019).  A survey on multimodal data-driven smart healthcare systems: Approaches and applications.  IEEE Access, 7:133583–133599.
  * Carrette et al., (2003) Carrette, O., Demalte, I., Scherl, A., Yalkinoglu, O., Corthals, G., Burkhard, P., Hochstrasser, D. F., and Sanchez, J.-C. (2003).  A panel of cerebrospinal fluid potential biomarkers for the diagnosis of Alzheimer’s disease.  PROTEOMICS: International Edition, 3(8):1486–1494.
  * 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.
  * Doshi et al., (2016) Doshi, J., Erus, G., Ou, Y., Resnick, S. M., Gur, R. C., Gur, R. E., Satterthwaite, T. D., Furth, S., Davatzikos, C., and Alzheimer’s Neuroimaging Initiative (2016).  MUSE: MUlti-atlas region Segmentation utilizing Ensembles of registration algorithms and parameters, and locally optimal atlas selection.  Neuroimage, 127:186–195.
  * Duits et al., (2018) Duits, F. H., Brinkmalm, G., Teunissen, C. E., Brinkmalm, A., Scheltens, P., Van der Flier, W. M., Zetterberg, H., and Blennow, K. (2018).  Synaptic proteins in CSF as potential novel biomarkers for prognosis in prodromal Alzheimer’s disease.  Alzheimer’s Research & Therapy, 10(1):1–9.
  * Frisoni et al., (2009) Frisoni, G. B., Prestia, A., Rasser, P. E., Bonetti, M., and Thompson, P. M. (2009).  In vivo mapping of incremental cortical atrophy from incipient to overt Alzheimer’s disease.  Journal of Neurology, 256(6):916.
  * Goldhardt et al., (2019) Goldhardt, O., Warnhoff, I., Yakushev, I., Begcevic, I., Förstl, H., Magdolen, V., Soosaipillai, A., Diamandis, E., Alexopoulos, P., and Grimmer, T. (2019).  Kallikrein-related peptidases 6 and 10 are elevated in cerebrospinal fluid of patients with Alzheimer’s disease and associated with CSF-TAU and FDG-PET.  Translational Neurodegeneration, 8(1):1–13.
  * Higgins et al., (2018) Higgins, I. A., Kundu, S., and Guo, Y. (2018).  Integrative Bayesian analysis of brain functional networks incorporating anatomical knowledge.  Neuroimage, 181:263–278.
  * Hol et al., (2003) Hol, E., Roelofs, R., Moraal, E., Sonnemans, M., Sluijs, J., Proper, E., De Graan, P., Fischer, D., and Van Leeuwen, F. (2003).  Neuronal expression of GFAP in patients with Alzheimer pathology and identification of novel GFAP splice forms.  Molecular Psychiatry, 8(9):786–796.
  * Insausti et al., (1995) Insausti, R., Tunon, T., Sobreviela, T., Insausti, A., and Gonzalo, L. (1995).  The human entorhinal cortex: a cytoarchitectonic analysis.  Journal of Comparative Neurology, 355(2):171–198.
  * Jack Jr et al., (2011) Jack Jr, C. R., Barkhof, F., Bernstein, M. A., Cantillon, M., Cole, P. E., DeCarli, C., Dubois, B., Duchesne, S., Fox, N. C., Frisoni, G. B., et al. (2011).  Steps to standardization and validation of hippocampal volumetry as a biomarker in clinical trials and diagnostic criterion for Alzheimer’s disease.  Alzheimer’s & Dementia, 7(4):474–485.
  * Jacobs et al., (2018) Jacobs, H. I., Hopkins, D. A., Mayrhofer, H. C., Bruner, E., van Leeuwen, F. W., Raaijmakers, W., and Schmahmann, J. D. (2018).  The cerebellum in Alzheimer’s disease: evaluating its role in cognitive decline.  Brain, 141(1):37–47.
  * Jagust, (2018) Jagust, W. (2018).  Imaging the evolution and pathophysiology of alzheimer disease.  Nature reviews. Neuroscience, pages 687–700.
  * Kanekiyo et al., (2007) Kanekiyo, T., Ban, T., Aritake, K., Huang, Z.-L., Qu, W.-M., Okazaki, I., Mohri, I., Murayama, S., Ozono, K., Taniike, M., et al. (2007).  Lipocalin-type prostaglandin D synthase/$\beta$-trace is a major amyloid $\beta$-chaperone in human cerebrospinal fluid.  Proceedings of the National Academy of Sciences, 104(15):6412–6417.
  * Khan et al., (2015) Khan, W., Aguilar, C., Kiddle, S. J., Doyle, O., Thambisetty, M., Muehlboeck, S., Sattlecker, M., Newhouse, S., Lovestone, S., Dobson, R., et al. (2015).  A subset of cerebrospinal fluid proteins from a multi-analyte panel associated with brain atrophy, disease classification and prediction in Alzheimer’s disease.  PloS one, 10(8).
  * Kruthika et al., (2019) Kruthika, K., Maheshappa, H., Initiative, A. D. N., et al. (2019).  Multistage classifier-based approach for Alzheimer’s disease prediction and retrieval.  Informatics in Medicine Unlocked, 14:34–42.
  * Lin et al., (2017) Lin, F., Ren, P., Lo, R. Y., Chapman, B. P., Jacobs, A., Baran, T. M., Porsteinsson, A. P., and Foxe, J. J. (2017).  Insula and inferior frontal gyrus’ activities protect memory performance against Alzheimer’s disease pathology in old age.  Journal of Alzheimer’s Disease, 55(2):669–678.
  * Liu et al., (2015) Liu, S., Cai, W., Liu, S., Zhang, F., Fulham, M., Feng, D., Pujol, S., and Kikinis, R. (2015).  Multimodal neuroimaging computing: a review of the applications in neuropsychiatric disorders.  Brain informatics, 2(3):167–180.
  * Liu et al., (2018) Liu, Z., Dai, X., Tao, W., Liu, H., Li, H., Yang, C., Zhang, J., Li, X., Chen, Y., Ma, C., et al. (2018).  APOE influences working memory in non-demented elderly through an interaction with SPON1 rs2618516.  Human Brain Mapping, 39(7):2859–2867.
  * Long et al., (2020) Long, J. P., Irajizad, E., Doecke, J. D., Do, K.-A., and Ha, M. J. (2020).  A framework for mediation analysis with multiple exposures, multivariate mediators, and non-linear response models.  arXiv preprint arXiv:2011.06061.
  * Ma et al., (2020) Ma, D., Fetahu, I. S., Wang, M., Fang, R., Li, J., Liu, H., Gramyk, T., Iwanicki, I., Gu, S., Xu, W., et al. (2020).  The fusiform gyrus exhibits an epigenetic signature for Alzheimer’s disease.  Clinical Epigenetics, 12(1):1–16.
  * Mattsson et al., (2014) Mattsson, N., Insel, P., Nosheny, R., Trojanowski, J. Q., Shaw, L. M., Jack Jr, C. R., Tosun, D., Weiner, M., Initiative, A. D. N., et al. (2014).  Effects of cerebrospinal fluid proteins on brain atrophy rates in cognitively healthy older adults.  Neurobiology of Aging, 35(3):614–622.
  * Mattsson et al., (2013) Mattsson, N., Insel, P., Nosheny, R., Zetterberg, H., Trojanowski, J., Shaw, L., Tosun, D., and Weiner, M. (2013).  CSF protein biomarkers predicting longitudinal reduction of CSF $\beta$-amyloid42 in cognitively healthy elders.  Translational psychiatry, 3(8):e293–e293.
  * Mormino et al., (2009) Mormino, E., Kluth, J., Madison, C., Rabinovici, G., Baker, S., Miller, B., Koeppe, R., Mathis, C., Weiner, M., Jagust, W., et al. (2009).  Episodic memory loss is related to hippocampal-mediated $\beta$-amyloid deposition in elderly subjects.  Brain, 132(5):1310–1323.
  * Murakami et al., (2011) Murakami, K., Murata, N., Noda, Y., Tahara, S., Kaneko, T., Kinoshita, N., Hatsuta, H., Murayama, S., Barnham, K. J., Irie, K., et al. (2011).  SOD1 (copper/zinc superoxide dismutase) deficiency drives amyloid $\beta$ protein oligomerization and memory loss in mouse model of Alzheimer disease.  Journal of Biological Chemistry, 286(52):44557–44568.
  * Nadel and Hardt, (2011) Nadel, L. and Hardt, O. (2011).  Update on memory systems and processes.  Neuropsychopharmacology, 36(1):251–273.
  * Nathoo et al., (2019) Nathoo, F. S., Kong, L., Zhu, H., and Alzheimer’s Disease Neuroimaging Initiative (2019).  A review of statistical methods in imaging genetics.  Canadian Journal of Statistics, 47(1):108–131.
  * O’Connor et al., (1993) O’Connor, D. T., Kailasam, M. T., and Thal, L. J. (1993).  Cerebrospinal fluid chromogranin A is unchanged in Alzheimer dementia.  Neurobiology of Aging, 14(3):267–269.
  * Park et al., (2020) Park, S. Y., Kang, J. Y., Lee, T., Nam, D., Jeon, C.-J., and Kim, J. B. (2020).  SPON1 can reduce amyloid beta and reverse cognitive impairment and memory dysfunction in Alzheimer’s disease mouse model.  Cells, 9(5):1275.
  * Parker et al., (2018) Parker, T. D., Slattery, C. F., Zhang, J., Nicholas, J. M., Paterson, R. W., Foulkes, A. J., Malone, I. B., Thomas, D. L., Modat, M., Cash, D. M., et al. (2018).  Cortical microstructure in young onset Alzheimer’s disease using neurite orientation dispersion and density imaging.  Human Brain Mapping, 39(7):3005–3017.
  * Paterson et al., (2014) Paterson, R., Bartlett, J., Blennow, K., Fox, N., Shaw, L., Trojanowski, J., Zetterberg, H., Schott, J., Initiative, A. D. N., et al. (2014).  Cerebrospinal fluid markers including trefoil factor 3 are associated with neurodegeneration in amyloid-positive individuals.  Translational Psychiatry, 4(7):e419.
  * Pini et al., (2016) Pini, L., Pievani, M., Bocchetta, M., Altomare, D., Bosco, P., Cavedo, E., Galluzzi, S., Marizzoni, M., and Frisoni, G. B. (2016).  Brain atrophy in Alzheimer’s disease and aging.  Ageing Research Reviews, 30:25–48.
  * Reitz et al., (2011) Reitz, C., Tokuhiro, S., Clark, L. N., Conrad, C., Vonsattel, J.-P., Hazrati, L.-N., Palotás, A., Lantigua, R., Medrano, M., Z. Jiménez-Velázquez, I., et al. (2011).  SORCS1 alters amyloid precursor protein processing and variants may increase Alzheimer’s disease risk.  Annals of Neurology, 69(1):47–64.
  * Richardson et al., (2016) Richardson, S., Tseng, G. C., and Sun, W. (2016).  Statistical methods in integrative genomics.  Annual review of statistics and its application, pages 181–209.
  * Rosenbaum, (2002) Rosenbaum, P. R. (2002).  Covariance adjustment in randomized experiments and observational studies.  Statistical Science, 17(3):286–327.
  * Sadigh-Eteghad et al., (2014) Sadigh-Eteghad, S., Majdi, A., Farhoudi, M., Talebi, M., and Mahmoudi, J. (2014).  Different patterns of brain activation in normal aging and Alzheimer’s disease from cognitional sight: meta analysis using activation likelihood estimation.  Journal of the Neurological Sciences, 343(1-2):159–166.
  * Schultz et al., (2018) Schultz, S. A., Gordon, B. A., Mishra, S., Su, Y., Perrin, R. J., Cairns, N. J., Morris, J. C., Ances, B. M., and Benzinger, T. L. (2018).  Widespread distribution of tauopathy in preclinical Alzheimer’s disease.  Neurobiology of Aging, 72:177–185.
  * Simon et al., (2013) Simon, N., Friedman, J., Hastie, T., and Tibshirani, R. (2013).  A sparse-group lasso.  Journal of Computational and Graphical Statistics, 22(2):231–245.
  * Singh et al., (2019) Singh, K., Jayaram, M., Kaare, M., Leidmaa, E., Jagomäe, T., Heinla, I., Hickey, M. A., Kaasik, A., Schäfer, M. K., Innos, J., et al. (2019).  Neural cell adhesion molecule Negr1 deficiency in mouse results in structural brain endophenotypes and behavioral deviations related to psychiatric disorders.  Scientific Reports, 9(1):1–15.
  * Small et al., (1993) Small, D. H., Nurcombe, V., Clarris, H., Beyreuther, K., and Masters, C. L. (1993).  The role of extracellular matrix in the processing of the amyloid protein precursor of Alzheimer’s disease.  Annals of the New York Academy of Sciences, 695(1):169–174.
  * Smith et al., (2015) Smith, L. K., He, Y., Park, J.-S., Bieri, G., Snethlage, C. E., Lin, K., Gontier, G., Wabl, R., Plambeck, K. E., Udeochu, J., et al. (2015).  $\beta$2-microglobulin is a systemic pro-aging factor that impairs cognitive function and neurogenesis.  Nature Medicine, 21(8):932–937.
  * Song et al., (2018) Song, Y., Zhou, X., Zhang, M., Zhao, W., Liu, Y., Kardia, S. L., Roux, A. V. D., Needham, B. L., Smith, J. A., and Mukherjee, B. (2018).  Bayesian shrinkage estimation of high dimensional causal mediation effects in omics studies.  Biometrics.
  * Sun et al., (2008) Sun, B., Zhou, Y., Halabisky, B., Lo, I., Cho, S.-H., Mueller-Steiner, S., Devidze, N., Wang, X., Grubb, A., and Gan, L. (2008).  Cystatin C-cathepsin B axis regulates amyloid beta levels and associated neuronal deficits in an animal model of Alzheimer’s disease.  Neuron, 60(2):247–257.
  * Swanson et al., (2016) Swanson, A., Willette, A., Initiative, A. D. N., et al. (2016).  Neuronal pentraxin 2 predicts medial temporal atrophy and memory decline across the Alzheimer’s disease spectrum.  Brain, Behavior, and Immunity, 58:201–208.
  * Thomas et al., (2003) Thomas, E. A., Laws, S. M., Sutcliffe, J. G., Harper, C., Dean, B., McClean, C., Masters, C., Lautenschlager, N., Gandy, S. E., and Martins, R. N. (2003).  Apolipoprotein D levels are elevated in prefrontal cortex of subjects with Alzheimer’s disease: no relation to apolipoprotein E expression or genotype.  Biological Psychiatry, 54(2):136–141.
  * VanderWeele, (2016) VanderWeele, T. J. (2016).  Mediation analysis: A practitioner’s guide.  Annual Review of Public Health, 37(1):17–32.
  * Vassar et al., (1999) Vassar, R., Bennett, B. D., Babu-Khan, S., Kahn, S., Mendiaz, E. A., Denis, P., Teplow, D. B., Ross, S., Amarante, P., Loeloff, R., et al. (1999).  $\beta$-secretase cleavage of Alzheimer’s amyloid precursor protein by the transmembrane aspartic protease BACE.  Science, 286(5440):735–741.
  * Wesenhagen et al., (2020) Wesenhagen, K. E., Teunissen, C. E., Visser, P. J., and Tijms, B. M. (2020).  Cerebrospinal fluid proteomics and biological heterogeneity in Alzheimer’s disease: A literature review.  Critical Reviews in Clinical Laboratory Sciences, 57(2):86–98.
  * Xiao et al., (2017) Xiao, M.-F., Xu, D., Craig, M. T., Pelkey, K. A., Chien, C.-C., Shi, Y., Zhang, J., Resnick, S., Pletnikova, O., Salmon, D., et al. (2017).  NPTX2 and cognitive dysfunction in Alzheimer’s disease.  Elife, 6:e23798.
  * Zhang, (2019) Zhang, Q. (2019).  High dimensional mediation analysis with applications to causal gene identification.  bioRxiv, page 497826.
  * Zhao et al., (2020) Zhao, Y., Li, L., and Caffo, B. S. (2020).  Multimodal neuroimaging data integration and pathway analysis.  Biometrics.
  * Zhao and Luo, (2016) Zhao, Y. and Luo, X. (2016).  Pathway lasso: Estimate and select sparse mediation pathways with high dimensional mediators.  arXiv preprint arXiv:1603.07749.
  * Zou and Hastie, (2005) Zou, H. and Hastie, T. (2005).  Regularization and variable selection via the elastic net.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67(2):301–320.