  We use the BayeSN hierarchical probabilistic SED model to analyse the
optical-NIR ($BVriYJH$) light curves of 86 Type Ia supernovae (SNe Ia) from the
Carnegie Supernova Project to investigate the SN Ia host galaxy dust law
distribution and correlations between SN Ia Hubble residuals and host mass. Our
Bayesian analysis simultaneously constrains the mass step and dust $R_V$
population distribution by leveraging optical-NIR colour information. We
demonstrate how a simplistic analysis where individual $R_V$ values are first
estimated for each SN separately, and then the sample variance of these point
estimates is computed, overestimates the $R_V$ population variance
$\sigma_R^2$. This bias is exacerbated when neglecting residual intrinsic
colour variation beyond that due to light curve shape. Instead, Bayesian
shrinkage estimates of $\sigma_R$ are more accurate, with fully hierarchical
analysis of the light curves being ideal. For the 75 SNe with low-to-moderate
reddening (peak apparent $B-V\leq0.3$), we estimate an $R_V$ distribution with
population mean $\mu_R=2.59\pm0.14$, and standard deviation
$\sigma_R=0.62\pm0.16$. Splitting this subsample at the median host galaxy mass
($10^{10.57}~\mathrm{M}_\odot$) yields consistent estimated $R_V$ distributions
between low- and high-mass galaxies, with $\mu_R=2.79\pm0.18$,
$\sigma_R=0.42\pm0.24$, and $\mu_R=2.35\pm0.27$, $\sigma_R=0.74\pm0.36$,
respectively. When estimating distances from the full optical-NIR light curves
while marginalising over various forms of the dust $R_V$ distribution, a mass
step of $\gtrsim0.06$ mag persists in the Hubble residuals at the median host
mass.
