A novel statistical analysis method for functional MRI to overcome the drawbacks of conventional data-driven methods such as the independent component analysis (ICA) is developed. Although ICA has been broadly applied to functional MRI due to its capacity to separate spatially or temporally independent components, the assumption of independence has been challenged by recent studies showing that ICA does not guarantee independence of simultaneously occurring distinct activity patterns in the brain. Instead, sparsity of the signal has been shown to be more promising. This coincides with biological findings such as sparse coding in V1 simple cells, electrophysiological experiment results in the human medial temporal lobe, and etc. The main contribution of this paper is, therefore, a new data driven fMRI analysis that is derived solely based upon the sparsity of the signals. A compressed sensing based data-driven sparse generalized linear model is proposed that enables estimation of spatially adaptive design matrix as well as sparse signal components that represent synchronous, functionally organized and integrated neural hemodynamics. Furthermore, an MDL based model order selection rule is shown to be essential in selecting unknown sparsity level for sparse dictionary learning. Multi-level analysis of fMRI data using data-driven sparse GLM is also investigated. Using simulation and real fMRI experiments, we show that the proposed method can adapt individual variation better compared to the conventional ICA methods.