For positive \(p\)-harmonic functions on Riemannian manifolds, we derive a gradient estimate and Harnack inequality with constants depending only on the lower bound of the Ricci curvature, the dimension \(n\), \(p\) and the radius of the ball on which the function is defined. Our approach is based on a careful application of the Moser iteration technique and is different from Cheng-Yau's method employed by Kostchwar and Ni, in which a gradient estimate for positive \(p\)-harmonic functions is derived under the assumption that the sectional curvature is bounded from below.