使用 R 进行 PCA 特征选择

Question

我是一名生物学家。我的实验输出包含大量特征（存储为列数和 563 行）。列是8603个相当高的特征。

所以，当我尝试在 R 中进行 PCA 分析时，它会出现“内存不足”错误。

我也尝试过分段进行 princomp，但它似乎不适用于我们的 接近。

我尝试使用链接中给出的脚本...

http://www.r-bloggers.com/introduction-to-feature-selection-for-bioinformaticians-using-r-correlation-matrix-filters-pca-backward-selection/

但还是不行:(

我正在尝试使用以下代码

bumpus <- read.table("http://www.ndsu.nodak.edu/ndsu/doetkott/introsas/rawdata/bumpus.html", 
                     skip=20, nrows=49, 
                     col.names=c("id","total","alar","head","humerus","sternum"))

boxplot(bumpus, main="Boxplot of Bumpus' data") ## in this step it is showing the ERROR

# we first standardize the data:
bumpus.scaled <- data.frame( apply(bumpus,2,scale) )
boxplot(bumpus.scaled, main="Boxplot of standardized Bumpus' data")

pca.res <- prcomp(bumpus.scaled, retx=TRUE)
pca.res

# note:
# PC.1 is some kind of average of all the measurements 
#    => measure of size of the bird
# PC.2 has a negative weight for 'sternum' 
#    and positive weights for 'alar', 'head' and 'humerus'
#    => measure of shape of the bird

# first two principal components:
pca.res$x[,1:2]
plot(pca.res$x[,1:2], pch="", main="PC.1 and PC.2 for Bumpus' data (blue=survived, red=died)")
text(pca.res$x[,1:2], labels=c(1:49), col=c(rep("blue",21),rep("red",28)))
abline(v=0, lty=2)
abline(h=0, lty=2)

# compare to segment plot:
windows()
palette(rainbow(12, s = 0.6, v = 0.75)) 
stars(bumpus, labels=c(1:49), nrow=6, key.loc=c(20,-1), 
      main="Segment plot of Bumpus' data", draw.segment=TRUE) 

# compare to biplot:
windows()
biplot(pca.res, scale=0)
# what do the arrows mean?
# consider the arrow for sternum:
abline(0, pca.res$rotation[5,2]/pca.res$rotation[5,1])
# consider the arrow for head:
abline(0, pca.res$rotation[3,2]/pca.res$rotation[3,1])

但是第二行

boxplot(bumpus, main="Bumpus' data的Boxplot") ##显示错误

错误是

Error: cannot allocate vector of size 1.4 Mb

In addition: There were 27 warnings (use warnings() to see them)

请帮忙！

Answer 1

在功能数量巨大或超过 观察，建议根据以下公式计算主成分 转置数据集。在您的情况下尤其如此，因为默认 意味着计算一个 8603 x 8603 的协方差矩阵，它本身已经 消耗大约 500 MB 内存（哦，好吧，这不算太多，但是嘿......）。

假设矩阵X 的行对应于观察值 和列对应于特征，将您的数据居中，然后在 居中的X 的转置。特征对的数量不会超过 反正观察。最后，将每个结果特征向量乘以X^T。你做 不需要对特征值做后者（详细解释见下文）：

你想要什么

这段代码演示了PCA在转置数据集上的实现，并比较了prcomp和“转置PCA”的结果：

pca.reduced <- function(X, center=TRUE, retX=TRUE) {
  # Note that the data must first be centered on the *original* dimensions
  # because the centering of the 'transposed covariance' is meaningless for
  # the dataset. This is also why Sigma must be computed dependent on N
  # instead of simply using cov().
  if (center) {
    mu <- colMeans(X)
    X <- sweep(X, 2, mu, `-`)
  }
  # From now on we're looking at the transpose of X:
  Xt <- t(X)
  aux <- svd(Xt)
  V <- Xt %*% aux$v
  # Normalize the columns of V.
  V <- apply(V, 2, function(x) x / sqrt(sum(x^2)))
  # Done.
  list(X = if (retX) X %*% V else NULL,
       V = V,
       sd = aux$d / sqrt(nrow(X)-1),
       mean = if (center) mu else NULL)
}

# Example data (low-dimensional, but sufficient for this example):
X <- cbind(rnorm(1000), rnorm(1000) * 5, rnorm(1000) * 3)

original   <- prcomp(X, scale=FALSE)
transposed <- pca.reduced(X)

# See what happens:    
> print(original$sdev)
[1] 4.6468136 2.9240382 0.9681769
> print(transposed$sd)
[1] 4.6468136 2.9240382 0.9681769
> 
> print(original$rotation)
               PC1           PC2          PC3
[1,] -0.0055505001  0.0067322416  0.999961934
[2,] -0.9999845292 -0.0004024287 -0.005547916
[3,]  0.0003650635 -0.9999772572  0.006734371
> print(transposed$V)
              [,1]          [,2]         [,3]
[1,]  0.0055505001  0.0067322416 -0.999961934
[2,]  0.9999845292 -0.0004024287  0.005547916
[3,] -0.0003650635 -0.9999772572 -0.006734371

详情

要了解为什么可以处理转置矩阵，请考虑 以下：

特征值方程的一般形式是

          A x = λ x                               (1)

不失一般性，让M 成为您原件的居中“副本” 数据集X。用M^T M 替换A 产生

          M^T M x = λ x                           (2)

这个等式乘以M 得到

          M M^T M x = λ M x                       (3)

y = M x 的后续替换产生

          M M^T y = λ y                           (4)

已经可以看到y 对应于“协方差”的一个特征向量 转置数据集的矩阵（注意M M^T实际上不是真实的 作为数据集X 的协方差矩阵沿其列而不是其列居中 行。此外，缩放必须通过样本数（M 的行）来完成 而不是特征的数量（M 的列和M^T 的行）。

还可以看出M M^T和M^T M的特征值相同。

最后，最后一次乘以M^T 得到

          (M^T M) M^T y = λ M^T y                 (5)

其中M^T M 是原始协方差矩阵。

从等式（5）可以看出M^T y 是M^T M 的特征向量 特征值λ.