This is not correct, p-element software is available to compute engineering quantities for structures with nonlinear materials (StressCheck, for instance). Might not do rubbers, but that's a limitation of the software, not the p-version. More than a decade ago I used p-version to model nonlinear viscoelastic materials, which worked quite well; I don't recall any h-version FE software being able to handle those material types, though I am sure they could, if the FE companies bothered to develop them.
I can imagine the numerical difficulties intrinsic to finite differences contribute significantly to their now being out of favor. Anybody who has ever run an experiment in which an output signal had to be differentiated I am sure is aware that integrating is intrinsically more stable than differentiating.
The 'proof' why 'hp' refinements converge more quickly (that is, for fewer degrees of freedom) can be found in the following references:
[1] Babuska, I, "The p- and hp-versions of the Finite Element Method: the state of the art." Finite Elements: Theory and Applications, ed. D.L. Dwoyer, M.Y. Hussaini, and R.G. Voigt, pp. 199-239, 1988
[2] Gui, W., and Babuska, I, "The h- p- and hp-versions of the finite element method in one dimension. Part I: The error analysis of the p-version. Part II: The error analysis of the h and hp-versions. Part III: The Adaptive hp-version," Numerische Mathematik, Vol. 49, pp. 577-612, 613-657, and 659-683, 1986.
[3] Babuska, I, and Suri, M. "The p and hp-versions of the finite element method. An Overview," Comp. Meth. in Appl. Mech. and Engng, Vol 80, pp. 5-26.
Warning! These are written in the language of a branch of high level mathematics called functional analysis. There'll be many terms unfamiliar to 99% of all engineers: Sobolev spaces, Lebesque integration, etc.
Those are the mathematics based proof of superiority of hp-version over h-version.