The second to eighth virial coefficients of methane were determined for temperatures up to 1200 K using an existing ab initio-based and empirically fine-tuned two-body potential combined with a new empirical nonadditive three-body potential. Nuclear quantum effects were accounted for by the semiclassical Feynman-Hibbs approach. The numerical evaluation of the high-dimensional integrals through which the virial coefficients are expressed was performed employing the Mayer-sampling Monte Carlo technique. By fitting suitable mathematical functions to the calculated virial coefficients, an analytical eighth-order virial equation of state (VEOS8) was obtained. Pressures p computed as a function of temperature T and density ρ using VEOS8 agree highly satisfactorily with p(ρ, T) values obtained with the experimentally based reference equation of state for methane of Setzmann and Wagner (SWEOS) at state points at which VEOS8 is sufficiently converged. It is shown that it is essential to account for nonadditive three-body interactions in the calculations in order to achieve good agreement with the SWEOS.