Over the last decade, significant amounts of high-spectral and time-resolution spectroscopic data have been acquired for a number of rapidly oscillating Ap stars (roAp). Progress in the understanding of the information held by these data requires the development of theoretical models that can be directly compared with them. In this work, we present a theoretical model for the radial velocities of roAp stars that takes full account of the coupling between the pulsations and the magnetic field. We explore the impact on the radial velocities of changing the position of the observer, the mode frequency, and angular degree, as well as of changing the region of the disc where the elements are concentrated. We find that for integrations over the full disc, in the outermost layers the radial velocity is generally dominated by the acoustic waves, showing a rapid increase in amplitude. The most significant depth-variations in the radial velocity phase are seen for observers directed towards the equator and for even degree modes with frequencies close to, or above the acoustic cutoff. Comparison between the radial velocities obtained for spots of elements located around the magnetic poles and around the magnetic equator, shows that these present distinct amplitude-phase relations, resembling some of the differences seen in the observations. Finally, we discuss the conditions under which one may expect to find false nodes in the pulsation radial velocity of roAp stars.