in the context of neuroimaging , machine learning approaches have been used so far to address diagnostic problems , where patients were classified into different groups based on anatomical or functional data . by contrast , in cognitive studies , the standard framework for functional or anatomical brain mapping was based on mass univariate inference procedures .
recently , a new way of analyzing functional neuroimaging data has emerged [ 2 , 3 ] , and it consists in assessing how well behavioral information or cognitive states can be predicted from brain activation images such as those obtained with functional magnetic resonance imaging ( fmri ) .
this approach opens new ways for understanding the mental representation of various perceptual and cognitive parameters , which can be regarded as the study of the corresponding neural code , albeit at a relatively low spatial resolution .
the accuracy of the prediction of the behavioral or cognitive target variable , as well as the spatial layout of predictive regions , can provide valuable information about functional brain organization ; in short , it helps to decode the brain system .
many different pattern recognition and machine leaning methods have been used to extract information from brain images and compare it to the corresponding target . among them , linear discriminant analysis ( lda ) [ 3 , 5 ] , support vector machine ( svm ) [ 69 ] , or regularized prediction [ 10 , 11 ] has been particularly used . the major bottleneck in this kind of analytical framework
is that there are far more features than samples , so that the problem is plagued by the curse of dimensionality , leading to overfitting .
dimension reduction can be used to extract relevant information from the data , the standard approach in functional neuroimaging being feature selection ( e.g. , anova ) [ 3 , 6 , 11 , 12 ] . however , by performing feature selection and parameter estimation separately , such approach is not optimal .
thus , a popular combined selection / estimation scheme , recursive feature elimination , may be used .
however , this approach relies on a specific heuristic , which does not guarantee the optimality of the solution and is particularly costly .
by contrast , there is great interest in sparsity - inducing regularizations , which optimize both simultaneously .   in this paper
, we assume that the code under investigation is about some scalar parameter that characterizes the stimuli , such as a scale / shape parameters but possibly also position , speed ( assuming a 1-d space ) , or cardinality .
thus , we focus on regression problems and defer the generalization to classification to future work .
let us introduce the following predictive linear model : 
						 ( 1)y = xw+b , 
					 where y represents the behavioral variable and ( w , b ) are the parameters to be estimated on a training set .
a vector w   can be seen as an image ; p is the number of features ( or voxels ) , and b   is called the intercept .
each row is a p - dimensional sample , that is , an activation map related to the observation . with
n  p , the estimation of w is ill posed . to cope with the high dimensionality of the data , one can penalize the estimation of w , for example , based on the 2 norm of the weights .
classical regularization schemes have been used in functional neuroimaging , such as the ridge regression , lasso , or elastic net regression . however
, these approaches require the amount of penalization to be fixed beforehand and possibly optimized by cross - validation . to deal with the choice of the amount of penalization
, one can use the bayesian regression techniques , which include the estimation of regularization parameters in the whole estimation procedure .
standard bayesian regularization schemes are based on the fact that a penalization by weighted 2 norm is equivalent to setting the gaussian priors on the weights w : 
						 ( 2)w~(0,a1 ) ,  a = diag  ( 1,  ,p ) ,  i[1,
 ,p ] ,  i+ , 
					 where  is the gaussian distribution and i the precision of the ith feature .
the first one is bayesian ridge regression ( brr ) , which corresponds to the particular case 1
=  = m . by regularizing all the features identically , brr is not well suited when only few features are relevant .
the second classical scheme is automatic relevance determination ( ard ) , which corresponds to the case i  j if i  j. the regularization performed by ard is very adaptive , as all the weights are regularized differently .
however , by regularizing each feature separately , ard is prone to underfitting when the model contains too many regressors   and also suffers from convergence issues .  
these classical bayesian regularization schemes have been used in fmri inverse inference studies [ 10 , 14 , 21 ] . however , these studies used only sparsity as built - in feature selection and do not consider neuroscientific assumptions for improving the regularization ( i.e. , within the design of the matrix a ) .
indeed , due to the intrinsic smoothness of functional neuroimaging data , predictive information is rather encoded in different groups of features sharing similar information .
a potentially more adapted approach is the bayesian regression scheme presented in , which regularizes patterns of voxels differently .
the weights of the model are defined by w = u , where u is a matrix defined as set of spatial patterns ( one pattern by column ) and  are the parameters of the decomposition of w in the basis defined by u. the regularization is controlled through the covariance of  , which is assumed to be diagonal with only m possible different values cov  (  ) = exp  ( 1)i + +exp  ( m)i . the matrices i are diagonal and defined subsets of columns of u sharing similar variance exp  ( i ) . due to its class - based model , this approach is similar to the one proposed in this paper , but the construction of i relies on ad hoc voxel selection steps , so that there is no proof that the solution is correct .
a contrario , the proposed approach jointly optimizes , within the same framework , the construction of the pattern of voxels and the regularization parameter of each pattern .   in this paper , we detail a model for the bayesian regression in which features are grouped into k different classes that are subject to different regularization penalties . the estimation of the penalty is performed in each class separately , leading to a stable and adaptive regularization .
the construction of the group of features and the estimation of the predictive function are performed jointly .
this approach , called multiclass sparse bayesian regression ( mcbr ) , is thus an intermediate solution between brr and ard .
it requires less parameters to estimate than ard and is far more adaptive than brr .
another asset of the proposed approach in fmri inverse inference is that it creates a clustering of the features and thus yields useful maps for brain mapping . after introducing our model and giving some details on the parameter estimation algorithms ( the variational bayes or gibbs sampling procedures ) ,
we show that the proposed algorithm yields better accuracy than reference methods , while providing more interpretable models .
we first detail the notations of the problem and describe the priors and parameters of the model .
we recall the linear model for regression : 
								 ( 3)y = f(x , w , b)=xw+b . 
							
we denote by y   the targets to be predicted and x   the set of activation images related to the presentation of different stimuli .
the integer p is the number of voxels and n the number of samples ( images ) .
typically , p ~ 10  to  10 ( for a whole volume ) , while n ~ 10  to  10 . 

priors on the noisewe use classical priors for regression , and we model the noise on y as an i.i.d .
gaussian variable : 
									 ( 4)~(0,1in),~(;1,2 ) , 
								 where  is the precision parameter and  stands for the gamma density with two hyperparameters 1 , 2 : 
									 ( 5)(x;1,2)=21x11exp  x2(1 ) . 
 
we use classical priors for regression , and we model the noise on y as an i.i.d .
gaussian variable : 
									 ( 4)~(0,1in),~(;1,2 ) , 
								 where  is the precision parameter and  stands for the gamma density with two hyperparameters 1 , 2 : 
									 ( 5)(x;1,2)=21x11exp  x2(1 ) . 

priors on the class assignmentin order to combine the sparsity of ard with the stability of brr , we introduce an intermediate representation , in which each feature j belongs to one class among k indexed by a discrete variable zj ( z = { z1 ,  , zp}).all the features within a class k  { 1 ,  , k } share the same precision parameter k , and we use the following prior on z : 
									 ( 6)z~j=1p  k=1kkjk , 
								 where  is kronecker 's  , defined as 
									 ( 7)jk={0if  zjk,1if  zj = k.we finally introduce an additional dirichlet prior   on  : 
									 ( 8)~dir( ) 
								 with a hyperparameter . by updating at each step the probability k of each class , it is possible to prune classes .
this model has no spatial constraints and thus is not spatially regularized . 
  in order to combine the sparsity of ard with the stability of brr
, we introduce an intermediate representation , in which each feature j belongs to one class among k indexed by a discrete variable zj ( z = { z1 ,  , zp}).all the features within a class k  { 1 ,  , k } share the same precision parameter k , and we use the following prior on z : 
									 ( 6)z~j=1p  k=1kkjk , 
								 where  is kronecker 's  , defined as 
									 ( 7)jk={0if  zjk,1if  zj = k .
we finally introduce an additional dirichlet prior   on  : 
									 ( 8)~dir( ) 
								 with a hyperparameter . by updating at each step the probability k of each class , it is possible to prune classes .
priors on the weightsas in ard , we make use of an independent gaussian prior for the weights : 
									 ( 9)w~(0,a1 )  with  diag  ( a)={z1,  ,zp } , 
								 where zj is the precision parameter of the jth feature , with zj  { 1 ,  , k}. we introduce the following prior on k : 
									 ( 10)k~(k;1,k,2,k ) 
								 with hyperparameters 1,k , 2,k .
, we make use of an independent gaussian prior for the weights : 
									 ( 9)w~(0,a1 )  with  diag  ( a)={z1,  ,zp } , 
								 where zj is the precision parameter of the jth feature , with zj  { 1 ,  , k}. we introduce the following prior on k : 
									 ( 10)k~(k;1,k,2,k ) 
								 with hyperparameters 1,k , 2,k .
the link between the proposed mcbr model and the other regularization methods , bayesian ridge regression and automatic relevance determination , is obvious . with k = 1 ,
=  = zp , we retrieve the brr model ,   with k = p , that is ,   zi  zj if i  j , and assigning each feature to a singleton class ( i.e. , zj = j ) , we retrieve the ard model .   moreover , the proposed approach is related to the one developed in . in this paper , the authors proposed , for the distribution of weights of the features , a binary mixture of gaussians with small and large precisions .
our work can be viewed as a generalization of this model to a number of classes k  2 . for models with latent variables , such as mcbr
for instance in a mixture of components , a singularity is a component with one single sample and thus zero variance . in such cases , maximizing the log likelihood yields flawed solutions , and one can use the posterior distribution of the latent variables p(z | x , y ) for this maximization .
however , the posterior distribution of the latent variables given the data does not have a closed - form expression , and some specific estimation methods , such as variational bayes or gibbs sampling , have to be used .
we first estimate the model by the variational bayes , and the resulting algorithm is thus called vb - mcbr .
we also detail an algorithm , called gibbs - mcbr , based on a gibbs sampling procedure .                                    
the variational bayes ( or vb ) approach provides an approximation q( ) of p( | y ) , where q( ) is taken in a given family of distributions and  = [ w ,  ,  , z ,  ] .
additionally , the variational bayes approach often uses the following mean field approximation , which allows the factorization between the approximate distribution of the latent variables and the approximate distributions of the parameters : 
 ( 11)q()=q(w)q()q()q(z)q( ) . 
  
we introduce the kullback - leibler divergence (q( ) ) that measures the similarity between the true posterior p( | y ) and the variational approximation q( ) .
( y ) as 
										 ( 12)log  p  ( y    )=(q())+(q( ) ) 
									 with 
										 ( 13)(q())=dq()log  p(y,)q(),(q())=  dq()log  q()p(    y ) , 
									 where (q( ) ) is called free energy and can be seen as the measure of the quality of the model . as (q( ) )  0 , the free energy is a lower bound on log  p  ( y ) with equality if and only if q( ) = p( | y ) .
so , inferring the density q( ) of the parameters corresponds to maximizing  with respect to the free distribution q( ) . in practice ,
the vb approach consists in maximizing the free energy  iteratively with respect to the approximate distribution q(z ) of the latent variables and with respect to the approximate distributions of the parameters of the model q(w ) , q( ) , q( ) , and q( ) .
the variational distributions and the pseudocode of the vb - mcbr algorithm are provided in appendix a. this algorithm maximizes the free energy . in practice , iterations are performed until convergence to a local maximum of . with an ard prior ( i.e. , k = p and fixing zj = j ) , we retrieve the same formulas as the ones found for variational ard .
the resulting algorithm is called gibbs - mcbr ; the pseudocode of the algorithm and the candidate distributions are provided in appendix b. the gibbs sampling algorithm is used for generating a sequence of samples from the joint distribution to approximate marginal distributions .
the main idea is to use conditional distributions that should be known and possibly easy to sample from , instead of directly computing the marginals from the joint law by integration ( the joint law may not be known or hard to sample from ) .
the sampling is done iteratively among the different parameters , and the final estimation of parameters is obtained by averaging the values of the different parameters across the different iterations ( one may not consider the first iterations , this is called the burn in ) . our model needs few hyperparameters ; we choose here to use slightly informative and class - specific hyperparameters in order to reflect a wide range of possible behaviors for the weight distribution .
this choice of priors is equivalent to setting heavy - tailed centered student 's t - distributions with variance at different scales , as priors on the weight parameters .
we set k = 9 , with weakly informative priors 1,k = 10 , k  [ 1 ,  , k ] and 2,k = 10 , k  [ 1 ,  , k ] .
moreover , we set 1 = 2 = 1 . starting with a given number of classes and letting the model automatically prune the classes can be seen as a means of avoiding costly model selection procedures .
the choice of class - specific priors is also useful to avoid label switching issues and thus speeds up convergence .
crucially , the priors used here can be used in any regression problem , provided that the target data is approximately scaled to the range of values used in our experiments . in that sense , the present choice of priors can be considered as universal .
we also randomly initialize q(z ) for vb - mcbr ( or z for gibbs - mcbr ) .  
our method is evaluated with a cross - validation procedure that splits the available data into training and validation sets . in the following , ( x , y ) are a learning set ( x , y ) is a test set , and y^t = f(xtw^ ) refers to the predicted target , where w^ is estimated from the training set . the performance of the different models is evaluated using  , the ratio of explained variance : 
										 ( 14)(yt , y^t)=var(yt)var(yty^t)var(yt ) . 
									 this is the amount of variability in the response that can be explained by the model ( perfect prediction yields  = 1 , while  < 0 if prediction is worse than chance ) .   in our experiments ,
the proposed algorithms are compared to different state - of - the - art regularization methods .  

elastic net   regression , which requires setting two parameters 1 and 2 . in our analyzes , a cross - validation procedure within the training set is used to optimize these parameters . here
, we use 1{0.2,0.1,0.05,0.01 } , where =||xty|| , and 2  { 0.1,0.5,1 . , 10 .
support vector regression ( svr ) with a linear kernel , which is the reference method in neuroimaging .
the c parameter is optimized by cross - validation in the range of 10 to 10 in multiplicative steps of 10 .  

bayesian ridge regression ( brr ) , which is equivalent to mcbr with k = 1 and 1 = 2 = 1 = 2 = 10 , that is , weakly informative priors . 
 automatic relevance determination ( ard ) , which is equivalent to mcbr with k = p and 1 = 2 = 1 = 2 = 10 , that is , weakly informative priors .   
all these methods are used after an anova - based feature selection as this maximizes their performance . indeed , irrelevant features and redundant information
the optimal number of voxels is selected within the range { 50,100,250,500 } , using a nested cross - validation within the training set .
we do not directly select a threshold on p value or cluster size , but rather a predefined number of features .
the estimation of the parameters of the learning function is also performed using a nested cross - validation within the training set , to ensure a correct validation and an unbiased comparison of the methods .
the implementation of elastic net is based on coordinate descent , while svr is based on libsvm .
methods are used from python via the scikit - learn open source package .             for vb - mcbr and gibbs - mcbr , in order to avoid a costly internal cross - validation ,
the number of iterations used is fixed to 5000 ( burn in of 4000 iterations ) for gibbs - mcbr and 500 for vb - mcbr .
preliminary results on both simulated and real data showed that these values are sufficient enough for an accurate inference of the model .
as explained previously , we set k = 9 , with weakly informative priors 1,k = 10 ,   k  [ 1 ,  , k ] and 2,k = 10 , k  [ 1 ,  , k ] .
moreover , we set 1 = 2 = 1 , and we randomly initialize q(z ) for vb - mcbr ( or z for gibbs - mcbr ) .
  we now evaluate and illustrate mcbr on two different sets of simulated data .  
we first test mcbr on a simulated data set , designed for the study of ill - posed regression problem , that is , n  p. data are simulated as follows : 
										 ( 15)x~(0,1 )  with  ~(0,1),y=2(x1+x2x3x4)+0.5(x5+x6x7x8)+. 
									 we have p = 200 features , n = 50 images for the training set , and n = 50 images for the test set .
we compare mcbr to the reference methods , but we do not use feature selection , as the number of features is not very high .
we average the results of 15 different trials , and the average explained variance is shown in table 1 .
gibbs - mcbr outperforms the other approaches , yielding higher prediction accuracy than the reference elastic net and ard methods .
vb - mcbr falls into the local maximum of  and does not yield an accurate prediction .
indeed , it can not finely adapt the weights of the relevant features , as these features are regularized similarly as the irrelevant ones .
svr has also low accuracy , due to the fact that we do not perform any feature selection .
thus , svr suffers from the curse of dimensionality , unlike other methods such as ard or elastic net , which performs feature selection and model estimation jointly . in figure 2
, we represent the probability density function of the distributions of the weights obtained with brr ( a ) , gibbs - mcbr ( b ) , and ard ( c ) . with brr
the gibbs - mcbr algorithm creates a multimodal distribution , lots of weights being highly regularized ( pink distributions ) , and informative features are allowed to have higher weights ( blue distributions ) . with mcbr ,
weights are clustered into different groups , depending on their predictive power , which is interesting in application such as fmri inverse inference , as it can yield more interpretable models .
indeed , the class to the features with higher weights ( { x1 , x2 , x3 , x4 } ) belong which is small ( average size of 6 features ) but has a high purity ( percentage of relevant features in the class ) of 74% .
we now look at the values of w1 and w2 for the different steps of the two algorithms ( see figure 3 ) .
we can see that vb - mcbr ( b ) quickly falls into a local maximum , while gibbs - mcbr ( a ) visits the space and reaches the region of the correct set of parameters ( red dot ) .
the simulated data set x consists of n = 100 images ( size 12  12  12 voxels ) with a set of four square regions of interest ( roi ) ( size 2  2  2 ) .
we call  the support of the roi ( i.e. , the 32 resulting voxels of interest ) .
each of the four rois has a fixed weight in { 0.5,0.5 , 0.5,0.5}. we call wi , j , k the weight of the ( i , j , k ) voxel .
the resulting images are smoothed with a gaussian kernel with a standard deviation of 2 voxels , to mimic the correlation structure observed in real fmri data . to simulate the spatial variability between images ( intersubject variability , movement artifacts in intrasubject variability ) , we define a new support of the rois , called  such that , for each image lth , 50% ( randomly chosen ) of the weights w are set to zero . thus , we have . we simulate the target y for the lth image as  
 ( 16)yl=(i , j , k)wi , j , kxi , j , k , l+l 
									 with the signal in the ( i , j , k ) voxel of the lth image simulated as 
 ( 17)xi , j , k , l~(0,1 ) , 
									 and l ~ (0 ,  ) is a gaussian noise with standard deviation  > 0 .
we choose  in order to have a signal - to - noise ratio of 5  db .
the resulting images of weights are given in figure 4 , with the true weights ( a ) and resulting anova f - scores ( b ) .
the reference methods can detect the truly informative regions ( rois ) , but elastic net ( f ) and ard ( h ) retrieve only part of the support of the weights .
moreover , elastic net yields an overly sparse solution . brr ( g ) also retrieves the rois but does not yield a sparse solution , as all the features are regularized in the same way .
we note that the weights in the feature space estimated by svr ( e ) are nonzero everywhere and do not outline the support of the ground truth .
vb - mcbr ( c ) converges to a local maximum similar to the solution found by brr ( g ) ; that is , it creates only one nonempty class , and thus regularizes all the features similarly .
we can thus clearly see that , in this model , the variational bayes approach is very sensitive to the initialization and can fall into nonoptimal local maxima , for very sparse support of the weights .
finally , gibbs - mcbr ( d ) retrieves most of the true support of the weights by performing an adapted regularization . in this section ,
we assess the performance of mcbr in an experiment on the mental representation of object size , where the aim is to predict the size of an object seen by the subject during the experiment , in both intrasubject and intersubject cases .
the size ( or scale parameter ) of the object will be the target variable y. we apply the different methods on a real fmri dataset related to an experiment studying the representation of objects , on ten subjects , as detailed in . during this experiment , ten healthy volunteers viewed objects of 4 shapes in 3 different sizes ( yielding 12 different experimental conditions ) , with 4 repetitions of each stimulus in each of the 6 sessions .
we pooled data from the 4 repetitions , resulting in a total of n = 72 images by subject ( one image of each stimulus by session ) .
functional images were acquired on a 3-t mr system with an eight - channel head coil ( siemens trio , erlangen , germany ) as t2 * -weighted echo - planar image ( epi ) volumes .
twenty transverse slices were obtained with a repetition time of 2  s ( echo time : 30  ms ; flip angle : 70 ; 2  2  2-mm voxels ; 0.5  mm gap ) .
realignment , normalization to mni space , and general linear model ( glm ) fit were performed with the spm5 software ( http://www.fil.ion.ucl.ac.uk/spm/software/spm5/ ) .
the normalization is the conventional method of spm ( implying affine and nonlinear transformations ) and not the one using unified segmentation .
the normalization parameters are estimated on the basis of a whole - head epi acquired in addition and are then applied to the partial epi volumes .
the effect of each of the 12 stimuli convolved with a standard hemodynamic response function was modeled separately , while accounting for serial autocorrelation with an ar(1 ) model and removing low - frequency drift terms using a high - pass filter with a cutoff of 128  s. the glm is fitted separately in each session for each subject , and we used in the present work the resulting session - wise parameter estimate images ( the -maps are used as rows of x ) .
the four different shapes of objects were pooled across for each one of the three sizes , and we are interested in finding discriminative information on sizes .
this reduces to a regression problem , in which our goal is to predict a simple scalar factor ( size of an object ) .
all the analyzes are performed without any prior selection of regions of interest and use the whole acquired volume . 

the dimensions of the real data set for one subject are p ~ 7  10 and n = 72 ( divided in 3 different sizes , 24 images per size ) . we evaluate the performance of the method by a leave - one - condition - out cross - validation ( i.e. , leave-6-image - out ) , and doing so the glm is performed separately for the training and test sets .
the parameters of the reference methods are optimized with a nested leave - one - condition - out cross - validation within the training set , in the ranges given before . 
 
the dimensions of the real data set for one subject are p ~ 7  10 and n = 72 ( divided in 3 different sizes , 24 images per size ) . we evaluate the performance of the method by a leave - one - condition - out cross - validation ( i.e. , leave-6-image - out ) , and doing so the glm is performed separately for the training and test sets .
the parameters of the reference methods are optimized with a nested leave - one - condition - out cross - validation within the training set , in the ranges given before . 

the intersubject analysis relies on subject - specific fixed - effect activations that is , for each condition , the 6 activation maps corresponding to the 6 sessions are averaged together .
this yields a total of 12 images per subject , one for each experimental condition .
the dimensions of the real data set are p ~ 7  10 and n = 120 ( divided into 3 different sizes ) .
we evaluate the performance of the method by cross - validation ( leave - one - subject - out ) .
the parameters of the reference methods are optimized with a nested leave - one - subject - out cross - validation within the training set , in the ranges given before . 
  additionally , we perform an intersubject regression analysis on the sizes .
the intersubject analysis relies on subject - specific fixed - effect activations that is , for each condition , the 6 activation maps corresponding to the 6 sessions are averaged together .
this yields a total of 12 images per subject , one for each experimental condition .
the dimensions of the real data set are p ~ 7  10 and n = 120 ( divided into 3 different sizes ) .
we evaluate the performance of the method by cross - validation ( leave - one - subject - out ) .
the parameters of the reference methods are optimized with a nested leave - one - subject - out cross - validation within the training set , in the ranges given before . 

compared to the results on simulated data , vb - mcbr still falls in a local maximum similar to the bayesian ridge regression which performs well in this experiment .
moreover , both gibbs - mcbr and vb - mcbr are more stable than the reference methods . 
 
compared to the results on simulated data , vb - mcbr still falls in a local maximum similar to the bayesian ridge regression which performs well in this experiment .
moreover , both gibbs - mcbr and vb - mcbr are more stable than the reference methods . 

as in the intrasubject analysis , both mcbr approaches outperform the reference methods , svr , brr , and ard .
however , the prediction accuracy is similar to that of elastic net . in this case ,
gibbs - mcbr performs slightly better than vb - mcbr , but the difference is not significant.one major asset of mcbr ( and more particularly gibbs - mcbr , as vb - mcbr often falls into a one - class local maximum ) is that it creates a clustering of the features , based on the relevance of the features in the predictive model .
this clustering can be accessed using the variable z , which is implied in the regularization performed on the different features . in figure 5
, we give the histogram of the weights of gibbs - mcbr for the intersubject analysis .
we keep the weights and the values of z of the last iteration ; the different classes are represented as dots of different colors and are superimposed on the histogram .
we can notice than the pink distribution represented at the bottom of the histogram corresponds to relevant features .
this cluster is very small ( 19  voxels ) , compared to the two blue classes represented at the top of the histogram that contain many voxels ( 746  voxels ) which are highly regularized , as they are noninformative.the maps of weights found by the different methods are detailed in figure 6 .
the methods are used combined with an anova - based univariate feature selection ( 2500  voxels selected , in order to have a good support of the weights ) . as elastic net , gibbs - mcbr yields a sparse solution but extracts a few more voxels .
the map found by elastic net is not easy to interpret , with very few informative voxels scattered in the whole occipital cortex .
the map found by svr is not sparse in the feature space and is thus difficult to interpret , as the spatial layout of the neural code is not clearly extracted .
vb - mcbr does not yield a sparse map either , all the features having nonnull weights 
  the results obtained with the different methods are given in table 3 .
as in the intrasubject analysis , both mcbr approaches outperform the reference methods , svr , brr , and ard . however , the prediction accuracy is similar to that of elastic net . in this case ,
gibbs - mcbr performs slightly better than vb - mcbr , but the difference is not significant
. one major asset of mcbr ( and more particularly gibbs - mcbr , as vb - mcbr often falls into a one - class local maximum ) is that it creates a clustering of the features , based on the relevance of the features in the predictive model .
this clustering can be accessed using the variable z , which is implied in the regularization performed on the different features . in figure 5 , we give the histogram of the weights of gibbs - mcbr for the intersubject analysis .
we keep the weights and the values of z of the last iteration ; the different classes are represented as dots of different colors and are superimposed on the histogram .
we can notice than the pink distribution represented at the bottom of the histogram corresponds to relevant features .
this cluster is very small ( 19  voxels ) , compared to the two blue classes represented at the top of the histogram that contain many voxels ( 746  voxels ) which are highly regularized , as they are noninformative . the maps of weights found by the different methods are detailed in figure 6 .
the methods are used combined with an anova - based univariate feature selection ( 2500  voxels selected , in order to have a good support of the weights ) . as elastic net , gibbs - mcbr yields a sparse solution but extracts a few more voxels .
the map found by elastic net is not easy to interpret , with very few informative voxels scattered in the whole occipital cortex .
the map found by svr is not sparse in the feature space and is thus difficult to interpret , as the spatial layout of the neural code is not clearly extracted .
vb - mcbr does not yield a sparse map either , all the features having nonnull weights
it is well known that in high - dimensional problems , regularization of feature loadings significantly increases the generalization ability of the predictive model . however , this regularization has to be adapted to each particular dataset . in place of costly cross - validation procedures
, we cast regularization in a bayesian framework and treat the regularization weights as hyperparameters .
the proposed approach yields an adaptive and efficient regularization and can be seen as a compromise between a global regularization ( bayesian ridge regression ) that does not take into account the sparse or focal distribution of the information and automatic relevance determination .
additionally , mcbr creates a clustering of the features based on their relevance and thus explicitly extracts groups of informative features .
ard is subject to an underfitting in the hyperparameter space that corresponds to an underfitting in model selection ( i.e. , on the features to be pruned ) . indeed , as ard is estimated by maximizing evidence , models with less selected features are preferred , as the integration is done on less dimensions , and thus evidence is higher .
a contrario , mcbr requires far less hyperparameter ( 2  k , with k  p ) and suffers less from this issue , as the sparsity of the model is defined by groups .
moreover , a full bayesian framework for estimating ard requires to set some priors on the hyperparameters ( e.g. , 1 and 2 ) , and it may be sensitive to specific choice of these hyperparameters . a solution is to use an internal cross - validation for optimizing these parameters , but this approach can be computationally expensive . in the case of mcbr , the distributions of the hyperparameters are bound to a class and not to each feature .
indeed , the choice of good hyperparameters for the features is dealt with at the class level . on simulated data ,
our approach performs better than other classical methods such as svr , brr , ard , and elastic net and yields a more stable prediction accuracy .
moreover , by adapting the regularization to different groups of voxels , mcbr retrieves the true support of the weights and recovers a sparse solution .
results on real data show that mcbr yields more accurate predictions than other regularization methods .
as it yields less sparse solution than elastic net , it gives access to more plausible loading maps which are necessary for understanding the spatial organization of brain activity , that is , retrieving the spatial layout of the neural coding . on real fmri data ,
the explicit clustering of gibbs - mcbr is also an interesting aspect of the model , as it can extract few groups of relevant features from many voxels . in some experiments ,
the variational bayes algorithm yields less accurate predictions than the gibbs sampling approach , which can be explained by the difficulty of initializing the different variables ( especially z ) when the support of the weight is overly sparse .
moreover , the vb - mcbr algorithm relies on a variational bayes approach , which may not be optimal , due to strong approximations in model inference .
a contrario gibbs - mcbr is more time consuming but yields a better model inference .
finally , the variability in the results may be explained by the difficulty to estimate the model ( optimality is not ensured ) . the question of model selection ( i.e. , the number of classes k ) has not been addressed in this paper .
one can use the free energy in order to select the best model , but due to the instability of vb - mcbr , this approach does not seem promising .
a more interesting method is the one detailed in , which can be used with the gibbs sampling algorithm . here
, model selection is performed implicitly by emptying classes that do not fit the data well . in that respect ,
the choice of heterogeneous priors for different classes is crucial : replacing our priors with class - independent priors ( i.e. , 1,k = 10 , k  [ 1 ,  , k ] ) in the intersubject analysis on size prediction leads gibbs - mcbr to a local maximum similar to vb - mcbr .
finally , this model is not restricted to the bayesian regularization and can be used for classification , within a probit or logit model [ 35 , 36 ] .
the proposed model may thus be used for diagnosis in medical imaging , for the prediction of both continuous or discrete variables .
in this paper , we have proposed a model for adaptive regression , called mcbr .
the proposed method integrates , in the same bayesian framework , brr and ard and performs a different regularization for relevant and irrelevant features .
it can tune the regularization to the possible different level of sparsity encountered in fmri data analysis , and it yields interpretable information for fmri inverse inference , namely , the z variable ( latent class variable ) .
experiments on both simulated and real data show that our approach is well suited for neuroimaging , as it yields accurate and stable predictions compared to the state - of - the - art methods .
the variational bayes approach yields the following variational distributions :  
 q(w ) ~ (w |  ,  ) with 
										 ( a.1)a=diag  ( l1,  ,lp )  withlj=k=1kq(zj = k)l1,kl2,k  j{1,  ,p},(a.2)=(a1a2xtx+a)1,(a.3)=a1a2xty ; 
 q(k ) ~ (l1,k , l2,k ) with 
										 ( a.4)l1,k=1,k+12j=1pq(zj = k),(a.5)l2,k=2,k+12j=1p(jj2+jj)q(zj = k ) ; 
 q( ) ~ (a1 , a2 ) with 
										 ( a.6)a1=1+n2,(a.7)a2=2 + 12(yx)t(yx)+12tr  ( xtx ) ; 
 q(zj = k ) ~ exp   with 
										 ( a.8)jk=12(j2+jj)l1,kl2,k+ln  ( k)+12((l1,k)log  ( l2,k)),(a.9)k = exp  { (dk)(k=1k = kdk)},(a.10)dk=k+j=1pq(zj = k ) , 
									 where  is the digamma function (x ) = (x)/(x ) . the vb - mcbr algorithm is provided in pseudo - code in algorithm 1 .
with  = [ w ,  ,  , z ,  ] , we have the following candidate distributions ( i.e. , the distributions used for the sampling of the different parameters ) : 
 p(w |   { w } )  (w |  ,  ) with 
										 ( b.1)=(xtx+a)1  with  a = diag  ( z1,  ,zp),(b.2)=xty ; 
 p( |   {  } )  k=1(k | l1,k , l2,k ) with 
										 ( b.3)l1,k=1,k+12j=1p(zj = k),(b.4)l2,k=2,k+12j=1p(zj = k)wj2 ; 
 p( |   {  } )  (a1 , a2 ) with 
										 ( b.5)a1=1+n2,(b.6)a2=2 + 12(yx)t(yx ) ; 
 p(zj |   { z } )  mult(exp  j,1 ,  , exp  j , k ) with 
										 ( b.7)jk=12wj2k+ln  ( k)+12log  k ; 
 p(k |   {  } )  dir(dk ) with 
										 ( b.8)dk=k+j=1p(zj = k ) . the algorithm is provided in pseudocode in algorithm 2 .
 
 inverse inference has recently become a popular approach for analyzing neuroimaging data , by quantifying the amount of information contained in brain images on perceptual , cognitive , and behavioral parameters . 
 as it outlines brain regions that convey information for an accurate prediction of the parameter of interest , it allows to understand how the corresponding information is encoded in the brain . 
 however , it relies on a prediction function that is plagued by the curse of dimensionality , as there are far more features ( voxels ) than samples ( images ) , and dimension reduction is thus a mandatory step . 
 we introduce in this paper a new model , called multiclass sparse bayesian regression ( mcbr ) , that , unlike classical alternatives , automatically adapts the amount of regularization to the available data . 
 mcbr consists in grouping features into several classes and then regularizing each class differently in order to apply an adaptive and efficient regularization . 
 we detail these framework and validate our algorithm on simulated and real neuroimaging data sets , showing that it performs better than reference methods while yielding interpretable clusters of features . 