0

我有一个包含多个物种和大约 400 个变量的数据集。我想对每个物种执行主成分分析 (PCA)，并返回每个物种具有最高负载值的变量。

制作我的数据的复制虚拟集：

set.seed(45)
pcadata <- data.frame(matrix(sample(10, 26746*400, TRUE), ncol=400))
cbind(pcadata,"Species")

我遇到的一个问题是给定物种的样本量不同。例如，我可能有 250 个物种 A 的样本和 520 个物种 B 的样本。因此我必须使用该prcomp函数，因为我的变量比样本多。因此，如果 Species A (spA) 在 data.frame 中，我首先必须对数据进行子集化：

pcadata.s<-pcadata[,2:401]

pca<-prcomp(pcadata.s,cor=T,scale=T)
al<-abs(pca$rotation)                    #Absolute value of the loading value
loads<-sweep(al,2,colSums(al),"/")       #Percentage contribution
loads.mtx<-as.data.frame(loads)
rownames(loads.mtx)[apply(loads.mtx,2,which.max)] #Return the Column-name with the max value

我想，不必每次都进行子采样，而是获取每个 Species 分组的列名，例如：

Species  PC1     PC2      PC3      PC4      PC5
 spA      V3     V100     V287     V2       V65
 spB      V78    V197     V310     V23      V333 
........

刚刚意识到：我需要选择我感兴趣的组件，最好是解释方差的 95%，也许我也会尝试 99% ......但我必须包含代码。

任何建议将不胜感激。

4

2 回答 2

2

如果我正确理解您的问题，您希望将该prcomp函数应用于数据的子集。没有本地的方式来处理这个（据我所知）。

您可以尝试以下方法：

species <- unique(colnames(pcadata))
pcaresults <- list()
for (sp in species) {
  spIndices <- which(colnames(pcadata) == sp)
  pcaresults[sp] <- prcomp(pcadata[,spIndices], cor=T,scale=T)
}

这将为您提供一个列表，其中每个元素都是该物种的 PCA 的返回结果。您可以更改循环或格式化返回列表，以仅获取您想要的数据。

于 2013-09-14T07:54:01.683 回答
2

如果您将物种名称作为变量保留在数据框中，则可以ddply在plyr包中使用。

library(plyr)
# create data with a species variable
set.seed(45)
df <- data.frame(matrix(sample(1:10, size = 50, replace = TRUE), ncol = 5))
df$species <- rep(1:2, each = 5)

# run pca and massage data per species
df2 <- ddply(.data = df, .variables = .(species), function(x){
  pca <- prcomp(x[ , 1:5], cor = TRUE, scale = TRUE)
  load <- abs(pca$rotation) 
  prop_load <- apply(load, 2, function(x) x/sum(x))
  max_load <- rownames(prop_load)[apply(prop_load, 2, function(x) which.max(x))]
  max_load2 <- data.frame(t(max_load))
  names(max_load2) <- colnames(load)
  return(max_load2)
}
)
df2

# species PC1 PC2 PC3 PC4 PC5
# 1       1  X1  X2  X4  X3  X5
# 2       2  X2  X1  X3  X2  X5
于 2013-09-14T09:53:50.917 回答