Solution 1:

Here's a proof using generating functions; I haven't given much thought yet to how it could be translated into a combinatorial argument.

The generating function for the number of partitions is

$$p(x)=(1+x+x^2+\dotso)(1+x^2+x^4+\dotso)(1+x^3+x^6+\dotso)\dots=\prod_m\frac1{1-x^m}\;.$$

To count the number of times a part $k$ occurs in the partitions of $n$ (the left-hand side of the equation), we can replace the $k$-th factor by one including this count:

$$ \begin{align} f_k(x) &=\left(0+1x^k+2(x^k)^2+3(x^k)^3+\dotso\right)\prod_{m\neq k}\frac1{1-x^m}\\ &=\left(x^k\frac{\mathrm d}{\mathrm d(x^k)}\left(1+x^k+(x^k)^2+(x^k)^3+\dotso\right)\right)\prod_{m\neq k}\frac1{1-x^m}\\ &=\left(x^k\frac{\mathrm d}{\mathrm d(x^k)}\frac1{1-x^k}\right)\prod_{m\neq k}\frac1{1-x^m}\\ &=\frac{x^k}{(1-x^k)^2}\prod_{m\neq k}\frac1{1-x^m}\\ &=\frac{x^k}{1-x^k}\prod_m\frac1{1-x^m}\\ &=p(x)\frac{x^k}{1-x^k}\;.\\ \end{align} $$

To count the number of at-least-$k$-fold occurrences of parts (the right-hand side of the equation), consider this generation function (which I'll first write out for $k=2$ to illustrate the idea and then generalize):

$$\begin{align} g_2(x,y) &=\left(1+x+y(x^2+x^3+\dotso)\right)\left(1+x^2+y((x^2)^2+(x^2)^3+\dotso)\right)\dots\\ &=\left(1+x+\frac{yx^2}{1-x}\right)\left(1+x^2+\frac{y(x^2)^2}{1-x^2}\right)\dots\\ &=\frac{1-x^2+yx^2}{1-x}\frac{1-(x^2)^2+y(x^2)^2}{1-x^2}\dots\\ &=p(x)\left(1-x^2+yx^2\right)\left(1-(x^2)^2+y(x^2)^2\right)\dots\;,\\ g_k(x,y) &=\prod_m\left(\sum_{l=0}^{k-1}(x^m)^l+y\sum_{l=k}^\infty(x^m)^l\right)\\ &=p(x)\prod_m\left(1-(x^k)^m+y(x^k)^m\right)\;. \end{align} $$

Every factor of $y$ tracks the at-least-$k$-fold use of a part. What we want is the total count of factors of $y$ included in the coefficient of $x^n$; that is, we want to count $j$ times the coefficient of $y^jx^n$. This we can get by differentiating with respect to $y$ and then setting $y$ to $1$; the coefficient of $x^n$ in the result will be the desired count. But

$$\begin{align} \left.\frac{\mathrm d}{\mathrm dy}g_k(x,y)\right|_{y=1} &=\left.\frac{\mathrm d}{\mathrm dy}p(x)\prod_m\left(1-(x^k)^m+y(x^k)^m\right)\right|_{y=1}\\ &=\left.p(x)\sum_{l=1}^\infty(x^k)^l\prod_{m\neq l}\left(1-(x^k)^m+y(x^k)^m\right)\right|_{y=1}\\ &=p(x)\sum_{l=1}^\infty(x^k)^l\\ &=p(x)\frac{x^k}{1-x^k}\;,\\ \end{align}$$

which coincides with the result for the left-hand side.

Solution 2:

You can establish a bijection in the partitions of $n$ between parts equal to $k$ and equal parts occurring at least $k$ times by ‘rotating rectangles’ in the diagrams of the partitions.

If $f_\lambda(i)\geqslant k$, we define $\rho_i(\lambda)$ to be the partition resulting from removing $k$ parts equal to $i$ and replacing them with $i$ parts equal to $k$ (i.e. ‘rotating’ an $i\times k$ rectangle in the diagram of $\lambda$, as shown below for $\rho_5$ with $k=3$). We could write $\rho_i(\lambda)=\lambda\:-<\!\!i^k\!\!>+<\!\!k^i\!\!>$.

$\qquad\qquad\qquad\qquad\qquad\qquad$ $\rho_5$

$\rho_i(\lambda)$ is a partition that contains at least $i$ elements equal to $k$. Indeed, if $\mu$ is a partition such that $f_\mu(k)=j$ then there are $j$ distinct partitions, $\lambda_1,\dots,\lambda_j$, with $\rho_i(\lambda_i)=\mu$.