I object to the term "Levi-Civita symbol", and will call them by a more descriptive, and less inaccurate, term "permutation symbol".
Your observation about the permutation symbol in 2-D is a good one.
I don't know what form an answer to your question would take.
(I would love to hear one!)
Of course, $-i$ is not the only antisymmetric (bi-)linear
transformation: $i$ does just as well, the difference
being described geometrically as sense of rotation. In the four-dimensional
case, there are more transformations that are like rotation.
The definition you referred to (in terms of the sign of permutations)
specifies individual elements of a matrix representation, which is good for a starting point.
The rotations and the permutation symbol are manifestations of rather
different things, geometrically: rotations are simple transformations of the
space, whereas the permutation symbol is related to measure (area, volume) of
regions in the space.
The former can always be represented by square matrices, whereas the
complexity of the latter grows much faster — the summation
convention was invented to deal with that greater complexity. Only in two
dimensions does the permutation symbol coincide with a rotation — but in all
dimensions, the two things are related.
The direct generalization from the 2D case to 4D is the
four-dimensional permutation symbol, which can be imagined as a
$4\times 4\times 4\times 4$ nested block of numbers.
It is related to the quaternion units $i$, $j$, and $k$, but the relation
takes a little work to describe.
Note that, in contrast to the case of complex numbers, there are two
linear transformations associated with a quaternion: one for multiplication on
the left by the quaternion, and one for the right, and that they are typically
different.
Thus, multiplication by a quaternion is not represented by a single
matrix, but two.
For the $4\times 4$ matrices representing
(with respect to the usual basis for quaternions $\{ 1, i, j, k \}$)
quaternion multiplication on the left and right by $i$, $j$, and $k$,
write
\begin{align*}
i_L &=
\left[\begin{matrix}
0 & -1 & & \\
1 & 0 & & \\
& & 0 & -1 \\
& & 1 & 0 \\
\end{matrix}\right] ,
&j_L &=
\left[\begin{matrix}
& & -1 & 0 \\
& & 0 & 1 \\
1 & 0 & & \\
0 & -1 & & \\
\end{matrix}\right] ,
&k_L &=
\left[\begin{matrix}
& & & -1 \\
& & -1 & \\
& 1 & & \\
1 & & & \\
\end{matrix}\right] ; \\
i_R &=
\left[\begin{matrix}
0 & -1 & & \\
1 & 0 & & \\
& & 0 & 1 \\
& & -1 & 0 \\
\end{matrix}\right] ,
&j_R &=
\left[\begin{matrix}
& & -1 & 0 \\
& & 0 & -1 \\
1 & 0 & & \\
0 & 1 & & \\
\end{matrix}\right] ,
&k_R &=
\left[\begin{matrix}
& & & -1 \\
& & 1 & \\
& -1 & & \\
1 & & & \\
\end{matrix}\right] .
\end{align*}
The symbol $\epsilon_{\mu\nu\rho\sigma}$ can be written
as a $4 \times 4$ matrix of $4 \times 4$ matrices, so:
$$
\left[\begin{matrix}
0
& \left[\begin{matrix}
0 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 1\\0 & 0 & -1 & 0
\end{matrix}\right]
& \left[\begin{matrix}
0 & 0 & 0 & 0\\0 & 0 & 0 & -1\\0 & 0 & 0 & 0\\0 & 1 & 0 & 0
\end{matrix}\right]
& \left[\begin{matrix}
0 & 0 & 0 & 0\\0 & 0 & 1 & 0\\0 & -1 & 0 & 0\\0 & 0 & 0 & 0
\end{matrix}\right]\\
\left[\begin{matrix}
0 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & -1\\0 & 0 & 1 & 0
\end{matrix}\right]
& 0
& \left[\begin{matrix}
0 & 0 & 0 & 1\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\\-1 & 0 & 0 & 0
\end{matrix}\right]
& \left[\begin{matrix}
0 & 0 & -1 & 0\\0 & 0 & 0 & 0\\1 & 0 & 0 & 0\\0 & 0 & 0 & 0
\end{matrix}\right]\\
\left[\begin{matrix}
0 & 0 & 0 & 0\\0 & 0 & 0 & 1\\0 & 0 & 0 & 0\\0 & -1 & 0 & 0
\end{matrix}\right]
& \left[\begin{matrix}
0 & 0 & 0 & -1\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\\1 & 0 & 0 & 0
\end{matrix}\right]
& 0
& \left[\begin{matrix}
0 & 1 & 0 & 0\\-1 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0
\end{matrix}\right]\\
\left[\begin{matrix}
0 & 0 & 0 & 0\\0 & 0 & -1 & 0\\0 & 1 & 0 & 0\\0 & 0 & 0 & 0
\end{matrix}\right]
& \left[\begin{matrix}
0 & 0 & 1 & 0\\0 & 0 & 0 & 0\\-1 & 0 & 0 & 0\\0 & 0 & 0 & 0
\end{matrix}\right]
& \left[\begin{matrix}
0 & -1 & 0 & 0\\1 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0
\end{matrix}\right]
& 0
\end{matrix}\right]
$$
(I generated this using sympy.)
Here, $\mu$ and $\nu$ index the containing matrix, while $\rho$ and $\sigma$
index the inner matrices, both in row-column order.
In this form, the matrices for multiplication by quaternion units
can be read off in terms of blocks of the matrix for the permutation symbol:
\begin{eqnarray*}
i_L &=\quad \epsilon_{4 3 \rho \sigma}
+ \epsilon_{3 1 \rho \sigma} \,, \qquad
i_R &=\quad \epsilon_{4 3 \rho \sigma}
+ \epsilon_{1 2 \rho \sigma} \,, \\
j_L &=\quad \epsilon_{2 4 \rho \sigma}
+ \epsilon_{2 1 \rho \sigma} \,, \qquad
j_R &=\quad \epsilon_{2 4 \rho \sigma}
+ \epsilon_{1 3 \rho \sigma} \,, \\
k_L &=\quad \epsilon_{3 2 \rho \sigma}
+ \epsilon_{4 1 \rho \sigma} \,, \qquad
k_R &=\quad \epsilon_{3 2 \rho \sigma}
+ \epsilon_{1 4 \rho \sigma} \,.
\end{eqnarray*}
To express the permutation symbol in terms of the quaternion units, arrange
the six matrices to form two antisymmetric $4 \times 4$ matrices of matrices,
so:
$$
U_L = \left[\begin{matrix}
0 & -i_L & -j_L & -k_L \\
i_L & 0 & -k_L & j_L \\
j_L & k_L & 0 & -i_L \\
k_L & -j_L & i_L & 0
\end{matrix}\right] , \quad
U_R = \left[\begin{matrix}
0 & -i_R & -j_R & -k_R \\
i_R & 0 & k_R & -j_R \\
j_R & -k_R & 0 & i_R \\
k_R & j_R & -i_R & 0
\end{matrix}\right] .
$$
From here, it is a matter of verification that
$$
\epsilon_{\mu\nu\rho\sigma} = \frac{1}{2} ( U_L - U_R ) .
$$
I am at a loss to answer "why", beyond a simple verification, but
I expect that a succinct, convincing proof, based on symmetries, is achievable.