We use a set of N-body simulations employing a modified gravity (MG) model with Vainshtein screening to study matter and halo hierarchical clustering. As test-case scenarios we consider two normal branch Dvali-Gabadadze-Porrati (nDGP) gravity models with mild and strong growth rate enhancement. We study higher-order correlation functions $xi_n(R)$ up to $n=9$ and associated hierarchical amplitudes $S_n(R)equivxi_n(R)/sigma(R)^{2n-2}$. We find that the matter PDFs are strongly affected by the fifth-force on scales up to $50h^{-1}$Mpc, and the deviations from GR are maximised at $z=0$. For reduced cumulants $S_n$, we find that at small scales $Rleq10h^{-1}$Mpc the MG is characterised by lower values, with the deviation growing from $7%$ in the reduced skewness up to even $40%$ in $S_5$. To study the halo clustering we use a simple abundance matching and divide haloes into thee fixed number density samples. The halo two-point functions are weakly affected, with a relative boost of the order of a few percent appearing only at the smallest pair separations ($rleq 5h^{-1}$Mpc). In contrast, we find a strong MG signal in $S_n(R)$s, which are enhanced compared to GR. The strong model exhibits a $>3sigma$ level signal at various scales for all halo samples and in all cumulants. In this context, we find that the reduced kurtosis to be an especially promising cosmological probe of MG. Even the mild nDGP model leaves a $3sigma$ imprint at small scales $Rleq3h^{-1}$Mpc, while the stronger model deviates from a GR-signature at nearly all scales with a significance of $>5sigma$. Since the signal is persistent in all halo samples and over a range of scales, we advocate that the reduced kurtosis estimated from galaxy catalogues can potentially constitute a strong MG-model discriminatory as well as GR self-consistency test.