  We develop a higher regularity theory for general quasilinear elliptic
equations and systems in divergence form with random coefficients. The main
result is a large-scale $L^\infty$-type estimate for the gradient of a
solution. The estimate is proved with optimal stochastic integrability under a
one-parameter family of mixing assumptions, allowing for very weak mixing with
non-integrable correlations to very strong mixing (e.g., finite range of
dependence). We also prove a quenched $L^2$ estimate for the error in
homogenization of Dirichlet problems. The approach is based on subadditive
arguments which rely on a variational formulation of general quasilinear
divergence-form equations.
