I claim first that
$$\binom{n}kk^r=\sum_{i=\max\{0,k-r\}}^k\binom{n}{i,k-i,n-k}{r\brace{k-i}}(k-i)!\;.\tag{1}$$
The lefthand side of $(1)$ clearly counts the ways to choose $K\subseteq[n]$ such that $|K|=k$ and then choose a function from $[r]$ to $K$. The $i$ term on the righthand side of $(1)$ counts the ways to choose a $k$-element subset $K$ of $[n]$, choose an $i$-element subset $I$ of $K$, and then choose a function from $[r]$ onto $K\setminus I$. Clearly this is possible if and only if $\max\{0,k-r\}\le i\le k$, so the two sides of $(1)$ count the same thing.
Now just rearrange the righthand side of the desired identity, expanding $(1+x)^{n-j}$ by the binomial theorem, making a change of index ($k=i+j$), and reversing the order of summation:
$$\begin{align*}
\sum_{j=0}^r\binom{n}j{r\brace j}j!(1+x)^{n-j}x^j&=\sum_{j=0}^r\binom{n}j{r\brace j}j!x^j\sum_{i=0}^{n-j}\binom{n-j}ix^i\\\\
&=\sum_{i=0}^n\sum_{j=0}^{\min\{r,n-i\}}\binom{n}j\binom{n-j}i{r\brace j}j!x^{i+j}\\\\
&=\sum_{i=0}^n\sum_{k=i}^{\min\{i+r,n\}}\binom{n}{k-i}\binom{n-k+i}i{r\brace{k-i}}(k-i)!x^k\\\\
&=\sum_{k=0}^n\sum_{i=\max\{0,k-r\}}^k\binom{n}{i,k-i,n-k}{r\brace{k-i}}(k-i)!x^k\;.
\end{align*}$$
Finally, apply $(1)$.
Here is an approach based upon generating functions. The following can be found in section II.3.1 in Analytic combinatorics by P. Flajolet and R. Sedgewick:
The class $S^{(A,B)}$ of set partitions with block sizes in $A\subseteq \mathbb{Z}_{\geq 1}$ and with a number of blocks that belongs to $B$ has exponential generating function
\begin{align*}
S^{(A,B)}(z)=\beta(\alpha(z))\qquad\text{where}\qquad \alpha(z)=\sum_{a\in A}\frac{z^a}{a!},\quad \beta(z)=\sum_{b\in B}\frac{z^b}{b!}\tag{1}
\end{align*}
We decompose the problem into two parts and use (1) to derive generating functions for each part.
First part: We consider $p$ partitions with size $r$.
\begin{align*}
A_1=\{r\}, B_1=\{p\}\qquad\text{where}\qquad\alpha_1(z)=\frac{z^r}{r!},\beta_1(z)=\frac{z^p}{p!}
\end{align*}
We obtain a generating function
\begin{align*}
\beta_1\left(\alpha_1\left(z\right)\right)=\frac{1}{p!}\left(\frac{z^r}{r!}\right)^p=\frac{z^{rp}}{p!\left(r!\right)^p}\tag{2}
\end{align*}
Second part: We consider $k-p$ partitions with size $\geq 1$.
\begin{align*}
A_2={\mathbb{Z}}_{\geq 1}, B_2=\{k-p\}\qquad\text{where}\qquad
\alpha_2(z)=\sum_{n=1}^\infty \frac{z^n}{n!},\beta_2(z)=\frac{z^{k-p}}{(k-p)!}
\end{align*}
We obtain a generating function
\begin{align*}
\beta_2\left(\alpha_2\left(z\right)\right)&=\frac{1}{(k-p)!}\left(e^z-1\right)^{k-p}
=\sum_{n=k-p}^\infty {n\brace k-p}\frac{z^n}{n!}\tag{3}
\end{align*}
Note the coefficients of (3) are the Stirling numbers of the second kind.
Since the original problem is the cross product of the two parts, the resulting generating functions is the product of the generating functions in (2) and (3). We denote with
\begin{align*}
a_{r,p}(n,k)
\end{align*}
the number of distributions of $n$ elements into $k$ partitions with at least $p$ partitions having exactly $r$ elements.
We obtain for $1\leq p\leq k\leq n, 1\leq r\leq \left\lfloor\frac{n}{r}\right\rfloor$
\begin{align*}
\beta_1\left(\alpha_1\left(z\right)\right)\beta_2\left(\alpha_2\left(z\right)\right)
&=\sum_{n=k-p}^\infty {n\brace k-p}\frac{z^n}{n!}\cdot\frac{z^{rp}}{p!\left(r!\right)^p}\\
&=\frac{1}{p!\left(r!\right)^p}\sum_{n=k-p+rp}{n-rp\brace k-p}\frac{z^n}{(n-rp)!}\\
&=\frac{(rp)!}{p!(r!)^p}\binom{n}{rp}\sum_{n=k+(r-1)p}{n-rp\brace k-p}\frac{z^n}{n!}\\
&=\sum_{n=k+(r-1)p}a_{r,p}(n,k)
\end{align*}
We conclude:
\begin{align*}
\color{blue}{a_{r,p}(n,k)=\frac{(rp)!}{p!(r!)^p}\binom{n}{rp}{n-rp\brace k-p}}\tag{4}
\end{align*}
Example: Calculating OPs example with $n=4,k=2,p=1$ and $r=1$ using (4) results in
\begin{align*}
a_{1,1}(4,2)=\binom{4}{1}{3\brace 1}=4
\end{align*}
in accordance with OPs calculation.
Best Answer
The answer will be a recurrence together with an explanation justifying that recurrence. For example, the binomial coefficients satisfy the two-variable recurrence $$\binom{k}n=\binom{k-1}{n-1}+\binom{k-1}n\tag{1}$$ (which you may recognize as the rule used to form Pascal’s triangle). The Stirling numbers of the second kind, which are also written $$\left\{\matrix{k\\n}\right\}$$ instead of $S(k,n)$, satisfy a two variable recurrence quite similar to $(1)$, though a little more complicated. One term counts the partitions of $[k]$ in which $\{k\}$ is a block by itself; the other counts the partitions of $[k]$ in which $k$ is in a block with something else.
To get you started, if $\{k\}$ is a block by itself, the remaining $k-1$ elements of $[k]$ must make up the other $n-1$ blocks. There are $S(k-1,n-1)$ ways in which they can do this, so there are how many partitions of $[k]$ in which $\{k\}$ is a block by itself?