In view of the expectation that the existence of complex poles is a signal of confinement, we investigate the analytic structure of gluon, quark, and ghost propagators in the Landau gauge QCD and QCD-like theories by employing an effective model of Yang-Mills theory with a gluon mass term, which we call the massive Yang-Mills model. In this model, we particularly investigate the number of complex poles in the parameter space of the model consisting of gauge coupling constant, gluon mass, and quark mass for the gauge group SU (3) and various numbers of quark flavors NF within the asymptotic free region. This investigation extends the previous result obtained for the pure Yang-Mills theory with no flavor of quarks NF = 0 that the gluon propagator has a pair of complex conjugate poles and the negative spectral function while the ghost propagator has no complex pole. The gluon and quark propagators at the best-fit parameters for NF = 2 QCD have one pair of complex conjugate poles as in the zero flavor case. By increasing quark flavors, we find a new region in which the gluon propagator has two pairs of complex conjugate poles for light quarks with the intermediate number of flavors 4 N f < 10. However, the gluon propagator has no complex poles if very light quarks have many flavors N f ≥ 10 or both of the gauge coupling and quark mass are small. In the other regions, the gluon propagator has one pair of complex conjugate poles. Moreover, as a general feature, we argue that the gluon spectral function of this model with nonzero quark mass is negative in the infrared limit. In sharp contrast to gluons, the quark and ghost propagators are insensitive to the number of quark flavors within the current approximations adopted in this paper. These results suggest that details of the confinement mechanism may depend on the number of quark flavors and quark mass.