1

我在 stats.stackexchange 中尝试了这个问题，有人建议我在这里尝试一下，所以这里是：

我已经用 VEGAN 包在 R 中完成了一些关于树木健康的生态数据的 PCA 分析。总共有 80 棵树（因此，80 个“地点”）分为四个处理类别。我已经用颜色编码点绘制了数据 - 根据治疗组的颜色。我不想在 PCA 双图上绘制单个站点/树，而是制作类似于盒须图的东西，它有四个“十字”，显示每个组的质心和两个 PCA 维度的 SE。我在论文中看到过这样的数字，但我似乎找不到这样的 R 脚本。有什么建议么？（我想在此处发布我正在寻找的示例图片，但我能找到的都是收费的，抱歉）。

我想另一种选择是只获取站点分数并手动查找方法和 SE 并创建我自己的情节，但如果可能的话，我宁愿为它找到一个脚本。

我一直在运行的代码非常简单：

p1<-princomp(scale(health, scale=T))
summary(p1)
scores(p1)
plot(p1)
loadings(p1)
biplot(p1, xlab = "PC 1 (38%)", ylab = "PC 2 (22%)",cex=0.6)
plot(p1$scores[,1],p1$scores[,2])
names(p1)

plot(p1$scores[,1],p1$scores[,2], type='n', xlab="PC I", ylab="PC II")
text(p1$scores[,1],p1$scores[,2] labels=Can$tree)
4

1 回答 1

0

我可能会先ordiellipse看看这是否适合您的需求。

### Reproducible example
require("vegan")
data(varespec)
data(varechem)

pca <- rda(varespec, scale = TRUE)
grp <- with(varechem, cut(Baresoil, 4, labels = 1:4))
cols <- c("red","orange","blue","forestgreen")
scl <- 1 ## scaling

plot(pca, display = "sites", scaling = scl, type = "n")
points(pca, display = "sites", scaling = scl, col = cols[grp], pch = 16)
lev <- levels(grp)
for (i in seq_along(lev)) { ## draw ellipse per group
  ordiellipse(pca, display = "sites", kind = "se", scaling = scl,
              groups = grp, col = cols[i], show.groups = lev[i])
}
## centroids
scrs <- as.data.frame(scores(pca, display = "sites", scaling = scl, 
                             choices = 1:2))
cent <- do.call(rbind, lapply(split(scrs, grp), colMeans))
points(cent, col = cols, pch = 3, cex = 1.1)

这产生

在此处输入图像描述

您可以points()从上面的代码中删除该行以阻止它绘制实际样本，但我认为它在理解ordiellipse正在做什么方面具有指导意义。

在图中，质心通过分组标记为每个轴上的站点得分的平均值grp。椭圆是一个连续区域（给定我在ordiellipse()调用中选择的设置）关于该质心的 1 个标准误差。您对每个方向上的误差线的建议是椭圆的特定情况，ordiellipse如果您要计算质心的标准误差，它们应该在水平和垂直方向上延伸到椭圆的极值点.

但是，这将无法考虑两个轴上分数的协方差。请注意，在下面的示例中，在那些与轴成角度的椭圆中，标准误差条不会在其极值点与椭圆相交。如果您要绘制一个包含由误差线定义的区域的框，它将包含椭圆，但它对质心的不确定性给出了非常不同的印象。

serrFun <- function(df) {
  apply(df, 2, function(x) sd(x) / sqrt(length(x)))
}
serr <- do.call(rbind, lapply(split(scrs, grp), serrFun))
for (i in seq_along(lev)) {
    arrows(cent[i, 1] - serr[i, 1], cent[i, 2],
           cent[i, 1] + serr[i, 1], cent[i, 2],
           col = cols[i], code = 3, angle = 90, length = 0.05)
    arrows(cent[i, 1], cent[i, 2] - serr[i, 2],
           cent[i, 1], cent[i, 2] + serr[i, 2], 
           col = cols[i], code = 3, angle = 90, length = 0.05)
}

在此处输入图像描述

于 2014-03-21T17:29:02.550 回答