Type Ia supernovae (SN Ia) are the most precise cosmological distance indicators, important for measuring cosmic acceleration and the properties of dark energy. Current and upcoming automated wide-field surveys will find large numbers of SN Ia, but cosmological inferences are already limited by systematic errors. Current cosmological analyses use optical light curves to estimate distances, whose accuracy is limited by the confounding effects of host galaxy dust extinction. The combination of optical and near infrared (NIR) light curves and spectroscopic data has the potential to improve inference and distance predictions in supernova cosmology. I have constructed a principled, hierarchical Bayesian framework, described by a directed acyclic graph, to coherently model the multiple random and uncertain effects underlying the observed SN Ia data, including measurement error, intrinsic supernova covariances, host galaxy dust extinction and reddening, and distances. An MCMC code, BayeSN, efficiently computes probabilistic inferences for the parameters of individual SN and the hyperparameters of the population. Application to optical, NIR, and spectroscopic data demonstrates that the combination of optical and NIR data approximately doubles the precision of cross-validated SN Ia distance predictions compared to using optical data alone, and estimates correlations between the intrinsic colors and characteristics of supernova spectral lines.