Is there anything wrong with the following argument?
First of all, by scaling all $x_i$ and all $y_i$ by a positive constant, we may safely assume that $\prod_i x_i = \prod_i y_i =1$.
The result would now follow from the following more general conjecture.
$\bf Conjecture:$ For ${\bf a}=(a_1,\ldots,a_{n-1})\in (\mathbb R_{>0})^{n-1}$ consider
$h_{\bf a}(z) = z^n+a_{n-1}z^{n-1}+\cdots+a_1z+1$. Define the function $f:(\mathbb R_{>0})^{n-1}\to {\mathbb R}$ by
$$
f(a_1,\ldots,a_{n-1}) = \sum_{z\vert h_{\bf a}(z)=0} (\log(-z))^2
$$
where the branch of $\log$ is picked as usual on the complement of ${\mathbb R}_{\leq 0}$ in $\mathbb C$. Then all partial derivatives $\frac{\partial f}{\partial a_k}$ are positive.
First, a few comments. The function $f$ is well-defined because none of the roots
of $h_a(z)$ are positive reals. In addition, $f$ is real-valued, because roots of $h_a(z)$ come in complex conjugate pairs for which values of $\log(-z)$ are complex conjugates.
The consequence of this conjecture is that if ${\bf a}\leq {\bf b}$ coordinate-wise,
then $f({\bf a})\leq f({\bf b})$. Applied to the case when roots of $h_{\bf a}(z)$
and $h_{\bf b}(z)$ are real numbers $-x_i$ and $-y_i$, we get the original conjecture.
Now, let me describe what I think is the proof of this new conjecture.
First of all, as standard, I can write the function $f(\bf a)$ by a contour integral. In a neighborhood of fixed ${\bf a}$ for any large enough $R$ and small enough $\epsilon>0$ there holds
$$
f({\bf a}) = \frac 1{2\pi i}\int_{C} (\log(-z))^2 \frac {h_{\bf a}'(z)}{h_{\bf a}(z)}\,dz
$$
over the contour $C$ on the union of $\mathbb C$ with two copies of $\mathbb R_{>0}$
which I will call "north" and "south" shores (so that $\log(-t)=\log t -\pi i$ on the north shore and $\log(-t)=\log t +\pi i$ on the south shore).
The contour $C$ is the union of the following four pieces $C_\epsilon$, $C_R$, $C_+$, $C_-$.
$C_\epsilon$ is a circle of radius $\epsilon$ around $z=0$ traveled from $\epsilon$-south
to $\epsilon$-north clockwise.
$C_R$ is a circle of radius $R$ around $z=0$ traveled from $R$-north to $R$-south counterclockwise.
$C_+$ is the line segment $[\epsilon,R]$-north.
$C_-$ is the line segment $[R,\epsilon]$-south.
The derivative $\frac{\partial f({\bf a})}{\partial a_k}$ is the integral of the
derivative, so we get:
$$
\frac{\partial f({\bf a})}{\partial a_k}= \frac 1{2\pi i}\int_{C} (\log(-z))^2 \frac \partial{\partial a_k}\frac {h_{\bf a}'(z)}{h_{\bf a}(z)}\,dz
$$
$$
=
\frac 1{2\pi i}\int_{C} (\log(-z))^2 \Big(\frac {z^k}
{h_{\bf a}(z)}
\Big)'\,dz
=
-\frac 1{2\pi i}\int_{C} \Big((\log(-z))^2 \Big)'\frac {z^k}
{h_{\bf a}(z)}
\,dz
$$
$$
=-\frac 1{\pi i}\int_{C} \log(-z)\frac {z^{k-1}}
{h_{\bf a}(z)}\,dz = -\frac 1\pi {\rm Im}\Big(\int_{C} \log(-z)\frac {z^{k-1}}
{h_{\bf a}(z)}\,dz\Big)
$$
We can take a limit as $R\to +\infty$ and $\epsilon\to 0$. Since $k\leq {n-1}$ the integral over $C_R$ goes to zero (the length is $2\pi R$ and the size of the function is
$O(R^{-2}\log R)$). The integral over $C_\epsilon$ also goes to zero, because the $k>=1$,
so the function is $O(\log \epsilon)$ and the length is $2\pi\epsilon$.
So we get
$$
\frac{\partial f({\bf a})}{\partial a_k}
= -\lim_{\epsilon\to 0^+}\frac 1\pi {\rm Im}\Big( \int_{[\epsilon,+\infty]-{\rm north}\cup [+\infty,\epsilon]-{\rm south}}
\log(-z)\frac {z^{k-1}}
{h_{\bf a}(z)}\,dz\Big)
$$
$$
=-\lim_{\epsilon\to 0^+}\frac 1\pi
\int_{\epsilon}^{+\infty}{\rm Im}\Big(
(\log(t)-\pi i)\frac {t^{k-1}}
{h_{\bf a}(t)} -(\log(t)+\pi i)\frac {t^{k-1}}
{h_{\bf a}(t)}\Big)\,dt
$$
$$
=2\lim_{\epsilon \to 0^+}
\int_{\epsilon}^{+\infty}
\frac {t^{k-1}}{h_{\bf a}(t)}\,dt
>0.$$
This finishes the proof of the new conjecture, and consequently of the old conjecture.
Remark:
I am guessing that there is a simpler argument for
$$
\frac{\partial f({\bf a})}{\partial a_k} = 2\int_{0}^{+\infty}
\frac {t^{k-1}}{h_{\bf a}(t)}\,dt
$$
but I am just writing the first thing that came to my mind.
Best Answer
Conditions (1,2,3) are not enough for the inequality to hold. Take $n=3$, $x_1=x_2=3,x_3=0$, $y_1=3,y_2=y_3=0$. Then we need the multiset $(3,3,0,1,1,1)$ (all three $x$'s and 3 times mean of $y$'s) to majorate $(3,0,0,2,2,2)$. But four largest elements of the first multiset have sum $3+3+1+1=8$, while in the second it is $3+2+2+2=9$. So, the claim does not hold in full generality. To be more explicit, we get an opposite inequality for $f(x)=\max(x-2,0)$.
Conditions (1) and
(2') $x_1+...+x_i \geqslant y_1+.....+y_i$ and $x_{i+1}+...+x_n \leqslant y_{i+1}+...+y_n$ for $i=1,....,n-1$
are enough. We prove it by verifying that the multiset of $2n$ numbers $A:=\{x_1,\dots,x_n,y,y,\dots,y\}$ majorates the multiset $B:=\{y_1,\dots,y_n,x,x,\dots,x\}$, where $x=\frac1n \sum x_i$, $y=\frac1n \sum y_i$. Without loss of generality $x\leqslant y$, else change signs of all $x$'s and $y$'s. We have to check that the sum $w_m$ of $m$ largest elements of $B$ does not exceed the sum of $m$ largest elements of $A$. Let $m$ largest elements of $B$ be $y_1,\dots,y_s$ and $(m-s)$ times $x$. Consider two cases.
1-st case. $s\leqslant n-1$. Then $w_m=y_1+\dots+y_s+(m-s)x\leqslant x_1+\dots+x_s+(m-s)y$.
2-nd case. $s=n$. Then $w_m=y_1+\dots+y_n+(m-n)x\leqslant n\cdot y+x_1+\dots+x_{m-n}$.
In both cases we found $m$ elements of $A$ with a sum at least $w_m$, as desired.