On a $d$-dimensional Riemannian, spin manifold $(M,g)$ we consider non-linear, stochastic partial differential equations for spinor fields, driven by a Dirac operator and coupled to an additive Gaussian, vector-valued white noise. We extend to the case in hand a procedure, introduced in [DDRZ20] for the scalar counterpart, which allows to compute at a perturbative level the expectation value of the solutions as well as the associated correlation functions accounting intrinsically for the underlying renormalization freedoms. This framework relies strongly on tools proper of microlocal analysis and it is inspired by the algebraic approach to quantum field theory. As a concrete example we apply it to a stochastic version of the Thirring model proving in particular that it lies in the subcritical regime if $d\leq 2$.