基于R的群落学多元统计分析

1、1基于 的群落多元统计分析 用vegan包进行排序分析赖江山（janson) 中国科学院植物研究所2010. 11.52vegan软件包简介软件包简介3vegan是是Vegetation analysis的缩写的缩写,是群落分析的是群落分析的package作者：作者：Jari Oksanen /web/packages/vegan/index.htmlhttp:/cc.oulu.fi/jarioksa/softhelp/vegan.htmllibrary(vegan)什么是排序什么是排序(ordination)?排序的过程是将样方或植物种排列在一

2、定的空间，排序的过程是将样方或植物种排列在一定的空间，使得排序轴能够反映一定的生态梯度，从而，能使得排序轴能够反映一定的生态梯度，从而，能够解释植被或植物种的分布与环境因子间的关系，够解释植被或植物种的分布与环境因子间的关系，也就是说排序是为了揭示植被也就是说排序是为了揭示植被-环境间的生态关系。环境间的生态关系。因此，排序也叫梯度分析（因此，排序也叫梯度分析（gradient analysis）。）。间接梯度分析间接梯度分析 (Indirect gradient analysis)直接梯度分析直接梯度分析 (direct gradient analysis)2个种的排序图个种的排序图3个种的

3、排序图个种的排序图4个种的排序图？个种的排序图？40个种排序图？个种排序图？排序的目标：排序的目标：1.降低维数，减少坐标轴的数目降低维数，减少坐标轴的数目 ；2. 由降低维数引起的信息损失尽量少，即发生最由降低维数引起的信息损失尽量少，即发生最小的畸变，也就是让新的坐标系第小的畸变，也就是让新的坐标系第1-3轴排序轴包轴排序轴包含大量的生态信息含大量的生态信息 。排序的目的：排序的目的：表示植被与环境之间的关系：表示植被与环境之间的关系：所有排序方法都反映植物种和环境之间的关系以所有排序方法都反映植物种和环境之间的关系以及在某一环境梯度上的种间关系。及在某一环境梯度上的种间关系。1.线形模型

4、（线形模型（linear model），短的梯度，主成分），短的梯度，主成分分析（分析（Principle component analysis)，需要对数据，需要对数据进行非线性转换，如取对数；进行非线性转换，如取对数；2.非线性模型（非线性模型（non-linear model）如高斯模型，）如高斯模型，长的梯度，对应分析长的梯度，对应分析 （Correspondence analysis) 群落数据输入群落数据输入gtsdata=read.table(gtsdata.txt,header=T)gtsdatadim(gtsdata)环境因子数据输入环境因子数据输入gtsenv=read.t

5、able(gtsenv.txt,header=T); gtsenvdim(gtsenv)数据的标准化数据的标准化1. decostand(x, method, MARGIN, )total: 除以行和或列和除以行和或列和 (default MARGIN = 1是是row);max:除以行或列的最大值除以行或列的最大值(default MARGIN = 2 是列是列);freq:除以行或列的最大值除以行或列的最大值,并乘以非零值的个数，非零值的平并乘以非零值的个数，非零值的平均值为均值为1 (default MARGIN=2);normalize:使行或列的平方和等于使行或列的平方和等于1 (d

6、efault MARGIN = 1);range: 标准化使行或列的值在标准化使行或列的值在0 . 1 (default MARGIN = 2).standardize:标准化使行或列的和为标准化使行或列的和为1且方差为且方差为1(default MARGIN = 2)；pa: 将数据转换为将数据转换为0、1数据；数据； chi.square: 除以行和及列和的平方根；除以行和及列和的平方根；hellinger: 采用采用total标准化以后再取平方根；标准化以后再取平方根；log：对数化，默认自然对数，对数化，默认自然对数，logbase参数是自选的参数是自选的base2. wisconsi

7、n(x) :除以列最大值，再除以行和。除以列最大值，再除以行和。排序类别排序类别(in CANOCO)间接梯度分析（间接梯度分析（Indirect Gradient Analysis） ：PCA (Principal components analysis)CA (Correspondence analysis)DCA (Detrended Correspondence Analysis)直接梯度分析（直接梯度分析（Direct Gradient Analysis） ：RDA (Redundance analysis)CCA (Canonical correspondence analysis

8、)DCCA (Detrended CCA )PCA RDA CA CCA DCA DCCA13决定排序的模型：单峰还是线性？决定排序的模型：单峰还是线性？decorana(gtsdata)Call:decorana(veg = gtsdata) Detrended correspondence analysis with 26 segments.Rescaling of axes with 4 iterations. DCA1 DCA2 DCA3 DCA4Eigenvalues 0.3939 0.2239 0.09555 0.06226Decorana values 0.5025 0.1756

9、 0.06712 0.03877Axis lengths 3.2595 2.5130 1.21445 1.00854如果这四个轴中梯度最长（最大值）超过4，选择单峰模型排序（CA、CCA、DCA）更合适。如果是小于3，选择线性模型（PCA、RDA）比较合理。如果介于3-4之间，单峰模型和线性模型结果差不多。间接梯度分析间接梯度分析（Indirect Gradient Analysis）PCA (Principal components analysis)CA (Correspondence analysis)DCA (Detrended Correspondence Analysis)主成分分

10、析主成分分析(Principle component analysis, PCA)主成分分析的主要原主成分分析的主要原理是：理是：使坐标旋转一定的角使坐标旋转一定的角度后，使第一轴表示数度后，使第一轴表示数据最大的方差，使第二据最大的方差，使第二轴表示数据第二的方差。轴表示数据第二的方差。而且轴与轴之间是正交而且轴与轴之间是正交的（的（orthogonal)。-4-2024-4-2024 PCA和和RDA都采用函数都采用函数rda实现：实现：在在vegan包中，包中，rda(formula, data, scale=FALSE, .) rda(X, Y, Z, scale=FALSE, .)

11、scores(x, choices, display=c(sites,species), .) 在在stat包中：包中：princomp(x, .) 主成分分析主成分分析princomp(formula, data = NULL, subset, na.action, .) gts.rda=rda(gtsdata)gts.rdaCall: rda(X = gtsdata) Inertia RankTotal 352.1 Unconstrained 352.1 22Inertia is variance Eigenvalues for unconstrained axes: PC1 PC2 PC

12、3 PC4 PC5 PC6 PC7 PC8 111.779 73.580 54.607 32.959 26.481 18.063 12.763 7.637 scores(gts.rda,choices=c(1:4) ,display=c(si,sp)summary(gts.rda) #类似类似Canoco的的log文件和文件和.sol文件的信息文件的信息plot(gts.rda,choices=c(1,2),display=c(sp,si)biplot(gts.rda,choices=c(1,2),display=c(sp,si)plot(rda(gtsdata,scale=T)plot(rd

13、a(gtsdata)!如果不对数据做标准化的话，丰富种的值就非常大，排序如果不对数据做标准化的话，丰富种的值就非常大，排序时就只能看清丰富种的位置，其它种就拥挤在一起。时就只能看清丰富种的位置，其它种就拥挤在一起。 如用如用x x1 1, ,x x2 2, ,x x3 3, ,x x4 4, ,x x5 5, ,x x66分别表示原先的变分别表示原先的变量，而用量，而用y y1 1, ,y y2 2, ,y y3 3, ,y y4 4, ,y y5 5, ,y y66表示新的主成表示新的主成分，那么，第一和第二主成分为分，那么，第一和第二主成分为11234562123456-0.806-0.6

14、74-0.6750.8930.8250.8360.3530.5310.5130.3060.4350.425yxxxxxxyxxxxxx 这些系数称为主成分载荷（这些系数称为主成分载荷（loading），它表示），它表示主成分和相应的原先变量的相关系数。比如主成分和相应的原先变量的相关系数。比如y1表表示式中示式中x1的系数为的系数为-0.806，这就是说第一主成分，这就是说第一主成分和的和的x x1 1变量的相关系数为变量的相关系数为-0.806。相关系数。相关系数(绝对绝对值）越大，主成分对该变量的代表性也越大。值）越大，主成分对该变量的代表性也越大。负荷负荷(loading)gts.pca

15、=princomp(gtsdata)gts.pca$loadingsgtsenv.pca=princomp(gtsdenv)gtsenv.pca$loadingsbiplot(gtsenv.pca) 第一主成分代表海拔高度，第二第一主成分代表海拔高度，第二主成分代表坡向主成分代表坡向对应分析对应分析(Correspondence analysis, CA)1.PCA在迭代运算过程是采用线性模型在迭代运算过程是采用线性模型2.CA在迭代运算过程采用单峰模型（加权平均法）在迭代运算过程采用单峰模型（加权平均法）CA在在vegan中也是用中也是用cca函数来实现：函数来实现：gts.ca=cca(g

16、tsdata)gts.casummary(gts.ca)Call:cca(X = gtsdata) Inertia RankTotal 1.424 Unconstrained 1.424 21Inertia is mean squared contingency coefficient Eigenvalues for unconstrained axes: CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8 0.50253 0.26564 0.14023 0.10502 0.09127 0.05540 0.05063 0.04204 26plot(gts.ca)CANOCO里面 s

17、caling of ordination scores27plot(gts.ca, scaling=1)用物种数据对样方坐标进行加权平均用物种数据对样方坐标进行加权平均,使样方坐使样方坐标在物种数据的中心，因此对样方感兴趣的话，标在物种数据的中心，因此对样方感兴趣的话，采用这种做图方法。采用这种做图方法。plot(gts.ca, scaling=2)plot(gts.ca, scaling=3)用样方数据对物种坐标进行加权平均用样方数据对物种坐标进行加权平均,使物种数使物种数据在样方数据的中心，因此对物种感兴趣的话，据在样方数据的中心，因此对物种感兴趣的话，采用这种做图方法。采用这种做图方法。

18、如果一个物种靠近某个样方，表明该物种可能对如果一个物种靠近某个样方，表明该物种可能对该样方的位置起很大的作用。特别是对于二元数该样方的位置起很大的作用。特别是对于二元数据的排序，这个样方可能就代表该物种。据的排序，这个样方可能就代表该物种。如图中，如图中，20号样方号样方与短柄枹与短柄枹(QUESER)靠得比较近，表明：靠得比较近，表明：短柄枹表征了短柄枹表征了20号号样方的特征，样方的特征，19号号样方与样方与20号样方距号样方距离近，生态关系也离近，生态关系也较近。较近。 只在少数样方出现的物只在少数样方出现的物种通常在排序空间的边缘，种通常在排序空间的边缘，表明它们只偶然发生，或表明它们

19、只偶然发生，或它们只在稀有生境（如米它们只在稀有生境（如米槠槠CASCAR）。）。 在排序空间中心的物种，在排序空间中心的物种，可能在取样区域是该物种可能在取样区域是该物种最优分布区，如甜槠，或最优分布区，如甜槠，或有两个或多个最优分布区，有两个或多个最优分布区，或与前两个轴不相关。或与前两个轴不相关。除趋势对应分析除趋势对应分析(Detrended correspondence analysis, DCA)：CA采用单峰曲线表示物种和环境关系采用单峰曲线表示物种和环境关系CA产生的弓形效应产生的弓形效应CA的第二排序轴在的第二排序轴在许多情况下是第一许多情况下是第一轴的二次变形，即轴的二次变

20、形，即所谓的所谓的“弓形效应弓形效应”（Arch effect）或）或者者“马蹄形效应马蹄形效应”（horse-shoe effect）。）。 （详见张金屯群落生态学168页） 和分别代表除趋势前和除趋势后的样方排序（引自Hill 和Gauch 1980）DCA在在R中的实现采用函数中的实现采用函数decorana。decorana(veg, iweigh=0, iresc=4, ira=0, mk=26, short=0, before=NULL, after=NULL) veg:群落数据群落数据;iweig:稀有物种的权重；稀有物种的权重；(稀有物种影响比较大）稀有物种影响比较大）ires

21、c：纠正弓形效应的次数；：纠正弓形效应的次数；ira:分析的类型分析的类型(DCA:0 , CA:1);mk：校正弓形效应轴的分段数：校正弓形效应轴的分段数;short: 需要校正的最短梯度。需要校正的最短梯度。plot(decorana(gtsdata)plot(cca(gtsdata)直接梯度分析（直接梯度分析（Direct Gradient Analysis） ：RDA (Redundance analysis)CCA (Canonical correspondence analysis)pRDA (partial RDA)pCCA (partial CCA)冗余分析冗余分析(redun

22、dancy analysis, RDA)及典范对应分析（及典范对应分析（Canonical correspondence analysis, CCA)1.通常采用通常采用PCA处理环境数据，采用处理环境数据，采用CA处理群落处理群落数据，但这些方法都只能处理一个数据表；数据，但这些方法都只能处理一个数据表；2.RDA和和CCA是多元分析（是多元分析（PCA，CA）和线性）和线性回归的结合，研究植被和环境之间的关系。回归的结合，研究植被和环境之间的关系。38 当我们在解释变量（环境因子数据）与响应变量（物种数据）之间建立预测模型的时候，经常会遇到这样的情况，往往我们仅仅考察解释变量中某几个环境因

23、子的对物种数据的影响，但剩下的环境因子也会对物种产生影响，这些剩余环境因子我们经常称为协变量(Covariables) 。在CANOCO中，协变量的影响可以用偏分析（partial analyze）剔除出来。 实际上，任何一个环境因子变量均可以成为协变量。例如，我们要研究管理模式对蝴蝶群落中组成的影响，我们可以在不同的海拔地点取样，海拔也许对群落物种组成影响很大，但此时我们感兴趣的是管理模式的影响，而非海拔梯度的影响。这个时候，如果能剔除出海拔的影响，我们能管理模型与蝴蝶种群之间更清晰的关系。 摘自 第一章6pPCA与环境因子结合是与环境因子结合是RDA，CA与环境因子结与环境因子结合是合是C

24、CA。RDA在在vegan中的实现：中的实现：rda(formula, data, .) rda(X, Y, Z, .) CCA在在vegan中的实现：中的实现：cca(formula, data, .) cca(X, Y, Z, .)gts.rda=rda(gtsdata,gtsenv);gts.rda=rda(gtsdataelev+convex+slope+aspect+N+K+P+pH,data=gtsenv);summary(gts.rda)plot(gts.rda)gts.prda=rda(gtsdata,gtsenv,1:4, gtsenv,5:8)summary(gts.prd

25、a)plot(gts.prda)plot(gts.rda)环境因子一般用箭头表环境因子一般用箭头表示，箭头所处的象限表示，箭头所处的象限表示环境因子与排序轴间示环境因子与排序轴间的正负相关性，箭头连的正负相关性，箭头连线的长度代表着某个环线的长度代表着某个环境因子与群落分布和种境因子与群落分布和种类分布间相关程度的大类分布间相关程度的大小，连线越长，说有相小，连线越长，说有相关性越大。反之越小。关性越大。反之越小。箭头连线和排序轴的夹箭头连线和排序轴的夹角代表着某个环境因子角代表着某个环境因子与排序轴的相关性大小，与排序轴的相关性大小，夹角越小，相关性越高；夹角越小，相关性越高；反之越低。反之

26、越低。 gts.cca=cca(gtsdata,gtsenv);gts.cca=cca(gtsdataelev+convex+slope+aspect+N+K+P+pH,data=gtsenv) summary(gts.cca)plot(gts.cca)gts.pcca=cca(gtsdata,gtsenv,1:4, gtsenv,5:8)summary(gts.pcca)plot(gts.pcca)plot(gts.cca)RDARDA和和CCACCA结果的检验：结果的检验：goodness(object, display = c(“species”, “sites”), choices,

27、model = c(“CCA”, “CA”), statistic = c(“explained”, “distance”), summarize = FALSE, .) ：物种或样方与轴累计解释量；：物种或样方与轴累计解释量；【Cumulative fit per species as fraction of variance of species】in CANOCOinertcomp(object, display = c(“species”, “sites”), statistic = c(“explained”, “distance”), proportional = FALSE) ：

28、物种或样方的方差分解分析；：物种或样方的方差分解分析；spenvcor(object) :物种和环境的相关分析物种和环境的相关分析;intersetcor(object) :环境因子和各轴的相关分析；环境因子和各轴的相关分析；最好使用：最好使用：envfit;vif.cca(object) ：方差膨胀因子分析；方差膨胀因子分析；RDARDA和和CCACCA结果的检验：结果的检验：envfit(X, P, permutations = 0, strata, choices=c(1,2), .) envfit(formula, data, .) gts.cca=cca(gtsdata,gtsenv

29、)ef=envfit(gts.cca,gtsenv,permu=1000)CCA1 CCA2 r2 Pr(r) elev -0.9999873 -0.0050479 0.5005 0.001 *convex -0.7746504 0.6323898 0.0844 0.28 slope 0.1526685 0.9882775 0.5296 0.001 *aspect -0.1364818 -0.9906426 0.0169 0.82 N 0.6878942 0.7258110 0.6342 0.001 *P 0.8101502 0.5862224 0.5856 0.001 *K 0.600498

30、4 0.7996260 0.3364 0.001 *pH -0.1929061 -0.9812172 0.6313 goodness(gts.cca,display=“sp”)inertcomp(gts.cca, proportional = T) spenvcor(gts.cca) intersetcor(gts.cca) vif.cca(gts.cca) RDARDA和和CCACCA模型的选择：模型的选择： step(mod1,scope)mod1=cca(gtsdata1,data=gtsenv);mod2=cca(gtsdataelev+convex+slope+aspect+P+N+

31、K+pH,data=gtsenv);mod=step(mod1,scope=(list(lower=formula(mod1),upper=formula(mod2);Start: AIC=180.2gtsdata 1 Df AIC+ P 1 174.66+ elev 1 174.67+ N 1 174.89+ pH 1 176.95+ slope 1 177.82+ K 1 178.14 180.20+ convex 1 180.85+ aspect 1 181.53RDARDA和和CCACCA结果的检验：结果的检验： gts.cca=cca(gtsdataelev+convex+slope

32、+aspect+P+N+K+pH,data=gtsenv)anova(gts.cca)anova(gts.cca,by=“axis”)anova(gts.cca,by=“terms”)Model: cca(formula = gtsdata elev + convex + slope + aspect + P + N + K + pH, data = gtsenv) Df Chisq F N.Perm Pr(F) elev 1 0.2446 7.0228 100 0.01 *convex 1 0.0624 1.7902 100 0.07 . slope 1 0.1356 3.8928 100

33、0.01 *aspect 1 0.0208 0.5981 100 0.55 P 1 0.0911 2.6159 100 modCall:cca(formula = gtsdata P + pH + elev + convex, data = gtsenv) Inertia RankTotal 1.4243 Constrained 0.5640 4Unconstrained 0.8603 21Inertia is mean squared contingency coefficient Eigenvalues for constrained axes: CCA1 CCA2 CCA3 CCA4 0

34、.31595 0.18030 0.04863 0.01913 Eigenvalues for unconstrained axes: CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8 0.27067 0.13398 0.08373 0.06294 0.05026 0.04852 0.04199 0.02758 作图函数：作图函数：plot(x, choices = c(1, 2), display = c(sp, wa, cn), scaling = 2, type, xlim, ylim, .) x: cca或或rda分析结果；分析结果；choices:选择的轴；选择的轴；di

35、splay:显示的类型，显示的类型，sp,物种，物种，wa,样方，样方，lc, 线性结合的样方坐标线性结合的样方坐标;作图函数：作图函数：text(x, display = sites, labels, choices = c(1, 2), scaling = 2, arrow.mul, head.arrow = 0.05, select, .)x: cca或或rda分析结果；分析结果；display:显示的类型，显示的类型，sp,wa, lc,bp,环境因环境因子子;cn,中心化后的环境因子；中心化后的环境因子；label:用来替代箭头的字符；用来替代箭头的字符；Arrow.mul:环境因子

36、放大倍数环境因子放大倍数;作图函数作图函数:points(x, display = sites, choices = c(1, 2), scaling = 2, arrow.mul, head.arrow = 0.05, select, .)plot(gts.cca,type=n);text(gts.cca,display=“cn”,arrow.mul=1.5,col=4); points(gts.cca,display=lc,pch=16,col=2,select=1:20)points(gts.cca,display=“lc”,pch=2,col=4,select=21:40)text(g

37、ts.cca,display=“lc”,col=4”);绘制三维排序图：绘制三维排序图：ordiplot3d(object, display = sites, choices = 1:3, ax.col = 2, arr.len = 0.1, arr.col = 4, envfit, xlab, ylab, zlab, .) ordirgl(object, display = sites, choices = 1:3, type = p, ax.col = red, arr.col = yellow, text, envfit, .) orglpoints(object, display =

38、sites, choices = 1:3, .) orgltext(object, text, display = sites, choices = 1:3, justify = center, adj = 0.5, .) ordiplot3d(gts.cca,type=“h”) ordirgl(gts.cca,type=p) 绘制环境因子梯度图：绘制环境因子梯度图：ordisurf(x, y, choices=c(1, 2), knots=10, family=gaussian, col=red, thinplate = TRUE, add = FALSE, display = sites, w = weights(x), main, nlevels = 10, l