Let me continue my comment with an illustration for the case where you're interested in existing R packages. There are several package in the Multivariate Task View that will provide enhanced method for PCA-related methods (as compared to R base
princomp), e.g. ade4 or FactorMineR. I personally like FactoMineR because of its simple syntax, and you check the associated website for more information on the available methods.
One can use supplementary categorical and/or numerical variables when applying a PCA. Those variables are not used to construct factor axes, but can be showed afterwards on the correlation circle (for numerical variables) or the individual map (for categorical variables). Here is a toy example of use (from the on-line help):
res.pca <- PCA(decathlon, quanti.sup = 11:12, quali.sup=13)
If you have multiple passive variables, you can select the one to display (with or without confidence ellipses) using the
keepvar= argument. Here is another picture with two illustrative variables.
Be careful with arguments that are a little bit non-standard if you are used to default plotting functions in R. The
plotellipses() function makes use of the helper function
ellipse::ellipse that you can use (or not) in any plot (look for
monpanel.ellipse subfunction in
plotellipses() to see how confidence lines are computed). That's what I did to build specific individual map (B&W, different plotting symbol, etc.). For example, the following snippet just plot all individuals with two different symbols depending on the type of sporting event (2004 Olympic Game or 2004 Decastar):
labs <- paste(round(res.pca$eig[1:2, 2], 2), "%", sep="")
xlab=paste("Dim. 1 (", labs, ")", sep=""),
ylab=paste("Dim. 2 (", labs, ")", sep=""))
abline(v=0, h=0, lty=2)
Besides, I would like to point you to @vqv's excellent ggbiplot package, available on GitHub, which follows from one of his answer. (It uses R base functions and ggplot2.)