Functional near-infrared spectroscopy (fNIRS) is an emerging technique for monitoring the concentration changes of oxy- and deoxy-hemoglobin (oxy-Hb and deoxy-Hb) in the brain. An important consideration in fNIRS-based neuroimaging Modality is to conduct group-level analysis from a set of time series measured from a group of subjects. We investigate the feasibility of multilevel statistical inference for NIRS. As a case study, we search for hemodynamic activations in the prefrontal cortex during Stroop interference. Hierarchical general linear model (GLM) is used for making this multilevel analysis. Activation patterns both at the subject and group level are investigated on a comparative basis using various classical and Bayesian inference methods. All methods showed consistent left lateral prefrontal cortex activation for oxy-Hb during interference condition, while the effects were much less pronounced for deoxy-Hb. Our analysis showed that mixed effects or Bayesian models are more convenient for faithful analysis of fNIRS data. We arrived at two important conclusions. First, fNIRS has the capability to identify activations at the group level, and second, the mixed effects or Bayesian model is the appropriate mechanism to pass from subject to group-level inference.