首先,请注意您对 PCA 的理解略有偏差。PCA 不会将变量“分组”到主成分中。相反,每个主要成分都是一个新变量(“新代谢物”),其长度与您的每个原始变量(在您的情况下为 16)相同。
的确,每个主成分都可以代表一组原始代谢物,但它比这更复杂和复杂。每种代谢物实际上可以与多个主要成分相关联,每个主要成分的量不同。一组代谢物完全由主成分 A 表示,而另一组完全由主成分 B 表示,这可能是真的,但更有可能的是,每个代谢物都与多个 PC 相关,每个 PC 的程度不同。
因此,将每台 PC 视为一种新的代谢物(用 PCA 的说法,我们称之为“本征代谢物”),它代表了整个代谢物的一般模式。一些代谢物高度代表这种模式,而另一些则不是,有些代谢物具有多种模式的某些特性。
现在我们来回答你的问题。我们引入响应变量的方法是查看每个主成分(每个特征代谢物)和响应变量之间的关系——例如,相关性或线性回归。这是 R 代码中的一个简洁的小示例。(请注意,我使用的是svd
函数而不是pls
包,正如您可能正在使用的那样)。
set.seed(1337)
center.rows = function(m) t(apply(m, 1, function(r) r - mean(r)))
response.var = rnorm(16, 0, 1)
relevant.metabolites = t(replicate(50, rnorm(16, response.var, 2)))
unrelated.metabolites = t(replicate(110, rnorm(16, 0, 2)))
metabolite.data = rbind(relevant.metabolites, unrelated.metabolites)
s = svd(center.rows(metabolite.data))
我所做的是创建一个模拟代谢物数据矩阵,其中包含 50 个与响应变量有些相关的基因,以及 110 个不相关的基因。我们可以看到只有一台 PC 是显着的,它解释了大约 33% 的方差:
var.explained = s$d^2 / sum(s$d^2)
plot(var.explained)
现在让我们将该 PC 与响应变量进行比较:
plot(response.var, s$v[, 1])
如您所见,第一台 PC 几乎完美地重新捕获了我们用来创建矩阵的响应变量(相关性为 -.993),尽管我们在其中添加了很多额外的噪声。(不要担心相关性是负的- 方向是任意的)。实际上,这种与“特征代谢物”的相关性远大于响应变量与任何单个重要基因的相关性:
hist(apply(relevant.metabolites, 1, function(r) cor(r, response.var)))
当然,这正是我们首先应用 PCA 的原因。虽然响应变量和任何单个代谢物之间的关系可能太弱而无法测量其余的噪音,但当您将许多基因组合在一起时,您可能会发现主成分很好地捕捉到了该响应!