Sum of inverse of multinomial coefficients

Your sum is a fancy way of writing $$ A_j=A_j^{(4)}=\sum_{\substack{a_i\ge\lambda_i, i=1,\dots,4\\a_1+\dots+a_4=j}}a_1!a_2!a_3!a_4! $$ (which, of course, can be treated as $A_j^{(k)}$ for any $k\ge2$). If I understand you correctly, your expectation is that $A_j\sim N\times\{\text{the maximum term of the sum}\}$ as $j\to\infty$, where $N$ is the number of hits of this maximum term. This is certainly wrong, as the terms in the neighbourhoods of those maximal entries contribute quite substantially to the sum. Some plausible asymptotics to consider here are $A_{j+1}/A_j$ or $A_j^{1/j}$ as $j\to\infty$, as in these cases one can indeed show that the leading term completely determines the growth. The related reference to study is the book "Asymptotic methods in analysis" by de Bruijn (1961), more specifically, Chapter 3 there. (I would recommend doing $k=2$ and $k=3$ first.) I really recommend this particular book, as the most accessible (and elementary enough) reference to the asymptotics of binomial sums.

Note that if $k_i\geq \lambda_i$ for all $i\in \{1,\dots,l\}$ the term $$ \binom{n}{k_1,\dots ,k_l}^{-1} $$ can only be maximal if for any $i,j\in \{1,\dots,l\}$ with $i\neq j$ we have $k_i=\lambda_i$ or $k_j=\lambda_j$. Hence to get the maximal value we need to put $k_i=\lambda_i$ for all but one $i\in \{1,\dots,l\}$. Hence we have for $n\geq \sum_{i=1}^l \lambda_i$ $$ \binom{n}{k_1,\dots ,k_l}^{-1}\leq \max_{j\in \{1,\dots,l\}} \frac{(n+\lambda_{j}-\sum_{i=1}^l \lambda_i)!\prod_{i=1}^l \lambda_i!}{\lambda_{j}!n!}\leq \frac{(n+\lambda_{max}-\sum_{i=1}^l \lambda_i)!\prod_{i=1}^l \lambda_i!}{\lambda_{max}!n!} $$ A double counting argument yields $$ \sum_{\substack{k_1+\dots+k_l=n\\k_i\geq \lambda_i}} 1=\binom{n-\sum_{l=1}^l\lambda_i+l-1}{l-1} $$ Hence we have $$ \sum_{\substack{k_1+\dots+k_l=n\\k_i\geq \lambda_i}} \binom{n}{k_1,\dots, k_l}^{-1}\leq \binom{n-\sum_{l=1}^l\lambda_i+l-1}{l-1}\frac{(n+\lambda_{max}-\sum_{i=1}^l \lambda_i)!\prod_{i=1}^l \lambda_i!}{\lambda_{max}!n!}\leq \frac{\prod_{i=1}^l \lambda_i!}{(l-1)!\lambda_{max}! }(n+l-1)^{(l-1)-\sum_{i=1}^l \lambda_i+\lambda_{max}}. $$ For example for $\lambda_1=\dots=\lambda_l=2$ we get $$ \sum_{\substack{k_1+\dots+k_l=n\\k_i\geq \lambda_i}} \binom{n}{k_1,\dots, k_l}^{-1}\leq \frac{2^{l-1}}{(l-1)!} (n+l-1)^{-(l-1)}. $$

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).