The problem of anomalous scaling in magnetohydrodynamics turbulence is considered within the framework of the kinematic approximation, in the presence of a large-scale background magnetic field. Field theoretic renormalization group methods are applied to the Kazantsev-Kraichnan model of a passive vector advected by the Gaussian velocity field with zero mean and correlation function $\propto \delta(t-t')/k^{d+\epsilon}$. Inertial-range anomalous scaling for the tensor pair correlators is established as a consequence of the existence in the corresponding operator product expansions of certain "dangerous" composite operators, whose negative critical dimensions determine the anomalous exponents. The main technical result is the calculation of the anomalous exponents in the order $\epsilon^2$ of the $\epsilon$ expansion (two-loop approximation).