Let $\phi$ be an $N$-function, (i.e. $\phi : \mathbb{R}_{\geq 0} \to \mathbb{R}_{\geq 0}$ is convex and satisfies $\lim_{t \to 0} \frac{\phi(t)}{t} = 0, \lim_{t\to \infty} \frac{\phi(t)}{t} = \infty$).
We can define the associated Luxemburg norm on the appropriate subspace of $\mathbb{R}$-valued random variables (Orlicz space $L_\phi$), by $\|Z\|_{\phi} := \inf \{ \lambda : \mathbb{E} \phi(|Z|/\lambda) \leq 1\}$. The space $L_\phi$ is just a space of all random variables for which this infimum is finite.
For example, when $\phi(t) = t^p$ for $p > 1$, this is just the $L_p$ norm, and when $\phi(t) = e^{t^2} -1$ the corresponding Orlicz space is the space of sub-gaussian random variables.
Question
Is the following statement true: Consider a sequence $Z_1, Z_2, \ldots Z_n$ of independent random variables with $\mathbb{E} Z_i=0$ and $\|Z_i\|_\phi \leq 1$ for all $i$. Then $$ \|\frac{1}{n} \sum_{i \leq n} Z_i\|_\phi \leq \lambda_\phi(n)$$ for some $\lambda_\phi$ s.t. $\lambda_\phi(n) \to 0$ as $n\to \infty$.
A bit more context
I am mostly interested in the scenario in which $\phi$ grows only slightly faster than linear, say $\phi(t) = t \ln(1+t)$.
Together with Markov inequality, this would imply that for any $\varepsilon, \delta$, there is some $n_\phi(\varepsilon, \delta)$ , s.t. $\mathbf{Pr}(|\sum_{i\leq n} \frac{1}{n} X_i| > \varepsilon) \leq \delta$. Law of large numbers asserts that as long as the first moment is bounded, the sample average $\sum \frac{1}{s} X_i$ converges to $0$ - and I hope to be able to quantify the speed of convergence under just slightly stronger assumption.
Note that this is true in $L_p$ spaces: for any $1< p \leq 2$, we have $\|\frac{1}{n}\sum_i Z_i\|_p \lesssim n^{1/p - 1}$, so the interesting case is $\phi$ growing slower than $t^{1+\gamma}$ for any $\gamma$.