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.