How can I prove the identity $2(n-1)n^{n-2}=\sum_k\binom{n}{k}k^{k-1}(n-k)^{n-k-1}$?

How can I prove the identity

$$2(n-1)n^{n-2}=\sum_k\binom{n}{k}k^{k-1}(n-k)^{n-k-1}?$$

I know that the number of trees on $n$ vertices is $n^{n-2}$, and that a tree with $n$ vertices has $n-1$ edges, but I don't know how to continue. Could you help me please?

The OP is in the correct direction. By Cayley's formula, the number of labeled trees on $n$ vertices is $n^{n-2}$. We will count the number of labelled trees in a different way and get the equality.

Take a tree $T$ and an edge $e \in T$. Then $T \setminus \{ e \}$ contains two components, say of sizes $k$ and $n-k$. Also, each of the components is a tree by itself. With this in mind, we build the tree $T$ as follows. Partition the vertices into two parts, of sizes $k$ and $n-k$ respectively. Draw trees $T_1$ and $T_2$ on the two parts. Then pick an arbitrary vertex from each of the parts and join the pair by an edge.

The number of ways of doing this procedure is: $$ \sum_{k=1}^{n-1} \binom{n}{k} \times k^{k-2} (n-k)^{n-k-2} \times (k (n-k)) = R, $$ the right hand side of the purported equation. The binomial coefficient arises from choosing $k$ vertices. The second factor is by applying the Cayley's formula for the number of trees on $k$ and $n-k$ vertices. The $k(n-k)$ factor is the number of ways to connect these components by a single edge.

Now, this procedure overcounts the number of trees. Specifically, each tree can be cut in any of its $n-1$ edges to give a pair of components. (That's a factor of $n-1$.) Further, in selecting $k$ vertices to make two components of sizes $k$ and $n-k$, we distinguish the two components (one that is selected, and one that is not). So we should further divide the number we got by $2$.

Thus the number of labelled trees on $n$ vertices is $$ \frac{R}{2(n-1)}. $$ Equating it to $n^{n-2}$, we get the claim.

The Lambert function has the following Maclaurin series:

$$-W(-x)=\sum_{k=1}^\infty \frac{k^{k-1}}{k!} x^k$$

(In some references, $-W(-x)$ is referred to as the "tree function", $T(x)$, as it is a generating function for rooted labeled trees.)

This series can be derived through Lagrangian inversion. In particular, the coefficient of $x^k$ in the power series of the tree function is given by the expression

$$\frac1{k!}\left.\dfrac{\mathrm d^{k-1}}{\mathrm dt^{k-1}}\left(\frac{t}{t\exp(-t)}\right)^k\right\vert_{t=0}=\frac{k^{k-1}}{k!}$$

The original equation can be re-expressed as

$$\frac{2(n-1)n^{n-2}}{n!}=\sum_{k=1}^{n-1} \frac{k^{k-1}}{k!}\frac{(n-k)^{n-k-1}}{(n-k)!}$$

The right hand side is the autoconvolution of the sequence $\dfrac{k^{k-1}}{k!}$; its ordinary generating function is thus the square of the generating function of $\dfrac{k^{k-1}}{k!}$, which is $(-W(-x))^2=W(-x)^2$.

To find an expression for the series coefficient of $W(-x)^2$, we consider first the related function $\left(\dfrac{W(-x)}{-x}\right)^2-1=\exp(-2\,W(-x))-1$. This function is the inverse of $\dfrac{\log\sqrt{1+x}}{\sqrt{1+x}}$. Applying Lagrangian inversion to $\dfrac{\log\sqrt{1+x}}{\sqrt{1+x}}$ and using the results obtained in this answer, we obtain the series

$$\left(\frac{W(-x)}{-x}\right)^2-1=\sum_{k=1}^\infty \frac{2(k+2)^{k-1}}{k!}x^k$$

which rearranges to

$$W(-x)^2=\sum_{k=0}^\infty \frac{2(k+2)^{k-1}}{k!}x^{k+2}$$

and can also be expressed as

$$W(-x)^2=\sum_{k=1}^\infty \frac{2(k-1) k^{k-2}}{k!} x^k$$

which proves the simpler expression for the autoconvolution.

The last series also appears as formula 11 of this paper by Corless, Jeffrey, and Knuth, as well as appearing in a disguised form as equation 5.60 in Concrete Mathematics.