The prediction of behavioral information or cognitive states from brain activation images such as those obtained with fMRI can be used to assess the specificity of several brain regions for certain cognitive or perceptual functions. This kind of analysis is implemented by learning a classifier or regression function that fits a giventarget variable given fMRI activations. The accuracy of this prediction depends on whether it uses the relevant variables i.e. the correct brain regions. Recovering the truly predictive pattern has proven to be challenging from a statistical point of view: the high dimensionality of the data together with the limited number of images makes the problem of brain pattern recovery an ill-posed problem.
So far, the approaches proposed to address this issue have relied on linear models, with univariate, i.e.
voxel-based, Anova (analysis of variance) for hypothesis testing, or, for predictive modeling, with the choice of a regularizer using a priori domain-specific knowledge, such as the-norm to promote sparsity [yamashita2008, Ryali_Supekar_Abrams_Menon_2010], total variation to promote spatial smoothness [Michel2011]
. Various data fit terms have been used, Logistic Regression (LR)[Ryali_Supekar_Abrams_Menon_2010], Linear SVM [mouraomiranda2005], Lasso [yamashita2008]
. While Linear SVM and LR cannot address the recovery problem for multiclass problems, linear regression models assume a linear relationship between the quantity to predict and the amplitude of the fMRI signals. If the linear relationship does not hold in practice, then the estimation of the predictive patterns may suffer from a loss of statistical power. This can be particularly relevant with Blood Oxygen-Level Dependent (BOLD) signals observed in fMRI, where a saturation effect is expected as the level of signal increases.
When targets to predict consist of ordered values, as in a parametric design, such as clinical scores, pain levels or the complexity of a cognitive task, the response to these different conditions can reflect the non-linearities in the data. In such situation, we propose to use a data fit term, known as loss function, not relying on an assumption of linearity but only of increasing response. We show on simulations that this new formulation opens the door to capturing the non-linearity and leads to better recovery of the predictive brain patterns. On an fMRI dataset we show that the new formulation leads to models with better recovery properties.
We write vectors in bold,, matrices with capital bold letters, . The dot product between two vectors is denoted . We denote by the norm of a vector.
Ii Learning a linear model from fMRI data
Following standard statistical learning notations we denote by , , the data and the target variables. In this paper, we aim at learning a weight vector such that the prediction of can be non-linearly related to the value of . The vector corresponds here to a brain map that can be represented in brain space as a volume for visualization of the predictive pattern of voxels. It is useful to rewrite these quantities in matrix form; more precisely, we denote by the design matrix assembled from fMRI volumes and by the corresponding targets. In other words, each row of is a -dimensional sample, i.e., an activation map of voxels related to one stimulus presentation.
A standard approach to perform the estimation of leads to the following minimization problem
where is the regularization term and is the loss function. The parameter balances the loss function and the penalty .
If the explained variable is a linear combination of the images, , we can estimate using the mean squared error loss function . However, with fMRI the linearity assumption may not be valid. Instead, we model our explained variable as , where is a non-decreasing function.
We introduce the use of pairwise loss functions. These loss functions only assume the target values to be a non-decreasing function of the data. They have been widely used in ranking, a type of supervised machine learning problem whose goal is to automatically construct an order from the training data. A pairwise loss function operates on pairs of images: given a pair of imagesand , we build a model that predicts the sign of .
Let denote the index set of all considered pairs. Note that in some settings it might be important to restrict ourselves to a selected subgroup of all pairs, e.g. to the pairs of images of a single subject or to the pairs of images corresponding to a single session. For this purpose we define to be a weight associated to each pair. We will now present the pairwise loss functions used in this article:
Pairwise hinge loss [Herbrich2000]
. This is a natural extension of the loss function used by Support Vector Machines and has been successfully used in information retrieval[Cao:2006:ARS:1148170.1148205].
Pairwise logistic loss [Dekel_Manning_Singer_2003]. This is the pairwise extension of the logistic regression loss function.
When noise is present in the model, the order of two samples might get inverted, a phenomenon known as label switching. Because this only affects labels that lie close, it is natural to penalize more the misclassification of distant labels. By setting the sample weights to to a value that increases as labels become more separated, we become more robust to label switching. In the case of hinge loss functions, several strategies for choosing the appropriate weights are discussed in [Cao:2006:ARS:1148170.1148205].
On the implementation side, both pairwise hinge loss and pairwise logistic loss can be implemented on top of existing SVM and Logistic Regression solvers, respectively, by taking the difference of pairs as input values. In practice, we used the liblinear [Fan_Chang_Hsieh_Wang_Lin_2008] library via the scikit-learn [scikit-learn] library.
The simulated data contains volumes of size and
, each one consisting of Gaussian white noise smoothed by a Gaussian kernel with standard deviation of 2 voxels. This mimics the spatial correlation structure observed in real fMRI data. The simulated vector of coefficientshas a support restricted to four cubic Regions of Interest (ROIs) of size (). The values of restricted to these ROIs are .
The target variable is simulated as a linear model:
where the noise
follows a uniform distribution.is chosen such that the signal-to-noise ratio verifies . We then define another target variable
to be a sigmoid function of, that is,
For each of the loss functions introduced earlier and mean squared error, we compute the correlation coefficient . This gives us the goodness of fit for the estimated . A method with correlation coefficient of 1 is able to recover perfectly the ground truth. Since we are interested in the estimation of , we restricted ourselves to linear models, leaving out models such as kernel ridge or support vector regression with Gaussian kernel.
For all models we added regularization in the form of and we cross-validated separately for each loss. In this setting, the model with mean squared error loss is equivalent to ridge regression. In order to focus on the effect of non-linearity, we first considered a noiseless simulation using trivial weights .
Once we learned a vector for each method we compute the correlation coefficient with the ground truth for different sizes of the training data. The experiment is repeated 10 times with different initialization of the pseudorandom number generator. We compute errorbars and show the results in figure 1. As the number of samples increases, the linear model stalls and pairwise loss functions outperform MSE on both and
dimension. As expected, the higher dimensionality of the second simulation makes the correlation coefficient decrease. However, unlike MSE, ranking tends to a perfect recovery as the number of samples increases. Both pairwise loss functions perform equivalently and have a significatively higher correlation coefficient than MSE. In the rest of the paper we will use pairwise logistic as loss function. As a result, pairwise loss functions should be preferred over MSE in situations where underlying model is non-linear. Notice that we fixed the non-linearity to be a sigmoid function, but the pairwise loss functions only assume that this function is non-decreasing. Unlike linear regression models, pairwise loss functions are indeed able to learn the structure on the non-linear transform.
We now consider the model with noise as in (4) and use non-trivial weights . To account for label switching, we set to zero for pairs with too similar labels:
In the case of discrete values, this would be equivalent to zeroing weights for which the labels are adjacent. We now compute the correlation coefficient for weighted and unweighted pairwise logistic models and linear ridge regression model. The result can be seen in figure 2 for dimension . The unweighted logistic model breaks down in presence of noise and performs worse than linear ridge regression. On the other hand, appropriately setting the weights has a major effect on robustness, where this model outperforms MSE in a noisy setting. Note also that weighted pairwise logistic has smaller variance than MSE.
Iv Results on fMRI data
This dataset, described in [Cauvet2012], consists of 34 healthy volunteers scanned while listening to 16 words sentences with five different levels of complexity. These were 1 word constituent phrases (the simplest), 2 words, 4 words, 8 words and 16 words respectively, corresponding to 5 levels of complexity which was used as class label in our experiments. To clarify, a sentence with 16 words using 2 words constituents is formed by a series of 8 pairs of words. Words in each pair have a common meaning but there is no meaning between each pair. A sentence has therefore the highest complexity when all the 16 words form a meaningful sentence.
The dataset consists of 8 manually labeled ROIs, some informative and some not. For each ROI separately, we split the data into 60% training samples, 20% for parameter selection and 20% for validation. We trained a pairwise logistic model and set the regularization by cross validation. We choose to be zero if classes are adjacent, i.e. if and if and do not belong to the same subject, in order to consider exclusively non-adjacent pairs of images from the same subject. In all other cases, was set to one. We computed the generalization score on the validation set as the mean number of inversions with respect to the order in labels, i.e. the sign flips for all pairs of images in the validation set.
We kept the four ROIs with highest scores (see figure 3). These are: Anterior Superior Temporal Sulcus (aSTS), Temporal Pole (TP), Inferior Frontal Gyrus Orbitalis (IFGorb) and Inferior Frontal Gyrus triangularis (IFG tri).
In order to inspect the properties of the estimated functions for each ROI, we estimated using a pairwise logistic model. We then projected our data along this vector and regularized the result using non parametric locally weighted scatterplot smoothing (LOWESS). Results in figure 4 show that the linearity varies in shape across ROIs which suggests that different regions exhibit different sensitivities to the complexity parameter under investigation. We see however a trend in the figures towards non-linear and non-decreasing functions with some saturation effect of the BOLD signal as in the temporal pole (TP).
In the case of the Temporal Pole (TP), which is the ROI revealing the highest saturation effect, an F-test on the data reveals that the quadratic polynomial model fits better the data than a linear model (p-value ). As shown in the simulations, in this particular case pairwise loss functions are likely to improve the recovery of active brain regions.
In this paper, we investigated the use of pairwise loss functions to improve the problem of brain pattern recovery with supervised learning applied to fMRI data. Through simulations, we showed the benefit of such loss functions when the target to predict is non-linearly related to the voxel amplitudes. Experimental results on fMRI data confirmed the presence of such non-linear effects in the data which suggest that the pairwise approach should improve the identification of predictive brain patterns on experimental data.
This work shows that improvements in recovery of brain activation patterns should not only rely on the choice of a particular regularizer, but also on an appropriate loss function. Here we have only considered -penalized models, but a natural extension to work with full brain data would be to consider pairwise loss functions combined with sparse structured penalizations which incorporate domain-specific knowledge. This opens the path to further improvements and refinements in the recovery of brain pattern activation via supervised learning.
This work was supported by the ViMAGINE ANR-08-BLAN-0250-02, IRMGroup ANR-10-BLAN-0126-02 and Construct ANR grants.