Parkinson’s disease (PD) is a progressive neurodegenerative movement disorder characterized by a large number of motor and non-motor deficits, which significantly contribute to reduced quality of life. Despite the definition of the broad spectrum of clinical characteristics, mechanisms triggering illness, the nature of its progression and a character of therapeutic effects still remain unknown. The enormous advances in magnetic resonance imaging (MRI) in the last decades have significantly affected the research attempts to uncover the functional and structural abnormalities in PD and have helped to develop and monitor various treatment strategies, of which dopamine replacement strategies, mainly in form of levodopa, has been the gold standard since the late seventies and eighties. Motor, task-related functional MRI (fMRI) has been extensively used to assess the pathological state of the motor circuitry in PD. Several studies employed motor paradigms and fMRI to review the functional brain responses of participants to levodopa treatment. Interestingly, they provided conflicting results. Wide spectrum of symptoms, variability and asymmetry of the disease presentation, several treatment approaches and their divergent outcomes make PD enormously heterogeneous. In this work we hypothesized that not considering the disease heterogeneity might have been an adequate cause for the discrepant results in aforementioned studies. We show that not accounting for the disease variability might indeed compromise the results and invalidate the consequent interpretations. Accordingly, we propose and formalize a statistical approach to account for the intra and inter subject variability. This might help to minimize this bias in future motor fMRI studies revealing the functional brain dysfunction and contribute to the understanding of still unknown pathophysiological mechanisms underlying PD.