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) 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" 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) and $p_n(t)$ are the refined Lah / Laguerre polynomials of A130561 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 (the Comtet A polynomials), which acting on $x$ gives A145271 (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.
Edit (Apr 9, 2018): Here's an excerpt from an email to me from Nigel Ray in 2014 concerning Rota's interest in this topic:
"I'm afraid that "Extensions of UC (I)" long predates latex, but here's a link to the Advances
http://www.sciencedirect.com/science/article/pii/0001870886900654
(so long as you have permission). It was Gian-Carlo Rota who suggested I write this up, when I bumped into him one day in Berkeley - he was intrigued by the connection with formal group laws. When I got the chance to explain this to him in more detail (I think it was in Boston) he homed in on the fact that it explained the unresolved question he had always had with binomial sequences - namely how to express their products as linear combinations of themselves. I still think of UC as "formal group laws via Hurwitz Series" (ie divided formal power series), and hope there remains scope for developing that viewpoint."