Antisymmetric wave functions instantly imply the Pauli exclusion principle – essentially because $\psi(x,x)=-\psi(x,x)=0$, to write the concept schematically – which implies that the occupation numbers are $N=0,1$ and statistical physics is therefore inevitably governed by the Fermi-Dirac statistics which may be derived from Boltzmann/statistical physics for these occupation numbers.
Similarly, symmetric wave functions imply that particles are indistinguishable but the occupation numbers may be $N=0,1,2,3,\dots$. That implies the Bose-Einstein distribution by applying the Boltzmann steps to the multiparticle states with these occupation numbers. For the proofs of the first two paragraphs of my answer, see
How to derive Fermi-Dirac and Bose-Einstein distribution using canonical ensemble?
The bulk of the spin-statistics theorem is to link the antisymmetric functions with the half-integer spin and symmetric functions with integer spin. It was proved by Pauli and all the evidence available to me suggests that you haven't seen a clear proof because you haven't tried. I won't reproduce the full proof here because I don't believe it would be a good investment of time but I will give a sketch. The Lagrangian for a spin-0 real field $\phi$ has to contain the kinetic term
$$\frac{1}{2} \partial_\mu \phi \partial^\mu \phi $$
which is dictated by the Lorentz symmetry etc. If $\phi$ were anticommuting, the object above would identically vanish and there would be no dynamics. For spin-1/2 fields, it's the other way around (the Majorana kinetic term would vanish if the fields were not anticommuting), and so on. For the wrong combinations of spin and statistics, Pauli actually showed one can't have positive norms and/or Hamiltonians bounded from below.
