  We provide an end-to-end exploration of a distinct modified gravitational
theory in Jordan-Brans-Dicke (JBD) gravity, from an analytical and numerical
description of the background expansion and linear perturbations, to the
nonlinear regime captured with a hybrid suite of $N$-body simulations, to the
parameter constraints from existing cosmological probes. The nonlinear
corrections to the matter power spectrum due to baryons, massive neutrinos, and
modified gravity are simultaneously modeled and propagated in the cosmological
analysis for the first time. In the combined analysis of the Planck CMB
temperature, polarization, and lensing reconstruction, Pantheon supernova
distances, BOSS measurements of BAO distances, the Alcock-Paczynski effect, and
the growth rate, along with the joint ($3\times2$pt) dataset of cosmic shear,
galaxy-galaxy lensing, and overlapping redshift-space galaxy clustering from
KiDS and 2dFLenS, we constrain the JBD coupling constant, $\omega_{\rm
BD}>1540$ (95% CL), the effective gravitational constant, $G_{\rm
matter}/G=0.997\pm0.029$, the sum of neutrino masses, $\sum m_{\nu}<0.12$ eV
(95% CL), and the baryonic feedback amplitude, $B<2.8$ (95% CL), all in
agreement with the standard model expectation. We show that the uncertainty in
the gravitational theory alleviates the tension between KiDS$\times$2dFLenS and
Planck to below $1\sigma$ and the tension in the Hubble constant between Planck
and the direct measurement of Riess et al. (2019) down to ~$3\sigma$; however,
we find no substantial model selection preference for JBD gravity relative to
$\Lambda$CDM. We further show that the neutrino mass bound degrades by up to a
factor of three as the $\omega_{\rm BD}$ parameterization becomes more
restrictive and that a positive shift in $G_{\rm matter}/G$ suppresses the CMB
damping tail in a way that might complicate future inferences of small-scale
physics. (Abridged)
