Abstract: The equation of state (EoS) of neutron star matter plays a key role in both the structure and evolution of a neutron star. However, as lattice QCD faces significant difficulties in simulating dense matter, only effective models can be relied on to unveil their properties, such as the microscopic many body theories and density functional theories. Due to the uncertainties reside in model assumptions and parameters, the matter state and composition of dense matter are still unclear and exhibit large ambiguities, leading to large uncertainties in the corresponding EoSs. The measurements of neutron star properties in multi-messenger era, on the other hand, have placed a stringent constraint on its EoS. For example, the measurements of2M_⊙ pulsars have ruled out various soft EOSs for dense stellar matter. Based on the multi-messenger observations of the binary neutron star merger event GW170817, the tidal deformability of 1.4M_⊙ neutron star is believed to lie within 70≤Lambda1.4≤580 with the corresponding radii R=（11.9±1.4）km, excluding the EoSs which are too stiff at small densities. The masses and radii of the two millisecond pulsars PSR J0030+0451 and PSR J0740+6620 were measured accurately via pulse-profile modeling with NICER and XMM-Newton data, where similar radii (about 12.4 km) were obtained despite their large differences in masses. This reduces the likelihood of a strong first-order phase transition inside neutron stars with masses M≤2M_⊙. Adopting various EoS parameterizations which are either phenomenological or nuclear physics-motivated, extensive investigations were carried out by applying those observational data to constraining the EoSs and neutron star properties with Bayesian parameter estimation. On the other hand, as proposed by Biswas, one can take a complementary approach to estimate statistically the most preferred EoS according to the current observations. In light of the recently obtained 10 unified neutron star EoSs using relativistic mean field models with nonlinear self couplings (NL3, PK1, TM1, GM1, and MTVTC) and density-dependent couplings (DD-LZ1, DDME-X, PKDD, DD-ME2, and TW99), combined with the other 6 unified EoSs (FSU2, DD2, NL3wr-L55, FSU2R, FSU2H, and TM1e) which were obtained in a similar manner, those EoS models then can be studied for Bayesian model selection, which is not included in Biswas’s original study. In this work, a Bayesian model selection on the 16 unified neutron star EoSs predicted by relativistic mean field models is thus performed, which are available on CompOSE. In particular, the tidal deformability measurements were used from the binary neutron star merger event GW170817 and the simultaneous mass-radius measurements of PSR J0030+0451 and PSR J0740+6620 by the NICER collaboration. The most preferred EoS model is DD2, which predicts the radius and tidal deformability of a 1.4M_⊙ neutron star to be 13.19 km and 687, respectively. The selected EoS models that are not ruled out by observation can then be ordered with respect to the relative odds ratio: DD2, TW99, DD-LZ1, DD-ME2, TM1e, FSU2H, DDME-X, PKDD, FSU2R, and MTVTC. These EoSs predict that the radii and tidal deform abilities of 1.4M_⊙ neutron stars should lie within 12.3-13.6 km and 400-784, while the maximum masses range from 2.02M_⊙ to 2.56M_⊙.