  We explore observational constraints on possible deviations from Newtonian
gravity by means of large-scale clustering of galaxies. We measure the power
spectrum and the bispectrum of Sloan Digital Sky Survey galaxies and compare
the result with predictions in an empirical model of modified gravity. Our
model assumes an additional Yukawa-like term with two parameters that
characterize the amplitude and the length scale of the modified gravity. The
model predictions are calculated using two methods; the second-order
perturbation theory and direct N-body simulations. These methods allow us to
study non-linear evolution of large-scale structure. Using the simulation
results, we find that perturbation theory provides reliable estimates for the
power spectrum and the bispectrum in the modified Newtonian model. We also
construct mock galaxy catalogues from the simulations, and derive constraints
on the amplitude and the length scale of deviations from Newtonian gravity. The
resulting constraints from power spectrum are consistent with those obtained in
our earlier work, indicating the validity of the previous empirical modeling of
gravitational nonlinearity in the modified Newtonian model. If linear biasing
is adopted, the bispectrum of the SDSS galaxies yields constraints very similar
to those from the power spectrum. If we allow for the nonlinear biasing
instead, we find that the ratio of the quadratic to linear biasing
coefficients, b_2/b_1, should satisfy -0.4 < b_2/b_1<0.3 in the modified
Newtonian model.
