用vegan进行排序分析 米湘成.

Presentation on theme: "用vegan进行排序分析 米湘成."— Presentation transcript:

1 用vegan进行排序分析 米湘成

2 一、vegan软件包 >library(vegan) 是Vegetation analysis的缩写。是一个群落分析的软件。
作者：Jari Oksanen >library(vegan)

3 >gtsdata=read.table(“gtsdata.txt”,header=T);
群落数据 >gtsdata=read.table(“gtsdata.txt”,header=T); >gtsdata >dim(gtsdata) CASEYR SCHSUP DAPOLD PINMAS CINSUB CYCGLA MACTHU LITGLA CASFAR MELOLD MYRRUB ADIMIL ALBKAL SORFOL CYCMYR TOXSUC CASTIB

4 elev convex slope aspect N P K pH
环境数据 >gtsenv=read.table(“gtsenv.txt”,header=T); >gtsenv >dim(gtsenv) elev convex slope aspect N P K pH

5 1. decostand(x, method, MARGIN, range.global, na.rm = FALSE)
数据的标准化 1. decostand(x, method, MARGIN, range.global, na.rm = FALSE) total: 除以行和或列和 (default MARGIN = 1); max:除以行或列的最大值(default MARGIN = 2); freq:除以行或列的最大值,并乘以非零值的个数(default MARGIN=2); normalize:使行或列的平方和等于1 (default MARGIN = 1); range: 标准化使行或列的值在 (default MARGIN = 2). standardize:标准化使行或列的和为1且方差为1(default MARGIN = 2)； pa: 将数据转换为0、1数据； chi.square: 除以行和及列和的平方根； hellinger: 采用total标准化以后再取平方根； 2. wisconsin(x) :除以列最大值，再除以行和。

6 一、什么是排序(ordination)? 间接梯度分析 (Indirect gradient analysis)

7 一、什么是排序(ordination)? 间接梯度分析 (Indirect gradient analysis)

8 3个种的排序图 2个种的排序图 4个种的排序图？？？ 40个种排序图？？？

9 排序的内容： 1.降低维数，减少坐标轴的数目 ；
2. 由降低维数引起的信息损失尽量少，即发生最小的畸变，也就是说它的第1-3轴排序轴包含大量的生态信息 。

10 排序的目标： 表示植被与环境之间的关系： 所有排序方法都反映植物种和环境之间的关系以及在某一环境梯度上的种间关系。
线形模型（linear model），短的梯度，主成分分析（Principle component analysis)，需要对数据进行非线性转换，如取对数； 非线性模型（non-linear model）如高斯模型，长的梯度，对应分析 （Correspondence analysis)

11 排序采用的距离方法： 欧氏距离： 2.卡方距离

12 二、主成分分析(Principle component analysis, PCA)
主成分分析的主要原理是： 使坐标旋转一定的角度后，使第一轴表示数据最大的方差，使第二轴表示数据第二的方差。而且轴与轴之间是正交的（orthogonal)。

13 rda(formula, data, scale=FALSE, ...) rda(X, Y, Z, scale=FALSE, ...)
PCA和RDA都采用函数rda实现： 在vegan包中， rda(formula, data, scale=FALSE, ...) rda(X, Y, Z, scale=FALSE, ...) scores(x, display=c("sites", "species"), choices, ...) 在stat包中： princomp(x, ...) princomp(formula, data = NULL, subset, na.action, ...)

14 主成分分析表示信息的方法： 特征值(eigenvalue)和Inertia 每一个轴都有一个特征值表示其信息量也就是方差的大小。

15 >sum(apply(gtsdata,2,var)) >summary(gts.rda)
>gts.rda=rda(gtsdata) >gts.rda Call: rda(X = gtsdata) Inertia Rank Total Unconstrained Inertia is variance Eigenvalues for unconstrained axes: PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 >sum(apply(gtsdata,2,var)) >summary(gts.rda) !!所以如果不对数据做标准化的话，丰富种的值就非常大，排序时就只能看清丰富种的位置，其它种就拥挤在一起。

16 plot(rda(gtsdata,scale=T))

17 2. 特征值向量(eigenvector)和负荷(loading)
为特征向量， 为特征根, 为负荷;

18 >gts.pca=princomp(gtsdata); >summary(gts.pca)
Importance of components: Comp.1 Comp.2 Comp.3 Comp.4 Comp Comp Comp Comp Comp.9 Standard deviation Proportion of Variance Cumulative Proportion Comp Comp Comp Comp Comp Comp Comp Comp.17 Standard deviation Proportion of Variance Cumulative Proportion Comp Comp Comp Comp Comp.22 Standard deviation e-01 Proportion of Variance e-05 Cumulative Proportion e+00

19 >gts.pca$loadings Loadings:
Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9 Comp.10 Comp.11 Comp.12 Comp.13 Comp.14 Comp.15 Belper Empnig Junbuf Junart Airpra Elepal Rumace Viclat Brarut Ranfla Cirarv Hyprad Leoaut Potpal Poapra Calcus Tripra Trirep Antodo Comp.16 Comp.17 Comp.18 Comp.19 Belper Empnig ……

20 Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9 Comp.10
>dune.pca$scores Comp Comp Comp Comp Comp Comp Comp Comp Comp Comp.10

21 >biplot(gts.pca) 箭头所指的方向就是旋转以前物种的轴，箭头的长短表示物种对排序轴的影响大小。

22 >princomp(gtsdata,cor=T) >rda(gtsdata, scale=T)
我们常用的PCA是采用协方差矩阵进行计算的，但是也可以采用相关矩阵进行计算。 >princomp(gtsdata,cor=T) >rda(gtsdata, scale=T) 采用相关矩阵进行计算要使各个种之间变得差异小一些，两者结果是不同的。在实际应用中，两个都进行计算，选更易于解释和表现生态关系的。如果要在排序图上表现聚类，可能要采用协方差计算结果；如果不同的种有类似的贡献时，可能要采用相关矩阵计算结果。

23 如果前三轴包含的信息量非常少，表明种间的信息相关性非常少，信息冗余非常少，不能用PCA分析。
3. 每个轴包含的信息量 如果前三轴包含的信息量非常少，表明种间的信息相关性非常少，信息冗余非常少，不能用PCA分析。 >summary(gts.rda) Eigenvalues, and their contribution to the variance PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 PC12 PC13 PC14 PC15 lambda accounted PC16 PC17 PC18 PC19 PC20 PC21 PC22 lambda accounted

24 二、对应分析(Correspondence analysis, CA)
1. CA采用单峰模型，因此结果一般优于PCA； 2. CA采用 距离，因此用于处理有很多0的数据，特别是处理“有”或“无”的数据； 3. 矩阵中不允许有负数，PCA分析矩阵可有负数。

25 4. CA除了采用卡方距离以外，其它计算步骤都与PCA一致（如旋转等)
5. CA在vegan中的实现： cca(formula, data, ...) cca(X, Y, Z, ...)

26 >gts.ca=cca(gtsdata) >gts.ca
4. CA除了采用卡方距离以外，其它计算步骤都与PCA一致（如旋转等） >gts.ca=cca(gtsdata) >gts.ca Call: cca(X = gtsdata) Inertia Rank Total Unconstrained Inertia is mean squared contingency coefficient Eigenvalues for unconstrained axes: CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8

27 >chisq.test(gtsdata/sum(gtsdata))
4. Inerta在这里是矩阵标准化后的卡方统计 >chisq.test(gtsdata/sum(gtsdata)) Pearson's Chi-squared test data: gtsdata/sum(gtsdata) X-squared = , df = 819, p-value = 1 >summary(gts.ca) 5. 对样方和物种的排序坐标进行加权平均。（1）用物种数据对样地坐标进行加权平均；（2）用样方数据对物种坐标进行加权平均；（3）对样方和物种坐标都进行加权平均；加权平均后的排序图会收缩，收缩的因子等于特征值。

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

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

30 6. 对CA结果的解释： （1）每个轴的特征值表示各轴样方坐标和物种坐标的相关性，如：例子中第一轴的特征值是： 。 (2) 采用物种数据对样方坐标进行加权平均时， A. 样方之间的距离就是它们之间的卡方距离; B. 由于加权平均时，采用样方内物种的相对多度进行加权，因此如果一个样方内没有某个物种，这个物种对样方的坐标没有贡献。

31 C. 因此如果一个物种靠近某个样方，表明该物种可能对该样方的位置起很大的作用。特别是对于二元数据的排序，这个样方可能就代表该物种。
如图中，20号样方与短柄枹(QUESER)靠得比较近，表明：短柄枹表征了20号样方的特征，19号样方与20号样方距离近，生态关系也较近。

32 （3）生态位的阐述： 物种是按单峰曲线沿生境梯度分布的， 物种在排序空间的分布表明物种和环境的关系，物种和样方的相对位置表明物种和样方的相对关系。 如图中，22号样方和青冈(CYCGLA)

33 B. 只在少数样方出现的物种通常在排序空间的边缘，表明它们只偶然发生，或它们只在稀有生境（如米槠CASCAR）。

34 除趋势对应分析(Detrended correspondence analysis, DCA)：

35 CA的第二排序轴在许多情况下是第一轴的二次变形，即所谓的“弓形效应”（Arch effect）或者“马蹄形效应”（horse-shoe effect）。

36 DCA在R中的实现采用函数decorana。
decorana(veg, iweigh=0, iresc=4, ira=0, mk=26, short=0, before=NULL, after=NULL) veg:群落数据; iweig:稀有物种的权重； iresc：纠正弓形效应的次数； ira:分析的类型(0: DCA, 1:CA); mk：校正弓形效应轴的分段数; short: 需要校正的最短梯度。

37 plot(decorana(gtsdata))
plot(cca(gtsdata)) plot(decorana(gtsdata))

38 冗余分析(redundancy analysis, RDA)及典范对应分析（Canonical correspondence analysis, CCA)
1.通常采用PCA处理环境数据，采用CA处理群落数据，但这些方法都只能处理一个数据表； 2.RDA和CCA是多元分析（PCA，CA）和线性回归的结合，研究植被和环境之间的关系。

39 CCA在vegan中的实现： cca(formula, data, ...) cca(X, Y, Z, ...)
PCA与环境因子结合是RDA，CA与环境因子结合是CCA。 CCA在vegan中的实现： cca(formula, data, ...) cca(X, Y, Z, ...) RDA在vegan中的实现： rda(formula, data, ...) rda(X, Y, Z, ...)

40 ＞gts.rda=rda(gtsdata,gtsenv);
>gts.rda=rda(gtsdata~elev+convex+slope+aspect+N+K+P+pH,data=gtsenv); scale=F scale=T

41 >gts.rda Call: rda(X = gtsdata, Y = gtsenv) Inertia Rank
Total Constrained Unconstrained Inertia is variance Eigenvalues for constrained axes: RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7 RDA8 Eigenvalues for unconstrained axes: PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8

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

43 >summary(gts.rda) Call: rda(X = gtsdata, Y = gtsenv) Inertia Rank
Total Constrained Unconstrained Inertia is variance Eigenvalues for constrained axes: RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7 RDA8 Eigenvalues for unconstrained axes: PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8

44 >gts.cca=cca(gtsdata,gtsenv);
>gts.cca=cca(gtsdata~elev+convex+slope+aspect+N+K+P+pH,data=gtsenv) plot(gts.ccca)

45 >gts.cca Call: cca(X = gtsdata, Y = gtsenv) Inertia Rank
Total Constrained Unconstrained Inertia is mean squared contingency coefficient Eigenvalues for constrained axes: CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7 CCA8 Eigenvalues for unconstrained axes: CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8

46 plot(gts.cca)

47 RDA和CCA结果的检验： goodness(object, display = c(“species”, “sites”), choices, model = c(“CCA”, “CA”), statistic = c(“explained”, “distance”), summarize = FALSE, ...) ：检验物种或样方与环境因子的相关性； inertcomp(object, display = c(“species”, “sites”), statistic = c(“explained”, “distance”), proportional = FALSE) ：物种或样方的方差分解分析； spenvcor(object) :物种和环境的相关分析; intersetcor(object) :物种和各轴的相关分析； 最好使用：envfit; vif.cca(object) ：方差膨胀因子分析；

48 如： >goodness(gts.cca,display=“sp”) >inertcomp(gts.cca, proportional = T) >spenvcor(gts.cca) >intersetcor(gts.cca) >vif.cca(gts.cca)

49 display:显示的类型，sp,物种，wa,样方，lc, 线性结合的样方坐标;
作图函数： plot(x, choices = c(1, 2), display = c("sp", "wa", "cn"), scaling = 2, type, xlim, ylim, ...) x: cca或rda分析结果； choices:选择的轴； display:显示的类型，sp,物种，wa,样方，lc, 线性结合的样方坐标;

50 display:显示的类型，sp，wa, lc，bp,环境因子;cn,中心化后的环境因子；
作图函数： 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:环境因子放大倍数;

51 >plot(gts.cca,type="n");
作图函数: 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(gts.cca,display=“lc”,col=4”);

52 envfit(X, P, permutations = 0, strata, choices=c(1,2), ...)
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) CCA CCA2 r2 Pr(>r) elev <0.001 *** convex slope <0.001 *** aspect N <0.001 *** P <0.001 *** K <0.001 *** pH <0.001 ***

53 RDA和CCA结果的检验： anova(object, alpha=0.05, beta=0.01, step=100, perm.max=10000, by = NULL, ...) permutest.cca(x, permutations=100, model=c("direct", "reduced","full"), first = FALSE, strata, ...) x, object:cca或rda的分析结果； by:为”axis”,或”terms”，分析各轴或因子的显著性；

54 Df Chisq F N.Perm Pr(>F)
RDA和CCA结果的检验： >gts.cca=cca(gtsdata~elev+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) Df Chisq F N.Perm Pr(>F) elev <0.01 *** convex slope <0.01 *** aspect P <0.01 *** N * K pH

55 >mod1=cca(gtsdata~1,data=gtsenv);
RDA和CCA模型的选择： step(mod1,scope) >mod1=cca(gtsdata~1,data=gtsenv); >mod2=cca(gtsdata~elev+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 + elev + N + pH + slope + K <none> + convex + aspect

56 RDA和CCA模型的选择： >mod Call:
cca(formula = gtsdata ~ P + pH + elev + convex, data = gtsenv) Inertia Rank Total Constrained Unconstrained Inertia is mean squared contingency coefficient Eigenvalues for constrained axes: CCA1 CCA2 CCA3 CCA4 Eigenvalues for unconstrained axes: CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8

57 绘制环境因子梯度图： x:cca或rda的结果， y:需要绘图的环境因子或物种;
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:需要绘图的环境因子或物种;

58 >ef=envfit(gts.cca,gtsenv,permutation=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))

59 绘制物种因子梯度图： x:cca或rda的结果， y:需要绘图的因子; add, 添加到别的图上;
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, ...) 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, 添加到别的图上;

60 >plot(gts.cca,type="p",display="sites")
> with(gtsdata,points(gts.cca,dis="sites", cex=sqrt(QUESER),pch=16,col=4)) > with(gtsdata,ordihull(gts.cca,dis="sites", QUESER>0,show=T)) >with(gtsdata,ordispider(gts.cca,QUESER>0, show=T,w=QUESER,col=4)) >with(gtsdata,ordisurf(gts.cca, QUESER,add=T)

61 绘制三维排序图： 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 = "sites", choices = 1:3, ...) orgltext(object, text, display = "sites", choices = 1:3, justify = "center", adj = 0.5, ...)

62 >ordiplot3d(gts.cca,type=“h”) > ordirgl(gts.cca,type="p")

63 排序结果的比较： procrustes(X, Y, scale = TRUE, symmetric = FALSE, scores = "sites", ...) protest(X, Y, scores = "sites", permutations = 1000, strata, ...) x,y: 两个排序结果；

64 排序结果的比较： >gts.cca=cca(gtsdata,gtsenv); >gts.rda=rda(gtsdata,gtsenv); >procrustes(gts.cca,gts.rda,scores=“sp”); >protest(gts.cca,gts.rda, scores=“sp”);

65 THANK YOU ! 谢谢！