The exact formula for $p_\Lambda$ and $p_S$ is given in Theorem 3.3 in D. Panyushev and O. Yakimova "The PRV-formula for tensor product decompositions and its applications", Funct. Anal. Appl., 42 (2008), 45-52. Setting $a$ the number of fundamental weights that appear with a nonzero coefficient in the highest weight of $R$, and $b$ the number of those that are moreover self-dual, we have:

$$p_\Lambda = (a + \varepsilon b)/2 \quad\text{and} \quad p_S = (a - \varepsilon b)/2,$$

where $\varepsilon$ denotes the Frobenius-Schur indicator of $R$, which is by definition $1$ if $R$ is real, $-1$ if $R$ is quaternionic.

From there, it is easy to deduce that the answer is "yes". Indeed the weak inequalities between $p_\Lambda$ and $p_S$ (depending on the sign of $\varepsilon$) follow immediately from the formula. Moreover the equality case can only occur for real representations, essentially because the sum of a fundamental weight with its dual is always of real type.

Bruno and I actually ended up writing down our own proof of these results, unaware that it had already been done :-/. Thanks to the anonymous referee for pointing out the proper reference to us!