我看着統計::: biplot.default和統計::: biplot.prcomp,我接近你想要什麼。您可以修改此代碼以滿足您的需求。 (我使用了虹膜數據集)
require(plyr)
data(iris)
pcr <- prcomp(~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, data = iris)
indiv <- data.frame(pcr$x[,1:2])
indiv$species <- iris$Species
column <- data.frame(pcr$rotation[ ,1:2])
n <- nrow(indiv)
eigenval <- pcr$sdev[1:2]
eigenval <- eigenval * sqrt(n)
indiv <- transform(indiv, pc1 = PC1/eigenval[1], pc2 = PC2/eigenval[2])
column <- transform(column, pc1 = PC1 * eigenval[1], pc2 = PC2 * eigenval[2])
### based on stats:::biplot.default
unsigned.range <- function(x) c(-abs(min(x, na.rm = TRUE)), abs(max(x, na.rm = TRUE)))
rangx1 <- unsigned.range(indiv[, 4])
rangx2 <- unsigned.range(indiv[, 5])
rangy1 <- unsigned.range(column[, 3])
rangy2 <- unsigned.range(column[, 4])
mylim <- range(rangx1, rangx2)
ratio <- max(rangy1/rangx1, rangy2/rangx2)
nspecies <- table(iris$Species)
# compute the convex hull for each species
hull <- dlply(indiv[,1:3], .(species), chull)
# get points connected
hull <- llply(hull, function(x) c(x, x[1]))
plot(pc2 ~ pc1, data = indiv, cex = 0.5, col = c("blue", "yellow", "green")[iris$Species], xlim = mylim, ylim = mylim)
lines(indiv$pc1[hull$setosa], indiv$pc2[hull$setosa] , col = "blue")
lines(indiv$pc1[cumsum(nspecies)[1] + hull$versicolor], indiv$pc2[cumsum(nspecies)[1] + hull$versicolor], col = "yellow")
lines(indiv$pc1[cumsum(nspecies)[2] + hull$virginica], indiv$pc2[cumsum(nspecies)[2] + hull$virginica], col = "green")
par(new = TRUE)
plot(pc1 ~ pc2, data = column, axes = FALSE, type = "n", xlim = mylim * ratio, ylim = mylim * ratio, xlab = "", ylab = "")
text(column$pc1, column$pc2, labels = rownames(column), cex = 0.5, col = "red")
arrows(0, 0, column$pc1 * 0.8, column$pc2 * 0.8, col = "red", length = 0.2)
axis(3, col = "red")
axis(4, col = "red")
如果您還沒有找到它,請查看'chull'(簡稱'凸包')。它位於基本的R軟件包'grDevices'中,所以只需輸入'?chull',以獲得如何從一組x-y座標計算和繪製最小凸多邊形的示例。 – 2012-01-04 16:58:50
另一個指針:要獲得PC空間中點的x-y座標,輸入'pcr $ x'。 – 2012-01-04 17:07:07
也可以看看來自素食包 – EDi 2012-01-04 17:37:35