Another reference which takes a very explicit point of view (i.e. entirely matrix calculations and explicit circle inversions) is Beardon's *The Geometry of Discrete Groups*, which deduces the invariance in $n$ dimensions from the following useful formula: if $\sigma$ is the reflection in the Euclidean sphere of radius $r$ and centre $a$, then

$$ |\sigma(y) - \sigma(x)| = \frac{r^2 |y-x|}{|x-a||y-a|} $$

(this is formula 3.1.5 on p.26). Observe now that the extension $\tilde{\phi}$ to $\mathbb{R}^{n+1} $ of any inversion $\phi$ in $ \mathbb{R}^n$ is a reflection in some circle of radius $r$ centred on $ a $ with $a_{n+1} = 0 $, and explicitly writing down the formula for circle inversion gives that the $(n+1)$th coordinate of $\tilde{\phi}(x) $ is

$$ \frac{r^2 x_{n+1}}{|x-a|^2}. $$

But this formula is just $ x_{n+1} $ times the first displayed formula above, and this similarity suggests

\begin{align*}
  \frac{|\tilde{\phi}(y)-\tilde{\phi}(x)|^2}{\tilde{\phi}(x)_{n+1} \tilde{\phi}(y)_{n+1}} &= \frac{r^4 |y-x|^2}{|x-a|^2|y-a|^2 \frac{r^2 x_{n+1}}{|x-a|^2} \frac{r^2 y_{n+1}}{|y-a|^2}}\\
  &= \frac{|y-x|^2}{x_{n+1} y_{n+1}}
\end{align*}

i.e. the form $ |y-x|^2/x_{n+1} y_{n+1} $ is invariant under Poincare extensions of sphere inversions. The same is true for extensions of plane reflections, and since Poincare extension is a homomorphism we get that the extension of every Moebius transformation preserves that form; but the form is just the hyperbolic metric:

$$ d(x,y) = \int_\gamma \frac{|y-x|^2}{x_{n+1} y_{n+1}} ds. $$

This calculation is found on pp.34--35 of Beardon (I have added a very small amount of detail). At this point of the book he has only done the structure of Moebius maps as products of reflections, i.e. it is very hands-on and requires essentially no "heavy" machinery. The downside of this approach is that the result appears magically out of thin air as opposed to other appproaches which come from symmetric spaces and Lie groups, e.g. showing that $ \mathrm{PSL}(2,\mathbb{C}) $ is naturally isomorphic to $ O(2,1) $ as a Lie group and hence admits an action by isometries on $ \mathbb{H}^3 $ which you then show is the usual action by fractional linear transformations.