In 3-dimensional space, we have an explicit formula for the rotation matrix which will rotate about a vector $\vec{a} = [a_x, a_y, a_z]$. This is given by:
$$ \begin{bmatrix} \cos\theta+a_x^2(1-\cos\theta) & a_xa_y-a_z\sin\theta & a_xa_z(1-\cos\theta)+a_y \sin\theta \\ a_xa_y(1-\cos\theta)+a_z\sin\theta & \cos\theta+a_y^2(1-\cos\theta) & a_ya_z(1-\cos\theta)-a_x\sin\theta\\ a_xa_z(1-\cos\theta)-a_y\sin\theta & a_y a_z(1-\cos\theta)+a_x\sin\theta& \cos\theta+a_z^2(1-\cos\theta) \end{bmatrix} $$
In 4-d space, rotations are defined in a plane (in 3-d space as well, the rotation is happening in the plane orthogonal to the vector, $\vec{a}$).
Is there a corresponding formula that takes two vectors spanning the space that rotation is going to happen in (or the spanning the space which is going to be unaffected by the rotation) and the angle, $\theta$ we are going to rotate by as input and returns a $4 \times 4$ matrix?
Very good question! You may want to check out Givens rotation matrix $G(i,j,\theta)$, which induces rotation of $\theta$ only in the plane defined by $i$ and $j$ axis of the space, with the orthongonal subspace unchanged.
If you do a basis change $P$ such that the first axes $1,2$ correspond to the vectors spanning the rotation plane you defined. Then this shall be the required rotation matrix $$ PG(1,2,\theta)P^{-1} $$