  The direct detection of gravitational waves by the LIGO-Virgo collaboration
has opened a new window with which to measure cosmological parameters such as
the Hubble constant $H_0$, and also probe general relativity on large scales.
In this paper we present a new phenomenological approach, together with its
inferencial implementation, for measuring deviations from general relativity
(GR) on cosmological scales concurrently with a determination of $H_0$. We
consider gravitational waves (GWs) propagating in an expanding homogeneous and
isotropic background, but with a modified friction term and dispersion relation
relative to that of GR. We find that a single binary neutron star GW detection
will poorly constrain the GW friction term. However, a joint analysis including
the GW phase and GW-GRB detection delay could improve constraints on some GW
dispersion relations provided the delay is measured with millisecond precision.
We also show that, for massive gravity, by combining 100 binary neutron stars
detections with observed electromagnetic counterparts and hosting galaxy
identification, we will be able to constrain the Hubble constant, the GW
damping term and the GW dispersion relation with 2\%, 15\% and 2 \% accuracy,
respectively. We emphasise that these three parameters should be measured
together in order avoid biases. Finally we apply the method to GW170817, and
demonstrate that for all the GW dispersions relations we consider, including
massive gravity, the GW must be emitted $\sim$ 1.74s before the Gamma-ray burst
(GRB). Furthermore, at the GW merger peak frequency, we show that the
fractional difference between the GW group velocity and $c$ is $\lesssim
10^{-17}$.
