R 语言｜主成分分析 PC1 和 PC2 得分有相关性吗？

目录

引言
1 使用真实数据
2 使用虚拟数据
  2.1 使用随机数据演示
  2.2 多次模拟
结语

引言

主成分分析（PCA）是一种统计方法，它通过正交变换将一组可能相关的变量转换为一组线性不相关的变量，这组变量称为主成分。

因此，在理想情况下，PCA 得到的各主成分轴之间应该是完全不相关的，即它们之间的相关系数应该为零。这是 PCA 的一个关键特性。

1 使用真实数据

实测的 iris 数据，

• 著名的（Fisher’s or Anderson’s）iris 数据集。给出了三种鸢尾花各 50 朵花的萼片长度和宽度，以及花瓣长度和宽度的测量值，单位为厘米。

# 使用 iris 数据演示
pca_fit <- prcomp(scale(iris[, -5]))

pca_dat <- data.frame(pca_fit$x)

# 拟合 PC1 和 PC2 关系

lm_fit <- lm(pca_dat$PC2 ~ pca_dat$PC1)

# 可视化

par(family = "sans", cex = 1.5, col.main = "green4")
plot(
  pca_dat$PC1,
  pca_dat$PC2,
  pch  = 19,
  cex  = 0.3,
  xlab = "PC1",
  ylab = "PC2",
  main = sprintf("slope = %.5f", coef(lm_fit)[2])
)
abline(lm_fit, col = "green4")

2 使用虚拟数据

2.1 使用随机数据演示

# 使用随机数据演示
set.seed(1)
dat2     <- matrix(rnorm(20 * 100), ncol = 20)
pca_fit2 <- prcomp(scale(dat2))
pca_dat2 <- data.frame(pca_fit2$x)

par(family = "sans", cex = 1.5, col.main = "green4")
plot(
  pca_dat2$PC1,
  pca_dat2$PC2,
  pch  = 19,
  cex  = 0.3,
  xlab = "PC1",
  ylab = "PC2",
  main = sprintf("slope = %.5f", coef(lm_fit)[2])
)
lm_fit <- lm(pca_dat2$PC2 ~ pca_dat2$PC1)
abline(lm_fit, col = "green4")

2.2 多次模拟

先模拟 10 次看看，

# 模拟 10 次看看相关不相关
sapply(
  1:10,
  function(i) {
    set.seed(i)
    d <- matrix(rnorm(20 * 100), ncol = 20)
    r <- data.frame(prcomp(scale(d))$x)
    f <- lm(r$PC2 ~ r$PC1)
    sprintf("slope = %.5f", coef(f)[2])
  }
)

 [1] "slope = 0.00000"  "slope = 0.00000"  "slope = 0.00000"  "slope = -0.00000"
 [5] "slope = 0.00000"  "slope = 0.00000"  "slope = 0.00000"  "slope = 0.00000" 
 [9] "slope = 0.00000"  "slope = -0.00000"

模拟 1000 次看看，

# 模拟 1000 次看看相关不相关
slope_df <- do.call(
  rbind.data.frame,
  lapply(
    1:1000,
    function(i) {
      # set.seed(i)
      r <- rnorm(20 * 100) |> 
        matrix(ncol = 20) |> 
        scale() |> 
        (function(i) prcomp(i)$x)() |> # 非匿名函数
        # (\(i) prcomp(i)$x)() |> # 匿名函数
        data.frame()
      f <- lm(r$PC2 ~ r$PC1)
      s <- data.frame(coef(f)[2])
      colnames(s) <- "slope"
      rownames(s) <- NULL
      return(s)
    }
  )
)

par(family = "sans", cex = 1.5, col.main = "green4")
hist(
  slope_df$slope,
  col    = "green4",
  border = "green4",
  xlim   = c(-0.0001, 0.0001), # 显示范围为万分之一斜率
  main   = "Histogram of 1000 slopes"
)

结语

与理论一致，PC1 和 PC2 得分没有相关性。