Magnetohydrodynamical simulations show that the magnetic field can drive molecular outflows during the formation of massive protostars. The best probe to observationally measure both the morphology and the strength of this magnetic field at scales of 10-100au is maser polarization. Measuring the direction of magnetic fields at milliarcsecond resolution around a sample of massive star forming regions to determine whether there exists a relation between the orientation of the magnetic field and of the outflows. In addition by estimating the magnetic field strength via the Zeeman splitting measurements, the role of magnetic field in the dynamics of the massive star-forming region is investigated. We selected a flux-limited sample of 31 massive star-forming regions to perform a statistical analysis of the magnetic field properties with respect to the molecular outflows characteristics. We report the linearly and circularly polarized emission of 6.7GHz CH_3_OH masers towards seven massive star-forming regions of the total sample with the European VLBI Network. The sources are: G23.44-0.18, G25.83-0.18, G25.71-0.04, G28.31-0.39, G28.83-0.25, G29.96-0.02, and G43.80-0.13. We identified a total of 219 CH_3_OH maser features, 47 and 2 of which showed linearly and circularly polarized emission, respectively. We measured well ordered linear polarization vectors around all the massive young stellar objects and Zeeman splitting towards G25.71-0.04 and G28.83-0.25. Thanks to recent theoretical results, we were able to provide lower limits to the magnetic field strength from our Zeeman splitting measurements. We further confirm (based on ~80% of the total flux-limited sample) that the magnetic field on scales of 10-100 au is preferentially oriented along the outflow axes. The estimated magnetic field strength of |B_||_|>61mG and >21mG towards G25.71-0.04 and G28.83-0.2, respectively, indicates that it dominates the dynamics of the gas in both regions.