We present a calculation of the hyperfine splittings in bottomonium using lattice Nonrelativistic QCD. The calculation includes spin-dependent relativistic corrections through O(v^6), radiative corrections to the leading spin-magnetic coupling and, for the first time, non-perturbative 4-quark interactions which enter at alpha_s^2 v^3. We also include the effect of u,d,s and c quark vacuum polarisation. Our result for the 1S hyperfine splitting is M(Upsilon,1S) - M(eta_b,1S)= 60.0(6.4) MeV. We find the ratio of 2S to 1S hyperfine splittings (M(Upsilon,2S) - M(eta_b,2S))/ (M(Upsilon,1S) - M(eta_b,1S)) = 0.445(28).