There is an involution of $\tilde M_k^{\le l}$, $$i : \; f(\tau) \mapsto \overline{f(-\overline{\tau})}, \quad \sum_{n=0}^{\infty} a_n q^n \mapsto \sum_{n=0}^{\infty} \overline{a_n} q^n$$ which commutes with $D$, whose fixed points are exactly the quasimodular forms with real Fourier coefficients. So if $f \in \tilde M_{k, \mathbb{R}}^{\ell}$ is decomposed $$f = \sum_{j=0}^{\ell - 2} D^j f_j + c \cdot D^{\ell - 1} E_2, \quad f_j \in M_{k - 2j}, \; c \in \mathbb{C}$$ then $$f = i(f) = \sum_{j=0}^{\ell - 2} D^j i(f_j) + \overline{c} \cdot D^{\ell - 1} E_2.$$
Since the sum $\tilde M_k^{\le \ell} = \bigoplus_{j = 0}^{\ell - 2} D^j M_{k - 2j} \oplus \mathbb{C} D^{\ell - 1} E_2$ is direct, we have $$f_j = i(f_j) \; \text{for all} \; j \; \text{and} \; c = \overline{c},$$ so $f_j \in M_{k - 2j, \mathbb{R}}$ and $c \in \mathbb{R}$.