In QFT, why do fermions have to anticommute in order to insure causality? I have seen this question and I believe I understand the answer to it. However, AFAIK, only for bosons the causality condition is a vanishing commutator. For fermions we expect the anticommutator $[\phi,\phi^\dagger]_+$ to turn zero. The answer given to the question above does not seem to address this.
 A: Here is the formal answer on your question based on particular result of Pauli theorem. Calculations are rather cumbersome, but they are general.
Arbitrary fermionic field (with invariance under discrete transformations of the Lorentz group)
It can be shown that each Poincaré-covariant fermion field with half-integer spin $s = n + \frac{1}{2}$ and mass $m$ can be represented as
$$
\tag 1 \hat{\varphi}_{a} = \sum_{\sigma = -s}^{s}\int \frac{d^{3}\mathbf p}{\sqrt{(2 \pi )^{3}2E_{\mathbf p}}}\left( u_{a}^{\sigma }(\mathbf p)\hat{a}_{\sigma}(\mathbf p)e^{-ipx} + v_{a}^{\sigma } (\mathbf p)\hat{b}^{\dagger}_{\sigma}(\mathbf p)e^{ipx}\right),
$$
where
$$
\tag 2 \hat{\varphi}_{a} = \hat{\varphi}_{\mu_{0}...\mu_{n}} = \begin{pmatrix}\hat{\psi}_{\mu_{0}...\mu_{n}b} \\ \hat{\kappa}_{\ \mu_{0}...\mu_{n}}^{\dot {b}} \end{pmatrix} \in \left( \frac{n + 1}{2}, \frac{n}{2}\right) \oplus \left( \frac{n}{2}, \frac{n + 1}{2}\right),
$$
$b, \dot{b}$ are spinor indices, and
$$
\tag 3 \gamma^{\mu_{j}}\hat{\varphi}_{\mu_{0}...\mu_{j}...\mu_{n}} = 0, \quad (i\gamma^{\mu}\partial_{\mu} - m)\hat{\varphi}_{\mu_{0}...\mu_{n}} = 0, 
$$
$$
\tag 4 \quad \partial^{\mu_{j}} \hat{\varphi}_{\ \ \mu_{0}...\mu_{j}...\mu_{n}} = 0, \quad u^{\sigma}_{a}(\mathbf p ) = (-1)^{s + \sigma}\gamma_{5}v^{-\sigma }_{a}(\mathbf p).
$$
As you can see, for $n = 0$ eqs. $(1)-(4)$ give the Dirac field. 
The massless case is the same, but $\sigma $ may take only following sets of values: $\{-s\}$, $\{s\}$, or $\{-s,s\}$.
You can read about general aspects of irreducible Poincaré representations which is realised as a sum of fields of creators and annihilators in Weinberg QFT Vol. 1 (chapter about general causal fields). Equations $(3)-(4)$ may be given as the requirement that field $(1)$ transforms under an irreducible representations of the Poincaré group with spin $s$ and mass $m$.
Causality for fermionic theories and anticommutator
From the causality principle we must have 
$$
\tag 5 [\hat{\varphi}_{a}(x), \hat{\bar{\varphi }}_{b}(y)]_{\pm} = 0, \quad (x - y)^{2} < 0 , \quad g_{00} = 1.
$$
In $3 + 1$-dimensional spacetime and for indistinguishable particles, the first homotopy group of the configuration space with 3 spatial dimensions is the permutation group $S_N$. This means that the clockwise and anticlockwise exchange of two particles is equal and hence there are only two possible statistics, Bose-Einstein or Fermi-Dirac,
$$
\tag 6 [\hat{a}_{\sigma}(\mathbf p ), \hat{a}^{\dagger}_{\sigma {'}}(\mathbf k) ]_{\pm} = \delta (\mathbf p - \mathbf k ) \delta_{\sigma \sigma {'}}.
$$ 
Let's use $(1)-(4)$ and $(6)$ for clarifying the statistic which is obeyed by field $(1)$.
We have by using these equations and the identity $[\gamma_{\mu} , \gamma_{5}]_{+} = 0$ that
$$
\tag 7 [\hat{\varphi}_{a}(x), \hat{\bar{\varphi }}_{b}(y)]_{\pm} = \sum_{\sigma}\int \frac{d^{3}\mathbf p}{(2 \pi )^{3}2E_{\mathbf p}}\left( u^{\sigma}_{a}(\mathbf p)\bar{u}^{\sigma}_{b}(\mathbf p) e^{-ipX} \mp \gamma_{5}u^{\sigma}_{a}(\mathbf p)\bar{u}_{b}^{\sigma}(\mathbf p)e^{ipX}\right), 
$$
where $X = x - y$ and $\hat{\bar {\varphi}} = \hat{\varphi}^{\dagger}\gamma_{0}$.
By using the second identity of $(3)$ and the requirement of Lorentz covariance, it can be shown that $\sum_{\sigma}u^{\sigma}_{a}(\mathbf p)\bar{u}_{b}^{\sigma}(\mathbf p) = R_{ab}(p)$ can be represented in a form
$$
R_{ab}(p) = (\gamma^{\mu}p_{\mu} + m)P_{ab}(p), 
$$
where $P_{ab}(p)$ is constructed from sum of summands with even number of momentums and even number of gamma-matrices and from summands with odd numbers of momentums and odd number of gamma-matrices. So by using $[\gamma_{\mu} , \gamma_{5}]_{+} = 0$ again, we may get
$$
[\hat{\varphi}_{a}(x), \hat{\bar{\varphi }}_{b}(y)]_{\pm} = \int \frac{d^{3}\mathbf p}{(2 \pi )^{3}2E_{\mathbf p}}\left( R_{ab}(p) e^{-ipX} \mp R_{ab}(-p)e^{ipX}\right) = 
$$
$$
\tag 8  = R_{ab}\left( i\partial_{x}\right)(D(X) \mp D(-X)), \quad D(X) = \int \frac{d^{3}\mathbf p }{(2 \pi )^{3}2E_{\mathbf p}}e^{-ipX}. 
$$
It can be shown that for $X^{2} < 0$, the function $D(X)$ satisfies the relation $D(X) = D(-X)$, so $(8)$ vanishes if and only if fermionic fields have Fermi-Dirac statistics.
