June 18, 2019 Edit: The identity is now proved and generalized. See below.
I did not have time to look at your identity for very long but I am pretty sure it follows from Lemma II.2 in my article "Trees, forests and jungles: a botanical garden for cluster expansions" with Vincent Rivasseau.
A starting point. Make a set of $n+1$ elements by adding a root numbered say $0$. Define the edge variables $u_{ij}$ $(i<j)$ to be zero except $u_{0,i}=z_i$. An ordered forest of length $\tau$, i.e., $(\{i_1,j_1\},\{i_2,j_2\},\ldots,\{i_{\tau},j_{\tau}\})$ which contributes in Lemma II.2 namely for which
$$
\prod_{\nu=1}^{\tau}u_{i_{\nu},j_{\nu}}\neq 0
$$
basically is sequence of distinct elements $(j_1,\ldots,j_{\tau})$ in $\{1,2,\ldots,n\}$. And for maximal length $\tau=n$, this is just a permutation $\sigma$.
Not quite there yet, although it looks promising.
I use formula (II.6) in Lemma II.2 of my article for the set $\{0,1,\ldots,n\}$ instead of $\{1,\ldots,n\}$. In the notations of the paper, I define the edge variables $u_{\{ij\}}$ by $u_{\{0i\}}=z_i$ for $1\le i\le n$ and zero otherwise. I also define the edge variables $v_{\{ij\}}$ by $v_{\{0i\}}=y_i$ for $1\le i\le n$ and $1$ otherwise. I will rename $\nu$ from the paper and now call it $j$. I will also rename $\mu$ from the paper and now call it $i$. An ordered forest of length $\tau$ is a sequence of edges $(l_1,\ldots,l_{\tau})$. In order to contribute it has to be of the form $l_a=\{0,\sigma(a)\}$ for $\sigma=(\sigma(1),\ldots,\sigma(\tau))$ a sequence of distinct elements of $\{1,2,\ldots,n\}$ of length $\tau$.
With this harmonization of notations, Lemma II.2 from the paper gives
$$
y_1\cdots y_n=\sum_{\tau=0}^{n}\sum_{\sigma}z_{\sigma(1)}\cdots z_{\sigma(\tau)}
\sum_{j=0}^{\tau}\frac{y_{\sigma(1)}\cdots y_{\sigma(j)}}{\prod_{i=0,i\neq j}^{\tau}
(Z_{j}^{\sigma}-Z_{i}^{\sigma})}
$$
where for $k=0,1,\ldots,\tau$, $Z_{k}^{\sigma}=z_{\sigma(1)}+\cdots+z_{\sigma(k)}$.
Might as well define also for $k=0,1,\ldots,\tau$, $W_{k}^{\sigma}=w_{\sigma(1)}+\cdots+w_{\sigma(k)}$.
I will then carefully separate the cases where $i$ or $j$ are zero from the rest.
$$
y_1\cdots y_n=\sum_{\tau=0}^{n}\sum_{\sigma}z_{\sigma(1)}\cdots z_{\sigma(\tau)}
\sum_{j=1}^{\tau}\frac{y_{\sigma(1)}\cdots y_{\sigma(j)}}{Z_{j}^{\sigma}\prod_{i=1,i\neq j}^{\tau}
(Z_{j}^{\sigma}-Z_{i}^{\sigma})}
$$
$$
+\sum_{\tau=0}^{n}\sum_{\sigma}z_{\sigma(1)}\cdots z_{\sigma(\tau)}
\frac{1}{\prod_{i=1}^{\tau}
(-Z_{i}^{\sigma})}
$$
Now I set $y_i=e^{sz_i+tw_i}$. I take $\frac{\partial^2}{\partial s\partial t}$ and set $s=t=0$.
This gives the identity
$$
(w_1+\cdots+w_n)(z_1+\cdots+z_n)=
\sum_{\tau=0}^{n}\sum_{\sigma}z_{\sigma(1)}\cdots z_{\sigma(\tau)}
\sum_{j=1}^{\tau}\frac{W_{j}^{\sigma}}{\prod_{i=1,i\neq j}^{\tau}
(Z_{j}^{\sigma}-Z_{i}^{\sigma})}
$$
Now I am stuck. I would like to "take the coefficient" of $z_1\cdots z_n$ which would somehow force $\tau=n$ and $\sigma$ to be a permutation but I don't see it yet.
It's a bit late for me so can't give details but I think I know how to peal off $\tau=n$. Use Möbius inversion in the (Boolean) poset of subsets $I=\{i_1<\cdots<i_k\}$ of $[n]=\{1,\ldots,n\}$.
Namely apply the identity I just proved not just for the whole set $\{1,\ldots,n\}$ and corresponding variables, but also for all the $I$'s.
Really really have to go but basically the wanted result amounts to showing
$$
\sum_{I\subset [n]} (-1)^{n-|I|}\left(\sum_{i\in I} w_i\right)
\left(\sum_{i\in I} z_i\right)=0
$$
when $n\ge 3$.
OK here is the endgame.
I will denote by $|\cdot|$ the cardinality of a finite set. For $I\subset[n]$, let
$$
B_I=\sum_{\sigma}z_{\sigma(1)}\cdots z_{\sigma(|I|)}
\sum_{j=1}^{|I|}\frac{W_{j}^{\sigma}}{\prod_{i=1,i\neq j}^{|I|}
(Z_{j}^{\sigma}-Z_{i}^{\sigma})}
$$
where the sum is over sequences $\sigma=(\sigma(1),\ldots,\sigma(|I|))$ which form a permutation of the elements of $I$.
Also define
$$
A_I=\sum_{J\subset I}B_J\ .
$$
Then the previous identity applied to the variables in $I$ instead of the whole set $[n]$, give
$$
A_I=\left(\sum_{a\in I} w_a\right)
\left(\sum_{b\in I} z_b\right)\ .
$$
By Möbius inversion in the Boolean lattice, we get
$$
B_I=\sum_{J\subset I}(-1)^{|I|-|J|}A_J\ .
$$
The LHS in the OP's identity is $B_{[n]}$ divided by $z_1\cdots z_n$.
However,
$$
B_{[n]}=\sum_{J\subset[n]}(-1)^{n-|J|}
\left(\sum_{a\in J} w_a\right)
\left(\sum_{b\in J} z_b\right)
$$
$$
=\sum_{a,b\in [n]}w_a z_b\sum_{J,\{a,b\}\subset J\subset[n]}
(-1)^{n-|J|}=\sum_{a,b\in [n]}w_a z_b(1-1)^{n-|\{a,b\}|}=0
$$
by Newton's binomial theorem and the hypothesis $n\ge 3$. QED
Note that the OP's identity admits a trivial generalization to families of weights $w^{(1)},\ldots,w^{(k)}$, namely
$$
\sum_{\sigma\in\mathfrak{S}_n}
\sum_{j=1}^{n}\frac{W_{j}^{(1),\sigma}\cdots W_{j}^{(k),\sigma}}{\prod_{i=1,i\neq j}^{n}
(Z_{j}^{\sigma}-Z_{i}^{\sigma})}\ =\ 0
$$
as soon as $n\ge k+2$. Here I used the generalized notation $W_{j}^{(r),\sigma}=w_{\sigma(1)}^{(r)}+\cdots+w_{\sigma(j)}^{(r)}$.
Finally more information about identities like in Lemma II.2 of my article with Rivasseau can be found in my dedicated webpage (in dire need of updates though).
