Only numerical relativity simulations can capture the full complexities of binary black hole mergers . These simulations , however , are prohibitively expensive for direct data analysis applications such as parameter estimation . We present two new fast and accurate surrogate models for the outputs of these simulations : the first model , NRSur7dq4 , predicts the gravitational waveform and the second model , NRSur7dq4Remnant , predicts the properties of the remnant black hole . These models extend previous 7-dimensional , non-eccentric precessing models to higher mass ratios , and have been trained against 1528 simulations with mass ratios q \leq 4 and spin magnitudes \chi _ { 1 } , \chi _ { 2 } \leq 0.8 , with generic spin directions . The waveform model , NRSur7dq4 , which begins about 20 orbits before merger , includes all \ell \leq 4 spin-weighted spherical harmonic modes , as well as the precession frame dynamics and spin evolution of the black holes . The final black hole model , NRSur7dq4Remnant , models the mass , spin , and recoil kick velocity of the remnant black hole . In their training parameter range , both models are shown to be more accurate than existing models by at least an order of magnitude , with errors comparable to the estimated errors in the numerical relativity simulations . We also show that the surrogate models work well even when extrapolated outside their training parameter space range , up to mass ratios q = 6 .