An equivalent problem is to show the positivity of the connection factors $c^1_{i,j}$ in the expansions

$$p_i(t)p_j(t) = \sum^{i+j}_{n=1}\; c^n_{i,j}p_n(t),$$

where $p_n(t)$ are cycle index partition polynomials of the symmetric groups ([A036039][1]) with the indeterminates $x_n = (-1)^{n-1}h_{n-1}t \;$  and $h_n$ are the complete homogeneous symmetric polynomials with all of their indeterminates positive. The $c^1_{i,j}$ are essentially the coefficients of the FGL expansion Strickland displays. With $(a.)^n = a_n = f^{(n)}(0)\; $ and $\phi_n= n!t_n = e_n$, the elementary symmetric polynomials,  

$$c^1_{j,k} = p_j(a.)p_k(a.)= p_j(t)p_k(t)|_{t^n=a_n}.$$

Jair, in your Sage computations, if the coefficients of f(x) are expressed as its Taylor series coefficients, i.e., they are normalized by the factorials, it is easier to recognize them as A145271, the refined Eulerian numbers. Same for your polynomials p, and if the coefficients of p are grouped together by powers of t and expressed as the elementary symmetric polynomials/functions $e_n=n!t_n=n!\phi_n \;$, it is easy to see they are signed A036039 with the appropriate determinates given above, e.g., $3! p_3 = 2(e_1^2-e_2)t - 3e_1t^2+t^3 = 2h_2t-3h_1t^2+t^3$. 

See "[Formal group laws and binomial Sheffer sequences][2]" for details and examples.

Edit (Feb. 8, 2018):

The computations seem better characterized in terms of $f^{-1}(x)= x - (c_2x^2+c_3x^3+\cdots)$. Then $f(x)=e^{a.x}, \;$ where $a_n/n!$ are the refined face polynomials of the Stasheff associahedra (positive coefficients of [A133437][3]) and $p_n(t)$ are the refined Lah / Laguerre polynomials of [A130561][4] with indeterminates $(x_1,x_2,x_3,..)= (t,-c_2t,-c_3t,..)\;,$ related to the elementary Schur polynomials. Then again $p_n(a.)=0$ and $c^1_{j,k}=p_j(a.)p_k(a.).$

Edit (Feb. 12, 2018):

Further to a conjecture by Majer,

$$f[f^{-1}(x)+f^{-1}(y)]=\exp[f^{-1}(y) \cdot \phi(x)D_x]x = \exp[y 
 \cdot p.(\phi(x)D_x)]x,$$

so

$$\psi_n(x) = p_n(\phi(x)D_x)x.$$

The iterated infinitesimal generator is given in terms of $t_n$ by [A139605][5] (the Comtet A polynomials), which acting on $x$ gives [A145271][6] (the refined Eulerian polynomials, call them $a_n(x)$) all with integer coefficients. The polynomials $p_n(x)$ may be expressed as the the refined Stirling polynomials of the first kind (cycle index polynomials for $S_n$), noted above, or as numerous other composition partition polynomials. Then indeed we have the integer coefficients as Majer feels for

$$\psi_n(x) = p_n(a.(x))$$

with $\psi_n(0)=\delta_{n-1}.$

These polynomials and operators can be represented as various combinatoric structures, so with nice combinatorial interpretations of compositions can the ultimate constructs.


  [1]: https://oeis.org/A036039
  [2]: https://tcjpn.wordpress.com/2018/01/23/formal-group-laws-and-binomial-sheffer-sequences/
  [3]: https://oeis.org/A133437
  [4]: https://oeis.org/A130561
  [5]: https://oeis.org/A139605
  [6]: https://oeis.org/A145271