Modification of General Relativity can be invoked to explain the accelerated expansion of the Universe. Ongoing and future surveys (e.g. Euclid and LSST) will probe any departure from the predictions of GR in exquisite detail. In this talk I plan to introduce two novel techniques, based on large scale structure, to constrain modified gravity (MG) theories. In particular, I will show how the future 3D weak lensing surveys as well as the all-sky spectroscopic galaxy redshift surveys can be used very effectively to test the MG theories. As a concrete example, I will take the F(R) theories and the generalised Dilaton scalar-tensor theories. Going beyond the usual power spectrum analysis, I will show how the recently introduced integrated bispectrum (IB)  statistics can be used to tighten the constraints.  The link between the IB statistics and its generalisations to the recently discovered consistency relations will be discussed. To probe higher order statistics, I will outline how the effective field theory(EFT) based approaches can be coupled to the existing generating function formalism to compute the entire PDF, one point cumulants and their two-point generalisations.  I will also discuss briefly the relevance of these statistics for future 21cm Intensity Mapping surveys.