It isn't so much that we want to compare the $p$-refinement and $h$-refinement errors directly, instead we want to compare the convergence properties (e.g. speed) of each refinement strategy. This requires more knowledge of the constant in the apriori error estimate. We'll illustrate by looking at the apriori error estimate of a discontinuous Galerkin spectral method (which can be interpreted as a high-order FEM).
It is true that the constant $C$ in the inequality $$ e \leq Ch^{p+1}, $$ depends on the approximation order $p$. Thus, it is not immediately clear what to expect when we increase the order $p$. However, the error inequality above hides a lot of information about the convergence behaviour in that constant. In $p$-refinement it becomes important how the $C$ depends on $p$.
In a 1D spectral element method the error estimate on a given element $k$ has a form (derived in a similar fashion as in Canuto et. al.)$$ e \leq \mathcal{C}h^{\min(m,p+1)}p^{\frac{3}{2}-m}\|q^k\|_{H^m(\Omega)}, $$
where $\mathcal{C}$ is independent of $h$ and $p$, $q^k$ is the true solution on element $k$, and $H^m$ is the $m$ order Sobolev space. In we apply the numerical method to a problem with a "smooth" solution it will generally be true that $p+1\leq m$. Now we can more readily see the difference in $p$ versus $h$ refinement.
If we fix $p$ then the terms $$ p^{\frac{3}{2}-m}\|q^k\|_{H^m(\Omega)} $$ are a constant that we absorb into $\mathcal{C}$ to obtain the $C$ in the first error estimate. Now we can perform a standard $h$-refinement argument and see that the convergence will be similar to that of a finite difference scheme, i.e., some polynomial order on $h$.
However, if we fix the grid $h$ then the term $h^{\min(m,p+1)}$ is constant and the error estimate looks like $$ e \leq \hat{\mathcal{C}}p^{\frac{3}{2}-m}\|q^k\|_{H^m(\Omega)},$$
where $\hat{\mathcal{C}} = \mathcal{C}h^{\min(m,p+1)}$. Now we immediately see that the rate of convergence changes in the $p$-refinement study. Provided $m>\frac{3}{2}$ the error decays as a power of $p$. For "very smooth" solutions, e.g. infinitely differentiable, the decay in the error is faster than any polynomial order. Sometimes this is referred to as exponential convergence because if you look at a semilog plot of the error against $p$ you see a straight line.
As a quick aside we note that in spectral methods the convergence rate depends on the smoothness of the solution. Heuristically, we can interpret this as "the larger $m$ is in the error estimate, the smoother the solution," because the true solution is in a higher order Sobolev space, i.e., the true solution is square integrable up to the $m^{\textrm{th}}$ order derivative. The convergence of a Sobolev norm is often viewed in terms of the decay of Fourier coefficients of the solution $q^k$ (particularly for fractional derivatives). This is one reason the family of methods are referred to as spectral methods.