As Ah Wong mentioned above, the justification is the Inclusion-Exclusion Principle. Without loss of generality, consider a system of $N$ subsystems such that subsystem $i$ is in the state $r_i$. The state of the full system can be regarded as the set $R = \{r_i\}_{i=1}^N$. While the specific mathematical description of the states $r_i$ is irrelevant to the considerations below, it is helpful for grounding one's intuition to consider a specific example. In the case of quantum-mechanical ground-state energy of a system of many molecules as a function of the positions and atomic numbers of the atoms within each molecule (in the Born-Oppenheimer approximation), the state of molecule $i$ is an element of the Cartesian product $\mathbb{R}^{3n_i} \times \mathbb{N}^{n_i}$, where $n_i$ is the number of atoms in molecule $i$. The component of the state from $\mathbb{R}^{3n_i}$ describes the positions of the atoms in the molecule, while the component from $\mathbb{N}^{n_i}$ supplies their atomic numbers.
The energy can be regarded as a well-defined real-valued function of any subset $S \subseteq R$, and so can be denoted as $E(S)$. There are $2^N$ such subsets, as each subset can either include or exclude each member of $R$. The subset that includes all members is, of course, just $R$ itself, and the subset that includes none is the empty set. These $2^N$ subsets comprise the elements of the power set of $R$, which is often denoted as $\mathcal{P}(R)$. Formally, the energy function can be described as the map $E: \mathcal{P}(R) \to \mathbb{R}$.
The ansatz of the many-body expansion is that there exists a second function $V: \mathcal{P}(R) \to \mathbb{R}$ defined such that for each $S \subseteq \mathcal{P}(R)$, of which there are $2^N$, the following relation holds: $$E(S) = \sum_{T \in \mathcal{P}(S)} V(T)$$
Note that the sum is over elements of $\mathcal{P}(S)$ (i.e. subsets of $S$), not over subsets of $R$. The simplest case is when $S = \emptyset$: $$E(\emptyset) = \sum_{T \in \mathcal{P}(\emptyset)} V(T) = V(\emptyset)$$
In the quantum-mechanical case, the energy of zero particles is zero, and so the zeroth-order many-body term must also be zero.
For the general sum, while any order of summation is valid, it is perhaps most intuitive to consider the sum as ascending through a partial ordering of $\mathcal{P}(R)$, which can be regarded as an indexed sequence of the elements $\{S_i \in \mathcal{P}(R)\}_{i = 1}^{2^N}$ such that: $$S_1 \subseteq S_2 \subseteq \cdots \subseteq S_{2^N - 1} \subseteq S_{2^N}$$
Any power set is guaranteed to have a partial ordering, and there is more than one partial ordering for the power set of any set with more than one element. For the purposes of this derivation, any choice of partial ordering will suffice. Thus, the substitutions $T \to S_i$ and $S \to S_j$ give: $$E(S_j) = \sum_{S_i \in \mathcal{P}(S_j)} V(S_i)$$
As the sum over $V(S_i)$ above is defined over $\mathcal{P}(S_j)$ instead of the full $\mathcal{P}(R)$, it will be helpful as a bookkeeping device to introduce the following indicator function: $$h(A, B) =\begin{cases}1 & \textrm{if} \ A \subseteq B \\ 0 & \textrm{if} \ B \subset A \\ \end{cases}$$
By the definition of power sets, if $S_i \in \mathcal{P}(S_j)$, it follows that $S_i \subseteq S_j$ and thus $h(S_i, S_j) = 1$. This enables the $2^N$ sums to all be written over the full set of indices as follows: $$E(S_j) = \sum_{i = 1}^{2^N} V(S_i) h(S_i, S_j)$$
The factor of $h(S_i, S_j)$ serves to effectively restrict the summation to the subsets of $S_j$. Due to the partial ordering, for every nonzero term in the sum, the size of the set $S_i$ is greater than or equal to the size of the set in the term that came before. The value of $V$ for the empty set (zero) comes first, followed by the values of $V$ for each subset of $S_j$ that contains only one element. These are the self-energies of each molecule within $S_j$. The values of $V$ for subsets containing two elements, or the pair-energies, come next, and so on. In this intuitive picture, each term contributes to $E(S_j)$ only what the terms corresponding to its own subsets have not already contributed. To verify that this is indeed the case, the existence and uniqueness of $V$ must be established. Then, closed forms can be found for each of the $2^N$ values of $V(S_j)$.
One easy way to find $V$ is to recognize that the $2^N$ equations described above comprise a linear system. For notational convenience, the following definitions can be made: $$e_j = E(S_j) \\ v_j = V(S_j) \\ H_{ij} = h(S_i, S_j)$$
These definitions lead directly to the following set of equations for row vector-matrix multiplication: $$e_j = \sum_{i = 1}^{2^N} v_i H_{ij}$$
$H$ is an upper-triangular matrix, because the partial order-preserving indexing scheme ensures that lower-triangular elements correspond to values of $h(S_i, S_j)$ for sets that do not satisfy the relation $S_i \subseteq S_j$. The determinant of an upper-triangular matrix is the product of its diagonal entries, which are identically $1$ for $H$ because $\forall S_i \in \mathcal{P}(S)$, it holds that $S_i \subseteq S_i$ and thus $h(S_i, S_i) = 1$. Since $H$ has a nonzero determinant, there exists a unique non-trivial solution to the system of linear equations it represents. The solutions to this linear system are precisely the values $V(S_i)$, and so the existence and uniqueness of $V$ are both proven.
Solving for the values of $V(S)$ in terms of the values of $E(S)$, then, amounts to inverting $H$ and taking the following row vector-matrix product: $$v_j = \sum_{i = 1}^{2^N} e_i H^{-1}_{ij}$$
To do so, consider the matrix $H$ as the sum of the $2^N \times 2^N$ identity matrix $I$ and a strictly upper-triangular matrix $U$: $$H = I + U$$
$U$ is a nilpotent matrix by virtue of being strictly triangular. In specific, $U^{2^N + 1} = 0$, where $0$ is the $0$-matrix. Given this fact, it can be verified via explicit multiplication by $H$ that: $$H^{-1} = \sum_{k = 0}^{2^N} (-1)^{k} U^k$$
The above equation can be considered in analogy to the geometric series for scalars. Each entry of $H^{-1}_{ij}$, then, can be computed straightforwardly by summing the entries of $(U^k)_{ij}$ with alternating signs. By the definition of $U$, the entries of $U$ are almost identical to those of $H$ except the diagonal elements are set equal to 0. This is equivalent to defining the entries of $U$ in terms of a slightly different indicator function than $h$ that does not tolerate set equality. In other words: $$U_{ij} = h'(S_i, S_j)$$ $$h'(A, B) =\begin{cases}1 & \textrm{if} \ A \subset B \\ 0 & \textrm{if} \ B \subseteq A \\ \end{cases}$$
The matrix entries $(U^k)_{ij}$, in turn, can be written as explicit matrix products as follows: $$(U^k)_{ij} = \sum_{l_1 = 1}^{2^N} \cdots \sum_{l_{k - 1} = 1}^{2^N} U_{il_1} U_{l_1l_2} \cdots U_{l_{k - 1}j} = \sum_{l_1 = 1}^{2^N} \cdots \sum_{l_{k - 1} = 1}^{2^N} h'(S_i, S_{l_1}) h'(S_{l_1}, S_{l_2}) \cdots h'(S_{l_{k - 1}}, S_j)$$
These matrix elements now have a clear interpretation: they are counts of
possible totally ordered chains of subsets $\{S_{l_1}, S_{l_2}, \cdots S_{l_{k - 1}}\}$ of all of which $S_i$ is a subset and all of which are subsets of $S_j$. The subset $S_{l_1}$ can be created from $S_i$ by adding in some of the elements of $S_{j}$ that are not in $S_i$ (i.e. elements of the set subtraction $S_j \backslash S_i$). Then $S_{l_2}$ can be created by adding more of those elements. This process is repeated a total of $k$ times to obtain $S_j$. Accordingly, the problem of determining the number of totally ordered chains between $S_i$ and $S_j$ is equivalent to the problem of counting the number of unique partitions of the set subtraction $S_j \backslash S_i$ into $k$ subsets, where the order of subsets matters. In addition to $k$, this count will depend on the size of the set substraction: $\lvert S_j \backslash S_i \rvert = \lvert S_j \rvert - \lvert S_i \rvert$. Clearly, if $\lvert S_j \rvert - \lvert S_i \rvert \leq 0$, this count of partitions must be 0 because the set subtraction is empty.
The case where the order does not matter is solved by the Stirling numbers of the second kind (https://en.wikipedia.org/wiki/Stirling_numbers_of_the_second_kind), which account for the redundancy of equivalent orderings of the partitions by dividing by a prefactor. The corresponding case for which order does matter is given as: $$(U^k)_{ij} = W(\lvert S_j \rvert - \lvert S_i \rvert, k) = \sum_{l = 0}^{k} (-1)^l {k \choose l} {(k - l)}^{\lvert S_j \rvert - \lvert S_i \rvert}$$
For $k > 0$ and $\lvert S_j \rvert - \lvert S_i \rvert > 0$, these counts are tabulated nicely in upper right triangle form at OEIS (https://oeis.org/A019538/table). Recalling that $H^{-1}_{ij} = \sum_{k = 0}^{2^N} (-1)^k (U^k)_{ij}$, it follows that: $$H^{-1}_{ij} = I_{ij} + \sum_{k = 1}^{2^N} (-1)^k W(\lvert S_j \rvert - \lvert S_i \rvert, k)$$
This equation is equivalent to summing over the ${(\lvert S_j \rvert - \lvert S_i \rvert)}^{\textrm{th}}$ column of the upper right triangle at OEIS (indexing columns from $0$) with alternating positive and negative signs and values of 0 for the lower triangular entries. This sum always gives $(-1)^{\lvert S_j \rvert - \lvert S_i \rvert}$ when $\lvert S_j \rvert - \lvert S_i \rvert \geq 0$.
Putting everything together, the inverse matrix can be written succinctly as:
$$H^{-1}_{ij} = H_{ij} (-1)^{\lvert S_j \rvert - \lvert S_i \rvert}$$
The introduction of $H_{ij}$ ensures that each pair of $S_i$ and $S_j$ satisfy the foundational assumption that $S_i \subseteq S_j$ for nonzero entries of $H^{-1}$. As a result: $$v_j = \sum_{i = 1}^{N^2} e_i H_{ij} (-1)^{\lvert S_j \rvert - \lvert S_i \rvert}$$ $$V(S_j) = \sum_{i = 1}^{N^2} E(S_i) h(S_i, S_j) (-1)^{\lvert S_j \rvert - \lvert S_i \rvert}$$ $$V(S) = \sum_{T \in \mathcal{P}(S)} E(T) (-1)^{\lvert S \rvert - \lvert T \rvert}$$
The final expression makes use of the properties of the indicator function $h$ to rewrite a restricted sum as a sum over the power set of $S_j$, then substitutes $S_j \to S$ and $S_i \to T$ as before. This final expression also justifies the functional forms provided in the original post. In the set-centric notation used throughout this derivation, the function $V_n$ is defined as follows: $$V_n(R) = \sum_{S \in \mathcal{P}(R) \\ \& \lvert S \rvert = n} V(S) = \sum_{S \in \mathcal{P}(R) \\ \& \lvert S \rvert = n} \sum_{T \in \mathcal{P}(S)} E(T) (-1)^{\lvert S \rvert - \lvert T \rvert}$$
