You can derive the probability in a manner similar to that for the usual derivation of the derangement probabilities (the probability that none of the $N$ men get their own hat back).
There are a total of $N!^n$ ways that all of the items can be distributed among the $N$ men so that each has exactly one of each type of item. Let $A_i$ denote the event that the $i$th man obtains the $n$ items that belong to him. The number of ways this can happen is $|A_i| = (N-1)!^n$, as this involves distributing all items but those belonging to the $i$th man among the other $N-1$ men. Similarly, $|A_i \cap A_j| = (N-2)!^n$, and, in general, $|A_{i_1} \cap A_{i_2} \cap \cdots \cap A_{i_j}| = (N-j)!^n$. There are also $\binom{N}{j}$ ways to choose which $j$ men will receive their own $n$ items.
Let $D(N,n,k)$ denote the number of ways that exactly $k$ of the $N$ men receive all $n$ of their items back. By the principle of inclusion/exclusion, $$D(N,n,0) = \sum_{j=0}^N (-1)^j \binom{N}{j} (N-j)!^n = N! \sum_{j=0}^N (-1)^j \frac{(N-j)!^{n-1}}{j!}.$$
Now, $D(N,n,k) = \binom{N}{k} D(N-k,n,0)$, as this is the number of ways of choosing $k$ of the $N$ men to receive all of their items back times the number of ways that none of the remaining $N-k$ men receive all of their items back.
Thus $$D(N,n,k) = \binom{N}{k} (N-k)! \sum_{j=0}^{N-k} (-1)^j \frac{(N-k-j)!^{n-1}}{j!} = \frac{N!}{k!} \sum_{j=0}^{N-k} (-1)^j \frac{(N-k-j)!^{n-1}}{j!}.$$
Dividing by $N!^n$, we have that the probability that exactly $k$ of the $N$ men receive all $n$ of their items back is
$$\frac{1}{k! N!^{n-1}} \sum_{j=0}^{N-k} (-1)^j \frac{(N-k-j)!^{n-1}}{j!}.$$
Note that this formula agrees with the values obtained by Henry for the case $N = 4$, $n=2$.
Added: In fact, the Poisson approximation suggested by Henry appears to match up well with the exact values provided by the formula given here for small values of $k$. The accuracy of the Poisson approximation appears to deteriorate, relatively speaking, as $k$ increases. However, the Poisson approach still gives a good absolute approximation when $k$ is large because the probabilities are extremely small.
Here are some observations. Start with the species $\mathcal{Q}$ of
permutations with fixed points marked, so that we have
$$\mathcal{Q} = \mathfrak{P}
\left(\mathcal{U}\mathcal{Z}
+ \mathfrak{C}_{=2}(\mathcal{Z})
+ \mathfrak{C}_{=3}(\mathcal{Z})
+ \mathfrak{C}_{=4}(\mathcal{Z})
+ \cdots\right).$$
Translating to generating functions we obtain
$$Q(z, u) = \exp\left(
uz + \frac{z^2}{2} + \frac{z^3}{3} + \frac{z^4}{4} + \cdots \right)
\\= \exp\left(uz -z + \log\frac{1}{1-z}\right)
= \frac{e^{-z}}{1-z} e^{uz}.$$
The probability that no one gets their hat is computed from the count
of permutations with no fixed points,
$$D_n = n! [z^n] [u^0] Q(z, u)
= n! [z^n] \frac{e^{-z}}{1-z}
= n! \sum_{k=0}^n \frac{(-1)^k}{k!} \approx \frac{n!}{e}.$$
It follows that the probability is $1/e.$
For the expected value of people getting their hat we obtain
$$[z^n] \left. \frac{d}{du} Q(z, u)\right|_{u=1}
= [z^n] \left. \frac{z e^{-z}}{1-z} e^{uz}\right|_{u=1}
= [z^n] \frac{z}{1-z} = 1$$
and thus there is an average of one fixed point per permutation.
As for the location of the smallest fixed point call it $q$ we can
reason as follows. Take any permutation $\sigma$ of the $n-q$ elements
that are larger than $q$ and factor it into cycles. Partition the
$q-1$ elements smaller than $q$ into a derangement containing $p$
elements and distribute the remaining $q-1-p$ elements onto the cycles
of the permutation $\sigma$ of $n-q$, which will assure that none of
them becomes a fixed point and $q$ is the smallest fixed point because
the elements less than $q$ are either in the derangement of $p$ elements
or on the cycles of $\sigma,$ where none of them is a singleton.
This gives for $p$ fixed the value
$$(n-q)! \times {q-1\choose p} \times D_p \times
(n-q)(n-q+1)(n-q+2)\\ \cdots (n-q+(q-1-p-1))
\\= (n-q)! \times {q-1\choose p} \times D_p \times
(n-q)(n-q+1)(n-q+2)\ldots (n-p-2)
\\= (n-q)! \times {q-1\choose p} \times D_p \times
{n-p-2\choose q-1-p} (q-1-p)!$$
When $q<n$ this can also be written as
$$(n-q)! \times {q-1\choose p} \times D_p \times
\frac{(n-p-2)!}{(n-q-1)!}$$
or
$$(n-q) \times {q-1\choose p} \times D_p \times (n-p-2)!$$
The fourth term in the product above on the first line i.e. the
product at the end represents the fact that we have $n-q$ choices when
placing the first element from the $q-1-p$ elements on the cycles of
$\sigma,$ and $n-q+1$ choices for the next one and so on.
The first binomial coeffcient represents the choice of $p$ elements for
the derangement.
The total contribution which is the desired sum of the smallest fixed
point over all permutations is thus given by
$$\sum_{q=1}^n q\times (n-q)! \times
\sum_{p=0}^{q-1} {q-1\choose p} \times D_p \times
{n-p-2\choose q-1-p} (q-1-p)!$$
or alternatively
$$n\times D_{n-1} + \sum_{q=1}^{n-1} q\times (n-q)\times
\sum_{p=0}^{q-1} {q-1\choose p} \times D_p \times (n-p-2)!$$
where we have used the fact that if $n$ is the first fixed point it must have been preceded by a derangement of the first $n-1$ elements.
This yields the sequence
$$1, 1, 7, 31, 191, 1331, 10655, 95887, 958879, 10547659, \ldots$$
which points us to OEIS A155521 where
additional material awaits and which would appear to confirm the above
calculation.
The following admittedly simple Maple program was used to confirm the
above formulae for small $n:$
P := proc(n)
local gf, p, pos;
option remember;
gf := 0;
for p in combinat:-permute(n) do
for pos to n do if p[pos] = pos then break end if end do;
if pos <= n then gf := gf + z^pos end if
end do;
gf
end proc
Addendum.
(Inspired by the work of @YuvalFilmus.)
We learn from the OEIS entry that by a combinatorial bijection the
sequence above is the number of permutations on $n+1$ elements having
at least two fixed points. This has EGF
$$Q_{\ge 2}(z) =
\left. (Q(z, u) - [u^0] Q(z, u) - [u^1] Q(z, u))\right|_{u=1}
= \frac{1}{1-z}
- \frac{e^{-z}}{1-z}
- \frac{ze^{-z}}{1-z}.$$
Extracting coefficients we get
$$n! [z^n] Q_{\ge 2}(z) =
n! \left(1 - \sum_{k=0}^n \frac{(-1)^k}{k!}
- \sum_{k=0}^{n-1} \frac{(-1)^k}{k!}\right)
\approx n! \left(1-\frac{2}{e}\right)$$
Dividing by the count of non-derangements we obtain
$$\frac{(n+1)!(1-2/e)}{n!(1-1/e)} =
(n+1) \times \frac{e-2}{e-1}.$$
Best Answer
You're right, the solution focuses only on matches for a particular subset of the hats. This is a general feature of the inclusion-exclusion method: The probability that all conditions in some set $S$ hold is calculated from the probabilities for each subset $S'$ of all conditions in $S'$ holding, and in calculating that probability for $S'$ the conditions not contained in $S'$ are irrelevant.