We present results obtained by applying our BAyesian HierArchical Modeling for the Analysis of Supernova cosmology (BAHAMAS) software package to the 740 spectroscopically confirmed supernovae type Ia (SNIa) from the "Joint Light-curve Analysis" (JLA) dataset. We simultaneously determine cosmological parameters and standardization parameters, including host galaxy mass corrections, residual scatter and object-by-object intrinsic magnitudes. Combining JLA and Planck Cosmic Microwave Background data, we find significant discrepancies in cosmological parameter constraints with respect to the standard analysis: we find Omega_M = 0.399+/-0.027, 2.8\sigma\ higher than previously reported and w = -0.910+/-0.045, 1.6\sigma\ higher than the standard analysis. We determine the residual scatter to be sigma_res = 0.104+/-0.005. We confirm (at the 95% probability level) the existence of two sub-populations segregated by host galaxy mass, separated at log_{10}(M/M_solar) = 10, differing in mean intrinsic magnitude by 0.055+/-0.022 mag, lower than previously reported. Cosmological parameter constraints are however unaffected by inclusion of host galaxy mass corrections. We find ~4\sigma\ evidence for a sharp drop in the value of the color correction parameter, beta(z), at a redshift z_trans = 0.662+/-0.055. We rule out some possible explanations for this behaviour, which remains unexplained.