I add two hopefully useful remarks.
I consider the general situation. In the sequel $k\geq 1$ and $\lambda=(\lambda_1,\ldots,\lambda_{k+1})$ are fixed, $s:=\sum_{i=1}^{k+1}\lambda_i$ and $n\geq s$.
Let $A_n^{(k+1)}(\lambda):=\sum\limits_{\stackrel{a_i\geq \lambda_i, i=1,\ldots,k+1}{a_1+\ldots +a_{k+1}=n}} \prod_{i=1}^{k+1} a_i!$
and let $\mathcal{S}_{k}:=\{ (x_1,\ldots,x_{k+1})\;:\;x_i\geq 0, \sum_{i=1}^{k+1} x_i=1\}$ denote the $k$-dimensional standard simplex.
(i) $A_n^{(k+1)}(\mathbf{\lambda})$ has a nice geometric representation:
$$ A_n^{(k+1)}(\lambda)=\frac{(n-s+k)!\,(n+k)!}{(n-s)!} \int_{\mathcal{S}_{k}} \int_{\mathcal{S}_{k}} \left(x_1y_1+\ldots+x_{k+1}y_{k+1}\right)^{n-s}\,\prod_{i=1}^{k+1} x_i^{\lambda_i}d\mathbf{x}\,\,d\mathbf{y} $$
Proof: recall
Dirichlet's integral: $$((\sum_{i=1}^{k+1}a_i)+k)!\, \int_{\mathcal{S}_{k}} \prod_{i=1}^{k+1} x_i^{a_i}\,d\mathbf{x}=\prod_{i=1}^{k+1} a_i!$$
Now expand the integrand using the multinomial theorem and use Dirichlet's integral repeatedly.
EDIT: I have replaced the previous (unnecessarily complicated) derivation and again corrected the factor (apologies). (The conclusions remain unchanged).
(ii) Thus the large $n$ asymptotic of $A_n^{(k+1)}(\lambda)$ can be investigated along the lines of the Laplace method (explained e.g. in chapter 4 of de Bruijn's book).
Remark: I have done that, and (unfortunately?) the result confirmed your original conjecture. Thus either my analysis or the accepted answer (or both) are incorrect (no offence).
EDIT:
The Laplace analysis follows the usual scheme:
(1) identify maxima
(2) cutoff tails, approximate integrand (3) complete tails
I sketch the main steps. Let $I$ denote the integral above.
Basic intuition:
(1) Consider the scalar product $\langle \mathbf{x},\mathbf{y}\rangle>=\sum_{i=1}^{k+1}x_iy_i\leq 1$. On the domain of integration $\frac{1}{k+1}\leq \langle\mathbf{x},\mathbf{y}\rangle\leq 1$, and $\langle\mathbf{x},\mathbf{y}\rangle$ can be 1 only in a "corner" $x_iy_i=1$ for a certain $i$.
(2) For any $0<\epsilon<1$ the part of $I$ where $\langle\mathbf{x},\mathbf{y}\rangle<1-\epsilon$ is asymptotically negligible (compared to its complement). Asymptotically thus only the parts around the corners (i.e. $x_iy_i\approx 1$) need to be considered.
With this in mind:
(3) choose $\epsilon$ so small that the corners $C_i:=\{ x\in \mathcal{S}_k\;:\; x_i\geq 1-\epsilon, y_i\geq 1-\epsilon\}$ are disjoint, and discuss their influence individually.
E.g. for corner $k+1$ choose $x_1,\ldots,x_k$ resp. $y_1,\ldots,y_k$ as independent coordinates.
Then $x_{k+1}=1-\sum_{i=1}^k x_i=:1-u$, $y_{k+1}=1-\sum_{i=1}^k y_i=:1-t$.
We have $\langle\mathbf{x},\mathbf{y}\rangle=(1-u)(1-t)+\sum_{i=1}x_iy_i\leq (1-u)(1-t)+ut$, thus $(\langle\mathbf{x},\mathbf{y}\rangle)^2\leq (1-2u(1-u))(1-2t(1-t))$.
The parts of the integral where $u\geq 1/n^{2/3}$ or $t\geq 1/n^{2/3}$
are negligible because $\langle\mathbf{x},\mathbf{y}\rangle^n=\mathcal{O}(\exp(-n^{1/3}))$ there.
On the remaining part $\langle\mathbf{x},\mathbf{y}\rangle^n =\exp(-n(u+t))(1+o(1))$ uniformly.
Therefore introduce new coordinates $x_1,\ldots,x_{k-1},u$, $y_1,\ldots,y_{k-1},t$, then
\begin{align*}
&\int_{{C}_{k+1}}\int \left(x_1y_1+\ldots+x_{k+1}y_{k+1}\right)^n\,\prod_{i=1}^{k+1} x_i^{\lambda_i}d\mathbf{x}\,\,d\mathbf{y} \\
&=\int\limits_{{x_i\geq 0, x_1+\ldots +x_k=u}\atop{{y_i\geq 0, y_1+\ldots +y_k=t}\atop {0\leq u,t\leq n^{-2/3}}}} e^{-n(u+t)}\prod_{i=1}^k x_i^{\lambda_i}(1-u)^{\lambda_{k+1}} dx_1\ldots,dx_k dy_1\ldots dy_k\,du\,dt\big(1+o(1)\big)\\
&= \int_0^{n^{-2/3}}\int_0^{n^{-2/3}} e^{-n(u+t)} (\frac{u}{n})^{s-\lambda_{k+1}+k-1} \frac{\prod_{i=1}^k\lambda_i!}{(s-\lambda_{k+1}+k-1)!} (1-u)^{\lambda_{k+1}}
(\frac{t}{n})^{k-1}\frac{1}{(k-1)!}\,du\,dt\big(1+o(1)\big)\\
&=\frac{\prod_{i=1}^k \lambda_i!}{n^{2k+s-\lambda_{k+1}}}\left(1+o(1)\right)
\end{align*}
where last line follows after substituting $z=nu,w=nt$ and completing the tails.
Summing over all corners one now gets your formula (after applying Stirling's formula to it, and up to the additional factors).