How does one go about constructing a mixed entangled state?

An example of a state which yeilds an entangled state when you trace one qubit out, is the 3-qubit "W" state: $$ lvert W_3 rangle = tfrac{1}{sqrt3} Bigl( lvert 100 rangle + lvert 010 rangle + lvert 001 rangle Bigr) $$ Taking the outer product with itself, we obtain $$ lvert W_3 rangle!langle W_3 rvert = tfrac{1}{3} Bigl( !begin{aligned}[t]& lvert 100 rangle!langle 100 rvert + lvert 100 rangle!langle 010 rvert + lvert 100 rangle!langle 001 rvert phantom{Big)} \&+ lvert 010 rangle!langle 100 rvert + lvert 010 rangle!langle 010 rvert + lvert 010 rangle!langle 001 rvert \&+ lvert 001 rangle!langle 100 rvert + lvert 001 rangle!langle 010 rvert + lvert 001 rangle!langle 001 rvert Bigr) end{aligned} $$ If we trace out the third qubit, we then obtain the state $$ begin{align} rho = mathrm{tr}_3Bigl(lvert W_3 rangle!langle W_3 rvertBigr) &= tfrac{1}{3} Bigl( lvert 10 rangle!langle 10 rvert + lvert 10 rangle!langle 01 rvert + lvert 01 rangle!langle 10 rvert + lvert 01 rangle!langle 01 rvert + lvert 00 rangle!langle 00 rvert Bigr) \&= tfrac{1}{3} lvert 00 rangle!langle 00 rvert + tfrac{2}{3} lvert Psi^+ rangle!langle Psi^+ rvert, end{align} $$ where in particular $lvert Psi^+ rangle = tfrac{1}{sqrt 2}bigl( lvert 01 rangle + lvert 10 rangle bigr)$ is a Bell state.

Note that $rho$ is a rank-2 operator with different eigenvalues: any decomposition of $rho$ as a convex combination of other density operators, can only involve terms whose eigenvectors are supported on $mathrm{span},{ lvert 00 rangle, lvert Psi^+ rangle }$. Any mixed tensor product density operator $rho_A otimes rho_B$ has an eigenbasis consisting of two or more product states; but the only product state contained in $mathrm{span},{ lvert 00 rangle, lvert Psi^+ rangle }$ is the state $lvert 00 rangle$ itself. It follows that $rho$ cannot be decomposed as a convex combination of products of possibly-mixed states, and is entangled.

More generally, if for $n > 1$ we define $lvert W_n rangle$ as the analogous $n$-term uniform superposition of standard basis states with a single 1, we may describe it as $$ lvert W_n rangle = tfrac{sqrt{n{-}1}}{sqrt n} lvert W_{n{-}1}ranglelvert0rangle + tfrac{1}{sqrt n} lvert00cdots0ranglelvert1rangle $$ so that $$ mathrm{tr}_nBigl(lvert W_n rangle!langle W_n rvertBigr) = tfrac{n{-}1}{n} lvert W_{n-1} rangle!langle W_{n-1} rvert + tfrac{1}{n} lvert 00cdots0rangle!langle 00cdots 0rvert $$ which for larger values of $n$ is closer and closer to being a pure entangled state, while still being mixed for any finite $n$.

A mixed separable state is written in the form $$ rho=sum_ip_isigma^A_iotimessigma^B_i $$ where the $sigma_i$ are valid density matrices on a single site.

The example you give, say $rho=frac12|00ranglelangle00|+frac12|11ranglelangle 11|$ is exactly of this form. Specifically, $$ p_0=p_1=frac12,qquad sigma^A_0=sigma^B_0=|0ranglelangle 0|,qquad sigma^A_1=sigma^B_1=|1ranglelangle 1|. $$

In general, it can be tricky to definitively prove that such a decomposition does not exist. However, in the case of two-qubits, there's an if and only if condition: the state is entangled if and only if its partial transpose is not non-negative (i.e. contains a negative eigenvalue).

The classic example is the Werner state, $$ rho=frac{1-p}{4}I+p|psiranglelanglepsi| $$ where $|psirangle$ is a two-qubit Bell state. This can be written out as a matrix $$ rho=left(begin{array}{cccc} frac{1-p}{4} & 0 & 0 & 0 \ 0 & frac{1+p}{4} & -frac{p}{2} & 0 \ 0 & -frac{p}{2} & frac{1+p}{4} & 0 \ 0 & 0 & 0 & frac{1-p}{4} end{array}right) $$ If we take the partial transpose of this, we get $$ rho^T=left(begin{array}{cccc} frac{1-p}{4} & 0 & 0 & -frac{p}{2} \ 0 & frac{1+p}{4} & 0 & 0 \ 0 & 0 & frac{1+p}{4} & 0 \ -frac{p}{2} & 0 & 0 & frac{1-p}{4} end{array}right) $$ This has a negative eigenvalue if $p/2>(1-p)/4$.