We analyze a set of 89 Type Ia supernovae (SN Ia) that have both optical and near-infrared (NIR) photometry to derive distances and construct low redshift (z ≤ 0.04) Hubble diagrams. We construct mean light curve (LC) templates using a hierarchical Bayesian model. We explore both Gaussian process (GP) and template methods for fitting the LCs and estimating distances, while including peculiar velocity and photometric uncertainties. For the 56 SN Ia with both optical and NIR observations near maximum light, the GP method yields a NIR-only Hubble-diagram with a RMS of 0.117 ± 0.014 mag when referenced to the NIR maxima. For each NIR band, a comparable GP method RMS is obtained when referencing to NIR-max or B-max. Using NIR LC templates referenced to B-max yields a larger RMS value of 0.138 ± 0.014 mag. Fitting the corresponding optical data using standard LC fitters that use LC shape and color corrections yields larger RMS values of 0.179 ± 0.018 mag with SALT2 and 0.174 ± 0.021 mag with SNooPy. Applying our GP method to subsets of SN Ia NIR LCs at NIR maximum light, even without corrections for LC shape, color, or host-galaxy dust reddening, provides smaller RMS in the inferred distances, at the ∼ 2.3-4.1σ level, than standard optical methods that do correct for those effects. Our ongoing RAISIN program on the Hubble Space Telescope will exploit this promising infrared approach to limit systematic errors when measuring the expansion history of the universe to constrain dark energy.