Binomial coefficients and derivatives of modular forms 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_8-2E_6E_2+E_4E_2^2)$$
$$E_4^{(3)}=\frac{5}{72}(-E_{10}+3E_8E_2-3E_6E_2^2+E_4E_2^3) $$
$$E_4^{(4)}=\frac{35}{864}(E_4^3-4E_{10}E_2+6E_8E_2^2-4E_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_2-3E_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+2n-2k}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
non-cusp forms.
 A: 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 Halperin-Fricke 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.
A: 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 Rankin--Cohen brakets; see Section 5.2 in [D. Zagier, Elliptic modular forms and their applications, The 1-2-3 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)_{n-r}}{(2\pi i)^{n-r}}
c^{n-r}(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)^{n-r}\binom{n}{r}E_{k+2n-2r}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\mapsto-1/\tau$
(that is, $a=d=0$, $b=-1$, and $c=1$). Then (we take $s=n-r$ in the above formula)
$$
D^nE_k(-1/\tau)
=\sum_{s=0}^n\binom ns\frac{(k+n-s) _ s}{(2\pi i)^s}\tau^{k+2n-s}D^{n-s}E_k(\tau)
$$
and
$$
\begin{aligned}
&
\frac{(k) _ n}{12^n}\sum_{r=0}^n(-1)^{n-r}\binom nrE_{k+2n-2r}E_2^r\bigg|_{\tau\mapsto-1/\tau}
\cr &\qquad
=\frac{(k) _ n}{12^n}\sum_{r=0}^n(-1)^{n-r}\binom nr\tau^{k+2n-r}E_{k+2n-2r}\biggl(\tau E_2+\frac{12}{2\pi i}\biggr)^r
\cr &\qquad
=\frac{(k) _ n}{12^n}\sum_{r=0}^n(-1)^{n-r}\binom nr\tau^{k+2n-r}E_{k+2n-2r}\sum_{s=0}^r\binom rs\frac{12^s}{(2\pi i)^s}\tau^{r-s}E_2^{r-s}
\cr &\qquad
=\frac{(k) _ n}{12^n}\sum_{s=0}^n\binom ns\tau^{r-s}\frac{12^s}{(2\pi i)^s}\tau^{k+2n-s}
\sum_{r=s}^n(-1)^{n-r}\binom{n-s}{r-s}E_{k+2n-2r}E_2^{r-s}.
\end{aligned}
$$
Subtracting the latter from the former we obtain
$$
\begin{aligned}
g_n(-1/\tau)
&=\sum_{s=0}^n\binom ns\frac{(k+n-s) _ s}{(2\pi i)^s}\tau^{k+2n-s}g_{n-s}(\tau)
\cr
&=\tau^{k+2n}g_n(\tau)+\sum_{s=1}^n\binom ns\frac{(k+n-s) _ s}{(2\pi i)^s}\tau^{k+2n-s}g_{n-s}(\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_{n-s}=0$ for $s=1,\dots,n$. 
The latter however does not happen when $k+2n>12$.
