Multivariate Analyses with Two-Step Dimension Reduction for an Association Study Between 11C-Pittsburgh Compound B and Magnetic Resonance Imaging in Alzheimer’s Disease

Abstract

The neuropathological diagnosis of Alzheimer’s disease (AD) relies on amyloid beta (Aβ) deposition in brain tissues. To study the relationship between Aβ deposition and brain structure, as determined using ¹¹C-Pittsburgh compound B (PiB) and magnetic resonance imaging (MRI), respectively, we developed a regression model with PiB and MRI data as the predictor and response variables, respectively, and proposed a regression method for studying the association between them based on a supervised sparse multivariate analysis with dimension reduction based on a composite paired basis function. By applying this method to imaging data of 61 patients with AD (age: 55–85), the first component showed the strongest correlation with the composite score, owing to the supervised feature. The spatial pattern included the hippocampal and parahippocampal regions for MRI. The peak value was observed in the posterior cingulate and precuneus for PiB. The differences in PiB scores among the diagnosis groups 12 months after PiB imaging were significant between the normal and AD groups (p = 0.0284), but not between the normal and mild cognitive impairment (MCI) groups or the MCI and AD groups (p = 0.3508). Our method may facilitate the development of a dementia biomarker from brain imaging data. Scoring imaging data allows for visualization and the application of traditional analysis, facilitating clinical analysis for better interpretation of results.

1. Introduction

Alzheimer disease (AD) is the most common form of dementia, and its incidence is growing worldwide [1]. The neuropathological diagnosis of AD relies on the deposition of amyloid beta (Aβ) in the brain tissue. The spatial pattern of Aβ deposition can be visualized by administering the radiotracer ¹¹C-Pittsburgh compound B (PiB) and performing subsequent positron emission tomography (PET) [2]. When performing regression analysis, response variables based on cognitive function tests are often used. Chapleaut et al. [3] provide a comprehensive review of the latest applications of amyloid PET in neurodegenerative diseases. Amyloid PET is important in the diagnosis of AD and can minimally invasively detect amyloid plaques. With the approval of anti-amyloid monoclonal antibodies, their use has increased, helping treat mild cognitive impairment (MCI) and mild dementia. It has also been applied to neurodegenerative diseases other than AD. In particular, it has been applied to the study of dementia with Lewy bodies (DLB), Parkinson’s disease dementia (PDD), and Parkinson’s disease (PD) and some specific studies are presented below. The study by Edison et al. [4] measured amyloid load in DLB and PDD patients using ¹¹C-PiB PET; they observed a marked increase in amyloid load in approximately 80% of DLB patients, while PDD patients showed less of an increase in amyloid. A study by Donaghy et al. [5] evaluated the utility of amyloid PET in DLB. In particular, DLB patients showed a greater accumulation of amyloid plaques than DLB patients. In particular, we found that DLB patients often had less amyloid accumulation; DLB patients tended to have higher Aβ ligand binding than PDD patients. PD without dementia does not show increased binding. The Walker et al. [6] study examined the prevalence, distribution, and severity of cerebral amyloid angiopathy (CAA) in DLB, PDD, and PD, showing that the distribution and severity of CAA in DLB is similar to that in AD. This may be related to the development of dementia. Amyloid PET images such as PiB image are interpreted visually and evaluated using quantification methods such as SUV ratios. It is also important in clinical trials and utility studies, where it is compared with other biomarkers and analyzed in detail. This approach is less informative because the diagnosis is a univariate, categorical variable. Indeed, cognitive function tests play an important role in the diagnosis of AD, such as the Mini-Mental State Examination (MMSE), a simple cognitive assessment tool that evaluates memory, attention, computation, language, and visuospatial cognition, the Montreal Cognitive Assessment, which provides a more detailed assessment than MMSE and is particularly superior in detecting MCI, and the clock-drawing test, a simple test that assesses visuospatial cognition and executive function. These tests may include information directly related to the illness, as well as other factors. This is because of the wide range of difficulties of the questions and their content. Additionally, patients’ cultural backgrounds and educational levels may influence test results. Thus, in this study, images that can assess pathophysiology more directly than cognitive test functions are used as the objective variable.

Magnetic resonance imaging (MRI) is a useful tool for examining the structural changes caused by Aβ deposition in the brain and provides more relevant response information than cognitive tests. This concept is especially useful in imaging genetics [3]. In this context, brain imaging data were used to identify the intermediate phenotypes. Multiple datasets have been classified into multimodal imaging analyses in neuroinformatics research and multiomics in bioinformatics research, which are among the most widely explored topics. These are being advanced through the development of new methodologies and technologies and have the potential for major breakthroughs in the understanding of human health and disease [7]. Consequently, the use of structural MRI (sMRI) as a response is expected to yield more information through regression analysis. However, sMRI provides high-dimensional data with millions of voxels as predictor variables, making statistical analysis very challenging unless dimensionality reduction is performed.

Generally, in multivariate analyses of brain imaging using response and predictor variables, partial least squares (PLS), canonical component analysis (CCA), reduced rank regression, and parallel independent component analysis (ICA) have been used. These approaches are common because the datasets are scored using a linear combination but differ in the objective function used to compute weights. In brain imaging analyses, methods other than parallel ICA can be considered equivariant [8]. Although parallel ICA offers high decomposability, its incorporation of sparsity is challenging, leading to increased use of other methods. Lin et al. [9] showed that sparse models are effective in the integrated analysis of high-dimensional and heterogeneous imaging and genetic data. In the present study, an extension of the PLS or CCA was considered. Including imaging data as a response variable presents challenges owing to its large dimensionality and high correlations. Because voxel-based imaging data serve as both predictors and responses in the model, a multivariate approach with dimensionality reduction would be ideal. Although PLS and CCA can deal with correlated data and reduce dimensions, they may lead to the loss of important information, such as spatial relationships. Moreover, the objective function for weights maximizes the correlation between scores, and the resulting scores are not always related to the outcomes. Kawaguchi and Yamashita [10] developed a supervised and multiblock version of the PLS or CCA to analyze multimodal imaging and genetic data. In this study, we aimed to evaluate the relationship between Aβ deposition and brain structure, as determined using PiB and MRI. Moreover, we identified the optimal brain regions for assessing this relationship. A regression model was used with PiB data as the predictor and MRI data as the response variable. Furthermore, we propose a regression method for studying the association between PIB and MRI based on a supervised sparse multivariate analysis (SSMA) method with dimension reduction based on a composite paired basis function. This is a modified method of that proposed by Yoshida et al. [11] (radial basis function-sparse PLS without supervision and with a single basis function) and that of Kawaguchi and Yamashita [10] (multiblock sparse multivariable analysis). Initially, spatial brain images are reduced in dimensionality via basis expansion, followed by further dimensionality reduction using regularized matrix factorization. This is supervised by a univariate composite score representing a linear combination of cognitive test scores, enabling the simultaneous data-driven selection of relevant brain regions. While previous studies take the means of condensing to anatomical regions due to the huge voxel values, the proposed method is able to capture voxel-level relationships among brain images, and furthermore, super vision improves interpretability. These are considered to be the main advantages of the proposed method. Because it is a flexible method that is not restricted to specific information, we believe it has the potential to be applicable to a wider variety of clinical applications.

2. Materials and Methods

2.1. Data

This method was applied to real data from the AD Neuroimaging Initiative, which is a collection of imaging data from 61 individuals aged 55–85. We used MRI data obtained 12 months after PiB imaging.

2.2. Preprocessing

For image preprocessing, three-dimensional (3D) T1-weighted images of the participants were first segmented into gray matter (GM), white matter, and cerebrospinal fluid space, followed by anatomical normalization to a template image using DARTEL [12]. For PiB, images were coregistered to the bias-corrected T1 image, and the same warping parameter as that for T1 images was applied following intensity normalization to the cerebellar GM for each participant.

2.3. Notation

${\mathbf{S}}_Y$ and ${\mathbf{S}}_X$ are the $n\times N$ matrices for PiB and sMRI data, respectively. Each row corresponds to vectorized 3D imaging data for each participant.

2.4. Dimension Reduction

In the first phase, the number of dimensions for the imaging data are reduced by applying basis expansion. Since sMRI and PiB data have the same dimensions, we apply a common basis function to reduce the number of parameters from $N$ to $q$. $\mathit{Y}={\mathbf{S}}_Y\mathbf{B}$ and $\mathit{X}={\mathbf{S}}_X\mathbf{B}$ are the $n\times q$ matrix, where ${\mathbf{S}}_m=\left\{s_{m\alpha }\right\}$ is an $n\times N$ matrix, and $\mathbf{B}={\left\{{\varphi }_k\left({\mathit{v}}_j\right)\right\}}_{j=1,\dots ,N,k=1,\dots ,q}$ is an $N\times q$ matrix. Each element in B is defined by a radial B-spline function, which provides a basis for reducing the dimensions of the original data.

${\varphi }_k\left(\mathit{w}\right)={\displaystyle \frac{1}{4h^2}}\left\{\begin{array}{c}{h}^3-3h^2{d}_k\left(\mathit{w}\right)+3h{d_k\left(\mathit{w}\right)}^2+3{d_k\left(\mathit{w}\right)}^3\\ {\left(h-d_k\left(\mathit{w}\right)\right)}^3\\ 0\end{array}\begin{array}{c}\left(d_k\left(\mathit{w}\right)\le 0\right)\\ \left(0<d_k\left(\mathit{w}\right)\le h\right)\\ \left(d_k\left(\mathit{w}\right)>h\right)\end{array}\right.,$

where $d_k\left(\mathit{w}\right)=\parallel \mathit{w}-{\kappa }_i\parallel -h$, ${\kappa }_i\in {\mathbb{Z}}^3$ is the pre-specified knot, $h>0$ is the distance between knots, and $‖\mathit{a}‖=\sqrt{a_1^2+a_2^2+a_3^2}$ for 3D vector $\mathit{a}$. Notably, basis expansion was useful for 3D neuroimaging analysis, as also reported by Araki et al. [13,14], Yoshida et al. [11,15], Kawaguchi and Yamashita [10], and Kawaguchi [16].

The number of basis functions, $q$, is determined by the number of voxels and the distance between the prespecified knots. In this study, we employed four-voxel equal-spacing knots, resulting in a distance $h=\sqrt{3\times 4^2}=6.93$. This choice was influenced by our previous simulation study [8], which indicated that a smaller distance between knots could improve accuracy.

Dimension reduction based on the basis function was followed by the application of the SSMA method. Thus, the result is obtained as a linear combination of the basis functions ${\psi }_X\left({\mathit{v}}_j\right)={\sum }_{k=1}^q{\varphi }_k\left({\mathit{v}}_j\right)w_{X,m}$ and ${\psi }_Y\left({\mathit{v}}_j\right)={\sum }_{k=1}^q{\varphi }_k\left({\mathit{v}}_j\right)w_{Y,m}$, which is referred to as a composite function (Figure 1). ${\psi }_X\left(\cdot \right)$ and ${\psi }_Y\left(\cdot \right)$ are paired because of the hierarchical structure. This can be regarded as a paired composite function, which represents related brain regions across modalities. The composite basis $\psi \left(\mathit{v}\right)$ can exhibit a flexible shape, whereas the original basis function ${\varphi }_k\left(\mathit{v}\right)$ maintains a spherical shape.

2.5. Supervised Sparse Multivariate Analysis

In this section, we present the proposed SSMA method. Consider $n$ subject, with $\mathit{X}$ as the $n\times p$ predictor matrix and $\mathit{Y}$ as the $n\times q$ response matrix. Each subject also has a univariate outcome represented by the $n$-dimensional vector $\mathit{Z}$. We define scores $\mathit{t}$ for the predictor and $\mathit{u}$ for the response, with the following structure:

(1)$\mathit{t}=\mathit{X}{\mathit{w}}_X, \mathit{u}=\mathit{Y}{\mathit{w}}_Y,$

where ${\mathit{w}}_X$ and ${\mathit{w}}_Y$ are the weight vectors for $\mathit{X}$ and $\mathit{Y}$, respectively. Subsequently, matrices $\mathit{X}$ and $\mathit{Y}$ are normalized by their columns. The association between the two datasets (images) was evaluated by maximizing the correlation between the super scores $\mathit{t}$ and $\mathit{u}$ supervised by the outcomes. Weights ${\mathit{w}}_X$ and ${\mathit{w}}_Y$ were estimated by maximizing the following function:

(2)$\mathrm{L}\left({\mathit{w}}_X,{\mathit{w}}_Y\right)={\mu }_{XY}{\mathit{t}}^{\mathrm{\top }}\mathit{u}+{\mu }_{XZ}{\mathit{t}}^{\mathrm{\top }}\mathit{Z}+{\mu }_{YZ}{\mathit{u}}^{\mathrm{\top }}\mathit{Z}-{\mathrm{P}}_{{\lambda }_X}\left({\mathit{w}}_X\right)-{\mathrm{P}}_{{\lambda }_Y}\left({\mathit{w}}_Y\right)$

$\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t} \mathrm{t}\mathrm{o} \Vert {\mathit{w}}_X{\Vert }^2=1 \mathrm{a}\mathrm{n}\mathrm{d} \Vert {\mathit{w}}_Y{\Vert }^2=1$, where $0\le {\mu }_{XY},{\mu }_{XZ},{\mu }_{YZ}\le 1$ and ${\mu }_{XY}+{\mu }_{XZ}+{\mu }_{YZ}=1$ are proportions of the supervision. The penalty function ${\mathrm{P}}_{\lambda }\left(x\right)$ is defined as ${\mathrm{P}}_{\lambda }\left(x\right)=2\lambda \left|x\right|$ in this paper, where $\lambda >0$ is the regularization parameter that controls sparsity, thereby facilitating the detection of associated regions in the brain or genome. It is possible for the penalty function to have other structures; however, these structures are beyond the scope of this study. When the proportions of supervision for the predictor-response pairs (${\mu }_{XZ}$ and ${\mu }_{YZ}$) are both set to zero, the proposed method simplifies to either multiblock PLS or CCA.

We estimate the weights in (1) by stating the following proposition.

Proposition 1 (coordinate updates).

The solution of optimization problem (2) satisfies

${\tilde{\mathit{w}}}_X=h_{{\lambda }_X}\left({\mathit{X}}^{\mathrm{\top }}\left\{{\mu }_{XY}\mathit{u}+{\mu }_{XZ}\mathit{Z}\right\}\right), {\tilde{\mathit{w}}}_Y=h_{{\lambda }_Y}\left(\left\{{\mu }_{XY}{\mathit{t}}^{\mathrm{\top }}+{\mu }_{YZ}{\mathit{Z}}^{\mathrm{\top }}\right\}\mathit{Y}\right)$

where $h_{\lambda }\left(y\right)=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(y\right){\left(\left|y\right|-\lambda \right)}_{+}$ is the sparse function.

The proof can be obtained in a manner similar to that reported by Kawaguchi and Yamashita [10]. This leads to the following algorithm for estimating the weights in (1) by maximizing L in (2).

1. Initialize $\mathit{t}$ and $\mathit{u}$ and normalize the super scores as follows:

   $\mathit{t}\leftarrow \mathit{t}/{‖\mathit{t}‖}_2, \mathit{u}\leftarrow \mathit{u}/{‖\mathit{u}‖}_2$

   where $\leftarrow$ means “replaced with”.

2. Repeat until convergence.

   • 2.1. For fixed $\mathit{u}$, ${\tilde{\mathit{w}}}_X=h_{{\lambda }_X}\left({\mathit{X}}^{\mathrm{\top }}\left\{{\mu }_{XY}\mathit{u}+{\mu }_{XZ}\mathit{Z}\right\}\right)$ and normalize ${\widehat{\mathit{w}}}_X$ = ${\tilde{\mathit{w}}}_X/{‖{\tilde{\mathit{w}}}_X‖}_2$.

   • 2.2. Putting $\mathit{t}=\mathit{X}{\widehat{\mathit{w}}}_X$, ${\tilde{\mathit{w}}}_Y=h_{{\lambda }_Y}\left(b_Y\left\{{\mu }_{XY}{\mathit{t}}^{\mathrm{\top }}+{\mu }_{YZ}{\mathit{Z}}^{\mathrm{\top }}\right\}\mathit{Y}\right)$ then normalize ${\widehat{\mathit{w}}}_Y$ = ${\tilde{\mathit{w}}}_Y/{‖{\tilde{\mathit{w}}}_Y‖}_2$.

   • 2.3. Set $\mathit{u}=\mathit{Y}{\widehat{\mathit{w}}}_Y$,

3. Set $\mathit{p}={\mathit{X}}^{\mathrm{\top }}\mathit{t}/{\mathit{t}}^{\mathrm{\top }}\mathit{t}$ and $\mathit{q}={\mathit{Y}}^{\mathrm{\top }}\mathit{u}/{\mathit{u}}^{\mathrm{\top }}\mathit{u}$, $\mathit{X}\leftarrow \mathit{X}-\mathit{t}{\mathit{p}}^{\top }$ and $\mathit{Y}\leftarrow \mathit{Y}-\mathit{u}{\mathit{q}}^{\top }.$

Note that by iterating Steps 1–3, the deflation step (Step 3) yields multiple components ($k=\mathrm{1,2},\dots ,K$), the new $\mathit{X}$ and $\mathit{Y}$ obtained in Step 3. The convergence of this algorithm is supported by the proposition provided by Kawaguchi and Yamashita [10]. The method for selecting the regularization parameter was consistent with the approach reported by Kawaguchi and Yamashita [10], which is based on the information criterion. The R package “msma” was used to implement the method described in this section and is available from the Comprehensive R Archive Network at http://CRAN.R-project.org/package=msma (accessed on 8 November 2024).

2.6. Statistical Analysis

The relationships between sMRI and PiB scores estimated using SSMA were evaluated using Pearson’s correlation coefficients. The scores are calculated as standardized continuous values and outliers are rarely observed. Therefore, Pearson’s correlation coefficient is considered to be applicable. The resulting PiB scores in the SSMA were tested using Tukey’s multiple comparison test. This validation was performed by using non-imaging data of future diagnoses (normal, MCI, and AD groups) and comparing them between groups. The comparison was limited to cases with diagnostic information. Statistical significance was set at p < 0.05.

2.7. Simulation Study Setting

In this section, the characteristics of the proposed method are presented based on synthetic data. Specifically, the aim is to demonstrate how the supervised feature functions work, mainly by comparing the original method with the effects of varying ${\mu }_{YZ}$ and ${\mu }_{XZ}$. Notably, as described in Step 2, the case of ${\mu }_{XZ}={\mu }_{YZ}=0$ corresponds to the original method (unsupervised version).

We considered a sample size of $n=50$ subjects and set the number of components to $K=2$. The true images were 2D with $100\times 100$ pixels ($q=\mathrm{10,000}$), as shown in Figure 2. The observed data consisted of $\left(\mathit{X},\mathit{Y},\mathit{Z}\right),$ where $\mathit{X}$ is the predictor matrix, $\mathit{Y}$ is the response matrix, and $\mathit{Z}$ is the outcome vector. Both $\mathit{X}$ and $\mathit{Y}$ are $n\times \mathrm{10,000}$ matrices generated from the $n\times K$ (latent) scores $\mathit{t}$ and $\mathit{u}$, respectively, with the true images serving as weights. The number of columns in $\mathit{t}$ and $\mathit{u}$ corresponds to the number of components. The outcome measured by $\mathit{Z}$ is influenced by the latent factors represented in $\mathit{u}$.

First, $\mathit{t}$ is obtained by generating a multivariate normal random number from $N\left(\mathbf{0},{\mathit{I}}_2\right)$. Subsequently, $\mathit{u}$ is generated from the relationship between $\mathit{t}$ and $\mathit{u}$, represented as

$\mathit{u}=\mathit{t}{\mathit{w}}_u^{\mathrm{\top }}+\mathsf{\epsilon },$

where ${\mathit{w}}_u={\left(1,w_{u,0}\right)}^{\mathrm{\top }}$, $w_{u,0}=0.2 \mathrm{o}\mathrm{r} 0.8$, and $\mathsf{\epsilon }\sim N\left(\mathbf{0},{\mathit{I}}_2\right)$. In the second component, the relationships between scores are considered to be either weak or strong.

The predictor and response matrices were generated as follows. The predictor is given by

$\mathit{X}=\mathit{t}{\mathit{w}}_X^{-}$

where ${\mathit{x}}^{-}$ is the generalized inverse for $\mathit{x}$. The weight for images was ${\mathit{w}}_X={\mathit{s}}_X{w}_X$ with the $k$-th column of ${\mathit{s}}_X$ as ${\mathit{s}}_{X,k}={\left(s_{X,k}\left({\mathit{v}}_1\right),s_{X,k}\left({\mathit{v}}_2\right),\dots ,s_{X,k}\left({\mathit{v}}_q\right)\right)}^{\mathrm{\top }}$ is the vectorized true images and $w_X=$ 5. The response is given by

$\mathit{Y}=\mathit{u}{\mathit{w}}_Y^{-}+\mathit{E}$

where ${\mathit{w}}_Y={\mathit{s}}_Y{w}_Y$ with the $k$-th column of ${\mathit{s}}_Y$ as ${\mathit{s}}_{Y,k}={\left(s_{Y,k}\left({\mathit{v}}_1\right),s_{Y,k}\left({\mathit{v}}_2\right),\dots ,s_{Y,k}\left({\mathit{v}}_q\right)\right)}^{\mathrm{\top }}$ is the vectorized true images, $\mathit{E}=e_0{\mathit{E}}_0$ by considering the noise level $e_0=0.005$ and $n\times q$ matrix ${\mathit{E}}_0$ with $(i,j)$-th element $e_{0,i,j}\sim N\left(0,1\right)$. The outcome is given by $\mathit{Z}=\mathit{u}{\mathit{w}}_z^{\mathrm{\top }}$ with ${\mathit{w}}_z={\left(1,0\right)}^{\mathrm{\top }}$; briefly, only the first component is related with the outcome.

Thus, the case $w_Y=1$ represents a strong effect of the score on the observed data, and the case $w_Y=20$ ($w_Y^{-1}=0.05$) represents a weak effect. Thus, only the Y-score for the second component is associated with the outcome.

The SSMA method is applied to observed data $\left(\mathit{X},\mathit{Y},\mathit{Z}\right)$ with parameter combinations with ${\mu }_{XZ}=0,0.25,0.5,0.75$ and ${\mu }_{YZ}=0,0.25,0.5,0.75$ satisfying ${\mu }_{XZ}+{\mu }_{YZ}\le 1$ by using two components. The iterations for generating data and applying the method were conducted 100 times for each parameter pair. The radial B-spline function with 2-pixel equal spacing knots is used.

3. Results

3.1. Application

Figure 3 illustrates the spatial pattern and correlation among scores by components for cases $\left({\mu }_{XZ}, {\mu }_{YZ}\right)=\left(\mathrm{0,0}\right)$ and $\left({\mu }_{XZ}, {\mu }_{YZ}\right)=\left(\mathrm{0,0.25}\right)$. Four components were selected for both cases. The first component exhibited the highest correlation with a composite score, influenced by the supervised nature of the analysis. In the second component, a relationship that makes sense from a biological perspective was found, suggesting that the data reflect known anatomical or functional relationships in the brain. The spatial pattern included the hippocampal and parahippocampal regions for MRI, while peak values for PiB data were observed in the posterior cingulate and precuneus. These are considered reasonable areas based on previous studies, and further description will be given in the discussion section. For all components, except for the second one in the case of $\left({\mu }_{XZ}, {\mu }_{YZ}\right)=\left(\mathrm{0,0.25}\right),$ the relationships among scores were stronger than those observed in the case of $\left({\mu }_{XZ}, {\mu }_{YZ}\right)=\left(\mathrm{0,0}\right)$.

The differences in PiB scores among the diagnosis groups 12 months after PiB imaging are presented in Figure 4. When comparing the three groups—normal, MCI, and AD—statistical significance was observed only in the first component of the analysis. This means that the differences among these groups were notable in that specific component. In contrast, in the second component, a significance level of approximately 10% was noted, which suggests that the differences among the groups were less pronounced or only marginally significant in this context. The differences in the PiB first component scores among the diagnosis groups 12 months after PiB imaging were significant between the normal and AD groups (p = 0.0284), but not between the normal and MCI groups or the MCI and AD groups (p = 0.3508).

3.2. Simulation Study

All resulting weights from the SSMA method were binarized, categorizing weights as either zero (being equal, unselected) or non-zero (selected). The c-index was then computed using the true images and variables. The c-index is defined as

c-index = sensitivity − (1 − specificity).

The average values of the c-index were derived from 100 simulations, using various parameter combinations, specifically $\left(w_{u,0},w_Y\right)=\left(\mathrm{0.2,20}\right), \left(0.8, 1\right),$ along with ${\mu }_{XZ}$, and ${\mu }_{YZ}$. These parameters influence how the model behaves and the results it generates. The results are summarized in Table 1 and Table 2 for $e_0=0.005$. Notably, to facilitate interpretation, $w_Y$ is displayed as the inverse $w_Y^{-1}$. The results were divided into X and Y, with the c-index for each component listed. The average values from the two components are presented in the third and sixth columns, while the last column displays the average value for both X and Y across the two components. The highest mean value among ${\mu }_{XZ}$ and ${\mu }_{YZ}$ for each $w_{u,0}$ and $w_Y$ is highlighted in bold. Figure 5 illustrates the 2D displays of the resulting probabilities for Y and X for $e_0=0.005$. These probabilities are computed from the pixel-wise mean value of the binary images, indicating whether a pixel is selected in the estimation.

In general, the supervised feature with ${\mu }_{XZ}>0$ or ${\mu }_{YZ}>0$ consistently resulted in better predictive accuracy, as evidenced by higher c-index values, than the unsupervised cases with ${\mu }_{XZ}=0$ and ${\mu }_{YZ}=0$.

In the case of $w_{u,0}=$ 0.2, which represents a lower association between the scores u and t of X and Y, respectively, the difference in ${\mu }_{YZ}$ may not be significant. Conversely, when $w_{u,0}=0.8$, which indicates a higher association between the scores of X and Y, the supervised method with ${\mu }_{YZ}>0$ demonstrated superior performance to the original method with ${\mu }_{YZ}=0$. Furthermore, in the case of $w_Y=$ 1, the supervised method with ${\mu }_{XZ}+{\mu }_{YZ}<1$ outperformed the original method with ${\mu }_{XZ}+{\mu }_{YZ}=1$. In contrast, when $w_Y=20$, the supervised method with ${\mu }_{XZ}+{\mu }_{YZ}=1$ outperformed the original method with ${\mu }_{XZ}+{\mu }_{YZ}>1$.

In the case of $w_Y=20$, detecting the true images in Figure 5 was challenging. In contrast, for $w_Y=1$, the true images were successfully identified. The best performance was observed for $\left({\mu }_{XZ}, {\mu }_{YZ}\right)=\left(\mathrm{0.25,0.5}\right)$. This was mainly owing to the data structure of the simulation, where only the second component was associated with the outcome. However, the over-supervised (${\mu }_{XZ}+{\mu }_{YZ}=1$) resulted in poor reconstruction of the true image. Other parameter-setting patterns ($\left(w_{u,0},w_Y\right)=\left(\mathrm{0.2,1}\right), \left(0.8, 20\right)$) are presented in Table A1 and Table A2 of Appendix A, and the findings are consistent with the previous results.

4. Discussion

This paper presents a novel method for analyzing high-dimensional datasets. In practical applications, the scores calculated by the proposed method allow the use of traditional statistical analysis methods such as analysis of variance, which is usually used in clinical data analysis. It is also possible to perform logistic regression analysis and calculate the probability of disease diagnosis according to the value of the score. From a methodological perspective, we developed an SSMA method that can determine the association between multivariable datasets incorporating clinical outcomes. From an application perspective, we developed a method to predict dementia from brain images and to identify associations with other brain images using the following key features: (1) a paired composite basis function with a data-driven shape, (2) data-driven brain regions, and (3) availability for brain-wide analysis. Many studies focus on the regions specified by previous studies before analyzing them. For example, in studies on AD, the hippocampal region may be focused on. Also, in analyses of the whole brain, anatomically segmented regions of interests such as the representative Automated Anatomical Labelling Atlas may be used. In contrast, the proposed method does not require a predefined segmentation for analysis on a voxel-by-voxel basis, and thus can be flexibly applied to other diseases. In this study, we adopt an application in well-studied AD that also allows the evaluation of regions. These features of our method allow us to derive associations that are easier to interpret and can be applied to other medical big data analyses.

CCA and PLS are traditional methods that are highly effective for regression analysis of multivariate data. In this study, we extend this method to facilitate regression analysis involving both explanatory and objective variables represented as image data. While CCA and PLS are dimensionality reduction methods, further dimensionality reduction is required to analyze images with large amounts of information. The final result achieved dimensionality reduction using a flexible representation of the brain region. Additionally, we incorporated disease-related variables beyond imaging data as “teacher” data during parameter estimation to obtain results that better reflect the pathological condition. Simulation studies and real data analyses demonstrated the usefulness of this method. The choice of teacher data is determined on a study-by-study basis, and the results may vary depending on this choice.

The spatial patterns identified in this study highlighted the hippocampal and parahippocampal regions as significant spatial patterns on MRI. In contrast, PiB-PET scans frequently revealed peak values in the posterior cingulate gyrus and precuneus. These regions are essential indicators for early diagnosis and progression assessment of AD. Yokoi et al. [17] used PiB-PET and THK5351-PET to examine the relationship between Aβ accumulation and functional networks in the posterior cingulate gyrus and precuneus of patients with AD. Results indicate that Aβ accumulation in these regions is strongly associated with cognitive decline. For instance, a meta-analysis by Ballmaier et al. [18] used MRI to assess atrophy of the hippocampus and posterior cingulate gyrus in patients with AD, demonstrating marked atrophy associated with the disease. Teipel et al. [19] further reported hypometabolism and hypoperfusion of the posterior cingulate gyrus and precuneus in MCI and AD. Additionally, Storandt et al. [20] employed PiB-PET scans to investigate the association between Aβaccumulation, cognitive decline, and brain volume loss, supporting the link between these pathological features and disease progression. Based on this, early detection of Aβ accumulation will allow for early intervention to slow the progression of AD. Specific treatment recommendations include anti-amyloid medications (e.g., aducanumab) and lifestyle modifications. A systematic review by Teipel et al. [19] evaluated different neuroimaging correlates in patients with MCI and AD. Hypometabolism and hypoperfusion of the posterior cingulate gyrus and precuneus emerged as key indicators for MCI and AD, underscoring these regions’ role in disease progression.

This study has some limitations. The model’s performance relies on the choice of “teacher” data during parameter estimation. As different studies may select different clinical variables as “teacher” data, results may vary, potentially limiting reproducibility and generalizability across studies with different data sources. Additionally, this study primarily focuses on data associated with AD. Expanding validation across other neurodegenerative conditions or broader, non-pathological datasets could reveal limitations or areas where the model may not generalize well. Future studies could investigate criteria for selecting “teacher” data, optimizing variable choices to enhance model robustness and transferability across different datasets and study designs. Furthermore, additional imaging modalities beyond MRI and PET, such as functional MRI, diffusion tensor imaging, or quantitative susceptibility mapping, could be incorporated by extending the proposed method to multi-block methods, improving the detection of biomarkers across a broader spectrum of brain function and structure. These verification studies will be the subject of future research.

The proposed method facilitates constructing a prediction model that uses images as the response variable, offering greater predictive objectives than other univariate variables, such as cognitive test results. This approach is expected to yield novel insights through exploratory analysis.

5. Conclusions

The reported method for studying the association between the two imaging modalities will facilitate the development of a biomarker for dementia using brain imaging data.

Author Contributions

Conceptualization, A.K. and F.Y.; methodology, A.K.; software, A.K.; formal analysis, A.K.; writing—original draft preparation, A.K.; writing—review and editing, F.Y.; visualization, A.K.; project administration, A.K.; funding acquisition, A.K. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

The ADNI database used in this study was approved by the Institutional Review Boards of all participating institutions, and written informed consent has been obtained from all participants.

Data Availability Statement

The data utilized in this study were obtained from the Alzheimer’s Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu).

Acknowledgments

The authors are grateful to the referees for their constructive comments. For this study, we used the ACCMS supercomputer at Kyoto University. Data used in the 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 the analysis or writing of this report. A complete listing of ADNI investigators can be found at: https://adni.loni.usc.edu/wp-content/uploads/how_to_apply/ADNI_Acknowledgement_List.pdf (accessed on 8 November 2024).

Conflicts of Interest

The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Composite paired basis functions. The heat color represents the shape of the function.

Figure 2. True images for the simulation study.

Figure 3. Spatial pattern results in the SSMA model.

Figure 4. PiB scores in the SSMA model.

Figure 5. Probability 2D maps for the result of the simulation study.

Table 1

Simulation results for $e_0=$ 0.005, $w_X$ = 5, $w_Y$ = 20, $w_{u,0}$ = 0.2.

${\mathit{\mu }}_{\mathit{X}\mathit{Z}}$ $, {\mathit{\mu }}_{\mathit{Y}\mathit{Z}}$	comp1	comp2	Xmean	comp1	comp2	Ymean	XYmean
0.00, 0	0.348	0.351	0.350	0.021	0.028	0.024	0.187
0.00, 0.25	0.356	0.358	0.357	0.020	0.027	0.023	0.190
0.00, 0.5	0.355	0.311	0.333	0.019	0.028	0.023	0.178
0.00, 0.75	0.364	0.283	0.324	0.021	0.023	0.022	0.173
0.25, 0	0.349	0.287	0.318	0.018	0.026	0.022	0.170
0.25, 0.25	0.349	0.280	0.314	0.018	0.023	0.021	0.168
0.25, 0.5	0.337	0.283	0.310	0.019	0.026	0.022	0.166
0.25, 0.75	0.381	0.389	0.385	0.020	0.013	0.017	0.201
0.50, 0	0.348	0.294	0.321	0.017	0.028	0.022	0.172
0.50, 0.25	0.350	0.292	0.321	0.017	0.021	0.019	0.170
0.50, 0.5	0.381	0.389	0.385	0.020	0.013	0.017	0.201
0.75, 0	0.348	0.280	0.314	0.018	0.028	0.023	0.168
0.75, 0.25	0.381	0.389	0.385	0.020	0.013	0.017	0.201

Table 2

Simulation results for $e_0=$ 0.005, $w_X$ = 5, $w_Y$ = 1, $w_{u,0}$ = 0.8.

${\mathit{\mu }}_{\mathit{X}\mathit{Z}}$ $, {\mathit{\mu }}_{\mathit{Y}\mathit{Z}}$	comp1	comp2	Xmean	comp1	comp2	Ymean	XYmean
0.00, 0	0.510	0.870	0.690	0.466	0.690	0.578	0.634
0.00, 0.25	0.504	0.862	0.683	0.463	0.690	0.577	0.630
0.00, 0.5	0.504	0.853	0.678	0.470	0.693	0.581	0.630
0.00, 0.75	0.499	0.853	0.676	0.473	0.686	0.580	0.628
0.25, 0	0.506	0.863	0.684	0.464	0.705	0.584	0.634
0.25, 0.25	0.499	0.864	0.681	0.460	0.704	0.582	0.632
0.25, 0.5	0.493	0.874	0.684	0.472	0.708	0.590	0.637
0.25, 0.75	0.375	0.602	0.489	0.362	0.008	0.185	0.337
0.50, 0	0.506	0.863	0.685	0.465	0.704	0.584	0.634
0.50, 0.25	0.489	0.868	0.679	0.457	0.701	0.579	0.629
0.50, 0.5	0.375	0.602	0.489	0.362	0.008	0.185	0.337
0.75, 0	0.505	0.874	0.690	0.465	0.704	0.584	0.637
0.75, 0.25	0.375	0.602	0.489	0.362	0.008	0.185	0.337

Appendix A

Table A1 and Table A2 list additional results for other parameter-setting patterns in the simulation study.

Table A1

Simulation results for $e_0$ = 0.005, $w_X$ = 5, $w_Y$ = 1, $w_{u,0}$ = 0.2.

${\mathit{\mu }}_{\mathit{X}\mathit{Z}}$ $, {\mathit{\mu }}_{\mathit{Y}\mathit{Z}}$	comp1	comp2	Xmean	comp1	comp2	Ymean	XYmean
0.00, 0	0.502	0.844	0.673	0.432	0.140	0.286	0.480
0.00, 0.25	0.504	0.829	0.666	0.439	0.139	0.289	0.478
0.00, 0.5	0.501	0.841	0.671	0.440	0.136	0.288	0.479
0.00, 0.75	0.504	0.800	0.652	0.437	0.132	0.284	0.468
0.25, 0	0.506	0.848	0.677	0.437	0.139	0.288	0.482
0.25, 0.25	0.503	0.841	0.672	0.439	0.136	0.288	0.480
0.25, 0.5	0.502	0.818	0.660	0.435	0.130	0.282	0.471
0.25, 0.75	0.336	0.634	0.485	0.297	−0.024	0.137	0.311
0.50, 0	0.505	0.842	0.674	0.433	0.139	0.286	0.480
0.50, 0.25	0.503	0.829	0.666	0.436	0.139	0.287	0.477
0.50, 0.5	0.336	0.634	0.485	0.297	−0.024	0.137	0.311
0.75, 0	0.506	0.840	0.673	0.439	0.140	0.289	0.481
0.75, 0.25	0.336	0.634	0.485	0.297	−0.024	0.137	0.311

Table A2

Simulation results for $e_0$ = 0.005, $w_X$ = 5, $w_Y$ = 20, $w_{u,0}$ = 0.8.

${\mathit{\mu }}_{\mathit{X}\mathit{Z}}$ $, {\mathit{\mu }}_{\mathit{Y}\mathit{Z}}$	comp1	comp2	Xmean	comp1	comp2	Ymean	XYmean
0.00, 0	0.344	0.383	0.364	0.021	0.026	0.024	0.194
0.00, 0.25	0.367	0.366	0.367	0.018	0.028	0.023	0.195
0.00, 0.5	0.343	0.360	0.352	0.019	0.030	0.024	0.188
0.00, 0.75	0.335	0.344	0.339	0.022	0.029	0.025	0.182
0.25, 0	0.354	0.304	0.329	0.019	0.030	0.024	0.177
0.25, 0.25	0.354	0.313	0.333	0.019	0.026	0.022	0.178
0.25, 0.5	0.346	0.330	0.338	0.018	0.029	0.024	0.181
0.25, 0.75	0.383	0.394	0.388	0.020	0.015	0.018	0.203
0.50, 0	0.343	0.295	0.319	0.017	0.031	0.024	0.171
0.50, 0.25	0.356	0.300	0.328	0.016	0.028	0.022	0.175
0.50, 0.5	0.383	0.394	0.388	0.020	0.015	0.018	0.203
0.75, 0	0.346	0.307	0.327	0.017	0.032	0.024	0.176
0.75, 0.25	0.383	0.394	0.388	0.020	0.015	0.018	0.203

References

1. Nichols, E.; Steinmetz, J.D.; Vollset, S.E.; Fukutaki, K.; Chalek, J.; Abd-Allah, F.; Abdoli, A.; Abualhasan, A.; Abu-Gharbieh, E.; Akram, T.T. et al. Estimation of the global prevalence of dementia in 2019 and forecasted prevalence in 2050: An analysis for the Global Burden of Disease Study 2019. Lancet Public Health; 2022; 7, pp. e105-e125. [DOI: https://dx.doi.org/10.1016/S2468-2667(21)00249-8] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34998485]

2. Adlard, P.A.; Tran, B.A.; Finkelstein, D.I.; Desmond, P.M.; Johnston, L.A.; Bush, A.I.; Egan, G.F. A review of β-amyloid neuroimaging in Alzheimer’s disease. Front. Neurosci.; 2014; 8, 327. [DOI: https://dx.doi.org/10.3389/fnins.2014.00327] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25400539]

3. Chapleau, M.; Iaccarino, L.; Soleimani-Meigooni, D.; Rabinovici, G.D. The role of amyloid PET in imaging neurodegenerative disorders: A review. J. Nucl. Med.; 2022; 63, (Suppl. 1), pp. 13S-19S. [DOI: https://dx.doi.org/10.2967/jnumed.121.263195] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35649652]

4. Edison, P.; Rowe, C.C.; Rinne, J.O.; Ng, S.; Ahmed, I.; Kemppainen, N.; Villemagne, V.L.; O’Keefe, G.; Någren, K.; Chaudhury, K.R. et al. Amyloid load in Parkinson’s disease dementia and Lewy body dementia measured with [¹¹C] PIB positron emission tomography. J. Neurol. Neurosurg. Psychiatry; 2008; 79, pp. 1331-1338. [DOI: https://dx.doi.org/10.1136/jnnp.2007.127878] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18653550]

5. Donaghy, P.; Thomas, A.J.; O’Brien, J.T. Amyloid PET imaging in Lewy body disorders. Am. J. Geriatr. Psychiatry; 2015; 23, pp. 23-37. [DOI: https://dx.doi.org/10.1016/j.jagp.2013.03.001] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23831180]

6. Walker, L.; Simpson, H.; Thomas, A.J.; Attems, J. Prevalence, distribution, and severity of cerebral amyloid angiopathy differ between Lewy body diseases and Alzheimer’s disease. Acta Neuropathol. Commun.; 2024; 12, 28. [DOI: https://dx.doi.org/10.1186/s40478-023-01714-7] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/38360761]

7. Xin, Y.; Sheng, J.; Miao, M.; Wang, L.; Yang, Z.; Huang, H. A review of imaging genetics in Alzheimer’s disease. J. Clin. Neurosci.; 2022; 100, pp. 155-163. [DOI: https://dx.doi.org/10.1016/j.jocn.2022.04.017] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35487021]

8. Liu, J.; Calhoun, V.D. A review of multivariate analyses in imaging genetics. Front. Neuroinform.; 2014; 8, 29. [DOI: https://dx.doi.org/10.3389/fninf.2014.00029] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24723883]

9. Lin, D.; Cao, H.; Calhoun, V.D.; Wang, Y.P. Sparse models for correlative and integrative analysis of imaging and genetic data. J. Neurosci. Methods; 2014; 237, pp. 69-78. [DOI: https://dx.doi.org/10.1016/j.jneumeth.2014.09.001] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25218561]

10. Kawaguchi, A.; Yamashita, F. Supervised multiblock sparse multivariable analysis with application to multimodal brain imaging genetics. Biostatistics; 2017; 18, pp. 651-665. [DOI: https://dx.doi.org/10.1093/biostatistics/kxx011]

11. Yoshida, H.; Kawaguchi, A.; Tsuruya, K. Radial basis function-sparse partial least squares for application to brain imaging data. Comput. Math. Methods Med.; 2013; 2013, 591032. [DOI: https://dx.doi.org/10.1155/2013/591032] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23762188]

12. Ashburner, J.; Friston, K.J. Voxel-based morphometry--The methods. Neuroimage; 2000; 11, pp. 805-821. [DOI: https://dx.doi.org/10.1006/nimg.2000.0582] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/10860804]

13. Araki, Y.; Kawaguchi, A.; Yamashita, F. Regularized logistic discrimination with basis expansions for the early detection of Alzheimer’s disease based on three-dimensional MRI data. Adv. Data Anal. Classif.; 2013; 7, pp. 109-119. [DOI: https://dx.doi.org/10.1007/s11634-013-0127-5]

14. Araki, Y.; Kawaguchi, A. Functional logistic discrimination with sparse PCA and its application to the structural MRI. Behaviormetrika; 2019; 46, pp. 147-162. [DOI: https://dx.doi.org/10.1007/s41237-019-00079-3]

15. Yoshida, H.; Kawaguchi, A.; Yamashita, F.; Tsuruya, K. The utility of a network–based clustering method for dimension reduction of imaging and non-imaging biomarkers predictive of Alzheimer’s disease. Sci. Rep.; 2018; 8, 2807. [DOI: https://dx.doi.org/10.1038/s41598-018-21118-1]

16. Kawaguchi, A. Multivariate Analysis for Neuroimaging Data; CRC Press: Boca Raton, FL, USA, 2021.

17. Yokoi, T.; Watanabe, H.; Yamaguchi, H.; Bagarinao, E.; Masuda, M.; Imai, K.; Ogura, A.; Ohdake, R.; Kawabata, K.; Hara, K. et al. Involvement of the precuneus/posterior cingulate cortex is significant for the development of Alzheimer’s disease: A PET (THK5351, PiB) and resting fMRI study. Front. Aging Neurosci.; 2018; 10, 304. [DOI: https://dx.doi.org/10.3389/fnagi.2018.00304]

18. Boccia, M.; Acierno, M.; Piccardi, L. Neuroanatomy of Alzheimer’s Disease and Late-Life Depression: A Coordinate-Based Meta-Analysis of MRI Studies. J. Alzheimers Dis.; 2015; 46, pp. 963-970. [DOI: https://dx.doi.org/10.3233/JAD-142955] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25869784]

19. Talwar, P.; Kushwaha, S.; Chaturvedi, M.; Mahajan, V. Systematic Review of Different Neuroimaging Correlates in Mild Cognitive Impairment and Alzheimer’s Disease. Clin. Neuroradiol.; 2021; 31, pp. 953-967. [DOI: https://dx.doi.org/10.1007/s00062-021-01057-7] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34297137]

20. Storandt, M.; Mintun, M.A.; Head, D.; Morris, J.C. Cognitive Decline and Brain Volume Loss as Signatures of Cerebral Amyloid-β Peptide Deposition Identified With Pittsburgh Compound B. Arch. Neurol.; 2009; 66, pp. 1476-1481. [DOI: https://dx.doi.org/10.1001/archneurol.2009.272] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20008651]