  Bayesian graphical models are an efficient tool for modelling complex data
and derive self-consistent expressions of the posterior distribution of model
parameters. We apply Bayesian graphs to perform statistical analyses of Type Ia
supernova (SN Ia) luminosity distance measurements from the joint light-curve
analysis (JLA) data set. In contrast to the $\chi^2$ approach used in previous
studies, the Bayesian inference allows us to fully account for the
standard-candle parameter dependence of the data covariance matrix. Comparing
with $\chi^2$ analysis results, we find a systematic offset of the marginal
model parameter bounds. We demonstrate that the bias is statistically
significant in the case of the SN Ia standardization parameters with a maximal
6 $\sigma$ shift of the SN light-curve colour correction. In addition, we find
that the evidence for a host galaxy correction is now only 2.4 $\sigma$.
Systematic offsets on the cosmological parameters remain small, but may
increase by combining constraints from complementary cosmological probes. The
bias of the $\chi^2$ analysis is due to neglecting the parameter-dependent
log-determinant of the data covariance, which gives more statistical weight to
larger values of the standardization parameters. We find a similar effect on
compressed distance modulus data. To this end, we implement a fully consistent
compression method of the JLA data set that uses a Gaussian approximation of
the posterior distribution for fast generation of compressed data. Overall, the
results of our analysis emphasize the need for a fully consistent Bayesian
statistical approach in the analysis of future large SN Ia data sets.
