您实际上在问两个不同的问题:
- 如何在 PCA 投影后对数据进行聚类。
- 如何获得上述图。
然而,在讨论这些之前,我想补充一点,如果您的样本位于列中,那么您没有正确执行 PCA。您应该在转置数据集上执行此操作,如下所示:
model <- prcomp(t(d), scale=TRUE)
但要使其发挥作用,您必须删除数据中的所有常量行。
现在我假设您已按照您想要的方式完成了 PCA 步骤。
当您指定 retX=TRUE 时,prcomp 返回旋转矩阵(默认情况下为 true)。所以你会想使用model$x
.
下一步是根据主成分对数据进行聚类。这可以通过多种方式来完成。一是层次聚类。如果你最终想要 5 组,这是一种方法:
fit <- hclust(dist(model$x[,1:3]), method="complete") # 1:3 -> based on 3 components
groups <- cutree(fit, k=5) # k=5 -> 5 groups
此步骤将为您提供稍后用于着色的组。
最后一步是绘图。在这里,我编写了一个简单的函数来一次完成所有操作:
library(rgl)
plotPCA <- function(x, nGroup) {
n <- ncol(x)
if(!(n %in% c(2,3))) { # check if 2d or 3d
stop("x must have either 2 or 3 columns")
}
fit <- hclust(dist(x), method="complete") # cluster
groups <- cutree(fit, k=nGroup)
if(n == 3) { # 3d plot
plot3d(x, col=groups, type="s", size=1, axes=F)
axes3d(edges=c("x--", "y--", "z"), lwd=3, axes.len=2, labels=FALSE)
grid3d("x")
grid3d("y")
grid3d("z")
} else { # 2d plot
maxes <- apply(abs(x), 2, max)
rangeX <- c(-maxes[1], maxes[1])
rangeY <- c(-maxes[2], maxes[2])
plot(x, col=groups, pch=19, xlab=colnames(x)[1], ylab=colnames(x)[2], xlim=rangeX, ylim=rangeY)
lines(c(0,0), rangeX*2)
lines(rangeY*2, c(0,0))
}
}
这个函数很简单:它需要两个参数:1)分数矩阵,其中主成分在列中,样本在行中。如果你想要(例如)第一个、第二个和第四个组件,你基本上可以使用 model$x[c(1,2,4)] 。 2)聚类的组数。
然后,它根据传递的主成分和图对数据进行聚类(2D 或 3D,具体取决于传递的列数)
以下是一些例子:
plotPCA(model$x[,1:2], 5)
3D 示例(基于 3 个第一主成分):
plotPCA(model$x[,1:3], 5)
最后一个图将是交互式的,因此您可以将其旋转或放大/缩小。
希望这可以帮助。