Is there some well established way on how to quantify rotations in $\mathbb{R}^n$? To say which rotation is greater and which is smaller?
In $\mathbb{R}^2$ the rotation is characterized by a single angle. This angle can be interpreted as the magnitude of the rotation.
In $\mathbb{R}^3$ however we have three possible axes around which to rotate. Furthermore we can get a zero rotation by rotating around each axis by $\pi$ radians. Therefore devices such as the sum of the three angles do not make sense.
In $\mathbb{R}^n$ the situation is even more complicated.
Let's say that the rotation is characterized by matrix $R$. Is the distance of $R$ from the identity matrix $I$ of the appropriate dimension a reasonable choice? For example $$ m_{R} = || R - I|| $$ where $||.||$ is some matrix norm. If so, which norm would be appropriate?
Every $n\times n$ rotation matrix is diagonalizable over $\mathbb C$, yielding $n$ complex eigenvalues.
In $\mathbb R^2$ the two eigenvalues are $e^{i\theta}$ and $e^{-i\theta}$. Take logs and you obtain $\pm i\theta$, which corresponds to the magnitude of the rotation being $\theta$.
In $\mathbb R^3$ the three eigenvalues are $e^{i\theta}$, $e^{-i\theta}$, and $1$. Take logs and you obtain $\pm i\theta$ and $0$, so again it makes sense to define the magnitude of the rotation as $\theta$.
In $\mathbb R^n$ there are $\lfloor n/2\rfloor$ pairs of conjugate eigenvalues, and possibly a $1$ if $n$ is odd. So you have $k=\lfloor n/2\rfloor$ different rotation angles, and it is up to you to combine them. You could take $\sqrt{\theta_1^2+\cdots+\theta_k^2}$ or $|\theta_1|+\cdots+|\theta_k|$ or $\max(\theta_1,\ldots,\theta_n)$ or something similar. Any of these would be consistent with the well-known low-dimensional cases, as long as it is symmetric in the $\theta_i$ and reduces to $\theta_1$ when $\theta_2=\cdots=\theta_k=0$.
That said, the $\|R-I\|$ idea you had can also be a reasonable choice. In the low-dimensional cases, it has a simple relationship with rotation angle.
In $\mathbb R^2$, if $R$ is a rotation by $\theta$, then $\|R-I\|$ in the Frobenius norm is $2\sqrt2|\sin(\theta/2)|$.
In $\mathbb R^3$, we can pick a basis so that the $z$-axis lies along the axis of rotation, so that $R$ becomes a block diagonal matrix consisting of a $2\times2$ rotation by $\theta$ and a $1\times1$ identity matrix. Since the Frobenius norm is invariant to orthogonal changes of basis, again $\|R-I\|_F = 2\sqrt2|\sin(\theta/2)|$.
In $\mathbb R^n$, if $R$ rotates a single $2$-dimensional subspace by $\theta$ and leaves the orthogonal subspace unchanged, we still have $\|R-I\|_F = 2\sqrt2|\sin(\theta/2)|$. In all the cases above, we can recover the angle of rotation as $\theta = 2\arcsin\big(\|R-I\|_F/2\sqrt2\big)$. Unfortunately, for a general rotation matrix, $\|R-I\|_F$ can be greater than $2\sqrt2$, so we can't recover a generalized rotation magnitude consistent with the low-dimensional cases by the above formula. Nevertheless, we can just work with $\|R-I\|_F$.
P.S. I think in general we have $\|R-I\|_F=2\sqrt2\sqrt{\sin^2(\theta_1/2) + \cdots + \sin^2(\theta_k/2)}$, but I'm not positive.
P.P.S. The overall idea is analogous to the fact that the geodesic distance between two points $x$ and $y$ on a unit sphere in $\mathbb R^n$ can be easily computed as $2\arcsin(\|x-y\|/2)$.