This equation is a particular case of the so-called Underwood equation
$$\sum_{i=1}^n \frac{\alpha_i\, z_i}{\alpha_i- \theta}=1-q$$ where the $\alpha_i> 0$ and $z_i >0$ and $n$ can be very large (potentially up to thousands) and $q$ is given.
With my former research group we spent decades to find efficient numerical methods to solve it since, in chemical engineering, it has to be solved zillions of times in a single simulation.
For sure, as soon as $n>4$, there is no nalytical solutions and numerical methods are required. To me, the key point is to avoid its transform into a polynomial.
Our most recent work was published in $2014$ in this paper (you can also find it here) where we proposed rapid and robust solution methods using convex transformations. Beside, and this is a key point, for any root, we proposed simple and efficient starting guesses which,typically, make that very few iterations are required (this is illustrated in the first figure showing that the starting guess is almost the solution).
If you are not too concerned by computing time, there are simple things you could do. Using you equation, for the root between, say, $n_1$ and $n_2$ which correspond to the vertical asymptotes, transform
$$f(x)=\sum_{k=1}^p\frac{n_k}{x + n_k} - 1$$ as
$$g(x)=(x+n_1)(x+n_2)f(x)$$ which is the most basic form of the so-called Leibovici & Neoschil method which has been widely used for this class of problems during the last $27$ years.
$$g(x)=n_1(x+n_2)+n_2(x+n_1)-(x+n_1)(x+n_2)+(x+n_1)(x+n_2)\sum_{k=3}^p\frac{n_k}{x + n_k}$$ which gives
$$g(-n_1)=n_1(n_2-n_1)\qquad \text{and}\qquad g(-n_2)=-n_2(n_2-n_1)$$ Then, a linear interpolation gives, as an estimate,
$$x_0=-\frac {2\,n_1\,n_2}{n_1+n_2}$$
The problem is that, in the range, $g(x)$ go through an extremum and using directly Newton method is dangerous. A very good way to avoid any problem is to use a method which combine Newton steps and bisection steps (when Newton method tends to lead outside the range).
I would recommend for example subroutine rtsafe from "Numerical Recipes" which does it. The code is here.
This works quite well without any convergence problem (but it is much less efficient than what was proposed in our paper). For $99.9$% of the zillions of cases we worked, only one bisection step is used.
Edit
Revisiting all the things we tried, another possibility is, for each interval, to change variable :
$$x=-n_{1}+\frac{n_{1}-n_2} {1+e^{-t}}$$ and use Newton method for $f(t)$ using
$$t_0=\log \left(\frac{n_2}{n_1}\right)$$
For example, using $p=5$ and $n_i=p_{i+3}$ Newton iterates will be
$$\left(
\begin{array}{cc}
n & t_n \\
0 & 0.4519851 \\
1 & 1.4953866 \\
2 & 1.3605234 \\
3 & 1.3532976 \\
4 & 1.3532792
\end{array}
\right)$$
