The stellar initial mass function (IMF) is a fundamental property in the measurement of stellar masses and galaxy star formation histories. In this work we focus on the most massive galaxies in the nearby universe log(M /M ) > 11.2. We obtain high quality Magellan/LDSS-3 long slit spectroscopy with a wide wavelength coverage of 0.4µm − 1.01µm for 41 early-type galaxies (ETGs) in the MASSIVE survey, and derive high S/N spectra within an aperture of R e /8. Using detailed stellar synthesis models, we constrain the elemental abundances and stellar IMF of each galaxy through full spectral modeling. All the ETGs in our sample have an IMF that is steeper than a Milky Way (Kroupa) IMF. The best-fit IMF mismatch parameter, α IMF = (M/L)/(M/L) MW , ranges from 1.12 to 3.05, with an average of α IMF = 1.84, suggesting that on average, the IMF is more bottomheavy than Salpeter. Comparing the estimated stellar mass with the dynamical mass, we find that most galaxies have stellar masses smaller than their dynamical masses within the 1σ uncertainty. We complement our sample with lower-mass galaxies from the literature, and confirm that log(α IMF ) is positively correlated with log(σ), log(M ), and log(M dyn ). The IMF in the centers of more massive ETGs is more bottom-heavy. In addition, we find that log(α IMF ) is positively correlated with both [Mg/Fe] and the estimated total metallicity [Z/H]. We find suggestive evidence that the effective stellar surface density Σ Kroupa might be responsible for the variation of α IMF . We conclude that σ, [Mg/Fe] and [Z/H] are the primary drivers of the global stellar IMF variation.