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

1,基于的群落多元统计分析用vegan包进行排序分析,赖江山（janson)lai中国科学院植物研究所2010.11.5,2,3,vegan软件包简介,vegan是Vegetationanalysis的缩写,是群落分析的package作者：JariOksanen/web/packages/vegan/index.htmlhttp:/cc.oulu.fi/jarioksa/softhelp/vegan.htmllibrary(vegan),什么是排序(ordination)?排序的过程是将样方或植物种排列在一定的空间，使得排序轴能够反映一定的生态梯度，从而，能够解释植被或植物种的分布与环境因子间的关系，也就是说排序是为了揭示植被-环境间的生态关系。因此，排序也叫梯度分析（gradientanalysis）。间接梯度分析(Indirectgradientanalysis)直接梯度分析(directgradientanalysis),2个种的排序图,3个种的排序图,4个种的排序图？,40个种排序图？,排序的目标：1.降低维数，减少坐标轴的数目；2.由降低维数引起的信息损失尽量少，即发生最小的畸变，也就是让新的坐标系第1-3轴排序轴包含大量的生态信息。,排序的目的：表示植被与环境之间的关系：所有排序方法都反映植物种和环境之间的关系以及在某一环境梯度上的种间关系。线形模型（linearmodel），短的梯度，主成分分析（Principlecomponentanalysis)，需要对数据进行非线性转换，如取对数；非线性模型（non-linearmodel）如高斯模型，长的梯度，对应分析（Correspondenceanalysis),群落数据输入,gtsdata=read.table(gtsdata.txt,header=T)gtsdatadim(gtsdata),环境因子数据输入,gtsenv=read.table(gtsenv.txt,header=T);gtsenvdim(gtsenv),数据的标准化1.decostand(x,method,MARGIN,)total:除以行和或列和(defaultMARGIN=1是row);max:除以行或列的最大值(defaultMARGIN=2是列);freq:除以行或列的最大值,并乘以非零值的个数，非零值的平均值为1(defaultMARGIN=2);normalize:使行或列的平方和等于1(defaultMARGIN=1);range:标准化使行或列的值在0.1(defaultMARGIN=2).standardize:标准化使行或列的和为1且方差为1(defaultMARGIN=2)；pa:将数据转换为0、1数据；chi.square:除以行和及列和的平方根；hellinger:采用total标准化以后再取平方根；log：对数化，默认自然对数，logbase参数是自选的base2.wisconsin(x):除以列最大值，再除以行和。,排序类别(inCANOCO)间接梯度分析（IndirectGradientAnalysis）：PCA(Principalcomponentsanalysis)CA(Correspondenceanalysis)DCA(DetrendedCorrespondenceAnalysis)直接梯度分析（DirectGradientAnalysis）：RDA(Redundanceanalysis)CCA(Canonicalcorrespondenceanalysis)DCCA(DetrendedCCA),PCARDA,CACCADCADCCA,13,决定排序的模型：单峰还是线性？decorana(gtsdata)Call:decorana(veg=gtsdata)Detrendedcorrespondenceanalysiswith26segments.Rescalingofaxeswith4iterations.DCA1DCA2DCA3DCA4Eigenvalues0.39390.22390.095550.06226Decoranavalues0.50250.17560.067120.03877Axislengths3.25952.51301.214451.00854,如果这四个轴中梯度最长（最大值）超过4，选择单峰模型排序（CA、CCA、DCA）更合适。如果是小于3，选择线性模型（PCA、RDA）比较合理。如果介于3-4之间，单峰模型和线性模型结果差不多。,间接梯度分析（IndirectGradientAnalysis）PCA(Principalcomponentsanalysis)CA(Correspondenceanalysis)DCA(DetrendedCorrespondenceAnalysis),主成分分析(Principlecomponentanalysis,PCA),主成分分析的主要原理是：使坐标旋转一定的角度后，使第一轴表示数据最大的方差，使第二轴表示数据第二的方差。而且轴与轴之间是正交的（orthogonal)。,PCA和RDA都采用函数rda实现：在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.rdaCall:rda(X=gtsdata)InertiaRankTotal352.1Unconstrained352.122InertiaisvarianceEigenvaluesforunconstrainedaxes:PC1PC2PC3PC4PC5PC6PC7PC8111.77973.58054.60732.95926.48118.06312.7637.637scores(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(rda(gtsdata),!如果不对数据做标准化的话，丰富种的值就非常大，排序时就只能看清丰富种的位置，其它种就拥挤在一起。,如用x1,x2,x3,x4,x5,x6分别表示原先的变量，而用y1,y2,y3,y4,y5,y6表示新的主成分，那么，第一和第二主成分为,这些系数称为主成分载荷（loading），它表示主成分和相应的原先变量的相关系数。比如y1表示式中x1的系数为-0.806，这就是说第一主成分和的x1变量的相关系数为-0.806。相关系数(绝对值）越大，主成分对该变量的代表性也越大。,负荷(loading),gts.pca=princomp(gtsdata)gts.pca$loadings,gtsenv.pca=princomp(gtsdenv)gtsenv.pca$loadings,biplot(gtsenv.pca),第一主成分代表海拔高度，第二主成分代表坡向,对应分析(Correspondenceanalysis,CA)PCA在迭代运算过程是采用线性模型CA在迭代运算过程采用单峰模型（加权平均法）,CA在vegan中也是用cca函数来实现：gts.ca=cca(gtsdata)gts.casummary(gts.ca)Call:cca(X=gtsdata)InertiaRankTotal1.424Unconstrained1.42421InertiaismeansquaredcontingencycoefficientEigenvaluesforunconstrainedaxes:CA1CA2CA3CA4CA5CA6CA7CA80.502530.265640.140230.105020.091270.055400.050630.04204,26,plot(gts.ca),27,CANOCO里面scalingofordinationscores,plot(gts.ca,scaling=1),用物种数据对样方坐标进行加权平均,使样方坐标在物种数据的中心，因此对样方感兴趣的话，采用这种做图方法。,plot(gts.ca,scaling=2),plot(gts.ca,scaling=3),用样方数据对物种坐标进行加权平均,使物种数据在样方数据的中心，因此对物种感兴趣的话，采用这种做图方法。,如果一个物种靠近某个样方，表明该物种可能对该样方的位置起很大的作用。特别是对于二元数据的排序，这个样方可能就代表该物种。,如图中，20号样方与短柄枹(QUESER)靠得比较近，表明：短柄枹表征了20号样方的特征，19号样方与20号样方距离近，生态关系也较近。,只在少数样方出现的物种通常在排序空间的边缘，表明它们只偶然发生，或它们只在稀有生境（如米槠CASCAR）。在排序空间中心的物种，可能在取样区域是该物种最优分布区，如甜槠，或有两个或多个最优分布区，或与前两个轴不相关。,除趋势对应分析(Detrendedcorrespondenceanalysis,DCA)：,CA采用单峰曲线表示物种和环境关系,CA产生的弓形效应,CA的第二排序轴在许多情况下是第一轴的二次变形，即所谓的“弓形效应”（Archeffect）或者“马蹄形效应”（horse-shoeeffect）。（详见张金屯群落生态学168页）,和分别代表除趋势前和除趋势后的样方排序（引自Hill和Gauch1980）,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(gtsdata),直接梯度分析（DirectGradientAnalysis）：RDA(Redundanceanalysis)CCA(Canonicalcorrespondenceanalysis)pRDA(partialRDA)pCCA(partialCCA),冗余分析(redundancyanalysis,RDA)及典范对应分析（Canonicalcorrespondenceanalysis,CCA)1.通常采用PCA处理环境数据，采用CA处理群落数据，但这些方法都只能处理一个数据表；2.RDA和CCA是多元分析（PCA，CA）和线性回归的结合，研究植被和环境之间的关系。,38,当我们在解释变量（环境因子数据）与响应变量（物种数据）之间建立预测模型的时候，经常会遇到这样的情况，往往我们仅仅考察解释变量中某几个环境因子的对物种数据的影响，但剩下的环境因子也会对物种产生影响，这些剩余环境因子我们经常称为协变量(Covariables)。在CANOCO中，协变量的影响可以用偏分析（partialanalyze）剔除出来。实际上，任何一个环境因子变量均可以成为协变量。例如，我们要研究管理模式对蝴蝶群落中组成的影响，我们可以在不同的海拔地点取样，海拔也许对群落物种组成影响很大，但此时我们感兴趣的是管理模式的影响，而非海拔梯度的影响。这个时候，如果能剔除出海拔的影响，我们能管理模型与蝴蝶种群之间更清晰的关系。摘自第一章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(gts.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+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),RDA和CCA结果的检验：goodness(object,display=c(“species”,“sites”),choices,model=c(“CCA”,“CA”),statistic=c(“explained”,“distance”),summarize=FALSE,.)：物种或样方与轴累计解释量；【Cumulativefitperspeciesasfractionofvarianceofspecies】inCANOCOinertcomp(object,display=c(“species”,“sites”),statistic=c(“explained”,“distance”),proportional=FALSE)：物种或样方的方差分解分析；spenvcor(object):物种和环境的相关分析;intersetcor(object):环境因子和各轴的相关分析；最好使用：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)CCA1CCA2r2Pr(r)elev-0.9999873-0.00504790.50050.001*convex-0.77465040.63238980.08440.28slope0.15266850.98827750.52960.001*aspect-0.1364818-0.99064260.01690.82N0.68789420.72581100.63420.001*P0.81015020.58622240.58560.001*K0.60049840.79962600.3364inertcomp(gts.cca,proportional=T)spenvcor(gts.cca)intersetcor(gts.cca)vif.cca(gts.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.2gtsdata1DfAIC+P1174.66+elev1174.67+N1174.89+pH1176.95+slope1177.82+K1178.14180.20+convex1180.85+aspect1181.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=gtsdataelev+convex+slope+aspect+P+N+K+pH,data=gtsenv)DfChisqFN.PermPr(F)elev10.24467.02281000.01*convex10.06241.79021000.07.slope10.13563.8928100plot(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”);,绘制三维排序图：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,.),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,levels,labcex=0.6,.)x:cca或rda的结果，y:需要绘图的环境因子或物种;,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),绘制物种因子梯度图：ordisur