General method of deriving the mean field theory of a microscopic theory

In second quantization, the mean field approximation consists in approximating some combination of operators $A$ by a $c$-number $\langle A\rangle$. For example, for Bose-Einstein condenstate $A=b_0$, for Cooper pairing $A= a_{\mathbf{p}\uparrow}a_{-\mathbf{p}\downarrow}$, in the Hartree-Fock approximation $A=a_{\mathbf{p}}^+a_{\mathbf{p}}$, in the charge density wave state $A=a_{\mathbf{p}+\mathbf{q}}^+a_{\mathbf{p}}$. Here $a_i$ and $b_i$ are fermionic and bosonic annihilation operators.

The same in the Hubbard-Stratonovich transformation: we can couple any desired combination $A$ of field operators to auxiliary field.

In all mentioned cases, $A$ is indeed bosonic. Since bosonic operators can have large occupation numbers $\langle A\rangle\gg1$, deviations of $A$ from $\langle A\rangle$ can be relatively small, and also noncommutativity of $A$ and $A^+$ can be neglected (if we have the right choice of $A$): $$ \frac{\sqrt{\langle (A-\langle A\rangle)^2\rangle}}{\langle A\rangle}\rightarrow0,\qquad\frac{\langle[A,A^+]\rangle}{\langle A\rangle}\rightarrow0. $$ This makes the approximation accurate.

For fermionic quantities $A$, the mean field approximation does not have much sense because the occupation number is limited, $\langle A\rangle\sim 1$. Therefore, $A$ will be strongly fluctuating with respect to $\langle A\rangle$, and the approximation will be inaccurate.

Update number 2 following the discussion in comments

Average value of the fermionic opearator can be nonzero, $\langle A\rangle\neq0$ only at coherent mixtires of different numbers of fermions in the system, e.g. $|\psi\rangle=\alpha|N\rangle+\beta|N+1\rangle$. Generally, such superpositions of integer- and half-integer angular momentum states are forbidden by superselection rules.

In any case, the average $\langle A\rangle$ of a fermionic quantity $A$ is either zero or remains small in the macroscopic limit, so it is a bad candidate for a mean field order parameter.

Actually Wikipedia has an answer for your question,

https://en.wikipedia.org/wiki/Mean_field_theory

which will tell you how to bulid a mean field approximation self-consistently based on the Bogoliubov inequality.

If you want to know more details about the fundamental inequality,you can go through the book,Statistical Mechanics: A Set Of Lectures,written by Feynman.

Hope it helps.