Pain is subjective, while pain neuroimaging analysis methods need to be as objective as possible. Most neuroimaging studies use a general linear model (GLM) approach which heavily relies on a number of key assumptions. Independent component analysis (ICA), on the other hand, is a data-driven approach and hence significantly reduces need for specific assumptions such as Gaussian distributed residuals and the definition of user-specified design matrix, both required by the GLM. In this paper, we propose a two-level ICA-based method as an attractive alternative, and apply it to the analysis of functional magnetic resonance imaging (fMRI) data acquired during thermal pain stimuli. We identify distinct female brain responses in parts of the "pain matrix", operculum (secondary somatosensory cortex, SII), anterior insular (AI), dorsal anterior cingular cortex (dACC), and default mode network (DMN). We also show that this pain and DMN network have significant correlation with behavioral measures.