Reliability and comparability of human brain structural covariance networks

Structural covariance analysis is a promising and increasingly used structural Magnetic Resonance Imaging (MRI) analysis method which characterises the co-relations of morphology between brain regions over a group of subjects. However, to our knowledge, little has been investigated in terms of the comparability of results between different data sets of healthy human subjects as well as the reliability of results over different rescan sessions and FreeSurfer versions. Our results show (significant) differences of structural covariance between two data sets of age- and sex-matched healthy human adults. This low comparability is unaltered by site correction, and is most severe when using average cortical thickness as a measure of cortical morphology. The low comparability further extends to significant differences in graph theoretic measures. Similar results of low reliability were found when comparing rescan sessions of the same subjects, and even between different FreeSurfer versions of the same scans. To understand our observations of low reliability, we estimated a measurement error covariance, and we show that it is largest in cortical thickness. On simulated data we further demonstrated that stronger errors cause stronger attenuation of correlations and proportionally decreases reliability, which agrees with associations we observed in the real data. To our knowledge this study is the first highlight the problem of reliability and comparability in structural covariance. Practically, we suggest that (1) combining data from different scanners and sites for structural covariance analysis in a naive manner should be avoided, (2) surface area and volume should be preferred as morphological measures of structural covariance over cortical thickness, and (3) some analysis of robustness of the results should be performed for every structural covariance study.


Abstract
Structural covariance analysis is a promising and increasingly used structural Magnetic Resonance Imaging (MRI) analysis method which characterises the co-relations of morphology between brain regions over a group of subjects. However, to our knowledge, little has been investigated in terms of the comparability of results between different data sets of healthy human subjects as well as the reliability of results over different rescan sessions and FreeSurfer versions.
Our results show (significant) differences of structural covariance between two data sets of age-and sex-matched healthy human adults. This low comparability is unaltered by site correction, and is most severe when using average cortical thickness as a measure of cortical morphology. The low comparability further extends to significant differences in graph theoretic measures. Similar results of low reliability were found when comparing rescan sessions of the same subjects, and even between different FreeSurfer versions of the same scans. To understand our observations of low reliability, we estimated a measurement error covariance, and we show that it is largest in cortical thickness. On simulated data we further demonstrated that stronger errors cause stronger attenuation of correlations and proportionally decreases reliability, which agrees with associations we observed in the real data.
To our knowledge this study is the first highlight the problem of reliability and comparability in structural covariance. Practically, we suggest that (1) combining data from different scanners and sites for structural covariance analysis in a naive manner should be avoided, (2) surface area and volume should be preferred as morphological measures of structural covariance over cortical thickness, and (3) some analysis of robustness of the results should be performed for 3 for volume and 1.1-7.7% for average cortical thickness [18]. However, despite a number of studies on reliability and comparability of the raw cortical morphology measures, little has been investigated on the reliability and comparability of structural covariance. The comparability aspect is of especially high relevance for studies that combine multisite data for their analysis.
To investigate comparability in this study, we compared the structural covariance for data sets of healthy human subjects with comparable demographics. For our analysis of reliability, we compared the structural covariance of scan and rescan data as well as the structural covariance of different FreeSurfer versions applied to the same dataset. We additionally investigated if, and why, different cortical morphology measures are differentially comparable and reliable.

Data selection and processing
To determine the comparability of structural covariances across scan sites we analysed structural T1-weighted MRI brain scans from two publicly available datasets: The Human Connectome Project (HCP) can be found under https://db.humanconnectome.org/ [20]. The Cambridge Centre for Ageing and Neuroscience (Cam-CAN) data is available at http://www.mrc-cbu.cam.ac.uk/datasets/camcan/ [21,22]. The MRI scanners of both data sets have 3T field strength. HCP uses a customized Siemens Skyra scanner with 0.7 mm isotropic voxel size located at Washington University, St. Louis, MO, USA. Cam-CAN uses a Siemens TIM Trio System with 1 mm isotropic voxel size located at the Medical Research Council Cognition and Brain Sciences Unit, Cambridge, UK. We selected 100 unrelated subjects from the HCP data set and used the preprocessed FreeSurfer folders provided by HCP (using a modified version of Freesurfer 5.3) [23]. The structural MRI images of Cam-CAN were pre-processed with FreeSurfer version 6.0 using the 'recon-all' script for segmentation, surface reconstruction and parcellation. We retained 644 subjects that successfully completed recon-all without errors. As regions of interest (ROI) we chose the parcellation of the Desikan-Killiany brain atlas [24] comprising 34 brain regions in each hemisphere, which is commonly chosen in covariance analysis [25,26]. We removed one subject from the HCP data as an outlier (subject ID: 414229). To reduce confounding effects of age and sex the subjects were strictly matched: we selected 43 males and 43 females in the age range of 22-34 from both data sets (age 28.26 ± 3.53 in Cam-CAN and 28.7 ± 3.49 in HCP).
To analyse reliability of structural covariances between two scan sessions we used structural T1-weighted images of healthy human subjects collected from the HCP scan and rescan data set. The data set comprises 45 subjects which were scanned at two different scanning sessions typically no more than two years apart. Again, we used the preprocessed FreeSurfer folders provided by HCP, and we also selected the Desikan-Killiany atlas.
For analysis of reliability of structural covariances between different FreeSurfer versions, we used the above-mentioned Cam-CAN data set pre-processed in  FreeSurfer version 6.0 and separately pre-processed in version 5.3. We retained 637 subjects that successfully completed recon-all without errors in both FreeSurfer versions. We again used the standard Desikan-Killiany atlas for parcellation, and selected the same 86 subjects as for our analysis of comparability. In Fig. 6 we use the above data, but instead of the Desikan Killiany atlas we select the whole left and the whole right hemisphere as ROIs. For average cortical thickness the estimates of the left and the right hemisphere are directly provided by FreeSurfer. For volume and surface area we summed the ROIs of the Desikan Killiany atlas to compute the estimates for the left and the right hemisphere. Table 1 provides an overview of all the subject numbers and demographics used.

Data analysis and visualisation 2.2.1. Statistical analysis
We standardised each measure across all subjects, if not explicitly stated otherwise, which normalizes the mean of each ROI to 0 and the standard deviation of each ROI to 1. This step has no effect on the calculation of correlations, but permits later investigation of the relative error variance (relative to the variance of 1 of the measurement). To compute differences in correlation, we Fisher transformed the correlation values first.

Covariance ellipse for visualisation purposes
To visualise some correlations/covariances of our study, we display them as an ellipse. The ellipse represents the region that contains 95% of all samples that can be drawn from the underlying covarying Gaussian distribution. We calculated the 95% confidence interval with the Chi-squared distribution. The alignment of the error ellipse is computed upon the eigenvalues of the covariance matrix of the respective ROIs. The eigenvector with the largest eigenvalue determines the direction of the the ellipsoid longer axis. Since eigenvectors are orthogonal, the shorter axis of the ellipsoid corresponds to the direction of the smaller eigenvector. Our code was inspired by Matlab source code provided by [27,28]. To compute ellipsoids for the estimated error structure and estimated underlying correlation, we use our estimated variance and covariance (see section 2.2.5).

Network measures
We computed several common network measures on the structural covariance matrix. As a first step, we thresholded (and binarised) the structural covariance matrix. We performed thresholding as a percentage of the network density. I.e. a threshold of 0.1, (10%) indicates that the top 10% of strongest edges (in absolute value) are retained. We computed different brain network measures with the brain connectivity toolbox [29]. We investigate node strength (or node degree for the binarised matrices), characteristic path and clustering coefficient. These network measures are typically used in downstream analysis of structural covariance. In the supplementary material we additionally include the global efficiency, eigenvector centrality, assortativity and k-core centrality.

Permutation test
To compare the structural covariance matrices and network measures between data sets for statistical significant differences, we apply a permutation test (Fig. 2). E.g. to compute the difference of the two matrices, we calculated the L 1 distance (sum of the absolute differences). To obtain the reference distribution, we then computed the L 1 distance of 1000 randomly mixed data splits. The p-value of the actual difference is estimated from this distribution. Such permutation tests are common in studies that are e.g. comparing a patient vs. a control group. We demonstrate here the effect of such a test on our comparisons of, for example, two healthy subject data sets. We do not use the concept of statistical significance in further downstream analyses (e.g. for the selection of specific ROIs).

Corr
Corr _ ∑ ∑ data set 1 data set 2 L distance between matrices We use the Pearson correlation coefficient Observed L 1 distance Permuted L 1 distance Figure 2: Schematic illustration of the permutation test on the example of L1 distance between the structural covariance matrices. We computed the L 1 distance (red marker) from the structural covariance matrices of the respective data sets. We repeated this process 1000 times with randomly mixed data splits to gain a distribution which serves as reference. With this distribution we can compute the p-value of the actual data split.

Estimation of error structure and underlying correlation
To analyse reliability of repeated measurements of the same subjects, we estimated an error covariance structure of the data. Such repeated measurements were for example scan vs. rescan of the same subject, or applying FreeSurfer 5.3 vs. 6.0 in the same subject and scan. Here we explain how we estimated the error and the underlying correlation for the left and right hemisphere. The estimations for all other pairs of brain regions are analogous. We denote LH 1 = X + E 1 for the first measurement of the left hemisphere and LH 2 = X + E 2 for the second measurement of the left hemisphere. Respectively we denote RH 1 = Y + D 1 for the first measurement of the right hemisphere and RH 2 = Y +D 2 for the second measurement of the right hemisphere. E 1 , E 2 , D 1 and D 2 are errors introduced at this measurement. X and Y capture the actual biological measure of interest (in reality, this will also include systematic errors of this measurement as well as random errors of other processing steps). Note that X, Y , E 1 ,E 2 ,D 1 and D 2 are vectors of the same length, where each entry in the vectors corresponds to one subject.
For our derivation we make the following independence assumptions: In other words, we assume the actual measurements (X and Y ) are independent of the errors, and the errors between the repeat measurements are independent of each other.

Calculation of the variances.
With the following quantities which can be directly calculated from the empirically measured data: we then get Calculation of the covariances.
Similarly the following covariances can be directly calculated from the empirical data: Therefore, for a given cortical morphology measure, say cortical thickness, we can infer its error variance and covariance, and hence we can visualise the error structure with an ellipse (see section 2.2.2). We show this in Fig. 6(B) and (C) in the middle column. From Var(X), Var(Y), and Cov(X,Y) we can then estimate the 'true' underlying correlation in the absence of error. This is displayed in the right column in Fig. 6(B) and (C).
We will also use the term 'attenuation', by which we refer to the difference between the measured correlations and the estimated true correlation. Usually the latter is larger than the former, hence the term attenuation of correlation.

Simulating the effect of measurement error
For our simulation in Fig. 7(A), we artificially generated data to demonstrate the effect of measurement error. We computed two sets of 1000 random data points sampled from a bivariate normal distribution with zero mean and unit variance. These data vectors represent the ground truth morphological data of two ROIs. We can generate them with any pre-defined correlation of the ROIs, which we use as the true underlying correlation.
We then added a random normal error vector (zero mean, and ν standard deviation, where ν represents the relative error/noise magnitude to the measurement) to each of the ROIs to simulate noisy measurements of the morphological data. The correlation between these simulated measurements corresponds to the measured correlations in the real data. We repeated this process 10000 times to compute the variability of the simulated measured correlations.
As a measure of reliability we computed the standard deviation over the 10000 simulated measured correlations, where a high standard deviation indicates low reliability. Further, we computed the difference between the true correlation and the mean correlation over the 10000 simulation runs as the attenuation.
We simulated two example levels of true correlation (0.7 and 0.4) each with two different example levels of relative noise (25% and 50%) in Fig. 7.

Data availability statement
Analysis code and data of all morphological measures (selected subjects of HCP, HCP scan and rescan, selected subjects of Cam-CAN preprocessed in all FreeSurfer versions) are available at Github (https://github.com/cnnp-lab/ 2019Carmon-ReliabityComparabilityStructuralCovariance).

Results
3.1. Differences in correlation structure of cortical thickness between different sites, scan-rescan, and different FreeSurfer versions.
To investigate the comparability of structural covariance we compared two data sets of healthy human subjects (HCP and Cam-CAN data restricted to the same number of males and females within a narrow age range of 22-34 year old).
In Fig. 3 we used average cortical thickness as a cortical measure, as it is widely chosen for structural covariance analysis, and standardised the data for each site in each ROI before comparison (effectively applying a site correction). In Fig. 3(A) we see that after site correction the thickness of the example ROIs are comparable between HCP and Cam-CAN. However, the correlation between the ROIs remains significantly different between the sites. We demonstrate with this example that, as expected, site correction of the univariate measures does not correct for differences in correlation. In Fig. 3(B) we visualised the comparability of datasets as differences in correlation between HCP and Cam-CAN for the whole correlation matrix. Mean differences for each ROI are shown as a heatmap on the brain.
To investigate the reliability of scan sessions, we used the HCP scan-rescan data. Similar to Fig. 3(B), we again show differences in the correlation of ROI pairs in Fig. 3(C), albeit less strong in effect than the site comparison. We show the effect of FreeSurfer version in Fig. 3(D) using the same selected subjects from the Cam-CAN data set. The results show again that there are substantial differences in correlation between the Freesurfer versions ( Fig. 3(D)).
In summary, we detected (significant) differences of structural covariance between two demographically comparable data sets of healthy human subjects. These differences are prominent despite site correction. Similarly, the structural covariance across scan sessions and FreeSurfer versions (Fig. 3(C) and (D)) also show differences, which are slightly less pronounced in effect than between sites. (A) After standardisation (µ = 0, σ = 1), two example ROIs are comparable in their cortical thickness distribution between Cam-CAN (mint green) vs. HCP (lime green). Nevertheless, the correlation of the ROIs still differ significantly between the two data sets (right panel, p = 0.003 **). Ellipsoids in the right panel visualise the correlation structure. (B) The first two columns show the full structural covariance matrix based on the Desikan-Killiany atlas. The third column displays the absolute difference between these structural covariance matrices in each matrix entry. The last column shows the average difference on the brain as a heatmap. (C) Same as (B) but for HCP scan data (blue) and rescan data (purple). (D) Also same as (B) but for Cam-CAN processed in FreeSurfer version 6.0 (red) and FreeSurfer version 5.3 (orange). Note that the colour code used for each data set will be maintained throughout the remaining figures.

Low comparability between sites in structural covariance of cortical thickness
In Fig. 3 we showed that some correlations of ROI pairs have low comparability which are substantial and statistically significant. Next, we investigate if these differences also apply to summary measures of the correlation matrix, e.g. as measured using graph theoretic network measures. Specifically, we use the data and thickness covariance matrix as in Fig. 3(B). From the thresholded (and binarised) correlation matrices we computed the L 1 distance, node strength/degree, characteristic path and clustering coefficient. Fig. 4 shows the p-values of the comparison between sites. Significant differences are seen in all network measures in both the thresholded and binarised cases for a range of thresholds. In supplementary Fig. S1 we also compare additional network measures. In our analysis we test comparability for various thresholds of the network (ranging from 2.5% to 35%).
In summary, we show that low comparability also extends to the level of the whole matrix as well as various network measures. This significant low comparability of healthy human subjects demonstrates the importance of addressing the question of comparability in structural covariance analysis.

Average cortical thickness shows largest difference in structural covariance
Average cortical thickness, surface area, and volume are three of the most popular measures to analyse structural covariance. Thus far we only investigated cortical thickness derived networks. Next, we investigated if the structural covariance of other cortical measures are differentially comparable and reliable. Fig. 5 shows the average differences of each ROI in their structural covariance using different measures of cortical morphology. The different panels display the differences between site, scan/rescan, and FreeSurfer versions, respectively. In all panels the differences are most pronounced for average cortical thickness, where more ROIs display larger differences compared to the other cortical morphology measures. We can also note, as in Fig. 3, that generally the differences between sites are larger than the differences between scan sessions and FreeSurfer versions.
In summary, the choice of the cortical measure also affects the comparability and reliability of the structural covariance. In our study average cortical thickness is the least comparable and reliable of all brain measures.

Lowest correlations, largest estimated error and strongest estimated attenuation in average cortical thickness
To investigate the cause of low comparability and low reliability we used a single correlation as an initial starting point: namely the correlation between the left and right hemispheres. On the level of a single correlation we can visualise important aspects why cortical thickness shows the prominent differences we found in Fig. 5 as opposed to area and volume.
We visualise in Fig. 6(A) the correlations of the left and right hemisphere. Each dot corresponds to a subject in one of the data sets. Each row depicts  14 one cortical measure and each column one data pair (HCP and CamCan, scan and rescan, and FreeSurfer version 5.3 and version 6.0). As to be expected from Fig. 5 we see that the difference between the data sets are largest in the correlations of average cortical thickness. We also observe that the correlations of average cortical thickness are lowest in each data set. This finding extends also to the correlations of ROIs of the Desikan Killiany atlas (see supplementary Fig. S2).
For the scan-rescan data and the FreeSurfer 5.3 and 6.0 data it is possible to estimate 1) the measurement error and 2) the underlying correlation. The rationale behind this estimation is to utilise the two measures from the same subject. From these two measures we analytically estimated the variance and covariance of the error as well as the underlying correlation (Fig. 6(B) and (C)) given several assumptions of independence -see section 2.2.5. Note that the estimated error structure is largest for average cortical thickness and smallest for surface area. Fig. 6(B) and (C) also shows that all measured correlations were attenuated by the error structure, such that the estimated underlying correlation is always larger than the empirically measured correlation. Note that the error structure estimated here is relative to unit variance of the measurement (see Methods section 2.2.1).
The findings of Fig. 6(B) and (C) also extend to most estimations of underlying correlations of ROI pairs of the Desikan Killiany atlas (see supplementary Fig. S3).
In summary, we observed that average cortical thickness is least correlated, the error of average cortical thickness is largest, and the estimated underlying correlations of average cortical thickness are most attenuated. (B and C) For subjects with repeat measurements, i.e. in the question of reliability, we could estimate the underlying relative error structure, and 'true correlation'. In both panels the left column shows the measured correlation, the middle column the estimated relative error and the right column the estimated underlying 'actual' correlation.

A larger error leads to greater levels of attenuation and unreliability
We demonstrated so far that thickness is less reliable than volume and area, and we estimated a larger error structure and a stronger attenuation of correlations for thickness. Here, we relate these observations to each other using the full Desikan-Killiany ROI correlation matrices.
In Fig. 7(A) we investigated the effect of error on reliability and attenuation for simulated data. Using simulated data with predefined true correlation levels (0.7 and 0.4), we added relative error of two different magnitudes (25% and 50% error). We repeated the simulation 10,000 times to compute the variation of the measured correlations. A smaller variation between the different simulations suggests greater reliability. In supplementary Fig. S8 we show that, in our simulations, the true correlations and the strength of the added error are wellcaptured by our estimations.
In the simulations, we note two key observations. First, stronger errors decrease reliability and proportionally increase the attenuation, regardless of the strength of true correlation. Second, for the same error magnitude, high correlations are more attenuated than low correlations. The correlation strength minimally affects the effect of error on reliability. In empirical data, and if our estimates of error are correct, we would thus expect to see error variance and covariance to be associated with the attenuation and the level of unreliability. In turn, we expect an association between unreliability and attenuation to arise. This is indeed observed in Fig. 7(B)-(D), where we used all entries of the Desikan-Killiany atlas matrix to demonstrate these associations. Fig. 7(B) shows a very high association between the estimated error covariance of the HCP scan data and the level of unreliability (difference between the measured correlations of HCP scan and rescan). In absolute terms, correlations with a stronger absolute error covariance show more difference between scan/rescan. This agrees with our simulations, where a larger error magnitude leads to a lower level of reliability. Note that unlike the estimated error covariance, the estimated error variance is the same between different measurements for standardised data (since 1 = V ar[X] = V ar[Y ], α, β and γ of section 2.2.5 are the same for X and Y ) and therefore does not contribute to the low reliability.
In Fig. 7(C), left column, we can see the negative associations between the estimated attenuation and the estimated error covariance. More negative covariances correspond to a larger decrease in the measured correlation compared to the estimated actual correlation, and vice versa. For the error variance, a larger error variance corresponds to a larger decrease in the measured correlation compared to the estimated actual correlation. In summary, more attenuated correlations are associated with stronger absolute error (co)variance, which also corresponds to our simulation.
Most correlations are positively attenuated (corrected correlation − measured correlation > 0). In Fig. 7(D) we excluded the small proportion of correlations that show a negative attenuation. The plot associates more attenuated correlations with a stronger absolute unreliability (higher absolute differences between the measured correlations of scan/rescan). In other words, more attenuated correlations tend to be less reliable. This agrees with our observation in the simulated data that reliability is proportional to attenuation.
In summary, we showed on simulated data that a stronger error causes stronger attenuation and decreases reliability proportionally. The empirical data supported this finding. We can thus link the two observations of attenuation and reliability: more strongly attenuated correlations tend to be less reliable.

Discussion
To our knowledge, this is the first study to investigate the comparability and reliability of structural covariance analysis. We analysed the comparability of data sets across sites, as well as the reliability of scan sessions and FreeSurfer versions. We showed that site differences in structural covariance are naturally not accounted for with site correction of the univariate brain morphology measures. These (significant) differences extend to common down-stream analyses of different network measures. They also persisted in scan vs. rescan data, and even in the same data processed by different FreeSurfer versions. Interestingly, we also found that the severity of the differences vary between different measures of cortical morphology. In our analysis, the structural covariance of average cortical thickness is least reliable and comparable. Structural covariance of surface area and volume are more reliable and comparable. By estimating an error covariance structure, and the underlying 'true' correlation we showed that the relative error is largest in thickness, where the attenuation of the correlation is also the strongest. Finally, by using simulated data in combination with empirical observations, we argue that it is large measurement errors that drives attenuation of correlations and low reliability and comparability, particularly for cortical thickness.
Commonly the structural covariance is computed with Pearson's correlation coefficient, to which also our main results are restricted to. An alternative correlation coefficient is the Spearman's rank correlation, known to be less prone to outliers. However, we do not find major changes to our results using Spearman's correlation (see supplementary Fig. S11). An important question to address is if there are other correlation types which are less prone to errors in covariance structure. For example, Geerligs et al. show in their study that multivariate distance correlation is more reliable than Pearson's correlation under certain circumstances [30]. Structural covariance is also measured in other ways e.g. partial least squares regression [8]. Future studies should investigate the comparability and reliability of different measures of structural covariance.
In this context of covariance estimation, it is also important to acknowledge the importance of sample size. The estimation of covariance is inherently noisy, especially with low sample sizes. Although some existing studies use relatively large sample sizes, many of those are pooled data sets across different scanners and sites, which may not be comparable in terms of their structural covariance. Large and homogeneous datasets from a single site are still rare. This essentially means that estimators of covariance that perform well with low sample sizes are needed. Future studies could investigate the effect of e.g. shrinkage estimators [31] on the reliability and comparability of structural covariance.
One way of improving reliability may be achieved by estimating the true underlying correlation from repeated measurements. Our approach could be used for such a purpose, and in artificial data we could show that the corrected correlations are close to the true underlying correlations (see supplementary Fig.  S9). However, we had to make several assumptions to arrive at the estimated corrected correlation. Future work should investigate the validity of these assumptions further in specific contexts such as scan and rescan. Especially where there are several repeated measurements the true correlation can be estimated with other frameworks (e.g. Bayesian [32]). Importantly, if the error covariance structure is derived once, then it can be applied to correct structural covariance matrices, even across scanners and sites. Theoretically, all that is required is that a group of subjects are scanned on both scanners/sites to infer the error covariance structure. Once inferred, the correction can be applied to any future subject groups to make them comparable across sites.
Additionally, reliability could also be improved by the choice of ROIs. Previous research showed that different ROIs are differently reliable in their univariate cortical morphology measures [19]. This may also extend to their covariance structure. There are many potential criteria in restricting ROIs which could improve reliability and comparability. One criterion is to select ROIs by their size. Indeed, we could find a small association of the ROI volume with the p-value of the comparability and reliability (see supplementary Fig. S10). ROIs could also be restricted by their distance (e.g. movement artifacts would cause covarying errors between ROIs). Future investigations of the source of low reliability will determin if ROI selection may be able to improve reliability.
Previously, studies reported distinctions between the structural covariance of different cortical measures [33,34]. Yang et al. found statistically significant differences between the structural covariance of cortical thickness, volume and surface area [34]. Although that observation was not made in the context of reliability and comparability, it agrees with our observation that there are some inherent differences between thickness and surface area covariances. It is alarming that in our study the structural covariance of average cortical thickness is the least reliable and comparable of the cortical morphology measures, as it is the most commonly used measure for structural covariance. It is known that cortical thickness has a lower level of biological variance compared to surface area and volume (also see supplementary Fig. S12). In relative terms, the same magnitude of error would thus have a stronger effect. In agreement, we indeed show a higher relative error for cortical thickness compared to surface area, or volume.
The results of our study are of particular interest for the comparison of different clinical conditions, as some studies combine sites to increase their sample size. In our investigation we could see significant differences between sites, even when using comparable healthy subjects with similar demographics. Thus, ideally for the comparison of clinical conditions we recommend performing the analysis for each site separately, and only test for agreement across sites. Advanced hierarchical modelling may help in estimating a more reliable and comparable covariance structure across sites in future work. This might be of particular interest for ongoing studies, especially if thickness has been used as cortical measure.
From a biological perspective, the argument for analysing structural covariance is that a strong covariance indicates co-regulation, or co-development [1]. It is however unclear if simply obtaining the covariance of a cortical measure across several ROIs is indeed the best way to capture and analyse such hypothesised co-development, or co-regulation patterns. For example, if there is no biological co-relation between two ROIs, then all that will be measured in experiments is the error and correlation in the errors. The study of covariance is a way of dimensionality reduction, and more robust method of dimensionality reduction inherently utilizing reliability (see e.g. Sotiras et al. [35] ) may be better suited for a reliable and comparable data-driven approach.
Alternatively, hypothesis-driven approaches can be taken, where a specific covariance between specific morphological variables and/or ROIs is predicted by theory. If the theory is correct, the data should support the theory in a comparable and reliable manner. One example of such a predicted covariance structure is the recently-described 'universal scaling law of cortical folding' [36,37,38]. According to this law, brain surface area, cortical thickness, and exposed surface area are linked by a strong covariance structure, which has been confirmed across species [38], within human populations [37], and even between different ROIs of the same brain [36]. Importantly, the scaling law has been shown to be comparable across sites and reliable within the same data set. Future work could take a combined data-driven and hypothesis-driven approach to discover true covariance structures in brain morphology that are reliable and comparable.
In summary, our analyses show that the question of comparability and reliability is crucial in the study and interpretation of structural covariance. Reliable or comparable univariate measures of ROIs do not imply reliable and comparable correlations between them. Practically, we make the following recommendations (1) combining sites for structural covariance analysis in a naive manner should be avoided, (2) surface area and volume should be preferred as morphological measures of structural covariance over cortical thickness, and (3) some analysis of robustness of the results (testing for reliability of the processing steps at least) should be performed. Future work should establish a pipeline of reliable and comparable structural covariance analysis based on robust data-driven dimensionality reduction and hypothesis-driven discovery of existing covariance structures.