Abstract Cosmological constraints of Modified Gravity (MG) models are seldom carried out rigorously. First, even though general MG models evolve differently (i.e., background and perturbations) to the standard cosmological model, it is usual to assume a ΛCDM background. This treatment is not correct and in the era of precision cosmology could induce undesired biases in cosmological parameters. Second, neutrino mass is usually held fixed in the analyses which could obscure its relation to MG parameters. In a couple of previous papers we showed that by using the Effective Fluid Approach we can accurately compute observables in fairly general MG models. An appealing advantage of our approach is that it allows a pretty easy implementation of this kinds of models in Boltzmann solvers (i.e., less error-prone) while having a useful analytical description of the effective fluid to understand the underlying physics. This paper illustrates how an effective fluid approach can be used to carry out proper analyses of cosmological constraints in MG models. We investigated three MG models including the sum of neutrino masses as a varying parameter in our Markov Chain Monte Carlo analyses. Two models (i.e., Designer f ( R ) [DES-fR] and Designer Horndeski [HDES]) have a background matching ΛCDM, while in a third model (i.e., Hu & Sawicki f ( R ) model [HS]) the background differs from the standard model. In this way we estimate how relevant the background is when constraining MG parameters along with neutrinos' masses. We implement the models in the popular Boltzmann solver CLASS and use recent, available data (i.e., Planck 2018, CMB lensing, BAO, SNIa Pantheon compilation, H 0 from SHOES, and RSD Gold-18 compilation) to compute tight cosmological constraints in the MG parameters that account for deviation from the ΛCDM model. For both the DES-fR and the HS model we obtain log 10 b < -8 at 68% confidence when all data are included. In the case of the HDES model we find a somewhat weaker value of log 10 J c > -5 at 68% confidence. We also find that constraints on MG parameters are a bit weakened when compared to the case where neutrinos' masses are held fixed in the analysis.