A die is rolled 10 times, what is the probability that 2 numbers will occur exactly 4 times? So I know that the probability of rolling any one number exactly 4 times would be $$\frac{10!}{4!6!}\cdot \left(\frac{1}{6}\right)^4 \cdot\left(\frac{5}{6}\right)^6,$$
so would the probability that two numbers occur exactly four times be?
$$\frac{10!}{8!2!}\cdot\left(\frac{1}{6}\right)^4\cdot\left(\frac{1}{6}\right)^4\cdot\left(\frac{4}{6}\right)^2$$
No, you have ignored two effects. The first is that you seem to be trying to calculate the probability that two specific numbers occur four times each, and the problem asks for any two numbers. The second is that once one number is given to have occured exactly four times, the odds of the other number also doing so are different that the original odds.
The easy way to do this is to ask how many ways to arrange to get two four-timers are there, and divide by $6^{10}$. So:
There are $\binom{6}{2}$ ways to choose the two numbers.
There are $\binom{10}{4}$ ways to choose which of the dice landed on the specific first number.
There are then $\binom{6}{4}$ ways to choose which of the dice landed on the specific second number.
There are then $4^2$ ways to choose the rolls of the remaining two dice.
The net probability is $$ \frac{1}{6^{10}} \cdot \frac{6!}{2!4!}\cdot\frac{10!}{6!4!}\cdot \frac{6!}{4!2!}\cdot 4^2 $$
After your clarifications, here's an approach to compute favorable ways:
[ Choose faces to come up ] $\times$ [ Permute ]
There are two possibilities:
$4-4-2\;of\;a\;kind: [\binom62\binom41] \times [\frac{10!}{4!4!2!}]$
$4-4-1-1\;of\;a\;kind: [ \binom62\binom42] \times [\frac{10!}{4!4!1!1!}]$
Add up to get favorable ways, and divide by $6^{10}$ to get the final result
We can ask about the more general problem of rolling an $n$-sided die $m$ times and observing exactly $p$ values that occur $q$ times. This corresponds to the marked combinatorial species
$$\mathfrak{S}_{=n}(\mathfrak{P}_{=0}(\mathcal{Z}) + \mathfrak{P}_{=1}(\mathcal{Z}) + \cdots + \mathcal{U} \mathfrak{P}_{=q}(\mathcal{Z}) + \mathfrak{P}_{=q+1}(\mathcal{Z}) + \cdots).$$
What we have here is a sequence of sets containing the positions from among the $m$ rolls where a value from the $n$ possible values occured, where instances of sets of $q$ positions are marked. We then obtain the mixed generating function
$$G(z, u) = \left(\exp(z) - \frac{z^q}{q!} + u\frac{z^q}{q!}\right)^n.$$
On extracting the coefficient on $[u^p]$ we get
$$m! [z^m] [u^p] G(z, u) = m! [z^m] {n\choose p} \frac{z^{pq}}{(q!)^p} \left(\exp(z) - \frac{z^q}{q!}\right)^{n-p} \\ = m! [z^m] {n\choose p} \frac{z^{pq}}{(q!)^p} \frac{z^{(n-p)q}}{(q!)^{n-p}} \left(\frac{q!\exp(z)}{z^q} -1\right)^{n-p} \\ = m! [z^m] {n\choose p} \frac{z^{nq}}{(q!)^n} \sum_{k=0}^{n-p} {n-p\choose k} (-1)^{n-p-k} \frac{(q!)^k \exp(kz)}{z^{kq}} \\ = m! [z^m] {n\choose p} \frac{1}{(q!)^n} \sum_{k=0}^{n-p} {n-p\choose k} (-1)^{n-p-k} (q!)^k \exp(kz) z^{nq-kq} \\ = m! {n\choose p} \frac{1}{(q!)^n} \sum_{k=0}^{n-p} {n-p\choose k} (-1)^{n-p-k} (q!)^k [z^{m+kq-nq}] \exp(kz) \\ = m! {n\choose p} \frac{1}{(q!)^n} \sum_{k=0}^{n-p} {n-p\choose k} (-1)^{n-p-k} (q!)^k \frac{k^{m+kq-nq}}{(m+kq-nq)!}.$$
Taking into account that we must have $m+kq-nq\ge 0$ we obtain that $k\ge n - m/q.$ The problem definition tells us that $m\ge pq$ (zero possibilities otherwise) and hence $m/q \ge p.$ Thus we finally obtain the probability
$$\bbox[5px,border:2px solid #00A000]{ \frac{m!}{n^m} {n\choose p} \sum_{k=\lceil n - m/q\rceil}^{n-p} {n-p\choose k} \frac{(-1)^{n-p-k}}{(q!)^{n-k}} \frac{k^{m+kq-nq}}{(m+kq-nq)!}.}$$
Observe that this will yield $${\frac {875}{69984}}$$ for the question posed by the OP ($n=6, m=10, p=2, q=4$) which incidentally matches the answer by @MarkFischler. Here we have used that ${n-p\choose k}$ will be zero should $\lceil n - m/q\rceil \lt 0$ and $k$ is negative at the lower end of the range.
There is some very basic Maple code to verify the correctness of the formula.
PROB :=
proc(n, m, p, q)
option remember;
local ind, digits, rolls, res, valset, match;
if m < p*q then return 0 fi;
res := 0;
for ind from n^m to 2*n^m-1 do
digits := convert(ind, base, n);
rolls := digits[1..m];
valset := convert(rolls, `multiset`);
match := select(ent->ent[2] = q, valset);
if nops(match) = p then
res := res + 1;
fi;
od;
res/n^m;
end;
PROB2 :=
proc(n, m, p, q)
option remember;
local ind, digits, rolls, res, valset, match;
if m < p*q then return 0 fi;
res := 0;
for ind from n^m to 2*n^m-1 do
digits := convert(ind, base, n);
rolls := digits[1..m];
valset := convert(rolls, `multiset`);
match :=
add(`if`(ent[2] = q, 1, 0), ent in valset);
if match = p then
res := res + 1;
fi;
od;
res/n^m;
end;
X := (n, m, p, q)->
m!/n^m*binomial(n,p)
*add(binomial(n-p,k)*(-1)^(n-p-k)/(q!)^(n-k)
*k^(m+k*q-n*q)/(m+k*q-n*q)!, k=ceil(n-m/q)..n-p);
