Q:

Calculate the sum of first $n$ natural numbers taken $k$ at a time

So sum of first $n$ natural numbers taken $1$ at a time is $n\cdot(n+1)/2$
but what about $2,3,\dots,k$ at a time?
Is there a general formula?
For example, taking 1 at a time
$$\sum_{i = 1}^{n} i = \frac{n(n+1)}{2}$$
taking 2 at a time
$$\sum_{i = 1}^{n}\sum_{j = 1}^{n - i} i\cdot(i+j) = \frac{1}{24} n(n+1)(n-1)(3n+2)$$
taking 3 at a time
$$\sum_{i = 1}^{n}\sum_{j = 1}^{n - i}\sum_{h = 1}^{n - j - i} i\cdot(i+j)\cdot(i + j +h) = \frac{1}{48} n^2(n+1)^2(n-1)(n-2)$$
taking $k$ at a time?

A:

Overview: This answer consists of three parts

At first we develop summation formulae for slightly different $k$-fold sums for $k=1,2,3$. We take indices starting from $0$ instead of $1$, which makes calculation somewhat more convenient.
We obtain for general $k\geq 1$ a nice expression of the $k$-fold sum
  \begin{align*}
\sum_{i_1=0}^{n}\sum_{i_2=0}^{n-i_1}&\ldots\sum_{i_k=0}^{n-i_1-i_2-\ldots-i_{k-1}}i_1(i_1+i_2)\cdots(i_1+i_2+\ldots+i_k)\\
\end{align*}
  as sum of binomial coefficients .
We look at the difference of the $k$-fold sums when starting from $0$ and when starting from $1$.

Hint: If you are not familiar with generating functions you might want to look at the backstage info at the end of this answer.

Part 1: Summation formula for $k=1,2,3$ and start index $0$
The following is valid
\begin{align*}
\sum_{i=0}^ni&=\binom{n+1}{2}=\frac{1}{2}n(n+1)\tag{1}\\
\sum_{i=0}^n\sum_{j=0}^{n-i}i(i+j)&=\binom{n+3}{4}+2\binom{n+2}{4}\\
&=\frac{1}{24}n(n+1)(n+2)(3n+1)\tag{2}\\
\sum_{i=0}^n\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}i(i+j)(i+j+k)&=\binom{n+5}{6}+8\binom{n+4}{6}+6\binom{n+3}{6}\\
&=\frac{1}{48}n^2(n+1)^2(n+2)(n+3)\tag{3}\\
\end{align*}

Note, the representation of binomial coefficients looks promising, cause a pattern for $k>3$ can easily be derived. We will see in Part 2, that the multiplicities of the binomial coefficients also follow a well-known pattern.
In order to show these formulae we need a repertoire of $\left(x\frac{d}{dx}\right)^k\frac{1}{1-x}$ for $k=1,\ldots,3$.
We obtain
\begin{align*}
\left(x\frac{x}{dx}\right)\frac{1}{1-x}&=\frac{x}{(1-x)^2}
=\sum_{n=0}^{\infty}nx^n\\
\left(x\frac{x}{dx}\right)^2\frac{1}{1-x}&=\frac{x(1+x)}{(1-x)^3}
=\sum_{n=0}^{\infty}n^2x^n\\
\left(x\frac{x}{dx}\right)^3\frac{1}{1-x}&=\frac{x(1+4x+x^2)}{(1-x)^4}
=\sum_{n=0}^{\infty}n^3x^n\\
\end{align*}
It's also convenient to use the coefficient of operator $[x^n]$ to denote the coefficient of $x^n$ in a generating series.

We obtain for $k=1$
  \begin{align*}
\sum_{i=0}^{n}i&=\sum_{i=0}^{n}i\cdot1\\
&=\sum_{i=0}^n[x^i]\left(x\frac{x}{dx}\right)\frac{1}{1-x}[x^{n-i}]\frac{1}{1-x}\tag{4}\\
&=\sum_{i=0}^n[x^i]\frac{x}{(1-x)^2}[x^{n-i}]\frac{1}{1-x}\\
&=[x^n]\frac{x}{(1-x)^3}\tag{5}\\
&=[x^n]x\sum_{k=0}^{\infty}\binom{-3}{k}(-x)^k\\
&=[x^{n-1}]\sum_{k=0}^{\infty}\binom{k+2}{2}x^k\tag{6}\\
&=\binom{n+1}{2}\\
&=\frac{1}{2}n(n+1)
\end{align*}
  and (1) follows.

Comment:

In (4) we use $$i=[x^i]\sum_{k=0}^{\infty}kx^x=[x^i]\left(x\frac{x}{dx}\right)\frac{1}{1-x}$$ and $$1=[x^{n-i}]\sum_{k=0}^{\infty}x^k=[x^{n-i}]\frac{1}{1-x}$$
In (5) we use $[x^n]A(x)B(x)=\sum_{i=0}^{n}\left([x^i]A(x)\right)\left([x^{n-i}]B(x)\right)$
In (6) we use the formula $\binom{-n}{k}=\binom{n+k-1}{k}(-1)^k=\binom{n+k-1}{n-1}(-1)^k$

We obtain for $k=2$
  \begin{align*}
\sum_{i=0}^{n}&\sum_{j=0}^{n-i}i(i+j)\\
&=\sum_{i=0}^ni^2\sum_{j=0}^{n-i}1+\sum_{i=0}^ni\sum_{j=0}^{n-i}j\\
&= \sum_{i=0}^n[x^i]\left(x\frac{x}{dx}\right)^2\frac{1}{1-x}[x^{n-i}]\frac{1}{(1-x)^2}\\
&\qquad+\sum_{i=0}^n[x^i]\left(x\frac{x}{dx}\right)\frac{1}{1-x}[x^{n-i}]\frac{1}{1-x}\left(x\frac{x}{dx}\right)\frac{1}{1-x}\\
&= \sum_{i=0}^n[x^i]\frac{x(1+x)}{(1-x)^3}[x^{n-i}]\frac{1}{(1-x)^2}+\sum_{i=0}^n[x^i]\frac{x}{(1-x)^2}[x^{n-i}]\frac{x}{(1-x)^3}\\
&=[x^n]\frac{x+2x^2}{(1-x)^5}\\
&=\left([x^{n-1}]+2[x^{n-2}]\right)\sum_{k=0}^{\infty}\binom{-5}{k}(-x)^k\\
&=\left([x^{n-1}]+2[x^{n-2}]\right)\sum_{k=0}^{\infty}\binom{k+4}{4}x^k\\
&=\binom{n+3}{4}+2\binom{n+2}{4}\\
&=\frac{1}{24}n(n+1)(n+2)(3n+1)
\end{align*}
  and (2) follows.

$$ $$

We obtain for $k=3$
  \begin{align*}
\sum_{i=0}^{n}&\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}i(i+j)(i+j+k)\\
&=\sum_{i=0}^{n}\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}(i^3+2i^2j+i^2k+ij^2+ijk)\\
&=\sum_{i=0}^{n}\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}(i^3+4i^2j+ijk)\tag{7}\\
&=\sum_{i=0}^{n}i^3\sum_{j=0}^{n-i}1\sum_{k=0}^{n-i-j}1
+4\sum_{i=0}^{n}i^2\sum_{j=0}^{n-i}j\sum_{k=0}^{n-i-j}1\tag{8}\\
&\qquad+\sum_{i=0}^{n}i\sum_{j=0}^{n-i}j\sum_{k=0}^{n-i-j}k\\
\end{align*}

Note, in (7) we use the symmetry 
\begin{align*}
\sum_{i=0}^{n}\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}i^2j
=\sum_{i=0}^{n}\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}i^2k
=\sum_{i=0}^{n}\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}ij^2
\end{align*}
The calculation of the three sums in (8) is straight forward and can be done similarly to $k=1,2$.

We obtain
\begin{align*}
\sum_{i=0}^{n}i^3\sum_{j=0}^{n-i}1\sum_{k=0}^{n-i-j}1&=[x^n]\frac{x(1+4x+x^2)}{(1-x)^7}
=\binom{n+3}{6}+4\binom{n+4}{6}+\binom{n+5}{6}\\
\sum_{i=0}^{n}i^2\sum_{j=0}^{n-i}j\sum_{k=0}^{n-i-j}1&=[x^n]\frac{x^2(1+x)}{(1-x)^7}
=\binom{n+4}{6}+\binom{n+3}{6}\\
\sum_{i=0}^{n}i\sum_{j=0}^{n-i}j\sum_{k=0}^{n-i-j}k&=[x^n]\frac{x^3}{(1-x)^7}
=\binom{n+3}{6}\\
\end{align*}
Combining  the three sums according to (7) results in
\begin{align*}
\sum_{i=0}^{n}\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}(i^3+4i^2j+ijk)
&=6\binom{n+3}{6}+8\binom{n+4}{6}+\binom{n+5}{6}\\
&=\frac{1}{48}n^2(n+1)^2(n+2)(n+3)
\end{align*}
  and the claim (3) follows.

$$ $$

Part 2: Summation formula for all $k\geq 1$ and startindex $0$.
In the following we don't give a proof for general $k$, but we provide some aspects which give strong evidence for the correctness of the claim.
When looking at 
\begin{align*}
\sum_{i=0}^ni&=\color{blue}{1}\binom{n+1}{2}\\
\sum_{i=0}^n\sum_{j=0}^{n-i}i(i+j)&=\color{blue}{1}\binom{n+3}{4}+\color{blue}{2}\binom{n+2}{4}\\
\sum_{i=0}^n\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}i(i+j)(i+j+k)&=\color{blue}{1}\binom{n+5}{6}+\color{blue}{8}\binom{n+4}{6}+\color{blue}{6}\binom{n+3}{6}\\
\end{align*}
  the shape of the binomial coefficients can be easily generalized. The coefficients 
  \begin{align*}
\color{blue}{1};\quad \color{blue}{1},\color{blue}{2};\quad \color{blue}{1},\color{blue}{8},\color{blue}{6}
\end{align*}
  are part of the OEIS sequence A008517 and give the values of the Second order Eulerian Triangle $T(k,l),1\leq l\leq k$.

The values $T(4,l), 1\leq l\leq 4$ are $\color{blue}{1},\color{blue}{22},\color{blue}{58},\color{blue}{24}$ and indeed, it is easy to verify that
\begin{align*}
\sum_{i=0}^n&\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}\sum_{l=0}^{n-i-j-l}i(i+j)(i+j+k)(i+j+k+l)\\
&=
\color{blue}{1}\binom{n+7}{8}+\color{blue}{22}\binom{n+6}{8}+\color{blue}{58}\binom{n+5}{8}+\color{blue}{24}\binom{n+4}{8}
\end{align*}

We are now in the position to
Claim: The following formula of the $k$-fold sum is valid
  \begin{align*}
\sum_{i_1=0}^{n}\sum_{i_2=0}^{n-i_1}&\ldots\sum_{i_k=0}^{n-i_1-i_2-\ldots-i_{k-1}}i_1(i_1+i_2)\cdots(i_1+i_2+\ldots+i_k)\\
&=\sum_{l=1}^{k}\color{blue}{T(k,l)}\binom{n+2k-l}{2k}\tag{9}
\end{align*}
  with $T(k,l)$ the numbers of the second order Eulerian Triangle (A008517).

$$ $$

Part 3: Summation formula for all $k\geq 1$ and startindex $1$.
Finally some aspects about the summation formula with indices starting with $1$.
When looking at the differences of
  \begin{align*}
\sum_{i=0}^ni&-\sum_{i=1}^ni=0\\
\sum_{i=0}^n\sum_{j=0}^{n-i}i(i+j)&-\sum_{i=1}^n\sum_{j=1}^{n-i}i(i+j)=\frac{1}{6}n(n+1)(2n+1)\\
\sum_{i=0}^n\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}i(i+j)(i+j+k)&-
\sum_{i=1}^n\sum_{j=1}^{n-i}\sum_{k=1}^{n-i-j}i(i+j)(i+j+k)\\
&\qquad=\frac{1}{12}(n+1)^2(n+2)^2(2n+3)\\
\end{align*}
we can find them as A000330 ($k=2)$ and as A108674 ($k=3)$ in the OEIS database.
Some further elaboration could be to look for a difference formula for general $k$ and combine it with the result of (9) to obtain a nice expression for OPs $k$-fold sum with indices starting from $1$.

Backstage info:
You might want to skip this section, if generating functions are already known.
When looking at the $2$-fold sum
\begin{align*}
\sum_{i=1}^n\sum_{j=1}^{n-i}i(i+j)=\sum_{i=1}^n\left(i^2\sum_{j=1}^{n-i}1\right)+\sum_{i=1}^n\left(i\sum_{j=1}^{n-i}j\right)
\end{align*}
the outer sum has the shape of a Cauchy product
\begin{align*}
\sum_{i=0}^{n}a_ib_{n-i}
\end{align*}
with index $i$ starting at $1$ instead of $0$.
Since we want to use generating functions $A(x)=\sum_{n= 0}^{\infty}a_nx^n$ to derive the summation formulae, and the product of generating function 
\begin{align*}
A(x)B(x)=\sum_{k= 0}^{\infty}a_kx^k\sum_{l= 0}^{\infty}b_lx^l
= \sum_{n= 0}^{\infty}\left(\sum_{k=0}^{n-k}a_kb_{n-k}\right)x^n
\end{align*}
 gives Cauchy products, we consider $k$-fold sums with index starting from $0$ instead.
This is interesting by its own and later we can look at the difference to $k$-fold sums with indices starting from $1$.
We can successively apply the $\left(x\frac{d}{dx}\right)$-operator to a generating function $A(x)$ to obtain
\begin{align*}
  \left(x\frac{d}{dx}\right)A(x)&=\sum_{n=0}^{\infty}na_nx^n\\
  \left(x\frac{d}{dx}\right)^2A(x)&=\sum_{n=0}^{\infty}n^2a_nx^n
  \end{align*}
Multiplication of $A(x)$ with $\frac{1}{1-x}$ results in summing up the coefficients $a_n$
\begin{array}{crl}
  (a_n)_{n\geq 0}\qquad &\qquad A(x)=&\sum_{n=0}^{\infty}a_nx^n\\
  \left(\sum_{k=0}^{n}a_k\right)_{n\geq 0}\qquad&\qquad\frac{1}{1-x}A(x)=&\sum_{n=0}^{\infty}\left(\sum_{k=0}^{n}a_k\right)x^n
  \end{array}

