Another way is to compute manually, but first you need to extract all factors (e.g., all axes)
You should specify the nfactors
as the number of variables that you have, which in your example is 4. So it should be like like this:
pca <- principal(iris[1:4],nfactors = 4,rotate="varimax",scores=TRUE)
Then, extract the pca$loadings
after which you can compute the proportional variance by getting the sum of squared loadings per component (RC in this case) and divide them by the total sum squared loadings, which you can do by this:
colSums(pca$loadings[ , ]^2)/sum(pca$loadings[ , ]^2)
This should give you the same information as the proportion variance in the pca$loadings
, albeit for all the components (RC1 to RC4 in this case).