做微生物组、做植被调查、做水生态监测的朋友大概率都经历过这种时刻测序跑完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时需要计算四件事全模型X1X2的总解释率、仅X1模型的解释率、仅X2模型的解释率然后通过嵌套模型的减法得到各分块。全模型的总解释率中包含了两组变量单独和重叠的所有贡献而单组模型的解释率里则天然混入了另一组变量能解释的公共部分所以必须相减才能剥离出所谓“纯效应”纯X1贡献 E(X1X2) − E(X2)纯X2贡献 E(X1X2) − E(X1)共同贡献 E(X1) E(X2) − E(X1X2)残差 1 − E(X1X2)。这四步就是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用CCA3到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]残差。我把各分块的含义整理成一个对照表方便你贴进自己的实验记录本输出标签含义解读建议[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 样本量太小还硬跑VPAVPA的自由度消耗比普通回归大得多样本量建议至少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。这种“隔夜复核”听着土但确实帮我躲过至少两次返工你也可以试试。