  Despite its elegance, the theory of General Relativity is subject to
experimental, observational, and theoretical scrutiny to arrive at tighter
constraints or an alternative, more preferred theory. In alternative gravity
theories, the macroscopic properties of neutron stars, such as mass, radius,
tidal deformability, etc. are modified. This creates a degeneracy between the
uncertainties in the equation of state (EoS) and gravity since assuming a
different EoS can be mimicked by changing to a different theory of gravity. We
formulate a hierarchical Bayesian framework to simultaneously infer the EoS and
gravity parameters by combining multiple astrophysical observations. We test
this framework for a particular 4D Horndeski scalar-tensor theory originating
from higher-dimensional Einstein-Gauss-Bonnet gravity and a set of 20 realistic
EoS and place improved constraints on the coupling constant of the theory with
current observations. Assuming a large number of observations with upgraded or
third-generation detectors, we find that the $A+$ upgrade could place
interesting bounds on the coupling constant of the theory, whereas with the
LIGO Voyager upgrade or the third-generation detectors (Einstein Telescope and
Cosmic Explorer), the degeneracy between EoS and gravity could be resolved with
high confidence, even for small deviations from GR.
