Let us call $P$ the expression to be computed. As this is a symmetric polynomial, we will use the notation $[M]$, where $M$ is a monomial, to represent the sum of the distinct monomials obtained from $M$ by permuting the indeterminates in all possible ways (this is called an orbit sum). We choose as representative of each class the monomial that comes first in lexicographic order. In a polynomial, we order the orbit sums in lexicographic order.
The advantage of this technique is that, when doing polynomial arithmetic, we only need to compute the coefficient of the representative of each class.
As the polynomial is homogeneous of degree $4$, the only orbit sums involved are $[a^4]$ (which does not actually appear), $[a^3b]$, $[a^2b^2]$, $[a^2bc]$, and $[abcd]$.
We will express the polynomial in terms of the elementary symmetric polynomials:
\begin{align*}
s_1 &= a + b + c + d = [a] = 7\\
s_2 &= ab + ac + ad + bc + bd + cd = [ab] = 3\\
s_3 &= abc + abd + acd + bcd = [abc] = 21\\
s_4 &= abcd = [abcd] = 1
\end{align*}
where the values are obtained using Viète's relations.
The given polynomial is $P = [a^3b] + 2[a^2b^2] + 4[a^2bc] + 9[abcd]\ $; the coefficients are found by computing the number of possible ways to obtain the representative of each orbit sum.
Since the leading term $a^3b$ only contains two variables, we can produce it using $s_1$ and $s_2$. Specifically, we have:
\begin{align*}
s_1^2s_2 &= (a + b + c + d)^2(ab + ac + ad + bc + bd + cd)\\
&= [a^3b] + 2[a^2b^2] + 5[a^2bc] + 12[abcd]
\end{align*}
We obtain, after simplification:
$$
P - s_1^2s_2 = -[a^2bc] - 3[abcd]
$$
Since the leading term $-a^2bc$ only uses three variables, we use $s_1$, $s_2$, and $s_3$ to cancel it. We obtain, after some work:
\begin{align*}
s_1s_3 &= [a^2 b c] + 4[a b c d]\\
P - s_1^2s_2 + s_1s_3 &= [abcd] = s_4
\end{align*}
and this gives:
\begin{align*}
P &= s_1^2s_2 - s_1s_3 + s_4\\
&= 7^2\times 3 - 7\times 21 + 1\\
&= \mathbf{1}
\end{align*}