For simplicity assume $f$ is smooth. It turns out you can always approximate a subharmonic function by smooth subharmonic functions, so you aren't really losing any generality by doing so.
$(iv)\iff(iii)$: This is more or less by definition. Suppose $\phi\in C_c^\infty(U)$ is such that $\phi\geq 0$. Integration by parts says $$\int_Uf\Delta\phi = \int_U\phi\Delta f.$$ Suppose $(iv)$ holds, i.e., $\Delta f\geq 0$. Since $\phi\geq 0$, the right hand side of the equality is $\geq 0$, hence so is the left. Now assume $(iii)$. We know the left hand side is always $\geq 0$, and hence $\int \phi\Delta f\geq 0$ whenever $\phi\geq 0$ is compactly supported smooth. This means $\Delta f\geq 0$, as otherwise you could choose an appropriately constructed $\phi$ for which the right hand side would be negative.
You say you have shown $(i)\iff (iv)$, so this shows $(i)\iff (iv)\iff (iii)$.
The proof of $(ii)\implies (i)$ that I know uses a bit more -- the solution of the Dirichlet problem. Suppose $\mathbb{B}_r(x)\subset U$, and let $u$ be a harmonic function on $\mathbb{B}_r(x)$ that agrees with $f$ on the boundary of $\mathbb{B}_r(x)$. Because of $(ii)$, one has $f\leq u$ on all of $\mathbb{B}_r(x)$. In particular, $$f(x)\leq u(x) = \frac{1}{n\alpha(n)r^{n-1}}\int_{\partial\mathbb{B}_r(x)} u(y)\,dS(y) = \frac{1}{n\alpha(n)r^{n-1}}\int_{\partial\mathbb{B}_r(x)} f(y)\,dS(y).$$
It only remains to show that something implies $(ii)$. We'll do $(i)\implies (ii)$. Suppose that $(ii)$ is false. Let $D\Subset U$ be a domain and $u\colon D\to \mathbb{R}$ a harmonic function on $D$ which agrees with $f$ on $\partial D$, but for which there is a point $x\in D$ such that $f(x)>u(x)$. Let $V$ be the collection of points $\{x\in D : f(x)>u(x)$}. This is an open set which is nonempty by assumption. Fix an $x_0\in V$ for which $f(x_0) - u(x_0)$ is maximal. Then for small $r>0$, one has that $\mathbb{B}_r(x_0)\subset V$, and thus $$f(x_0) - u(x_0)\leq \frac{1}{n\alpha(n)r^{n-1}}\int_{\partial\mathbb{B}_r(x_0)} f(y) - u(y)\,dS(y).$$ On the other hand, the integrand $f(y)-u(y)$ is $\leq f(x_0) - u(x_0)$ by maximality, so the only way for this inequality to hold is if $f(y) - u(y)\equiv f(x_0) - u(x_0)$ for $y\in\partial\mathbb{B}_r(x_0)$. Since $r>0$ could have been anything (small), it follows that $f(x) - u(x)$ is constant in a neighborhood of $x_0$. By a connectedness argument, we see that $f(x) - u(x)$ is constant in $D$. However, $f$ and $u$ agree on $\partial D$, so it must be that $f\equiv u$ on $D$, contradicting our assumption that there is some $x$ with $f(x)>u(x)$. We conclude that $(i)\implies (ii)$.
In order to invoke the mean-value theorem, you need to assume that $B_{R}(x_{0}) \subset U$.
What you have shown is that for any $\mathbf{x}_{0} \in V$, the $B_{R}(x_{0})$ consists entirely of points where $u(\mathbf{x})=u(\mathbf{x}_{0})$. This means that $\mathbf{x} \in V$ for all $\mathbf{x} \in B_{R}(x_{0})$. Hence, $B_{R}(x_{0}) \subseteq V$.
You have hence shown that for any $\mathbf{x}_{0} \in V$ there exists an $R>0$ for which $B_{R}(x_{0}) \subseteq V$. By definition, this means that $V$ is an open set.
Best Answer
The proof is correct. The strong maximum principle "a nonconstant subharmonic function attains its global maximum only on the boundary" does hold for subharmonic functions.
I guess your source was referring to "a nonconstant harmonic function has no local maxima", which indeed does not generalize to subharmonic functions. E.g., $u(x) = \max(x_1, 0)$ is a nonconstant subharmonic function on the unit ball which has a local maximum at every point of the left half of the ball.