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

《基于R的群落学多元统计分析.ppt》由会员分享，可在线阅读，更多相关《基于R的群落学多元统计分析.ppt（65页珍藏版）》请在三一文库上搜索。

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

2、轴能够反映一定的生态梯度，从而，能够解释植被或植物种的分布与环境因子间的关系，也就是说排序是为了揭示植被-环境间的生态关系。因此，排序也叫梯度分析（gradient analysis）。 间接梯度分析 (Indirect gradient analysis) 直接梯度分析 (direct gradient analysis),2个种的排序图,3个种的排序图,4个种的排序图？,40个种排序图？,基础知识,、三相交流电： 由三个频率相同、电势振幅相等、相位差互差 120 角的交流电路组成的电力系统，叫三相交流电。 、一次设备： 直接与生产电能和输配电有关的设备称为一次设备。包括各种高压断路器、隔离

3、开关、母线、电力电缆、电压互感器、电流互感器、电抗器、避雷器、消弧线圈、并联电容器及高压熔断器等。,基础知识,、二次设备： 对一次设备进行监视、测量、操纵控制和保护作用的辅助设备。如各种继电器、信号装置、测量仪表、录波记录装置以及遥测、遥信装置和各种控制电缆、小母线等。 4、负荷开关： 负荷开关的构造与隔离开关相似，只是加装了简单的灭弧装置。它也是有一个明显的断开点，有一定的断流能力，可以带负荷操作，但不能直接断开短路电流，如果需要，要依靠与它串接的高压熔断器来实现。,排序的目标： 1.降低维数，减少坐标轴的数目 ； 2. 由降低维数引起的信息损失尽量少，即发生最小的畸变，也就是让新的坐标系第

4、1-3轴排序轴包含大量的生态信息 。,排序的目的： 表示植被与环境之间的关系： 所有排序方法都反映植物种和环境之间的关系以及在某一环境梯度上的种间关系。 线形模型（linear model），短的梯度，主成分分析（Principle component analysis)，需要对数据进行非线性转换，如取对数； 非线性模型（non-linear model）如高斯模型，长的梯度，对应分析 （Correspondence analysis),群落数据输入,gtsdata=read.table(“gtsdata.txt“,header=T) gtsdata dim(gtsdata),环境因子数据输入

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

6、. 1 (default MARGIN = 2). standardize:标准化使行或列的和为1且方差为1(default MARGIN = 2)； pa: 将数据转换为0、1数据； chi.square: 除以行和及列和的平方根； hellinger: 采用total标准化以后再取平方根； log：对数化，默认自然对数，logbase参数是自选的base 2. wisconsin(x) :除以列最大值，再除以行和。,排序类别(in CANOCO) 间接梯度分析（Indirect Gradient Analysis） ： PCA (Principal components analysis)

7、 CA (Correspondence analysis) DCA (Detrended Correspondence Analysis) 直接梯度分析（Direct Gradient Analysis） ： RDA (Redundance analysis) CCA (Canonical correspondence analysis) DCCA (Detrended CCA ),PCA RDA,CA CCA DCA DCCA,15,决定排序的模型：单峰还是线性？ decorana(gtsdata) Call: decorana(veg = gtsdata) Detrended corres

8、pondence analysis with 26 segments. Rescaling of axes with 4 iterations. DCA1 DCA2 DCA3 DCA4 Eigenvalues 0.3939 0.2239 0.09555 0.06226 Decorana values 0.5025 0.1756 0.06712 0.03877 Axis lengths 3.2595 2.5130 1.21445 1.00854,如果这四个轴中梯度最长（最大值）超过4，选择单峰模型排序（CA、CCA、DCA）更合适。如果是小于3，选择线性模型（PCA、RDA）比较合理。如果介于3

9、-4之间，单峰模型和线性模型结果差不多。,间接梯度分析 （Indirect Gradient Analysis） PCA (Principal components analysis) CA (Correspondence analysis) DCA (Detrended Correspondence Analysis),主成分分析(Principle component analysis, PCA),主成分分析的主要原理是： 使坐标旋转一定的角度后，使第一轴表示数据最大的方差，使第二轴表示数据第二的方差。而且轴与轴之间是正交的（orthogonal)。,PCA和RDA都采用函数rda实现：

10、在vegan包中， rda(formula, data, scale=FALSE, .) rda(X, Y, Z, scale=FALSE, .) scores(x, choices, display=c(“sites“,“species“), .) 在stat包中： princomp(x, .) 主成分分析 princomp(formula, data = NULL, subset, na.action, .),gts.rda=rda(gtsdata) gts.rda Call: rda(X = gtsdata) Inertia Rank Total 352.1 Unconstrained

11、352.1 22 Inertia is variance Eigenvalues for unconstrained axes: PC1 PC2 PC3 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“) bip

12、lot(gts.rda,choices=c(1,2),display=c(“sp“,“si“),plot(rda(gtsdata,scale=T),plot(rda(gtsdata),!如果不对数据做标准化的话，丰富种的值就非常大，排序时就只能看清丰富种的位置，其它种就拥挤在一起。,如用x1,x2,x3,x4,x5,x6分别表示原先的变量，而用y1,y2,y3,y4,y5,y6表示新的主成分，那么，第一和第二主成分为,这些系数称为主成分载荷（loading），它表示主成分和相应的原先变量的相关系数。比如y1表示式中x1的系数为-0.806，这就是说第一主成分和的x1变量的相关系数为-0.806

13、。相关系数(绝对值）越大，主成分对该变量的代表性也越大。,负荷(loading),gts.pca=princomp(gtsdata) gts.pca$loadings,gtsenv.pca=princomp(gtsdenv) gtsenv.pca$loadings,biplot(gtsenv.pca),第一主成分代表海拔高度，第二主成分代表坡向,对应分析(Correspondence analysis, CA) PCA在迭代运算过程是采用线性模型 CA在迭代运算过程采用单峰模型（加权平均法）,CA在vegan中也是用cca函数来实现： gts.ca=cca(gtsdata) gts.ca su

14、mmary(gts.ca) Call: cca(X = gtsdata) Inertia Rank Total 1.424 Unconstrained 1.424 21 Inertia 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,28,plot(gts.ca),CANOCO里面 scaling of

15、 ordination scores,29,plot(gts.ca, scaling=1),用物种数据对样方坐标进行加权平均,使样方坐标在物种数据的中心，因此对样方感兴趣的话，采用这种做图方法。,plot(gts.ca, scaling=2),plot(gts.ca, scaling=3),用样方数据对物种坐标进行加权平均,使物种数据在样方数据的中心，因此对物种感兴趣的话，采用这种做图方法。,如果一个物种靠近某个样方，表明该物种可能对该样方的位置起很大的作用。特别是对于二元数据的排序，这个样方可能就代表该物种。,如图中，20号样方与短柄枹(QUESER)靠得比较近，表明：短柄枹表征了20号样方

16、的特征，19号样方与20号样方距离近，生态关系也较近。,只在少数样方出现的物种通常在排序空间的边缘，表明它们只偶然发生，或它们只在稀有生境（如米槠CASCAR）。 在排序空间中心的物种，可能在取样区域是该物种最优分布区，如甜槠，或有两个或多个最优分布区，或与前两个轴不相关。,除趋势对应分析(Detrended correspondence analysis, DCA)：,CA采用单峰曲线表示物种和环境关系,CA产生的弓形效应,CA的第二排序轴在许多情况下是第一轴的二次变形，即所谓的“弓形效应”（Arch effect）或者“马蹄形效应”（horse-shoe effect）。 （详见张金屯群落

17、生态学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:稀有物种的权重；(稀有物种影响比较大） iresc：纠正弓形效应的次数； ira:分析的类型(DCA:0 , CA:1); mk：校正弓形效应轴的分段数; short: 需要校正的最短梯度。,plot(decorana(gtsdata),plot(cca(g

18、tsdata),直接梯度分析（Direct Gradient Analysis） ： RDA (Redundance analysis) CCA (Canonical correspondence analysis) pRDA (partial RDA) pCCA (partial CCA),冗余分析(redundancy analysis, RDA)及典范对应分析（Canonical correspondence analysis, CCA) 1.通常采用PCA处理环境数据，采用CA处理群落数据，但这些方法都只能处理一个数据表； 2.RDA和CCA是多元分析（PCA，CA）和线性回归的结合，

19、研究植被和环境之间的关系。,40,当我们在解释变量（环境因子数据）与响应变量（物种数据）之间建立预测模型的时候，经常会遇到这样的情况，往往我们仅仅考察解释变量中某几个环境因子的对物种数据的影响，但剩下的环境因子也会对物种产生影响，这些剩余环境因子我们经常称为协变量(Covariables) 。在CANOCO中，协变量的影响可以用偏分析（partial analyze）剔除出来。 实际上，任何一个环境因子变量均可以成为协变量。例如，我们要研究管理模式对蝴蝶群落中组成的影响，我们可以在不同的海拔地点取样，海拔也许对群落物种组成影响很大，但此时我们感兴趣的是管理模式的影响，而非海拔梯度的影响。这个时

20、候，如果能剔除出海拔的影响，我们能管理模型与蝴蝶种群之间更清晰的关系。 摘自 第一章6p,PCA与环境因子结合是RDA，CA与环境因子结合是CCA。 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(g

21、ts.rda) gts.prda=rda(gtsdata,gtsenv,1:4, gtsenv,5:8) summary(gts.prda) plot(gts.prda),plot(gts.rda),环境因子一般用箭头表示，箭头所处的象限表示环境因子与排序轴间的正负相关性，箭头连线的长度代表着某个环境因子与群落分布和种类分布间相关程度的大小，连线越长，说有相关性越大。反之越小。箭头连线和排序轴的夹角代表着某个环境因子与排序轴的相关性大小，夹角越小，相关性越高；反之越低。,gts.cca=cca(gtsdata,gtsenv); gts.cca=cca(gtsdataelev+convex+sl

22、ope+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),RDA和CCA结果的检验： goodness(object, display = c(“species”, “sites”), choices, model = c(“CCA”, “CA”), statistic = c(“explained”, “distance”), summari

23、ze = FALSE, .) ：物种或样方与轴累计解释量； 【Cumulative fit per species as fraction of variance of species】in CANOCO inertcomp(object, display = c(“species”, “sites”), statistic = c(“explained”, “distance”), proportional = FALSE) ：物种或样方的方差分解分析； spenvcor(object) :物种和环境的相关分析; intersetcor(object) :环境因子和各轴的相关分析； 最好使用

24、：envfit; vif.cca(object) ：方差膨胀因子分析；,RDA和CCA结果的检验： envfit(X, P, permutations = 0, strata, choices=c(1,2), .) envfit(formula, data, .) gts.cca=cca(gtsdata,gtsenv) 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 s

25、lope 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.6004984 0.7996260 0.3364 0.001 * pH -0.1929061 -0.9812172 0.6313 0.001 *,RDA和CCA结果的检验： anova(object, alpha=0.05, beta=0.01, step=100, perm.ma

26、x=10000, by = NULL, .) permutest.cca(x, permutations=100, model=c(“direct“, “reduced“,“full“), first = FALSE, strata, .) x, object:cca或rda的分析结果； by:为”axis”,或”terms”，分析各轴或因子的显著性；,如： goodness(gts.cca,display=“sp”) inertcomp(gts.cca, proportional = T) spenvcor(gts.cca) intersetcor(gts.cca) vif.cca(gts.

27、cca),RDA和CCA模型的选择： step(mod1,scope) mod1=cca(gtsdata1,data=gtsenv); mod2=cca(gtsdataelev+convex+slope+aspect+P+N+K+pH,data=gtsenv); mod=step(mod1,scope=(list(lower=formula(mod1),upper=formula(mod2); Start: AIC=180.2 gtsdata 1 Df AIC + P 1 174.66 + elev 1 174.67 + N 1 174.89 + pH 1 176.95 + slope 1 1

28、77.82 + K 1 178.14 180.20 + convex 1 180.85 + aspect 1 181.53,RDA和CCA结果的检验： gts.cca=cca(gtsdataelev+convex+slope+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)

29、 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 0.01 * aspect 1 0.0208 0.5981 100 0.55 P 1 0.0911 2.6159 100 0.01 * N 1 0.0591 1.6966 100 0.04 * K 1 0.0484 1.3906 100 0.07 . pH 1 0.0307 0.8814 100 0.39,RDA和CCA模型的选择： mod Call: cca(f

30、ormula = gtsdata P + pH + elev + convex, data = gtsenv) Inertia Rank Total 1.4243 Constrained 0.5640 4 Unconstrained 0.8603 21 Inertia is mean squared contingency coefficient Eigenvalues for constrained axes: CCA1 CCA2 CCA3 CCA4 0.31595 0.18030 0.04863 0.01913 Eigenvalues for unconstrained axes: CA1

31、 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:选择的轴； display:显示的类型，sp,物种，wa,样方，lc, 线性结合的样方坐标;,作图函数： text(x, display = “sites“, label

32、s, 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:环境因子放大倍数;,作图函数: points(x, display = “sites“, choices = c(1, 2), scaling = 2, arrow.mul, head.arrow = 0.05, select, .) plot(gts.cca,type

33、=“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(gts.cca,display=“lc”,col=4”);,绘制三维排序图： ordiplot3d(object, display = “sites“, choices = 1:3, ax.col = 2, arr.len = 0.1, arr.col = 4,

34、envfit, xlab, ylab, zlab, .) ordirgl(object, display = “sites“, choices = 1:3, type = “p“, ax.col = “red“, arr.col = “yellow“, text, envfit, .) orglpoints(object, display = “sites“, choices = 1:3, .) orgltext(object, text, display = “sites“, choices = 1:3, justify = “center“, adj = 0.5, .),ordiplot3

35、d(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, levels, labcex = 0.6, .) x:cca或rda的结果， y:需要绘图的环境因子或物种;,ef=envfit(gts.cca,gtsenv,permu

36、tation=1000) plot(gts.cca,display=“sites“) plot(ef) tmp=with(gtsenv,ordisurf(gts.cca,K,add=T) with(gtsenv,ordisurf(gts.cca,pH,add=T,col=6),绘制物种因子梯度图： ordisurf(x, y, choices=c(1, 2), knots=10, family=“gaussian“, col=“red“, thinplate = TRUE, add = FALSE, display = “sites“, w = weights(x), main, nlevel

37、s = 10, levels, labcex = 0.6, .) ordihull(ord, groups, display = “sites“, draw = c(“lines“,“polygon“), show.groups, .) ordispider(ord, groups, display=“sites“, w = weights(ord, display), show.groups, .) x:cca或rda的结果， y:需要绘图的因子; add, 添加到别的图上;,plot(gts.cca,type=“p“,display=“sites“) with(gtsdata,points

38、(gts.cca,dis=“sites“, cex=sqrt(QUESER),pch=16,col=4) with(gtsdata,ordihull(gts.cca,dis=“sites“, QUESER0,show=T) with(gtsdata,ordispider(gts.cca,QUESER0, show=T,w=QUESER,col=4) with(gtsdata,ordisurf(gts.cca, QUESER,add=T),排序结果的比较： procrustes(X, Y, scale = TRUE, symmetric = FALSE, scores = “sites“, .) protest(X, Y, scores = “sites“, permutations = 1000, strata, .) x,y: 两个排序结果；,排序结果的比较： gts.cca=cca(gtsdata,gtsenv); gts.rda=rda(gtsdata,gtsenv); procrustes(gts.cca,gts.rda,scores=“sp”); protest(gts.cca,gts.rda, scores=“sp”);,致谢,本课件基于米湘成博士R培训ppt基础上修改，在此表示致谢！,持之以恒=R的高手！,谢谢，请批评指正！,