资料来源：https://www.freesion.com/article/2497946072/

这篇微文写的非常棒，因此转发至自己的智库中，以便日后学习。

若有侵权，请联系我。

衷心的感谢每一位知识的奉献者。。。。

VPA，全称Variance Partitioning Analysis，中文成为方差分解分析，该分析的目的是确定指定的环境因子对群落结构变化的解释比例。

我们使用CCA/RDA的排序分析方法可以得到所有参与分析的环境因子对群落变化的解释比例。

那么在进行VPA时，首先就要对这些环境因子进行一个分类，然后在约束其它类环境因子的情况下，对某一类环境因子进行排序分析，这种分析也成为偏分析，即partial CCA/RDA。

在对每一类环境因子均进行偏分析之后，即可计算出每一个环境因子单独以及不同环境因子相互作用分别对生物群落变化的贡献。

分析实战

这里使用R语言vegan包的varpart()函数进行VPA分析，之后使用plot函数对结果进行可视化。

VPA是确定不同类型环境因子对群落变化的解释，那么首先就要对环境因子进行一个分类，这个类怎么分呢？

简单的说就是你自己想怎么分就怎么分，根据你研究的实际情况自己确定怎么分类。

两种环境因子分类

当分析的环境因子只有两类时，可以将两类环境因子放在不同的数据框中进行分析。

首先我们导入示例数据。

data(mite)

data(mite.env)

data(mite.pcnm)

这里mite为群落丰度数据表格，行为样本，列为物种；mite.env和mite.pcnm分别为两个环境因子的数据表格，同样行为样本，列为环境因子。

mod <- varpart(mite, mite.env, mite.pcnm, transfo="hel")

mod

进行VPA的时候，第一个数据框为群落数据，之后两个数据框分别代表两类环境因子，transfo时对数据进行转换，hel为hellinger转换，可以避免分析的“弓形效应”。

在结果中我们看Individual fractions部分即可。

a为X1也就是命令中第二个数据框单独对群落变化的贡献。

c为X2也就是命令中第三个数据框单独对群落变化的贡献。

b为X1和X2的相互作用对群落变化的贡献。

d为X1和X2无法解释的群落变化。

使用plot()函数对结果进行可视化。

plot(mod, bg = c("hotpink","skyblue"))

三种及以上的分类

对与将环境因子分为3类或4类的情况，可以将环境因子放在同一个数据框中，之后使用formula的形式指定不同的分类因子。

mod <- varpart(mite, ~ SubsDens + WatrCont, ~ Substrate + Shrub + Topo,

mite.pcnm, data=mite.env, transfo="hel")

可以看到3种分类的结果就相对复杂一些，a-h的含义需要根据Partition table中的对应情况确定一下，这个我就不从头捋一遍了，感兴趣的朋友可以自行画个交集和并集的图分解一下。

为什么不捋一遍结果呢，是因为可视化之后解释的比例就直接给出了????

plot(mod, bg=2:4)

VPA最多能将环境因子分为4组，再多就不行了，不过我想也几乎不会遇到能够将环境因子分为很多组的情况。

下面是干货

做VPA分析对样本的数量有一定的要求，记得最开始使用Canoco分析的时候，如果样本数小于环境因子数目减2，软件就会报错。

也就是说样本数目至少要比因子数目多2个，不要问我为什么，我也不知道为什么，没研究过具体的算法，我只是只要如果环境因子比样本数还多的话，就算做完了，结果也很奇葩，根本没法解释。

这对于一些大样本量的研究项目当然不成问题，但是对于一些经费有效的研究，可能样本数目就会是使用VPA的一个限制因素。

这里我有一个变通的方法，就是先对不同分类环境因子做降维分析，比如说PCoA，之后使用主要的PC的结果替代环境因子，从而达到降低实际使用因子数目的目的。

但是具体使用前几个PC就要不断的尝试，根据结果进行调整了。

对于结果的可视化，R默认的结果图确实不是很好看，关键还不太好调整，大家可以先默认出一个图，然后把每一部分的解释比例记下来，之后使用在线的Venn图绘制工具，画一个自己满意的只有圆圈的Venn图，再手动把解释比例和环境因子的名称给P上去