  Within the framework of $f(R)=R+\alpha R^2$ gravity, we study realistic
models of neutron stars, using equations of state compatible with the LIGO
constraints. i.e. APR4, MPA1, SLy, and WW1. By numerically solving modified
Tolman-Oppenheimer-Volkoff equations, we investigate the Mass--Radius relation
in both metric and torsional $f(R)=R+\alpha R^2$ gravity models. In particular,
we observe that torsion effects decrease the compactness and total mass of
neutron star with respect to the General Relativity predictions, therefore
mimicking the effects of a repulsive massive field. The opposite occurs in the
metric theory, where mass and compactness increase with $\alpha$, thus inducing
an excess of mass that overtakes the standard General Relativity limit. We also
find that the sign of $\alpha$ must be reversed whether one considers the
metric theory (positive) or torsion (negative) to avoid blowing up solutions.
This could draw an easy test to either confirm or discard one or the other
theory by determining the sign of parameter $\alpha$.
