  A joint hierarchical Bayesian analysis of the binary black hole (BBH) mass
function, merger rate evolution and cosmological parameters can be used to
extract information on both the cosmological and population parameters. We
extend this technique to include the effect of modified gravitational wave (GW)
propagation. We discuss the constraints on the parameter $\Xi_0$ that describes
this phenomenon (with $\Xi_0=1$ in General Relativity, GR) using the data from
the GWTC-3 catalog. We find the constraints $\Xi_0 = 1.2^{+0.7}_{-0.7}$ with a
flat prior on $\Xi_0$, and $\Xi_0 = 1.0^{+0.4}_{-0.8}$ with a prior uniform in
$\log\Xi_0$ ($68\%$ C.L., maximum posterior and HDI), which only rely on the
presence of a feature in the BBH mass distribution around $\sim 30-45
M_{\odot}$, and are robust to whether or not the event GW190521 is considered
an outlier of the population. We then study in more detail the effects of
modified GW propagation on population and cosmological analyses for LIGO/Virgo
at design sensitivity. For a given data-taking period, the relative error
$\Delta\Xi_0/\Xi_0$ has a significant dependence on the fiducial value of
$\Xi_0$, since the latter has a strong influence on the detection rate. For
five years of data, the accuracy ranges from $\sim 10\%$ on $\Xi_0$ when
$\Xi_0=1$ to $\Delta\Xi_0/\Xi_0\sim 20\%$ for $\Xi_0=1.8$ - a large deviation
from GR, still consistent with current limits and predicted by viable
cosmological models. For the Hubble parameter, we forecast an accuracy of
$\Delta H_0/H_0 \sim 20\%$, and an accuracy on $H(z)$ of $\sim7\%$ at a pivot
redshift $z_*\sim 0.8$. We finally show that, if Nature is described by a
modified gravity theory with a large deviation from the GR value $\Xi_0=1$,
such as $\Xi_0=1.8$, analysing the data assuming GR produces a significant bias
in the inferred values of the mass scales, Hubble constant, and BBH merger
rate.
