Total number of plane partitions for $4$ or more dimensions According to MacMahon formula the total number $P_3(r, s, t)$ of plane partitions that fit in the $r \times s \times t$ box $\mathcal{B}(r,s,t)$ is equal to the following product formula:
$$
P_3(r,s,t)=\prod_{(i,j,k)\in \mathcal{B}(r,s,t)}\frac{i+j+k-1}{i+j+k-2}=\prod_{i=1}^{r}\prod_{j=1}^{s}\frac{i+j+t-1}{i+j-1}=\frac{H(r-1)H(s-1)H(t-1)H(r+s+t-1)}{H(r+s-1)H(r+t-1)H(s+t-1)}
$$
where $H(n) = 1! 2! \cdots n!$ is the superfactorial.
In case of two dimensions and the $r \times s$ rectangle $\mathcal{B}(r,s)$ the above reduces to:
$$P_2(r,s) = \binom{r+s}{r}=\frac{(r+s)!}{r!s!}$$
Is there a formula known for $4$ or higher dimensions, e.g. $P_4(r,s,t,u)$? Maybe using a "super-superfactorial" in a similar way?
 A: My comments above were a bit condensed, so let me spell things out in a little more detail here.
First of all, there is a question of what "dimension" one should consider a plane partition to be. When we say a plane partition "fits inside an $r\times s \times t$ box" we are treating it as a 3-dimensional object. But the word "plane" suggest two dimensions. And indeed, the original way MacMahon envisioned plane partitions was as 2-dimensional arrays of nonnegative integers, weakly decreasing along rows and down columns, which have only finitely many nonzero entries. This was supposed to be a two-dimensional generalization of (integer) partitions, which are 1-dimensional arrays of weakly decreasing nonnegative integers, eventually all zero. Hence the term "plane partition." To get back to the 3D picture we take the 2D array and view it as a floor plan for how to stack unit boxes shoved into a corner.
Thus, if you are interested in the next dimension up, the term to search for is "solid partitions": see, e.g., https://en.wikipedia.org/wiki/Solid_partition. There are some things known about these (including nontrivial results about their asymptotics), but they are much less tractable than plane partitions. In particular, no formulas like the beautiful ones you mention are known for these. In fact, as I suggested above, MacMahon guessed something wrong about the generating function of solid partitions.
To be more precise, we need to work with $q$-analogs of the quantities you are considering. Also, I am going to reindex so that plane partitions correspond to the more conventional dimension $d=2$. Namely, let
$$ P_2(q; r,s,t) := \sum_{\pi \subseteq r\times s \times t}q^{|\pi|}$$
be the generating function of all plane partitions $\pi$ in an $r\times s \times t$ box, according to sum of entries $|\pi|$ (when viewed as a 2d-array; or in the 3d picture, total number of unit boxes). MacMahon in fact proved a $q$-analog of the product formula you mentioned:
$$ P_2(q; r,s,t) = \prod_{i=1}^{r}\prod_{j=1}^{s}\prod_{k=1}^{t} \frac{[i+j+k-1]_q}{[i+j+k-2]_q} = \prod_{i=1}^{r}\prod_{j=1}^{s} \frac{[i+j+t-1]_q}{[i+j-1]_q},$$
where $[k]_q = (1-q^k)/(1-q) = 1+q+\cdots+q^{k-1}$ is the usual $q$-number. Note that the analogous $P_1(q;r,s)$ is then just the usual $q$-binomial coefficient. By taking the limit $r,s,t\to \infty$ we get the generating function of all plane partitions:
$$ P_2(q) := \sum_{\pi} q^{|\pi|} = \prod_{i=1}^{\infty} \frac{1}{(1-q^i)^i}.$$
Compare this to the generating function of all (1-dimensional) partitions:
$$ P_1(q) :=\sum_{\lambda} q^{|\lambda|} = \prod_{i=1}^{\infty} \frac{1}{1-q^i}.$$
Apparently MacMahon suggested that for $d$-dimensional partitions the corresponding generating function might satisfy
$$ P_d(q) = \prod_{i=1}^{\infty} (1-q^i)^{-\binom{i+d-2}{d-1}}$$
(see Stanley, Enumerative Combinatorics, Vol. 2, equation (7.122) on pg. 402). However, this fails already for $d=3$ at the coefficient of $q^6$ (see https://en.wikipedia.org/wiki/Solid_partition#Generating_function).
Obtaining the generating function for all $d$-dimensional partitions should be "easier" than counting $d$-dimensional partitions in a box, so this suggests the latter problem is also hopeless.
Then the question arises, what is special about the cases $d \leq 2$? Since the smaller dimensions are included in the case $d=2$, we may as well ask: what is special about two dimensions? There are a couple of possible answers. One is that, as 2D arrays of numbers, plane partitions are very close to Young diagrams, which govern the representation theory of the general linear group and hence enjoy nice formulas coming from algebra. Another answer is that plane partitions are related to exactly solvable models in statistical mechanics, like the dimer model, and these models exist only in two dimensions. Even more concretely, nonintersecting lattice paths can be efficiently enumerated by determinants using the Lindström-Gessel-Viennot lemma, but only in the special case of planar networks.

In the comments, you mention that the Dedekind numbers, counting the number of order ideals in the Boolean lattice, can also be thought of as a special case of higher-dimensional partition enumeration. This is true. But since to get the Dedekind numbers we need to be increasing the dimension $d$, the fact that Dedekind numbers are intractable does not in and of itself imply that for fixed $d$, enumeration of $d$-dimensional partitions must be intractable.
