Harnack estimates for degenerate parabolic equations modeled on the subelliptic p-Laplacian