Let $x^4+c_3x^3+c_2x^2+c_1x+c_0$ be a real polynomial with no real root. Then there are two pairs of conjugate complex roots, $a_1\pm b_1 i$ and $a_2\pm b_2 i$, and one has the identity
$$ c_1^2-c_1c_2c_3+c_0c_3^2=-4a_1a_2\bigg((a_1+a_2)^2+(b_1+b_2)^2\bigg) \bigg((a_1+a_2)^2+(b_1-b_2)^2\bigg) \tag{1} $$
(1) can be easily confirmed by a CAS, and even the natural, "greedy" method consisting in replacing the $c_k$ by their values in terms of $a_k$ and $b_k$ in the LHS and expanding, can be carried out in a not too unreasonable way. But I wonder if the computation can be simplified by some extra insight ?
The pairwise sums of the roots are $2a_1, 2a_2,$ and $(a_1+a_2)\pm(b_1 \pm b_2)i$, so the right-hand side is the product of these six factors and an extra minus sign. So your expression could be written as $-\prod_{i<j}(r_i + r_j)$.
So, since the right-hand side is symmetric in the four roots, it's automatic it can be written as a polynomial expression in the $c_i$'s. This much is hopefully clear.
In fact, though this is not exactly obvious at first, it is the Schur polynomial $s_{3,2,1}(a_1,a_2,a_3,a_4)$. More generally, if the roots are $r_1, \ldots, r_n$, the product of the $r_i+r_j$ is the Schur polynomial $s_{n-1, n-2, \ldots, 1}(x_1, \ldots, x_n)$. This follows from the formula
$s_\lambda = \frac{\det(r_i^{\lambda_j+j-1})}{\det(r_i^{j-1})}$ (see the Schur Polynomial Wikipedia page for this formula and the one below.)
The denominator is the Vandermonde determinant and the numerator is just like a Vandermonde determinant, but the exponent in the $j$-th column of the matrix is $j-1 + \lambda_j$ instead of just $j-1$.
In this case $\lambda = ((n-1),(n-2),\ldots,1)$ so the exponents are $2(j-1)$, so the numerator is actually just another Vandermonde determinant, instead using the squares of the roots, $r_i^2$. So we have
$s_\lambda = \frac{\prod_{i<j}(r_i^2 - r_j^2)}{\prod_{i<j}(r_i-r_j)} = \frac{\prod_{i<j}(r_i+r_j)(r_i - r_j)}{\prod_{i<j}(r_i-r_j)} = \prod_{i<j}(r_i+r_j)$.
You can convert this to an expression in the $c_i$'s, that is, the elementary symmetric polynomials, via the Jacobi-Trudi formula (see the Wiki page)
$s_\lambda = \det(e_{\lambda_i' + j - i}).$
Here $e_k = c_{n-k}$ (the usual indexing goes the other way). Anyway, in the case $n=4$, this gives
$s_\lambda = \det \begin{pmatrix}e_3 & e_4 & 0 \\ e_1 & e_2 & e_3 \\ 0 & e_0 & e_1 \end{pmatrix} = e_3e_2e_1 - e_4 e_1^2 - e_3^2e_0 = c_1c_2c_3 - c_0c_3^2 - c_1^2.$
This matches with your formula (with the extra minus sign).