R统计绘图 |
您所在的位置：网站首页 › VPA是什么药 › R统计绘图 |
R统计绘图
|
变差分解分析(Variance Partitioning Analysis)可用于确定指定环境因子对微生物(原生生物/植物/动物等等)群落结构变化的解释比例。要计算指定环境因子与群落结构的相关性，就需要约束非指定环境因子的同时，对指定环境因子做排序分析。其实就是相当于做partial排序分析。公众号文章《R统计-PCA/PCoA/db-RDA/NMDS/CA/CCA/DCA等排序分析教程》写过如何使用vegan包进行偏分析。本文记录一下使用vegan包进行VPA分析的两种方法。 一、 数据准备 # 1. 设置工作路径，老生常谈，设置工作路径和查看路径内容一定是第一步 setwd("D:\\VPAanalysis") getwd() dir() # 2. 调用R包 install.packages("ecodist") library(vegan) browseVignettes("vegan") # 3. 读入数据 spe=read.csv("spe.csv",header=TRUE,row.names=1,sep=",", colClasses = c("character",rep("numeric",times=8)) )# 群落组成数据,物种组成数据是整数，为了使读入数据格式为数值型，设置colClasses，7为物种种类。 ## ?rep # 查看函数用法和函数中参数意义 spe group=read.csv("group.csv",header=TRUE,row.names=1,sep=",")#分类数据,包含两种类别，grazing和soil depth group # 这里不用group数据 env=read.csv("env.csv",header=TRUE,row.names=1,sep=",") #环境数据：土壤理化因子与气候因子。虚构数据仅用作教程 env
图1|物种组成数据(spe.csv)。
图2|环境因子数据(env.csv)。
图3|分类数据(group.csv)。 二、变差分解分析(Variation partitioning analysis,VPA) 2.1 方法一：RDA及VPA分析 使用rda()进行偏RDA分析，然后自行计算指定环境因子解释率以及不同环境因子方差解释率重叠部分。此部分也有两种计算方式。R统计-VPA分析(RDA/CCA)文中记录了先计算指定环境因子的partial effect(局部效应)，再使用总体环境因子效应-指定环境环境因子效应，得到共有方差解释率。此处记录计算指定环境因子的marginal effect(边际效应)，然后所有边际效应相加-环境因子总体方差解释率，得到环境因子对微生物群落结构的解释率的共有部分。 # 2.1.1 RDA-全部环境因子。 RDA=rda(decostand(spe, "hel"),env) RDA sumr=summary(RDA) # 2.1.2 RDA-WC、pH和TC ## RDA-WC、pH和TC的局部效应：单独只由WC、pH和TC解释的方差 RDA2=rda(decostand(spe, "hel"),env[c(1,2,5)],env[c(3,4,6)]) RDA2 sumr2=summary(RDA2) anova(RDA2) # 显著性检验 #permutest(RDA2,permu=999) # anova()作用一样 RsquareAdj(RDA2) # 指定环境因子相关性结果校正，环境因子分类中不止一个环境因子，需要对R^2结果进行校正。 ## RDA-WC、pH和TC的边际效应：只由WC、pH和TC解释的方差+WC、pH和TC与其它环境因子共同解释的方差 RDA4=rda(decostand(spe, "hel"),env[c(1,2,5)]) RDA4 sumr4=summary(RDA4) # 2.1.3 RDA-氮形态环境因子 ## 氮形态环境因子的局部效应 RDA3=rda(decostand(spe, "hel"),env[c(3,4,6)],env[c(1,2,5)]) RDA3 sumr3=summary(RDA3) ## 氮形态环境因子的边际效应 RDA5=rda(decostand(spe, "hel"),env[c(3,4,6)]) RDA5 sumr5=summary(RDA5) # 2.1.4 计算WC、pH和TC与氮形态环境因子对微生物群落变化的解释度的共有部分 ## 环境因子间的共线性，导致对群落结构方差贡献率的重叠。 ## 局部效应计算方式 con=sumr$constr.chi/sumr$tot.chi-(sumr2$constr.chi/sumr2$tot.chi+sumr3$constr.chi/sumr3$tot.chi) con ## 边际效应计算方式 con2=(sumr4$constr.chi/sumr4$tot.chi+sumr5$constr.chi/sumr5$tot.chi)-sumr$constr.chi/sumr$tot.chi con2 # 两种计算方法的结果一致 # 2.1.5 VPA分析结果 ## 局部效应计算方式 data = data.frame(RDA2=sumr2$constr.chi/sumr2$tot.chi,con=con,RDA3=sumr3$constr.chi/sumr3$tot.chi,un=sumr$unconst.chi/sumr$tot.chi) data ## 边际效应计算方式 data2 = data.frame(`WC+pH+TC`=sumr4$constr.chi/sumr4$tot.chi-con,con=con,`N form`=sumr5$constr.chi/sumr5$tot.chi-con,un=sumr$unconst.chi/sumr$tot.chi) data2注：如果群落结构数据是高通量测序数据，存在很多物种为0的情况，可选择对数据进行hellinger转换，避免“弓形效应”。“弓形效应”就是CA/RA的第二排序轴在许多情况下是第一轴的二次变形。此部分的输出结果在R统计-VPA分析(RDA/CCA)文中，已经解释过了。这里不再赘述。 2.2 方法二：VPA分析 直接使用varpart()函数进行VPA分析。 # 2.2.1 两个环境因子分类VPA分析 ## transfo参数定义微生物群落结构数据的转换方法："hellinger", "chi.square", "total", "norm"可选，当已经对微生物群落结构数据进行了矩阵计算，则不设置此参数。 ## chisquare = FALSE表示进行partial RDA分析，chisquare = TRUE则进行partial CCA分析。 rda.vpa |
【本文地址】
今日新闻 |
推荐新闻 |