Beyond Scalar Treatment: A Causal Analysis of Hippocampal Atrophy on Behavioral Deficits in Alzheimer’s Studies

Alzheimer’s disease is a progressive form of dementia that results in problems with memory, thinking and behavior. It often starts with abnormal aggregation and deposition of β amyloid and tau, followed by neuronal damage such as atrophy of the hippocampi, and finally leads to behavioral deficits. Despite significant progress in finding biomarkers associated with behavioral deficits, the underlying causal mechanism remains largely unknown. Here we investigate whether and how hippocampal atrophy contributes to behavioral deficits based on a large-scale observational study conducted by the Alzheimers Disease Neuroimaging Initiative (ADNI). As a key novelty, we use 2D representations of the hippocampi, which allows us to better understand atrophy associated with different subregions. It, however, introduces methodological challenges as existing causal inference methods are not well suited for exploiting structural information embedded in the 2D exposures. Moreover, our data contain more than 6 million clinical and genetic covariates, necessitating appropriate confounder selection methods. We hence develop a novel two-step causal inference approach tailored for our ADNI data application. Analysis results suggest that atrophy of CA1 and subiculum subregions may cause more severe behavioral deficits compared to CA2 and CA3 subregions. We further evaluate our method using simulations and provide theoretical guarantees.


Introduction
Alzheimer's disease (AD) is an irreversible, progressive 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 (Morishima-Kawashima and Ihara, 2002) that may cause the abnormal aggregation and deposition of the amyloid plaques. 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 β amyloid (Aβ) fibrils, as evidenced by a corresponding drop in the levels of soluble Aβ-42 in cerebrospinal fluid (CSF) (Shaw et al., 2009;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β and tau impact the hippocampus (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 (Braak et al., 1993;Fox et al., 1996), becomes detectable by magnetic resonance imaging (MRI). 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. The aim of the paper is to answer whether (and how) brain changes, more specifically hippocampal changes, affect the behavioral deficits in Alzheimer's studies. 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). Specifically, pathological and MRI studies show that hippocampal formation is severely atrophied in established Alzheimer's disease (Fox et al., 1996). Previous studies have found strong associations between hippocampal atrophy and behavioral deficits (Whitwell et al., 2007). Recently, there has been a growing interest in going beyond association and assessing causation to better understand the causes of the behavioral deficits. One particular question of interest is whether hippocampal atrophy is a cause of behavioral deficits. Some neuroscientists conjecture that the association between hippocampal atrophy and behavioral deficits is causal, because the former destroys the connections that help the neuron communicate and results in a loss of function (www.healthline.com/health/brain-atrophy).
In this article, we will develop new statistical tools to assess the causal relationship between hippocampal atrophy and behavioral deficits in an Alzheimer's study. A randomized trial to alter the brain size is clearly unethical. Instead, we based our analysis on a large-scale observational study by the Alzheimer's Disease Neuroimaging Initiative (ADNI). The ADNI study is a world-wide research project aiming to improve the prevention and treatment of Alzheimer's disease. We use the Alzheimer's Disease Assessment Scale cognitive score to quantify behavioral deficits, where higher score indicates more severe behavioral deficits. Our exposure of interest is hippocampal morphometry surface measure, summarized as two 100 × 150 matrices corresponding to the left/right hippocampus. Each element of the matrices is a continuousvalued 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 atrophy (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;Zhang et al., 2016). For example, Li et al. (2007) showed that the surface measures of the hippocampus can 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 identify how local changes of hippocampal subfields affect the behavioral deficits. In addition, the ADNI study measures ultra-high dimensional baseline covariates, including more than 6 million genetic variants on several hundreds of subjects. The special structure of the ADNI data application presents new challenges for statistical analysis. First, unlike conventional causal analysis that deals with a scalar treatment, 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 in the case of univariate binary or continuous treatment. The conventional approach is to adjust for all baseline covariates (Rubin, 2009). Recent studies, however, have shown that inclusion of the variables only associated with the exposure but not the outcome (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 Richardson et al. (2018) and references therein for an overview. Therefore, when a large number of covariates are available, our goal is to include all confounders and precision variables, while excluding all instrumental variables and irrelevant variables (not related to both outcome and exposure). To accomplish this goal, Shortreed and Ertefaie (2017) proposed an outcome-adaptive lasso approach that selects appropriate covariates for inclusion in propensity score models to account for confounding bias and maintain statistical efficiency. Ertefaie et al. (2018) proposed a variable selection method using a penalized objective function that is based on both the outcome and treatment assignment models, facilitating confounder selection in high-dimensional settings. The confounder selection problem was also investigated in VanderWeele and Shpitser (2011); Wang et al. (2012); Vansteelandt et al. (2012); Wilson and Reich (2014); Schnitzer et al. (2016); Antonelli et al. (2019), among others. However, these methods are not directly applicable to our ADNI data for two reasons. First, all the methods are developed for univariate binary or continuous treatment, and thus they are not directly applicable to our 2D high-dimensional continuous imaging treatment. Second, we have ultra-high dimensional baseline covariates (more than 6 million covariates), which makes the current penalization or Bayesian based confounder selection procedures computationally infeasible.
To address these challenges, we develop a novel two-step approach to estimate the causal effect of high-dimensional 2D hippocampal surface exposure on the Alzheimer's behavior score, while taking into account the ultra-high dimensional baseline covariates. The first step is a fast screening procedure based on both outcome and treatment generating 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 Section 3.3. The second step is a penalized regression procedure for the outcome generating model to further exclude instrumental and irrelevant variables, and simultaneously estimate the treatment effect. 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;Khadka et al., 2016). Although we consider a similar set of variables and causal structure among these variables as illustrated later in Figure  2, our goal is to estimate the causal effect of hippocampal atrophy on behavioral deficits. In contrast, in mediation analysis (Figure 4), researchers are often interested in the causal effects of genetic factors on behavioral deficits, and how those are mediated through hippocampal atrophy. Direct application of methods developed for imaging genetics mediation analysis to our problem may select genetic factors that are a subset of confounders having effects on both hippocampal atrophy and behavioral deficits. As discussed above, omitting the precision variables may inflate the variance of our causal effect estimate, and hence these precision variables should in fact be included in the adjustment set.
We summarize our main contributions as follows. Our work holds the potential to be the very first study of its kind investigating the causal relationship between complex neuroimaging data and human behavior using observational data. By developing a two-step causal inference procedure for high dimensional confounder selection and 2D treatment effect estimation, we are able to identify important biomarkers that need to be adjusted in making causal inference and locate subregions of the hippocampus that affect the behavioral deficits. To the best of our knowledge, this is the first time that 2D continuous treatment is studied in causal inference literature. Our result promises to unravel the underlying causal mechanism of behavioral deficits and help early Alzheimer's disease diagnosis.

5
The rest of the article is organized as follows. Section 2 includes detailed data and problem description. We introduce our model and two-step causal inference approach in Section 3. We analyze the ADNI data and estimate the causal effect of hippocampal atrophy on behavioral deficits in Section 4. Simulations are conducted in Section 5 to evaluate the finite-sample performance of the proposed method. Section 6 investigates the theoretical properties of our procedure. We finish with a discussion in Section 7.

Data and problem description
Understanding how human brains work and how they connect to human behavior is a central goal in medical studies. Traditionally, cognitive neuroscientists have been probing the relationship between brain atrophy and behavior using neuroimaging data collected by means of Magnetic Resonance Imaging (MRI) and Functional Magnetic Resonance Imaging (fMRI). The results usually reflect a correlation rather than a causal relation. Cognitive neuroscientists habitually deny that imaging data can provide causal information about the relationship between brain and behavior. However, recent studies have found evidence for the causal relationship between atrophy of certain brain regions and behavioral deficits (Godefroy et al., 1998;Rigoux and Daunizeau, 2015;Weber and Thompson-Schill, 2010).
In this paper, we are interested in studying whether (and how) hippocampal atrophy affects 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 24 months after baseline measurements. The Alzheimer's Disease Assessment Scale cognitive score, often referred to as the ADAS-11 (Rosen et al., 1984), consists of 11 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. A higher ADAS-11 score indicates more severe behavioral deficits. The ADAS-11 score is a widely used measure of cognitive behavior for Alzheimer's disease, and has been shown to have high test-retest reliability (Rosen et al., 1984) and moderately high internal consistency (Weyer et al., 1997).
The exposure of interest is the baseline 2D surface data obtained from left/right hippocampus. The hippocampus surface data were preprocessed from the raw MRI data, which were collected across a variety of 1.5 Tesla MRI scanners with protocols individualized for each scanner. Standard T1-weighted images were obtained by using volumetric 3-dimensional sagittal MPRAGE or equivalent protocols with varying resolutions. The typical protocol includes: inversion time (TI) = 1000 ms, flip angle = 8 o , repetition time (TR) = 2400 ms, and field of view (FOV) = 24 cm with a 256 × 256 × 170 acquisition matrix in the x−, y−, and z−dimensions yielding a voxel size of 1.25 × 1.26 × 1.2 mm 3 . We adopted a surface fluid registration based hippocampal subregional analysis package (Shi et al., 2013), which uses isothermal coordinates and fluid registration to generate one-to-one hippocampal surface registration for surface statistics computation. It introduced two cuts on a hippocampal surface to convert it into a genus zero surface with two open boundaries. The locations of the two cuts were at the front and back of the hippocampal surface. By using conformal parameterization, it essentially converts a 3D surface registration problem into a 2D image registration problem. The flow induced in the parameter domain establishes high-order correspondences between 3D surfaces. Finally, the radial distance was computed on the registered surface. This software package and associated image processing methods have been adopted and described by various studies (Wang et al., 2011;Shi et al., 2014;Monje et al., 2013). After preprocessing, we obtained left and right hippocampal shape representations as two 100 × 150 matrices.
In the ADNI data, there are millions of potential confounders that one may need to adjust for, including the whole genome sequencing data from all of the 22 autosomes. For the genetic data, we applied the following preprocessing technique to the 756 subjects in ADNI1 study. The first line quality control steps include (i) call rate check per subject and per single nucleotide polymorphism (SNP) marker, (ii) gender check, (iii) sibling pair identification, (iv) the Hardy-Weinberg equilibrium test, (v) marker removal by the minor allele frequency, and (vi) population stratification. The second line preprocessing steps include removal of SNPs with (i) more than 5% missing values, (ii) minor allele frequency smaller than 10%, and (iii) Hardy-Weinberg equilibrium p-value < 10 −6 . The 503,892 SNPs obtained from 22 autosomes were included for further processing. MACH-Admix software (http://www.unc.edu/∼yunmli/MaCH-Admix/)  is applied on all the subjects to perform genotype imputation, using 1000G Phase I Integrated Release Version 3 haplotypes (http://www.1000genomes.org) (Consortium et al., 2012) as reference panel. Quality control was also conducted after imputation, excluding markers with (i) low imputation accuracy (based on imputation output R 2 ), (ii) Hardy-Weinberg equilibrium p-value 10 −6 , and (iii) minor allele frequency < 5%. Finally, 6, 087, 205 bi-allelic markers (including SNPs and indels) of 756 subjects were retained in the data analysis.
We deleted those subjects that have missingness of hippocampus shape representations or ADAS-11 score observed at Month 24, ending up with 584 subjects for the analysis. Our aim is to estimate the causal effect of the hippocampal surface exposure on the ADAS-11 score by adjusting for clinical measures including age, gender and length of education, and a suitable subset of the 6, 087, 205 bi-allelic markers using the method developed in the next section.

Basic set-up
Suppose we observe independent and identically distributed samples is a continuous outcome of interest, and X ∈ X ⊆ R s denotes a vector of ultra-high dimensional genetic (and clinical) covariates, where we assume s >> n. We are interested in characterizing the causal effects of the 2D exposure Z on the outcome Y .
Under the potential outcome framework, we define Y (z) to be the potential outcome that would have been observed under treatment level z ∈ Z ⊆ R p×q . For each unit, Y i = Y i (Z i ) corresponds to the potential outcome under Z i , the level of treatment received. We are interested in the mean potential outcome µ(z) = E[Y (z)], which we call effect surface.
We make the following assumptions for identifiability of effect surface, which are standard in the causal inference literature: and there is only one version of treatment z.
A directed acyclic graph showing the causal structure is given in Figure 2. Under Assumptions 1-3, the effect surface can be identified with observed data as Denote X i = (X i1 , . . . , X is ) T . Without loss of generality, we assume that X il has been standardized for every 1 ≤ l ≤ s, and Z i and Y i have been centered. To estimate E[Y | Z = z, X], we assume the following linear regression models: (1) In (1), β l is the effect of the l-th potential confounder X il on the behavioral outcome Y i , B ∈ R p×q is the coefficient matrix, representing effect of the 2D imaging treatment Z i on the behavioral outcome Y i , and i is the random error. The inner product between two matrices is defined as Z i , B = vec(Z i ), vec(B) , where vec(·) is a vectorization operator that stacks the columns of a matrix into a vector. Model (2), previously introduced in Kong et al. (2020), builds the relationship between the 2D imaging treatment and the potential confounders. The C l is a p × q coefficient matrix characterizing the effect of the lth covariate X il on the 2D imaging treatment Z i and E i is a p × q matrix of random errors with mean 0. The symbol " * " denotes element-wise multiplication. Define where we assume |M 1 | < n and |M 2 | < n; here |M 1 | and |M 2 | represent the number of elements in M 1 and M 2 respectively. Under models (1) and (2), the effect surface can be expressed as µ(z) = z, B . Therefore, the estimation of E[Y | Z = z, X] boils down to estimation of the matrix coefficient B.
To obtain an unbiased estimate of the effect surface, the first step is to perform variable selection in models (1) and (2). For all the covariates X l , we can group them into four categories. Let A = {1, . . . , s}, and denote C the indices of confounders, i.e. variables associated with both the outcome and the exposure; P the indices of precision variables, i.e. predictors of outcome, but not exposure; I the indices of instrumental variables, i.e. covariates that predict the exposure, but not the outcome; S the indices of irrelevant variables, i.e. covariates that are not related to both the outcome and the exposure. These sets can be written as C = {l ∈ A|β l = 0 and C l = 0}, P = {l ∈ A|β l = 0 and C l = 0}, I = {l ∈ A|β l = 0 and C l = 0} and S = {l ∈ A|β l = 0 and C l = 0}.
The ideal adjustment set includes all confounders to avoid bias and all precision variables to increase statistical efficiency, while excluding instrumental variables and irrelevant variables (Brookhart et al., 2006;Shortreed and Ertefaie, 2017). In other words, we aim to retain all covariates from M 1 = C ∪ P = {l ∈ A|β l = 0}, while excluding covariates from I ∪ S = {l ∈ A|β l = 0}.

Naive screening methods
To find the nonzero β l 's, a straightforward idea is to consider a penalized estimator obtained from the outcome generating model (1), where one imposes, say Lasso penalties, on β l 's. However, this is computationally infeasible in our ADNI data application as the number of baseline covariates s is over 6 million.
Consequently, it is important to employ a screening procedure (e.g. Fan and Lv, 2008) to reduce the model size. To find covariates X l 's that are associated with the outcome Y conditional on the treatment Z, one might consider a conditional screening procedure for model (1) (Barut et al., 2016). Specifically, one can fit the model Y i = X il β l + Z i , B + i for each 1 ≤ l ≤ s, obtain marginal estimates of β M Z l 's and then sort the | β M Z l |'s for screening. This procedure works well if the treatment variable Z is of low dimension as one only needs to fit low dimensional ordinary least squares (OLS) s times. However, in our ADNI data application, the imaging treatment Z is of dimension pq = 15, 000, and one cannot obtain an OLS as n < pq. Thus, to apply the conditional sure independence screening procedure to our application, one may need to solve a penalized regression problem for each 1 ≤ l ≤ s, such as arg min where P λ (B) is a penalty of B. One can then obtain the estimates β M Z l,λ , and then rank the | β M Z l,λ |'s. The procedure will be repeated for s > 6, 000, 000 times in ADNI data, and is computationally prohibitive for two reasons. First, the penalized regression problem is much slower to solve compared to the OLS. Second, selection of the tuning parameter λ based on grid search substantially increases the computational burden.
Alternatively one may apply the marginal screening procedure of Fan and Lv (2008) to model (1). Specifically, one may solve the following marginal OLS on each X il by ignoring the exposure Z i : The marginal OLS estimate has a closed form β M l = n −1 n i=1 X il Y i , and one can rank | β M l |'s for screening. Specifically, the selected sub-model is defined as where γ 1,n is a threshold. Computationally, it is much faster than conditional screening for model (1) as we only need to fit one dimensional OLS for s > 6, 000, 000 times. However, the procedure is likely to miss some important confounders. To see this, plugging model (2) into (1) yields Even in the ideal case when X il 's are orthogonal for 1 ≤ l ≤ s, β M l is not a good estimate of β l because of the bias term C l , B . Thus, we may miss some nonzero β l 's in the screening step if the effects of β l and C l , B are of similar magnitudes but different signs. We illustrate this point in Figures 6-7 in the simulations, in which cases the conventional marginal screening on (1) fails to capture some of the true confounders.

Joint screening
To overcome the drawbacks of the estimation methods discussed in Section 3.2, we develop a joint screening procedure, specifically for our ADNI data application. The procedure is not only computationally efficient, but can also select all the confounders and precision variables with high probability. The key insight here is that although we are interested in selecting important variables in the outcome generating model, this can be done much more efficiently by incorporating information from the treatment generating model. Specifically, let C M l = n −1 n i=1 X il * Z i ∈ R p×q be the marginal OLS estimate in model (2) for l = 1, . . . , s. Following Kong et al. (2020), the important covariates in model (2) can be selected by where || · || op is the operator norm of a matrix and γ 2,n is a threshold. We then take the union of covariates sets in (4) and (7) to define a new set Intuitively, M * 1 includes all the important confounders and precision variables except for X l 's, for which both | β M l | and || C l || op are large. These X l 's, however, will be included in M 2 and hence also M. In Section 6, we show that with properly chosen γ 1,n and γ 2,n , P (M 1 ⊂ M) → 1 as n → ∞. In practice, we recommend choosing γ 1,n and γ 2,n such that | M * 1 | = | M 2 | = k, where k is the smallest integer such that | M| ≥ n/ log(n) . The threshold n/ log(n) was suggested by Fan and Lv (2008).
In general the set M includes not only confounders and precision variables in M 1 = C P, but also instrumental variables in I and a small subset of the irrelevant variables S. Nevertheless, the size | M| is greatly reduced compared to all the baseline covariates. For example, in our data application, we reduce the s > 6, 000, 000 potential confounders to around 100 covariates through our joint screening procedure. This makes it feasible to perform the second step in our procedure, a refined penalized estimation of B based on the covariates {X il , l ∈ M}.

Estimation of causal effects
In this step, we aim to estimate the treatment effect by excluding the instrument variables I and irrelevant variables in S from M and keeping the remaining covariates. This can be done by solving the following optimization problem arg min Denote ( B, β) the solution to the above optimization problem. The estimated effect surface is given by µ(z) = z, B . The Lasso penalty on β l is used to exclude instrumental and irrelevant variables in M, whose corresponding coefficients β l 's are zero. The nuclear norm penalty ||·|| * , defined as the sum of all the singular value of a matrix, is used to achieve a low-rank estimate of B, where the low-rank assumption in estimating 2D structural coefficients is commonly used in the literature (Chen et al., 2013;Zhou and Li, 2014;Kong et al., 2020). To solve the minimization problem (8), 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).
Algorithm 1 Shrinkage thresholding algorithm to solve (8) 2. Repeat (a) to (f) until the objective function Q(β, B) converges: A sufficient condition for the convergence of {β (t) } t≥1 and {B (t) } t≥1 is that the step size δ should be smaller than or equal to 1/L f , where L f is the smallest Lipschitz constant of the function l(β, B) (Beck and Teboulle, 2009). In our case, L f is equal to λ max (X T new X new )/n, where λ max (·) denotes the largest eigenvalue of a matrix.
For the tuning parameters, we use five-fold cross validation to select λ 1 and λ 2 based on a two-dimensional grid search.

ADNI data appliations
We use the data obtained from the ADNI study (adni.loni.usc.edu). The data usage acknowledgement is included in Section ?? of the supplement material. As described in Section 2, our data set includes 584 subjects from the ADNI1 study. The treatment of interest is the baseline 2D hippocampal surface radial distance measures, which can be represented as a 100 × 150 matrix for each part of the hippocampus. The outcome of interest is the ADAS-11 score observed at Month 24. The average ADAS-11 score is 14.4 with standard deviation 11.2. The potential confounders to adjust for include 6, 087, 205 bi-allelic markers as well clinical covariates: age, gender and education length. The average age is 75.5 years with standard deviation 6.7, and the average education length is 15.7 years with standard deviation 2.9. Among all the 584 subjects, 58.2% is female. Studies have shown that age and gender are main risk factors for Alzheimer's disease (Vina and Lloret, 2010;Guerreiro and Bras, 2015) with older people and female 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 hippocampal atrophy ( Van de Pol et al., 2006;Jack et al., 2000;Noble et al., 2012). Therefore, we include age, gender and education length as the confounders, and adjust for them in our analysis. In addition, we also adjust for the 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 left and right hippocampi separately, i.e. treat each part as a treatment and apply our procedure twice.
In particular, for each part of hippocampal outcome Y i , we regress both of our outcome Y i and exposure Z i on the three clinical measures (age, gender and education length) as well as the top five principal components of the whole genome data, to obtain the residual Y * i and Z * i . We then apply our joint screening procedure on those SNPs by treating Y * i and Z * i as outcome and exposure, respectively. In other words, the first step screening is only performed on the 6, 087, 205 SNPs. We choose the thresholds γ 1,n and γ 2,n such that | M| = n/ log(n) = 91. In Table 1, we list the top 20 SNPs corresponding to left and right hippocampi respectively.
From Table 1, one can see that there are quite a few important SNPs found by analysis results 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, which is 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). The SNP has also been found to be the most significant SNP associated with Alzheimer's disease (Kamboh et al., 2012). In addition, a great portion of the SNPs in Table 1 have 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;Herold et al., 2016), an A/G snv variation in the APOC1 gene; rs111789331 (Gao et al., 2016;Lutz et al., 2019), an A/T snv variation; and rs4420638 (Coon et al., 2007), an A/G snv variation in the APOC1 gene.  2  rs11694935  19  rs56131196  19  rs56131196  22  rs134828  2  rs11696076  19  rs4420638  19  rs4420638  16  rs194326  2  rs11692218  19  rs484195  19  rs484195  16  rs60922721  2  rs10495737  19  rs111789331  19  rs111789331  13  rs342705   Table 1: The top 20 SNPs selected through our joint screening procedure. The left two columns correspond to results from the left hippocampus, and the right two columns correspond to results from the right hipppocampus.
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 has 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 atrophy rate.
In addition, there are also a few SNPs found to be associated with hippocampi, though no evidence in the literature shows that they are strongly related to Alzheimer's. For example, the SNP rs484195, a G/A/T/C snv variation in the APOC1 gene, was identified to strongly associated to the volumes of both left hippocampus and right hippocampus . The SNP rs10495737, a C/T snv variation on chromosome 2, was among the top SNPs that are associated with the left hippocampal surface (Huang et al., 2017).
We then perform our second-step estimation procedure for each part of the hippocampi. In particular, we let X M be those selected SNPs, the population stratification (top five principal components of the whole genome data) and the clinical measures age, gender and education and Z be the left/right hippocampal surface image matrix. We use five fold cross validation to select λ 1 and λ 2 . The estimates B's corresponding to each hippocampus are mapped onto a representative hippocampal surface and plotted in Figure 3(a). We have also plotted the hippocampal subfield (Apostolova et al., 2006) in Figure 3(b). From the result, one can see that the atrophies of CA1 as well as the subiculum subregions increase the ADAS-11 score, which translates to more severe behavioral deficits. Existing literature (Schönheit et al., 2004;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 atrophy of CA1 and subiculum may cause more severe behavioral deficits compared to CA2 and CA3 subregions.
For comparison, we conduct mediation analysis (Bi et al., 2017) using the same dataset. The mediation analysis considers the same variables, but takes the genetic covariates as the treatment, and studies causal effects of genetic factors on behavioral deficits, disentangling the total effect of genetic factors on the behavioral outcome into a direct effect and an indirect effect operating through the hippocampal imaging mediator. A directed acyclic graph showing the mediation structure is given in Figure 4.
We take each element of imaging treatment Z as an intermediate imaging measure, and in total we have 30, 000 scans for left and right hippocampi (15, 000 each). Following Bi et al. (2017), linear structural equation models are assumed for the relationships among each of the genetic variants, each of the 30, 000 scans and the behavior score, while conditioning on participants age, gender, education length and the top five principal components of the whole genome data corresponding to population stratification. The p-value for each SNP effect was determined by the Wald test. A threshold of p-value is considered for significance after Bonferroni correction, controlling for the 6, 087, 205 SNPs.
Similar to Bi et al. (2017), we identified the pairs of SNP and imaging measure, for which the direct effect of SNP on behavioral outcome, the effect of SNP on imaging measure and the effect of imaging measure on the behavioral outcome are all significant. We list those SNPs with at least one paired imaging measure (i.e. hippocampal imaging pixel) in Table 2.
One can see most of the SNPs in Table 2 are also selected by our procedure, including SNPs rs429358, rs4420638, rs56131196, rs73052335, rs10414043 and rs7256200. As the mediation   Chromesome number SNP name  19  rs429358  19  rs4420638  19  rs56131196  19  rs6857  19  rs12721051  19  rs73052335  19  rs10414043  19  rs7256200  19 rs59007384 analysis can be used to find a subset of confounders, while our method aims to select the set of true confounders and precision variables, these SNPs selected by both procedures are likely to be confounders.

Simulation for screening
We perform our screening procedure (denoted by "union") and report the coverage proportion of M 1 , which is defined as | M∩M 1 | |M 1 | , where the size of the selected set | M| changes from 1 to 100. In addition, we report the coverage proportion for each of the confounding and precision variables, i.e. each of the j's in the set M 1 = {1, 2, 3, 104, 105, 106}. All the coverage proportions are averaged over 100 Monte Carlo runs.
To control the changing size of the | M|, we first set | M * 1 | = | M 2 | = 1 by specifying appropriate γ 1,n and γ 2,n . Then we sequentially add two variables, one to M * 1 by increasing γ 1,n and one to M 2 by increasing γ 2,n , until | M| reaches 100. Note that we always keep | M * 1 | = | M 2 | in the procedure. We may not obtain all the sizes between 1 to 100 because | M| may increase by at most 2. Therefore, for those sizes that cannot be reached, we use a linear interpolation to estimate the coverage proportion of M by using the closest two end points.
From the plots, one can see that both the "intersection" and "outcome" methods miss the confounder X 3 with a very high probability even as the size of the selected set approaches 100. In contrast, our method can select X 3 with high probability when | M| is relatively small. For confounders X 1 and X 2 , all three methods perform similarly. For the precision variables, the "outcome" method and our "union" method performs similarly in covering these variables, while the "intersection" performs badly. Combining the results, one can see that our method performs the best as our method selects all the confounders and precision variables with high probabilities. In addition, we find that the coverage proportion of our method increases when the sample size increases, which validates the sure independence screening property developed  Figure 6: Simulation results for the case (n, s, σ 2 ) = (200, 5000, 1): Panels (a) -(f) plot the average coverage proportion for X l , where l = 1, 2, 3, 104, 105 and 106. Panels (a)-(c) correspond to strong outcome and weak exposure predictor, moderate outcome and moderate exposure predictor and weak outcome and strong exposure predictor; Panels (d)-(f) correspond to strong, moderate and weak predictors of outcome only. Panel (g) plots the average coverage proportion for the index set M 1 = {1, 2, 3, 104, 105, 106}. The x-axis represents the size of M, while yaxis denotes the average proportion. The green solid, the red dashed and the black dash dotted lines denote our union screening method, the outcome screening method, and the intersection screening method, respectively. 19 Figure 7: Simulation results for the case (n, s, σ 2 ) = (1000, 5000, 1): Panels (a) -(f) plot the average coverage proportion for X l , where l = 1, 2, 3, 104, 105 and 106. Panels (a)-(c) correspond to strong outcome and weak exposure predictor, moderate outcome and moderate exposure predictor and weak outcome and strong exposure predictor; Panels (d)-(f) correspond to strong, moderate and weak predictors of outcome only. Panel (g) plots the average coverage proportion for the index set M 1 = {1, 2, 3, 104, 105, 106}. The x-axis represents the size of M, while y-axis denotes the average proportion. The green solid, the red dashed and the black dash dotted lines denote our union screening method, the outcome screening method, and the intersection screening method, respectively. 20 in Section 6.

Simulation for estimation
In this part, we evaluate the performance of our estimation procedure after the first-step screening. For the size of M in the screening step, we set | M| = n/ log(n) , so that | M| = 37, 80 and 144 for sample sizes n = 200, 500 and 1000, respectively. We report the mean squared errors (MSEs) for β and B defined as ||β − β|| 2 2 and Ḃ − B 2 F , respectively. Table 3 summarizes the average MSEs for β and B among 100 Monte Carlo runs. We can see that the MSE decreases with the sample size, which validates the consistency results developed in Section 6. In this subsection, we study theoretical properties for our screening procedure. We let M 1 = {1 ≤ l ≤ s n :β * l = 0} =Ċ ∪Ṗ, whereĊ = {1 ≤ l ≤ s n :Ċ l = 0 andβ * l = 0} anḋ P = {1 ≤ l ≤ s n :Ċ l = 0 andβ * l = 0}. Hereβ * l andĊ l are the true values for β l and C l , respectively, andḂ is the true value of B.
We have the following theorems, where the assumptions needed are included in Section ?? of the supplementary material.
Since the screening procedure automatically includes all the significant covariates for small value of γ 1,n and γ 2,n , it is necessary to consider the size of M, which we quantify in Theorem 2.
Theorem 1 shows that if γ 1,n and γ 2,n are chosen properly, our screening procedure will include all significant variables with a high probability. Theorem 2 guarantees that the size of selected model from the screening procedure is only of a polynomial order of n even though the original model size is of an exponential order of n. Therefore, the false selection rate of our screening procedure vanishes as n → ∞, while the size of M grows in a polynomial order of n, where the order depends on two constants κ and τ defined in the supplementary material. The proofs of these two theorems are collected in the supplementary material as well.

Theory for two-step estimator
In this section, we develop a unified theory for our two-step estimator. In particular, we derive a non-asymptotic bound for the final estimates. We first introduce some notation.
The bound in Theorem 3 implies that the convergence rate of the proposed estimator θ λ is O(max{ n 2κ+τ −1 , n 1−2τ }). Here ι is a positive constant as defined in Assumption (A6) in the supplementary material, and r is the rank ofḂ. The convergence rate is controlled by κ and τ , where κ controls the exponential rates of model complexity that can diverge and τ controls the rate of largest eigenvalue of population covariance matrix that can grow. The proof of the theorem is deferred to the supplementary material.

Discussion
This paper answers the question whether (and how) hippocampal atrophy affects behavorial deficits in Alzheimer's studies. With the hippocampal morphemetry surface measure as an exposure, the proposed causal inference procedure is successfully applied to the ADNI data analysis and identifies important genetic biomarkers that need to be adjusted for in causal inference and located subregions of the hippocampus that may be causally associated with behavioral deficits. To the best of our knowledge, this is the first causal inference approach for 2D high-dimensional continuous treatment with ultra-high dimensional potential confounders. The unique features of the data motivate us to develop a computationally efficient two-step screening and estimation procedure, which can select biomarkers among more than 6 million potential confounders and estimate the causal effect simultaneously.
There are a number of important directions for future work. First, in this article, we assume that the potential genetic confounders only have main effects on both the treatment and the outcome. Further investigation is warranted to extend the proposed approach to the case with interaction effects, which is of great scientific interest in genome-wide association studies (Cordell, 2009). Second, it is an interesting topic to study the statistical inference procedure for our causal estimates and characterize the uncertainty of the causal estimand. One may consider either bootstrap or debiased approaches (Chatterjee and Lahiri, 2011;Van de Geer et al., 2014). Note that existing approaches focus on the 1 penalization, whereas our estimation procedure is a combination of 1 sparse and nuclear norm low-rank penalization. Extending existing approaches to the low-rank case is highly non-trivial. Third, in the ADNI dataset, there are longitudinal ADAS-11 scores observed at different months and also 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 requires substantial effort for further invesitgation. Fourth, this study simply removes observations that have missingness. Accommodation of missing treatment, confounders and outcome under the proposed model framework is of great practical value, and we leave it for future research. Finally, it would be interesting to extend the current approach to MRI or fMRI imaging exposure, where one can combine the joint screening procedure with a second-step estimation procedure for tensorvalued images (Feng et al., 2020).

Supplementary Material
Supplementary Material available online contains detailed derivation and explantion of the Algorithm 1, ADNI data usage acknowledgement, additional simulation results, assumptions needed for our main theorems, and proofs of auxiliary lemmas and main theorems.