R语言中如何提取PCA分析95%置信椭圆外的异常值数据？

假设你已经完成PCA分析，且有包含PCA坐标（PC1、PC2）和分组信息的数据框（示例命名为pca_df），以下是具体实现步骤：

1. 计算95%置信椭圆的边界坐标 #

ggplot2的stat_ellipse()默认用MASS包的稳健协方差估计生成椭圆，我们复用该逻辑提取椭圆边界点：

# 加载依赖包
library(ggplot2)
library(MASS)
library(ellipse)

# 定义函数：按分组生成椭圆坐标
get_ellipse <- function(data, group_col, conf_level = 0.95) {
  ellipse_list <- lapply(unique(data[[group_col]]), function(g) {
    group_data <- subset(data, data[[group_col]] == g)
    # 计算稳健协方差（与stat_ellipse默认逻辑一致）
    cov_result <- cov.trob(group_data[, c("PC1", "PC2")])
    # 生成椭圆坐标
    ellipse_coords <- ellipse(cov_result$cov, centre = cov_result$center, level = conf_level)
    # 转成带分组的数据集
    data.frame(PC1 = ellipse_coords[,1], PC2 = ellipse_coords[,2], group = g)
  })
  do.call(rbind, ellipse_list)
}

# 生成椭圆边界数据集
ellipse_df <- get_ellipse(pca_df, group_col = "你的分组列名")

2. 识别并提取椭圆外的样本 #

使用sp包的point.in.polygon()判断样本点是否在对应分组的椭圆内：

library(sp)

# 给原数据集添加"是否在椭圆内"的标记
pca_df$in_ellipse <- sapply(1:nrow(pca_df), function(i) {
  current_group <- pca_df$你的分组列名[i]
  current_ellipse <- subset(ellipse_df, group == current_group)
  # 返回值：1=内部，2=边界，0=外部
  point_status <- point.in.polygon(pca_df$PC1[i], pca_df$PC2[i], 
                                   current_ellipse$PC1, current_ellipse$PC2)
  point_status != 0
})

# 提取椭圆外的异常值
outliers <- subset(pca_df, in_ellipse == FALSE)

3. 结果验证（可选） #

将异常值标在原图上确认：

ggplot(pca_df, aes(x = PC1, y = PC2, color = 你的分组列名)) +
  geom_point() +
  stat_ellipse(level = 0.95) +
  # 用红色高亮异常值
  geom_point(data = outliers, aes(x = PC1, y = PC2), color = "red", size = 3)

关键说明 #

• 替换代码中的你的分组列名为实际列名（如pop、population）
• 若PCA结果来自poppr的pca()函数，可通过pca_scores <- as.data.frame(scores(你的PCA对象))提取坐标，再与样本信息合并成pca_df
• 若需使用普通协方差（对应stat_ellipse(type="norm")），将代码中的cov.trob替换为cov.wt即可

内容的提问来源于stack exchange，提问作者jkr68