Define
$$A_n:=\frac{(a_1)_n\cdot\ldots\cdot(a_{q+1})_n}{(b_1)_n\cdot\ldots\cdot(b_q)_n}\frac{1}{n!}$$
Notice that
\begin{align}
\Big(\frac{|A_{n+1}|}{|A_n|}\Big)^n&=\Big(\frac{\prod^{q+1}_{k=1}|a_k+n|}{(n+1)\prod^p_{j=1}|b_j+n|}\Big)^n\\
&=\Big|1+\frac{a_{q+1}-1}{n+1}\Big|^n\prod^q_{k=1}\Big|1+\frac{a_k-b_k}{n+b_k}\Big|^n\xrightarrow{n\rightarrow\infty}e^{\mathfrak{R}(a_{q+1}-1)}\prod^q_{k=1}e^{\mathfrak{R}(a_k-b_k)}\\
&=\exp\Big(\mathfrak{R}\big(\sum^{q+1}_{k=1}a_j-\sum^q_{k=1}b_k -1\Big)\Big)
\end{align}
Consequently, $\sum_n|A_n|$ converges if $\mathfrak{R}\big(\sum^{q+1}_{k=1}a_j-\sum^q_{k=1}b_k\big)<0$; $\sum_n|A_n|$ diverges if $\mathfrak{R}\big(\sum^{q+1}_{k=1}a_j-\sum^q_{k=1}b_k\big)>0$; the test is inconclusive if $\mathfrak{R}\big(\sum^{q+1}_{k=1}a_j-\sum^q_{k=1}b_k\big)=0$.
Thus is based on the the following not well known convergence test:
Theorem: Suppose $a_n>0$ for all $n$.
- If $\limsup_n\Big(\frac{a_{n+1}}{a_n}\Big)^n<\frac1e$, then the series $\sum_na_n$ converges.
- If there is $N$ such that for all $n\geq N$, $\Big(\frac{a_{n+1}}{a_n}\Big)^n\geq\frac1e$, then the series $\sum_na_n$ diverges.
An alternative, a proof similar to the one used for $p=2$here works. Notice that
\begin{align}
\frac{A_{n+1}}{A_n}&=\frac{\prod^{p+1}_{k=1}(a_k+n)}{(n+1)\prod^p_{j=1}(b_j+n)}=1+\frac{\prod^{p+1}_{k=1}(a_k+n)-(n+1)\prod^p_{j=1}(b_j+n)}{\prod^{p+1}_{k=1}(a_k+n)}\\
&=1+\frac{\big(\sum^{p+1}_ja_k - \sum^p_kb_k-1)\big)\big)n^p+P(n)}{Q(n)}\tag{1}\label{one}
\end{align}
where $P$ is a polynomial in $n$ of degree at most $p-1$ and $Q$ is a monomial in $n$ of degree $p+1$. It follows immediately that the series $f(z)=\sum_nA_nz^n$ has radius of convergence $1$, hence $f$ is analytic in $B(0;1)$.
Now for $|z|=1$, a further inspection of \eqref{one} yields
$$\frac{A_{n+1}}{A_n}=1+\frac{\xi}{n z_n}+r_n$$
where $\xi=\sum^{p+1}_ja_k - \sum^p_kb_k-1$, $z_n\xrightarrow{n\rightarrow\infty}1$, and $nr_n\xrightarrow{n\rightarrow\infty}0$. Then
\begin{align}
n\Big(1-\Big|\frac{A_{n+1}}{A_n}\Big|\Big)&=n\Big(1-\big|1+\frac{\xi}{n z_n}+r_n\big|\Big)\\
&=n\left(\frac{1- \big|1+\frac{\xi}{n z_n}+r_n\big|^2}{1+ \big|1+\frac{\xi}{n z_n}+r_n\big|}\right)\\
&=n\left(\frac{-2\mathfrak{R}(\frac{\xi}{n z_n}+r_n)-\tfrac{1}{n^2}\mathfrak{R}^2(\frac{\xi}{z_n}+nr_n) +\frac{1}{n^2}\mathfrak{I}^2(\frac{\xi}{z_n}+nr_n)}{1+ \big|1+\frac{\xi}{n z_n}+r_n\big|} \right)
\end{align}
Hence
$$ n\Big(1-\Big|\frac{A_{n+1}}{A_n}\Big|\Big)\xrightarrow{n\rightarrow\infty}-\mathfrak{R}(\xi)=-\mathfrak{R}\Big(\sum^{p+1}_ja_k - \sum^p_kb_k-1\Big)=1+\mathfrak{R}\big(\sum^p_kb_k-\sum^{p+1}_ja_k\big)
$$
The conclusion then follows by the Raabe's theorem: if $-\mathfrak{R}(\xi)>1$, $\sum_nA_ne^{in \theta}$ converges absolutely; when $-\mathfrak{R}(\xi)<1$, $\sum_n|A_n|$ diverges; when $\mathfrak{R}(\xi)=1$ the test is inconclusive.
