We present a framework for self-consistent cosmological analyses of the full-shape anisotropic bispectrum, including the quadrupole ( = 2) and hexadecapole ( = 4) moments. This features a novel window-free algorithm for extracting the latter quantities from data, derived using a maximum-likelihood prescription. Furthermore, we introduce a theoretical model for the bispectrum multipoles (which does not introduce new free parameters), and test both aspects of the pipeline on several high-fidelity mocks, including the PT Challenge suite of gigantic cumulative volume. This establishes that the systematic error is significantly below the statistical threshold, both for the measurement and modeling. As a realistic example, we extract the large-scale bispectrum multipoles from BOSS DR12 and analyze them in combination with the power spectrum data. Assuming a minimal ΛCDM model, with a BBN prior on the baryon density and a Planck prior on n s , we can extract the remaining cosmological parameters directly from the clustering data. The inclusion of the unwindowed higher-order ( > 0) large-scale bispectrum multipoles is found to moderately improve one-dimensional cosmological parameter posteriors (at the 5% − 10% level), though these multipoles are detected only in three out of four BOSS data segments at ≈ 5σ. Combining information from the power spectrum and bispectrum multipoles, the real space power spectrum, and the post-reconstructed BAO data, we find H 0 = 68.2 ± 0.8 km s −1 Mpc −1 , Ω m = 0.33 ± 0.01 and σ 8 = 0.736 ± 0.033 (the tightest yet found in perturbative full-shape analyses). Our estimate of the growth parameter S 8 = 0.77 ± 0.04 agrees with both weak lensing and CMB results. The estimators and data used in this work have been made publicly available.