企业官网建设流程全解析

做微生物组、做植被调查、做水生态监测的朋友，大概率都经历过这种时刻：测序跑完，OTU表整理好，环境指标也测了一堆——pH、全氮、有机碳、降水量、温度、海拔，数据全在手里，却回答不了老板或审稿人那句最核心的问题：“到底是哪几个环境因子在驱动群落结构的变化？”排序图可以给出定性的方向，但对方要的是数字：你百分之多少的变化，由哪个因子解释。这时候需要用到的就是方差分解分析（Variance Partitioning Analysis，简称VPA）。

VPA不是新算法，它建立在约束排序（RDA、CCA）的基础上，把总方差切成几块，定量回答环境因子组的独立解释率、共同解释率和未解释率。无论你是做16S/ITS扩增子测序，还是做植物群落样方调查、浮游生物监测，只要手里有群落数据矩阵加环境变量矩阵，VPA都能派上用场。这篇内容我会从原理、R语言实操、结果解读到常见的坑，一步一步讲清楚，保证你拿到手就能跑，跑完能看懂，看得懂还敢写进文章里。

说明：文中代码基于R语言vegan包，个人在R 4.3.2 + vegan 2.6-4环境下实测通过。函数接口在不同大版本之间基本稳定，低版本用户注意检查vegan版本即可。

1. 项目概述：VPA到底在解决什么问题

1.1 一个真实的研究痛点

先还原一下我第一次做土壤微生物多样性项目时的场景。样品来自20个采样点，每点测了pH、水分、有机质、全氮等十来个指标，OTU表有几百万条序列。一开始我试图从热图和相关矩阵里找规律，发现指标之间互相纠缠：pH高的地方有机质也高，水分和海拔又是相关的。我根本说不清楚到底哪个因子起了主要作用，更别提解释比例了。

VPA要解决的正是这个问题。它不把每个环境因子单独拿出来比相关性，而是把环境变量分成若干组，比如气候组、土壤组、空间变量组，然后通过约束排序模型，把群落总方差的来源拆开：土壤组的独立贡献、气候组的独立贡献、两组重叠的共同贡献，以及模型无法解释的残差。输出是一张方差分配表和一个Venn图，谁贡献大谁贡献小一目了然。

1.2 VPA适用的典型场景

这套方法适用面很广，核心条件只有一个：你有响应变量矩阵（通常是物种多度或出现/不出现数据）以及解释变量矩阵（可以是连续性或分类环境因子）。常见场景包括：

• 扩增子测序项目：评估pH、土壤养分、气候等对细菌/真菌群落结构变化的解释程度。
• 植物生态学：量化地形、土壤、干扰史对植被样方组成的相对贡献。
• 淡水与海洋生态：比较水体理化参数、空间距离和生物相互作用对浮游/底栖群落的影响。
• 宏观生态与生物地理：区分环境筛选（niche）和空间过程（如扩散限制）的相对作用。

在所有这些场景里，VPA的核心价值是“定量归因”。认知上我们可以讨论很多生态过程机制，但落到文章里，审稿人希望看到量化的数字：环境解释了多少、空间解释了多少、残差还有多少，这直接决定了你的结论站不站得住脚。

2. 方法原理：VPA的数学本质与关键设计

2.1 它不是一个独立的算法，而是一次“组合式拆解”

很多人第一次听说VPA，以为它是一个单独打包好的统计算法。其实VPA是建立在约束排序之上的方差拆解策略。最底层的基础是RDA（冗余分析）或CCA（典范对应分析）。你可以把RDA理解成“多元线性回归的矩阵版本”：响应矩阵Y（物种）×解释矩阵X（环境因子），模型算出被X解释的方差占总方差的比例，这个比例就是R²。VPA只是把这个R²的计算逻辑进一步细化，用“全模型—子模型”的相减，把不同解释变量组各自贡献的方差剥离出来。

这背后的数学逻辑并不复杂。有两组变量X1和X2时，需要计算四件事：全模型（X1+X2）的总解释率、仅X1模型的解释率、仅X2模型的解释率，然后通过嵌套模型的减法得到各分块。全模型的总解释率中包含了两组变量单独和重叠的所有贡献，而单组模型的解释率里则天然混入了另一组变量能解释的公共部分，所以必须相减才能剥离出所谓“纯效应”：

1. 纯X1贡献 = E(X1+X2) − E(X2)；
2. 纯X2贡献 = E(X1+X2) − E(X1)；
3. 共同贡献 = E(X1) + E(X2) − E(X1+X2)；
4. 残差 = 1 − E(X1+X2)。

这四步就是VPA的全部核心。用Venn图表示，你会看到两个重叠的圆（X1、X2），中间重叠的“眼睛”区域是共同贡献，两侧月牙是纯贡献，外圈空白是残差。注意共同贡献说明的是“两组变量无法区分的重合解释部分”，并不代表一个独立的生态机制，这一点后面我会重点展开。

2.2 为什么大家都用调整R²而不是原始R²

使用原始R²会有一个严重问题：解释变量个数越多，R²天然越高，哪怕这些变量全是随机噪声。你可能听说过这个现象，放到多元回归里叫“过拟合”，放到RDA里也一样：往模型里加一个没意义的变量，总能多解释一点点方差。VPA的结果里，各组纯效应和共同效应都是通过模型相减得到的，如果直接用原始R²，变量多的组必然虚高，比较就没意义了。

所以vegan的varpart()同时返回Raw和Adjusted两套结果，你把数字写进文章时，一定要用Adjusted那一列。调整R²（由Legendre和Anderson提出，常写作R²adj）会按自由度对解释量打折扣，变量每多一个、样本每少一个，折扣就越大。我和同行交流的共识是：文章里凡是出现VPA结果，一律写调整后的解释百分比，否则审稿人大概率会打回来追问一句“用的是raw还是adjusted”。

2.3 线性RDA，还是单峰CCA，怎么选

VPA底下是RDA还是CCA，取决于你对物种-环境响应关系的假设。RDA假设物种多度沿环境梯度线性变化，适合数据梯度短、物种关系以线性为主的情况，比如很多土壤细菌群落和环境pH的关系，在几个pH单位跨度内基本可以按线性处理。CCA假设单峰响应，即每个物种在环境梯度的某个最适点达到最多，适合梯度长、有明显生态位分化的情况，典型如山地植物群落沿海拔梯度的分布。

实操中我自己的倾向：先用DCA（去趋势对应分析）看一下第一轴的梯度长度。梯度长度低于3，用RDA；高于4，用CCA；3到4之间两者都行，但为了稳妥我一般选CCA。这个规矩来自《Numerical Ecology》的建议，已经用了很多年，新手直接照搬即可。另外需要注意，vegan的varpart()默认走RDA路径，做CCA版本的方差分解需要手动基于cca模型写循环，非必要不建议新手折腾。

3. 实操指南：基于R的VPA完整流程

3.1 数据准备：三张表的规范格式

跑VPA之前，先把数据整理成三张表。第一张是群落数据表（物种表），行是样本，列是物种/OTU，值是多度、丰度或0/1数据；第二张是环境变量表，行与物种表完全一致，列是各个环境因子；第三张是空间变量表（可选，但强烈建议准备），行同样是样本，列是采样点坐标或由坐标生成的空间特征向量。

有一个最容易被新手忽略的坑：行名必须是一一对应的。我见过太多次OTU表和环境表读进来后发现匹配不上，原因只是样本ID的格式不一致，比如一个用“S01”，一个用“Sample_01”。读取数据后第一件该做的事就是核对行名，用identical(rownames(otu), rownames(env))检查，返回FALSE就赶紧统一格式。这一步不做好，后面所有结果都是空谈。

群落数据建议做Hellinger转化。因为OTU表经常以绝对丰度呈现，包含大量零值和巨大多度值，直接丢进RDA会严重受限于“物种总丰度”造成的虚假关联。Hellinger转化把绝对丰度转成相对丰度后再开平方，能把样本间的欧氏距离和生态学上常用的Bray-Curtis距离拉近，是当前群落约束排序的主流预处理。环境变量则做标准化（均值0、标准差1），让量纲不同的指标（pH是0-14，全氮可能是mg/kg上千）在模型里公平竞争。

library(vegan) # 读取数据（示例，实际请按自己的文件调整） otu <- read.delim("otu_table.txt", row.names = 1, check.names = FALSE) env <- read.delim("env.txt", row.names = 1, check.names = FALSE) coords <- read.delim("coordinates.txt", row.names = 1) # 行名核对 identical(rownames(otu), rownames(env)) # 应为 TRUE # 群落数据Hellinger转化 otu_hel <- decostand(otu, method = "hellinger") # 环境数据标准化 env_std <- decostand(env, method = "standardize")

3.2 环境变量预处理：必须做，不能省

把环境变量一股脑全塞进VPA是我见过的第二大坑。我当年就干过这事，把11个土壤指标全放到一组，结果调整R²低得可怜，还因为变量间共线性导致结果完全没法看。后来才知道，约束排序的多重共线性问题和回归一样严重，VIF可以轻松飙到几十。

推荐的流程分两步。第一步，用方差膨胀因子（VIF）过滤原变量，一般认为VIF大于10的变量需要剔除，或者两两相关性高于0.8的只保留一个。第二步，用前向选择（forward selection）挑选在解释群落变化中贡献显著的变量，保证每组最终只保留3~6个有效变量。这一步在vegan里用ordiR2step()，它基于调整R²做每一步判断，并且用置换检验控制显著性：

# 先把环境变量框拆成两组，例如气候组和土壤组 clim <- env_std[, c("temp", "precip", "humidity")] soil <- env_std[, c("pH", "TN", "SOC", "moisture")] # 土壤组的前向选择示例 mod0 <- rda(otu_hel ~ 1, data = soil) # 空模型 mod1 <- rda(otu_hel ~ ., data = soil) # 全模型 sel <- ordiR2step(mod0, scope = mod1, perm.max = 999) # 前向选择 sel$anova # 查看每一步的显著性 # 提取被保留的变量名 kept_vars <- labels(terms(sel)) kept_vars

前向选择的结果就是进入正式VPA的变量集合。写文章的时候，方法部分通常要写清楚：环境变量经VIF筛选后，每组保留哪些变量，随后进行前向选择，置换次数是多少。这一句话虽然朴素，却是审稿人判断你结果可信度的重要依据，别省略。

3.3 vegan::varpart 核心代码与绘图

变量筛选完之后，VPA本体其实只有两三行代码。把两组或三组解释变量矩阵准备好，直接调用varpart()：

# 两因子组VPA：气候组 vs 土壤组 vp12 <- varpart(otu_hel, clim_selected, soil_selected) vp12 plot(vp12, bg = c("steelblue", "orange"), Xnames = c("气候因子", "土壤因子")) # 三因子组VPA：气候组 vs 土壤组 vs 空间变量组 # 先构造空间变量，见3.4节 spa_selected # 前向选择后保留的空间变量矩阵 vp123 <- varpart(otu_hel, clim_selected, soil_selected, spa_selected) vp123 plot(vp123, bg = c("steelblue", "orange", "forestgreen"), Xnames = c("气候因子", "土壤因子", "空间变量"))

plot()会画出一个Venn图，但圆的大小并不严格按解释率比例缩放，所以阅读时以数字为准，图只是直观展示。打印的varpart结果包含Raw和Adjusted两行，Adjusted行下面就是拆好的各分块数字。需要提醒的是，如果解释变量矩阵的变量个数接近甚至超过样本量，自由度会不够，调整R²可能直接变成负数，说明模型已经严重过拟合，结果基本没法用。

另外一个老生常谈但必须做的事：显著性检验。varpart()本身不返回p值，要检验每个分数是否显著，得对对应的偏RDA模型做置换检验：

# 检验全模型的显著性 anova(rda(otu_hel, cbind(clim_selected, soil_selected)), perm.max = 999) # 检验气候组纯效应的显著性（以土壤组为协变量） anova(rda(otu_hel, clim_selected, soil_selected), perm.max = 999) # 检验土壤组纯效应的显著性（以气候组为协变量） anova(rda(otu_hel, soil_selected, clim_selected), perm.max = 999)

关于置换次数，投稿级别的工作我习惯用9999次，既能保证稳定，又不至于慢到等凉一杯咖啡。别用默认的199次，审稿人确实会挑这个细节。

3.4 空间变量怎么构造

如果你研究的是跨区域采样（比如沿着一条山脉、一条河流采样），生物群落天然有空间距离带来的相似性衰减：距离近的样点，就算环境因子完全一样，群落也会更相似。这种空间结构如果不从环境效应里剥离，环境因子组的解释率就会被系统性高估。所以建议把空间变量作为一组参与VPA。

构造空间变量的通行方法有两个。简单的办法是用采样点坐标生成多项式项，直接放一次项、二次项和交叉项：

coords <- read.delim("coordinates.txt", row.names = 1) spa_mat <- cbind(coords$x, coords$y, coords$x * coords$y, coords$x^2, coords$y^2) colnames(spa_mat) <- c("x", "y", "xy", "x2", "y2") spa_std <- decostand(spa_mat, method = "standardize")

如果采样点较多（超过30个），我更推荐使用vegan里的pcnm()生成PCNM特征向量，它能捕捉更精细的多尺度空间结构，结果比直接放多项式稳定。不过特征向量可能生成很多，同样需要经过前向选择，只保留显著的前几根作为spa_selected。

4. 结果解读：VPA图表与数字的深层逻辑

4.1 读懂Venn图中的数字

varpart()的打印结果里会用字母标出各分块。以三组变量为例，输出包括[a]、[b]、[c]（各组纯效应），[d]、[e]、[f]（两两共同），[g]（三组共同），[h]（残差）。我把各分块的含义整理成一个对照表，方便你贴进自己的实验记录本：

举个例子。假设结果里气候纯效应是9%，土壤纯效应是21%，共同解释率是6%，残差是64%。那么可以写结论：土壤性质对群落变化的独立解释力显著大于气候因子，两者共同影响仅占6%，说明土壤pH和营养等核心指标可能是直接驱动因素，气候更多通过影响土壤性质间接起作用。完整写法是：环境与空间因子共解释了36%的群落变异，其中土壤独立解释21%，气候独立解释9%，残差占64%，表明还存在大量未测环境因子或随机过程的影响。这个写法既完整又不夸大。

4.2 共享分数：最大方也最容易翻车的数字

共同解释率是最容易被过度解读的区域。很多人看到气候和土壤有10%的共同解释，就写“气候和土壤存在强烈的交互作用”。这是非常危险的说法。共同解释率本质上来源于解释变量之间的相关性（共线性），它只说明这部分变异的归属在两组变量之间无法区分，并不等同于生态过程中的真正联合效应。

打个生活化类比：一个人中午吃了火锅又喝了奶茶，夜里胃不舒服，你能分清是火锅的“独立责任”还是奶茶的“独立责任”？那一部分没法归因的不舒服程度，就是“共同解释率”。你可以说这个症状与两者都有关联，但不能确证哪个是主犯。所以写论文时的标准姿势是：报告独立解释率，指出共同解释率的存在，并说明由于环境因子本身的相关性，这部分变异无法被唯一归因。想进一步拆解谁才是主犯，需要更精细的实验设计或梯度分析，这不是VPA的局限，而是观察性数据共线性天然带来的问题。

4.3 残差很大，不是模型失败，而是一个重要信号

刚开始跑VPA的人看到残差经常心里一凉：三组变量加起来只解释了不到三成，剩下七成都不知道为什么。这是生态学研究的常态，别慌。微生物群落尤其如此，解释率动辄只有个位数，残差六成以上稀松平常。残差大说明三件事：要么有重要的环境变量没测到（比如微量养分、生物间相互作用），要么群落本身就包含大量随机过程和中性漂变，再要么数据噪音偏高。

残差大反而值得写进讨论部分，它提示了后续研究方向。我在一篇底泥微生物的文章里就明确写过：VPA显示理化因子仅解释18.2%的群落变异，残差达七成以上，提示竞争、捕食和扩散过程可能在塑造群落中占据主要地位。这样的表述是审稿人乐于看到的——说明你不只是跑了一个模板，而是理解了结果背后的生态学含义。

5. 常见错误与避坑指南

5.1 用原始R²而不看调整R²

这是VPA最经典的一个错，犯错的人特别多，原因是vegan输出里Raw就在第一行，位置醒目，很多人顺手就抄了。我踩过一次后养成了习惯：任何模型比较、任何方差分配百分比，都只从Adjusted那行取数。如果你担心自己在结果表里抄错，可以在代码里直接提取调整后的分块数字，或者打印时加一行注释提醒自己用调整值，宁可多一步也不能抄错。

5.2 变量组塞太满，VIF不检查

如果一组变量里塞了十多个高度相关的指标，“纯效应”会被严重低估，“共同效应”虚高，整个结果看起来像是随机数。前面已经给了前向选择流程，这里再补充一个速查标准：每组变量控制在3~6个，任何变量的VIF不超过10，两两Pearson相关系数不超过0.8。超过这个红线就砍变量，换来的是更干净、更可重复的结果。砍完变量后记得重新算一遍VIF确认，因为前向选择后的变量组合VIF可能和初筛时不一样。

5.3 忽视空间变量导致环境解释虚高

在没有控制空间自相关时，环境组解释率常常虚高。想象一下采样点沿着一条河流从上游到下游，水温、pH、溶氧都和距离相关，最终群落也随距离渐变。如果模型里只有理化变量，它们会把空间变化的大部分功劳都“认领”走。所以只要采样范围跨了地理尺度，我建议无条件放一组空间变量进入VPA。有些审稿人就是盯着这一点，没放空间变量往往会收到“是否考虑空间自相关？”的尖锐问题，提前放进去能省一整个返修周期。

5.4 把置换检验当摆设

显著性检验的逻辑和线性回归的F检验一样，检验的是“解释率是否显著大于0”。很多人跑完varpart()、Venn图画完，收工写论文，完全没有做anova()。结果是文章里写“气候组显著解释X%”，但没有任何显著性证据支撑。正确的做法是对全模型和每个纯分数对应的偏RDA都做置换检验，并在方法部分写明检验方式、置换次数和错误率控制方法。我也习惯对共享分数做检验，虽然它难以解释，但至少可以说明它是否在统计上存在。

5.5 样本量太小还硬跑VPA

VPA的自由度消耗比普通回归大得多，样本量建议至少30个，组数越多要求越高。如果只有12个样点却拆分三组变量，调整R²很容易出现负数，这时候任何解释都没有意义。碰到样本量不足的项目，我的建议是退一步，改用单组的RDA加变量贡献排序，或者用Mantel检验做初步探索，宁可少讲归因，也比给审稿人递刀强。样本量这件事没有办法靠数据变换技巧补救，唯一的出路是补采样，没有捷径。

6. 扩展与进阶：VPA之外还能怎么用

6.1 RDA、CCA、db-RDA怎么选

VPA之外还有一个变体叫db-RDA（基于距离的冗余分析），它允许你用Bray-Curtis、Jaccard等任何相异度矩阵作为对象，再做约束排序。这在分析微生物群落时很有用，因为Bray-Curtis是生态学里最符合直觉的距离指标。R语言里做db-RDA的方差分解一般用adonis2()或capscale()配合方差分解循环。但就我个人的使用体感，经典RDA的varpart()在文章里依然是接受度最高的方案，db-RDA版本的VPA结果在比较时更容易被质疑，能用经典RDA就先用经典RDA。

6.2 和随机森林、Mantel检验的搭配

VPA擅长回答“变量组的解释比例”，但不擅长回答“单个变量的重要性排序”。想要后者，可以和随机森林回归搭配，把每个环境变量对群落整体属性的重要性排个序，再和VPA的组间结果互相印证。Mantel检验则适合做初步筛选：跑VPA之前，用mantel()批量检验每个环境因子与群落距离的相关性，能帮你在前向选择之前先摸个底。组合打法常见于高分论文：先用Mantel检验找候选因子，再用前向选择定变量，最后用VPA给组间解释率定案。

6.3 三组、四组变量的VPA怎么玩

vegan::varpart()支持最多四组变量。四组时Venn图变成四圆交叠，输出分块会飙升到15个，肉眼已经很难直接解读，通常还是看数字表格。实际操作中，最常见的是三组组合：环境因子一组、空间变量一组、干扰或土地利用归到第三组。超过四组我并不推荐，因为每个分块的自由度被切得太碎，估计会很不可靠。只保留最核心的两到三组，结论往往更清楚，也更容易被读者理解。

最后再分享一个我自己的习惯：VPA做完之后，我会顺手把调整R²的三张关键数字（纯效应、共同效应、残差）写进实验记录本，并附上一句“为什么”的备注。第二天再回来看，往往能挑出昨天忽略的问题——比如某个共同解释率高得可疑，就回去检查VIF。这种“隔夜复核”听着土，但确实帮我躲过至少两次返工，你也可以试试。