By rotational symmetry we may assume that $Du(x_0)$ points in the direction of the first basis vector $e_1$. We must prove that
$${\int\!\!\!\!\!\!-}_{\partial B_R} [u-u(x_0)]\nu_1\le \sup u-u(x_0)$$
where $\nu_1$ is the first component of the unit normal vector. By the mean value property $u-u(x_0)$ has zero mean on $\partial B_R$. Thus, we can add any number to $\nu_1$ without changing the integral. Let's add $1$:
$${\int\!\!\!\!\!\!-}_{\partial B_R} [u-u(x_0)]\nu_1={\int\!\!\!\!\!\!-}_{\partial B_R} [u-u(x_0)](\nu_1+1)$$
Now that the factor $\nu_1+1$ is nonnegative, we use a one-sided bound on $u-u(x_0)$:
$${\int\!\!\!\!\!\!-}_{\partial B_R} [u-u(x_0)](\nu_1+1)\le (\sup u-u(x_0)){\int\!\!\!\!\!\!-}_{\partial B_R}(\nu_1+1)$$
Finally, $${\int\!\!\!\!\!\!-}_{\partial B_R}(\nu_1+1)=1$$ because $\nu_1$ has zero mean.
This
$$|D_iu(y)| \leq \frac nR \sup_{\partial\Omega}|u| \leq \frac n {d(y,\partial\Omega)} \sup_\Omega |u|\tag{$\ast$}$$
is wrong for $R < d(y,\partial\Omega)$ and harmonic $u \not\equiv 0$. The first of the two inequalities is okay if $u$ has decent boundary values on $\partial\Omega$ or we replace $\partial\Omega$ with $\partial B_R(y)$ or $\Omega$ there, but not the second. Since $u$ is harmonic, if it has nice (for example continuous) boundary values, we have
$$\sup_{\Omega} \lvert u\rvert = \sup_{\partial\Omega} \lvert u\rvert,$$
and since $R < d(y,\partial\Omega)$, the second inequality can only hold if the factor $\sup\limits_\Omega \lvert u\rvert$ is $0$. If we don't have good boundary values, the $\sup\limits_{\partial\Omega}$ would anyway have to be replaced with $\sup\limits_\Omega$ or $\sup\limits_{\partial B_R(y)}$. In either case, there are harmonic $u$ (constants, for example) for which the second inequality in $(\ast)$ doesn't hold.
However, the first inequality of $(\ast)$ [with the supremum taken over $\Omega$ to avoid demanding boundary regularity] holds for all $R < d(y,\partial\Omega)$, and hence
$$\lvert D_i u(y)\rvert \leqslant \inf_{R < d(y,\partial\Omega)} \frac{n}{R}\sup_{\Omega} \lvert u\rvert = \frac{n}{d(y,\partial\Omega)} \sup_{\Omega} \lvert u\rvert$$
gives us the estimate we are interested in.
Best Answer
This follows from the scalar version of the estimate:
noting that, since $u$ is harmonic, $D_iu$ is harmonic too.
The proof relies on Caccioppolli inequality and the estimate
Proof of Theorem: If $4r\geq s$, then we have simply $$ \int_{B_r(0)}[f-f(0)]^2\leq\left(\frac{4r}s\right)^{n+2}\int_{B_r(0)}[f-f(0)]^2\leq 4^{n+2}\left(\frac rs\right)^{n+2}\int_{B_s(0)}[f-f(0)]^2. $$ Otherwise, we have $r\leq s\leq 4r$ and we need to do a little more work. Separately bounding each component of $\nabla f$ using Lemma, we have $$ \sup_{B_r(0)}\lvert\nabla f\rvert^2\leq2^n-\!\!\!\!\!\!\!\int_{B_{2r}(0)} \lvert\nabla f\rvert^2 $$ and so $$ \begin{align*} -\!\!\!\!\!\!\int_{B_r(0)}[f-f(0)]^2 &\leq -\!\!\!\!\!\!\int_{B_r(0)}\left(r\sup_{B_r(0)}\lvert\nabla f\rvert\right)^2\\ &=r^2\sup_{B_r(0)}\lvert\nabla f\rvert^2\\ &\leq r^2\sup_{B_{s/4}(0)}\lvert\nabla f\rvert^2\\ &\leq r^2 2^n-\!\!\!\!\!\!\!\int_{B_{s/2}(0)}\lvert\nabla f\rvert^2\\ &\leq r^2 2^n \frac{2^n}{s^2}-\!\!\!\!\!\!\!\int_{B_s(0)}[f-f(0)]^2 \end{align*} $$ where the last inequality is by Caccioppolli. Rearranging, we obtain $$ \int_{B_r(0)}[f-f(0)]^2\leq 4^n\left(\frac rs\right)^{n+2}\int_{B_s(0)}[f-f(0)]^2 $$