When $Z$ is a $\mathcal{N}(0,1)$ random variable, $f$ smooth from $\mathbb{R} \to \mathbb{R}$ we have the Gaussian integration by parts formula $$ \mathbb{E}(Zf(Z)) = \mathbb{E} f'(Z). $$
One analog of the above in the multivariate case is when $Z$ is $\mathcal{N}(0, \Sigma)$, $f: \mathbb{R}^n \to \mathbb{R}$ in which case $$ \mathbb{E}[Z_i f(Z)] = \sum_j \Sigma_{ij} \mathbb{E}\left[ \frac{\partial f}{\partial x_j}(Z)\right]. $$
Now suppose $Z$ is an arbitrary mean 0 variance 1 random variable. In the one-dimensional case we now have the "approximate" integration by parts formula $$ |\mathbb{E}(Zf(Z) - f'(Z))| < C \|f''\|_{\infty} \mathbb{E}(|Z|^3). $$
Is there any known multivariate generalization of this formula? Crucially, in the proof I know of the multivariate gaussian case we use the fact that if the covariance of a Gaussian random vector is diagonal the components are independent to reduce to the 1-D case, but this doesn't seem possible in the general case.
