We compare GW150914 directly to simulations of coalescing binary black holes in full general relativity, including several performed specifically to reproduce this event. Our calculations go beyond existing semianalytic models, because for all simulations—including sources with two independent, precessing spins—we perform comparisons which account for all the spin-weighted quadrupolar modes, and separately which account for all the quadrupolar and octopolar modes. Consistent with the posterior distributions reported by Abbott et al. [Phys. Rev. Lett. 116, 241102 (2016)] (at the 90% credible level), we find the data are compatible with a wide range of nonprecessing and precessing simulations. Follow-up simulations performed using previously estimated binary parameters most resemble the data, even when all quadrupolar and octopolar modes are included. Comparisons including only the quadrupolar modes constrain the total redshifted mass $M_z\in [64M_{\odot }-82M_{\odot }]$, mass ratio $1/q=m_2/m_1\in [0.6,1]$, and effective aligned spin ${\chi }_{\mathrm{eff}}\in [-0.3,0.2]$, where ${\chi }_{\mathrm{eff}}=(S_1/m_1+S_2/m_2)\cdot ^L/M$. Including both quadrupolar and octopolar modes, we find the mass ratio is even more tightly constrained. Even accounting for precession, simulations with extreme mass ratios and effective spins are highly inconsistent with the data, at any mass. Several nonprecessing and precessing simulations with similar mass ratio and ${\chi }_{\mathrm{eff}}$ are consistent with the data. Though correlated, the components’ spins (both in magnitude and directions) are not significantly constrained by the data: the data is consistent with simulations with component spin magnitudes $a_{1,2}$ up to at least 0.8, with random orientations. Further detailed follow-up calculations are needed to determine if the data contain a weak imprint from transverse (precessing) spins. For nonprecessing binaries, interpolating between simulations, we reconstruct a posterior distribution consistent with previous results. The final black hole’s redshifted mass is consistent with $M_{f,z}$ in the range $64.0M_{\odot }-73.5M_{\odot }$ and the final black hole’s dimensionless spin parameter is consistent with $a_f=0.62–0.73$. As our approach invokes no intermediate approximations to general relativity and can strongly reject binaries whose radiation is inconsistent with the data, our analysis provides a valuable complement to Abbott et al. [Phys. Rev. Lett. 116, 241102 (2016)].