"Simple English":
I do not understand why the intrinsic angular momentum of a particle has anything to do with whether it can coexist with another identical particle
The coexistence of a particle with another (identical) particle is quantified by the value of the commutator. The commutator depends on the particular representation (how the particle transforms under Poincaré tranformations, that is Lorentz boosts, rotations and translations). The representation of a particle is uniquely labelled by its spin and mass.
So this is how the spin enters the picture.
And I really don't understand why the distinction would be
half-integer vs integer values
The integer and half-integer values comes out from (simple) angular momentum algebra, see later.
How does it enter the particle statistics business? If you do the maths (see later), you'll see that the spin enters as $2j$ and not just $j$.
$\newcommand{\mb}{\mathbf}$
$\newcommand{\im}{\mathrm{i}}$
Let's work in 3D and assume you know that there are only two particle statistics, such that:
$$ [a, a^\dagger ] = 0 \quad \text{for bosons}, \tag{1}$$
$$ \{a, a^\dagger \} = 0 \quad \text{for fermions}, \tag{2}$$
where $a^\dagger$ ($a$) is the creation (annihilation) operator for the quanta of the field.
There is a mathematical (topological) way of proving this, and an intuitive picture for visualising the topologically distinct paths in $n$ dimensions. But that's another matter.
How does spin come in?
The above distribution statistics for quantum particles stems from the requirement of indistinguishability between identical particles. In the true sense o the word, identical particles transform in the same irreducible representation of the Poincaré group, uniquely identified by their mass and spin. Hence, quantum mechanics is not enough and we need to do quantum field theory.
Quantum field theory's main addition to the pot is that it allows for the theory to be causal.
For a theory to be causal, the time ordering of physical events affecting the system's evolution cannot be reversed. This is especially problematic for space-like separations where a Lorentz boost may reverse the chronological order $t_{\mathrm{final}} - t_{\mathrm{initial}} < 0$. For causality to hold, any two space-like separated operators are required to commute:
\begin{equation}
[\mathcal{O}_1(x), \mathcal{O}_2(y)] = 0 \quad \text{if} \quad (x-y)^2 < 0, \quad g_{\mu \nu} = (+,-,-,-),
\end{equation}
to ensure their time ordering is irrelevant and not resulting in any physical consequence.
The relativistic counterparts of the commutation relations in eqs. $1$ and $2$ are:
\begin{equation}
\begin{gathered}
\text{bosons: }\qquad [a(\mb{p}, s), a^{\dagger}(\mb{p}', s')] = 2E_{\mb{p}} (2\pi)^3 \delta^{(3)}(\mb{p} - \mb{p}') \delta_{s, s'}, \\
\text{fermions: } \qquad \{a(\mb{p}, s), a^{\dagger}(\mb{p}', s') \} = 2E_{\mb{p}} (2\pi)^3 \delta^{(3)}(\mb{p} - \mb{p}') \delta_{s, s'} .
\end{gathered}
\end{equation}
A general field operator $\phi$ for non-interacting particles $A$ and $B$ is given by:
\begin{equation} \label{eq:fieldeq}
\begin{gathered}
\phi_{A, B}(x) = \int \frac{\mathrm{d}^3 \mb{p}}{(2\pi)^3} \frac{1}{2 E_{\mb{p}}} \sum_s \left [ \mathrm{e}^{-\im px}f_{A, B}(\mb{p},s)a(\mb{p},s) + \mathrm{e}^{\im px}h_{A, B}(\mb{p},s)a^{\dagger}(\mb{p},s) \right ]_{p^0 = E_{\mb{p}}}.
\end{gathered}
\end{equation}
Because operators $\mathcal{O}(x)$ are usually just a product of $\prod_i \phi_i(x)$, requiring $[\mathcal{O}_1(x), \mathcal{O}_2(y)]=0$ is the same as requiring $\left [ \phi_A(x), \phi_B(y)\right ] = 0$.
This results in:
\begin{equation} \label{eq:commutators}
\begin{gathered}
\text{for bosons: } \quad \left [ \phi_A(x), \phi_B(y)\right ] \propto 1 - (-1)^{2j} , \\
\text{for fermions: } \quad \left [ \phi_A(x), \phi_B(y)\right ] \propto 1 + (-1)^{2j},
\end{gathered}
\end{equation}
where $j$, in the absence of orbital angular momentum, is the particle's spin. See later for a bit of a proof.
The ladder operators used in angular momentum algebra result in a quantised spectrum bounded by $j$ and $-j$, so that $-j + n = j$ with $n \in \mathbb{N}$. This requires $j = n/2$, i.e. either an integer or a half integer.
Hence, in conclusion:
\begin{equation}
\begin{gathered}
\text{bosons} \Leftrightarrow \left [ a,a^\dagger \right ] = 0 \Leftrightarrow \text{integer spin}, \\
\text{fermions} \Leftrightarrow \left \{ a,a^\dagger \right \} = 0 \Leftrightarrow \text{half-integer spin}.
\end{gathered}
\end{equation}
Sketch of a proof
The foundation of the proof lies in the commutator $[\phi_A, \phi_B]$. The full, rigorous and index-full proof is given in The Quantum Field of Fields I by Weinberg, in the section about General causal fields. I am just going to sketch the proof here, as a guide for the argument presented in the reference.
The commutator $[\phi_A, \phi_B]$ depends on the commutators $[f_A, f_B]$ and $[h_A, h_B]$ of the functions that define the field operator. Hence we will turn our attention to the coefficients $f$ and $h$.
The function $f(\mb{p})$ for a generic momentum $\mb{p}$ can be generated from one at rest $f(\mb{0})$ by a Lorentz boost $\Lambda$: $f(\mb{p}) = D(\Lambda)f(\mb{0})$, with $D$ being some irreducible representation of the Lorentz group. The latter affects both space and time, and therefore obeys $SO(4)$ algebra -- technically $SO(1,3)$ -- which is isomorphic to $SO(3)\times SO(3)$. This can be easily understood by considering the same isomorphism in one fewer dimension: $SO(3)$ is isomorphic to $SO(2)\times SO(2)$, because any position on a 3D sphere can be written as a linear superposition of two orthogonal 2D circles. Hence, we can replace the four-dimensional vector $\mathcal{J}\in SO(4)$ by two three-dimensional vectors $\mb{J}, \mb{K}\in SO(3)$, with quantum numbers $j$ and $k$ obeying the usual angular momentum algebra and hence $2j+1, 2k+1 \in \mathbb{N}$.
Luckily we already know that the spherical harmonics $Y^m_{\ell}$ are the basis functions of the irreducible $2\ell+1$ space, and they transform according to $Y^m_{\ell}(-\mathbf{r}) = (-1)^{\ell} Y^m_{\ell}(\mathbf{r})$. Extending this to the generic representation $D(\Lambda)$, we expect a transformation $D(\Lambda) \rightarrow (-1)^j (-1)^k = (-1)^{j+k}$.
Identical particles transform according to the same irreducible representation of the Poincaré group, uniquely labelled by their mass $m = \sqrt{\frac{1}{c^2}\,p^\mu p_\mu}$ and spin. If $A$ and $B$ are identical, that is with equal mass and spin, then:
\begin{equation}
[f_A, f_B] \propto D(\Lambda)D(\Lambda) \propto (-1)^{j+k} (-1)^{j+k} = (-1)^{2(j+k)}.
\end{equation}
Hence, the total angular momentum $j+k$ determines the sign of the commutator and the statistics obeyed by the particles, resulting in the expressions above.
The treatment in this section was simplified by only considering real scalar fields, though it can be easily extended to complex fields representing charged particles. More importantly, the fields were assumed to be non-interacting and hence expressed in terms of single mode $\mb{p}$ creation and annihilation operators. Proofs of the spin-statistics theorem for interacting fields are much harder, and to be quite frank I do not understand them. However, it is argued in Weinberg that the cluster decomposition theorem ensures interacting fields can always be constructed from the non-interacting components considered here.
