Let $E_2$, $E_4$, and $E_6$ denote the standard Eisenstein series. The usual variables $q=e^{2\pi i\tau}$ allow us to regard the $E_n$'s as functions on either the upper half plane or the unit disk and we can define $E_n'=\frac{1}{2\pi i}\frac{d}{d\tau}E_n(\tau)=q\frac{d}{dq}E_n(q)$. I had cause to calculate a few of these and saw $$ E_4'=\frac{1}{3}(E_6+E_4E_2) $$ $$ E_4''=\frac{5}{36}(E_82E_6E_2+E_4E_2^2)$$ $$E_4^{(3)}=\frac{5}{72}(E_{10}+3E_8E_23E_6E_2^2+E_4E_2^3) $$ $$E_4^{(4)}=\frac{35}{864}(E_4^34E_{10}E_2+6E_8E_2^24E_6E_2^3+E_4E_2^4)40\Delta $$ and $$E_6'=\frac{1}{2}(E_8+E_6E_2) $$ $$E_6''=\frac{7}{24}(E_{10}2E_8E_2+E_6E_2^2) $$ $$E_6^{(3)}=\frac{7}{36}(E_4^3+3E_{10}E_23E_8E_2^2+E_6E_2^3)+168\Delta $$ It's a standard fact that the derivative of a modular form is quasimodular, so it's not surprising that we have polynomials in $E_2$. I am surprised about the appearance of the binomial coefficients though. Is there a deeper reason for their appearance? Also, I wonder if the/a pattern continues. For instance, it would be interesting if it happens that there always is some $\alpha \in \mathbb{Q}$ so that $$E_4^{(n)}\alpha \sum_{k=0}^{n} (1)^{k+n}\binom {n}{k}E_{4+2n2k}E_2^{k}$$ is modular (and similarly for $E_6$). The other direction you could ask if the pattern extends is for other modular forms besides $E_4$ and $E_6$. I've taken a handful of derivatives of other Eisenstein series and saw similar results. You don't get the binomial coefficients though when you take derivatives of $\Delta$, so maybe at most something general can be said is for noncusp forms.

The constant $\alpha$ in your question can be in fact written explicitly as $(k)_n/12^n$, where $(a)_n=\Gamma(a+n)/\Gamma(n)$ is the Pochhammer symbol (shifted factorial) and $k$ denotes the (even) weight of the corresponding Eisenstein series. Your observation is indeed related to the RankinCohen brakets; see Section 5.2 in [D. Zagier, Elliptic modular forms and their applications, The 123 of modular forms, Universitext (Springer, Berlin, 2008), pp. 1–103]. Preserving the notation $D$ of Zagier's lectures for your differential operator and picking a modular form $f$ of weight $k$, one can show that $D^nf$ transforms under the modular group as $$ D^nf\biggl(\frac{a\tau+b}{c\tau+d}\biggr) =\sum_{r=0}^n\binom{n}{r}\frac{(k+r)_{nr}}{(2\pi i)^{nr}} c^{nr}(c\tau+d)^{k+n+r}D^rf(\tau), $$ by the induction on $n\ge 0$. In addition, the function $E_2(\tau)$ transforms as $$ E_2\biggl(\frac{a\tau+b}{c\tau+d}\biggr) =\frac{12c(c\tau+d)}{2\pi i}+(c\tau+d)^2E_2(\tau). $$ Therefore, it remains to verify that the difference $$ g_n=D^nE_k\frac{(k) _ n}{12^n}\sum_{r=0}^n(1)^{nr}\binom{n}{r}E_{k+2n2r}E_2^r $$ satisfies $$ g_n\biggl(\frac{a\tau+b}{c\tau+d}\biggr)=(c\tau+d)^{k+2n}g_n(\tau). $$ Technicalities. Indeed, I found the remaining details quite boring, but going through my yesterday writing I have realised that your expectation fails already for $D^5E_4$ and $D^4E_6$. Here is my explanation why. Because $g_n(\tau)$ is a $q$series, so it is invariant under $\tau\mapsto\tau+1$, we can restrict to verifying the claim under the transformation $\tau\mapsto1/\tau$ (that is, $a=d=0$, $b=1$, and $c=1$). Then (we take $s=nr$ in the above formula) $$ D^nE_k(1/\tau) =\sum_{s=0}^n\binom ns\frac{(k+ns) _ s}{(2\pi i)^s}\tau^{k+2ns}D^{ns}E_k(\tau) $$ and $$ \begin{aligned} & \frac{(k) _ n}{12^n}\sum_{r=0}^n(1)^{nr}\binom nrE_{k+2n2r}E_2^r\bigg_{\tau\mapsto1/\tau} \cr &\qquad =\frac{(k) _ n}{12^n}\sum_{r=0}^n(1)^{nr}\binom nr\tau^{k+2nr}E_{k+2n2r}\biggl(\tau E_2+\frac{12}{2\pi i}\biggr)^r \cr &\qquad =\frac{(k) _ n}{12^n}\sum_{r=0}^n(1)^{nr}\binom nr\tau^{k+2nr}E_{k+2n2r}\sum_{s=0}^r\binom rs\frac{12^s}{(2\pi i)^s}\tau^{rs}E_2^{rs} \cr &\qquad =\frac{(k) _ n}{12^n}\sum_{s=0}^n\binom ns\tau^{rs}\frac{12^s}{(2\pi i)^s}\tau^{k+2ns} \sum_{r=s}^n(1)^{nr}\binom{ns}{rs}E_{k+2n2r}E_2^{rs}. \end{aligned} $$ Subtracting the latter from the former we obtain $$ \begin{aligned} g_n(1/\tau) &=\sum_{s=0}^n\binom ns\frac{(k+ns) _ s}{(2\pi i)^s}\tau^{k+2ns}g_{ns}(\tau) \cr &=\tau^{k+2n}g_n(\tau)+\sum_{s=1}^n\binom ns\frac{(k+ns) _ s}{(2\pi i)^s}\tau^{k+2ns}g_{ns}(\tau). \end{aligned} $$ Therefore, $g_n(1/\tau)=\tau^{k+2n}g_n(\tau)$, hence $g_n(\tau)$ is a modular form (of weight $k+2n$), if and only if the additional sum over $s$ vanishes, that is, $g_{ns}=0$ for $s=1,\dots,n$. The latter however does not happen when $k+2n>12$. 


This isn't really an answer, but a long comment with a bit of LaTeX that thought would render poorly in the comment box. The following fact may be lurking in the background here: While the derivative of a modular form of weight $k$ is not generally modular, the map $D$ on modular forms of weight $k$ defined by $$D(f) = q\frac{df}{dq}  \frac{k}{12}f\cdot E_2$$ actually does preserve modularity. I think this is sometimes called the HalperinFricke operator or something like this. It's also a derivation, for what it's worth. It certainly directly explains the first of your equations above, and I wonder if some cleverness iterating it would yield your more general observations. 

