Inverse of a skew-symmetric matrix For $a, x, y, z \in \mathbb R$, let $$M= \left(
  \begin{array}{cccc}
    \cos(a)  & \sin(a) \, x & \sin(a)\,  y & \sin(a) \, z \\
    -\sin(a) \, x & \cos(a) & \sin(a) \,z & -\sin(a)\, y \\
   -\sin(a) \, y & -\sin(a) \, z & \cos(a) & \sin(a) \, x \\
  -\sin(a) \, z & \sin(a) \, y & -\sin(a)\, x & \cos(a)
  \end{array}
\right).$$ 
Without doing much calculation, why the matrix $M-I_4$ is invertible and why its inverse is given $$(M-I_4)^{-1}=\,\frac{-1}{2} I_{4} - \frac{\cot(\sqrt{x^2+y^2+z^2} )}{2\sqrt{x^2+y^2+z^2}} A,$$
where $A$ is the skew-symmetric matrix given by 
$$A=\left(
  \begin{array}{cccc}
    0 & x & y &  z \\
    -x & 0 & z & -y \\
    -y & -z & 0 & x \\
    -z & y & -x & 0
  \end{array}
\right).$$
Thank you in advance

Here is a proof of the first part (inversibility of $M-I_4$)
and a computation of the inverse (though I do not get a result in the form of yours).
One obtains an efficient simplification by using half-angle formulas (https://en.wikipedia.org/wiki/Tangent_half-angle_substitution) :
$$\cos(a)=\dfrac{1-t^2}{1+t^2} \ \  \text{and} \ \ \sin(a)=\dfrac{2t}{1+t^2}$$ 
Let $N=M-I_4$.
All entries in matrix $N$ have a common factor $\dfrac{2t}{1+t^2}$. Therefore we can set
$$N=\dfrac{2t}{1+t^2}P \ \ \text{with} \ \  
P:=\begin{bmatrix}-t&x&y&z\\-x&-t&z&-y\\-y&-z&-t&x\\-z&y&-x&-t\end{bmatrix}$$
Therefore, an equivalent issue is to show the invertibility of $P$.
Please note that $P=A-tI_4$ with matrix $A$ as you have defined it.
It turns out that the determinant of $P$ is very compact :
$$\det(P)=(t^2 + x^2 + y^2 + z^2)^2\tag{1}$$
which is non-zero, excepted... exceptional cases $x=y=z=0$ and $t=0$ (the last case corresponding to angles $a=k\pi, k \in \mathbb{Z}$). Set apart these cases, $P$ is always invertible.
Now, what is the inverse of $M-I$ ? The inverse of $P$, considered as a quaternion (see remark 1 below) is
$$P^{-1}=\dfrac{1}{n}\begin{bmatrix}-t&-x&-y&-z\\x&-t&-z&y\\y&z&-t&-x\\z&-y&x&-t\end{bmatrix}=\dfrac{1}{n}(-A-tI_4)\tag{2}$$
$$\text{with }n:=t^2+x^2+y^2+z^2\tag{3}$$ 
(one can easily check (2) by multiplying $P$ and $P^{-1}$).
Therefore, as $M-I_4=\dfrac{2t}{1+t^2}P=\sin(a)P$, using (2) and (3), we get
$$(M-I_4)^{-1}=\dfrac{1}{\sin(a)}P^{-1}=\dfrac{1}{n\sin(a)}(A-tI_4)\tag{4}$$
which is different from your formula (as it has been remarked, this formula needs to depend on $a$).
Remarks : 
1) It is not surprizing to get formula (1) : it is to be related to the $4 \times 4$ matrix representation of quaternions (see paragraph "Matrix representation" in https://en.wikipedia.org/wiki/Quaternion) and their so-called "norm". By the way, I wouldn't be that surprised that your issue comes from quaternionic representations...
2) Using (1) with $t=0$, we get 
$$\det(A)=(x^2 + y^2 + z^2)^2\tag{2}$$
3) Formula (1) has to be connected to the so-called Pfaffian. If we have had to compute the determinant of a skew-symmetric matrix of even order, it is good to know that there exists a formula as a square of a polynomial expression in its entries called its associated Pfaffian. See https://en.wikipedia.org/wiki/Pfaffian where you will find the following formula for order $4$ :
$$\det\begin{bmatrix}0&a&b&c\\-a&0&d&e\\-b&-d&0&f\\-c&-e&-f&0\end{bmatrix}=(af-be+dc)^2$$
