Soodeh Moallemian is affiliated with GIGA-CRC Human Imaging at the University of Liège, Belgium, and the Center for Molecular & Behavioral Neuroscience at Rutgers University–Newark in the United States. Her academic work is positioned at the intersection of human neuroimaging, molecular neuroscience, and behavioral research, reflecting an interdisciplinary approach to brain science. Through her involvement in the referenced preprint, she contributes expertise in imaging methodologies and neuroscience-informed analysis relevant to the study’s clinical or translational focus. Her dual institutional affiliations underscore engagement with collaborative research spanning European and U.S.–based neuroscience centers.
Soodeh Moallemian12, Christine Bastin1, Martina F. Callaghan3†* and Christophe Phillips1*
1GIGA-CRC Human Imaging, University of Liège, Belgium
2Center for Molecular & Behavioral Neuroscience, Rutgers University–Newark, Newark, NJ 07102, United States of America
3Wellcome Centre for Human Neuroimaging, UCL Queen Square Institute of Neurology, University College London, UK
†Corresponding author. m.callaghan{at}ucl.ac.uk
medRxiv preprint DOI: https://doi.org/10.1101/2023.10.19.23297253
Posted: November 07, 2025, Version 2
Copyright: * Authors equally contributed.
Abstract
This study applied multivariate ANOVA to investigate age-related microstructural changes in the brain tissues driven primarily by myelin, iron, and water content, as observed in MRI (semi-)quantitative R1, R2*, MTsat and PD maps. This is effectively a re-analysis of the data analyzed in a univariate way by Callaghan et al., 2014. Voxel-wise analyses were performed on gray matter (GM) and white matter (WM), in addition to region of interest (ROI) analyses. The multivariate approach identified brain regions showing coordinated alterations in multiple tissue properties and demonstrated bidirectional correlations between age and all examined modalities in various brain regions, including the caudate nucleus, putamen, insula, cerebellum, lingual gyri, hippocampus, and olfactory bulb. The multivariate model was more sensitive than univariate analyses, as evidenced by detecting a larger number of significant voxels within clusters in the supplementary motor area, frontal cortex, hippocampus, amygdala, occipital cortex, and cerebellum bilaterally. The examination of normalized, smoothed, and z-transformed maps within the ROIs revealed concurrent age-dependent alterations in myelin, iron, and water content. These findings contribute to our understanding of age-related brain differences and provide insights into the underlying mechanisms of aging. The study emphasizes the importance of multivariate analysis for detecting subtle microstructural changes associated with aging when dealing with multiple quantitative MRI parameter maps.
1. Introduction
Aging is an inevitable part of our lifecycle associated with the physical deterioration of different organs. Various hallmarks of aging have been identified over the past years that associate with neurodegenerative pathological changes in the brain, making age a primary risk factor for most neurodegenerative diseases, including Alzheimer’s disease (AD), Parkinson’s disease (PD), and frontotemporal lobar dementia (FTD) (Azam et al., 2021; Jeremic et al., 2021).
Quantitative MR imaging (qMRI) enables us to extract sensitive and specific information about the microstructural properties of the brain tissue in-vivo, such as axon, myelin, iron, and water concentration (Weiskopf et al., 2021). The estimation of (semi-)quantitative metrics normally includes effective longitudinal relaxation rate (R1), which is sensitive to iron, myelin and water content, transverse relaxation rates (R2*), which is primarily sensitive to iron, proton density (PD) indicative of free water content, and magnetization transfer saturation (MTsat) associated with macromolecular content, predominantly myelin (Tabelow et al., 2019a; Taubert et al., 2020a).
Many studies of aging and neurodegenerative diseases focus on alterations in the nervous system, such as (de-)myelination or iron accumulation (Moallemian et al., 2023; Peters, 2002; Tian et al., 2022). Callaghan et al. (2014) investigated age-related differences of biologically relevant in-vivo measures over the course of normal aging using quantitative multiparameter mapping (MPM) and showed significant demyelination in white matter (WM), concurrent with an increase in iron levels in the basal ganglia, red nucleus, and extensive cortical regions, but decreases along the superior occipitofrontal fascicle and optic radiation (Callaghan et al., 2014). Steiger et al. (2016) investigated the difference in iron and myelin levels between two groups of young and older participants using MPM qMRI, and showed that age-related higher levels of iron are accompanied by a negative correlation of iron and myelin in the ventral striatum (Steiger et al., 2016). Although demyelination primarily affects the WM of the brain, recent research shows that it can also occur in gray matter (GM), which is made up of the cell bodies of neurons. Studies have shown that gray matter myelination can occur during development, particularly in the prefrontal cortex, and may continue throughout late adulthood in response to learning and experience (Fields, 2008; Timmler and Simons, 2019). Khattar and colleagues assessed the association of myelination and iron accumulation in the brain of a cohort of 21-94 year-old healthy controls and found a negative correlation between whole brain myelin water fraction and iron content in most brain regions; they also highlighted that the myelination continues until middle age in the brain (Khattar et al., 2021). Another study reported robust evidence for spatial overlap between volume, myelination, and iron decomposition changes in aging that affect predominantly motor and executive networks under a modified normal probability curve approach from the Permutation Analysis of Linear Models (PALM) toolbox (Taubert et al., 2020b; Winkler et al., 2016, 2014).
In this technical note, we re-analyze the data from (Callaghan et al., 2014), now available as an open dataset on OpenNeuro (doi:10.18112/openneuro.ds005851.v1.0.0) with a multivariate statistical modeling approach as implemented in the MSPM toolbox (Gyger et al., 2021), available on GitHub (https://github.com/LREN-CHUV/MSPM), which allows a single inference over a combination of modalities. Our aim is to assess the advantages of such an approach, in terms of sensitivity and specificity, over the conventional multiple univariate analysis.
2. Method
2.1. Participants and data preprocessing
We took advantage of the processed data from (Callaghan et al., 2014; Callaghan and Phillips, 2025; Phillips et al., 2025) which include 138 healthy participants aged 19-75 years (35.5% male, mean = 46.64, s.d = 21). Quantitative multiparameter maps (R1, R2*, PD, and MTsat) were reconstructed with the VBQ toolbox, a preliminary version of the hMRI toolbox (Draganski et al., 2011; Tabelow et al., 2019b). Processing steps included segmentation using the “unified segmentation” approach (Ashburner and Friston, 2005), and diffeomorphic morphing to MNI space using DARTEL (Ashburner, 2007). Tissue-weighted smoothing (for GM and WM separately) with a 3mm FHWM isotropic kernel was applied to account for residual misalignment while preserving the quantitative nature of the data. These steps produced 8 smoothed normalized qMRI maps, one for each parameter (R1, R2*, PD, and MTsat) and tissue class (GM and WM). Finally, group-level GM and WM masks were created to be further used as exclusive masks in the statistical analysis. For full details see (Callaghan et al., 2014; Callaghan and Phillips, 2025).
2.2. Univariate GLM analyses
In this re-analysis, as recommended for MSPM (Gyger et al., 2021), the 8 resulting sets of qMRI maps were z-transformed per modality and across subjects to ensure comparability between maps, utilizing the mean and variance over each voxel. This procedure ensured comparability of different modalities for our multivariate analysis. The univariate general linear models were performed on each set of z-transformed maps using age, gender, total intracranial volume (TIV), and scanner as regressors.
2.3. Multivariate GLM (mGLM)
While a standard multivariate regression model is useful for predicting the values of the dependent variables, it doesn’t directly address the question of how the entire set of dependent variables relates to the entire set of independent variables. In the context of this study, it is well-established that the brain undergoes simultaneous changes as we age. Therefore, a single independent variable may not adequately predict a single dependent variable; instead, the relationship might be more complex, where a combination of predictors may explain variance across multiple outcomes. Hence, a multivariate approach appears as a method of choice.
Here, we describe the multivariate GLM (mGLM) and its principles, including model and inference, as implemented in the MSPM toolbox are described here after (Gyger et al., 2021). The multivariate GLM (mGLM) is specified using the design matrices of the 4 univariate models and the multivariate observations at each voxel as Y = XB + E, where Y138×4 = [Y1, Y2, Y3, Y4] is the multi-modal data matrix, each row of Y represents one participant, and each column represents one map; and X138×5 = [X1, X2, X3, X4, X5] is the design matrix, the first column represents the mean over subjects, the rest of the columns represent mean-centered regressors (age, TIV, gender, scanner) respectively. Thus, B is a 5 × 4 matrix of regression coefficients estimated using an ordinary least-square method, as B^ = (XTX)−1XTX. And Ê = Y − Ŷ = Y − XB^ is the residual matrix of size 138 × 4. Importantly, the residuals are estimated on a per-voxel basis that allows for a straightforward determination of a unique covariance structure for each voxel. This feature is a significant advantage of mass multivariate approaches when dealing with dependent neuroimaging data. However, it is important to note that in this framework, the assumption of normality for the residuals implies that the covariance structure is assumed to be the same across different groups or conditions. There is an assumed degree of correlation between the columns of Y, this correlation is expressed by estimation of variance-covariance matrix 
, where n is the number of subjects and k is number of covariates.
The total sum-of-square-and-cross-product (SSCP) matrix of the model is SSCPTotal = Y’Y − Ny–y–’ and can be partitioned into the regression and residual SSCP:
Hypothesis testing
Hypothesis testing in mGLM, relies on testing the linear contrast CBL = 0. This extension of the univariate scheme combines standard hypotheses on the rows of matrix B, i.e. on the regressor, represented by matrix C, with hypotheses on the columns of B, i.e. on the different modalities, represented by matrix L. In the context of multivariate ANOVA (MANOVA) models, contrasts of main effects and interactions are formulated by setting L = It, the t × t identity matrix, as the dependent variables in multimodal neuroimaging applications are not assumed to be directly proportional. This method is the most suitable for conducting hypothesis testing on multimodal neuroimaging applications.
The SSCP matrices are defined as:
And an SSCP matrix associated with the residuals:
Test statistics in mGLM
The global significance of the GLM, i.e. H0: B = 0, one can simply be estimated as:
Then we can calculate the eigenvalues of the resulting matrix by solving to get the m = 4 eigenvalues.
In hypothesis testing using the multivariate GLM, there are four standard test statistics available, which can be constructed based on the calculation of the SSCP matrices: Pillai’s trace (Pillai, 1955), Wilks’ lambda (Wilks, 1932), Hotelling-Lawley trace (Hotelling, 1951), and Roy’s largest root (Roy, 1945). The Wilks Lambda statistic can be calculated based on the calculated m eigenvalues as . This has an approximate F distribution with degrees of freedoms a and b:
Where a = lq and b = rt − 2u.
If the minimum value between l = 4 and q = 1 is ≤ 2, the distribution is exactly F.
The hypothesis H0: CBL = 0 can be tested to assess the different potential co-occurrence of change between the modalities. To test the joint effect on the 4 quantitative maps in a specific tissue type, L4×4 is defined as an identity matrix corresponding to the number of dependent variables (semi-quantitative maps). As explained before, each column of L will perform a univariate analysis on each column of B. Here, C was defined as [0 1 0 0 0] to only see the correlation between age and maps. In this case, Λ has an exact F distribution with a = 4 and b = 130 degrees of freedom.
These matrices are generalizations of the numerator and denominator sums-of-squares from the univariate GLM hypothesis-testing approach. When L is an identity matrix, the main diagonal of SSCPhyppo contains the sums of squares for the hypothesis in C as applied to the estimated parameters for each dependent variable separately. And the SSCPresidual matrix is an unscaled form of the estimated covariance matrix Σ^.
Construction of the test statistics rely on some linear combination of m eigenvalues (λ1, …, λq) of SSCP−1 SSCPhyppo. Here, we will only focus on Wilks’ lambda test statistics. It quantifies the proportion of variance not accounted for by the hypothesis compared to the total variance in the data.
Canonical Correlation Analysis
Canonical vectors are calculated under the assumption that matrix L involves multiple dependent variables (l > 1), to extract the contribution of each dependent variable to the test statistics Λ. This contribution corresponds to the eigenvectors of the eigen decomposition of SSCP−1 SSCPhyppo(Tabachnick and Fidell, 2007). Once linear combinations of canonical variates for the dependent (U1 = a1y1 + a2y2 + a3y3 + a4y4) and independent (V1 = b1x1 + b2x2 + b3x3 + b4x4 + b5x5) variables are formed, the first canonical correlation is calculated by solving for λ1. After finding the first pair of canonical variates, the process continues. The analysis finds a second pair of variates (U2and V2) that are also maximally correlated with each other, but with the added constraint that they are uncorrelated with the first pair of variates. This process continues until a certain number of canonical variate pairs are extracted. The number of possible pairs is equal to the number of variables in the smaller of the two sets (in our case, 4).
2.4. Statistical inference and multiple comparison
All statistical analyses were performed on standardized z-transformed data. To assess the correlation between the individual tissue property maps and age, an F-test was performed on each map at voxel-level over GM and WM separately. All univariate statistical analyses were performed under the general linear model framework in SPM12, considering two p-values, .05 and .0125, family-wise error rate corrected (FWER) thresholds. The latter threshold accounts for the fact that for both tissue classes, 4 similar inferences are performed (one per parametric map); thus, applying a Bonferroni correction, the applied threshold is divided by 4, i.e., .0125= .05/4. The former threshold, on the other hand, aims at replicating the previously published results (Callaghan et al., 2014).
With the multivariate GLM, one could simply rely on a .05 FWER threshold, as only one inference is performed for all 4 maps together.
2.5. ROI analyses
The age-related parameter differences in selected regions, including caudate, cerebellum, hippocampus, middle frontal gyrus, pallidum, putamen, superior motor cortex, heschl and precentral gyri, and thalamus, will be further described. ROIs were selected based on identified significant areas in our mGLM results in (Callaghan et al., 2014). The ROIs were defined according to the Neuro-morphometrics AAL3 atlas (Rolls et al., 2020). As this atlas contains multiple subregions of thalamus and cerebellum, to create their ROI masks, we took the union of all thalamus and cerebellum subregions. For each participant and semi-quantitative map, the measure from each voxel in an ROI was extracted. The regressors of no interest (gender, TIV and scanner) were regressed out, and for each selected ROI, the relation between age and median value across subjects will be examined post hoc.
3. Results
The analyses were conducted on both gray matter (GM) and white matter (WM). Consequently, explicit non-overlapping masks for GM and WM were applied to each analysis. Since each parametric map has its specific unit, such as Hertz for R1 and R2* maps and p.u. for MTsat maps, their values cannot be directly compared.
3.1. uGLM vs mGLM: Voxel level analyses
The individual GM and WM analyses on R1, R2*, PD, and MTsat maps at a corrected threshold of p<.05 FWER, concur with those in (Callaghan et al., 2014). The statistical parametric maps for age-related changes in the microstructure of the brain in GM and WM thresholded at a standard p<.05 FWER are depicted in Figure S1 and Figure S2 of the supplementary data. The statistical parametric maps for the same analyses thresholded at a corrected p<.0125 FWER are presented in Figures Figure 1 and Figure 2 for GM and WM, respectively.
Figure 1.
Statistical parametric maps of uGLMs in GM; showing all the voxels with significant correlation with age, as detected by uGLMs for PD, MTsat, R1, and R2* maps. The F-tests were thresholded at p<.0125 FWER corrected at voxel-level. The SPMs were overlayed on the mean MTsat map for the cohort in the MNI space. Abbreviation: GLM, general linear model; uGLM, univariate GLM; GM, gray matter; FWER, family-wise error rate; SPM, statistical parametric map.
Figure 2.
Statistical parametric maps of uGLMs in WM; showing all the voxels with significant correlation with age, as detected by uGLMs for PD, MTsat, R1, and R2* maps. The F-tests were thresholded at p<.0125 FWER corrected at voxel-level. The SPMs were overlayed on the mean MTsat map for the cohort in the MNI space. Abbreviation: GLM, general linear model, uGLM, univariate GLM, mGLM, multivariate GLM, WM, white matter, FWER, family-wise error rate, SPM, statistical parametric map.
Voxel-wise mGLM results presented in Figure 3 indicate bidirectional correlation between all modalities and age at p<.05 FWER corrected. The correlation is observed bilaterally in caudate nucleus, putamen, insula, cerebellum, lingual gyri, hippocampus, and olfactory bulb which matches the regions observed with the uGLMs.
Figure 3.
Statistical parametric maps of mGLMs in GM and WM; showing all the voxels with significant correlation with age, as detected by the multivariate model. The F-tests were thresholded at p<.05 FWER corrected at voxel-level. The SPMs were overlayed on the mean MTsat map for the cohort in the MNI space. Abbreviation: GLM, general linear model; mGLM, multivariate GLM; WM, white matter, FWER, family-wise error rate, SPM, statistical parametric map.
Table 1 provides a summary of the F-tests results (thresholded at voxel-level p<.05 FWER) from both the univariate general linear models (uGLMs) and mGLM for all maps within the two tissue classes. Additionally, Table 2 presents the summary of the uGLM F-tests results with Bonferroni-adjusted threshold. Comparing the spatial extent of clusters of significant voxels between uGLMs and mGLM in Tables 1 and 2 reveals that the multivariate model identifies a larger number of significant voxels compared to the union of all uGLMs with both p<.05 FWER and p<.0125 Bonferroni corrections.
Table 1.Summary statistics for significant voxels in uGLMs and mGLM.
“United” rows show the union of significant voxels in SPMs for all modalities. Univariate GLMs were thresholded at p<0.05 FWER corrected per tissue class (GM or WM). Abbreviation: GLM, general linear model; uGLM, univariate GLM; mGLM, multivariate GLM; GM, gray matter; WM, white matter.
Table 2.Summary statistics for significant voxels in uGLMs and mGLM.
“United” rows show the union of significant voxels in SPMs for all modalities. Univariate GLMs were thresholded at p<0.0125=.05/4 FWER corrected, to account for the multiplicity of maps tested (4) per tissue class (GM or WM). The mGLM was thresholded at p<0.05 FWER corrected. Abbreviation: GLM, general linear model, uGLM, univariate GLM, mGLM, multivariate GLM, SPM, statistical parametric maps, GM, gray matter, WM, white matter.
To illustrate the voxels affected by age at a threshold of p<.0125 and p<.05 in at least one of the (semi-)quantitative maps, the union of all statistical parametric maps derived from uGLMs in GM was binarized and depicted in Figure 4.
Figure 4.
The union of all uGLM statistical parametric masks in GM. On the left side, the F-tests for the uGLMs were thresholded at at p<.05 FWER corrected at voxel-level. On the right side, the F-tests for the uGLMs were thresholded at P<.0125 after correction for FWER at voxel-level. The masks were overlayed on the mean MTsat map for the cohort in the MNI space. The right side of the brain is illustrated here. Abbreviation: GLM, general linear model; uGLM, univariate GLM; GM, gray matter; FWER, family-wise error rate.
To quantitatively and visually compare the maps of significant voxels in GM from the uGLMs and mGLM, we built the union of binarized statistical maps from the uGLMs, thresholded at p<.0125 FWER, and a binarized statistical map from the mGLM, thresholded at p<.05 FWER. In Figure 5, we then show in red the voxels that are uniquely detected by the mGLM, and in green, those uniquely detected by any of the uGLMs.
Figure 5.mGLM vs multiple uGLMs in GM.
This figure shows the voxels detected either in the multivariate GLM model (thresholded at P<.05 FWER corrected at voxel-level in green), or in any of the univariate GLMs (thresholded at P<.0125 after correction for FWER at voxel-level in red). The mask is overlaid on the MNI152 template image in MRIcroGL. Abbreviation: GLM, general linear model; uGLM, univariate GLM; mGLM, multivariate GLM; GM, gray matter; FWER, family-wise error rate.
Among the models used, only the mGLM detected an age effect in certain regions, including portions of the superior medial frontal lobe, supplementary motor area, paracentral lobule, middle and anterior cingulum, parts of the precuneus, cuneus, calcarine, lingual gyrus, cerebellum, hippocampus, and parahippocampal gyrus bilaterally, as well as the left fusiform gyrus.
3.2. ROI analyses
We decided to perform Region of interest (ROI) analyses on the regions of interest listed in (Callaghan et al., 2014) to replicate previous findings. Moreover, we compared the canonical vectors withing the listed ROIs to investigate the contribution of each semi-quantitative map in those critical regions.
Median values
The median values of normalized, smoothed, and z-transformed R1, R2*, MTsat, and PD maps within the Putamen and Hippocampus ROIs with respect to age are illustrated in Figure 6 and Figure 7. ROI-based regression analysis for the adjusted medians and age, within each cluster, is also provided, along with the respective linear age dependence observed from the uGLM analyses in the selected regions.
Figure 6.
Median voxel values bilaterally in the Putamen. The red lines depict the linear model fit. These data are shown for illustration purposes only and were not used for any additional analyses.
Figure 7.
Median voxel values within the hippocampi. The red lines depict the linear model fit. These data are shown for illustration purposes only and were not used for any additional analyses.
From the ROI-based univariate regression analysis, we observe a bilateral decrease in the normalized MTsat and PD values in GM with respect to aging. While the PD and MTsat median values decrease with age in all regions of the brain, median R2* values show an increase in most regions of the brain except for the thalamus, where there is a significant negative age-related variation. These results are consistent with previously published results for the mean values in the thalamus for R2* (Taubert et al., 2020) but were not present in the uGLM results. R1 median signals show a weak positive correlation with age. The alteration in median values as a function of age concurs with the magnitude of the associated canonical vector. The bivariate correlation analysis indicates the strongest correlations for PD in most selected regions, see Table 3.
Table 3.Pearson Partial Correlations on the median voxel values within different regions of interest.
Canonical vectors
To further investigate the impact of aging on myelination, iron content, and water concentration in various brain regions, quantitative MR parameters were selected from bilateral regions including the putamen, thalamus, hippocampus, cerebellum, caudate, middle frontal gyrus, precentral gyrus, heschl gyrus, supplementary motor area, caudate, and pallidum (Callaghan et al., 2014; Darnai et al., 2017; Tian et al., 2022; Wang et al., 2020).
For each ROI, the canonical vectors corresponding to each quantitative map within the gray matter (GM) were estimated, as shown in the upper panel of Figure 8. It is important to note that the canonical vectors displayed in the lower part of Figure 8 represent the contribution of each map in the peak voxel of the respective ROI and do not represent the contribution factor for the entire ROI. Additionally, the direction of these vectors cannot be interpreted as they are derived from F-tests; therefore, we have displayed the absolute weights. As illustrated, MTsat signals exhibit the highest contribution, while R1 and PD signals contribute the least in the mGLM as indicated by the canonical vectors.
Figure 8.
Canonical vectors for different modalities from the mGLM model, representing the contribution of each modality in each voxel. The color bars show the vector sizes [0,1]. The vectors correspond to the peak voxels at the selected ROIs. Detailed vector sizes are reported in Table sup-1 of the supplementary data. Abbreviation: mGLM, multivariate general linear model; ROI, region of interest.
4. Discussion
The multivariate approach used in this study to investigate age-related changes in the microstructural tissue properties of the brain, incorporating image-derived quantitative maps for myelin, iron, and free water content, enables the identification of regions that are influenced by the simultaneous occurrence of various parameter differences. By considering multiple quantitative maps simultaneously, the multivariate approach provides a comprehensive description of how these different tissue properties interact and contribute to age-related differences. This method enables the identification of specific brain regions that exhibit coordinated alterations in myelin, iron, and water content, offering insights into the underlying mechanisms of aging.
The observed association with age in the quantitative MR parameters align with findings from ex vivo histologic studies and demonstrate a high specificity for tissue properties, including myelin content, iron content, and free water content. Using voxel-wise analysis with the multivariate GLM (mGLM), a bidirectional correlation between age and all the examined modalities was observed bilaterally in multiple brain regions. These regions encompassed the caudate nucleus, putamen, insula, cerebellum, lingual gyri, hippocampus, and olfactory bulb. Importantly, the multivariate approach demonstrated advantages over univariate analyses that focus on individual tissue parameters separately.
While examining individual tissue properties in isolation may provide insights into specific aspects of brain aging, the multivariate model revealed large clusters in the brain that could not be detected by analyzing each property individually. This indicates that the combined examination of multiple tissue properties enables the detection of additional regions associated with aging despite the presence of opposite changes in different properties. As observed in Table 2, mGLM showed improved sensitivity compared to the individual univariate GLMs (uGLMs) by detecting a larger number of significant voxels within clusters that cover the supplementary motor area, frontal cortex, hippocampus, amygdala, occipital cortex, and cerebellum bilaterally. This finding suggests that mGLM is a more effective and sensitive technique for detecting age-related differences in the brain, which should extend to other qMRI studies.
In this study, the application of a multivariate model such as MANOVA proves advantageous due to the well-established correlations among brain tissue characteristics and the interrelated nature of quantitative map values. By accounting for this inherent correlation, MANOVA effectively reduces potential biases in the results that could arise from using multiple ANOVAs. MANOVA’s ability to consider the complex relationships between multiple variables allows it to detect effects that might be smaller than those detectable by ANOVA, providing a more comprehensive understanding of the data.
Additionally, independent covariates can affect the relationship between the brain microstructure characteristics rather than influencing only a single variable; such patterns cannot be detected by the univariate models. Moreover, MANOVA offers a convenient means to manage the family-wise error rate when simultaneously analyzing multiple dependent variables, effectively reducing type-1 errors. In this study, where four different variables were examined concurrently, better results were achieved under the same p-value threshold compared to multiple univariate analyses on the same GM and WM maps. This underscores the robustness and appropriateness of a multivariate approach for understanding the intricate relationships within the dataset. Here, we utilized Bonferroni threshold for combining the p-values to build up the same power for univariate and multivariate techniques. However, there are other methods for combining p-values for detection of partial association such ad Fisher and weighted Fisher methods, which are more relevant in case of many comparisons (Yoon et al., 2021).
The nonlinear relationship between MTsat values and age, as visible in Figures Figure 6 and Figure 7, could be explained by the fact that myelination is not limited to early development, but can also occur throughout adulthood, with the pattern of myelination depending on the hierarchy of connections between different brain regions (Peters, 2002; Snaidero and Simons, 2014; Timmler and Simons, 2019). Potential nonlinear age dependency is also seen in the R2* profile in some regions such as amygdala and putamen (Figure 7) which could be due to slower accumulation of iron in these regions. Hagiwara et al. showed the nonlinear behavior of different brain tissue properties (Hagiwara et al., 2021). However, within the scope of this study, we limited our examination to linear age-related variations, as in the original analysis of these data.
The canonical vectors shown in Figure 8 indicate the contribution of each modality in the multivariate model at the peak voxels in the selected ROIs. Of note, the canonical vectors can only highlight the contribution by their size and cannot inform us about the direction of the effect.
The ROI-based partial Pearson’s correlations, presented in Table 3, are evidence of stronger correlations between the GM PD maps and age, which concur with the canonical weights regarding the contribution of PD maps within the GM.
Limitations and conclusion
It is important to note that MANOVA differs from multiple ANOVA analyses, as it does not focus on the signal-to-noise ratio of independent variable effects on each dependent variable individually. Instead, it tests for effects of interest on a combination of outcome variables. For an assessment of the former, a return to univariate analyses is necessary. Aging involves not only microstructural changes in the brain but also macrostructural and functional alterations, such as changes in GM and WM volume and cognitive behavior changes. Therefore, there remains a need to investigate aging by considering different combinations of changes together. Furthermore, longitudinal studies are required to elucidate how normal aging deviates from pathological aging within a multivariate system.
In summary, the findings of this study underscore the importance of employing advanced statistical models like the mGLM to detect subtle microstructural changes associated with aging, when using multiple maps or multimodal data. The results highlight the significance of ROI analyses in identifying specific brain regions affected by aging and their relationships with different modalities. This study provides valuable insights into the neural mechanisms underlying age-related differences in brain structure, which may have implications for developing effective interventions to slow down or prevent cognitive decline in older adults.
Data & Code Availability
The data supporting this study is available on Open Neuro: https://openneuro.org/datasets/ds005851/versions/1.0.0. Moreover, the scripts used to generate the results presented in this study are provided on GitHub (https://github.com/CyclotronResearchCentre/qMRI_MSPM_AgeFx) as part of the replication package.
Supporting information
Supplementary results[supplements/297253_file02.pdf]
Data Availability
The processed data used in this study are available online as well as the code to reproduce the results.
Funding
This work was supported by the Walloon Region in the framework of the PIT program PROTHER-WAL under grant agreement No. 7289 and ULiege Research Concerted Action (SLEEPDEM, grant 17/2109). CB and CP are Research Directors at the F.R.S.-FNRS. This research was also funded in part by the Wellcome Trust [203147/Z/16/Z]. For the purpose of Open Access, the author has applied a CC BY public copyright licence to any Author Accepted Manuscript version arising from this submission.
Disclosure statement
The authors declare that this research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary
Figure S1, Figure S2, Table S1.
Acknowledgments
The authors would like to thank E. Anderson, M. Cappelletti, R. Chowdhury, J. Diedirchsen, T.H.B. Fitzgerald, and P. Smittenaar, who took part in data acquisition as part of multiple cognitive neuroimaging studies performed at the WCHN.
Footnotes
- The overall presentation of the manuscript has been reworked, including the introduction, methods, results and discussion sections. In particular, the data as well as the scripts to reproduce the results are now publicly available. Most figures and numbers in the tables have been updated accordingly.
This pre-print is available under a Creative Commons License (Attribution 4.0 International), CC BY 4.0, as described at http://creativecommons.org/licenses/by/4.0/
References
Yoon, S., Baik, B., Park, T., Nam, D., 2021. Powerful p-value combination methods to detect incomplete association. Sci. Rep. 11, 6980. doi:10.1038/s41598-021-86465-y
Ashburner, J., 2007. A fast diffeomorphic image registration algorithm. NeuroImage 38, 95–113. doi:10.1016/j.neuroimage.2007.07.007
Ashburner, J., Friston, K.J., 2005. Unified segmentation. NeuroImage 26, 839–851. doi:10.1016/j.neuroimage.2005.02.018
Azam, S., Haque, M.E., Balakrishnan, R., Kim, I.-S., Choi, D.-K., 2021. The Ageing Brain: Molecular and Cellular Basis of Neurodegeneration. Front. Cell Dev. Biol. 9, 683459. doi:10.3389/fcell.2021.683459
Callaghan, M.F., Freund, P., Draganski, B., Anderson, E., Cappelletti, M., Chowdhury, R., Diedrichsen, J., FitzGerald, T.H.B., Smittenaar, P., Helms, G., Lutti, A., Weiskopf, N., 2014. Widespread age-related differences in the human brain microstructure revealed by quantitative magnetic resonance imaging. Neurobiol. Aging 35, 1862–1872. doi:10.1016/j.neurobiolaging.2014.02.008
Callaghan, M.F., Phillips, C., 2025. Processed MPM qMRI aging data. doi:10.18112/openneuro.ds005851.v1.0.0
Darnai, G., Nagy, S.A., Horváth, R., Ács, P., Perlaki, G., Orsi, G., Kovács, N., Altbäcker, A., Plózer, E., Tényi, D., Weintraut, R., Schwarcz, A., John, F., Varga, E., Bereczkei, T., Clemens, Z., Komoly, S., Janszky, J., 2017. Iron Concentration in Deep Gray Matter Structures is Associated with Worse Visual Memory Performance in Healthy Young Adults. J. Alzheimers Dis. 59, 675–681. doi:10.3233/JAD-170118
Draganski, B., Ashburner, J., Hutton, C., Kherif, F., Frackowiak, R.S.J., Helms, G., Weiskopf, N., 2011. Regional specificity of MRI contrast parameter changes in normal ageing revealed by voxel-based quantification (VBQ). NeuroImage 55, 1423–1434. doi:10.1016/j.neuroimage.2011.01.052
Fields, R.D., 2008. White matter in learning, cognition and psychiatric disorders. Trends Neurosci. 31, 361–370. doi:10.1016/j.tins.2008.04.001
Gyger, L., Ramponi, C., Mall, J.F., Swierkosz-Lenart, K., Stoyanov, D., Lutti, A., von Gunten, A., Kherif, F., Draganski, B., 2021. Temporal trajectory of brain tissue property changes induced by electroconvulsive therapy. NeuroImage 232, 117895. doi:10.1016/j.neuroimage.2021.117895
Hagiwara, A., Fujimoto, K., Kamagata, K., Murata, S., Irie, R., Kaga, H., Someya, Y., Andica, C., Fujita, S., Kato, S., Fukunaga, I., Wada, A., Hori, M., Tamura, Y., Kawamori, R., Watada, H., Aoki, S., 2021. Age-Related Changes in Relaxation Times, Proton Density, Myelin, and Tissue Volumes in Adult Brain Analyzed by 2-Dimensional Quantitative Synthetic Magnetic Resonance Imaging. Invest. Radiol. 56, 163–172. doi:10.1097/RLI.0000000000000720
Hotelling, H., 1951. A generalized T test and measure of multivariate dispersion. Proc. Second Berkeley Symp. Math. Stat. Probab. 23–41.
Jeremic, D., Jiménez-Díaz, L., Navarro-López, J.D., 2021. Past, present and future of therapeutic strategies against amyloid-β peptides in Alzheimer’s disease: a systematic review. Ageing Res. Rev. 72, 101496. doi:10.1016/j.arr.2021.101496
Khattar, N., Triebswetter, C., Kiely, M., Ferrucci, L., Resnick, S.M., Spencer, R.G., Bouhrara, M., 2021. Investigation of the association between cerebral iron content and myelin content in normative aging using quantitative magnetic resonance neuroimaging. NeuroImage 239, 118267. doi:10.1016/j.neuroimage.2021.118267
Moallemian, S., Salmon, E., Bahri, M.A., Beliy, N., Delhaye, E., Balteau, E., Degueldre, C., Phillips, C., Bastin, C., 2023. Multimodal imaging of microstructural cerebral alterations and loss of synaptic density in Alzheimer’s disease. Neurobiol. Aging 132, 24–35. doi:10.1016/j.neurobiolaging.2023.08.001
Peters, A., 2002. The effects of normal aging on myelin and nerve fibers: A review. J. Neurocytol. 13.
Phillips, C., Moallemian, S., Callaghan, M.F., 2025. Lifespan quantitative MR images from 138 subjects, an open and spatially preprocessed dataset. Presented at the Organization for Human Brain Mapping (OHBM).
Pillai, K.C.S., 1955. Some New Test Criteria in Multivariate Analysis. Ann. Math. Stat. 26, 117–121. doi:10.1214/aoms/1177728599
Rolls, E.T., Huang, C.-C., Lin, C.-P., Feng, J., Joliot, M., 2020. Automated anatomical labelling atlas 3. NeuroImage 206, 116189. doi:10.1016/j.neuroimage.2019.116189
Roy, S.N., 1945. The Individual Sampling Distribution of the Maximum, the Minimum and Any Intermediate of the p-Statistics on the Null-Hypothesis on JSTOR. Indian J. Stat. 1933–1960 7, 133–158.
Snaidero, N., Simons, M., 2014. Myelination at a glance. J. Cell Sci. 127, 2999–3004. doi:10.1242/jcs.151043
Steiger, T.K., Weiskopf, N., Bunzeck, N., 2016. Iron Level and Myelin Content in the Ventral Striatum Predict Memory Performance in the Aging Brain. J. Neurosci. 36, 3552–3558. doi:10.1523/JNEUROSCI.3617-15.2016
Tabachnick, B.G., Fidell, L.S., 2007. Using multivariate statistics, 5th ed. ed. Pearson/Allyn & Bacon, Boston.
Tabelow, K., Balteau, E., Ashburner, J., Callaghan, M.F., Draganski, B., Helms, G., Kherif, F., Leutritz, T., Lutti, A., Phillips, C., Reimer, E., Ruthotto, L., Seif, M., Weiskopf, N., Ziegler, G., Mohammadi, S., 2019a. hMRI – A toolbox for quantitative MRI in neuroscience and clinical research. NeuroImage 194, 191–210. doi:10.1016/j.neuroimage.2019.01.029
Tabelow, K., Balteau, E., Ashburner, J., Callaghan, M.F., Draganski, B., Helms, G., Kherif, F., Leutritz, T., Lutti, A., Phillips, C., Reimer, E., Ruthotto, L., Seif, M., Weiskopf, N., Ziegler, G., Mohammadi, S., 2019b. hMRI – A toolbox for quantitative MRI in neuroscience and clinical research. Neuroimage 194, 191–210. doi:10.1016/j.neuroimage.2019.01.029
Taubert, M., Roggenhofer, E., Melie-Garcia, L., Muller, S., Lehmann, N., Preisig, M., Vollenweider, P., Marques-Vidal, P., Lutti, A., Kherif, F., Draganski, B., 2020a. Converging patterns of aging-associated brain volume loss and tissue microstructure differences. Neurobiol. Aging 88, 108–118. doi:10.1016/j.neurobiolaging.2020.01.006
Taubert, M., Roggenhofer, E., Melie-Garcia, L., Muller, S., Lehmann, N., Preisig, M., Vollenweider, P., Marques-Vidal, P., Lutti, A., Kherif, F., Draganski, B., 2020b. Converging patterns of aging-associated brain volume loss and tissue microstructure differences. Neurobiol. Aging 88, 108–118. doi:10.1016/j.neurobiolaging.2020.01.006
Tian, Yao, Tian Yuanliangzi, Yuan, Z., Zeng, Y., Wang, S., Fan, X., Yang, D., Yang, M., 2022. Iron Metabolism in Aging and Age-Related Diseases. Int. J. Mol. Sci. 23, 3612. doi:10.3390/ijms23073612
Timmler, S., Simons, M., 2019. Grey matter myelination. Glia 67, 2063–2070. doi:10.1002/glia.23614
Wang, F., Ren, S.-Y., Chen, J.-F., Liu, K., Li, R.-X., Li, Z.-F., Hu, B., Niu, J.-Q., Xiao, L., Chan, J.R., Mei, F., 2020. Myelin degeneration and diminished myelin renewal contribute to age-related deficits in memory. Nat. Neurosci. 23, 481–486. doi:10.1038/s41593-020-0588-8
Weiskopf, N., Edwards, L.J., Helms, G., Mohammadi, S., Kirilina, E., 2021. Quantitative magnetic resonance imaging of brain anatomy and in vivo histology. Nat. Rev. Phys. 3, 570–588. doi:10.1038/s42254-021-00326-1
Wilks, S.S., 1932. Certain Generalizations in the Analysis of Variance. Biometrika 24, 471–494. doi:10.2307/2331979
Winkler, A.M., Ridgway, G.R., Webster, M.A., Smith, S.M., Nichols, T.E., 2014. Permutation inference for the general linear model. NeuroImage 92, 381–397. doi:10.1016/j.neuroimage.2014.01.060
Winkler, A.M., Webster, M.A., Brooks, J.C., Tracey, I., Smith, S.M., Nichols, T.E., 2016. Non-parametric combination and related permutation tests for neuroimaging. Hum. Brain Mapp. 37, 1486–1511. doi:10.1002/hbm.23115