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];

1. Converting sums of powers of $x$ into the power sum variables $s_i$:

        sum[x] := Subscript[s, 1]; sum[x^p_] := Subscript[s, p];

2. Additivity:

        sum[a_ + b_] := sum[a] + sum[b];

3. 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]];

4. 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}$.)