Let $P(n)=(n+r_1)(n+r_2)...(n+r_k)$ be a polynomial with simple, rational, negative roots (i.e. $r_i>0$) and degree $k\geq 2$ (I stick with negative roots as I don't have to worry about dividing by $0$). My question is: What can we say about the dimension of the set
$$S=\text{span}\left\{\sum_{n=0}^\infty \frac{1}{P(n)},\sum_{n=0}^\infty \frac{n}{P(n)},...,\sum_{n=0}^\infty \frac{n^{k-2}}{P(n)}\right\}$$
over $\mathbb{Q}$ in terms of the roots $r_1,r_2,...,r_k$?
Work so far: For the case $k=2$, it is not hard to show that
$$\sum_{n=0}^\infty \frac{1}{(n+a)(n+b)}=\frac{\psi(a)-\psi(b)}{a-b}$$
(where $\psi(x)$ is the digamma function). Since this is an increasing function for $x>0$ we know
$$\sum_{n=0}^\infty\frac{1}{(n+a)(n+b)}=\frac{\psi(a)-\psi(b)}{a-b}\neq 0$$
and therefore the dimension can only be $1$. In fact, we can use this fact to prove the dimension is never $0$. Assume by way of contradiction that there exists a polynomial $P(n)=(n+r_1)(n+r_2)+...+(n+r_k)$ such that
$$\dim\left(\text{span}\left\{\sum_{n=0}^\infty \frac{1}{P(n)},\sum_{n=0}^\infty \frac{n}{P(n)},...,\sum_{n=0}^\infty \frac{n^{k-2}}{P(n)}\right\}\right)=0$$
(i.e. all of the elements are $0$). Then pick rationals $a_i$ such that
$$Q(n)=\sum_{i=0}^{k-2}a_i n^i=(n+r_3)(n+r_4)+...+(n+r_k)=\frac{P(n)}{(n+r_1)(n+r_2)}$$
Since the dimension of the space is $0$, we know
$$0=\sum_{n=0}^\infty\frac{1}{P(n)}\sum_{i=0}^{k-2}a_in^i=\sum_{n=0}^\infty\frac{Q(n)}{P(n)}=\sum_{n=0}^\infty\frac{1}{(n+r_1)(n+r_2)}\neq 0$$
(which follows from the $k=2$ case). And finally, for the $k=3$ case, a paper by Saradha and Tijdeman proves the dimension is $1$ if and only if all the roots are integer lengths apart (otherwise it is $2$). This pattern extends to higher $k$ as well: For $k\geq 3$, the dimension of $S$ is $1$ if and only if all roots of $P(n)$ are integer lengths apart (this takes more effort to prove but it is indeed true).
For all other $k\geq 4$, the question seems much more difficult.
Motivation: This question is deeply connected to a previous (false) conjecture that I had (found here) which concerned the rationality of infinite sums of rational functions. This question could be thought of as a modification of that original question as I tried to nail down the important behavior of these objects.
No idea for the dimension but each of your series has a closed form. Doing the partial fraction decomposition of $\frac{x^d}{P(x)}$ which has only simple poles and vanishes at $\infty$ we get
$$\sum_{n=1}^\infty \frac{n^d}{P(n)} = \sum_{n=1}^\infty\sum_{m=1}^k \frac{r_m^d/P'(r_m)}{n+r_m}= \sum_{n=1}^\infty\sum_{m=1}^k \frac{B r_m^d/P'(r_m)}{B n+B r_m}$$ $$= -\sum_{m=1}^k \sum_{n,1\le Bn+Br_m\le Br_m }\frac{B r_m^d/P'(r_m)}{B n+B r_m}+ \sum_{n=1}^\infty\sum_{l=1}^{B-1} c_{l,d} \frac{ e^{2i\pi ln/B}}{ n}$$ $$=-\sum_{m=1}^k \sum_{n,1\le Bn+Br_m\le Br_m}\frac{B r_m^d/P'(r_m)}{B n+B r_m}- \sum_{l=1}^{B-1} c_{l,d}\log(1-e^{2i\pi l/ B})$$ where $B$ is the product of the denominators of the $r_m$, $\sum_{m=1}^k Br_m^k/P'(r_m) =0$ and $$c_{l,d} = \frac1B \sum_{m=1}^k e^{-2i\pi l r_m} \frac{B r_m^d}{P'(r_m)} $$ are the FFT coefficients.