Magnetic resonance imaging (MRI) plays a vital role in the scientific investigation and clinical management of multiple sclerosis. Analyses of binary multiple sclerosis lesion maps are typically "mass univariate" and conducted with standard linear models that are ill suited to the binary nature of the data and ignore the spatial dependence between nearby voxels (volume elements). Smoothing the lesion maps does not entirely eliminate the non-Gaussian nature of the data and requires an arbitrary choice of the smoothing parameter. Here we present a Bayesian spatial model to accurately model binary lesion maps and to determine if there is spatial dependence between lesion location and subject specific covariates such as MS subtype, age, gender, disease duration and disease severity measures. We apply our model to binary lesion maps derived from T2-weighted MRI images from 250 multiple sclerosis patients classified into five clinical subtypes, and demonstrate unique modeling and predictive capabilities over existing methods.