This really depends where you want to start from. For instance, you can construct the Choi state of $\mathcal E$, i.e.,
$$
\sigma = (\mathcal E \otimes \mathbb I)(|\Omega\rangle\langle\Omega|)\ ,
$$
with $\Omega = \tfrac{1}{\sqrt{D}}\sum_{i=1}^D |i,i\rangle$, and then extract the Kraus operators of $\mathcal E(\rho)=\sum M_i\rho M_i^\dagger$ by taking any decomposition
$$
\sigma = \sum |\psi_i\rangle\langle\psi_i|\ ,\tag{*}
$$
and writing $|\psi_i\rangle = (M_i\otimes\mathbb I)|\Omega\rangle$ (which is always possible).
Note that the decomposition $(*)$ is highly non-unique (any $|\phi_j\rangle = \sum V_{ij} |\psi_i\rangle$, with $V$ an isometry, is also a valid decomposition), which relates to the fact that the Kraus decomposition is equally non-unique. Obviously, the eigenvalue decomposition is a simple choice (which, moreover, minimizes the number of Kraus operators).
Let's look at your example in a bit more detail. Here, $D=2$. You have that
$$
\mathcal E(X)=p\mathrm{tr}(X)\,\frac{\mathbb I}{2}+(1-p)X
$$
for any $X$ (due to linearity) -- the $\mathrm{tr}(X)$ is required to make this trace-preserving for general $X$.
We now have that
\begin{align}
\sigma &= (\mathcal E \otimes \mathbb I)(|\Omega\rangle\langle \Omega|)
\\
& = \tfrac1D \sum_{ij} \mathcal E(|i\rangle\langle j|)\otimes |i\rangle\langle j|\
\end{align}
inserting the definition of $|\Omega\rangle$ and using linearity.
This yields
$$
\sigma = \frac{p}{2D}\mathbb I\otimes \sum_{i}|i\rangle\langle i| +
(1-p)\frac1D \sum_{ij}|i\rangle\langle j|\otimes |i\rangle\langle j|\ .
$$
The second term is just $(1-p)|\Omega\rangle\langle\Omega|$, and the first term is
$\frac{p}{2D}\mathbb I\otimes\mathbb I$.
You can now see that one possible eigenvalue decomposition of $\sigma$ is given by the four Bell states (I leave it to you to work out the weights), and it is well known and easy to check that that the four Bell states can be written as
$$
(\sigma_k\otimes \mathbb I)|\Omega\rangle\ ,
$$
where $\sigma_k$ are the three Pauli matrices or the identity.
Thus, you get that the $M_i$ in the Kraus representation are the Paulis and the identity, with the weight given by the eigenvalue decomposition of $\sigma$.