我有不同的变量在人群中相互作用。基本上我一直在做千足虫的清单并测量地形的其他一些值，比如：
- 采集标本的种类和数量
- 动物所在的不同环境
- 酸碱度
- 有机物百分比
- P、K、Mg、Ca、Mn、Fe、Zn、Cu的含量
- Ca + Mg/K 关系
基本上我想使用 PCA 来确定哪些变量驱动样本的可变性并使森林（环境）不同；我应该将哪些变量用于“变量”，哪些变量用于“个人”？
使用什么标准将变量分离为解释变量和生态学排序方法的响应？
机器算法验证
主成分分析
多元分析
生态
对应分析
2022-03-14 07:50:32
1个回答
正如@amoeba 在评论中提到的那样，PCA 只会查看一组数据，它将向您展示这些变量的主要（线性）变化模式、这些变量之间的相关性或协方差，以及样本之间的关系（行) 在您的数据集中。
通常对物种数据集和一组潜在解释变量所做的事情是适应受约束的排序。在 PCA 中，主成分，即 PCA 双图上的轴，是作为所有变量的最佳线性组合推导出来的。 、TotalCarbon的土壤化学数据集上运行此程序，您可能会发现第一个成分是
和第二个组件
这些组件可以从测量的变量中自由选择，并且被选择的是那些依次解释数据集中最大变化量的变量，并且每个线性组合都是正交的（与其他组合不相关）。
在一个受约束的排序中，我们有两个数据集，但我们不能自由地选择我们想要的第一个数据集（上面的土壤化学数据）的任何线性组合。相反，我们必须选择第二个数据集中变量的线性组合，以最好地解释第一个数据集中的变化。此外，在 PCA 的情况下，一个数据集是响应矩阵并且没有预测变量（您可以将响应视为预测自身）。在受约束的情况下，我们有一个响应数据集，我们希望用一组解释变量来解释它。
尽管您没有解释哪些变量是响应，但通常人们希望使用环境解释变量来解释这些物种的丰度或组成（即响应）的变化。
PCA 的约束版本是生态圈中称为冗余分析（RDA）的东西。这假设了物种的潜在线性响应模型，这要么不合适，要么仅在物种响应的梯度较短时才合适。
PCA 的替代方法是一种称为对应分析 (CA) 的东西。这是不受约束的，但它确实有一个潜在的单峰响应模型，这在物种如何沿较长梯度响应方面更加现实。另请注意，CA 模拟相对丰度或组成，PCA 模拟原始丰度。
CA 有一个约束版本，称为约束或规范对应分析(CCA) - 不要与称为规范相关分析的更正式的统计模型相混淆。
在 RDA 和 CCA 中，目的是将物种丰度或组成的变化建模为解释变量的一系列线性组合。
从描述中听起来你的妻子想根据测量的其他变量来解释千足虫物种组成（或丰度）的变化。
一些警告的话；RDA 和 CCA 只是多元回归；CCA 只是一个加权多元回归。你学到的关于回归的任何东西都适用，还有一些其他的陷阱：
- 当您增加解释变量的数量时，约束实际上变得越来越少，并且您并没有真正提取出最佳解释物种组成的组件/轴，并且
- 使用 CCA，当您增加解释因素的数量时，您就有可能将曲线的伪影引入 CCA 图中的点配置中。
- 与更正式的统计方法相比，RDA 和 CCA 的理论基础欠发达。我们只能合理地选择继续使用逐步选择的解释变量（由于我们不喜欢将其作为回归中的选择方法的所有原因，这并不理想），我们必须使用置换测试来做到这一点。
所以我的建议与回归相同；提前考虑你的假设是什么，并包括反映这些假设的变量。不要只是将所有解释变量都混在一起。
例子
不受约束的排序
主成分分析
我将展示一个使用我帮助维护的 R 的vegan包比较 PCA、CA 和 CCA 的示例，该包旨在适合这些排序方法：
library("vegan") # load the package
data(varespec) # load example data
## PCA
pcfit <- rda(varespec)
## could add `scale = TRUE` if variables in different units
pcfit
> pcfit
Call: rda(X = varespec)
Inertia Rank
Total 1826
Unconstrained 1826 23
Inertia is variance
Eigenvalues for unconstrained axes:
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
983.0 464.3 132.3 73.9 48.4 37.0 25.7 19.7
(Showed only 8 of all 23 unconstrained eigenvalues)
与 Canoco 不同，素食主义者不标准化惯性，因此总方差为 1826，特征值采用相同的单位，总和为 1826
> cumsum(eigenvals(pcfit))
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
982.9788 1447.2829 1579.5334 1653.4670 1701.8853 1738.8947 1764.6209 1784.3265
PC9 PC10 PC11 PC12 PC13 PC14 PC15 PC16
1796.6007 1807.0361 1816.3869 1819.1853 1821.5128 1822.9045 1824.1103 1824.9250
PC17 PC18 PC19 PC20 PC21 PC22 PC23
1825.2563 1825.4429 1825.5495 1825.6131 1825.6383 1825.6548 1825.6594
我们还看到第一个特征值大约是方差的一半，而前两个轴我们已经解释了总方差的约 80%
> head(cumsum(eigenvals(pcfit)) / pcfit$tot.chi)
PC1 PC2 PC3 PC4 PC5 PC6
0.5384240 0.7927453 0.8651851 0.9056821 0.9322031 0.9524749
可以从样本和物种在前两个主成分上的得分绘制双图
> plot(pcfit)
这里有两个问题
- 排序基本上由三个物种主导——这些物种离起源最远——因为它们是数据集中最丰富的分类群
- 排序中有一个强烈的曲线拱形，表明一个长的或占主导地位的单一梯度已被分解为两个主要的主要成分，以保持排序的度量属性。
加州
由于单峰响应模型，CA 可能会在这两个方面提供帮助，因为它可以更好地处理长梯度，并且它可以模拟物种的相对组成而不是原始丰度。
执行此操作的vegan / R代码类似于上面使用的 PCA 代码
cafit <- cca(varespec)
cafit
> cafit <- cca(varespec)
> cafit
Call: cca(X = varespec)
Inertia Rank
Total 2.083
Unconstrained 2.083 23
Inertia is mean squared contingency coefficient
Eigenvalues for unconstrained axes:
CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
0.5249 0.3568 0.2344 0.1955 0.1776 0.1216 0.1155 0.0889
(Showed only 8 of all 23 unconstrained eigenvalues)
在这里，我们解释了大约 40% 的站点之间相对组成的变化
> head(cumsum(eigenvals(cafit)) / cafit$tot.chi)
CA1 CA2 CA3 CA4 CA5 CA6
0.2519837 0.4232578 0.5357951 0.6296236 0.7148866 0.7732393
物种和地点分数的联合图现在较少由少数物种主导
> plot(cafit)
您选择哪一个 PCA 或 CA 应由您希望对数据提出的问题决定。通常对于物种数据，我们更经常对物种套件的差异感兴趣，因此 CA 是一种流行的选择。如果我们有环境变量的数据集，比如水或土壤化学，我们不会期望它们以单峰方式沿梯度响应，因此 CA 将是不合适的，而 PCA（相关矩阵，scale = TRUE在rda()调用中使用）将是更合适。
受约束的协调；CCA
现在，如果我们希望使用第二组数据来解释第一个物种数据集中的模式，我们必须使用受约束的排序。通常这里的选择是 CCA，但 RDA 是一种替代方案，RDA 也是在数据转换后使其能够更好地处理物种数据。
data(varechem) # load explanatory example data
我们重复使用该cca()函数，但我们要么提供两个数据框（X用于物种，以及Y用于解释/预测变量），要么提供一个模型公式，列出我们希望拟合的模型的形式。
要包含所有变量，我们可以将varechem ~ ., data = varechem其用作包含所有变量的公式——但正如我上面所说，这通常不是一个好主意
ccafit <- cca(varespec ~ ., data = varechem)
> ccafit
Call: cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe + Mn +
Zn + Mo + Baresoil + Humdepth + pH, data = varechem)
Inertia Proportion Rank
Total 2.0832 1.0000
Constrained 1.4415 0.6920 14
Unconstrained 0.6417 0.3080 9
Inertia is mean squared contingency coefficient
Eigenvalues for constrained axes:
CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7 CCA8 CCA9 CCA10 CCA11
0.4389 0.2918 0.1628 0.1421 0.1180 0.0890 0.0703 0.0584 0.0311 0.0133 0.0084
CCA12 CCA13 CCA14
0.0065 0.0062 0.0047
Eigenvalues for unconstrained axes:
CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8 CA9
0.19776 0.14193 0.10117 0.07079 0.05330 0.03330 0.01887 0.01510 0.00949
上述排序的三元组是使用该plot()方法产生的
> plot(ccafit)
当然，现在的任务是找出其中哪些变量实际上是重要的。另请注意，我们仅使用 13 个变量解释了大约 2/3 的物种差异。在此排序中使用所有变量的问题之一是我们在样本和物种分数中创建了一个拱形配置，这纯粹是使用太多相关变量的人工制品。
如果您想了解更多相关信息，请查看vegan文档或一本关于多元生态数据分析的好书。
与回归的关系
用 RDA 来说明联系是最简单的，但 CCA 是一样的，除了所有内容都涉及行和列双向表边际总和作为权重。
从本质上讲，RDA 相当于将 PCA 应用于拟合值矩阵，该矩阵来自拟合到每个物种（响应）值（例如丰度）的多元线性回归，预测变量由解释变量矩阵给出。
在 R 中，我们可以这样做
## centre the responses
spp <- scale(data.matrix(varespec), center = TRUE, scale = FALSE)
## ...and the predictors
env <- as.data.frame(scale(varechem, center = TRUE, scale = FALSE))
## fit a linear model to each column (species) in spp.
## Suppress intercept as we've centred everything
fit <- lm(spp ~ . - 1, data = env)
## Collect fitted values for each species and do a PCA of that
## matrix
pclmfit <- prcomp(fitted(fit))
这两种方法的特征值相等：
> (eig1 <- unclass(unname(eigenvals(pclmfit)[1:14])))
[1] 820.1042107 399.2847431 102.5616781 47.6316940 26.8382218 24.0480875
[7] 19.0643756 10.1669954 4.4287860 2.2720357 1.5353257 0.9255277
[13] 0.7155102 0.3118612
> (eig2 <- unclass(unname(eigenvals(rdafit, constrained = TRUE))))
[1] 820.1042107 399.2847431 102.5616781 47.6316940 26.8382218 24.0480875
[7] 19.0643756 10.1669954 4.4287860 2.2720357 1.5353257 0.9255277
[13] 0.7155102 0.3118612
> all.equal(eig1, eig2)
[1] TRUE
出于某种原因，我无法使轴分数（负载）匹配，但这些分数总是被缩放（或不被缩放），所以我需要看看这些是如何在这里完成的。
我们没有rda()像我展示的那样通过lm()等方式进行 RDA，但我们对线性模型部分使用 QR 分解，然后对 PCA 部分使用 SVD。但基本步骤是相同的。