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?
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}&amp;\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&amp;=\binom{n+1}{2}=\frac{1}{2}n(n+1)\tag{1}\\
\sum_{i=0}^n\sum_{j=0}^{n-i}i(i+j)&amp;=\binom{n+3}{4}+2\binom{n+2}{4}\\
&amp;=\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)&amp;=\binom{n+5}{6}+8\binom{n+4}{6}+6\binom{n+3}{6}\\
&amp;=\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&gt;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}&amp;=\frac{x}{(1-x)^2}
=\sum_{n=0}^{\infty}nx^n\\
\left(x\frac{x}{dx}\right)^2\frac{1}{1-x}&amp;=\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}&amp;=\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&amp;=\sum_{i=0}^{n}i\cdot1\\
&amp;=\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}\\
&amp;=\sum_{i=0}^n[x^i]\frac{x}{(1-x)^2}[x^{n-i}]\frac{1}{1-x}\\
&amp;=[x^n]\frac{x}{(1-x)^3}\tag{5}\\
&amp;=[x^n]x\sum_{k=0}^{\infty}\binom{-3}{k}(-x)^k\\
&amp;=[x^{n-1}]\sum_{k=0}^{\infty}\binom{k+2}{2}x^k\tag{6}\\
&amp;=\binom{n+1}{2}\\
&amp;=\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}&amp;\sum_{j=0}^{n-i}i(i+j)\\
&amp;=\sum_{i=0}^ni^2\sum_{j=0}^{n-i}1+\sum_{i=0}^ni\sum_{j=0}^{n-i}j\\
&amp;= \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}\\
&amp;\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}\\
&amp;= \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}\\
&amp;=[x^n]\frac{x+2x^2}{(1-x)^5}\\
&amp;=\left([x^{n-1}]+2[x^{n-2}]\right)\sum_{k=0}^{\infty}\binom{-5}{k}(-x)^k\\
&amp;=\left([x^{n-1}]+2[x^{n-2}]\right)\sum_{k=0}^{\infty}\binom{k+4}{4}x^k\\
&amp;=\binom{n+3}{4}+2\binom{n+2}{4}\\
&amp;=\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}&amp;\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}i(i+j)(i+j+k)\\
&amp;=\sum_{i=0}^{n}\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}(i^3+2i^2j+i^2k+ij^2+ijk)\\
&amp;=\sum_{i=0}^{n}\sum_{j=0}^{n-i}\sum_{k=0}^{n-i-j}(i^3+4i^2j+ijk)\tag{7}\\
&amp;=\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}\\
&amp;\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&amp;=[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&amp;=[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&amp;=[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)
&amp;=6\binom{n+3}{6}+8\binom{n+4}{6}+\binom{n+5}{6}\\
&amp;=\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&amp;=\color{blue}{1}\binom{n+1}{2}\\
\sum_{i=0}^n\sum_{j=0}^{n-i}i(i+j)&amp;=\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)&amp;=\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&amp;\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)\\
&amp;=
\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}&amp;\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)\\
&amp;=\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&amp;-\sum_{i=1}^ni=0\\
\sum_{i=0}^n\sum_{j=0}^{n-i}i(i+j)&amp;-\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)&amp;-
\sum_{i=1}^n\sum_{j=1}^{n-i}\sum_{k=1}^{n-i-j}i(i+j)(i+j+k)\\
&amp;\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)&amp;=\sum_{n=0}^{\infty}na_nx^n\\
  \left(x\frac{d}{dx}\right)^2A(x)&amp;=\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 &amp;\qquad A(x)=&amp;\sum_{n=0}^{\infty}a_nx^n\\
  \left(\sum_{k=0}^{n}a_k\right)_{n\geq 0}\qquad&amp;\qquad\frac{1}{1-x}A(x)=&amp;\sum_{n=0}^{\infty}\left(\sum_{k=0}^{n}a_k\right)x^n
  \end{array}
This is quite easy if you keep one term constant then its simply $$\sum c.\sum \frac{(n)(n+1)}{2}-(\sum c\sum k)$$  where $c$ is the constant term and $n=c+1$ and $k=1..c$ for you example $c=1..4,n=1..5\ \text{and}\ k=i$ and if its random I hardly think there exists a method.