In the literature there exist many theories that provide predictions for the rheological properties of concentrated colloidal suspensions, some in excellent agreement with existing experimental data. However, the manner in which hydrodynamic interactions are included differs greatly among the various approaches. Here we incorporate hydrodynamic interactions into the formalism developed previously to account for many body thermodynamic interactions in concentrated suspensions. A conservation equation involving two particle conditional averages of the hydrodynamic functions is derived and solved with different approximations for the functions, including the dilute limit, the rescaling due to Brady, and a lubrication approximation from our previous work. Through our previous calculations without hydrodynamic interactions we have a well characterized approximation for the thermodynamic couplings in concentrated suspensions. The addition of approximations for the hydrodynamic interactions allows prediction for hard sphere suspensions of the low shear viscosity, linear viscoelasticity, long-time self-diffusion coefficient, dichroism, and the nonequilibrium structure for comparison with extensive experimental data. We demonstrate that the form of the nonequilibrium structure is relatively insensitive to the model for the hydrodynamic interactions, while its magnitude is increased by the overall slowing down of relaxations. The calculation of the stress is dominated by the structure near contact at volume fractions approaching random close packing while at intermediate volume fractions it is more sensitive to the model for the interaction functions over the entire range of separations. Comparison of the predictions with a wide range of experimental measurements and simulation results suggests that existing hydrodynamic approximations are of the correct magnitude, while existing thermodynamic approximations are correct in form but underestimate the magnitude of the nonequilibrium structure.