I'm looking for combinatorial proofs (using, e.g., trees) of the following particular matrix equation $(I)$ and also combinatoric operational analogs of its solution via matrix inversion and/or Cramer's method with determinants;
$$(I): \; [B][\bar{a}] = [a] $$
and, conversely,
$$[\bar{a}] = [\hat{B}][a] = [B^{(-1)}][a]= [B]^{-1}[a]$$
stemming from the compositional relation
$$\bar{A}(t) = A(\hat{B}(t)) = A(B^{(-1)}(t))= \frac{\hat{B}(t)}{t}= \frac{B^{(-1)}(t)}{t},$$
where
$$\hat{B}(t)= B^{(-1)}(t)$$
is the compositional inverse of the analytic function or formal power series $B(t)$, $[B]$ is the lower triangular matrix with the row polynomials given by the e.g.f. $e^{B(t)x}$, and $[a]$ and $[\bar{a}]$ are column vectors containing the coefficients of the formal Taylor series
$$A(t) = e^{a.t} = \frac{t}{B(t)}$$
and
$$\bar{A}(t) = e^{\bar{a}.t} = \frac{B^{(-1)}(t)}{t}.$$
With $B(t) = e^t-1$, these dual matrix equations relate the Bernoulli numbers to the sequence $(-1)^n \;n!/(n+1)$ via the Stirling numbers of the second and first kinds with e.g.f.s $e^{(e^t-1)x}$ and $e^{\ln(1+t)x}=(1+t)^x$, respectively.
I use the notation $[R]$ and $[r]$ to denote, respectively, a square matrix and a column vector of finite or infinite extent; $R(t)$, an associated function analytic at the origin with $t$ an independent variable (real or complex), from which $[R]$ or $[r]$ is generated; $R_n(x)$, a polynomial (or polynomials) of degree $n$ in the independent variable $x$; $R_{n,k} = [R]_{n,k}$, the coefficient of $x^k$ of $R_n(x)$; and $r_n$, the coefficient of the (possibly formal) Taylor series for $R(t)$ and the $n$'th component of $[r]$. To make formulas and derivations brief, readable, and intuitive, I use a dot as a subscript to denote an umbral character $U.$ so that in an identity reduced ultimately to an analytic series the umbral variable $(U.)^n$ can be reduced to $U_n$ regardless of the nature of $U_n$ as an indeterminate, polynomial, etc.
Background analysis:
This section is somewhat motivational, providing a more general result, and might be helpful, but you can skip it and go to the next section for a briefer derivation of the matrix equation $(I)$.
Given the parent function (or formal Taylor series), with $t$ a real or complex independent variable and $(b.)^n = b_n$ the indeterminates with the dot as a subscript indicating an umbral character,
$$B(t) = e^{b.t} = \sum_{n \ge 0}\frac{(b.t)^n}{n!} = \sum_{n \ge 0}b_n \frac{t^n}{n!} $$
with $b_0=0$ and $b_1 \ne 0$ and its formal compositional inverse about the origin
$$B^{(-1)}(t) = \hat{B}(t) = e^{\hat{b}.t},$$
construct the following Sheffer polynomial sequences;
I) the parent binomial Sheffer sequence
$$e^{B.(x)t} = e^{B(t)x}$$
II) its umbral inverse binomial Sheffer sequence
$$e^{\hat{B}.(x)t} = e^{\hat{B}(t)x} = e^{B^{(-1)}(t)x}$$
III) the Appell Sheffer sequence associated to the parent function
$$e^{A.(x)t} = \frac{t}{B(t)} \; e^{xt} = A(t) \; e^{xt} = e^{a.t}\; e^{xt} = e^{(a.+x)t}$$
IV) the conjugate Appell Sheffer sequence
$$e^{\bar{A}.(x)t} = \frac{\hat{B}(t)}{t} \; e^{xt} = \frac{B^{(-1)}(t)}{t} \; e^{xt}= \bar{A}(t) \; e^{xt} = e^{\bar{a}.t}\; e^{xt} = e^{(\bar{a}.+x)t}.$$
Then
$$ e^{\hat{B}.(A.(B.(x)))t} = e^{A.(B.(x))\hat{B}(t)}= e^{(a.+B.(x))\hat{B}(t)}$$
$$ = e^{a.\hat{B}(t)} \; e^{B.(x)\hat{B}(t)} = A(\hat{B}(t)) \; e^{xB(B^{(-1)}(t))}$$
$$ = \frac{t}{B(t)}\; |_{t \rightarrow \hat{B}(t)=B^{(-1)}(t)} \; e^{xt} =\frac{B^{(-1)}(t)}{t} \; e^{xt}$$
$$=\frac{\hat{B}(t)}{t} \; e^{xt} = \bar{A}(t) \; e^{xt},$$
so stemming from the compositional relation
$$\bar{A}(t) = A(\hat{B}(t)) = A(B^{(-1)}(t))= \frac{\hat{B}(t)}{t}$$
is the umbral compositional conjugation
$$\hat{B}_n(A.(B.(x))) = \bar{A}_n(x),$$
which is equivalent to the matrix conjugation
$$[\hat{B}][A][B] =[B^{(-1)}][A][B] = [B]^{-1}[A][B] = [\bar{A}].$$
Since $B_n(0) = [B]_{n,0} = \delta_n$, i.e., the constant terms of the binomial Sheffer polynomials are zero except for $B_0(x)=1$,
$$[\bar{a}] = [\hat{B}][a] = [B^{(-1)}][a]= [B]^{-1}[a]$$
and, conversely,
$$[a] = [B][\bar{a}] .$$
Given the parent function $B(t)$, the matrix $[B]$ and the column vector $[a]$ can be constructed, and then Cramer's rule can be used to construct the column vector $[\bar{a}]$, containing the normalized coefficients of the formal Taylor series of the inverse function $B^{(-1)}(t)$.
These last two matrix equations are the matrix reps for a Faa di Bruno-type algebra dealing with the compositions
$$\bar{A}(t) = A(\hat{B}(t))$$
and
$$A(t) = \bar{A}(B(t)).$$
Caveat: In order for the Appell sequences to have the canonical form with $A_0(x) =a_0= 1$, the imposed condition $A(t)=\frac{t}{B(t)}$ forces $B'(0)=b_1 =1$, but these unity conditions are not necessary in establishing the conjugation relations and $B'(0)=b_1$ is assumed to be only finite and non-zero below.
Composition partition polynomials
The matrix equation
$$[B][\bar{a}] = [a]$$
represents the functional composition (with $D = d/dz$)
$$A(t) = \bar{A}(B(t)) = e^{B(t)D_{z=0}}\bar{A}(z) = e^{tB.(D_{z=0})}\bar{A}(z)$$
$$=e^{tB.(D_{z=0})}e^{\bar{a}.z} = e^{tB.(\bar{a}.)} = \sum_{n \ge 0} \frac{t^n}{n!} \;B_n(\bar{a}.),$$
so the row polynomials associated with the matrix $[B]$ are the compositional partition polynomials $Prt_n(b_1,b_2,b_3,...,b_n;x)$ with the e.g.f.
$$e^{Prt.(b_1,...;x)t} = e^{B(t)x}.$$
Now switch reps but use the same symbol for the indeterminates and express $B(t)$ as the formal power series, or o.g.f.,
$$B(t) =\frac{b.t}{1-b.t}= b_1 t + b_2 t^2 + b_3t^3 ...\;,$$
giving matrices with integral factors. (Of course, for a given fixed function the indeterminates evaluate differently here from the indeterminates with the same symbols in the earlier section.) Then the refined Lah partition polynomials (cf. OEIS A130561 and link therein to "Lagrange a la Lah") represent the composition of the exponential function with $B(t)$ and the composition partition polynomials have the e.g.f.
$$e^{B(t)x} =\exp[\frac{b.t}{1-b.t}x] = e^{Lah.(b_1,b_2,...;x)t}.$$
The column vector for the conjugate slightly generalized Appell sequence $\bar{A}_n(x)$ should now be associated with the Lagrange inversion partition polynomials of A133437 for compositional inversion of formal o.g.f.s, which are intimately related to the refined Euler characteristic polynomials of the associahedra. With the additional scale factor,
$$\bar{a}_n = \frac{\hat{b}_{n+1}}{n+1} .$$
The associated slightly generalized Appell partition polynomials in the appropriate indeterminates stem from
$$A(t) = e^{a.t} = \frac{t}{B(t)} = \frac{t}{b_1t+b_2t^2+b3t^3+...}.$$
These partition polynomials for inverting a formal o.g.f. are explored in A263633 (with $x[n]$ in the OEIS replaced by $b_{n+1}$ and $b_1=1$) and are a refinement of the Pascal matrix, but they must be multiplied by factorials to give $a_n$, which are coefficients for an e.g.f. (With the indeterminates simply scaled, the refined Pascal polynomials are associated with the refined Euler characteristic partition polynomials, or signed refined face polynomials, of the permutahedra.)
As an illustration and for numerical spot checks, the first few rows of the matrices for the matrix equation
$$[B][\bar{a}] = [a]$$
are
$$[B] \rightarrow\begin{bmatrix} 1 & & & &\\ 0 & b_1 & & &\\ 0 & 2b_2 & b_1^2 & &\\ 0 & 6b_3 & 6b_1 b_2 & b_1^3 &\\ 0 & 24b_4 & 24b_1 b_3+12b_2^2 & 12b_1^2 b_2 & b_1^4 \end{bmatrix}.$$
$$[\bar{a}] \rightarrow\begin{bmatrix} \frac{1}{b_1}(1) & \\ \frac{1}{2b_1^3}(-2b_2) & \\ \frac{1}{3b_1^5}(12b_2^2-6b_1b_3) & \\ \frac{1}{4b_1^7}( -120 b_2^3 + 120 b_1b_2b_3 - 24b_1^2b_4) & \\ ... & \end{bmatrix}$$
$$=\begin{bmatrix} \frac{1}{b_1}(1) & \\ \frac{1}{b_1^3}(-b_2) & \\ \frac{1}{b_1^5}(4b_2^2-2b_1b_3) & \\ \frac{1}{b_1^7}( -30 b_2^3 + 30 b_1b_2b_3 - 6b_1^2b_4) & \\ ... & \end{bmatrix}$$
$$[a]\rightarrow \begin{bmatrix} \frac{1}{b_1}(1) & \\ \frac{1}{b_1^2}(-b_2) & \\ \frac{2!}{b_1^3} (b_2^2-b_1b_3) & \\ \frac{3!}{b_1^4}( -b_2^3 + 2 b_1b_2b_3 - b_1^2b_4) & \\ ... & \end{bmatrix}$$
$$= \begin{bmatrix} \frac{1}{b_1}(1) & \\ \frac{1}{b_1^2}(-b_2) & \\ \frac{1}{b_1^3} (2b_2^2-2b_1b_3) & \\ \frac{1}{b_1^4}( -6b_2^3 + 12 b_1b_2b_3 - 6b_1^2b_4) & \\ ... & \end{bmatrix}.$$
Various associations among these elements follow from the application of matrix inversion and Cramer's rule for solving matrix equations with determinants, operations that have the potential for combinatorial interpretations.
The use of Cramer's method to solve for $\bar{a}_n$ using $a_n$ and $[B]$ is equivalent to the use of the noncrossing partitions A134264 to solve for $B^{(-1)}(t)$ using the coefficients of the shifted reciprocal $t/B(t) = e^{a.t}$ expressed as an o.g.f., i.e., using $a_n/n!$.
