The detection and characterization of exoplanets—planets that orbit stars other than the Sun—is one of the most active areas of research in modern astronomy. Many exoplanets are discovered using the radial velocity technique, which involves detecting the Doppler shift in a star’s spectral lines resulting from the gravitational effects of an orbiting planet. A challenge to this approach is that measured radial velocity signals are often corrupted by stellar activity such as spots rotating across the star’s surface. A principled method for recovering the underlying planetary radial velocity signal was proposed by Rajpaul et al. (2015), which uses dependent Gaussian processes to jointly model the corrupted radial velocity signal and proxies for stellar activity. Our approach extends Rajpaul et al. (2015) by (i) incorporating science- and data-driven dimension reduction techniques to automatically extract informative stellar activity proxies, and (ii) introducing a model comparison procedure to select the best model for the stellar activity proxies at hand from a larger class of models.