Another proof of the identity for the Hilbert transform can be found by means of Fourier transform. If $\sigma(\xi) = i \operatorname{sign} \xi$, then
$$ \mathcal{F}(H(f g) - f H g - g H f)(\xi) = \int \mathcal{F} f(\eta) \mathcal{F} g(\xi - \eta) (\sigma(\xi) - \sigma(\xi - \eta) - \sigma(\eta)) d\eta . $$
It turns out that
$$ \sigma(\xi) - \sigma(\xi - \eta) - \sigma(\eta) = \sigma(\xi)\sigma(\xi - \eta) \sigma(\eta) , $$
and so $H(f g) - f H g - g H f = H(Hf Hg)$.
The same calculation for the Riesz transform, that is, $\sigma(\xi) = \xi_j / |\xi|$, leads to something much more complicated: $\sigma(\xi) - \sigma(\xi - \eta) - \sigma(\eta)$ is not a product of three functions as it was for the Hilbert transform. Therefore, $R_j(f g) - f R_j g - g R_j f$ is not of the form $A(Bf Bg)$ for whatever Fourier multipliers $A$ and $B$.
However, the expression for $\sigma(\xi) - \sigma(\xi - \eta) - \sigma(\eta)$ can be be slightly simplified: if I am not mistaken,
$$ \sigma(\xi) - \sigma(\xi - \eta) - \sigma(\eta) = \sigma(\xi) (1 - (2 + 2 |\xi - \eta|^{-1} |\eta|^{-1} \langle \xi - \eta, \eta\rangle)^{1/2}).$$
It follows that indeed $R_j(f g) - f R_j g - g R_j f$ is the Riesz transform of something rather regular; namely, it is equal to $R_j(B(f,g))$, where $$B(f, g)(z) = \text{p.v.}\iint f(x) g(y) k(z - x, z - y) dx dy$$ with singular kernel $k$ given by $$\mathcal{F} k(\xi, \eta) = (1 - (2 + 2 |\xi|^{-1} |\eta|^{-1} \langle \xi, \eta\rangle)^{1/2}).$$
