我去年分析一组土壤微生物数据时碰到一个特别常见的场景冗余分析RDA跑出来pH、含水量、总氮几个箭头指向清清楚楚p值也都小于0.001看起来结果很完美。可到了写文章的时候审稿人问了一句总变异中环境因子到底解释了多少比例空间距离又解释了多少它们之间重叠的部分占多少我当场愣住了。RDA只告诉我这些因子显著但并没有告诉我每个部分贡献了多少解释比例。于是我去翻方法学文献找到了方差分解分析VPAVariation Partitioning Analysis这个工具它刚好解决这个问题。这篇博文就把我从数据准备到结果解读的完整流程、踩过的坑以及为什么VPA是群落生态学中归因分析标配的方法一次讲清楚。无论你是做微生物群落、浮游生物还是植物生态这个方法大概率能用在你的数据上。1. 单跑RDA的局限你只有显著没有比例很多刚接触群落数据分析的同学会把RDA的结果当成终点图里画出样方和环境因子的箭头算出一个整体解释量再标上p值就完事了。但实际上这个整体解释量是一个大杂烩里面混着所有环境因子的共同贡献、独立贡献甚至还有与空间结构纠缠不清的部分。VPA存在的意义就是把这块大杂烩按因子组拆开让解释比例落到具体的生态过程上。1.1 一个典型场景p值很好看但回答不了多少假设你在研究某个湖泊浮游植物群落的空间分布。你采集了20个样点记录了物种多度同时测了水温、pH、总氮、总磷还记录了每个样点的经纬度坐标。你用RDA分析得到总解释量为38.6%p0.001图也漂亮——样点沿着温度和营养盐梯度分开。但问题来了这38.6%里有多少是环境因子单独贡献的有多少是空间位置单独贡献的又有多少是环境因子与空间位置无法区分的因为水温和总氮本身就存在空间分布趋势——湖心温度高、总氮高岸边温度低、总氮低——环境因子和空间坐标天然纠缠在一起。普通RDA无法把这两者拆开。如果你想在文章里说环境异质性解释了群落变异的X%空间过程如扩散限制解释了Y%那你手里只有整体38.6%这个数字根本说不出口。VPA就是为这种归因需求设计的。1.2 VPA的切蛋糕思路VPA的核心思路是把响应变量物种多度或群落组成的总变异看作一整块蛋糕然后按因子组切分成不同部分。最简单的情形是把解释变量分成两组比如环境因子和空间变量于是得到四块[a]环境因子单独解释的部分纯环境效应[c]空间变量单独解释的部分纯空间效应[b]环境与空间共同解释、无法区分归属的部分共享效应[d]完全未解释的残差这里有个容易混淆的点VPA不是替代RDA/CCA的独立方法而是在它们基础上的二次加工。每个组分的计算都要用到偏RDA或偏CCApartial RDA / partial CCA也就是在约束排序模型中加入协变量先控制住另一组因子的影响再看当前因子组的纯贡献。我习惯把VPA理解为正交化实验设计的思想搬到观测数据上。实验里我们通过随机分配来控制无关变量但观测数据没法随机分配环境因子只能用统计手段——把空间坐标当协变量控住或者把环境因子当协变量控住——来模拟控制变量的效果。当然它的效果不如真实验设计那么干净但已经是观测数据里能做的最好的归因手段之一。2. 变差分解的数学逻辑全模型、分组模型与减法理解VPA的数学逻辑不需要高深的统计背景核心就是先算总数再算纯贡献最后用减法求共享。我最初看varpart()输出时一脸懵直到自己手推了一遍才彻底明白。2.1 三个模型算出核心组分假设有两组解释变量X环境和W空间响应矩阵为Y。VPA的核心是三个模型全模型用X和W一起约束Y得到总解释量AdjR²_all。用RDA表达就是rda(Y ~ X W)。偏模型1在控制W的前提下看X对Y的纯解释量即rda(Y ~ X Condition(W))得到[a]。偏模型2在控制X的前提下看W对Y的纯解释量即rda(Y ~ W Condition(X))得到[c]。然后共享部分[b] AdjR²_all - [a] - [c]残差[d] 1 - AdjR²_all。为什么能这样减因为全模型的解释量包含X的贡献、W的贡献、以及X和W重叠部分的贡献。两个偏模型分别扣掉了另一组的影响所以各自只留下纯贡献重叠的部分自然就等于总体减去两个纯贡献的差值。这个减法逻辑很朴素但只有在使用校正R²时才成立下面细说。2.2 为什么必须用校正R²而不是原始R²如果你尝试用原始R²去做上面的加减法会发现一个致命问题原始R²只要往模型里加变量就必然增大哪怕是加一串随机数也会增大。这会导致全模型的R²虚高偏模型的R²也是虚高减法算出来的共享部分[b]可能变成负数或者超出合理范围。校正R²Adjusted R²会对模型中变量的个数做惩罚。每增加一个解释变量如果它不能带来足够的解释增量R²adj不仅不升还可能下降。这就保证了在不同模型之间做加减运算时用的是一把公平的尺子。vegan包的varpart()函数默认输出的就是校正后的R²表格里标注得很清楚。这一点特别重要如果你自己写代码手动计算一定要用RsquareAdj()函数拿校正R²不要直接取summary(rda_model)$r.squared。2.3 分组模型与生态学含义的对应很多教科书会把[a]、[b]、[c]、[d]直接对应到生态过程纯环境组分代表环境过滤的贡献纯空间组分代表空间过程如扩散限制、历史因素的贡献共享部分代表环境因子本身具有空间结构导致的不可区分贡献残差则代表随机过程、未测量的关键因子、以及物种间的相互作用等。不过这里要提醒一句这种对应关系是生态学解释上的约定俗成不是严格的因果推断。比如纯空间组分显著可能确实是扩散限制但也可能是因为你没测到的某个环境因子恰好呈空间分布。VPA能帮你量化可解释的比例但为什么是它解释的需要结合专业知识去论证。引用结果时最好写成空间变量独立解释了12%的群落变异这可能与扩散限制有关而不是斩钉截铁地说扩散限制贡献了12%。3. 完整实操流程从数据清洗到Venn图这一节是全文最能抄作业的部分。我用R的vegan包走一遍VPA的标准流程从数据准备到出图每一步都说明为什么这么做。3.1 数据准备的三张表做VPA需要三类数据建议在分析前整理成三个独立对象响应矩阵spe样方为行、物种为列的丰度矩阵。如果是微生物OTU表、浮游植物计数表、植物盖度表都可以但需要注意数据量纲和稀疏程度。环境因子框env样方为行、环境变量为列的数据框。pH、温度、养分含量等量纲各异通常需要标准化。空间变量spa/xy样方坐标或者由坐标生成的趋势面变量、距离矩阵特征向量。我见过不少人在数据准备上偷懒直接把原始经纬度扔进VPA。经纬度如果跨度很大比如几十公里两个坐标的数值量级和变异尺度差异可能影响分解结果。更规范的做法是用坐标生成空间变量常见有两种趋势面分析直接用坐标和高次项x, y, x², y², xy作为空间变量适合刻画大尺度的空间趋势。MEMMorans Eigenvector Maps对样点间的距离矩阵做特征分解得到多尺度的空间特征向量能同时捕捉大尺度和局部尺度的空间结构。R里adespatial包的dbmem()函数可以直接算。从原理上说MEM比趋势面更灵活但会生成很多特征向量需要筛选。实操中如果只想快速看一下空间效应直接标准化坐标也能用但写文章时建议用MEM审稿人观感更好。3.2 数据预处理的三步第一步物种数据做Hellinger转化。decostand(spe, method hellinger)。为什么要做因为原始多度数据里优势种和稀有种的数值差异极大直接做RDA会让少数优势种主导结果而Hellinger转化能降低这种偏差。如果你的数据是0/1有无数据或者丰度跨度不大这一步也可以省但多数情况下做了更稳妥。第二步环境因子标准化。scale(env)把每个变量变成均值0、方差1。这步是为了消除量纲影响——pH的数值范围是6到9总氮可能是0.1到5不标准化的话数值大的变量会在距离计算中占主导。第三步共线性筛选。用vif.cca()或普通的car::vif()检查环境因子的方差膨胀因子。VIF大于10的因子建议剔除严格一点大于5就剔除。共线性会在VPA里制造出虚高的共享部分后面第4节我会专门讲这个坑。3.3 核心代码varpart()与显著性检验下面是完整可运行的核心代码我加了一行行注释library(vegan) library(adespatial) # spe : 样方×物种丰度矩阵 # env : 样方×环境因子数据框 # xy : 样方坐标矩阵两列 # 1. 响应变量Hellinger转化 spe.hel - decostand(spe, method hellinger) # 2. 环境因子标准化 env.scale - as.data.frame(scale(env)) # 3. 由坐标生成空间特征向量MEM # 这一步会返回若干列空间变量代表不同尺度的空间结构 mem - dbmem(xy, thresh NULL) # 4. 如果MEM列数太多可以做前向选择精简 # library(adespatial) # mem.sel - forward.sel(spe.hel, mem, adjR2thresh 0.01) # mem - as.data.frame(mem[, mem.sel$order]) # 5. VPA核心分解环境组和空间组 vp - varpart(spe.hel, env.scale, mem, scale FALSE) vp # 直接打印分解表 # 6. 检验各组分显著性 # [a] 纯环境组分控制空间变量后做RDA test_env - rda(spe.hel, env.scale, mem) anova(test_env, permutations 999) # [c] 纯空间组分控制环境变量后做RDA test_spa - rda(spe.hel, mem, env.scale) anova(test_spa, permutations 999) # [ab] 环境总组分不加空间协变量 test_env_total - rda(spe.hel, env.scale) anova(test_env_total, permutations 999) # [cb] 空间总组分不加环境协变量 test_spa_total - rda(spe.hel, mem) anova(test_spa_total, permutations 999) # 7. 画Venn图 plot(vp, digits 2, bg c(#E69F00, #56B4E9, #009E73))这里有个细节varpart()里我写了scale FALSE是因为环境因子已经手动标准化了避免重复标准化。如果环境因子没做标准化建议scale TRUE让函数自动处理。3.4 结果表与Venn图的阅读方法运行vp后控制台会打印出类似这样的分解表组分含义解释量校正R²[a]环境因子纯解释量0.182[c]空间变量纯解释量0.073[b]环境与空间共享0.051[d]残差/未解释量0.694[ab]环境因子总解释量0.233[bc]空间变量总解释量0.124[abc]环境空间总解释量0.306Venn图里两个椭圆的重叠区域面积表示[b]。如果重叠区域很大说明环境因子与空间变量高度耦合比如环境因子本身有强烈的空间分布格局这时候讨论共享部分的生态学含义要格外谨慎它既可能是环境因子驱动了群落、又随空间变化也可能是空间过程影响了环境再影响群落因果关系并不明确。4. 实操中甩不掉的坑共享虚高、负解释量与因子组问题VPA看起来就是一个函数、一张图、几个数字实际操作里至少有三个坑是必然要踩的提前知道能省很多返工时间。4.1 共享部分爆炸共线性让环境与空间分不开我最开始做VPA时环境组里同时放了pH、总氮、总磷、电导率空间组放了MEM向量结果共享部分[b]占总解释量的一大半纯环境组分只有可怜的几个百分点。当时第一反应是程序算错了仔细排查后发现是共线性问题电导率和总离子浓度高度相关而它们又受地理位置的强烈影响相当于环境组里塞了太多重复信息这些信息又和空间组重叠共享部分自然膨胀。应对方式有两个层面。一是分析前做VIF筛选把VIF大于10的变量剔除或者用主成分分析把共线环境因子压缩成几个正交的主成分。二是调整因子组的划分逻辑把高度空间自相关的环境因子和代表纯空间过程的MEM向量视作两组时一定要清楚共享解释量在这个设定下就是这个结果不是bug。4.2 负解释量不是出错是没有纯贡献varpart()的表格里偶尔会出现负号比如[a] -0.021。很多人看到负数就以为模型跑错了。实际上校正R²本身允许负值它表示该因子组在扣除另一组协变量之后连解释量的基准线都没达到相当于纯贡献为零。换句话说这组因子单独拿出来在统计上没有任何解释能力。我处理负值的原则是报告中直接写该组分未检测到独立的解释贡献或解释量为0不要保留负号。Venn图里也会出现负数标注最好把负值统一改成0或者用不显示负值的绘图方式否则审稿人会追问。4.3 三组及以上因子组的划分学问VPA不只是两组vegan包支持三组甚至更多组。比如把解释变量分成环境因子空间变量土地利用变量三组分解表会变成八块Venn图变成三个圆。组数增加后共享部分的解释变得更复杂比如三组间的三重共享项[abc]三组共同解释的部分生态学含义往往非常难说清。我的建议是除非确有理论支撑否则尽量保持两组设计。三组VPA虽然看起来更精细但共享项太多写文章时会陷入这段共享究竟算谁的的泥潭。如果确实要做三组可以聚焦在某两组间的二元共享和该组的纯贡献上三组共享部分保守处理为不可分辨的联合效应。5. 一个完整实例复盘从原始数据到文章结论最后用一个我实际处理过、改写过结构的例子完整演示数据分析→结果解读→论文写作的链路让你看到VPA在真实研究里到底怎么落地。5.1 案例背景与数据研究对象是某流域的河流底栖动物群落。我们在30个样点定量采集了底栖动物同时记录了水温、溶解氧、pH、电导率、总氮、总磷六个环境因子以及每个样点的GPS坐标。科学问题是环境过滤和空间过程谁对底栖动物群落组成的影响更大物种数据是典型的多度数据存在大量零值优势种和稀有种差异大所以做Hellinger转化。环境因子经过VIF筛选后保留水温、pH、总氮、总磷四个变量电导率和总离子浓度VIF偏高剔除。空间变量用dbmem()生成了5个显著MEM向量。5.2 分析结果与解读varpart()的结果如下[a] 纯环境组分0.132p0.008[c] 纯空间组分0.064p0.012[b] 共享组分0.048[d] 残差0.756总解释量abc0.244结合生态学知识我的解读是环境因子独立解释了13.2%的群落变异说明环境过滤确实在起作用空间变量独立解释了6.4%暗示存在一定的扩散限制或空间自相关过程共享的4.8%说明部分环境因子本身有空间分布导致两者无法完全区分。残差75.6%看起来很大但在生态学观测数据里相当正常——随机过程、未测量的生物相互作用以及噪声都会归入这里。后来我在文章里特意强调了一个要点不要把残差大等同于模型失败。残差大只能说明我们手里的变量解释不了那么多变异不代表群落变化没有规律更不代表VPA没用。恰恰相反残差大小可以启发下一步研究——如果残差很大是否遗漏了关键环境因子是否应该做多营养级互作分析或功能基因层面研究5.3 文章里怎么写VPA结果给一段我实际写过的结果描述作为参考变差分解分析VPA结果表明环境因子与空间变量共同解释了底栖动物群落组成的24.4%的变异图X。其中环境因子独立解释了13.2%p0.008空间变量独立解释了6.4%p0.012两者重叠部分为4.8%。环境因子的独立贡献大于空间变量的独立贡献表明在流域尺度上环境过滤对底栖动物群落结构的塑造作用强于扩散限制这与流水中底栖动物普遍具有较强扩散能力的特点一致。这段写法遵循一个逻辑先给数字再给统计检验最后结合生态学机制解释。审稿人看到这段一般不会再有疑问。我在实际使用VPA时还有一个习惯就是把它和Mantel检验、PCNM空间尺度分析搭配着看。如果VPA算出来的纯空间组分显著我紧接着会去看不同尺度MEM向量的解释占比判断空间效应主要来自大尺度趋势还是局部斑块。这些信息拼在一起能让空间过程这个模糊的概念更具体。VPA不是一套孤立的方法它适合嵌入到环境-空间-群落这个三角关系的完整证据链里分析结论才站得住脚。