A plane partition is a subset of $\mathbb Z_{\geqslant0}^3$ s.t. if it contains $(i+1,j,k)$ or $(i,j+1,k)$ or $(i,j,k+1)$ it also contains $(i,j,k)$.

What is the generating function $R(q)$ of (volumes of) plane partitions not containing the cell $(1,1,1)$?

Since a plane partition containing $(1,1,1)$ contains a $2\times2\times2$ cube, mod $q^8$ this g.f. coincides with MacMahon’s $$ M(q)=\frac1{\prod(1-q^i)^i}=1+q+3q^2+6q^3+13q^4+24q^5+48q^6+86q^7+160q^8+282q^9+500q^{10}+\ldots $$ and computing a couple more terms $$ R(q)=1+q+3q^2+6q^3+13q^4+24q^5+48q^6+86q^7+159q^8+279q^9+488q^{10}+\ldots $$ OEIS doesn’t know this sequence, but in the form $$ R(q)=1+\frac{q-3q^3+6q^6-10q^{10}+\ldots}{[(1-q)(1-q^2)(1-q^3)\ldots]^3}. $$ a pattern is evident.

So actually the question is

How to prove that $R(q)=1+\frac1{(q)_\infty^3}\sum\limits_i(-1)^i\frac{i(i-1)}2q^{\tfrac{i(i-1)}2}$?

AFAIK the statement can be found somewhere in the literature on representation theory of quantum toroidal something. But surely there is a well-known combinatorial (in the broad sense, not necessarily bijective) proof?

Context 1. G.f. for plane partitions not containing $(1,1,0)$ is $1/(q)_\infty^2\cdot\sum(-1)^iq^{\frac{i(i+1)}2}$. See the subsection about «V-shaped partitions» in Stanley's Enumerative Combinatorics for a (simple) combinatorial proof. Perhaps there is something like that proof in «(1,1,1)-case».

Context 2. G.f. for plane partitions not containing $(1,1,1)$ but intersecting axes at $(i,0,0)$, $(0,j,0)$ and $(0,0,k)$ resp. is $\genfrac[]{0pt}{}{i+j}i\genfrac[]{0pt}{}{j+k}j\genfrac[]{0pt}{}{k+i}kq^{i+j+k+1}$. So maybe it’s possible to compute $R(q)$ as $1+\sum_n q^{n+1}\sum_{i+j+k=n}\genfrac[]{0pt}{}{i+j}i\genfrac[]{0pt}{}{j+k}j\genfrac[]{0pt}{}{k+j}k$. And for the internal sum I know the answer at least at $q=1$: https://math.stackexchange.com/q/177209/