Let $f(x) = |x|^{p-1}x$ for some $2 \leq p \leq 3$. I've seen the pointwise estimate
$$\left| f(x+y+z)-[f(y)+f(z)+f'(y)x + f'(y)z] \right| \lesssim |f(x)| + |f'(x)y| + |f'(z)x| + |f'(z)y|$$
in Lemma 3.12 from this paper https://arxiv.org/abs/1503.05024.
By homogeneity considerations, we can reduce to the case $x^2 + y^2 + z^2 = 1$ and then it suffices to check the inequality when the RHS is zero, i.e. for $(x,y,z)$ in a neighborhood of $(0,\pm 1,0)$ or $(0,0,\pm 1)$. Then, the author states that the inequality follows from the Taylor expansion of $f(x)$ but I cannot see how it is applied.
I think I have found a proof by myself. It's done by comparing the order at which each side vanishes. Assume that $2 \leq p \leq 3$.
Case 1 : Neighborhood of $(0,\pm 1, 0)$
For the LHS, we use a Taylor approximation around $y$ : $$LHS = |f(x+y+z) - f(y) - f'(y)(x+z) - f(z)| \lesssim |x+z|^2 + |z|^{p} \lesssim x^2 + z^2$$
where we used that $p \geq 2$.
The RHS behaves like : $$RHS \sim |x|^{p} + |x|^{p-1} + |z|^{p-1}|x| + |z|^{p-1} \gtrsim |x|^{p-1} + |z|^{p-1}$$
so we need $2 \geq p-1$.
Case 2 : Neighborhood of $(0,0,\pm1)$
For the LHS, we use a Taylor approximation around $z$ :
$$LHS = |f(x+y+z) - f(z) - f(y) - f'(y)(x+z)| \lesssim |x+y| + |y|^{p} + |y|^{p-1} \lesssim |x| + |y|$$
where we used that $p-1 \geq 1$.
The RHS behaves like : $$RHS \sim |x|^p + |x|^{p-1}|y| + |x| + |y| \gtrsim |x| + |y|$$