Successful application of two-dimensional transition metal dichalcogenides in optoelectronic, catalytic, or sensing devices heavily relies on the materials' quality, that is, the thickness uniformity, presence of grain boundaries, and the types and concentrations of point defects. Raman spectroscopy is a powerful and nondestructive tool to probe these factors but the interpretation of the spectra, especially the separation of different contributions, is not straightforward. Comparison to simulated spectra is beneficial, but for defective systems first-principles simulations are often computationally too expensive due to the large sizes of the systems involved. Here, we present a combined first-principles and empirical potential method for simulating Raman spectra of defective materials and apply it to monolayer MoS 2 with random distributions of Mo and S vacancies. We study to what extent the types of vacancies can be distinguished and provide insight into the origin of different evolutions of Raman spectra upon increasing defect concentration. We apply to our simulated spectra the phonon confinement model used in previous experiments to assess defect concentrations, and show that the simplest form of the model is insufficient to fully capture peak shapes, but a good match is obtained when the type of phonon confinement and the full phonon dispersion relation are accounted for.