Numerical methods based on the Magnus expansion are an efficient class of integrators for Schrödinger equations with time-dependent Hamiltonian. Though their derivation assumes an unreasonably small time step size as would be required for a standard explicit integrator, the methods perform well even for much larger step sizes. This favorable behavior is explained, and optimal-order error bounds are derived which require no or only mild restrictions of the step size. In contrast to standard integrators, the error does not depend on higher time derivatives of the solution, which is in general highly oscillatory. Even more interesting than the error bounds themselves are the mechanisms which lead to these bounds and which make Magnus methods perform so well for Schrödinger equations. Therefore, we will not only present the results but also outline a general procedure for obtaining such error bounds.