最近处理数据，需要知道环境因子与生物群落的关系，在看文献的过程中，发现了一个VPA的分析，自己想要复现文章中的每个环境因子单独的解释效应（图4D）。研究了一下，特地写成推文一方面帮助需要的人，一方面记录一下自己的学习过程。

首先先浅浅的介绍一下VPA分析，如下图所示：

图1   以3组解释变量X1、X2和X3为例, 使用韦恩图概括了响应变量Y的总变差组成。A, 基于变量的顺序R2的方法, 根据分配的顺序, [a]、[b]和[c]分别为X1、X2和X3的条件效应。B, 基于变差分解和层次分割法, [a]、[b]和[c]分别为X1、X2和X3的边际效应, [d]、[e]和[f]分别为X1、X2和X3两两之间的共同效应, [g]为三者之间的共同效应, X1、X2和X3的单独效应分别为[a] + [d]/2 + [f]/2 + [g]/3、[b] + [d]/2 + [e]/2 + [g]/3以及[c] + [e]/2 + [f]/2 + [g]/3。残差[h]代表了未被X1、X2和X3解释的Y的变差部分。（Liu et al.,2023）

变差分解广泛应用于生态学的分析，但是传统的vegan包中的varpart函数只能处理不超过4组的解释变量，并且解释变量多的话，看起来会很复杂，共同的解释效应也没有办法量化，也不能很方便的看到每个环境因子的单独贡献。刚好赖江山老师新近发表的rdacca.hp函数包（LAI et al.,2022）可以解决传统的varpart函数进行变差分解的缺点，但是我在尝试重现论文中目标图的结果时也遇到了问题，那就是我按照varpart函数进行VPA分析与使用rdacca.hp函数进行VPA分析计算的环境因子边际效应不一样，我现在没有想清楚这个问题，不知道是不是采取的算法不一致导致的。

以下是代码部分：

##RDA 的变差分解

#使用数据为发表论文的附件数据，大家如果需要可以直接去下载这是链接：https://www.sciencedirect.com/science/article/pii/S0043135423004839?via%3Dihub

#设置工作目录

setwd("按照你自己文件存放的路径填进来")

#安装rdacca.hp,vegan包

install.packages("rdacca.hp")

install.packages("vegan")

#加载包

library(rdacca.hp)

library(vegan)

#读取物种/ASV数据，输入文件格式为行为样本列为物种/ASV，如果不是则需要转置

spe <- read.csv('asv.csv', row.names = 1, header = T)

spe <- data.frame(t(phylum)) #转置数据

#推荐在对纯物种多度数据执行 PCA、RDA 等线性排序方法前，执行 Hellinger 转化

#示例的物种数据应该是已经执行过转化，因此不需要转化

#spe_hel <- decostand(spe, method = 'hellinger')

#读取环境变量

env <- read.csv('env.csv', row.names = 1, header = T)

#运行变差分解，首先使用vegan包中的varpart函数进行VPA分析，比较异同。

#以两组（以TN为一组，'pH', 'EC','WT','DO','TUR','WS','NH3-N','COD','TP','NO3-N'指标为另一组）环境变量为例

rda_vp <- varpart(spe, env['TN'], env[c('pH', 'EC','WT','DO','TUR','WS','NH3.N','COD','TP','NO3.N')])

rda_vp

plot(rda_vp, digits = 2, Xnames = c('TN', 'pH + EC + WT + DO + TUR + WS + NH3.N + COD + TP + NO3.N'), bg = c('blue', 'red'))

这里面TN的边际效应，是与文章中不符合的，并且与rdacca.hp计算出来的也不符合，而rdacca.hp计算出来的边际效应与文章中的结果吻合。

#视自己的环境数据而定，选择是否要对数转化

env.log <- log1p(env) 

#使用rdacca.hp进行VPA分析

cap.hp = rdacca.hp(spe, env, method = 'RDA', type = 'R2', scale = FALSE, var.part = TRUE) 

cap.hp$Total_explained_variation #提取所有环境因子对群落结构变化的总解释率，代表所有环境因子能对群落结构解释的方差部分

0.48

cap.hp$Hier.part #每个环境因子对群落结构的解释率，所有环境因子的解释率之和等于“cap.hp$Total_explained_variation”即所有环境因子对群落结构变化的总解释率，即层次分割

cap.hp$Var.part #如果Var.part=TRUE，则包含所有解释变量共同的解释率和百分比的矩阵(N个预测器或矩阵为2^N-1)。即方差分解（VPA）。

这是所有的结果后续可以把结果提取出来，用来绘制韦恩图或者柱形图

参考文献：

LIU Yao, YU Xin, YU Yang, HU Wen-Hao, LAI Jiang-Shan. Application of “rdacca.hp” R package in ecological data analysis: case and progress[J]. Chin J Plant Ecol, 2023, 47(1): 134-144.

Lai JS, Zou Y, Zhang JL, Peres-Neto PR.Generalizing hierarchical and variation partitioning in multiple regression and canonical analyses using the rdacca.hp R package Methods in Ecology and Evolution, 2022，13, 782-788.

转载本文请联系原作者获取授权，同时请注明本文来自王超凡科学网博客。
链接地址：https://blog.sciencenet.cn/blog-3436070-1409048.html

上一篇：高通量测序原理及特点——以Illumina测序平台为例-王超凡的博文