Let n ≥ 2 and Ω be a bounded Lipschitz domain in R n . In this article, the authors investigate global (weighted) estimates for the gradient of solutions to Robin boundary value problems of second order elliptic equations of divergence form with real-valued, bounded, measurable coefficients in Ω. More precisely, let p ∈ (n/(n − 1), ∞). Using a real-variable argument, the authors obtain two necessary and sufficient conditions for W 1,p estimates of solutions to Robin boundary value problems, respectively, in terms of a weak reverse Hölder inequality with exponent p or weighted W 1,q estimates of solutions with q ∈ (n/(n − 1), p] and some Muckenhoupt weights. As applications, the authors establish some global regularity estimates for solutions to Robin boundary value problems of second order elliptic equations of divergence form with small BMO coefficients, respectively, on bounded Lipschitz domains, C 1 domains or (semi-)convex domains, in the scale of weighted Lebesgue spaces, via some quite subtle approach which is different from the existing ones and, even when n = 3 in case of bounded C 1 domains, also gives an alternative correct proof of some know result. By this and some technique from harmonic analysis, the authors further obtain the global regularity estimates, respectively, in Morrey spaces, (Musielak-)Orlicz spaces and variable Lebesgue spaces.