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
$\begin{Bmatrix} n\\ k \end{Bmatrix}$ is the number of ways to distribute $n$ distinct objects into $k$ identical boxes with no empty boxes allowed/each box containing at least $1$ object. We will first consider distributions into $k$ distinct boxes.
The total number of ways to distribute $n$ distinct objects into $k$ distinct boxes with empty boxes allowed is: $$k^n$$
Let $A_i$ be the number of ways to perform the distribution with the condition that box $i$ must be an empty box. There are no restrictions on the other boxes. (This is equivalent to distributing the $n$ objects into the remaining $k-1$ boxes, empty boxes allowed.)
For example $A_1$ is the number of ways to distribute the objects into the $k$ boxes with box 1 being empty and is equivalent to distributing the objects into the remaining $k-1$ boxes, empty boxes allowed.
$$\implies |A_1| = (k-1)^n $$
Now we notice that $|A_1|=|A_2|=|A_3|=\dots$ since it doesn't matter which box is empty. (We are always distributing into the remaining $k-1$ boxes).
We know there are $\binom{k}{1}$ possible $A_i $
$$\sum_{i=1}^k|A_i| = \binom{k}{1}(k-1)^n$$
Using a similar argument, $|A_i\cap A_j|=(k-2)^n$ because it involves distribution into $k-2$ remaining boxes. Also, using the fact that $|A_1\cap A_2|=|A_1\cap A_3|=|A_1\cap A_4|...$ and that there are $\binom{k}{2}$ such terms,
$$\sum_{1\leq i,j\leq k}|A_i\cap A_j| = \binom{k}{2}(k-2)^n$$
On a similar note, we try extending this to when $l$ boxes must be empty:
$$\sum_{1\leq i,j\dots l\leq k}|A_i\cap A_j\cap \dots\cap A_l| = \binom{k}{l}(k-l)^n$$
$|A_1\cup A_2\cup \dots \cup A_k|$ is the number of ways to distribute $n$ distinct objects into $k$ distinct boxes with at least one empty box.
We are trying to find the distribution with no empty boxes: $$|\overline{A_1\cup A_2\cup \dots \cup A_k}|=k^n -|A_1\cup A_2\cup \dots \cup A_k|$$
By the principle of inclusion and exclusion,
$$=k^n-(|A_1|+|A_2|+|A_3|+...)+(|A_1\cap A_2|+|A_1\cap A_3|+|A_1\cap A_4|+...)-...$$
$$=k^n-\binom{k}{1}(k-1)^n+\binom{k}{2}(k-2)^n-\binom{k}{3}(k-3)^n+\dots +(-1)^l\binom{k}{j}(k-j)^n$$
$$=\sum_{l=0}^{k}(-1)^l\binom{k}{l}(k-l)^n$$
Let $j=k-l$, when $l=0$, $j=k$. When $l=k$, $j=0$.
$$=\sum_{j=k}^{0}(-1)^{k-j}\binom{k}{k-(k-j)}j^n (\text{this is technically wrong notation})$$
This is just the sum backwards. So we have:
$$=\sum_{j=0}^{k}(-1)^{k-j}\binom{k}{j}j^n$$
Now we divide by $k!$ since we are dealing with identical boxes:
$$\begin{Bmatrix} n\\ k \end{Bmatrix} =\frac{1}{k!}\sum_{j=0}^{k}(-1)^{k-j}\binom{k}{j}j^n$$