We analyze the mean square error of stochastic Runge--Kutta integrators for overdamped Langevin dynamics whose potential is strongly convex outside a bounded region. A decomposition splits the local error into a mean-zero term and a smaller remainder, and the discrete Poisson equation turns their moments into a bound on the error of a time average. We carry this out for two stochastic Runge--Kutta integrators proposed by Yang & Wang (2026), both of strong order $frac32$ and weak order $2$, which evaluate the gradient of the potential and no higher derivative. We show that their laws approach the law of the exact solution at second order in Wasserstein-1 distance, up to a logarithm, uniformly in the number of steps. For a test function with bounded derivatives up to third order, we prove that the mean square error over $N$ steps with step size $h$ is $mathcal O ( frac{1}{Nh} + h^4 )$, which is the optimal order in the discretization. Experiments measure the strong and weak orders and the sampling bias on a nonconvex potential in $mathbb R^2$, and compare the integrators on a diffusion model of CIFAR-10.