Mapping the Genetic-Imaging-Clinical Pathway with Applications to Alzheimer’s Disease

Alzheimer’s disease (AD) is an irreversible brain disorder that slowly destroys memory and thinking skills. According to World Alzheimer Reports (Gaugler et al., 2019), there are around $55$ million people worldwide living with Alzheimer’s disease and related dementia. The total global cost of Alzheimer’s disease and related dementia was estimated to be a trillion US dollars, equivalent to $1.1\%$ of global gross domestic product. Alzheimer’s patients often suffer from behavioral deficits including memory loss and difficulty of thinking, reasoning and decision making.

In the current model of AD pathogenesis, it is well established that deposition of amyloid plaques is an early event that, in conjunction with tau pathology, causes neuronal damage. Scientists have identified risk genes that may cause the abnormal aggregation and deposition of the amyloid plaques (e.g. Morishima-Kawashima and Ihara, 2002). The neuronal damage typically starts from the hippocampus and results in the first clinical manifestations of the disease in the form of episodic memory deficits (Weiner et al., 2013). Specifically, Jack Jr et al. (2010) presented a hypothetical model for biomarker dynamics in AD pathogenesis, which has been empirically and collectively supported by many works in the literature. The model begins with the abnormal deposition of $\beta$ amyloid (A$\beta$) fibrils, as evidenced by a corresponding drop in the levels of soluble A$\beta$-42 in cerebrospinal fluid (CSF) (Aizenstein et al., 2008). After that, neuronal damage begins to occur, as evidenced by increased levels of CSF tau protein (Hesse et al., 2001). Numerous studies have investigated how A$\beta$ and tau impact the hippocampus (e.g. Ferreira and Klein, 2011), known to be fundamentally involved in acquisition, consolidation, and recollection of new episodic memories (Frozza et al., 2018). In particular, as neuronal degeneration progresses, brain atrophy, which starts with hippocampal atrophy (Fox et al., 1996), becomes detectable by magnetic resonance imaging (MRI). Studies from recent years also found other important CSF proteins that may be related to hippocampal atrophy. For instance, the low levels of chromogranin A (CgA) and trefoil factor 3 (TFF3), and high level of cystatin C (CysC) are evidently associated with hippocampal atrophy (Khan et al., 2015; Paterson et al., 2014). Indeed, the impact of protein concentration on behavior can also be through atrophy of other brain regions. For example, there exists potential entorhinal tau pathology on episodic memory decline (Maass et al., 2018). As sufficient brain atrophy accumulates, it results in cognitive symptoms and impairment. This process of AD pathogenesis is summarized by the flow chart in Figure 1. Note that it is still debatable how A$\beta$ and tau interact with each other as mentioned by Jack Jr et al. (2013); Majdi et al. (2020); however, it is evident that A$\beta$ may still hit a biomarker detection threshold earlier than tau (Jack Jr et al., 2013). In addition, as noted by Hampel et al. (2018), it is likely that highly complex interactions exist between A$\beta$ as well as tau, and the cholinergic system. For instance, the association has been found between CSF biomarkers of amyloid and tau pathology in AD (Remnestål et al., 2021). It has also been found that other factors, such as dysregulation and dysfunction of the Wnt signaling pathway, may also contribute to A$\beta$ and tau pathologies (Ferrari et al., 2014). In addition, the M1 and M3 subtypes of muscarinic receptors increase amyloid precursor protein production via the induction of the phospholipase C/protein kinase C pathway and increase BACE expression in AD brains (Nitsch et al., 1992).

The aim of this paper is to map the genetic-imaging-clinical (GIC) pathway for AD, which is the most important part of the hypothetical model of AD pathogenesis in Figure 1. Histological studies have shown that the hippocampus is particularly vulnerable to Alzheimer’s disease pathology and has already been considerably damaged at the first occurrence of clinical symptoms (Braak and Braak, 1998). Therefore, the hippocampus has become a major focus in Alzheimer’s studies (De Leon et al., 1989). Some neuroscientists even conjecture that the association between hippocampal atrophy and behavioral deficits may be causal, because the former destroys the connections that help the neuron communicate and results in a loss of function (Watson, 2019). We are interested in delineating the genetically-regulated hippocampal shape that drives AD related behavioral deficits and disease progression.

To map the GIC pathway, we extract clinical, imaging, and genetic variables from the Alzheimer’s Disease Neuroimaging Initiative (ADNI) study as follows. First, we use the Alzheimer’s Disease Assessment Scale (ADAS) cognitive score to quantify behavioral deficits, for which a higher score indicates more severe behavioral deficits. Second, we characterize the exposure of interest, hippocampal shape, by the hippocampal morphometry surface measure, summarized as two $100\times 150$ matrices corresponding to the left/right hippocampi. Each element of the matrices is a continuous-valued variable, representing the radial distance from the corresponding coordinate on the hipppocampal surface to the medial core of the hippocampus. Compared with the conventional scalar measure of hippocampus shape (Jack et al., 2003), recent studies show that the additional information contained in the hippocampal morphometry surface measure is valuable for Alzheimer’s diagnosis (Thompson et al., 2004). For example, Li et al. (2007) showed that the surface measures of the hippocampus could provide more subtle indexes compared with the volume differences in discriminating between patients with Alzheimer’s and healthy control subjects. In our case, with the 2D matrix radial distance measure, one may investigate how local shapes of hippocampal subfields are associated with the behavioral deficits. Third, the ADNI study measures ultra-high dimensional genetic covariates and other demongraphic covariates at baseline. There are more than $6$ million genetic variants per subject.

The special data structure of the ADNI data application presents new challenges for statistically mapping the GIC pathway. First, unlike conventional statistical analysis that deals with scalar exposure, our exposure of interest is represented by high-dimensional 2D hippocampal imaging measures. Second, the dimension of baseline covariates, which are also potential confounders, is much larger than the sample size. Recently there have been many developments for confounder selection, most of which are in the causal inference literature. Studies show inclusion of the variables only associated with exposure but not directly with the outcome except through the exposure (known as instrumental variables) may result in loss of efficiency in the causal effect estimate (e.g. Schisterman et al., 2009), while inclusion of variables only related to the outcome but not the exposure (known as precision variables) may provide efficiency gains (e.g. Brookhart et al., 2006); see Shortreed and Ertefaie (2017); Richardson et al. (2018); Tang et al. (2021) and references therein for an overview.

When a large number of covariates are available, the primary difficulty for mapping the GIC pathway is how to include all the confounders and precision variables, while excluding all the instrumental variables and irrelevant variables (not related to either outcome or exposure). We develop a novel two-step approach to estimate the conditional association between the high-dimensional 2D hippocampal surface exposure and the Alzheimer’s behaviorial score, while taking into account the ultra-high dimensional baseline covariates. The first step is a fast screening procedure based on both the outcome and exposure models to rule out most of the irrelevant variables. The use of both models in screening is crucial for both computational efficiency and selection accuracy, as we will show in detail in Sections 3.3 and 3.4. The second step is a penalized regression procedure for the outcome generating model to further exclude instrumental and irrelevant variables, and simultaneously estimate the conditional association. Our simulations and ADNI data application demonstrate the effectiveness of the proposed procedure.

Our analysis represents a novel inferential target compared to recent developments in imaging genetics mediation analysis (Bi et al., 2017). Although we consider a similar set of variables and structure among these variables as illustrated later in Figure 2, our goal is to estimate the conditional association of hippocampal shape with behavioral deficits. In contrast, in mediation analysis, researchers are often interested in the effects of genetic factors on behavioral deficits, and how those are mediated through hippocampus. Direct application of methods developed for imaging genetics mediation analysis to our problem may select genetic factors that are confounders affecting both hippocampal shape and behavioral deficits. In comparison, we aim to include precision variables in the adjustment set as they may improve efficiency.

The rest of the article is organized as follows. Section 2 includes a detailed data and problem description. We introduce our models and a two-step variable selection procedure in Section 3. We analyze the ADNI data and estimate the association between hippocampal shape and behavioral deficits conditional on observed clinical and genetic covariates in Section 4. Simulations are conducted in Section 5 to evaluate the finite-sample performance of the proposed method. We finish with a discussion in Section 6. The theoretical properties of our procedure are included in Section 15 in the supplementary material.

Understanding how human brains work and how they connect to human behavior is a central goal in medical studies. In this paper, we are interested in studying whether and how hippocampal shape is associated with behavioral deficits in Alzheimer’s studies. We consider the clinical, genetic, imaging and behavioral measures in the ADNI dataset. The outcome of interest is the Alzheimer’s Disease Assessment Scale cognitive score observed at 24th month after baseline measurements. The Alzheimer’s Disease Assessment Scale cognitive 13 items score (ADAS-13) (Mohs et al., 1997) includes 13 items: word recall task, naming objects and fingers, following commands, constructional praxis, ideational praxis, orientation, word recognition task, remembering test directions, spoken language, comprehension, and word-finding difficulty, delayed word recall and a number cancellation or maze task. A higher ADAS score indicates more severe behavioral deficits.

The exposure of interest is the baseline 2D surface data obtained from the left/right hippocampi. The hippocampus surface data were preprocessed from the raw MRI data, where the detailed preprocessing steps are included in Section 9.1 of the supplementary material. After preprocessing, we obtained left and right hippocampal shape representations as two $100\times 150$ matrices. The imaging measurement at each pixel is an absolute metric, representing the radial distance from the pixel to the medial core of the hippocampus. The unit for the measurement is in millimeters.

In the ADNI data, there are millions of observed covariates that one may need to adjust for, including the whole genome sequencing data from all of the 22 autosomes. We have included detailed genetic preprocessing techniques in Section 9.2 of the supplementary material. After preprocessing, $6,087,205$ bi-allelic markers (including SNPs and indels) of $756$ subjects were retained in the data analysis.

We excluded those subjects with missing hippocampus shape representations, baseline intracranial volume (ICV) information or ADAS-13 score observed at Month 24, after which there are $566$ subjects left. Our aim is to estimate the association between the hippocampal surface exposure and the ADAS-13 score conditional on clinical measures including age, gender and length of education, ICV, diagnosis status, and $6,087,205$ bi-allelic markers.

We use the data obtained from the ADNI study (adni.loni.usc.edu). The data usage acknowledgement is included in Section 8 of the supplement material. As described in Section 2, we include $566$ subjects from the ADNI1 study. The exposure of interest is the baseline 2D hippocampal surface radial distance measures, which can be represented as a $100\times 150$ matrix for each part of the hippocampus. The outcome of interest is the ADAS-13 score observed at Month 24. The average ADAS-13 score is $20.8$ with standard deviation $14.1$. The covariates to adjust for include $6,087,205$ bi-allelic markers as well as clinical covariates, including age, gender and education length, baseline intracranial volume (ICV), and baseline diagnosis status. The average age is $75.5$ years old with standard deviation $6.6$ years, and the average education length is $15.6$ years with standard deviation $2.9$ years. Among all the $566$ subjects, $58.1\%$ were female. The average ICV was $1.28\times 10^{6}$ $\rm{mm}^{3}$ with standard deviation $1.35\times 10^{5}$ $\rm{mm}^{3}$. There were $175$ ($184$ at Month 24) cognitive normal patients, $268$ ($157$ at Month 24) patients with mild cognitive impairment (MCI), and $123$ ($225$ at Month 24) patients with AD at the baseline. Studies have shown that age and gender are the main risk factors for Alzheimer’s disease (Vina and Lloret, 2010; Guerreiro and Bras, 2015) with older people and females more likely to develop Alzheimer’s disease. Multiple studies have also shown that prevalence of dementia is greater among those with low or no education (Zhang et al., 1990). On the other hand, age, gender and length of education have been found to be strongly associated with the hippocampus (Van de Pol et al., 2006; Jack et al., 2000; Noble et al., 2012). Previous studies (Sargolzaei et al., 2015) suggest that the ICV is an important measure that needs to be adjusted for in studies of brain change and AD. In addition, the baseline diagnosis status may help explain the baseline hippocampal shape and the AD status at Month 24. Therefore, we consider age, gender, education length, baseline ICV and baseline diagnosis status as part of the confounders, and adjust for them in our analysis. In addition, we also adjust for population stratification, for which we use the top five principal components of the whole genome data. As both left and right hippocampi have 2D radial distance measures and the two parts of hippocampi have been found to be asymmetric (Pedraza et al., 2004), we apply our method to the left and right hippocampi separately.

We use the default method (Gabriel et al., 2002) of Haploview (Barrett et al., 2005) and PLINK (Purcell et al., 2007) to form linkage disequilibrium (LD) blocks. Previous studies report that about 50 genetic variants are associated with AD; see the review in Sims et al. (2020) for details. This provides support for our assumption that $|\mathcal{M}_{1}|<n$ ($n=566$). On the other hand, a genome-wide association analysis of 19,629 individuals by Zhao et al. (2019) shows that $57$ genetic variants are associated with the left hippocampal volumes and $54$ are associated with the right hippocampal volumes. This provides support for our assumption that $|\mathcal{M}_{2}|<n$. Therefore, we apply our blockwise joint screening procedure on those SNPs on each part of hippocampal outcome $Y_{i}$ and exposure $\bm{Z}_{i}$ marginally. We choose the thresholds $\gamma_{1,n}$, $\gamma_{2,n}$,$\gamma_{3,n}$ and $\gamma_{4,n}$ such that $|\widehat{\mathcal{M}}^{block}|=2\lfloor n/\log(n)\rfloor=178$. In Table 3 of the supplementary material, we list the top $20$ SNPs corresponding to left and right hippocampi, respectively. As suggested by one referee, we plot similar figures as the Manhattan plot for $\widehat{\mathcal{M}}^{*}_{1}$, $\widehat{\mathcal{M}}^{block,*}_{1}$, $\widehat{\mathcal{M}}^{\textbf{}}_{2}$ and $\widehat{\mathcal{M}}^{block}_{2}$ in Figure 7 of the supplementary material, where genomic coordinates are displayed along the x-axis, the y-axis represents the magnitude of $|\widehat{\beta}^{M}_{l}|$, $\widehat{\beta}^{block,M}_{l}$, $\|\widehat{\bm{C}}_{l}^{M}\|_{op}$ and $\widehat{C}_{l}^{block,M}$ and the horizontal dashed line represents the threshold values $\gamma_{1,n}$, $\gamma_{2,n}$, $\gamma_{3,n}$, and $\gamma_{4,n}$, respectively.

From Table 3 and Figure 7, one can see that there are quite a few important SNPs for both hippocampi. For example, the top SNP is the rs429358 from the 19th chromosome. This SNP is a C/T single-nucleotide variant (snv) variation in the APOE gene. It is also one of the two SNPs that define the well-known APOE alleles, the major genetic risk factor for Alzheimer’s disease (Kim et al., 2009). In addition, a great portion of the SNPs in Table 3 has been found to be strongly associated with Alzheimer’s. These include rs10414043 (Du et al., 2018), an A/G snv variation in the APOC1 gene; rs7256200 (Takei et al., 2009), an A/G snv variation in the APOC1 gene; rs73052335 (Zhou et al., 2018), an A/C snv variation in the APOC1 gene; rs769449 (Chung et al., 2014), an A/G snv variation in the APOE gene; rs157594 (Hao et al., 2017), a G/T snv variation; rs56131196 (Gao et al., 2016), an A/G snv variation in the APOC1 gene; rs111789331 (Gao et al., 2016), an A/T snv variation; and rs4420638 (Coon et al., 2007), an A/G snv variation in the APOC1 gene.

Among those SNPs that have been found to be associated with Alzheimer’s, some of them are also directly associated with hippocampi. For example, Zhou et al. (2020) revealed that the SNPs rs10414043, rs73052335 and rs769449 are among the top SNPs that have significant genetic effects on the volumes of both left and right hippocampi. Guo et al. (2019) identified the SNP rs56131196 to be associated with hippocampal shape.

We then perform our second-step estimation procedure for each part of the hippocampi. Here $X_{\widehat{\mathcal{M}}}$ denotes the SNPs selected in the screening step, the population stratification (top five principal components of the whole genome data) and the five clinical measures (age, gender, education length, baseline ICV and baseline diagnosis status), and ${\bm{Z}}$ denotes the left/right hippocampal surface image matrix. To visualize the results, we map the estimates $\widehat{\bm{B}}$ corresponding to each part of the hippocampus onto a representative hippocampal surface and plot it in Figure 3(a). We have also plotted the hippocampal subfield (Apostolova et al., 2006) in Figure 3(b). Here Cornu Ammonis region 1 (CA1), Cornu Ammonis region 2 (CA2) and Cornu Ammonis region 3 (CA3) are a strip of pyramidal neurons within the hippocampus proper. CA1 is the top portion, named as “regio superior of Cajal” (Blackstad et al., 1970), which consists of small pyramidal neurons. Within the lower portion (regio inferior of Cajal), which consists of larger pyramidal neurons, there are a smaller area called CA2 and a larger area called CA3. The cytoarchitecture and connectivity of CA2 and CA3 are different. The subiculum is a pivotal structure of the hippocampal formation, positioned between the entorhinal cortex and the CA1 subfield of the hippocampus proper (for a complete review, see Dudek et al. 2016).

From the plots, we can see that all the $15,000$ entries of $\widehat{\bm{B}}$ corresponding to both hippocampi are negative. This implies that the radial distance of each pixel of both hippocampi is negatively associated with the ADAS-13 score, which depicts the severity of behavioral deficits. The subfields with strongest associations are mostly CA1 and subiculum. Existing literature (Apostolova et al., 2010) has found that as Alzheimer’s disease progresses, it first affects CA1 and subiculum subregions and later CA2 and CA3 subregions. This can partially explain why the shapes of CA1 and subiculum may have stronger associations with ADAS-13 scores compared to CA2 and CA3 subregions.

We examine the effect size of the whole hippocampal shape by evaluating the term $\langle\bm{Z}_{i},\widehat{\bm{B}}\rangle$. Specifically, we calculate the proportion of variance explained by imaging covariates as follows:

Our results show that the shape of the left hippocampi and that of the right one account for 5.83% and 4.71% of the total variations in behavior deficits, respectively. Such effect sizes are quite large compared with those for polygenic risk scores in genetics. In addition, we perform permutation test to test whether the $R^{2}$ statistic is significant. In particular, we randomly permutate the $\{Y_{1},\ldots,Y_{n}\}$, denoted by $\{Y_{i}^{*},\ldots,Y_{n}^{*}\}$, and we then apply our estimation procedure based on $(X_{i},\bm{Z}_{i},Y_{i}^{*})$, obtain $\widehat{\bm{B}}^{*}$, and calculate $(R^{2})^{*}$. We repeat this for 1,000 times and and we obtain the $\{(R_{(k)}^{2})^{*},1\leq k\leq 1000\}$, which mimics the null distribution. Finally, the p-value can be calculated as $\frac{1}{1000}\sum_{k=1}^{1000}1\{(R_{(k)}^{2})^{*}\geq R^{2}\}$. The p-values for both hippocampi are less than $0.001$, suggesting that the $R^{2}$’s are significantly different from zero.

We also conduct sensitivity analysis by varying the relative sizes of $\widehat{\mathcal{M}}_{1}^{*}$ and $\widehat{\mathcal{M}}_{2}$ in the joint screening procedure. The estimates $\widehat{\bm{B}}$s are similar among different choices of $|\widehat{\mathcal{M}}_{1}^{*}|$ and $|\widehat{\mathcal{M}}_{2}|$; see supplementary material Section 11 for details. In addition, we repeated our analysis on the $391$ MCI and AD subjects. We have similar findings that the radial distances of each pixel of both hippocampi are mostly negatively associated with the ADAS-13 score. And the subfields with strongest associations are mostly CA1 and subiculum; see supplementary material Section 12 for details. As suggested by one referee, we have performed the SNP-imaging-outcome mediation analysis proposed by Bi et al. (2017); see Section 13 in the supplementary material for the detailed procedure. There is no evidence for the mediating relationship of SNP-imaging-outcome from our analysis.

In this section, we perform simulation studies to evaluate the finite sample performance of the proposed method. The dimension of covariates is set as $s=5000$, and the exposure is a $64\times 64$ matrix. The $X_{i}\in\mathbb{R}^{s}$ is independently generated from $N(\bm{0},\bm{\Sigma}_{x})$, where $\bm{\Sigma}_{x}=(\sigma_{x,ll^{\prime}})$ has an autoregressive structure such that $\sigma_{x,ll^{\prime}}=\rho_{1}^{|l-l^{\prime}|}$ holds for $1\leq l,l^{\prime}\leq s$ with $\rho_{1}=0.5$. Define ${\bm{B}}$ as a $64\times 64$ image shown in Figure 4(a), and ${\bm{C}}$ a $64\times 64$ image shown in Figure 4(b). For ${\bm{B}}$, the black regions of interest (ROIs) are assigned value $0.0408$ and white ROIs are assigned value $0$. For ${\bm{C}}$, the black ROIs are assigned value $0.0335$ and white ROIs are assigned value $0$. Further we set ${\bm{C}}_{l}=v_{l}*\bm{C}$, where $v_{1}=-1/3$, $v_{2}=-1$, $v_{3}=-3$, $v_{207}=-3$, $v_{208}=-1$, $v_{209}=-1/3$, and $v_{l}=0$ for $4\leq l\leq 206$ and $210\leq l\leq s$. We set the parameters $\beta_{1}=3$, $\beta_{2}=1$, $\beta_{3}=1/3$, $\beta_{104}=3$, $\beta_{105}=1$, $\beta_{106}=1/3$, and $\beta_{l}=0$ for $4\leq l\leq 103$ and $107\leq l\leq s$. In this setting, we have $\mathcal{C}=\{1,2,3\}$, $\mathcal{P}=\{104,105,106\}$, $\mathcal{I}=\{207,208,209\}$ and $\mathcal{S}=\{1\ldots,5000\}\backslash\{1,2,3,104,105,106,207,208,209\}$.

The random error $\mathrm{vec}(\bm{E}_{i})$ is independently generated from $N(\bm{0},\bm{\Sigma}_{e})$, where we set the standard deviations of all elements in $\bm{E}_{i}$ to be $\sigma_{e}=0.2$ and the correlation between $\bm{E}_{i,jk}$ and $\bm{E}_{i,j^{\prime}k^{\prime}}$ to be $\rho_{2}^{|j-j^{\prime}|+|k-k^{\prime}|}$ for $1\leq j,k,j^{\prime},k^{\prime}\leq 64$ with $\rho_{2}=0.5$. The random error $\epsilon_{i}$ is generated independently from $N(0,\sigma^{2})$, where we consider $\sigma=1$ or $0.5$. The $Y_{i}$’s and ${\bm{Z}}_{i}$’s are generated from models (1) and (2). We consider three different sample sizes $n=200,500$ and $1000$.

This paper aims at mapping the complex GIC pathway for AD. The unique features of the hippocampal morphometry surface measure data motivate us to develop a computationally efficient two-step screening and estimation procedure, which can select biomarkers among more than $6$ million observed covariates and estimate the conditional association simultaneously. If there was no unmeasured confounding, then the conditional association we estimate corresponds to the causal effect. This is, however, not the case in the ADNI study because we have unmeasured confounders such as A$\beta$ and tau protein levels. To control for unmeasured confounding and estimate the causal effect, one possible approach is to use the generic variants as potential instrumental variables (e.g. Lin et al., 2015).

There are a number of other important directions for future work. Firstly, the vast majority of AD, known as “polygenic AD”, is influenced by the actions and interactions of multiple genetic variants simultaneously, likely in concert with non-genetic factors, such as environmental exposures and lifestyle choices among many others (Bertram and Tanzi, 2020). Therefore, various types of interaction effects, such as genetic-genetic and imaging-genetic, could be incorporated into the outcome generating model (1). However, this may significantly increase the computation as the dimension of genetic relevant covariates will increase from $6,087,205$ to more than $90$ billion covariates, if we add all the possible imaging-genetics interaction terms. One may consider interaction screening procedures (Hao and Zhang, 2014) as the first-step. Secondly, this study simply removes observations with missingness. Accommodation of missing exposure, confounders and outcome under the proposed model framework is of great practical value and worth further investigation. Thirdly, baseline diagnosis status is an important effect modifier, as the effect of hippocampus shape on behavioral measures can be different across the CN/MCI/AD groups. However, the relatively small sample size in the ADNI study does not allow us to conduct a reliable subgroup analysis. The subgroup analyses are pertinent for further exploration when a larger sample size is available. Fourthly, in the ADNI dataset, there are longitudinal ADAS-13 scores observed at different months as well as other longitudinal behavioral scores obtained from Mini-Mental State Examination and Rey Auditory Verbal Learning Test, which can provide a more comprehensive characterization of the behavioral deficits. Integrating these different scores as a multivariate longitudinal outcome to improve the estimation of the conditional association requires substantial effort for future research. Lastly, one could consider incorporating information from other brain regions. For instance, an entorhinal tau may exist on episodic memory decline through other brain regions, such as the medial temporal lobe (Maass et al., 2018).

Supplementary material available online contains detailed derivations and explanations of the main algorithm, ADNI data usage acknowledgement, image and genetic data preprocessing steps, screening results and sensitivity analyses of the ADNI data application with a subgroup analysis including only MCI and AD patients, detailed procedure and results for the SNP-imaging-outcome mediation analyses, additional simulation results, theoretical properties of the proposed procedure including the main theorems, assumptions needed for our main theorems, and proofs of auxiliary lemmas and main theorems.

The authors thank the editor, associate editor and referees for their constructive comments, which have substantially improved the paper. Yu was partially supported by the Canadian Statistical Sciences Institute (CANSSI) postdoctoral fellowship and the start-up fund of University of Texas at Arlington. Wang and Kong were partially supported by the Natural Science and Engineering Research Council of Canada and the CANSSI Collaborative Research Team Grant. Zhu was partially supported by the NIH-R01-MH116527.

The supplementary file is organized as follows. The detailed description of the main algorithm is included in Section 7. We include the ADNI data usage acknowledgement in Section 8 and image and genetic data preprocessing steps in Section 9. In Section 10, we list the screening results of ADNI data applications. Section 11 examines the sensitivity of the estimate $\widehat{\bm{B}}$ from the ADNI data application by varying the relative sizes of $\widehat{\mathcal{M}}_{1}^{*}$ and $\widehat{\mathcal{M}}_{2}$. Section 12 includes a subgroup analysis including only the mild cognitive impairment (MCI) and Alzheimer’s disease (AD) patients. We include the detailed SNP-imaging-outcome mediation analysis procedure and results in Section 13. In Section 14, we list additional simulation results. The theoretical properties including the main theorems of our procedure are included in Section 15. We state the assumptions needed for the main theorems in Section 16. In Section 17, we include the auxiliary lemmas needed for the theorems and their proofs. We give the detailed proofs of our main theorems in Section 18.

To solve the minimization problem (3) of the main paper, we utilize the Nesterov optimal gradient method (Nesterov, 1998), which has been widely used in solving optimization problems for non-smooth and non-convex objective functions (Beck and Teboulle, 2009; Zhou and Li, 2014).

Before we introduce Nesterov’s gradient algorithm, we first state two propositions on shrinkage thresholding formulas (Beck and Teboulle, 2009; Cai et al., 2010).

The Nestrov’s gradient method utilizes the first-order gradient of the objective function to produce the next iterate based on the current search point. Differed from the standard gradient descent algorithm, the Nesterov’s gradient algorithm uses two previous iterates to generate the next search point by extrapolating, which can dramatically improve the convergence rate. The Nesterov’s gradient algorithm for problem (3) is presented as follows. Denote $l(\bm{\beta},\bm{B})=\frac{1}{2n}\sum_{i=1}^{n}\left(Y_{i}-\langle\bm{\beta},X_{i}\rangle-\langle\bm{Z}_{i},\bm{B}\rangle\right)^{2}$ and $P(\bm{\beta},\bm{B})=P_{1}(\bm{\beta})+P_{2}(\bm{B})$, where $P_{1}(\bm{\beta})=\lambda_{1}\sum_{l}|\beta_{l}|$ and $P_{2}(\bm{B})=\lambda_{2}||\bm{B}||_{*}$. We also define

where $\nabla l(\bm{\beta},\bm{B})=[(\partial_{\bm{\beta}}l)^{\mathrm{\scriptscriptstyle T}},\{\partial_{\mathrm{vec}(\bm{B})}l\}^{\mathrm{\scriptscriptstyle T}}]^{\mathrm{\scriptscriptstyle T}}\in\mathbb{R}^{|\widehat{\mathcal{M}}|+pq}$ denotes the first-order gradient of $l(\bm{\beta},\bm{B})$ with respect to $\left[\bm{\beta}^{\mathrm{\scriptscriptstyle T}},\{\mathrm{vec}(\bm{B})\}^{\mathrm{\scriptscriptstyle T}}\right]^{\mathrm{\scriptscriptstyle T}}\in\mathbb{R}^{|\widehat{\mathcal{M}}|+pq}$. We define