  Here we investigate the $q$-series \begin{align*}
\mathcal{U}_a(q)&=\sum_{n=0}^{\infty} MO(a;n)q^n&:=\sum_{0< k_1<k_2<\cdots<k_a}
\frac{q^{k_1+k_2+\cdots+k_a}}{(1-q^{k_1})^2(1-q^{k_2})^2\cdots(1-q^{k_a})^2},\\
\mathcal{U}_a^{\star}(q)&=\sum_{n=0}^{\infty}M(a;n)q^n&:=\sum_{1\leq k_1\leq
k_2\leq\cdots\leq k_a}
\frac{q^{k_1+k_2+\cdots+k_a}}{(1-q^{k_1})^2(1-q^{k_2})^2\cdots(1-q^{k_a})^2}.
\end{align*} MacMahon introduced the $\mathcal{U}_a(q)$ in his seminal work on
partitions and divisor functions. Recent works show that these series are sums
of quasimodular forms with weights $\leq 2a.$ We make this explicit by
describing them in terms of Eisenstein series. We use these formulas to obtain
explicit and general congruences for the coefficients $MO(a;n)$ and $M(a;n).$
Notably, we prove the conjecture of Amdeberhan-Andrews-Tauraso as the $m=0$
special case of the infinite family of congruences $$ MO(11m+10; 11n+7)\equiv
0\pmod{11}, $$ and we prove that $$ MO(17m+16; 17n+15)\equiv 0\pmod{17}. $$ We
obtain further formulae using the limiting behavior of these series. For $n\leq
a+\binom{a+1}2,$ we obtain a ``hook length'' formulae for $MO(a;n)$, and for
$n\leq 2a$, we find that $M(a;n)=\binom{a+n-1}{n-a}+\binom{a+n-2}{n-a-1}.$
