This paper is about the homogenization of linear elliptic operators in divergence form with stationary random coefficients that have only slowly decaying correlations. It deduces optimal estimates of the homogenization error from optimal growth estimates of the (extended) corrector. In line with the heuristics, there are transitions at dimension $d=2$, and for a correlation-decay exponent $β=2$; we capture the correct power of logarithms coming from these two sources of criticality. The decay of correlations is sharply encoded in terms of a multiscale logarithmic Sobolev inequality (LSI) for the ensemble under consideration --- the results would fail if correlation decay were encoded in terms of an $α$-mixing condition. Among other ensembles popular in modelling of random media, this class includes coefficient fields that are local transformations of stationary Gaussian fields. The optimal growth of the corrector $φ$ is derived from bounding the size of spatial averages $F=\int g\cdot\nablaφ$ of its gradient. This in turn is done by a (deterministic) sensitivity estimate of $F$, that is, by estimating the functional derivative $\frac{\partial F}{\partial a}$ of $F$ w.~r.~t.~the coefficient field $a$. Appealing to the LSI in form of concentration of measure yields a stochastic estimate on $F$. The sensitivity argument relies on a large-scale Schauder theory for the heterogeneous elliptic operator $-\nabla\cdot a\nabla$. The treatment allows for non-symmetric $a$ and for systems like linear elasticity.