Revamped Feb. 12, 2022:

I posted an answer to this (perennial) question in detail in the old MO-Q "Formula for n-th iteration of dx/dt=B(x)" and pointed out a common conflation of related but distinct number arrays, all related to 'natural growth' of rooted trees (and therefore Lehmer codes, see Adler ref below, p. 12). With decreasing order of refinement, they are

**1)** OEIS A139002: the Connes-Moscovici weights of the Connes-Moscovici Hopf algebra, enumerating forests of 'naturally-grown' rooted trees (see Hivert et al. in Taylor's answer, eqns. 34 and 36-40)

**2)** OEIS A139605: coefficients for the Scherk-Comtet partition polynomials, the normal-ordered operator expansion of the diff op $(f(z)\partial_z)^n$ in terms of the derivatives of $f$ and the derivative operator

**3)** OEIS A145271: the refined Eulerian numbers, coefficients for the expansion of $(f(z)\partial_z)^n \; f(z)$, or, equivalently, $(f(z)\partial_z)^{n+1} \; z$, as a partition polynomial in terms of the derivatives of $f(z)$ (see Hivert et al., eqns. 15-20). **This is the expansion the OP, M.G., is addressing**.

The paper by Hivert et al. presented by Peter Taylor addresses item 1 and illustrates item 3; however, there is no discussion of a multinomial-type formula for the refined Eulerian numbers--only eqns. 28 and 35 are given, essentially stating that item 3 is a coarsening/reduction of item 1. As far as I can tell, beyond eqn. 40, Hivert et al. deal with trees in relation to generic statistics generated by various codes and don't come back to the illustration of item 3. The Findstat entry created by Hivert that Taylor links to does give examples of the refined Eulerian numbers generated by the Lehmer code algorithm. This is an algorithm, rather than a formula of the multinomial-type (correct me if I'm mistaken)--formulas that are available for item 1 and for specialized Lagrange inversion partition polynomials (see notes III and IV below) generated from item 3, but not for item 3 directly. Giving the algorithm rather than a formula is analogous to saying a number array counts the number of perfect matchings of the vertices of the n-simplices and giving an algorithm that constructs the matchings and then counts them, or saying the number array enumerates a sum over the non-intersecting dissections of the convex polygons and giving a construction algorithm, without giving a numerical formula in terms of the number of vertices of the constructs. A direct formula for the refined Eulerian numbers is still an open question, I believe.

The paper by Domininci that skbmoore (any relation to eecummings?) cites deals with related but distinct expansions, giving as the most general example the associahedra partition polynomials I point out below. Although very important, this is not the expansion the OP is addressing. (Dominici doesn't address the connections to the associahedra nor flow equations nor does he go into any detail on trees or other combinatorial constructs, giving just a ref. on trees)

Additional refs and further notes:

I) "Set partitions and integrable hierarchies" by Adler presents the first few partition polynomials of the expansion the OP desires on p. 11 and discusses them in the context of 'natural growth' sequences $T_n$. He also alludes to Lehmer codes in the context of $T_n$. Cayley's rooted trees and the related Connes-Moscovici weights are other manifestations of 'natural growth'.

II) OEIS A145271: the partition polynomials for $(g(z)\partial_z)^n g(z)$ with $g(z) = 1/f'(z)$ (so just change notation). I called (unaware of the Hivert et al. paper until now) the coefficients of the polynomials the refined Eulerian numbers since they naturally reduce to the Eulerian numbers A008292, or A123125. My blog post "A Creation Op, Scaled Flows, and Operator Identities" contains detail on related flow functions, p.d.e.s, and more as a prelude to introducing the action of the iterated generalized Lie derivative $q(z)+g(z)\partial_z$. See also the recent MO-Q "How are Sheffer polynomials related to Lie theory?" on connections to the Sheffer polynomials, the core polynomials related to the umbral / finite operator calculus.

III) Particular series reps for $f(z)$ in $g(z) = 1/f'(z)$ of the refined Eulerian partition polynomials lead to the classic Lagrange inversion partition polynomials (LIPs) A134685 for compositional inversion of functions and formal Taylor series / e.g.f.s with $f(z)= a_1 z + a_2 \frac{z^2}{2!}+ a_3 \frac{z^3}{3!}+...$ ; the associahedra LIPs A111785 (renormalized A133437), for formal power series / o.g.f.s with $f(z)= b_1 z + b_2 z^2+ b_3 z^3+...$ ; the LIPs A133932 with the log rep $f(z) = c_1 z + c_2 \frac{z^2}{2}+ c_3 \frac{z^3}{3}+...$ ; and finally, but not least, the noncrossing partition LIPs / free cumulant partition polynomials A134264, so prominent in free probability and related quantum theory and random matrix theory, with $f(z) = z/h(z) = z/( h_1 z + h_2 z^2+ h_3 z^3+...)$. {For more info, see, e.g., my answer to the MO-Q "Important formulas in combinatorics" and the recent posts "Ruling the Inverse Universe, the inviscid Hopf-Burgers evolution equation ..." and "A Taste of Moonshine in Free Moments".}

IV) There are multinomial-type expressions for each numerical coefficient of the normal-ordering expansion of the operator $(g(z)\partial_z)^{n}$ (item 2), which are presented in the MO-Q "Differential operator power coefficients". In addition, direct simple multinomial-type expressions exist for the coefficient of any given partition monomial for all the LIPs listed in III. However, as far as I know, no such multinomial-type formula currently exists for the refined Eulerian numbers even though each full partition polynomial can be calculated independently of the others using the corresponding partition polynomials for the other LIPs in III and the permutahedra partition polynomials of A133314 (or the o.g.f. version A263633) or calculated via a matrix computation presented in the MO-Qs "Сlosed formula for $(g\partial)^n$" and "Expansions of iterated, or nested, derivatives, or vectors--conjectured matrix computation".

V) A noncommutative version of A139605 (link therein) was presented by Kentaro Ihara in "Derivations and automorphisms on non-commutative power series".

VI) Around 1853, when Lie was about ten, Charles Graves published the elegant generalized Taylor shift formula

$$ e^{t \; g(z)\partial_z} H(z) = H[f^{(-1)}(f(z)+t)]$$

and, in the 1850s also, published the operator commutator

$$[L,R] =LR-RL= 1$$ (e.g., $L=\partial_z$ and $R=z$) from which the Graves-Lie-Pincherle commutator $$[h(L),R] = \frac{dh(L)}{dL} = h'(L)$$

can be inferred and the dual

$$[L,h(R)] = \frac{dh(R)}{dR} = h'(R).$$

For Sheffer polynomial sequences, $p_n(z)$, the lowering/annihilation/destruction op $L$ and the raising/creation op $R$ are defined by $L \; p_n(z) = n \; p_{n-1}(z)$ and $R \; p_n(z) = p_{n+1}(z)$.

Obviously, $R= g(z)\partial_z$ is the raising op for the partition polynomials formed from $(g(z)\partial_z)^n g(z)$, and, with $g(z)=1/f'(z)$,

$$[g(z)\partial_z,f(z)] = 1,$$

implying, similar to the commutator identities above, that

$$[(g(z)\partial_z)^n,f(z)] = n \; (g(z)\partial_z)^{n-1}.$$

Then we have the functional identity

$$[(g(z)\partial_z)^n,f(z)]g(z) =(g(z)\partial_z)^n f(z)g(z) - f(z)(g(z)\partial_z)^n g(z) = n \; (g(z)\partial_z)^{n-1}g(z).$$

Recall the commutator acting on a function as in

$$[\partial_z,H(z)]K(z) = \partial_zH(z)K(z) - H(z)\partial_zK(z)$$

is the Newton-Leibniz formula in disguise

$$\partial_zH(z)K(z) = H(z)\partial_zK(z) + [\partial_z,H(z)]K(z) = H(z)\partial_zK(z) + K(z)\partial_zH(z).$$

3more comments