As bound variables, $r$ and $i$ must play no role in the expansion and so they shouldn't even appear in our solution. It should be equally evident that $n$ is just along for the ride as a placeholder for the upper limit; we could call it anything, and therefore we may call it nothing and ignore it except in expressions where it will necessarily appear in the output. Notice, too, that $x$ disappears under the summation: it really just serves as a placeholder for arguments of other expressions being summed.
Consequently, expressions like "$\sum_{i=1}^n f(x_i)$" should be understood as operations on $f$ itself. What properties should they have? Only the obvious ones implied by the ring operations $+$ and $\times$ in the algebra of polynomials in the $x_i$. This gives just four simple rules, after we establish the unique roles of the symbols "$x$", "$s$", and "$n$":
ClearAll[sum, x, s, n];
Converting sums of powers of $x$ into the power sum variables $s_i$:
sum[x] := Subscript[s, 1]; sum[x^p_] := Subscript[s, p];
Additivity:
sum[a_ + b_] := sum[a] + sum[b];
Linearity (with respect to "scalars" which do not have any $x_i$ in them):
sum[Times[a_, y__]] /; FreeQ[a, x] := a sum[Times[y]];
The effect of summing constant values (this is the only place $n$ need appear):
sum[a_] /; FreeQ[a, x] := n a;
That should do it. But to perform the algebra, we need to expand algebraic combinations of everything possible into powers of $x$ (and otherwise leave everything else alone):
expand[a_] /; ! FreeQ[a, sum] := Map[Expand[#, x] &, a, Infinity];
expand[a_] := a;
Examples
The question:
sum[x ((sum[x] - x)^2 - sum[x^2])] // expand
$s_1^3-3 s_1 s_2+s_3$
Variance:
sum[(x - sum[x]/n)^2] // expand
$-\frac{s_1^2}{n}+s_2$
Skewness:
sum[(x - sum[x]/n)^3] / sum[(x - sum[x]/n)^2] ^(3/2) // expand
$\frac{\frac{2 s_1^3}{n^2}-\frac{3 s_1 s_2}{n}+s_3}{\left(-\frac{s_1^2}{n}+s_2\right){}^{3/2}}$
etc.
We might worry about misapplications, such as to non-polynomial functions. Not to fear:
sum[Exp[(x + 1)^2]] // expand
$\text{sum}\left[e^{1+2 x+x^2}\right]$
The expansion proceeds insofar as it can, but sum does not know how to go any further, and so stops. sum is also ignorant of other non-polynomial objects, but when they can be converted to polynomials, it works:
sum[Series[Exp[x], {x, 0, 4}] // Normal] // expand
$n+s_1+\frac{s_2}{2}+\frac{s_3}{6}+\frac{s_4}{24}$
(This is a series expansion through order four of $\sum_{i=1}^n e^{x_i}$.)
Edit
Michael E2, in an extension to his answer, reminds us of the value of simplifying multiple sums. His generalization further clarifies the nature of these operations. Because the variable name "$x$" does not explicitly appear in the $s_i$ notation, the variable name is immaterial. What matters are the index names, only insofar as they are used to connect an $x$ with its enclosing summation: as I remarked at the outset, as bound variables they must disappear at the end.
Thus, for example, we could write sum[i^3, i] for $\sum_{i=1}^n x_i^3$: the first argument to sum is a polynomial (or rational function, even) and the second one is a symbol to indicate what is varying over the summation.
A comparable extension is trivial to make: just include the second argument explicitly in the definition of sum. So that we can reproduce the previous solution (where x was the only variable), we can have the second argument default to x if it's missing. Here is the entire (generalized) solution:
ClearAll[sum, s, n];
sum[a_] := sum[a, x];
sum[x_, x_] := Subscript[s, 1]; sum[x_^p_, x_] := Subscript[s, p];
sum[Times[a_, y__], x_] /; FreeQ[a, x] := a sum[Times[y], x];
sum[a_ + b_, x_] := sum[a, x] + sum[b, x];
sum[a_, x_] /; FreeQ[a, x] := n a;
expand[a_] /; ! FreeQ[a, sum] := Map[Expand, a, Infinity]; expand[a_] := a;
For example, compute $\sum_{k=1}^n\sum_{j=1}^n\sum_{i=1}^n (x_i+x_j+x_k - 3\bar{x})^2$, where $\bar{x} = \sum_{i=1}^n x_i / n$ is the average:
sum[sum[sum[(i + j + k - 3 sum[i, i]/n)^3, i], j], k]
$6 s_1^3-9 n s_1 s_2+3 n^2 s_3$
(Notice that it was no problem for i to appear in two summations.) You do have to be a little careful not to go overboard and stretch the use of sum too far:
sum[x^2, x^2] // expand
$s_1$
This surprising result is correct: the symbol x^2 in the second argument is being used to represent the variable summed over; in terms of this pattern, the first argument x^2 is the first power, whence the result is indeed $s_1$.
One can rather easily go further with generalizations, depending on the need. The next big step would be to implement multivariate sums, such as $\sum_{i=1}^n\sum_{j=1}^n x_i y_j^2$, perhaps to be written $s_1(x)s_2(y)$ or--assuming a specific ordering $(x,y)$ of the variables--even as $s_{1,2}$. As we appear to have gone beyond the scope of the present question, I leave these generalizations to interested readers. I hope that the examples given here have shown how stripping the notation down to its essentials can lead to simple and clear solutions.
