How to derive the Lambert W function series expansion? How do you use the Lagrange inversion theorem to derive the Taylor Series expansion of W(x)? How else can you derive a series expansion?
 A: While the common way to derive it is by using the Lagrange Inverse Theorem, there technically isn't anything stopping us from making a Taylor Series for it as you would with any other function. As always, we're going to need a list of derivatives. The first one can be found pretty easily via implicit differentiation as follows: 
$$y = W(x)$$
$$ye^y = x$$
$$\frac{d}{dx}\left(ye^y=x\right)$$
$$\left(y+1\right)e^y\cdot\frac{dy}{dx}=1$$
$$\therefore \frac{d}{dx}W\left(x\right)=\frac{e^{-W\left(x\right)}}{\left(W\left(x\right)+1\right)}$$
$$\frac{d}{dx}W\left(x\right)=\frac{W\left(x\right)}{x\left(W\left(x\right)+1\right)}$$
Also, to anyone unfamiliar with the last step, it is one of the main identities of the Lambert W function. 
Now that we have the first derivative, we can simply differentiate as many times as we want to get all subsequent derivatives. The important thing to note, however, is that all subsequent derivatives will only require the W function to be evaluated at x. Therefore, if we know the value of W(x) we can theoretically calculate the value of any nth derivative of W(x) at that x. We can use this to our advantage by considering a value of W which is easy to calculate, such as W(e). This can be calculated as follows: 
$$y=W\left(e\right)$$
$$ye^y=e=1e^1$$
$$\therefore y=1 \Rightarrow W(e)=1 $$
We now have everything we need to calculate a Taylor Series centered at x=e as usual.  
$$\sum_{n=0}^{\infty}\frac{W^{\left(n\right)}\left(e\right)}{n!}\left(x-e\right)^n$$
$$=1+\frac{1\left(x-e\right)}{2\cdot e\cdot1!}-\frac{3\left(x-e\right)^2}{2^3e^2\cdot2!}+\frac{19\left(x-e\right)^3}{2^5e^3\cdot3!}-\frac{185\left(x-e\right)^4}{2^7e^4\cdot4!}+\frac{2437\left(x-e\right)^5}{2^9e^5\cdot5!}...$$
I agree that this isn't as useful as the Lagrange Inverse method (mostly due to the lack of an explicit definition for the coefficients), however I felt it was worth noting as the question did ask if there were any other methods. And, as shown by the plot of the 5th degree polynomial, this does 
indeed work. 
5th Degree Polynomial Approximation of W(x)
A: I used robjohn's answer to produce an ad hoc solution that didn't involve inequalities or real numbers in particular, but instead produce an identity on formal power series so that it could be used for p-adic numbers for investigating things like $x^x=c$. I chose do that by demonstrating that the series for $W$ that he gives satisfies $W(xe^x)=x$.
$$W(x) = \sum_{k \ge 1} \frac{(-k)^{k-1}}{k!}x^k$$
Composing the power series gets us
$$W(xe^x) = \sum_{k \ge 1} \frac{(-k)^{k-1}}{k!}\left(\sum_{n \ge 1} \frac{1}{(n-1)!} x^n\right) ^k$$
The $n$th power on the series for $xe^x$ corresponds to a convolution of the coefficients with itself $n$ times,
$$W(xe^x) = \sum_{k \ge 1} \frac{(-k)^{k-1}}{k!}\sum_{n \ge 1}\left(\sum_{i_1+\cdots+i_k = n} \frac{1}{(i_1-1)!\cdots(i_k-1)!} \right) x^n$$
Recognizing the coefficients resemble the multinomial expansion, let's focus on it and rewrite it as such
$$\frac{1}{n!} \sum_{i_1+\cdots+i_k = n} i_1 \cdots i_k \binom{n}{i_1,\dots, i_k} $$
The indices multiplying the multinomial coefficient can be seen as coming from differentiating the expansion and then evaluating at 1.
$$\frac{1}{n!} \frac{\partial^k}{\partial x_1 \cdots \partial x_k}\left. \sum_{i_1+\cdots+i_k = n} \binom{n}{i_1,\dots, i_k} x_1^{i_1} \cdots x_k^{i_k} \right|_{x_1=\cdots =x_k = 1}$$
Now we have the multinomial expansion in its entirety, so we can simply write it out,
$$\frac{1}{n!} \left.\frac{\partial^k}{\partial x_1 \cdots \partial x_k} \left( x_1 + \cdots + x_k \right)^n \right|_{x_1=\cdots =x_k = 1}$$
Evaluating these derivatives is straight forward,
$$\left. \frac{1}{n!} \frac{n!}{(n-k)!} \left( x_1 + \cdots + x_k \right)^{n-k} \right|_{x_1=\cdots =x_k = 1}$$
$$\frac{k^{n-k}}{(n-k)!}$$
Let's now replace this more compact form of the convolution back into our series.
$$W(xe^x) = \sum_{k \ge 1} \frac{(-k)^{k-1}}{k!}\sum_{n \ge 1} \frac{k^{n-k}}{(n-k)!} x^n$$
Change order of summation and write in the binomial coefficient and clean up terms slightly,
$$W(xe^x) = \sum_{n \ge 1} -\frac{x^n}{n!} \sum_{k \ge 1} \binom{n}{k} (-1)^k k^{n-1} $$
Notice that because the binomial coefficient appears, the inner series is really only a finite sum since it forces the terms to 0 for $k>n$. Let's focus on this.
$$\sum_{k=1}^n \binom{n}{k} (-1)^k k^{n-1}$$
Similar to how we brought in the derivative for the multinomial series, we can think of the $k^{n-1}$ coming from an exponent on a variable $y$ that has been evaluated at 1.
$$\left. \sum_{k=1}^n \binom{n}{k} (-1)^k k^{n-1}y^k \right|_{y=1}$$
This time we need to multiply by $y$ after differentiating to keep the power at $k$.
$$\left. \left(y \frac{d}{dy} \right)^{n-1} \sum_{k=1}^n \binom{n}{k} (-1)^ky^k \right|_{y=1}$$
When $n=1$ there is no derivative taken, so the result is just $-1$ for this sum. For $n>1$ we can safely put in the $k=0$ constant term which gets removed by the derivative so that we can rewrite it now as a binomial expansion,
$$\left. \left(y \frac{d}{dy} \right)^{n-1} (1-y)^n \right|_{y=1}$$
Because we are taking only up to the $n-1$ derivative of $(1-y)^n$ there will always be a $(1-y)$ term in the product rule expansion. This means evaluating at $y=1$ forces all of these terms to $0$. This means the only surviving term was the $n=1$ term which gave $-1$, so we return to our series and find that it simplifies to,
$$W(xe^x) = x$$
A: It can be easily derived by applying the Lagrange inversion formula to $ze^z$ as it satisfies all criterias of the Lagrange inversion theorem, i.e. $f(0)=0$ and $f'(0)$ is non zero.
