We consider the p and hp versions of the virtual element method ([1,2]), for which the available theoretical error estimates, which are significantly sub-optimal with respect to the order of the method, are known to be quite pessimistic, as shown by the numerical evidence. In this work we look at the approximation properties of the polynomial part of the discrete solution. By carefully taking into account the role played by the stabilization term, we obtain sharper error bounds, which come much closer to the numerical results presented in [1]. We also propose a new, easily computable, design for the stabilization term, also leading to an improved theoretical error bound. For the different stabilization terms, we then compare, both theoretically and numerically, the condition number of the stiffness matrix, which is also a criterion to be considered in the choice of the stabilization strategy. This work is realized with the support of the Italian Ministry of Research, under the complementary action NRRP “D34Health - Digital Driven Diagnostics, prognostics and therapeutics for sustainable Health care” (Grant #PNC0000001)