2.1 Linear latent variable models for functional connectivity: PCA and its extensions

In this section we outline the linear latent variable models employed in the unsupervised learning stage of the proposed framework. We begin by discussing principal component analysis (PCA), a well-established technique for dimensionality reduction [23]. The common derivation for PCA poses it as an optimization problem seeking to learn the linear projection which maximizes explained variance within the projected space [24]. However, PCA can also be derived as inference under a simple linear latent variable model, which posits that observations are generated as a linear projection from low-dimensional latent variables, [25]. When both observations and latent variables are taken to follow a multivariate Gaussian distributions we obtain the following generative model for observed data: (1) (2) where is a diagonal matrix and denotes measurement noise. Eqs (1) and (2) serve to highlight how PCA can be seen as a low-rank model for the covariance matrix; by marginalizing over latent variables we obtain: (3) implying that the loading matrix, W, captures low-rank covariance structure. Learning the associated loading matrix, W, proceeds via maximizing the log-likelihood over observations across all N subjects: (4) where Σ⁽ⁱ⁾ is as defined in Eq (3) and K⁽ⁱ⁾ denotes the sample covariance matrix for the ith subject. In the context of PCA, the maximization is performed subject to the constraint that W be orthonormal, (5) and a closed-form solution is obtained via eigendecomposition.

Following [21] it is possible to interpret each column of W as encoding functional networks or “eigenconnectivities”. While the loading matrix, W, is shared across all subjects, each diagonal entry of G⁽ⁱ⁾ denotes the extent to which the associated network is expressed in subject i. This allows us to study connectivity as being composed of various distinct networks, resulting in significant benefits from the perspective of interpretability. We can further unpack Eq (3) as follows (see also Fig 2 below): (6) where W_j denotes the jth column of W and we write to denote the jth diagonal entry of the matrix . As such, we may interpret each W_j as encoding the jth network and as a measure of activity within the corresponding network in the ith subject.

There exist several extensions to the model described in Eqs (1) and (2), the prime example being factor analysis which allows the variances in Eq (2) to vary across dimensions. Recently, several extensions have been proposed where constraints such as non-negativity are introduced with the goal of improving the interpretability of results [26, 27, 28]. The motivation behind such methods stems from the fact that interpreting and visualizing PCA-based networks becomes very challenging, particularly in high-dimensions. Challenges arise from the fact that each principal component will correspond to a weighted sum of BOLD activities across all observed regions. As such, it is often difficult to identify which regions are the principal contributors to a certain principal component (and hence functional network) without applying ad-hoc post analysis. Furthermore, it is possible that some entries in the principal components may be negative, which further complicates the interpretation from the perspective of functional connectivity analysis.

The aforementioned issues can be mitigated via the introduction of non-negativity constraints on the loading matrix, W. This ensures that each principal component corresponds only to a weighted positive sum of activity over all brain regions. As such, the principal component can be directly interpreted as the contribution of each region to each functional network. Furthermore, the introduction of non-negativity will often yield sparsity in the sense that many of the entries of the principal components will be exactly zero [27]. It follows that such sparsity further facilitates the interpretation of the corresponding networks. From an optimization perspective, the loading matrix is inferred by maximizing the original log-likelihood objective, with the additional non-negativity constraint: (7)

It is important to note that the orthonormality constraint has been dropped in Eq (7), making the associated optimization problem less challenging. However, the combination of non-negativity and orthonormality, as enforced in [29], leads to several desirable properties. First, the loading matrix W has at most one non-zero entry per row. This implies that we may interpret the columns of W as encoding membership to k non-overlapping networks or clusters. Another very important benefit of introducing non-negativity and orthonormality constraints is that the matrix W is uniquely defined and identifiable. This is not the case in standard factor analytic models, where W is only identifiable up to an arbitrary rotation [30, 25]. Given that throughout this work we will directly interpret the columns of the loading matrix, W, as encoding functional connectivity networks, the lack of identifiability in PCA and factor analysis models is a significant limitation. We refer to the model presented in [29] as Modular Hierarchical Analysis (MHA). The associated optimization problem therefore becomes: (8)

MHA can therefore been seen to address the two important limitations of traditional models such as PCA and factor analysis; first that the presence of negative values in the loading matrix complicates the interpretation of such matrices (addressed via the use of non-negativity constraints) and second is the fact that the latent variables are rotationally invariant (addressed via the further introduction of orthogonality). A further limitation of models such as PCA and factor analysis is that they implicitly assume latent variables must be uncorrelated. In many cases, especially when such models are applied on data relating to a cohort of subjects, such an assumption will not be valid, implying the associated generated models are misspecified. In contrast, MHA is able to identify and recover components even when they are uncorrelated. This is an important theoretical advantage, as MHA continues to enjoy the same identifiability properties even in the presence of correlated latent variables, and practical advantage, as we demonstrate in this work. Finally, we note that in the context of fMRI data, MHA corresponds to an intuitive generative model whereby latent variables capture the activity within each functional network. The optimization of Eqs (5), (7) and (8) is discussed in S1 Appendix. Furthermore, we provide both Python and R code to implement MHA in S1 Code.

Moreover, we note that model introduced by [26], termed Modular Connectivity Factorization (MCF), shares many similarities with MHA. In fact, both methods introduce non-negativity and orthonormality over the loading matrix, W. The fundamental difference, however, is that MCF is not associated with a linear latent variable model, and instead parameters are inferred as follows: (9) where Σ⁽ⁱ⁾ is defined as in Eq (6) and K⁽ⁱ⁾ is the empirical covariance for the ith subject. A related approach was also proposed by [31].

Finally, it is important to note that whilst identifiability can be obtained via the combination of non-negativity and orthonormality, as is the case with the MHA model, it can also be obtained by relaxing the assumed distribution over latent variables, as is the case with independent component analysis (ICA) models. Formally, ICA is also a linear latent variable model, however, latent variables are no longer assumed to follow a Gaussian distribution [32]. While the relaxation of the Gaussianity assumption complicates the associated optimization, which must now be solved using gradient descent methods and accounting for the presence of multiple local optima due to the non-convex objective function [33], ICA has been widely employed in the study of functional connectivity [34, 35]. Moreover, we note that the “spatial” version of ICA used in fMRI reverses the roles of latent variables and loadings, which means that it is actually looking at the non-Gaussianity or sparsity of what we call here the loadings, corresponding to spatial patterns. Fig 2 provides a visualization of the benefits obtained by introducing each of the aforementioned constraints. In particular, we note that it is the combination of non-negativity together with orthonormality which yields interpretable and identifiable networks. We empirically validate such claims by applying all of the aforementioned models to synthetic and real fMRI datasets below.

2.2 Predicting brain age using functional network activity

The previous section outlined the various flavours of latent variable models which can be employed in order to learn functional networks across a cohort of N subjects. The aforementioned models allow us to decompose observed functional connectivity patterns as a linear sum of networks encoded by the columns of the loading matrix, W. While the loading matrix is shared across all subjects (indicating the same networks are present across all subjects), the extent to which they contribute to the observed covariance of the ith subject is denoted by the diagonal entries of G⁽ⁱ⁾, as stated in Eq (6).

We now consider the task of predicting the biological brain age, a⁽ⁱ⁾, using inferred functional connectivity networks as features. In the interest of interpretability we limit ourselves to linear regression models of the form: (10) Recall that corresponds to the jth diagonal entry of the matrix G⁽ⁱ⁾. As such, the proposed models will essentially seek to predict the biological age of subjects by considering activity within each inferred functional network. In the case of the ith subject, the observed activity in network j is quantified by . In practice, we will seek to quantify the activity of various functional networks on unseen subjects, defined to be subjects whose data was not employed to estimate loading matrix, W. We note that due to the orthonormality of W, together with Eq (6), we may estimate for data from unseen subjects, denoted by i*, as follows: (11) We note that Eq (11) requires the observation noise, v^(i*). This is not a concern for all subjects whose data is employed during the unsupervised learning of the latent variables, as parameters v⁽ⁱ⁾ are inferred alongside loading matrix, W. However, the primary goal of this work is to build predictive models which can generalize to unseen subjects. In this context, an estimate of the observation noise, v^(i*), can be obtained as follows: (12)

Although the class of models considered in Eq (10) may be considered amongst the simplest supervised regression models, they yield several important benefits when seeking to understand both the estimated parameters as well as the contribution of each of the features. In particular, each β_j corresponds to the regression coefficient summarizing the (linear) relationship between the activity of the jth network and biological age, conditional on all remaining networks. As such, if certain regression coefficients are deemed to be insignificant, we may conclude that the associated network is invariant during healthy ageing.

2.3 Hyper-parameter selection

The proposed two-stage estimation framework requires the input of only one hyper-parameter: the dimensionality of latent variables k. In the context of PCA and factor analysis, this hyper-parameter directly corresponds to the number of principal components or factors inferred, and a wide literature exists for tuning such a parameter [23]. One of the advantages of the latent variable models presented in Section 2.1 is that they each correspond to probabilistic models whose likelihood can be directly evaluated. As such, a logical choice to tuning hyper-parameter k is to directly maximize the log-likelihood over held out data.

In order to effectively perform hyper-parameter tuning as well as quantify the generalization performance of the proposed method, data was split into training, validation and test datasets as follows:

• First, a subset of subjects were held out as test data. As such, we obtain two datasets: where S_train, S_test ⊂ {1, …, N} denote the non-overlapping sets of training and test subjects respectively. Recall N is the number of subjects present and we write to denote the n observations available for the ith subject.
• Training data is further split into training and validation datasets on a subject-by-subject basis.

Splitting the data in this manner allows for effective hyper-parameter tuning, using training and validation datasets, as well as for generalization performance to be measured using test dataset which corresponds to unseen subjects.

2.4 Experimental data

The data employed in this manuscript corresponds to resting-state fMRI data taken from three distinct open-access repositories. There were small variations in the resting-state functional MR image acquisition for each of the repositories considered: CamCAN [38], Human Connectome Project [37], and the ATR Wide Age Range [38]. The pre-processing employed on each dataset was as follows:

• CamCAN: This dataset was pre-processed by us. Data was motion corrected, spatially smoothed with a 5mm FWHM Gaussian kernel, registered into MNI152 standard space using FLIRT [39] via a skull-stripped high-resolution T1 image and resampled to 4x4x4mm voxel sizes. Each high resolution T1 image was segmented into grey and white matter and cerebrospinal fluid using SPM Dartel [40]. Mean timecourses for cerebrospinal fluid and white matter as well as 6 motion parameters were linearly filtered from each voxel to reduce non-neural noise.
• HCP: We used the pre-processed resting-state fMRI data from a random subset of healthy participants. Notably, the pipeline involved FIX ICA-based noise reduction process [40], to remove individual sources of physiological, non-physiological and motion related noise. Full details of the pre-processing pipeline can be found at https://www.humanconnectome.org/study/hcp-young-adult/document/extensively-processed-fmri-data-documentation.
• ATR: We used the preprocessed data. The pre-processing pipeline notably included regressing out the global grey matter signal as well as signals from cerebrospinal fluid and white matter, to remove sources of spurious variation

All three pre-processed fMRI datasets were subsequently processed as follows: a cortical parcellation based on resting-state functional connectivity analyses [42] was used to define 264 distinct 10mm diameter regions of interest (ROIs). The fMRI time course averaging across all voxels within each ROI was extracted. These 264 average time courses were then used in subsequent analyses. Full details are provided here https://bicr-resource.atr.jp/var/www/webapp/bicrresource/bicrresource/staticfiles/pdf/Methods.pdf.

3.1 Synthetic data experiments

In this section we evaluate the performance of the proposed two-stage estimation framework using synthetic data. To this end, we generate artificial data whose properties approximately match those which are frequently reported in fMRI studies. The objective is then to quantify which of the linear latent variable models presented in Section 2.1 are able to both robustly recover the associated loading matrix, W, as well as learn the relevant factors which serve as accurate predictors of brain age on unseen subjects. Synthetic data was then generated in order to satisfy Eqs (1), (2) and (10). This is achieved as follows:

• First, we randomly generated a factor loading matrix, , which satisfied the constraints of both non-negativity and orthonormality. The reason for introducing both constraints is that we will seek to quantify how reliably each latent variable model can recover W, and it is therefore imperative to ensure we generate W from an identifiable model (see discussion in Section 2.1). In order to achieve this a dense matrix, W, was sampled with each entry following a uniform distribution over the interval [0, 1]. Subsequently, for each row only the entry with the largest value was retained with all other entries set to zero. Finally, the norm of each column was set to one.
• Second, the factor loadings for the ith subject, , were randomly generated as follows: with all negative samples being discarded.
• The regression coefficients, , were drawn uniformly at random from the interval [0, 10].
• Finally, we are able to randomly generate observations and ages for each subject as follows: (13) (14) Recall that is a diagonal matrix consisting of entries .

We note that the choices for sampling distributions of both the factor loadings, g⁽ⁱ⁾, as well as the regression coefficients, β, are necessarily somewhat heuristic. However, care was taken to ensure the implied distributions over subject ages approximately matched the empirical distributions observed within the CamCAN repository.

We note that throughout experiments we consider the performance of each method whilst varying two distinct factors: the number of observations per subject, n, and the number of training subjects, N. Furthermore, throughout simulations we fix the dimensionality of observations to be p = 50 and the number latent factors to be k = 5.

Given artificial data generated as described above, we look to quantify the performance of each of the linear latent variable models using the following two metrics:

1. Accurate recovery of the loading matrix, W. This is quantified in terms of the squared error between the true loading matrix and the estimated loading matrix.
2. Accurate brain age prediction over unseen subjects. In line with other literature, this is quantified in terms of the mean absolute error between true and predicted brain ages [11, 8].

3.1.1 Synthetic data results.

We begin by considering the performance of each linear latent variable model as the number of observations per subject, n, increases for a fixed number of training subjects, N = 25. The results are presented in Fig 3. We note that both in terms of recovery of the loading matrix, W, as well as in terms predicting the ages over unseen subjects, the introduction of regularity constraints, be they in the form of non-negativity, orthonormality or non-Gaussianity or sparsity (as in ICA), leads to improvements.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3.

Simulation results for recovery of the true loading matrix (left panel) and prediction of brain age for unseen subjects (right panel) as the number of observations per subject, n, increases. We note that the introduction of regularity constraints (e.g., non-negativity or orthonormality) on the loading matrix leads to improvement in performance.

https://doi.org/10.1371/journal.pone.0232296.g003

We also study the performance of the various latent variable models when the number of training subjects, N, increases and the number of observations is fixed at n = 100 per subject. These results are presented in Fig 4. In terms of recovery of the loading matrix, W, we again observe that introducing regularity constraints leads to significant improvements. In terms of predictions over unseen subjects (as shown in the right panel of Fig 4), the improvements due to the introduction of regularity conditions begin to fade as the number of training subjects increases. In particular, beyond a certain number of training subjects (approximately 25 in the case of these experiments), the improvement in out-of-sample predictions begins to plateau.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4.

Simulation results for recovery of the true loading matrix (left panel) and prediction of brain age for unseen subjects (right panel) as the number of training subjects, N, increases. We note that the introduction of regularity constraints (e.g., non-negativity or orthonormality) on the loading matrix leads to improvement in performance.

https://doi.org/10.1371/journal.pone.0232296.g004

3.2 Resting-state fMRI data experiments

While the previous section presented results relating to synthetic data, here we present experimental results where the proposed two-step procedure is applied to three open-source resting-state fMRI datasets. The datasets considered correspond to the Cambridge Center for Ageing and Neuroscience (CamCAN) repository, the Human Connectome Project (HCP) repository, and the ATR Wide-Age-Range repository. The purpose of employing three distinct datasets is to effectively measure the generalization performance of the proposed approach on unseen data. As such, data from the HCP and Wide-Age-Range repositories was not employed during any of the model training and instead used exclusively as unseen test data. It is important to note that in addition to significant inter-subject variability [45], fMRI data also suffers from the presence of several other well-documented issues such as variable scanner performance or noise [46, 47, 48]. As such, validating the performance of the proposed brain age prediction models in this way will provide a more realistic measure of their generalization performance.

3.2.1 CamCAN repository results.

Resting-state fMRI data was collected from a total of 647 subjects from the CamCAN repository. Subject ages ranged from 18 to 88 years of age (average age of 54.31±18.56, 318 males and 329 females). The CamCAN dataset was employed as the principal dataset in the proposed two-step procedure, implying that it was employed to learn both the functional network structure in the unsupervised learning stage and the linear regression models in the supervised learning stage. As such, the data was split into training, validation and test subsets as described in Section 2.3.

Step 1: Unsupervised functional network inference. The first stage of the proposed framework involves the estimation of reproducible functional connectivity networks via the use of the various linear latent variable models discussed in Section 2.1. The number of functional networks inferred corresponds directly to the dimensionality of latent variables, which is determined by hyper-parameter k. As each linear latent variable model can be interpreted as a probabilistic model, we select hyper-parameter k by maximizing the log-likelihood over the validation dataset. This resulted in the choice of k = 5 when the loading matrix was restricted to be both non-negative and orthonormal, as proposed by [26] and [29]. While it is possible that the choice of hyper-parameter may vary across distinct latent variable models (e.g., for PCA or factor analysis), we choose to keep the choice of k fixed across all models as this facilitates model comparison and interpretation of results.

The left panel of Fig 5 visualizes the results when the MHA linear latent variable model was employed (Figures produced using the plot glass brain function from the nilearn python module [49]). We note that, as discussed in Section 2.1, the MHA linear latent variable model effectively clusters regions into sub-networks via the introduction of non-negativity and orthonormality constraints. As such, each plot in the left panel of Fig 5 visualizes spatially remote brain regions which have been clustered together, indicating that these regions share strong positive correlations. We note that these correlations (i.e., edges in a network) are omitted for clarity in Fig 5. The results demonstrate that the inferred networks are spatially homogeneous and symmetric across both hemispheres. Furthermore, many of the inferred networks correspond to widely reported networks and regions: network 1 captures the default model network (DMN) and network 2 overlaps with the salience network, while networks 3 and 4 correspond to a higher-level visual network and the somatomotor network respectively. For comparison, we include equivalent plots for all other latent variable models considered in visualized in Fig 6, presented in the Supplementary Material. We note that alternative methods, such as PCA, which did not enforce the combination of both non-negativity and orthonormality, yielded results which were visibly less clustered and more difficult to interpret.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5.

Left panel: Inferred networks as recovered when non-negativity and orthonormality constraints are introduced over the loading matrix, W. Networks are spatially consistent and symmetric. Right panel: visualization of network activities against subject age demonstrating (mostly negative) linear trends with healthy aging.

https://doi.org/10.1371/journal.pone.0232296.g005

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Inferred networks using alternative linear latent variable models.

In the case of models such as PCA and factor analysis, networks were obtained by thresholding entries of W so only non-negative entries considered.

https://doi.org/10.1371/journal.pone.0232296.g006

The right panel of Fig 5 visualizes the correlation between the activity of each network (as defined in Eq (11)) with the age of each subject. For networks 1-3 we observe a significant negative correlation between the activity and age, suggesting that ageing induces a drop in activity of such networks. These results are in line with related research on ageing induced differences in functional connectivity. In particular, the decrease in activity of the DMN (network 1), has been widely reported [19, 50, 51].

Step 2: Supervised training of brain age prediction models. Recall that the overall objective of the proposed framework was build interpretable models of biological brain age. To this end, the features recovered from linear latent variable models where employed as features in a linear regression framework to predict the brain age of each subject. In particular, the five distinct the linear latent variable models detailed in Section 2.1 where employed to learn reproducible sub-networks parameterized by a loading matrix, . The activity within each functional network, defined as in Eq (11), was subsequently employed as features to predict biological age using linear regression.

We note that the CamCAN repository, as well as HCP and ATR repositories, each contained over a hundred subjects each. This is in contrast to typical fMRI studies, where the sample size is often in the range of 20 to 30 subjects [52, 48]. Furthermore, recall that the goal of experiments presented are to quantify performance on unseen resting-state fMRI data with a view to providing an indication of how each of the linear latent variable models employed would perform in a typical fMRI study. As such, throughout the remainder of this section we report the performance, in terms of mean absolute error, over random subsets of 30 subjects from each repository. This corresponds to a form of bootstrapping, where we average results over a random sample of possible cohorts. In practice, we report results over 1000 random subsets of 30 subjects for each of the three repositories considered.

Fig 7 visualizes the mean absolute error on unseen test data for various choices of k ∈ {2, …, 10}. We note that the combination of linear regression with the use of non-negativity and orthonormality constraints, as advocated by both the MCF and MHA models, leads to competitive performance over a range of choices of k. In particular, such algorithms out-perform both non-negative PCA and PCA, suggesting that the introduction of such constraints serves to improve the predictive properties of the model. Moreover, we note that Fig 7 indicates the presence of a bias-variance trade-off that is often encountered in supervised learning whereby performance on unseen test data begins to deteriorate as the number of parameters (in our case k) increases beyond a certain value.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Mean Absolute Error (MAE) performance for a varying number of networks, as determined by k (x-axis), on unseen test data from CamCAN.

We note that the combination of non-negativity and othonormality (MHA and MCF) yields competitive results across a wide range of k.

https://doi.org/10.1371/journal.pone.0232296.g007

As mentioned previously, the choice of k = 5 was selected in by maximizing log-likelihood over a validation dataset (i.e., in an entirely unsupervised manner—data regarding subject ages was not considered). Fig 8 visualizes the performance on the unseen test dataset for the specific choice of k = 5, for all possible choices of linear latent variable models. The results indicate that as additional constraints are introduced to the loading matrix, the generalization capabilities of the models also improve. As such, MCF and MHA, which introduce the most stringent constraints corresponding to both non-negativity and orthonormality, obtain the best generalization performance. We also note that ICA is also competitive. Moreover, non-negative PCA, which relaxes the requirement for orthonormality, is the next most competitive latent variable model. Finally, PCA and factor analysis, which relax all the aforementioned constraints, obtain the worst generalization performance.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Mean Absolute Error (MAE) performance on unseen testing data from CamCAN repository when the dimensionality of latent variables is fixed to k = 5 (implying we infer 5 networks).

We note that as regularity constraints are introduced, in particular non-negativity and orthonormality, predictive performance improves.

https://doi.org/10.1371/journal.pone.0232296.g008

3.2.2 Transfer onto HCP and ATR Wide-Age-Range repositories.

The results of Section 3.2.1 provide a measure of performance, in terms mean absolute error in predicted brain age, within a large-scale resting-state fMRI dataset. However, it is widely accepted that in addition subject-specific noise, there are several other significant contributors to noise in fMRI data: these include issues related to scanner noise and frequency of acquisition of images [46, 47, 48]. As a result, in order to thoroughly verify the generalization performance of the proposed methods, we employ resting-state fMRI data from the HCP and ATR Wide-Age-Range repositories. We note that data from the aforementioned repositories was employed only for testing purposes, as such it was not employed to learn the network structure across subjects, nor to tune the parameters of the linear regression models. For a summary of the characteristics of HCP and ATR Wide-Age-Range datasets see Fig 9 and S1 Table in the Supplementary Material.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Histogram visualizing age distribution for each of the repositories employed.

We note that the CamCAN dataset has the widest range of all repositories considered, validating its use as a the primary dataset in our study.

https://doi.org/10.1371/journal.pone.0232296.g009

Prediction of biological age on both the HCP and ATR Wide-Age-Range repositories was performed as follows: First, the loading matrix, was employed to obtain estimated activity within each network, as detailed in Eqs (11) and (12). Subsequently, predictions of biological age were obtained using Eq (10). At each stage both and are the parameters inferred using the CamCAN dataset (i.e., there was no fine-tuning of parameters). As a result, performance on both HCP and ATR Wide-Age-Range datasets provide a robust measure of generalization performance to entirely unseen data.

Results on the HCP data are provided in Fig 10. As expected, the mean absolute errors are larger for each of the distinct latent variable models when compared to the results of on the CamCAN dataset (Fig 8), which will be partially the result of varying scanner noise and image acquisition properties. Importantly we note that, as with the CamCAN dataset, there once again a relationship between the introduction of additional constraints (in the form of non-negativity, orthonormality or non-Gaussianity) and generalization performance. As before, methods such as PCA and factor analysis which do not introduce any constraints had the weakest performance as well as the largest drop in performance.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Mean Absolute Eerror (MAE) performance on unseen data from HCP repository.

Results are broadly consistent with performance on the CamCAN data, indicating good generalization. We note that the introduction of non-negativity or orthogonality constraints leads to improved generalization. The number of functional networks was k = 5.

https://doi.org/10.1371/journal.pone.0232296.g010

The HCP results presented above serve to partially validate the predictive models trained using the CamCAN dataset. However, one significant limitation of the HCP dataset is that subject ages only range from 22 to 37 years of age. This is particularly relevant in the context of brain age biomarkers, as many neurodegenerative diseases of interest will be associated with advanced ages. As a result, we further validated the generalization capabilities of the proposed brain age prediction models on the ATR Wide-Age-Range dataset, which had subjects ranging from 20 to 70 years of age. Results, presented in Fig 11 are consistent with results on the CamCAN and HCP datasets, again indicating that the introduction of constraints non-negativity and orthonormality constraints improves generalization performance.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. Mean Absolute Error (MAE) performance on unseen data from ATR Wide-Age-Range repository.

Results are broadly consistent with performance on the CamCAN data, indicating good generalization. Further, as with the HCP data, we note that the introduction of non-negativity or orthogonality constraints leads to improved generalization. The number of functional networks considered was k = 5.

https://doi.org/10.1371/journal.pone.0232296.g011

3.3 Extension to non-independent latent variable models

The results presented above employ linear latent variable models where the inferred latents are assumed to be independent. This is clearly stated in the generative model considered in Eq (1) where the covariance of latent variables, G⁽ⁱ⁾, is assumed to be diagonal. Note that in the case of PCA, factor analysis and MHA, since latent variables are assumed to be multivariate Gaussian, the fact the covariance is diagonal implies the latent variables are independent. However, such an assumption will often fail in practice, implying that the empirical covariance structure over latent variables will not be diagonal. In this section we seek to exploit this by directly introducing the off-diagonal entries of the latent variable covariances, G⁽ⁱ⁾, as features in our linear regression models for biological age. As such, whilst Eq (10) considered a linear model where only the diagonal entries of each G⁽ⁱ⁾ were employed to predict biological ages of each subject, we now consider linear regression models of the following form: (15) Note that in Eq (15) we employ the full upper triangular entries of the covariance matrix as features. This is equivalent to vectorizing the covariance matrix and removing duplicate entries due to symmetry. As such, whilst k features were employed in Eq (10), we now consider a linear models with features; many of which will seek to predict the biological age of individuals based on the off diagonal entries of each G⁽ⁱ⁾. It is important to note that the model presented in Eq (10) is a special case of Eq (15).

As in Section 3.2, we proceed in a two-stage approach whereby we first estimate the loading matrices for the various linear latent variable models employed and subsequently train linear regression models using the full vectorized covariance matrix as features.

Fig 12 visualizes the MAE error on unseen test data as a function of the dimensionality of latent variables, k. We note that for all choices of k the reported errors are smaller than those reported in Fig 7. This provides empirical evidence that the off-diagonal entries of the latent variable covariances are discriminative features for brain age prediction, and therefore can be seen as evidence that models which assume diagonal covariance structure over latents are misspecified. Fig 13 provides further visualizations in the case where k = 5. We note that the MHA model performs competitively, this is to be expected as this model directly accommodates the possibility of non-independent latent variables [29]. Moreover, we note that MHA performs particularly well when the number of networks is small (when dimension of latent variables, k, is less than or equal to 5), which is useful when we wish to prioritize the interpretability of results. Finally, the performance of various methods, as depicted in Fig 12, shows similar trends as in Fig 7; there is once again a bias-variance trade-off associated with the choice of k and the introduction of non-negativity or non-Gaussianity constraints (as in MCF or ICA) leads to improved generalization performance. Finally, whilst Fig 12 only shows generalization performance to unseen subjects from the CamCAN cohort, we also present results for generalization performance to brain age prediction on the HCP and ATR Wide-Age-Range datasets in Figs 14 and 15 of the Supplementary Material.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 12. Mean Absolute Error (MAE) performance for a varying number of networks, as determined by k (x-axis), on unseen test data from CamCAN when latent variables are no longer assumed to have an isotropic covariance structure and the full vectorized covariance is employed as features in the linear regression models.

We note that MHA is able to directly accommodate such a scenario and hence is competitive for all choices of latent variable dimension, k.

https://doi.org/10.1371/journal.pone.0232296.g012

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 13. Mean absolute error (MAE) performance on unseen testing data from CamCAN repository when the dimensionality of latent variables is fixed to k = 5 (implying we infer 5 networks).

Note that latent variables are no longer assumed to have an isotropic covariance structure and the full vectorized covariance is employed as features in the linear regression models.

https://doi.org/10.1371/journal.pone.0232296.g013

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 14. Mean Absolute Error (MAE) performance on unseen data from HCP repository.

Results are broadly consistent with performance on the CamCAN data, indicating good generalization. We note that the introduction of non-negativity or orthogonality constraints leads to improved generalization. The number of functional networks was k = 5.

https://doi.org/10.1371/journal.pone.0232296.g014

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 15. Mean absolute error (MAE) performance on unseen data from ATR Wide-Age-Range repository.

Results are broadly consistent with performance on the CamCAN data, indicating good generalization. Further, as with the HCP data, we note that the introduction of non-negativity or orthogonality constraints leads to improved generalization. The number of functional networks considered was k = 5.

https://doi.org/10.1371/journal.pone.0232296.g015