Multiple Sclerosis (MS) is a chronic neuroinflammatory disorder of the brain and spinal cord, widely studied nowadays, due to its relevant prevalence in the population. Even though no cure exists, an earlier and more adequate choice of treatment could delay its evolution and prevent irreversible sequelae and disability progression. Currently, Magnetic Resonance Imaging (MRI) represents an essential nonclinical tool for the detection of a hallmark of the disease, i.e. the presence of demyelinating lesions within cerebral white matter (WM), and, consequently, for the diagnosis of MS early within its course. However, errors in estimating lesions can contribute to a wrong diagnosis, if only the WM lesion load is taken into account, with a more relevant impact in individuals with a reduced lesion load at an initial clinical event, delaying the start of a treatment until a second clinical relapse or after confirming, successively, dissemination of lesions in time. In this context, this work proposes an innovative system, employing a multivariate analysis approach, with the aim of mining and integrating multiple sensitive neuroimaging markers, including but not limited to the WM lesion load, into classification models for supporting a more robust diagnosis of Relapsing-Remitting-MS (RR-MS) already at an initial clinical event. To this aim, a retrospective study of 81 patients with diagnosis of RR-MS (39 males and 42 females, 37.3 ± 8.1 years old, age range 20 - 58) and 29 healthy people of comparable age and gender (15 males and 14 females, 39.7 ± 11.1 years old, age range 22 - 57) is used. All the individuals were enrolled at Multiple Sclerosis Centre of the "Federico II" University Hospital (Naples, Italy). A machine learning method based on both statistics and Fuzzy Logic, already validated for its desirable characteristics, is applied to volumetric and relaxometric data extracted from brain MRI by means of a multiparametric relaxometric segmentation method. As results of this work, one- and multi-dimensional models are generated and compared with those obtained with other stateof-the-art methods, confirming their goodness in terms of performance, interpretability and robustness, even if their generalizability is limited by the dependence of measurements on MRI scanner. More in particular, the most promising model is two-dimensional and presents a high value of accuracy (89.1%), a particularly high sensitivity (93.8%) and satisfactory negative predictive value (81.5%) and positive predictive value (91.6%). This model identifies a further imaging marker related to grey matter alterations that, combined to the WM lesion load, allows to discriminate more accurately between RR-MS patients and healthy subjects, perform a more robust detection of RR-MS patients at the disease onset and, contextually, confidently exclude the disease presence in subjects classified as healthy. In addition, due to its interpretability, this model is able to clarify how and how much WM lesion load and grey matter alterations are correlated to RR-MS. Finally, it allows personalized diagnosis support, by calculating reliable classification confidence associated to each particular case.