Any finite subgroup of $GL_n(\Bbb Z)$ is isomorphic to a subgroup of $GL_n(\Bbb F_p)$ for $p \ge 3$.
I have two questions related to this problem:
The statement is much more precise: if $G$ is a finite subgroup of $GL_n(\mathbb{Z})$ then the quotient map to $GL_n(\mathbb{F}_p)$ for an odd prime $p$ is an embedding. This follows from the claim that the kernel of this quotient map, namely the congruence subgroup consisting of invertible matrices congruent to $I \bmod p$, is torsion-free (and so does not contain any elements of finite order, and so intersects trivially with $G$), which is essentially what the post you linked proves.
I think the proof can be cleaned up by working $p$-adically: more specifically, passing to $GL_n(\mathbb{Z}_p)$ I think we can show that for $p \ge 3$, $1$ is the unique root of unity (in the ring of algebraic integers over $\mathbb{Z}_p$) congruent to $1 \bmod p$, and so if an element of the congruence subgroup has finite order then its eigenvalues are roots of unity which must then all be equal to $1$.
As for a reference, you can consult Serre's Bounds for the order of the finite subgroups of $G(k)$ which contains information on this result (although he calls the torsion-freeness above an "easy exercise") and many generalizations of it, with references. It's not written for students but you may find it at least a bit helpful.
Edit: Here are some details on a somewhat different argument from the blog post, inspired by $p$-adic ideas but not explicitly working $p$-adically. Let $\Gamma(d) = \text{ker}(GL_n(\mathbb{Z}) \to GL_n(\mathbb{Z}/d))$ denote the congruence subgroup
$$\{ M \in GL_n(\mathbb{Z}) : M \equiv I \bmod d \}.$$
We want to show that for $p$ an odd prime, $\Gamma(p)$ is torsion-free. It suffices to show that there is no $q$-torsion where $q$ is prime. So suppose $M \in \Gamma(p)$ satisfies $M^q = I$. Then the eigenvalues of $M$ must be $q^{th}$ roots of unity, and in fact the characteristic polynomial of $M$ must be a product of copies of the cyclotomic polynomials $\Phi_1(x) = x - 1$ and
$$\Phi_q(x) = \frac{x^q - 1}{x - 1} = x^{q-1} + x^{q-2} + \dots + 1.$$
On the other hand, since $M \equiv I \bmod p$, the characteristic polynomial $\varphi(x)$ of $M$ must be congruent to $(x - 1)^n$ mod $p$. So a necessary condition is that $\Phi_q(x) \equiv (x - 1)^{q-1} \bmod p$, or
$$x^q - 1 \equiv (x - 1)^q \bmod p.$$
Taking the coefficient of $x^{q-1}$ on both sides gives $q \equiv 0 \bmod p$ so $q = p$; in other words, the only possible torsion is $p$-torsion so now it suffices to rule out this possibility. (Note that so far we have not used the hypothesis that $p$ is odd, and that $\Gamma(2)$ does have $2$-torsion, since it contains $-I$.)
Write $M = I + p^k N$ where $k \ge 1$ and $N \not\equiv 0 \bmod p$. Then $M^p - I$ is equal to
$$\sum_{i=1}^p {p \choose i} p^{ki} N^i.$$
Taking both sides $\bmod p^{k+2}$ gives $p^{k+1} N + p^{kp} N^p \bmod p^{k+2}$ (all the other terms are divisible at least by $p^{2k+1}$). If either $k \ge 2$ or $p \ge 3$ then $kp \ge k+2$ so this simplifies to $p^{k+1} N \bmod p^{k+2}$ which is nonzero by hypothesis. So $M^p \neq I$ and $\Gamma(p)$ has no $p$-torsion for odd $p$ as desired.
If $k = 1$ and $p = 2$ we get $p^2 N + p^2 N^2 \bmod p^3$ which could be zero, and in fact this happens if we take $N = -I$. What we learn from this argument is that $\Gamma(2)$ can have $2$-torsion but $\Gamma(4)$ cannot (this corresponds to $k \ge 2$), and so cannot have any torsion. This kind of exceptional behavior at $p = 2$ is common $p$-adically which is a useful hint here.
