It's worth noting that the 3-dimensional case is the most general case. If $x(s) = p + s u$ and $y(t) = q + t v$ generate the lines then the distance is the minimum norm of the residual $r(s, t) = x(s) - y(t) = w + s u - t v$ where $w = p-q$. The residual lies in the vector space spanned by $u$, $v$ and $w$. Within that space it occupies the 2-dimensional affine plane through $w$ with directions $u$ and $v$, assuming $u$ and $v$ are independent. We get the minimum norm residual by projecting parts of $w$ parallel to the $uv$ plane. If $u$ and $v$ are unit vectors, we get an orthogonal basis for the plane: $\hat{u} = u, \hat{v} = (1 - u u^T) v$. The distance is $\|(1 - \hat{u} \hat{u}^T - \hat{v} \hat{v}^T /(\hat{v}^T \hat{v})) w\|$.
Compared to cross product solution for the 3-dimensional case, you'll note that projecting onto the cross product of $u$ and $v$ is equivalent to the projection operator $1 - \hat{u} \hat{u}^T - \hat{v} \hat{v}^T /(\hat{v}^T \hat{v})$ when $u$ and $v$ are independent. More generally, if $P$ is an orthogonal projection onto a subspace then $1-P$ is the orthogonal projection onto the subspace's orthogonal complement.
What I went through is really just a special case of solving a least-squares problem with Gram-Schmidt. You can calculate the distance between arbitrary affine subspaces by a similar process. Let $A, B$ be linear maps into a vector space $V$ and let $a, b$ be points in $V$. Then $a + \text{im}(A)$ and $b + \text{im}(B)$ are affine subspaces of $V$. Their distance is the minimum norm of the residual $r(x, y) = a-b + Ax - By$. This is just a least-squares problem for a block matrix $C$ composed of $A$ and $-B$. Let $z = [x, y]$, $c = a-b$ and $C = [A, -B]^T$. Then $c + Cz = a-b + Ax - By$.