Classification of multi-site MR images in the presence of heterogeneity using multi-task learning.

With the advent of Big Data Imaging Analytics applied to neuroimaging, datasets from multiple sites need to be pooled into larger samples. However, heterogeneity across different scanners, protocols and populations, renders the task of finding underlying disease signatures challenging. The current work investigates the value of multi-task learning in finding disease signatures that generalize across studies and populations. Herein, we present a multi-task learning type of formulation, in which different tasks are from different studies and populations being pooled together. We test this approach in an MRI study of the neuroanatomy of schizophrenia (SCZ) by pooling data from 3 different sites and populations: Philadelphia, Sao Paulo and Tianjin (50 controls and 50 patients from each site), which posed integration challenges due to variability in disease chronicity, treatment exposure, and data collection. Some existing methods are also tested for comparison purposes. Experiments show that classification accuracy of multi-site data outperformed that of single-site data and pooled data using multi-task feature learning, and also outperformed other comparison methods. Several anatomical regions were identified to be common discriminant features across sites. These included prefrontal, superior temporal, insular, anterior cingulate cortex, temporo-limbic and striatal regions consistently implicated in the pathophysiology of schizophrenia, as well as the cerebellum, precuneus, and fusiform, middle temporal, inferior parietal, postcentral, angular, lingual and middle occipital gyri. These results indicate that the proposed multi-task learning method is robust in finding consistent and reliable structural brain abnormalities associated with SCZ across different sites, in the presence of multiple sources of heterogeneity.

Introduction

Neuroimaging studies have widely explored the clinical value of machine learning methods for differentiating psychiatric patients from healthy controls at the individual level. In addition to providing individualized indices for diagnostic purposes, machine learning methods may ultimately help identify brain regions affected by disease in subtle ways that can only be elucidated using multi-variate analysis. While most of these neuroimaging studies to date have been performed using single-site datasets, it is essential to integrate multi-site data for two reasons. First, multi-site data provide sufficient statistical power for detecting subtle, but informative patterns of brain structure and function (Brown et al., 2011; Friedman et al., 2006; Schnack et al., 2010), which may be difficult to unravel with the relative small sample sizes usually acquired in single centers (Pearlson, 2009; Segall et al., 2009). Second, large sample sizes enhance sample generalizability by pooling large patient populations with diverse demographic features and clinical characteristics including disease onset, symptom severity, and types and duration of treatment (Brown et al., 2011; Friedman et al., 2006; Glover et al., 2012; Pearlson, 2009; Sutton et al., 2008; Van Horn and Toga, 2009). Multi-site studies are therefore becoming increasingly the norm in neuroimaging research (Casey et al., 1998; Van Horn and Toga, 2009).

Multi-site data reflect a more comprehensive abnormal pattern of disease, and therefore may provide a richer understanding of disease signatures than single-site data. However, two recent studies found that simple pooling of multi-site data did not outperform single-site disease classification. Colby (Colby et al., 2012) and Nielsen (Nielsen et al., 2013) pooled multi-site data and trained a common classifier for all data, to identify attention deficit hyperactivity disorder (ADHD) and autism, respectively. Results in these two studies showed that the pooled dataset exhibited lower accuracy than each single-site datasets. Here we seek a new approach to synergistically integrating multi-site data, by emphasizing two points. First, since datasets collected in multiple imaging centers have a common disorder of interest (e.g. schizophrenia (SCZ), in our experiments herein), the abnormal patterns in each dataset are strongly related and thus, to some extent, may share a common imaging signature. The pattern reproducibility among multi-site data has been repeatedly demonstrated in several multi-site studies on functional MRI (fMRI) (Casey et al., 1998; Costafreda et al., 2007; Gee et al., 2015; Jovicich et al., 2015), morphometric MRI (Cannon et al., 2014; Schnack et al., 2010), and diffusion-tensor imaging (DTI) (Jovicich et al., 2014; Pfefferbaum et al., 2003). Most of these studies arrived at similar conclusions, namely that with appropriate multi-site data collection, different data sites shared highly consistent feature patterns. With this, the site-shared features reveal consistent brain abnormalities in multi-site data, which can lead to a more accurate neurobiological understanding of the psychiatric disorder under investigation. On the other hand, though the integration of multiple single-site data is advantageous, the unavoidable presence of site-specific features might decrease the accuracy of a classifier that merely pools data together across studies. Heterogeneity can emanate from multiple sources including scanner differences, differences in image acquisition protocols, or ethnic and treatment differences among participating patient populations (Jovicich et al., 2015; Schnack et al., 2010; Van Horn and Toga, 2009). Given such site-related heterogeneity, multi-site datasets should not simply be merged into larger cohorts for further machine learning investigation (Pan and Yang, 2010). In an attempt to eliminate or reduce the site-specific variability, studies have suggested same scanning protocols, consistent scanner parameters and etc. (Brown et al., 2011; Calhoun and Adali, 2009; Pearlson, 2009; You et al., 2011) in data collection, as well as the utilization of smoothness equalization (Friedman et al., 2006) and independent component analysis (Kim et al., 2009; Meda et al., 2008) in data preprocessing. Despite these efforts, site-specific heterogeneities still exist due to their complex causes (Pearlson, 2009; Segall et al., 2009).

The above considerations highlight the need for a feature-learning framework in multi-site disease classification that can extract the site-shared features, while also accounting for the site-specific features; this approach generally seeks an overarching signature of disease, whereas it accommodates potential sub-cohort and other differences to be taken into consideration. In recent years, multi-task learning has been successful in learning task-shared and task-specific features simultaneously, which effectively improves generalization compared with traditional machine learning methods. For example, support vector machine (SVM) with single-task learning (Caruana, 1997) learns a distinct feature pattern and finds a maximum margin hyperplane to classify two groups, which extracts information within a single learning task. In contrast, multi-task learning extracts a subset of task-shared features to generate more accurate models on multiple tasks, with the task-specific features learned simultaneously (Marquand et al., 2014). The basic assumption of multi-task learning is that the feature weights of different tasks share similar sparse patterns (Chen et al., 2012), which can be learned by imposing sparsity regularization penalties on the task weight matrix (Kumar and Daume III, 2012). l1-norm and l2,1-norm are two commonly used sparsity regulating terms in multi-task learning, which enforces the weight matrix of different tasks to be sparse across all tasks. Particularly, an l1-norm term highlights task-specific features by encouraging the weights of irrelevant features to be very small (Wang et al., 2015), while l2,1-norm introduces group sparsity and enforces task-shared features to have larger weights (Watanabe et al., 2014; Yan et al., 2015).

The advantage of multi-task learning makes it suitable for multi-site data learning, considering the site as task, and the site-shared and site-specific features as task-shared and task-specific features. Neuroimaging studies have shown the effectiveness of performing multi-task learning in the brain decoding and disease classification (Marquand et al., 2014; Obozinski et al., 2010; Rao et al., 2013; Wang et al., 2015; Watanabe et al., 2014). Specifically, multi-site fMRI data of ADHD was demonstrated better than single-site classification by learning site-shared and site-specific features using multi-task scheme (Watanabe et al., 2014). In this work, though the multi-task learning scheme successfully extracted site-shared and site-specific features in multi-site data, the form of the objective function was rather complex and specific as it included an l2,1-norm group sparsity regularization term and a 6-D spatial structure penalty (generated by the GraphNet, fused Lasso, or the isotropic total variation). In order to make the multi-task learning scheme more simple and applicable, an objective function including l1-norm and l2,1-norm penalty terms was used in the current study (Wang et al., 2011; Wang et al., 2015).

Building upon this emerging literature, we aim to distinguish SCZ patients from healthy controls across multiple-site MRI data using multi-task learning. We hypothesized that using multi-task learning framework on multi-site classification would not only have better performance than single-site data classification, but would also identify the abnormalities shared by all sites, and also specific to each site. These site-shared brain structural alterations should be consistent with the previously reported altered regions in SCZ, such as the brain regions involving prefrontal, superior temporal, insular, temporo-limbic regions, among others.

Materials and methods

Participants and MRI acquisitions

MRI data were collected by three academic centers, including locations in the United States (University of Pennsylvania; site A) (Davatzikos et al., 2005), Brazil (University of Sao Paulo; site B) (Schaufelberger et al., 2007; Zanetti et al., 2013), and China (Tianjin; site C). From each site, a balanced dataset was obtained with 50 normal controls (NCs) and 50 SCZ patients randomly selected from a larger pool of available subjects. In total, we had 150 NCs and 150 SCZs, which didn't differ by age and gender significantly (p > 0.05; see Table 1).

Table 1

Characteristics of the participants in this study.

Variable		Sample size	Gender (male/female)	Age (years) Mean ± SD (range)
Site A	SCZ	50	28/22	35.38 ± 11.78(19–60)
	NC	50	25/25	32.50 ± 12.96(15–65)
	p value		0.69a	0.25b
Site B	SCZ	50	34/16	27.48 ± 7.89(18–50)
	NC	50	29/21	30.60 ± 8.17(18–50)
	p value		0.30a	0.055b
Site C	SCZ	50	25/25	34.10 ± 8.44(16–56)
	NC	50	22/28	32.24 ± 11.42(21–57)
	p value		0.55a	0.36b
Site A + B + C	SCZ	150	87/63	32.32 ± 10.08(16–60)
	NC	150	76/74	31.78 ± 10.99(15–65)
	p value		0.20a	0.66b

Note: SCZ: schizophrenia; NC: normal controls.

Pearson Chi-square test.

Two-sample t-test.

All SCZ patients met DSM-IV criteria. Written informed consent was obtained from all participants before MRI scanning. In site A, the 50 SCZ patients had chronic symptoms and were receiving treatment with antipsychotics (mean duration of illness 16.2 ± 12.3 years). In site B, all SCZ subjects were recruited shortly after they made their first contact with mental health services due to psychotic symptoms, and their duration of illness was 1.0 ± 1.3 years; 31 patients had been on antipsychotic treatment within 3 weeks of MRI, while the remaining 19 patients were free of antipsychotics at the time of MRI scanning. Site C contributed 5 first-episode, never-treated SCZ patients and 45 chronic SCZ patients under antipsychotic treatment (mean duration of illness 10.5 ± 7.2 years).

In site A, the imaging data were acquired using a Siemens Trio 3-T scanner (Siemens Medical Systems, Erlangen Germany), with the following protocol: slice thickness = 1 mm, TE = 3.51 ms, TR = 18.1 ms, flip angle = 9°, acquisition matrix = 240 × 180, and slice number = 160, no gaps, 1-mm isotropic voxels. In site B, the T1 images were acquired using two identical 1.5-T GE Signa scanners (GE Medical Systems, Milwaukee WI, USA) with the following protocol: T1-SPGR sequence providing 124 contiguous slices, TE = 5.2 ms, TR = 21.7 ms, flip angle = 20°, FOV = 22 cm, acquisition matrix = 256 × 192. Details of the MRI scanning in site B can be found in Zanetti et al. (2013). In site C, all scans were obtained using a GE 3-T Signa scanner (GE Medical Systems, Milwaukee WI, USA) with the following protocol: slice thickness = 1 mm, TE = 3.2 ms, TR = 8.2 ms, flip angle = 12°, acquisition matrix = 256 × 256, FOV = 25.6 cm.

Preprocessing

The imaging data were preprocessed using the following steps as previously described (Davatzikos, 1998; Davatzikos et al., 2001; Goldszal et al., 1998; Shen and Davatzikos, 2003): (1) skull stripping, (2) bias field correction, (3) brain tissue segmentation into gray matter (GM), white matter (WM), cerebrospinal fluid (CSF), (4) spatial registration to Montreal Neurological Institute (MNI) template, (5) generation of the Regional Analysis of Volumes Examined in Normalized Space (RAVENS) (Davatzikos et al., 2001; Goldszal et al., 1998) maps of GM, WM, and CSF using deformable registration (DRAMMS) (Ou et al., 2011), and (6) spatial smoothing of RAVENS maps using a 6-mm Full Width at Half Maximum (FWHM) Gaussian filter. In our study, we investigated the RAVENS value changes in the GM images, which indicated regional GM volumetric alterations in the brain. Note that the RAVENS method has been demonstrated to provide more accurate and reliable quantification of volumes of different brain regions (Goldszal et al., 1998). In order to reduce computational requirements, and since SCZ is not known to be associated with very localized volumetric reductions, the GM RAVENS maps were resampled to have a spatial resolution of 3 × 3 × 3 mm3. After applying a GM mask, 66,732 voxel-wise features for each subject were extracted. These images yielded a feature matrix of dimensionality of 100 × 66,732 for each site dataset. In order to ensure that features in each site are within the same level and contribute equally in the process of feature weight learning, the feature matrices were then standardized to zero mean and unit variance across all subjects within each site, respectively.

Multi-site classification using multi-task learning framework

In our study, the multi-task learning framework was used to learn the site-shared and site-specific features in the three data sites. These selected features were then used to classify the different data sites respectively. The workflow of multi-site learning is shown in Fig. 1.

Fig. 1

Workflow of each iteration in multi-site classification (with total 20 repeated experiments). (1) There were three sites in our study (sites A, B, and C). After preprocessing, the N samples in each data site were randomly divided into 5 folds and designated to three sets: one fold for testing, 80% and 20% of the other 4 folds for training and validation, respectively. Each sample had the feature dimension of D. (2) Feature learning. The training sets from three sites were used in the multi-task learning framework, which generated the feature weights of the three sites W. The three column vectors in the weight matrix W were then sorted according to their absolute values respectively, while the feature weights of the site-shared features obtained according to the K top ranked feature weights were strengthened. The new feature weights matrix W′ was generated. (3) Parameter tuning. Using the feature weights learned in the last step, the K features and their corresponding weights in the validation samples from three sites were used to classify the validation sets. The parameter set which contributed to the best classification accuracy was selected. (4) Testing. The best parameter set was used in testing and the classification accuracy was obtained.

We model the classification between SCZs and NCs of different sites as a multi-task learning problem with t learning tasks (one task for each site). In the ith learning task, X = [x, x, …, x] ∈ R is features of training samples, where d is the number of features and n is number of training samples, and their class labels are Y = [y, y, …, y] ∈ R, where y ∈ {+1, −3} stands for NCs and SCZs, respectively. Let w ∈ R denotes a feature weight vector for the ith learning task, and W = [w1, w2, …, w] ∈ R denotes the feature weight matrix of all tasks. The columns of W (w ∈ R) are the feature weights in each task, and the rows of W (w ∈ R, where p = 1, 2,⋯, d) hold the weight vectors of each feature across all learning tasks. Our goal is to learn the feature weight matrix W by minimizing the following cost function:where is the dot product of the feature weights and feature values in each training sample, and the first term in the objective function is the hinge loss cost of all training samples in all tasks, the second term ‖W‖1 = ∑|w| is l1-norm of W, the last term ‖W‖2, 1 = ∑‖w‖2 is a l2,1-norm penalty. The l1-norm term encourages task-specific sparsity, while l2,1-norm term encourages group sparsity among all tasks to identify task-shared features. The l2,1 + l1 regularization terms consider both the task-specific and task-shared features, which make multi-task learning more suitable in classification of the heterogeneous data. The optimal stochastic alternating direction method of multipliers (SADMM) was adopted to solve the optimization function (Azadi and Sra, 2014; Wang et al., 2015). Details in the solution of the optimal problem using SADMM, and the converge prove can be found in Wang et al. (2015).

In our study, the learning of feature weights at each data site was considered as a different task. Due to the sparse penalty constrains in the object function of multi-task learning, many elements in the feature weight matrix W were forced to zero, or approximately equal to zero. Weight values reflect the importance levels of different features in learning the classification boundaries between two groups, while weight values that are positive or negative tend to classify the sample into either positive or negative classes, respectively. We assume that features with larger weight coefficients (in terms of absolute value) contribute more in discriminating SCZ from NCs. Moreover, features with zero and near-zero weights make small contributions in the classification. Thus, in order to select the more important features for further classification, we sorted the weight vector in each data site according to their absolute values after the weights W were learned.

The consensus features in a small portion of the top-ranked features-which were set to 5% across three data sites, were considered as the site-shared features. The threshold of 5% was chosen because in terms of classification performance, it outperformed others in a range of 5%, 10%, and 20% (Tan et al., 2013). For a single data site, the site-specific and site-shared features were mixed together, as well as some irrelevant features. Using the multi-task feature learning scheme, the feature weights of the site-specific and site-shared features were obviously greater than the irrelevant features. Though the site-shared and site-specific features were both very important, the site-shared features were considered underlying the pathology of SCZ and more critical in characterizing disease alterations. Thus, to further emphasize the influence of site-shared features in the classification, we increased the weights of site-shared features using the following equation,where γ ≥ 1.

Note that the γ was not set to be a free parameter, which was also tuned using the nested-loop cross-validation in our study.

Classification with top ranked features

The K top ranked features with the new feature weights were used to classify each dataset. The classifier equation is shown as follows.where w denotes the new weight value of the pth feature in the ith data site, x denotes the pth feature value, and y denotes the sample label. Given one testing sample in the ith data site, the dot product of the feature weights, and the feature values of the K selected features were summed up to get the sample label.

Validation experiments

To evaluate the proposed method, ten comparison experiments were designed in the consideration of three aspects. (1) To show the effectiveness of gathering multi-site data, we reported the classification results of single-site, pooled, and multi-site data, respectively. Single-site classification learns features of each single-site data separately for classification; pooled classification pools the three data sites together as a larger dataset regardless of the site differences, which contributed to a simple binary classification of 150 controls and 150 patients; and multi-site classification uses multi-task learning to learn the site-specific and site-shared features simultaneously, in the three data sites, and combine the two types of features to classify each dataset separately. (2) To show the effectiveness of using multi-task feature learning framework, two widely used feature learning methods including principal component analysis (PCA) and two-sample t-test (ttest2) were also tested. (3) The proposed classifier based on the learned feature weights was compared to the SVM classifier with a linear kernel. Differences between these two classifiers exist as the SVM focuses on the selected features and calculates new feature weights for classification, while our proposed classifier takes advantage of the feature weights obtained in the feature learning step, and then generate the classification label. Table 2 listed the ten comparison methods in the current study.

Table 2

Summary of the comparison methods.

Method	Description
l2,1 + l1 (MS)	The proposed multi-task classification framework on multi-site data, which contains l2,1 + l1 regularization terms in the feature learning step, and the classifier based on the learned feature weights.
l1 (SS)	Single-task classification on each single site data: the single-task feature learning step uses the l1 regularization term, and the classifier uses the same classifier formula as the multi-task classification framework.
l1 (PO)	Data from all the datasets are pooled together as a larger dataset. The classification framework is the same as the single-task classification.
l2,1 + l1 + SVM (MS)	Using multi-site data, the feature learning step is multi-task learning, and uses SVM classifier with the linear kernel.
l1 + SVM (SS)	Use single-task feature learning framework and a SVM classifier to classify each single site data.
l1 + SVM (PO)	Use single-task feature learning framework and a SVM classifier to classify the pooled data.
PCA + SVM (SS)	Use principal component analysis (PCA) to learn features in each single site data and a SVM classifier with the linear kernel.
PCA + SVM (PO)	Use principal component analysis (PCA) to learn features on the pooled data and a SVM classifier with the linear kernel.
ttest2 + SVM (SS)	Use two-sample t-test (ttest2) to learn features in each single site data and a SVM classifier with the linear kernel.
ttest2 + SVM (PO)	Using two-sample t-test (ttest2) to learn features on the pooled data and a SVM classifier with the linear kernel.

In the first three experiments, the same proposed feature learning and classifier framework was applied in the multi-site, single-site and pooled data classification, respectively. In the single-site data classification, for each data site, we used the l1-norm feature selection (Wang et al., 2015), which included the hinge loss function and the l1-norm sparsity regularization term in the cost function. In the pooled data classification, the feature learning cost function was similar to the single-site classification, which considered the pooled data as a single cohort. Note that the feature learning framework applied in the single-site data classification and the pooled data classification was a special case of the multi-task feature learning framework, with the l2,1 sparsity regularization term disabled. The site-specific features were learned in the same way in all three classification experiments, while the site-shared features were only emphasized in multi-site learning, which could show the influence of site-shared features more clearly. The classifier for single-site data and pooled data classification was the same as Eq. (3), however, increasing weights of site-shared features in Eq. (2) was not applied as we did in multi-site learning, since there were no site-shared features in the single-site and pooled data classification.

In other seven experiments, the SVM classifier with a linear kernel was used after different feature learning methods. As shown in Table 2, the 4th to 6th experiments used the proposed multi-task feature learning technique, taking the feature learning method in single-site and pooled data as a special case of multi-site feature learning. In the 7th and 8th experiments, the PCA method was used, while in the 9th and 10th experiments, the ttest2 method was used in the single-site and pooled data classification, respectively.

The nested loop 5-fold cross-validation was used in the ten experiments to evaluate the classification performance. In specific, the samples were randomly separated into five subsets: one (20% of all samples) was used for testing, and the other four (80% of all samples) were divided further into two parts with the proportion rates of 80% (64% of all samples) and 20% (16% of all samples), respectively. The larger of the two parts was used as a training set, and the smaller part was used as a validation set, as described previously (Koutsouleris et al., 2014). The training samples were used for feature learning and classifier training, while the validation samples were used to tune all parameters with the trained classifier, and the best parameter set which contributed to the best classification accuracy of the validation samples was selected. Finally, the best parameters in the validation step were used to classify the testing samples and obtain the classification result. The whole process was repeated for 20 repeated experiments. The average accuracies of the 20 iterations for each site was considered as the final classification result.

The parameters of the multi-task feature learning framework were tuned by grid searching the regularization factors in the range of α, β ∈ {10−8, 10−6, 10−4, 10−2, 10−1, 0.3, 0.5, 1, 5, 10, 102, 103}, while in the special case of using multi-task feature learning in the single-site and pooled data, the l2,1-norm term was disabled and β was set to 0. In the proposed classifier in the multi-site classification (the first experiment), the ratio to increase the site-shared feature weights γ was tuned in the range of {1, 2, …, 100}. The number of selected features K ranged from {500, …, 20,000}, with the step-size of 100. In the experiments with the SVM classifier, the linear kernel was used and the cost value C was set to 1.

The site-shared and site-specific features were obtained in the multi-site classification, which used the proposed multi-task learning framework. After the parameter tuning, the best parameter set of the regulation terms α, β, the number of selected features K, and the site-shared feature increasing ratio γ were obtained. In each iteration, the best classification accuracy in each site corresponded to a specific W vector. We first calculated the mean W vectors of the 20 iterations that were contributed to the best accuracy in all three sites respectively, then the mean W vector of sites A, B, and C was obtained. The consensus features with weights that ranked the top 5% of the three sites were shown as the site-shared features, and the remaining features were considered as the site-specific features for each site, respectively. The feature weights of the site-shared and site-specific features were displayed using the Caret package (Van Essen, 2005), in which the clusters comprising of >50 voxels were shown.

Results

The classification results of multi-site, single-site, and pooled data were listed in Table 3. Using the proposed multi-task feature learning framework and the classifier based on the learned feature weights, the average classification accuracies of sites A, B, and C, based on the multi-site data learning, were 78.5%, 64.5%, and 84.5%, which outperformed the accuracies of 76%, 56%, and 82.5% on single-site learning, as well as the accuracies of 70%, 57.5%, and 75.3% on the pooled classification.

Table 3

The average accuracy values of multi-site, single-site and pooling classification with 20 repetitions in 10 experiments.

Method	Site A	Site B	Site C	Average of 3 sites
l2,1 + l1 (MS)	0.785	0.645	0.845	0.758
l1 (SS)	0.76	0.56	0.825	0.713
l1 (PO)	0.70	0.575	0.753	0.676
l2,1 + l1 + SVM (MS)	0.713	0.648	0.658	0.673
l1 + SVM (SS)	0.728	0.598	0.825	0.690
l1 + SVM (PO)	0.610	0.573	0.668	0.617
PCA + SVM (SS)	0.668	0.474	0.724	0.622
PCA + SVM (PO)	0.646	0.494	0.732	0.624
ttest2 + SVM (SS)	0.562	0.468	0.680	0.570
ttest2 + SVM (PO)	0.672	0.530	0.646	0.616

MS: multi-site data; SS: single-site data; PO: pooled data.

Using the multi-task feature learning framework and the SVM classifier, the average accuracies in multi-site (67.3%), single-site (69.0%) and pooled (61.7%) data classification were lower than the results using the proposed classifier (75.8%, 71.3%, and 67.6%), respectively. However, the average accuracy of multi-site classification was poorer than the single-site classification using multi-task feature learning framework and SVM classifiers. In the experiments using PCA feature learning method and SVM classifiers, the average results of single-site and pooled data classification were 62.2% and 62.4%. In the experiments using ttest2 feature learning method, the classification accuracies were the lowest comparing to all other methods, which were 57.0% and 61.6%, respectively in single-site and pooled data classification.

The best average accuracy was 75.8% with our proposed multi-task feature learning classification method. The results indicated that multi-site learning captured more precise feature patterns in revealing the differences between SCZ patients and controls, than using only one single-site data. Simply pooling all data sites together did not improve classification performance compared with using single-site data.

The site-shared and site-specific features with the largest weights were mapped onto a brain image template. As shown in the left panel of Fig. 2, multi-site learning extracted a series of feature patterns as site-shared features, which revealed gray matter abnormalities in a variety of regions, including the fusiform, middle temporal, superior temporal, inferior parietal, postcentral, angular, inferior frontal, middle occipital and lingual gyri, insular, anterior cingulate cortex, precuneus, the dorsolateral prefrontal cortex, and the cerebellum. The most discriminative site-shared regions are all listed in Table 4. Detailed information including the BA area, side, cluster size and peak MNI coordinates is also shown in Table 4.

Fig. 2

Left panel: locations of the site-shared features Right panel: locations of the site-specific features corresponding to each site. These site-shared and site-specific features were shown with cluster size of >50 voxels obtained by multi-task learning on multi-site schizophrenia classification. The colorbar represents the weight values of features. The warm and hot colors corresponded to negative and positive weight values, respectively.

Table 4

The site-shared gray matter alteration features of brain regions.

Regions	Side	BA	Cluster size (voxels)	MNI coordinates
Cerebellum crus I, fusiform gyrus, lingual gyrus	L	18,19	146	−32,−86,−17/−29,−74,−3
Middle temporal gyrus	L	38,21,22	77	−44,12,−35/−59,−14,−8/−56,−59,−13/−65,−35,−1
Superior temporal gyrus, middle temporal gyrus, insula	R	22,41,13	77	52,−41,−10/49,−17,1
Anterior cingulate, superior medial frontal gyrus	L	32	51	−12,34,1/−12,34,25
Middle occipital gyrus, angular	L/R	39,40	111	45,−55,20/33,−76,32/−51,67,39
Insula	R	13	27	33,16,12
Middle frontal gyrus, inferior triangular frontal gyrus	L/R	9,10,44,45	162	−51,26,28/−33,52,16/−42,22,13/39,35,24
Postcentral	L	2	18	−49,−31,50
Precuneus	L/R	31	25	−6,−43,36
Inferior parietal gyrus	L	40	15	−40,−49,51
Superior parietal gyrus	L	7,5	42	−28,−51,75
Superior frontal gyrus, middle frontal gyrus	R	9,6	56	27,47,32/33,41,35/35,8,42
Superior frontal gyrus, middle frontal gyrus	L	32,6	63	−16,23,47/−34,11,52
dlPFC	L/R	9,6,32	119	−16,23,47/−34,11,52/27,47,32/33,41,35/35,8,42

Note: BA = Brodmann area, R = right, L = left.

The site-specific features were displayed in the right panel of Fig. 2. The site-specific features of site A included regions of the fusiform, lingual, inferior temporal, middle temporal, superior temporal, middle frontal, superior frontal, superior medial frontal, postcentral, superior parietal, precentral gyri, precuneus, putamen, supplemental motor area, and the cerebellum. The site-specific features of site B included regions of the fusiform, anterior cingulate cortex, posterior cingulate cortex, hippocampus, parahippocampus, angular, inferior frontal, middle frontal, postcentral, precentral, inferior parietal, middle occipital gyri, precuneus, and the cerebellum. The site-specific features of site C included regions of the fusiform, lingual, anterior cingulate cortex, angular, middle temporal, inferior temporal, superior temporal, middle frontal, postcentral, supramarginal, superior medial frontal gyri, precuneus, insular, and the cerebellum. As shown, the site-specific patterns in the three sites involved several regions that were located adjacent to the same site-shared regions, including the fusiform, lingual, postcentral, precentral, prefrontal gyri, and the cerebellum. Detailed information of the site-specific features are shown in Table 5, Table 6, Table 7.

Table 5

The site-specific gray matter alteration features of brain regions in site A.

Regions	Side	BA	Cluster size (voxels)	MNI coordinates
Fusiform gyrus, cerebellum VI, lingual gyrus	L	37,19,18	135	−20,−59,−15
Inferior temporal gyrus, middle temporal gyrus	L	37,20,19,21	50	−56,−56,−10
Superior temporal gyrus, middle temporal gyrus	L/R	38,21,41	148	−47,9,−32/55,−53,−4
Putamen	R		69	30,8,15
Middle frontal gyrus/superior frontal gyrus	R	9,8,10	114	29,38,38
Postcentral gyrus, superior parietal gyrus	R	3,2,5,40	236	−43,−37,65/26,−46,49
Superior medial frontal gyrus, supplemental motor area, superior frontal gyrus	L	8,6	203	−19,35,58
Superior parietal gyrus, postcentral gyrus, precuneus gyrus	L	7,5	92	−28,−54,75
Precentral gyrus	L	6	54	−34,−12,77

Note: BA = Brodmann area, R = right, L = left.

Table 6

The site-specific gray matter alteration features of brain regions in site B.

Regions	Side	BA	Cluster size (voxels)	MNI coordinates
Cerebellum IV, V, fusiform gyrus	R	19,36	51	22,−35,−20
Hippocampus, parahippocampus, fusiform gyrus	R	36	75	34,−35,−11
Anterior cingulate cortex	L	32,9,24,10	294	−3,38,23
Inferior parietal gyrus, middle occipital gyrus, angular	R	40,39,7	105	−46,−49,47
Posterior cingulum	L		65	−9,−43,36
Inferior frontal gyrus	L	13,44,47	38	−39,22,13
Middle frontal gyrus	L	8,9	53	−34,35,37
Precuneus	L	7	120	−13,−67,37
Postcentral gyrus, precentral gyrus	L/R	6,4,3,40	180	−52,−31,53

Note: BA = Brodmann area, R = right, L = left, dlPFC = dorsolateral prefrontal cortex.

Table 7

The site-specific gray matter alteration features of brain regions in site C.

Regions	Side	BA	Cluster size (voxels)	MNI coordinates
Cerebellum crus I, crus II, fusiform gyrus, lingual gyrus	L	18,19	209	−28,−84,−38
Middle temporal gyrus, inferior temporal gyrus	L	21,20,22	274	−59,−12,−26
Superior temporal gyrus, middle temporal gyrus, insula	R	22,13,41,21	112	46,−14,−2
Anterior cingulate cortex, superior medial frontal gyrus	L	10,32,9	122	0,46,14
Angular	R	39,22,13	62	45,−55,23
Middle frontal gyrus	L	9,10	81	−28,56,37
Supramarginal gyrus	L	40	80	−54,−37,29
Postcentral gyrus	R	2,40,3	62	44,−31,55
Precuneus	L/R	7	72	2,−66,67

Note: BA = Brodmann area, R = right, L = left.

Discussion

We present a method utilizing multi-task l2,1 + l1-norm learning to simultaneously learn the site-shared and site-specific features, for multi-site classification problems. The method was tested to differentiate SCZ patients from healthy controls in three data sites. Compared to single-site and pooled data classification, the multi-site classification performance was significantly higher, which demonstrated that multi-site data can enhance accuracy in identifying brain alterations in SCZ if multi-task learning is applied. The multi-task learning framework thus provides a new way to take advantage of multi-site data to better classify and understand brain diseases.

Reliable identification of multi-site schizophrenia data using multi-task learning

Single-site data contains information about pathophysiology, but may simultaneously include information related to different scanners or imaging protocols, as well as demographic differences that may influence the pattern investigations across individual subject samples. Using single-site data is limited by smaller sample sizes and smaller diversity in imaging features, both of which often lead to poor generalization performance of classifiers on new patient populations. Data pooled from multiple sites alleviate both of these problems at the expense of introducing confounding heterogeneity in the imaging data, which renders the pooled analysis challenging. In our study, using the proposed multi-task feature learning scheme, the classification accuracy of using pooled data was no better than using each single site and, in the case of site C data, pooled learning performed worse than single-site classification. This result might demonstrate that the site-specific features in each data site may not always benefit for the pooled classification. The advantage of the multi-task learning frame in our study is that it can simultaneously learn site-shared and site-specific features, which is suitable when combining multi-site data for disease classification. The l2,1 + l1-norm penalty terms focuses on the group sparsity, and sparsity to maximize weights of site-shared or site-specific features, while setting all other irrelevant features to be very small. In addition, the two penalty terms are very common and simple in sparsity regularization and can easily be generalized. The multi-task learning framework has attracted increasingly attention in other fields, such as in data mining, bioinformatics and computer vision (Heisele et al., 2001; Ji and Ye, 2009; Xue et al., 2007).

In the current study, the performance was poor in using the PCA and ttest2 feature learning methods in the single-site and pooled classification. PCA is a linear method for dimensionality reduction, which may not consider the nonlinear correlations of different features, while ttest2 is a univariate feature selection method, neglecting the relationship of multiple variables. Thus, the accuracies of using PCA and ttest2 were no better than using the proposed multi-task feature learning framework. Results also showed that accuracies using the proposed classifier were much better than using SVM. In the multi-task feature learning step, the feature weights were calculated based on all features and contained valuable information. Those features exhibited the largest weights were selected and used in the classification step. SVM classifier only used the sequencing information of the feature weights, which had lost some valuable information of these selected features. Therefore, better results using our proposed classifier showed that the feature weights learned by the multi-task feature learning method were accurate.

Site-shared and site-specific regions of gray matter alterations in schizophrenia

A pattern of gray matter intensity alterations implicating several brain regions in the SCZ sample was revealed by the site-shared feature analysis of multi-site MRI data in our study. This is consistent with a vast number of large-scale structural MRI investigations of SCZ, including both meta-analyses of MRI results (Chan et al., 2011; Ellison-Wright et al., 2008; Honea et al., 2005; Palaniyappan et al., 2012; Steen et al., 2006; Vita et al., 2012; Wright et al., 2000), and mega-group comparisons of mean indices combining data from multiple sites (Gupta et al., 2015; Segall et al., 2009; Torres et al., 2016). These previous studies have repeatedly highlighted the prefrontal cortex, lateral temporal cortex, insular, striatum, anterior cingulate cortex, hippocampus and other temporo-limbic structures as key components of the network of brain regions showing volumetric deficits relative to healthy controls from early stages of SCZ onwards. Most of those brain regions were implicated in our study, indicating a high degree of consistency between the spatial distribution of findings in our site-shared feature analysis and the previous MRI literature on SCZ. The brain regions involved in this network are critical to emotion processing and several cognitive domains known to be affected in schizophrenia, including language processing, working memory and other executive functioning operations, episodic memory, attention, self-monitoring, error detection, and attribution of salience to emotionally relevant stimuli, among others (Backasch et al., 2014; Belin et al., 2000; Cabeza and Nyberg, 2000; Chao et al., 1999; Jeong et al., 2009; Kim et al., 2003; Noppeney and Price, 2002; Ojemann et al., 2002; Palaniyappan and Liddle, 2012a; Palaniyappan and Liddle, 2012b; Potkin et al., 2009; Price, 2010; Tranel et al., 1997; Wright et al., 2003). Thus, the findings presented herein provide evidence that our method robustly identifies consistent and reliable structural brain abnormalities associated with SCZ across different sites.

Additional brain regions often implicated in the pathophysiology of SCZ were also highlighted in our site-shared analysis, such as the cerebellum (Andreasen and Pierson, 2008). Though the cerebellum is traditionally believed to participate in motor function, neuroimaging studies, in both nonhuman primates and humans, have demonstrated that it is also involved in brain cognitive and affective processing (Bernard and Mittal, 2015; Dum and Strick, 2003; Hu et al., 2008; Kelly and Strick, 2003). Several structural and functional neuroimaging studies have shown cerebellar abnormalities in SCZ subjects relative to controls (Bernard and Mittal, 2015; Gupta et al., 2015; Ha et al., 2004; Kim et al., 2014; Mothersill et al., 2016; Rusch et al., 2007; Shen et al., 2010; Shenton et al., 2001; Wagshal et al., 2015). Abnormalities in the cerebellum, through its connections to the prefrontal cortex via thalamus, have been proposed as critical to the state of “cognitive dysmetria” put forward by some authors to explain the behavior manifestations and cognitive deficits of schizophrenia (Andreasen et al., 1999; Rusch et al., 2007; Shen et al., 2010). However cerebellar findings should be interpreted cautiously, since segmentation between GM and WM in that region is difficult and prone to errors, due to poor contrast at this imaging resolution. Finally, other brain sites showing volume abnormalities in our study, including the precuneus and lateral parietal cortical regions, fusiform, lingual and middle occipital gyri, have been variably implicated in previous neuroimaging studies of SCZ (Borgwardt et al., 2010; Brunet-Gouet and Decety, 2006; Chou et al., 2014; Davatzikos et al., 2005; Di Rosa et al., 2009; Dong et al., 2006; Faget-Agius et al., 2012; Gaser et al., 1999; Herrmann et al., 2006; Johnston et al., 2005; Lee et al., 2002; Mashal et al., 2014; Nierenberg et al., 2005; Niznikiewicz et al., 2000; Onitsuka et al., 2005; Onitsuka et al., 2003; Quarantelli et al., 2014; Schultz et al., 2010; Seiferth et al., 2008; Shenton et al., 2001; Singh et al., 2014; Takahashi et al., 2011; Tanskanen et al., 2010; Walther et al., 2009; Whalley et al., 2006; Yoon et al., 2006). These brain regions are involved in a variety of cognitive operations of relevance to schizophrenia, including language, visual attention and emotional processing (Cohen et al., 1996; Dong et al., 2006; Herrmann et al., 2006; Johnston et al., 2005; Kaas, 1995; Mashal et al., 2014; Nierenberg et al., 2005; Niznikiewicz et al., 2000; Seiferth et al., 2008; Walther et al., 2009; Whalley et al., 2006; Yoon et al., 2006).

Findings of structural brain abnormalities in SCZ may be influenced by clinical variables not specifically related to the pathophysiology of the disorder, such as the use of antipsychotic medications. Several recent studies have shown that the continued use of antipsychotics may be significantly associated with some of the neuroanatomical changes detected in MRI studies of SCZ, even after taking into account the influence of other potential moderators such as illness severity and substance abuse (Fusar-Poli et al., 2013; Ho et al., 2011; Torres et al., 2016). In the present study, participants with SCZ from site B were all in their first psychotic episode and had a relatively shorter duration of illness in comparison to the samples contributed by sites A and C (both of which included a substantial proportion of chronic, medicated SCZ patients). Thus the overall SCZ sample from site B was less exposed to antipsychotic treatment in comparison to the two other sites. The selection of first-episode SCZ patients in site B may help to explain the more modest diagnostic accuracy figures obtained in the single-site analysis for that center in comparison to the two other sites, given the fact that medicated SCZ patients in chronic disease stages display more widespread brain volumetric abnormalities relative to healthy controls than first-episode SCZ patients (Meisenzahl et al., 2008; Torres et al., 2016). Since the three imaging sites contributed each an equal number of SCZ patients and controls to the present investigation, disease chronicity and antipsychotic exposure should be taken as clinical variables that increased inter-site variability in our study, instead of representing unifying features present in all sites. Therefore, the results of the site-shared analysis presented herein can be more safely assumed to reflect the commonality of SCZ neuropathology across the three sites, rather than the potential influence of confounding clinical variables ubiquitously present in the three centers. Disentangling the specific neurobiological features underlying psychiatric disorders from clinically relevant confounders is challenging in neuroimaging research, given the highly prevalent influence of factors such as use of nicotine and other drugs of abuse, cardiovascular risk factors, and exposure to pharmacological treatment, all of which can affect brain structure (Weinberger and Radulescu, 2016).

The site-specific regions of the three sites revealed in our study consisted of two parts. The larger part of the site-specific features was located in the same regions of the site-shared features, though the coordinates were slightly different. These regions involved most of the site-shared features, including the fusiform, lingual, postcentral, precuneus, prefrontal gyri and the cerebellum. The findings further implicated that these site-shared and also the site-specific regions revealed the abnormal feature patterns of the SCZ accurately, and multi-task feature learning framework provided a promising method in detecting these features. However, though these site-specific features were located in the same regions with the site-shared features, the feature distributed in adjacent areas. This was reasonable due to the data heterogeneity resulted from each site. Another part of the site-specific regions included some site-specific regions found exclusively in each site, such as putamen and supplemental motor area in site A, hippocampus, parahippocampus, posterior cingulate cortex in site B, and supramarginal gyrus in site C.

Studies have shown that supplemental motor area dysfunction was related to motor disturbances in SCZ patients (Schroder et al., 1995). The putamen was proposed to be related to the treatment response (Buchsbaum et al., 2003; Li et al., 2012; Mitelman et al., 2009), and also reflected the deficits of auditory verbal hallucinations of SCZ patients (Cui et al., 2016). The gyrus of hippocampus and parahippocampus were reduced (Razi et al., 1999), and hippocampus was demonstrated to be associated with the positive symptoms of SCZ patients (Duan et al., 2015). Meanwhile, the posterior cingulate cortex has recently been implicated in the pathophysiology of schizophrenia (Newell et al., 2006). These site-specific regions may provide potential evidences on interpreting the underlying multisite data heterogeneity introduced by different scanners, subjects, treatment exposure, or disease symptoms or etc.

Limitations and future directions

Although the classification result of multi-site MRI data using multi-task learning was encouraging, there are limitations in our work. First, due to the limited sample sizes, small site numbers, and the specific psychiatric disease type and data modality in our study, the effectiveness of using multi-task learning scheme on multi-site classification needs to be further validated in larger datasets. Moreover, an important challenge in multi-task learning is that when new data is analyzed, the multi-task learner would need to be recalculated. This would obviously limit clinical utility. One solution is to apply our approach to obtain a set of t discriminants, each associated with one of our training sites. It is possible that as new data and sites are added to our analyses, one would not need to solve the multi-task learning problem again, but rather could apply all these t vectors and pool the classification results via weighted consensus approaches.

Conclusion

This work proposes a multi-task learning framework for high-dimensional classification, which simultaneously learns the site-specific and site-shared features in multi-site MRI data of brain morphology. This method was tested on multi-site anatomical images of SCZ patients. The classification results of using multi-task learning outperformed single-task learning and the pooled data. The most significant site-shared features revealed brain structural alterations in SCZ consistent with those previously reported in the literature. In summary, multi-task learning provides a promising technique to learn the site-shared features in multi-site data study without neglecting site-specific features, which gives a new angle on integrating multiple single-site data together for big data studies.