We give various integral representation formulas simultaneously for a function and its derivatives in terms of vector field gradients of the function of appropriately high order. When the function has compact support, simpler formulas can be derived. Many of the results proved here appear to be new even in the special case of classical Euclidean space. For instance, Theorem 2.2 below reduces to the following result in the usual Euclidean case: Let B be a ball in RN with radius r(B), let m be a positive integer, and let f\in Wm,1(B). Then there is a polynomial P=π_m(B,f) of degree m-1 such that for any integers i, j with 0 ≤ j < i ≤ m and a.e. x\in B, |\nablaj(f-P)(x)| ≤ C\int_{B}\frac{|\nablaif(y)|}{|x-y|N-(i-j)} dy +Cr(B)i-j-N\intB|\nablaif(y)| dy. Moreover, if 0 < i-j ≤ N, then for a.e. x \in B we have the more refined formula |\nablaj(f-P)(x)| ≤ C\intB\frac{|\nablaif(y)|}{|x-y|N-i+j} dy.