Group Factor Analysis for Alzheimer’s Disease

Abstract

For any neuroimaging study in an institute, brain images are normally acquired from healthy controls and patients using a single track of protocol. Traditionally, the factor analysis procedure analyzes image data for healthy controls and patients either together or separately. The former unifies the factor pattern across subjects and the latter deals with measurement errors individually. This paper proposes a group factor analysis model for neuroimaging applications by assigning separate factor patterns to control and patient groups. The clinical diagnosis information is used for categorizing subjects into groups in the analysis procedure. The proposed method allows different groups of subjects to share a common covariance matrix of measurement errors. The empirical results show that the proposed method provides more reasonable factor scores and patterns and is more suitable for medical research based on image data as compared with the conventional factor analysis model.

1. Introduction

Modern medical imaging techniques are capable of measuring human brain in vivo [1]. For instance, magnetic resonance (MR) imaging measures nuclei of atoms, and positron emission tomography detects the positron-emitting radionuclides to construct three-dimensional images. The imaging procedures are designed and settled before medical or cognitive experiments. Once the protocol is established, the laboratory and the hospital begin to recruit a variety of subjects of interest into experimental sessions. Errors resulting from individual scans are actually generated from common sources, such as the scanner, protocol, and software. Initial classification of subjects into groups can be realized by using clinical diagnosis, which may be uncertain to some extent, provided by physicians along with subjects’ anamnesis.

Conventional factor analysis [2] models reduce high-dimensional data into a few latent variables and assume that data 𝐱 were generated by a set of unobserved independent unit-variance Gaussian source 𝐟 plus uncorrelated zero-mean Gaussian random noise 𝐮, 𝐱=𝐿𝐟+𝐮, where 𝐿 is the factor loading matrix. The sample covariance of 𝐱 can be expressed as 𝐿𝐿^𝑇+Ψ, where Ψ is a diagonal covariance matrix of random noises. The goal of factor analysis is to find 𝐿 and Ψ that maximally fit the sample covariance [3–5]. The EM algorithm was proposed to estimate the matrices [6]. Factor analysis is commonly applied to the dataset as a whole or to different groups of data separately, which may result in factor patterns hard to interpret and limit the potential use of the method in a wider range of medical applications. In this study, we propose a mixture factor analysis model (MFAM) to assign a common covariance matrix of noises or measurement errors to different groups of subjects but to allow individual groups having their own latent structures. In the empirical application, we analyzed an Alzheimer’s disease (AD) dataset by first extracting the volumetric information from MR anatomical images for both healthy controls and the patients suffering either AD or mild cognitive impairment, followed by applying the proposed MFAM to the volumetric data.

2. Material and Method

2.1. The Model

Let 𝑀 be the number of subject groups. To find multiple sets of factor loadings, {𝐿_𝑗;𝑗=1,…,𝑀}, with the 𝐟 scores distributed as Gaussian within each group, the data vector can be decomposed into a linear combination of factor loadings for each group [7, 8], that is, 𝐿_𝑗∈𝑅^𝐷×𝐾, 𝐱=_𝑗𝜋_𝑗𝝁_𝑗+𝐿_𝑗×𝐟∣𝑤_𝑗+𝐮, (1) where 𝐱 is 𝐷-dimensional and each factor scores 𝐟∣𝑤_𝑗 has 𝐾 variables, that is, 𝐟∈𝑅^𝐾. The parameter 𝝅 is associated with the proportion of subjects in the 𝑗th group, 𝜋_𝑗=𝑝(𝑤_𝑗). The indicator variable 𝑤 is one, 𝑤_𝑗=1, when the data belongs to 𝑗th group, otherwise 𝑤 is set to zero, 𝑤_𝑗=0. The formula (1) using 𝝅 introduces the main difference from previous mixture models of factor analysis. The data vector 𝐱 need not be centered and the mean of the 𝑗th group data is 𝝁_𝑗. The covariance matrix of residuals 𝐮 is a diagonal matrix Ψ=diag[Ψ₁,Ψ₂,…,Ψ_𝐷]. The data distribution can be expressed as 𝑃(𝐱)=^𝑀_𝑗=1𝑃𝐱∣𝐟,𝑤_𝑗𝑃𝐟∣𝑤_𝑗𝑝𝑤_𝑗𝑑𝐟. (2) In this work, capitalized 𝑃 denotes the probability function of a vector or a matrix and lowercase 𝑝 denotes the probability function of a scalar. The factor scores are assumed to be distributed as Gaussian 𝑃𝐟∣𝑤_𝑗=𝑁(0,𝐼),∀𝑗. (3) The notation 𝐼 is the identity matrix of order 𝐷. The distribution of data 𝐱 in each group is given by 𝑃𝐱∣𝐟,𝑤_𝑗=𝑁𝝁_𝑗+𝐿_𝑗𝐟_𝑗,Ψ. (4) Based on the MFAM (2), the likelihood function 𝑄 is as follows: 𝑄=𝐸^𝑁𝑀_𝑖=1𝑗=1(2𝜋)^−𝐷/2|Ψ|^−1/2exp−𝐱_𝑖−𝝁_𝑗−𝐿_𝑗𝐟_𝑖^𝑇Ψ⁻¹×𝐱_𝑖−𝝁_𝑗−𝐿_𝑗𝐟_𝑖^𝑤𝑗, (5) where 𝐸 denotes the expectation. The 𝑁 is the number of data vectors (subjects) with subscript 𝑖 for the 𝑖th subject. We need to compute the expectation of the variables, 𝐸𝑤_𝑗𝐟_𝑖∣𝐱_𝑖=𝐸𝑤_𝑗∣𝐱_𝑖𝐸𝐟_𝑖∣𝑤_𝑗,𝐱_𝑖. (6) To estimate 𝑄 in (5), the posterior probability of the 𝑗th group is calculated as 𝑃𝑤_𝑗=𝑃∣𝐱𝐱∣𝑤_𝑗𝑃𝑤_𝑗𝑃(𝐱)=𝜋_𝑗𝑁𝐱−𝝁_𝑗,𝐿_𝑗𝐿^𝑇_𝑗+Ψ∑_𝑢𝜋_𝑢𝑁𝐱−𝝁_𝑢,𝐿_𝑢𝐿^𝑇_𝑢,+Ψ (7) where the probability of 𝐱 given 𝑤_𝑗 is 𝑃𝐱∣𝑤_𝑗=𝑁𝐱−𝝁_𝑗,𝐿_𝑗𝐿^𝑇_𝑗+Ψ. (8) The parameters 𝝅 in (7) is the prior probability derived from the clinical diagnosis. Therefore, the expectation of 𝑤_𝑗 given 𝐱_𝑖 in (6) is proportional to the numerator in (7), ℎ_𝑖𝑗=𝐸𝑤_𝑗∣𝐱_𝑖∝𝜋_𝑗𝑁𝐱_𝑖−𝝁_𝑗,𝐿_𝑗𝐿^𝑇_𝑗+Ψ. (9) To calculate (6), we consider that the posterior probability of 𝐟 given 𝐱 is 𝑃()=𝑃()𝑃(𝐟)𝐟∣𝐱𝐱∣𝐟𝑃(𝐱)∝exp−𝐟^𝑇𝐿^𝑇Ψ⁻¹𝐿+𝐼𝐟−2𝐟^𝑇𝐿^𝑇Ψ⁻¹𝐱. (10) After some arithmetic calculation, 𝑃(𝐟∣𝐱) can be expressed as 𝑃(𝐟∣𝐱)∼𝑁𝑅⁻¹𝐿^𝑇Ψ⁻¹𝐱,𝑅, (11) where 𝑅=(𝐿^𝑇Ψ⁻¹𝐿+𝐼). Hence, the expectation of 𝐟 given 𝐱 is 𝐸[𝐟∣𝐱]=𝑅⁻¹𝐿^𝑇Ψ⁻¹𝐱. (12) From above, 𝐸(𝐟_𝑖∣𝑤_𝑗,𝐱_𝑖) in (6) is calculated as 𝐸𝐟_𝑖∣𝑤_𝑗,𝐱_𝑖=𝑅_𝑗⁻¹𝐿^𝑇_𝑗Ψ⁻¹𝐱_𝑖−𝝁_𝑗, (13) where 𝑅_𝑗=(𝐿^𝑇_𝑗Ψ⁻¹𝐿_𝑗+𝐼), according to (12).

There is no constraint on those factor loadings 𝐿_𝑗. The estimation of 𝐿_𝑗 is simply the maximum of 𝑄. A convenient way to express 𝑄 in (5) is achieved by setting 𝐟_𝑖̃=[𝐟^𝑇_𝑖1]^𝑇 and 𝐿_𝑗=[𝐿_𝑗𝝁_𝑗]. The expected log likelihood function can be expressed as 𝐸log𝑄=𝐸log^𝑁𝑀_𝑖=1𝑗=1(2𝜋)^−𝐷/2|Ψ|^−1/2−1exp2𝐱_𝑖−𝐿_𝑗𝐟̃_𝑖^𝑇Ψ⁻¹×𝐱_𝑖−𝐿_𝑗𝐟̃_𝑖^𝑤𝑗=−𝐷×𝑁2()−𝑁×log2𝜋2|Ψ|−log_𝑖,𝑗12ℎ_𝑖𝑗𝐱^𝑇_𝑖Ψ⁻¹𝐱_𝑖−ℎ_𝑖𝑗𝐱^𝑇_𝑖Ψ⁻¹𝐿_𝑗𝐸𝐟_𝑖̃∣𝐱_𝑖,𝑤_𝑗+12ℎ_𝑖𝑗×trace𝐿^𝑇_𝑗Ψ⁻¹𝐿_𝑗𝐸𝐟_𝑖̃𝐟^𝑇_𝑖̃∣𝐱_𝑖,𝑤_𝑗. (14) To maximize 𝑄 with respect to 𝐿_𝑗, we equate the derivative of (14) to zero, 𝜕log𝐸[𝑄]𝜕𝐿_𝑗=−_𝑖ℎ_𝑖𝑗Ψ⁻¹𝐱_𝑖𝐸𝐟_𝑖̃∣𝐱_𝑖,𝑤_𝑗^𝑇+ℎ_𝑖𝑗Ψ⁻¹𝐿_𝑗𝐸𝐟_𝑖̃𝐟^𝑇_𝑖̃∣𝐱_𝑖,𝑤_𝑗^𝑇⟹𝐿=0_𝑗=_𝑖ℎ_𝑖𝑗𝐱_𝑖𝐸𝐟_𝑖̃∣𝐱_𝑖,𝑤_𝑗^𝑇×_𝑠ℎ_𝑠𝑗𝐸𝐟_𝑖̃𝐟^𝑇_𝑖̃∣𝐱_𝑠,𝑤_𝑗^𝑇−1, (15) where 𝐸𝐟_𝑖̃∣𝑤_𝑗,𝐱_𝑖=𝐸𝐟_𝑖∣𝑤_𝑗,𝐱_𝑖^𝑇1^𝑇. (16) All the variables are estimated by the EM algorithm. In the E-step, the algorithm computes the expectation of the factor scores in (6) and the second moment of the scores, 𝐸𝑤_𝑗𝐟𝐟^𝑇∣𝐱=𝐸𝑤_𝑗∣𝐱𝐸𝐟𝐟^𝑇∣𝑤_𝑗,𝐱, (17) by 𝐸𝑤_𝑗𝐟_𝑖𝐟^𝑇_𝑖∣𝐱_𝑖=ℎ_𝑖𝑗Cov𝐟_𝑖∣𝑤_𝑗,𝐱_𝑖+ℎ_𝑖𝑗𝐸𝐟_𝑖∣𝑤_𝑗,𝐱_𝑖𝐸𝐟_𝑖∣𝑤_𝑗,𝐱_𝑖^𝑇. (18) The covariance matrix of residual, Ψ, can be estimated by its inverse matrix, 𝜕𝑄𝜕Ψ⁻¹=−𝑁21Ψ−2_𝑖,𝑗ℎ_𝑖𝑗𝐱_𝑖𝐱^𝑇_𝑖+ℎ_𝑖𝑗𝐱_𝑖𝐸𝐟_𝑖̃∣𝐱_𝑖,𝑤_𝑗^𝑇𝐿^𝑇_𝑗−12ℎ_𝑖𝑗𝐿_𝑗𝐸𝐟_𝑖̃𝐟^𝑇_𝑖̃∣𝐱_𝑖,𝑤_𝑗^𝑇𝐿^𝑇_𝑗=0. (19) Substituting (15) for 𝐿_𝑗 and making constraints on the diagonal of Ψ, we obtain 1Ψ=𝑁diag_𝑖,𝑗ℎ_𝑖𝑗𝐱_𝑖−𝐿_𝑗𝐸𝐟_𝑖̃∣𝐱_𝑖,𝑤_𝑗𝐱^𝑇_𝑖. (20) The prior probability 𝑝(𝑤_𝑗) should be proportional to the clinical diagnosis such that the estimation of the factor loadings and the factor scores can capture the latent factors of different disease groups. The proposed model also carries the same indeterminacy problem associated with factor patterns; that is there exist numerous orthogonal transformations to rotate the matrix of factor loadings without changing the maximum of 𝑄 [9]. Considering 𝐻 be any 𝐾×𝐾 orthogonal matrix, 𝐻𝐻^𝑇=𝐻^𝑇𝐻=𝐼. Equation (1) can be written 𝐱=_𝑗𝜋_𝑗𝝁_𝑗+𝐿_𝑗×𝐻𝐻^𝑇×𝐟∣𝑤_𝑗+𝐮=_𝑗𝜋_𝑗𝝁_𝑗+𝐿^∗_𝑗×𝐟^∗∣𝑤_𝑗+𝐮, (21) where 𝐿^∗_𝑗=𝐿_𝑗𝐻 and 𝐟^∗∣𝑤_𝑗=𝐻^𝑇×𝐟∣𝑤_𝑗. The assumption, 𝐟^∗∣𝑤_𝑗∼𝑁(0,𝐼), is kept. The covariance of 𝐱 is 𝐿^∗_𝑗(𝐿^∗_𝑗)^𝑇+Ψ=𝐿_𝑗𝐻𝐻^𝑇𝐿^𝑇_𝑗+Ψ=𝐿_𝑗𝐿^𝑇_𝑗+Ψ, which remains the same. Therefore, there are infinite equivalent solutions to satisfy the maximum of (5). Imposing reasonable constraints to identify a set of model parameters can make the factor loadings scientifically interpretable. A widely used approach for a simple factor structure is realized by setting some factor loadings to hypothetical values such as zeros.

The permutation and changing the sign of columns in the factor loading matrix with factor scores does not affect the model at all and the algorithm will yield the same solution. In order to realize consistent, interpretable, and comparable results, we suggest to recursively test all combinations to find the one of them that has the highest similarity among 𝑀 factor loading matrices so that we can find a coherent interpretation for different groups of subjects. Each pair of factor loading and factor scores can be multiplied by either +1 or −1. The 𝑀 sets of loadings has (2^𝐾)^𝑀 combinations. The possible permutation of the 𝑀 set of loadings is the factorial of 𝐾. The complexity of the recurrence is therefore 2^𝐾𝑀×(𝐾!)^𝑀. The problem can be formulated as a bipartite matching and the Hungarian algorithm can find the match in a lower complexity.

2.2. Data Description

The T1-weighted MR images of 416 subjects were downloaded from the Open Access Series of Imaging Studies [10], which is publically available for analysis. All the T1-weighted images were acquired on a 1.5-T Siemens Vision scanner. Among all 416 subjects, there are 316 normal subjects (average age: 45.09±23.90), 70 subjects who have been clinically diagnosed with very mild AD (average age: 76.21±7.19), and 30 are with moderate AD (average age: 78.03±6.91). The proportions of each type of subject are 𝝅=[75.96%,16.83%,7.21%]^𝑇. Multiple intrasession acquisitions provide extremely high signal-to-noise ratio, making the data amenable to our analysis. The available images were provided skull stripped, gain field corrected, and registered to the atlas space of Talairach and Tournoux [11] with a 12-parameter rigid affine transform. The resolution of the images is 176×208×176. The number of voxels, which is more than six million, is much larger than the number of subjects. We extracted the clinically and psychologically interested regions instead of processing whole voxels in the image. The subcortical structures are extracted by the segmentation method [12] which uses manually labeled image data as priori information for a Bayesian framework that utilizes the principles of the active shape and appearance models. The size of a subcortical region was calculated by multiplying the voxel size and the number of voxels in the region. Fifteen subcortical regions were successfully extracted.

According to a demographic study by the National Institute on Aging and Alzheimer’s Association based on the data collected in the Chicago Health and Aging Project, the prevalence of dementia among individuals aged 71 and older was 13.9%, and AD (Alzheimers disease) was 9.7% [13]. The study was based on a sample of 856 individuals. The 𝝅 was estimated to be [76.4%,13.9%,9.7%]^𝑇 which is close to the statistics in our empirical data. The data vector of each subject had fifteen dimensions, each corresponding to the volume size of a subcortical structure divided by the estimated total intracranial volume. The average size of all of the intracranial volume is 1480.5cm³. The intracranial volume is estimated by the linear registration from a manually measured intracranial volume of a standard brain to the individual brain [14]. The analysis of variance (ANOVA) of the data for each structure were calculated and shown in Table 1 and Figure 1. The smaller 𝑝 value indicates high probability of inequality of the structure size among the three groups.

Table 1 

The ANOVA results for the three groups in different subcortical structures.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

Figure 1 

The plots of means and standard deviations for the three groups in different subcortical structures.

We subtracted the mean from the data and used the remainder for analysis. Using the covariance matrix of the data to estimate the factor scores would cause that a few structures dominate the factor loadings; therefore, we divided each dimension by its standard deviation to compel each of them to have unit variance. After the algorithm converged, we used varimax rotation [15], which transforms the loadings into the space that maximizes the variance, to rotate the factor loadings. Given data, the expectation of its type was set to ℎ_𝑖1,ℎ_𝑖2,ℎ_𝑖3=[][][]1,0,0ifsubject𝑖belongstoNL0,1,0ifsubject𝑖belongstovAD0,0,1ifsubject𝑖belongstoAD. (22)

3. Results

Figure 3 shows the trend of the likelihood climbs as adding the number of factors in the analysis. In the scree plot in Figure 2, three eigenvalues of the covariance matrix of the whole dataset are greater than one and the cumulative percentage of variance from the largest three eigenvalues reaches 78%. Thus we set 𝐾=3 in this analysis.

Figure 2 

The scree plot for the ordered eigenvalues.

(a)

(b)

Figure 3 

The two curves record the approximate value of log likelihood for two methods.

The factor loadings for the three groups are shown in Figure 4, in which the vertical axis marks the fifteen regions. The vAD denotes the group of very mild AD. The log likelihood in (5) after the algorithm converges is −5475.814. The loading of structures has symmetric property and usually the right and left structures have similar loadings. Using the factor loadings to estimate 𝝅 and the expected group information given 𝐱 by (9), we obtain the adjusted and turned proportions as [82.89%,7.97%,9.14%]^𝑇. This may suggest the underlying variation among different groups of subject and need further investigation. Note that the reestimated proportion ℎ_𝑖𝑗 is not binary anymore.

Figure 4 

The factor loading matrices for the three groups.

We show the results of conventional factor analysis in Figure 5 as a comparison. The program run on the mild AD patients in the dataset cannot achieve reproducible results; therefore, the quantity of mild AD’s results in Figure 5 varies from time to time. The analysis for the AD group cannot converge, however the factor loadings are reproducible. The distance of whole factor loading matrices among the three groups for conventional factor analysis is 5.28 while the proposed method is 4.54. The correlation of the three-factor loading matrix estimated by conventional factor analysis methods is [𝐶₁₂,𝐶₁₃,𝐶₂₃]=[0.3879,0.1698,0.3633]. The correlation by proposed method is [0.5388,0.5564,0.4986]. Table 2 lists the 𝑝 values for all factors by the Kolmogorov-Smirnov test [16] on the factor scores against a Gaussian distribution. The test examines the difference between input distributions and a Gaussian distribution. The smaller the 𝑝 values, the more strongly the test rejects the Gaussian assumption. The algorithm tries to search for loads with normally distributed factor score, hence large 𝑝 indicates the factor fit well to Gaussian distribution.

Table 2 

The normality test of factor score by Kolmogorov-Smirnov test.

Figure 5 

The result of maximum likelihood method and rotated by varimax with Kaiser normalization. The results of healthy controls are not unique and unstable.

The means (centers) of the clusters are shown in Figure 6. The means are near the origin and include negative value because the data are standardized by the subtraction of the overall mean of the data in the preprocess stage. The yellow color in the first column indicates that healthy controls have larger sizes in subcortical structures, and the second and the third columns indicate that the patients have smaller sizes in different subcortical regions in general. The AD patient has very small thalamus, putamen, and hippocampus. The hippocampus is related to memory and learning. The putamen is a structure involved in the regulation of voluntary movement. The abnormal pallidum in Figure 4 can cause movement disorders. Figure 7 shows the associations of first factor scores with the score of minimental state examination (MMSE) by both methods.

Figure 6 

The cluster means of the three groups.

(a) Proposed method

(b) Traditional method

Figure 7 

Correlations of factor scores with the MMSE scores. The red color denotes healthy subjects; the green color denotes very mild AD patients; the blue color denotes the moderate AD patients. The regression lines of the three groups by the (a) proposed method and (b) conventional method are also plotted in the figures.

4. Conclusions

The proposed method finds closer and more correlated factor loadings than the conventional method because it considers the same error matrix for different groups of data. The result of conventional factor analysis having higher normality for AD patients than normal subjects is less convincing. Conventional factor analysis that decomposes the observed data together intermixes the latent factors. Taking the data apart will misseparate the noise. This work proposed using a mixture model of factor analysis method for neurodegenerative disease research by showing highly correlated factor loading across different groups of subjects and together with proper normality of the factor scores.

Acknowledgment

This work was supported by National Science Council under project NSC101-2811-M-001-082.