We study the propagation of cosmological gravitational wave (GW) backgrounds from the early radiation era until the present day in modified theories of gravity. Comparing to general relativity (GR), we study the effects that modified gravity parameters, such as the GW friction α M and the tensor speed excess α T, have on the present-day GW spectrum. We use both the WKB estimate, which provides an analytical description but fails at superhorizon scales, and numerical simulations that allow us to go beyond the WKB approximation. We show that a constant α T makes relatively insignificant changes to the GR solution, especially taking into account the constraints on its value from GW observations by the LIGO-Virgo collaboration, while α M can introduce modifications to the spectral slopes of the GW energy spectrum in the low-frequency regime depending on the considered time evolution of α M. The latter effect is additional to the damping or growth occurring equally at all scales that can be predicted by the WKB approximation. In light of the recent observations by pulsar timing array (PTA) collaborations, and the potential observations by future detectors such as SKA, LISA, DECIGO, BBO, or ET, we show that, in most of the cases, constraints cannot be placed on the effects of α M and the initial GW energy density ℰ* GW separately, but only on the combined effects of the two, unless the signal is observed at different frequency ranges. In particular, we provide some constraints on the combined effects from the reported PTA observations.