Bayesian approaches have been proposed by several functional magnetic resonance imaging (fMRI) researchers in order to overcome the fundamental limitations of the popular statistical parametric mapping method. However, the difficulties associated with subjective prior elicitation have prevented the widespread adoption of the Bayesian methodology by the neuroimaging community. In this paper, we present a Bayesian multilevel model for the analysis of brain fMRI data. The main idea is to consider that all the estimated group effects (fMRI activation patterns) are exchangeable. This means that all the collected voxel time series are considered manifestations of a few common underlying phenomena. In contradistinction to other Bayesian approaches, we think of the estimated activations as multivariate random draws from the same distribution without imposing specific prior spatial and/or temporal information for the interaction between voxels. Instead, a two-stage empirical Bayes prior approach is used to relate voxel regression equations through correlations between the regression coefficient vectors. The adaptive shrinkage properties of the Bayesian multilevel methodology are exploited to deal with spatial variations, and noise outliers. The characteristics of the proposed model are evaluated by considering its application to two real data sets.