Massive stars play an important role in many areas of astrophysics, but numerous details regarding their formation remain unclear. In this paper we present and analyse high resolution (R ~ 30,000) near-infrared 2.3 micron spectra of 20 massive young stellar objects from the RMS database, in the largest such study of CO first overtone bandhead emission to date. We fit the emission under the assumption it originates from a circumstellar disc in Keplerian rotation. We explore three approaches to modelling the physical conditions within the disc - a disc heated mainly via irradiation from the central star, a disc heated mainly via viscosity, and a disc in which the temperature and density are described analytically. We find that the models described by heating mechanisms are inappropriate because they do not provide good fits to the CO emission spectra. We therefore restrict our analysis to the analytic model, and obtain good fits to all objects that possess sufficiently strong CO emission, suggesting circumstellar discs are the source of this emission. On average, the temperature and density structure of the discs correspond to geometrically thin discs, spread across a wide range of inclinations. Essentially all the discs are located within the dust sublimation radius, providing strong evidence that the CO emission originates close to the central protostar, on astronomical unit scales. In addition, we show that the objects in our sample appear no different to the general population of MYSOs in the RMS database, based on their near- and mid-infrared colours. The combination of observations of a large sample of MYSOs with CO bandhead emission and our detailed modelling provide compelling evidence of the presence of small scale gaseous discs around such objects, supporting the scenario in which massive stars form via disc accretion.