I would like to suggest an interpretation using super symmetric functions. These are symmetric functions that are symmetric in two sets of variables $\{x_i\}$ and $\{y_j\}$ separately. They satisfy the property that setting a single $x_i$ variable to equal $z$ and a single $y_j$ variable equal to $z$ gives a polynomial independent of $z$, in addition to several other conditions (see Chapter 1, Section 3, Exercise 23 of Macdonald's book). Note that since it's independent of $z$, it's the same as just setting $x_i = y_j = 0$, which amounts to considering the same super symmetric function in a smaller set of variables.
A concrete description of supersymmetric functions is as follows. The supersymmetric functions are the image of the map $\Lambda \to \Lambda \otimes \Lambda$ given by $(\mathrm{Id} \otimes S) \circ \Delta$, where $\Lambda$ is the ring of symmetric functions, $\Delta$ is the usual comultiplication, and $S$ is the usual antipode (which sends the power sum $p_i(y)$ to $-p_i(y)$). The variables in the second factor (i.e. $\{y_j\}$) are the "super variables". Because the map is injective, the supersymmetric functions are abstractly isomorphic to $\Lambda$, but they are convenient for performing certain specialisation tricks. But let me give an example:
$$
\sum_{r=0}^\infty e_r(x/y) t^r = \frac{\prod_i(1+x_i t)}{\prod_j (1+y_j t)}.
$$
These are the elementary super symmetric functions. You can see that if you set a single $x_i = z$ and a single $y_j = z$, then a factor of $1+zt$ cancels in the numerator and denominator, giving an expression independent of $z$ as stated above. Another example is that $p_i(x/y) = p_i(x) - p_i(y)$.
Now we're in a situation where we can explain the identity. Clearly $s_{\mu}(x) = s_\mu(x,1/1)$ by the specialisation property of super symmetric functions. But the right hand side can be obtained by first setting an $x_i$ equal to 1, then passing from $\Lambda$ to super symmetric functions, then setting a $y_j=1$ and all other $y$'s to zero (that setting $x_i=1$ commutes with passing to super symmetric functions isn't totally obvious, but I'm omitting it for the sake of brevity). The upshot of this is that $s_\mu(x,1) = \sum_\lambda s_\lambda(x) s_{\mu/\lambda}(1)$, and $s_{\mu/\lambda}(1)$ is zero unless $\mu/\lambda$ is a horizontal strip. This gives the outer sum in your expression. Then we recognise that second part of the operation replaces $p_i(x)$ with $p_i(x/1) = p_i(x) - p_i(1) = p_i(x)-1$, which recovers your expression, where I have taken for granted that
$$
s_\lambda(x) = \sum_{\alpha \vdash |\lambda|} \frac{\chi_\alpha^\lambda}{z_\alpha}\prod_i p_i^{a_i}.
$$