R语言vegan包使用教程

更新时间：2023-07-18 22:13:01 阅读量：实用文档文档下载

说明：文章内容仅供预览，部分内容可能不全。下载后的文档，内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的，是否完整无缺。

R语言中的外在软件包“Vegan”是专门用于群落生态学分析的工具。Vegan能够提供所有基本的排序方法,同时具有生成精美排序图的功能,版本更新很快。我们认为Vegan包完全可以取代CANOCO,成为今后排序分析的首选统计工具。

MultivariateAnalysisofEcologicalCommunitiesinR:vegantutorial

JariOksanenFebruary8,2013

Abstract

ThistutorialdemostratestheuseofordinationmethodsinRpack-agevegan.ThetutorialassumesfamiliaritybothwithRandwithcommunityordination.Packagevegansupportsallbasicor-dinationmethods,includingnon-metricmultidimensionalscaling.Theconstrainedordinationmethodsincludeconstrainedanalysisofproximities,redundancyanalysisandconstrainedcorrespondenceanalysis.Packageveganalsohassupportfunctionsfor ttingen-vironmentalvariablesandforordinationgraphics.

Contents

1Introduction

2Ordination:basicmethod

2.1Non-metricMultidimensionalscaling............2.2Communitydissimilarities..................2.3Comparingordinations:Procrustesrotation........2.4Eigenvectormethods.....................2.5Detrendedcorrespondenceanalysis.............2.6Ordinationgraphics.............

........

3Environmentalinterpretation

3.1Vector tting.........................3.2Surface tting.........................3.3Factors.............................4Constrainedordination

4.1Modelspeci cation......................4.2Permutationtests.......................4.3Modelbuilding........................4.4Linearcombinationsandweightedaverages........4.5Biplotarrowsandenvironmentalcalibration........4.6Conditionedorpartialmodels........

........

23358811121414151618192123282930

1INTRODUCTION

5Dissimilaritiesandenvironment

5.1adonis:MultivariateANOVAbasedondissimilarities5.2Homogeneityofgroupsandbetadiversity......5.3Manteltest.......................5.4Protest:Procrustestest................

............

323233353636363839

6Classi cation

6.1Clusteranalysis........................6.2Displayandinterpretationofclasses............6.3Classi edcommunitytables.................

1Introduction

Thistutorialdemonstratestypicalwork owsinmultivariateordinationanalysisofbiologicalcommunities.Thetutorial rstdiscussesbasicun-constrainedanalysisandenvironmentalinterpretationoftheirresults.Thenitintroducesconstrainedordinationusingconstrainedcorrespon-denceanalysisasanexample:alternativemethodssuchasconstrainedanalysisofproximitiesandredundancyanalysiscanbeused(almost)similarly.Finallythetutorialdescribesanalysisofspecies–environmentrelationswithoutordination,andbrie ytouchesclassi cationofcommu-nities.

Theexamplesinthistutorialaretested:ThisisaSweavedocument.Theoriginalsource lecontainsonlytextandRcommands:theiroutputandgraphicsaregeneratedwhilerunningthesourcethroughSweave.However,youmayneedarecentversionofvegan.Thisdocumentwasgeneretatedusingveganversion2.0-6andRversion2.15.1(2012-06-22).

Themanualcoversordinationmethodsinvegan.Itdoesnotdis-cussmanyothermethodsinvegan.Forinstance,thereareseveralfunc-tionsforanalysisofbiodiversity:diversityindices(diversity,renyi,fisher.alpha),extrapolatedspeciesrichness(specpool,estimateR),speciesaccumulationcurves(specaccum),speciesabundancemodels(rad-fit,fisherfit,prestonfit)etc.NeitherisvegantheonlyRpack-ageforecologicalcommunityordination.BaseRhasstandardstatisticaltools,labdsvcomplementsveganwithsomeadvancedmethodsandpro-videsalternativeversionsofsomemethods,andade4providesanalter-nativeimplementationforthewholegammeofordinationmethods.

Thetutorialexplainsonlythemostimportantmethodsandshowstypicalwork ows.Iseeordinationprimarilyasagraphicaltool,andIdonotshowtoomuchexactnumericalresults.Instead,therearesmallvignettesofplottingresultsinthemarginsclosetotheplacewhereyouseeaplotcommand.Isuggestthatyourepeattheanalysis,trydi erentalternativesandinspecttheresultsmorethoroughlyatyourleisure.Thefunctionsareexplainedonlybrie y,anditisveryusefultocheckthecor-respondinghelppagesforamorethoroughexplanationofmethods.Themethodsalsoareonlybrie yexplained.Itisbesttoconsultatextbookonordinationmethods,ormylectures,for rmertheoreticalbackground.

2ORDINATION:BASICMETHOD

Ordination:basicmethod

2.1

Non-metricMultidimensionalscaling

Non-metricmultidimensionalscalingcanbeperformedusingisoMDSfunc-tionintheMASSpackage.Thisfunctionneedsdissimilaritiesasinput.Functionvegdistinvegancontainsdissimilaritieswhicharefoundgoodincommunityecology.ThedefaultisBray-Curtisdissimilarity,nowadaysoftenknownasSteinhausdissimilarity,orinFinlandasSørensenindex.Thebasicstepsare:

>library(vegan)>library(MASS)>data(varespec)

>vare.dis<-vegdist(varespec)>

vare.mds0<-isoMDS(vare.dis)

initialvalue18.026495iter5value10.095483finalvalue10.020469converged

Thedefaultisto ndtwodimensionsandusemetricscaling(cmdscale)asthestartingsolution.Thesolutionisiterative,ascanbeseenfromthetracinginformation(thiscanbesuppressedsettingtrace=F).

TheresultsofisoMDSisalist(itemspoints,stress)forthecon- gurationandthestress.StressSisastatisticofgoodnessof t,anditisafunctionofandnon-linearmonotonetransformationofobserved

dissimilaritiesθ(d)andordinationdistancesd

.nationNmdsspacemapsandobserveditcanhandlecommunitynonlineardissimilaritiesspeciesresponsesnonlinearlyofanyontoshape.ordi-WecaninspectthemappingusingfunctionShepardinMASSpackage,orasimplewrapperstressplotinvegan:

>stressplot(vare.mds0,vare.dis)

FunctionstressplotdrawsaShepardplotwhereordinationdistancesareplottedagainstcommunitydissimilarities,andthe tisshownasamonotonestepline.Inaddition,stressplotshowstwocorrelationlikestatisticsofgoodnessof t.ThecorrelationbasedonstressisR2=1 S2.The“ t-basedR2”isordinationdistancesd

thecorrelationbetweenthe ttedvaluesθ(d)and

,orbetweenthesteplineandthepoints.Thisshouldbelinearevenwhenthe tisstronglycurvedandisoftenknownasthe“linear t”.ThesetwocorrelationsarebothbasedontheresidualsintheShepardplot,buttheydi erintheirnullmodels.Inlinear t,thenullmodelisthatallordinationdistancesareequal,andthe tisa athorizontalline.Thissoundssensible,butyouneedN 1dimensionsforthenullmodelofNpoints,andthisnullmodelisgeometricallyimpossi-bleintheordinationspace.Thebasicstressusesthenullmodelwhereallobservationsareputinthesamepoint,whichisgeometricallypossible.Finallyawordofwarning:yousometimesseethatpeopleusecorrelationbetweencommunitydissimilaritiesandordinationdistances.Thisisdan-gerousandmisleadingsincenmdsisanonlinearmethod:animproved

S= i=j[θ(dij) d ij]2i=jd2

Non metric fit, R2 = 0.99 q

.Linear fit, R2 = 0.943

1qqqqqqq8

qqq

qqe

0qcqnqqqqqqqqqaqqqqqqqqqqqqqqqtsqqqiqqqqqqqqqqqqD6qqqq

qqqqqq .n0qqoqqqqqqqqqqqqqqiqqqqqtanqqqqqqqqqqqqqqqqqqqqqidqqqqqqqqqqqqqqqqqqqqqqrO4q.qq0qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq

qqqqqqqqqqqqqqqqqqqqqqqqqqqqqq2qqqqqqqqqqqqqqqqqqqqq.0qqqqqqqq

qqqqqqqqq

qqqqqqqqqqqqqqqqqqq

0.20.40.60.8

Observed Dissimilarity

2.1Non-metricMultidimensionalscaling2ORDINATION:BASICMETHOD

.014

.016

i.D023

.0 27

109

4.0 21

0.6 0.4 0.2

0.00.20.4

Dim1

Cla.phy.0Cet.islDic.pol

Cla.ste24

Cla.chl

Pin.syl12Bar.lyc

Cla.cer

10Poh.nut

Pti.cil

Cla.bot

Dic.sp

Emp.nig

Vac.vit0

.04Cla.ran19Pol.pilCet.eri6Cla.criCla.unc23Cla.graPel.aphCet.niv

Cal.vul13Cla.corCla.def20Ple.sch28Bet.pub2

S318Cla.fimCla.sp15

Vac.myrLed.pal16DCla.arbCla.coc22M7Pol.jun27

N14

Hyl.spl5Ste.sp

Dip.mon

Dic.fusDes.fle5

Vac.uli

.25

0 Cla.ama

Ich.eri

Nep.arc

0.50.00.51.0

NMDS1

ordinationwithmorenonlinearrelationshipwouldappearworsewiththiscriterion.

Functionsscoresandordiplotinvegancanbeusedtohandletheresultsofnmds:

>ordiplot(vare.mds0,type="t")

Onlysitescoreswereshown,becausedissimilaritiesdidnothaveinfor-mationaboutspecies.

Theiterativesearchisverydi cultinnmds,becauseofnonlinearre-lationshipbetweenordinationandoriginaldissimilarities.Theiterationeasilygetstrappedintolocaloptimuminsteadof ndingtheglobalop-timum.Thereforeitisrecommendedtouseseveralrandomstarts,andselectamongsimilarsolutionswithsmalleststresses.Thismaybete-dious,butveganhasfunctionmetaMDSwhichdoesthis,andmanymorethings.Thetracingoutputislong,andwesuppressitwithtrace=0,butnormallywewanttoseethatsomethinghappens,sincetheanalysiscantakealongtime:

>vare.mds<-metaMDS(varespec,trace=FALSE)>vare.mds

Call:

metaMDS(comm=varespec,trace=FALSE)

globalMultidimensionalScalingusingmonoMDSData:wisconsin(sqrt(varespec))Distance:bray

Dimensions:2Stress:0.1826

Stresstype1,weakties

Twoconvergentsolutionsfoundafter20tries

Scaling:centring,PCrotation,halfchangescaling

Species:expandedscoresbasedon‘wisconsin(sqrt(varespec))’>plot(vare.mds,type="t")

Wedidnotcalculatedissimilaritiesinaseparatestep,butwegavethe

originaldatamatrixasinput.Theresultismorecomplicatedthanpre-viously,andhasquiteafewcomponentsinadditiontothoseinisoMDSre-sults:nobj,nfix,ndim,ndis,ngrp,diss,iidx,jidx,xinit,is-tart,isform,ities,iregn,iscal,maxits,sratmx,strmin,sf-grmn,dist,dhat,points,stress,grstress,iters,icause,call,model,distmethod,distcall,data,distance,converged,tries,engine,species.Thefunctionwrapsrecommendedproceduresintoonecommand.Sowhathappenedhere?

1.Therangeofdatavalueswassolargethatthedataweresquareroottransformed,andthensubmittedtoWisconsindoublestandardiza-tion,orspeciesdividedbytheirmaxima,andstandsstandardizedtoequaltotals.Thesetwostandardizationsoftenimprovethequal-ityofordinations,butweforgottothinkaboutthemintheinitialanalysis.

2ORDINATION:BASICMETHOD2.FunctionusedBray–Curtisdissimilarities.

3.FunctionrunisoMDSwithseveralrandomstarts,andstoppedei-therafteracertainnumberoftries,orafter ndingtwosimilarcon gurationswithminimumstress.Inanycase,itreturnedthebestsolution.4.Functionrotatedthesolutionsothatthelargestvarianceofsitescoreswillbeonthe rstaxis.5.Functionscaledthesolutionsothatoneunitcorrespondstohalvingofcommunitysimilarityfromthereplicatesimilarity.6.Functionfoundspeciesscoresasweightedaveragesofsitescores,butexpandedthemsothatspeciesandsitescoreshaveequalvari-ances.Thisexpansioncanbeundoneusingshrink=TRUEindis-playcommands.ThehelppageformetaMDSwillgivemoredetails,andpointtoexplanationoffunctionsusedinthefunction.

2.2Communitydissimilarities

Non-metricmultidimensionalscalingisagoodordinationmethodbe-causeitcanuseecologicallymeaningfulwaysofmeasuringcommunitydissimilarities.Agooddissimilaritymeasurehasagoodrankorderrela-tiontodistancealongenvironmentalgradients.Becausenmdsonlyusesrankinformationandmapsranksnon-linearlyontoordinationspace,itcanhandlenon-linearspeciesresponsesofanyshapeande ectivelyandrobustly ndtheunderlyinggradients.

ThemostnaturaldissimilaritymeasureisEuclideandistancewhichisinherentlyusedbyeigenvectormethodsofordination.Itisthedistanceinspeciesspace.Speciesspacemeansthateachspeciesisanaxisorthogonaltoallotherspecies,andsitesarepointsinthismultidimensionalhyper-space.However,Euclideandistanceisbasedonsquareddi erencesandstronglydominatedbysinglelargedi erences.Mostecologicallymean-ingfuldissimilaritiesareofManhattantype,andusedi erencesinsteadofsquareddi erences.Anotherfeatureingooddissimilarityindicesisthattheyareproportional:iftwocommunitiessharenospecies,theyhaveamaximumdissimilarity=1.EuclideanandManhattandissimilaritieswillvaryaccordingtototalabundanceseventhoughtherearenosharedspecies.

PackageveganhasfunctionvegdistwithBray–Curtis,JaccardandKulczy´nskiindices.AlltheseareoftheManhattantypeanduseonly rstorderterms(sumsanddi erences),andallarerelativizedbysiteto-talandreachtheirmaximumvalue(1)whentherearenosharedspeciesbetweentwocomparedcommunities.Functionvegdistisadrop-inre-placementforstandardRfunctiondist,andeitherofthesefunctionscanbeusedinanalysesofdissimilarities.

Therearemanyconfusingaspectsindissimilarityindices.Oneisthatsameindicescanbewrittenwithverydi erentlookingequations:twoalternativeformulationsofManhattandissimilaritiesinthemarginserve

2.2Communitydissimilarities

djk

N(xij xik)2

Euclidean

i=1djk= N|xij xik|Manhattan

i=1

A= Nxij

Nxik

i=1i=1

Nmin(xij,xik)

i=1

djk=A+B 2JManhattan

dA+B 2Jjk=

A+BBraydA+B 2Jjk=

A+B J

Jaccardd1 1 JJ

jk=2A+

Kulczy´nski

2.2Communitydissimilarities2ORDINATION:BASICMETHOD

forA=BforA=B

dAB=0dAB>0dAB=dBAdAB≤dAx+dxB

asanexample.Anothercomplicationisnaming.Functionvegdistusescolloquialnameswhichmaynotbestrictlycorrect.ThedefaultindexinveganiscalledBray(orBray–Curtis),butitprobablyshouldbecalledSteinhausindex.Ontheotherhand,itscorrectnamewassupposedtobeCzekanowskiindexsomeyearsago(butnowthisisregardedasanotherindex),anditisalsoknownasSørensenindex(butusuallymisspelt).Strictlyspeaking,Jaccardindexisbinary,andthequantitativevariantinveganshouldbecalledRuˇziˇckaindex.However,vegan ndseitherquantitativeorbinaryvariantofanyindexunderthesamename.

Thesethreebasicindicesareregardedasgoodindetectinggradi-ents.Inaddition,vegdistfunctionhasindicesthatshouldsatisfyothercriteria.Morisita,Horn–Morisita,Raup–Cric,BinomialandMountfordindicesshouldbeabletocomparesamplingunitsofdi erentsizes.Eu-clidean,CanberraandGowerindicesshouldhavebettertheoreticalprop-erties.

FunctionmetaMDSusedBray-Curtisdissimilarityasdefault,whichusuallyisagoodchoice.Jaccard(Ruˇziˇcka)indexhasidenticalrankorder,buthasbettermetricproperties,andprobablyshouldbepreferred.Functionrankindexinvegancanbeusedtostudywhichoftheindicesbestseparatescommunitiesalongknowngradientsusingrankcorrelationasdefault.Thefollowingexampleusesallenvironmentalvariablesindatasetvarechem,butstandardizesthesetounitvariance:

>data(varechem)

>rankindex(scale(varechem),varespec,c("euc","man","bray","jac","kul"))eucmanbrayjackul0.23960.27350.28380.28380.2840

arenon-linearlyrelated,buttheyhaveidenticalrankorders,andtheirrankcorrelationsareidentical.Ingeneral,thethreerecommendedindicesarefairlyequal.

Itookaverypracticalapproachonindicesemphasizingtheirabilitytorecoverunderlyingenvironmentalgradients.Manytextbooksempha-sizemetricpropertiesofindices.Theseareimportantinsomemethods,butnotinnmdswhichonlyusesrankorderinformation.Themetricpropertiessimplysaythat

1.iftwositesareidentical,theirdistanceiszero,

2.iftwositesaredi erent,theirdistanceislargerthanzero,3.distancesaresymmetric,and

4.theshortestdistancebetweentwositesisaline,andyoucannotimprovebygoingthroughothersites.Theseallsoundverynaturalconditions,buttheyarenotful lledbyalldissimilarities.Actually,onlyEuclideandistances–andprobablyJaccardindex–ful llallconditionsamongthedissimilaritiesdiscussedhere,andaremetrics.Manyotherdissimilaritiesful llthree rstconditionsandaresemimetrics.

Thereisaschoolthatsaysthatweshouldusemetricindices,andmostnaturally,Euclideandistances.Oneoftheirdrawbackswasthat

2ORDINATION:BASICMETHODtheyhaveno xedlimit,buttwositeswithnosharedspeciescanvaryindissimilarities,andevenlookmoresimilarthantwositessharingsomespecies.Thiscanbecuredbystandardizingdata.SinceEuclideandis-tancesarebasedonsquareddi erences,anaturaltransformationistostandardizesitestoequalsumofsquares,ortotheirvectornormusingfunctiondecostand:

>dis<-vegdist(decostand(varespec,"norm"),"euclid")

Thisgiveschorddistanceswhichreachamaximumlimitof√

whentherearenosharedspeciesbetweentwosites.AnotherrecommendedalternativeisHellingerdistancewhichisbasedonsquarerootsofsitesstandardizedtounittotal:

>dis<-vegdist(decostand(varespec,"hell"),"euclidean")

Despitestandardization,thesestillareEuclideandistanceswithalltheirgoodproperties,butfortransformeddata.Actually,itisoftenusefultotransformorstandardizedataevenwithotherindices.Ifthereisalargedi erencebetweensmallestnon-zeroabundanceandlargestabundance,wewanttoreducethisdi uallysquareroottransformationissu cienttobalancethedata.Wisconsindoublestandardizationoftenimprovesthegradientdetectionabilityofdissimilarityindices;thiscanbeperformedusingcommandwisconsininvegan.Herewe rstdivideallspeciesbytheirmaxima,andthenstandardizesitestounittotals.Afterthisstandardization,manydissimilarityindicesbecomeidenticalinrankorderingandshouldgiveequalresultsinnmds.

Youarenotrestrictedtouseonlyvegdistindicesinvegan:vegdistreturnssimilardissimilaritystructureasstandardRfunctiondistwhichalsocanbeused,aswellasanyothercompatiblefunctioninanypackage.Somecompatiblefunctionsaredsvdis(labdsvpackage),daisy(clusterpackage),anddistance(analoguepackage),andbetadiversityindicesinbetadiverinvegan.Morever,veganhasfunctiondesigndistwhereyoucande neyourowndissimilarityindicesbywritingitsequationusingeitherthenotationforA,BandJabove,orwithbinarydata,the2×2contingencytablenotationwhereaisthenumberofspeciesfoundonbothcomparedsites,andbandcarenumbersofspeciesfoundonlyinoneofthesites.Thefollowingthreeequationsde nethesameSørensenindexwherethenumberofsharedspeciesisdividedbytheaveragespeciesrichnessofcomparedsites:

>d<-vegdist(varespec,"bray",binary=TRUE)>d<-designdist(varespec,"(A+B-2*J)/(A+B)")

>d<-designdist(varespec,"(b+c)/(2*a+b+c)",abcd=TRUE)

Functionbetadiverde nessomemorebinarydissimilarityindicesinvegan.

Mostpublisheddissimilarityindicescanbeexpressedasdesigndistformulae.However,itismucheasierandsafertousethecannedalter-nativesinexistingfunctions:itisveryeasytomakeerrorsinwritingthedissimilarityequations.

2.2Communitydissimilarities

Quadraticterms

J=N

=1xijxA= iNikB= i2N=1xiji=1x2ik

Minimumterms

J=N

=1min(xij,xik)A= iN=1xB= iNiji=1xikBinarytermsJ=SharedspeciesA=No.ofspeciesinjB=

No.ofspeciesink

Sitek

presentabsent

Sitej

presentababsent

J=aA=a+bB=a+c

2.3Comparingordinations:Procrustesrotation2ORDINATION:BASICMETHOD

2.3Comparingordinations:Procrustesrotation

Procrustes errors

Twoordinationscanbeverysimilar,butthismaybedi culttosee,becauseaxeshaveslightlydi erentorientationandscaling.Actually,innmdsthesign,orientation,scaleandlocationoftheaxesarenotde- ned,althoughmetaMDSusessimplemethodto xthelastthreecompo-nents.ThebestwaytocompareordinationsistouseProcrustesrotation.Procrustesrotationusesuniformscaling(expansionorcontraction)androtationtominimizethesquareddi erencesbetweentwoordinations.PackageveganhasfunctionprocrustestoperformProcrustesanalysis.

HowmuchdidwegainwithusingmetaMDSinsteadofdefaultisoMDS?

0.4

Dimension 2

0.0

>>>>>tmp<-wisconsin(sqrt(varespec))dis<-vegdist(tmp)

vare.mds0<-isoMDS(dis,trace=0)pro<-procrustes(vare.mds,vare.mds0)pro

0.2

Call:

procrustes(X=vare.mds,Y=vare.mds0)Procrustessumofsquares:0.156>plot(pro)

0.4 0.2

0.4 0.20.00.20.40.6

Dimension 1

Procrustes errors

0.30

Inthiscasethedi erenceswerefairlysmall,andmainlyconcernedtwopoints.Youcanuseidentifyfunctiontoidentifythosepointsinaninteractivesession,oryoucanaskaplotofresidualdi erencesonly:

>plot(pro,kind=

0.250.000.05

Thedescriptivestatisticis“Procrustessumofsquares”orthesumofsquaredarrowsintheProcrustesplot.Procrustesrotationisnonsym-metric,andthestatisticwouldchangewithreversingtheorderofordina-tionsinthecall.Withargumentsymmetric=TRUE,bothsolutionsare rstscaledtounitvariance,andamorescale-independentandsymmetricstatisticisfound(oftenknownasProcrustesm2).

Procrustes residual

0.100.150.20

2.4

Index

Eigenvectormethods

methodnmdsmdspcaca

metricanyanyEuclideanChi-square

mappingnonlinearlinearlinear

weightedlinear

djk

= (xij xik)2

i=1

Non-metricmultidimensionalscalingwasahardtask,becauseanykindofdissimilaritymeasurecouldbeusedanddissimilaritieswerenonlinearlymappedintoordination.Ifweacceptonlycertaintypesofdissimilaritiesandmakealinearmapping,theordinationbecomesasimpletaskofrotationandprojection.Inthatcasewecanuseeigenvectormethods.Principalcomponentsanalysis(pca)andcorrespondenceanalysis(ca)arethemostimportanteigenvectormethodsincommunityordination.Inaddition,principalcoordinatesanalysisa.k.a.metricscaling(mds)isusedoccasionally.PcaisbasedonEuclideandistances,caisbasedonChi-squaredistances,andprincipalcoordinatescanuseanydissimilarities(butwithEuclideandistancesitisequaltopca).

Pcaisastandardstatisticalmethod,andcanbeperformedwithbaseRfunctionsprcomporprincomp.Correspondenceanalysisisnotasubiq-uitous,butthereareseveralalternativeimplementationsforthatalso.In

2ORDINATION:BASICMETHODthistutorialIshowhowtoruntheseanalyseswithveganfunctionsrdaandccawhichactuallyweredesignedforconstrainedanalysis.

Principalcomponentsanalysiscanberunas:

>vare.pca<-rda(varespec)>vare.pca

Call:rda(X=varespec)InertiaRank

Total1826Unconstrained182623Inertiaisvariance

Eigenvaluesforunconstrainedaxes:PC1PC2PC3PC4PC5PC6PC7PC8983.0464.3132.373.948.437.025.719.7

(Showedonly8ofall23unconstrainedeigenvalues)>plot(vare.pca)

Theoutputtellsthatthetotalinertiais1826,andtheinertiaisvari-ance.Thesumofall23(rank)eigenvalueswouldbeequaltothetotalinertia.Inotherwords,thesolutiondecomposesthetotalvarianceintolinearcomponents.Wecaneasilyseethatthevarianceequalsinertia:

>sum(apply(varespec,2,var))[1]1826

Functionapplyappliesfunctionvarorvariancetodimension2orcolumns(species),andthensumtakesthesumofthesevalues.Inertiaisthesumofallspeciesvariances.Theeigenvaluessumuptototalinertia.Inotherwords,theyeach“explain”acertainproportionoftotalvariance.The rstaxis“explains”983/1826=53.8%oftotalvariance.

Thestandardordinationplotcommandusespointsorlabelsforspeciesandsites.Somepeopleprefertousebiplotarrowsforspeciesinpcaandpossiblyalsoforsites.Thereisaspecialbiplotfunctionforthispurpose:

>biplot(vare.pca,scaling=-1)

Forthisgraphwespeci edscaling=-1.Theresultsarescaledonlywhentheyareaccessed,andwecan exiblychangethescalinginplot,biplotandothercommands.Thenegativevaluesmeanthatspeciesscoresaredividedbythespeciesstandarddeviationssothatabundantandscarcespecieswillbeapproximatelyasfarawayfromtheorigin.

Thespeciesordinationlookssomewhatunsatisfactory:onlyreindeerlichens(Cladina)andPleuroziumschreberiarevisible,andallotherspeciesarecrowdedattheorigin.Thishappensbecauseinertiawasvari-ance,andonlyabundantspecieswithhighvariancesareworthexplaining(butwecouldhidethisinplotbysettingnegativescaling).Standard-izingallspeciestounitvariance,orusingcorrelationcoe cientsinsteadofcovarianceswillgiveamorebalancedordination:

>vare.pca<-rda(varespec,scale=TRUE)>vare.pca

2.4Eigenvectormethods

Cla.ran

Cla.arb

132

CP0

Cet.nivPol.pilCla.spBar.lycPin.sylCet.eriIch.eriDic.polCla.fimCla.criCla.cerCla.botDic.fusDic.spLed.palCla.corDes.fleCla.chlPoh.nutPol.junCet.islPti.cilEmp.nigVac.myrHyl.spl21Vac.vit2

152212

10Cla.ste

4 Ple.sch27

628

4 20246810

PC1

Cla.arbCla.ran

Ich.eri

Ste.sp62

Vac.uli

Cla.amaPol.pilCla.cocCla.gra18Cal.vul

13Dip.mon2

CCet.eriCet.nivP14Cla.fim430Cla.def20Cla.cer232121522Dic.fus25Pol.jun24

Cla.spPin.syl12Dic.spCet.islCla.chlCla.phy

Led.palCla.ste10

Emp.nig9

27Des.fleVac.vitPoh.nut

Vac.myr

Hyl.spl

Ple.sch

4 20

246

PC1

2.4Eigenvectormethods9

11Cla.phy

10Cla.coc

Poh.nut

Cla.stePin.sylCla.chl

Cla.spCet.isl

Cla.gra1223Vac.vit5

Cla.ran14Pol.pilCet.eriCla.cri

Cla.fimCla.arbSte.sp618Dic.pol21

Ich.eriCla.ama13Cla.uncDip.monCla.corCla.defPel.aph24Bet.pubCla.bot

Bar.lycPti.cilC4

7Cal.vulCet.niv1516Vac.uli20Dic.sp

Emp.nig

Led.palVac.myr

Des.fle

Ple.sch

Hyl.spl28

10123

PC1

.19Bar.lyc2128

Bet.pub

0.Pti.cil

Vac.myrHyl.spl

27Cla.ste

Cla.phyCla.chl

Cet.islDes.fle

Cla.bot

.0Cla.sp

Pin.sylPoh.nutPle.sch3

Emp.nigVac.vitNep.arc

Dic.sp24250

.Pol.jun

0Cla.cer

Cla.fim2

Cla.graCla.corPel.aph23A1522C5

Cla.coc

Cla.criCla.def20Cet.eri

Dic.fus

.0Cet.niv 4

Dip.monCla.ran

Cla.unc16

Pol.pil0

.Cla.amaCla.arb181Vac.uli

Cal.vul 14

Ste.sp

13.1Ich.eri

.2 5

1012

CA1

2ORDINATION:BASICMETHOD

Call:rda(X=varespec,scale=TRUE)InertiaRank

Total44Unconstrained4423Inertiaiscorrelations

Eigenvaluesforunconstrainedaxes:

PC1PC2PC3PC4PC5PC6PC7PC88.904.764.263.732.962.882.732.18

(Showedonly8ofall23unconstrainedeigenvalues)>plot(vare.pca,scaling=3)

Nowinertiaiscorrelation,andthecorrelationofavariablewithitselfisone.Thusthetotalinertiaisequaltothenumberofvariables(species).Therankorthetotalnumberofeigenvectorsisthesameaspreviously.Themaximumpossiblerankisde nedbythedimensionsofthedata:itisonelessthansmallerofnumberofspeciesornumberofsites:

>dim(varespec)[1]2444

Iftherearespeciesorsitessimilartoeachother,rankwillbereducedevenfromthis.

Thepercentageexplainedbythe rstaxisdecreasedfromthepreviouspcaabundant.Thisspeciesisnatural,withsincehighpreviouslyvariances,webutneedednowwetohave“explain”toexplainonlytheallspeciesequally.Weshouldnotlookblindlyatpercentages,buttheresultweget.

Correspondenceanalysisisverysimilartopca:

>vare.ca<-cca(varespec)>vare.ca

Call:cca(X=varespec)

InertiaRank

Total2.08Unconstrained2.0823

Inertiaismeansquaredcontingencycoefficient

Eigenvaluesforunconstrainedaxes:

CA1CA2CA3CA4CA5CA6CA7CA80.52490.35680.23440.19550.17760.12160.11550.0889(Showedonly8ofall23unconstrainedeigenvalues)>plot(vare.ca)

Nowtheinertiaiscalledmeansquaredcontingencycoe cient.Corre-spondenceanalysisisbasedonChi-squareddistance,andtheinertiais

theChi-squaredstatisticofadatamatrixstandardizedtounittotal:

>chisq.test(varespec/sum(varespec))

Pearson'sChi-squaredtest

data:varespec/sum(varespec)

X-squared=2.083,df=989,p-value=1

2ORDINATION:BASICMETHOD2.5Detrendedcorrespondenceanalysis

YoushouldnotpayanyattentiontoP-valueswhicharecertainlymis-leading,butnoticethatthereportedX-squaredisequaltotheinertiaabove.

Correspondenceanalysisisaweightedaveragingmethod.Inthegraphabovespeciesscoreswereweightedaveragesofsitescores.Withdi erentscalingofresults,wecoulddisplaythesitescoresasweightedaveragesofspeciesscores:

>plot(vare.ca,scaling=1)

Wealreadysawanexampleofscaling=3orsymmetricscalinginpca.Theothertwointegersmeanthateitherspeciesareweightedaveragesofsites(2)orsitesareweightedaveragesofspecies(1).Whenwetakeweightedaverages,therangeofaveragesshrinksfromtheoriginalval-ues.Theshrinkagefactorisequaltotheeigenvalueofca,whichhasatheoreticalmaximumof1.

2.5Detrendedcorrespondenceanalysis

Correspondenceanalysisisamuchbetterandmorerobustmethodfor

communityordinationthanprincipalcomponentsanalysis.However,withlongecologicalgradientsitsu ersfromsomedrawbacksor“faults”whichwerecorrectedindetrendedcorrespondenceanalysis(dca):

Singlelonggradientsappearascurvesorarcsinordination(arce ect):thesolutionistodetrendthelateraxesbymakingtheirmeansequalalongsegmentsofpreviousaxes.

Sitesarepackedmorecloselyatgradientextremesthanatthecen-tre:thesolutionistorescaletheaxestoequalvariancesofspeciesscores.

Rarespeciesseemtohaveanundulyhighin uenceontheresults:thesolutionisstodownweightrarespecies.

AllthesethreeseparatetricksareincorporatedinfunctiondecoranawhichisafaithfulportofMarkHill’soriginalprogrammewiththesamename.Theusageissimple:

>vare.dca<-decorana(varespec)>vare.dca

Call:

decorana(veg=varespec)

Detrendedcorrespondenceanalysiswith26segments.Rescalingofaxeswith4iterations.DCA1DCA2DCA3DCA4

Eigenvalues0.5240.3250.20010.1918Decoranavalues0.5250.1570.09670.0608Axislengths2.8162.2051.54651.6486>plot(vare.dca,display="sites")

Bar.lyc2Bet.pub

Hyl.spl

Pti.cil

Cla.ste

Cla.phyCla.chl

Des.fle

Cla.bot10

921

2Cla.sp

12Pin.sylPoh.nutEmp.nig27Ple.schVac.vit

Dic.spCla.cer3Pol.jun

2425Nep.arc2

ACCla.fimPel.aphCla.cor2015

224

Cla.gra16Cla.coc

18Cla.criCla.def1Cet.niv

613Cet.eriDic.fus14

Dip.monCla.ran7

Cla.uncPol.pil

5Cla.amaCla.arb

Vac.uli

Cal.vul2

Ste.sp

Ich.eri

2 10123

CA1

2.6Ordinationgraphics21

.1109

212

.28

027

CD0.020

2418

2515522

.0 4

616

713

0.5

114

1.0 0.50.00.51.01.5

DCA1

2ORDINATION:BASICMETHOD

Functiondecorana ndsonlyfouraxes.Eigenvaluesarede nedasshrinkagevaluesinweightedaverages,similarlyasinccaabove.The“Decoranavalues”arethenumbersthattheoriginalprogrammereturnsas“eigenvalues”—Ihavenoideaoftheirpossiblemeaning,andtheyshouldnotbeused.Mostoftenpeoplecommentonaxislengths,whichsometimesarecalled“gradientlengths”.Theetymologyisobscure:thesearenotgradients,butordinationaxes.Itisoftensaidthatiftheaxislengthisshorterthantwounits,thedataarelinear,andpcashouldbeused.Thisisonlyfolkloreandnotbasedonresearchwhichshowsthatcaisatleastasgoodaspcawithshortgradients,andusuallybetter.

Thecurrentdatasetishomogeneous,andthee ectsofdcaarenotverylarge.Inheterogeneousdatawithacleararce ectthechangesoftenaremoredramatic.Rescalingmayhavelargerin uencethandetrendinginmanycases.

Thedefaultanalysisiswithoutdownweightingofrarespecies:seehelppagesfortheneededarguments.Actually,downweightisanindependentfunctionthatcanbeusedwithccaaswell.

Thereisaschoolofthoughtthatregardsdcaasthemethodofchoiceinunconstrainedordination.However,itseemstobeafragileandvaguebackoftricksthatisbetteravoided.

2.6Ordinationgraphics

Wehavealreadyseenmanyordinationdiagramsinthistutorialwithonefeatureincommon:theyareclutteredandlabelsaredi culttoread.Ordinationdiagramsaredi culttodrawcleanlybecausewemustputalargenumberoflabelsinasmallplot,andoftenitisimpossibletodrawcleanplotswithallitemslabelled.Inthischapterwelookatproducingcleanerplots.Forthiswemustlookattheanatomyofplottingfunctionsinveganandseehowtogainabettercontrolofdefaultfunctions.

Ordinationfunctionsinveganhavetheirdedicatedplotfunctionswhichprovidesasimpleplot.Forinstance,theresultofdecoranaisdisplayedbyfunctionplot.decoranawhichbehindthescenesiscalledbyourplotfunction.Alternatively,wecanusefunctionordiplotwhichalsoworkswithmanynon-veganordinationfunctions,butusespointsinsteadoftextasdefault.Theplot.decoranafunction(orordiplot)actuallyworksinthreestages:

1.Itdrawsanemptyplotwithlabelledaxes,butwithnosymbolsforsitesorspecies.2.Itusesfunctionstextorpointstoaddspeciestotheemptyframe.Iftheuserdoesnotaskspeci cally,thefunctionwillusetextinsmalldatasetsandpointsinlargedatasets.3.Itaddsthesitessimilarly.

Forbettercontroloftheplotswemustrepeatthesestagesbyhand:drawanemptyplotandthenaddsitesand/orspeciesasdesired.

Inthischapterwestudyadi cultcase:plottingtheBarroColoradoIslandordinations.

2ORDINATION:BASICMETHOD2.6Ordinationgraphics

>data(BCI)

Thisisadi cultdatasetforplotting:ithas225speciesandthereisnowayoflabellingthemallcleanly–unlessweuseverylargeplottingareawithsmalltext.Wemustshowonlyaselectionofthespeciesorsmallpartsoftheplot.Firstanordinationwithdecoranaanditsdefaultplot:

>mod<-decorana(BCI)>plot(mod)

DCA2

Thereisanadditionalprobleminplottingspeciesordinationwiththesedata:

>names(BCI)[1:5]

[1]"Abarema.macradenium""Acacia.melanoceras"[3]"Acalypha.diversifolia""Acalypha.macrostachya"[5]"Adelia.triloba"

++++++++++

++++++++++++++++++++++++++++++++

+++++++++++++++++++++q+q++++q+q++++++++q+++q+++++q+++qqq++++++q+q++q++qqqq+++++++++qqq+q++++++q+qqqq++q+++++++qq+++++++++++qq++++++++++++++++++++++++++++++++++++

+++++++++++++++++++++++++

+++

++++

4 202

6 4 20DCA1

246

Thedatasetusesfullspeciesnames,andthereisnowayof ttingthoseinordinationgraphs.Thereisautilityfunctionmake.cepnamesinvegantoabbreviateLatinnames:

>shnam<-make.cepnames(names(BCI))>shnam[1:5]

[1]"Abarmacr""Acacmela""Acaldive""Acalmacr""Adeltril"

Alchlati

DCA2

>pl<-plot(mod,dis="sp")

Allveganordinationplotfunctionsreturninvisiblyanordiplotobjectwhichcontainsinformationonthepointsplotted.Thisinvisibleresultcanbecaughtandusedasinputtoidentify.Thefollowingselectivelylabelssomeextremespeciesasclicked:

>identify(pl,"sp",labels=shnam)

Theeasiestwaytoselectivelylabelspeciesistouseinteractiveiden-tifyfunction:whenyouclicknexttoapoint,itslabelwillappearonthesideyouclicked.Youcan nishlabellingclickingtherightmousebutton,orwithhandicappedone-buttonmouse,youcanhittheesckey.

AbarmacrSocrexor

++Pachquin

++++++++

+Entescho+Nectciss++++++++++++Pachsess+++++++++++++++++Ficuyopo+++++++++++++++++++++++++++Margnobi++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++Gustsupe+++++Poularma+++++++++++++++++++++++Ocotwhit++++++++++

+++

Macrrose+++++Tropcauc++++Cavaplat++++Sapibroa+++++++++

++Pourbico

++++

Senndari+

+Brosguia+

+Quasamar

6 4 20DCA1

246

Thereisan“ordinationtextorpoints”functionorditorpinvegan.Thisfunctionwilllabelanitemonlyifthiscanbedonewithoutover-writingpreviouslabels.Ifanitemcannotbelabelledwithtext,itwillbemarkedasapoint.Itemsareprocessedeitherfromthemargintowardthecentre,orindecreasingorderofpriority.Thefollowinggiveshigherprioritytothemostabundantspecies:

DCA2

CasecommAlchlati

>stems<-colSums(BCI)

>plot(mod,dis="sp",type="n")

>sel<-orditorp(mod,dis="sp",lab=shnam,priority=stems,pcol="gray",pch="+")

Wealsocanzoomintosomepartsoftheordinationdiagramsbysettingxlimandylim,andwecanseemoredetails.

Analternativetoorditorpisfunctionordilabelwhichdrawstextonopaquelabelsthatcoverotherlabelsbelowthem.Alllabelscannotbedisplayed,butatleasttheuppermostarereadable.Argumentpriorityworkssimilarlyasinorditorpandcanbeusedtoselectwhichofthelabelsaremostimportanttoshow:

++Socrexor+++Ingacocl+++++

+Entescho+++Cordbico+++++Pachsess++++++Ficuyopo++Cocccoro++++++++++++Unonpitt+++++Talinerv++++++++++++++++Tabeguay++++++++Heisacum+++++++++++++T+++Faraocci+++++++++++++++TFicuinsi+Sapiglan++++++++++++++++++Gustsupe+++++Garcinte+Poularma++++++++++++++++CaseaculSloatern

Ocotwhit++++++T+ripcumi

++Hirttria+MacrroseViromult++++Xylomacr+++Marilaxi++++Drypstan+++

4 202

6 4 20DCA1

246

3ENVIRONMENTALINTERPRETATION

>plot(mod,dis="sp",type="n")

ACD02 4 6 4 20246

DCA1

>ordilabel(mod,dis="sp",lab=shnam,priority=stems)

Finally,thereisfunctionordipointlabelwhichusesbothpointsand

labelstothesepoints.Thepointsarein xedpositions,butthelabelsareiterativelylocatedtominimizetheiroverlap.TheBarroColoradoIslanddatasethasmuchtoomanynamesfortheordipointlabelfunction,butitcanbeusefulinmanycases.

Inadditiontotheseautomaticfunctions,functionorditkplotallowseditingofplots.Ithaspointsin xedpositionswithlabelsthatcanbedraggedtobetterplaceswithamouse.Thefunctionusesdi erentgraphicaltoolset(Tcl/Tk)thanordinaryRgraphics,buttheresultscanbepassedtostandardRplotfunctionsforeditingordirectlysavedasgraphics les.Moreover,theordipointlabelouputcanbeeditedusingorditkplot.

Functionsidentify,orditorp,ordilabelandordipointlabelmayprovideaquickandeasywaytoinspectordinationresults.Oftenweneedabettercontrolofgraphics,andjudicuouslyselectthelabelledspecies.Inthatcasewecan rstdrawanemptyplot(withtype="n"),andthenuseselectargumentinordinationtextandpointsfunctions.Theselectargumentcanbeanumericvectorthatliststheindicesofselecteditems.Suchindicesaredisplayedfromidentifyfunctionswhichcanbeusedtohelpinselectingtheitems.Alternatively,selectcanbealogicalvectorwhichisTRUEtoselecteditems.Suchalistwasproducedinvisiblyfromorditorp.Youcannotseeinvisibleresultsdirectlyfromthemethod,butyoucancatchtheresultlikewedidaboveinthe rstorditorpcall,andusethisvectorasabasisforfullycontrolledgraphics.Inthiscasethe rstitemswere:

>sel[1:14]

AbarmacrAcacmelaAcaldiveAcalmacrAdeltrilAegipanaAlchcost

FALSEFALSEFALSEFALSEFALSEFALSEFALSEAlchlatiAlibedulAllopsilAlseblacAmaicoryAnacexceAndiiner

TRUEFALSEFALSEFALSETRUETRUEFALSE

3Environmentalinterpretation

Itisoftenpossibleto“explain”ordinationusingecologicalknowledgeonstudiedsites,uallyitispreferabletouseexternalenvironmentalvariablestointer-prettheordination.Therearemanywaysofoverlayingenvironmentalinformationontoordinationdiagrams.Oneofthesimplestistochangethesizeofplottingcharactersaccordingtoanenvironmentalvariables(argumentcexinplotfunctions).Theveganpackagehassomeusefulfunctionsfor ttingenvironmentalvariables.

3.1Vector tting

Themostcommonlyusedmethodofinterpretationisto tenvironmentalvectorsontoordination.The ttedvectorsarearrowswiththeinterpre-tation:

3ENVIRONMENTALINTERPRETATION Thearrowpointstothedirectionofmostrapidchangeinthetheenvironmentalvariable.Oftenthisiscalledthedirectionofthegradient.

Thelengthofthearrowisproportionaltothecorrelationbetweenordinationandenvironmentalvariable.Oftenthisiscalledthestrengthofthegradient.

Fittingenvironmentalvectorsiseasyusingfunctionenvfit.Theexampleusesthepreviousnmdsresultandenvironmentalvariablesinthedatasetvarechem:

>data(varechem)

>ef<-envfit(vare.mds,varechem,permu=999)>ef***VECTORS

NMDS1NMDS2r2Pr(>r)

N-0.0573-0.99840.250.036*P0.61970.78480.190.101K0.76650.64230.180.114Ca0.68520.72830.410.003**Mg0.63250.77450.430.004**S0.19140.98150.180.132Al-0.87160.49020.530.001***Fe-0.93600.35200.450.001***Mn0.7987-0.60170.520.001***Zn0.61760.78650.190.108Mo-0.90310.42940.060.512Baresoil0.9249-0.38030.250.055.Humdepth0.9328-0.36040.520.003**pH-0.64800.76170.230.060.---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1Pvaluesbasedon999permutations.

The rsttwocolumnsgivedirectioncosinesofthevectors,andr2givesthesquaredcorrelationcoe cient.Forplotting,theaxesshouldbescaledbythesquarerootofr2.Theplotfunctiondoesthisautomatically,andyoucanextractthescaledvalueswithscores(ef,"vectors").Thesigni cances(Pr>r),orP-valuesarebasedonrandompermutationsofthedata:ifyouoftengetasgoodorbetterR2withrandomlypermuteddata,yourvaluesareinsigni cant.

Youcanaddthe ttedvectorstoanordinationusingplotcommand.Youcanlimitplottingtomostsigni cantvariableswithargumentp.max.Asusual,moreoptionscanbefoundinthehelppages.

>plot(vare.mds,display="sites")>plot(ef,p.max=0.1)

3.2Surface tting

Vector ttingispopular,anditprovidesacompactwayofsimultaneouslydisplayingalargenumberofenvironmentalvariables.However,itimplies

3.2Surface tting

.0q

Mgq

AlpHq

.0Fe

qSq

M0Nq

Baresoil

.Humdepth

0 q

.0 q

0.4 0.20.00.20.40.6

NMDS1

3.3Factorsq

.0q

001

800

5506q65

.00 0

50 7 00

4005

0 q

750

D3.50q

M0 300 Nq

300q

250

.02 5q

0 200

150

.0 050

0.4 0.20.00.2

0.40.6

NMDS1

3ENVIRONMENTALINTERPRETATION

alinearrelationshipbetweenordinationandenvironment:directionandstrengthareallyouneedtoknow.Thismaynotalwaysbeappropriate.

Functionordisurf tssurfacesofenvironmentalvariablestoordi-nations.Itusesgeneralizedadditivemodelsinfunctiongamofpackagemgcv.Functiongamusesthinplatesplinesintwodimensions,andauto-maticallyselectsthedegreeofsmoothingbygeneralizedcross-validation.Iftheresponsereallyislinearandvectorsareappropriate,the ttedsur-faceisaplanewhosegradientisparalleltothearrow,andthe ttedcontoursareequallyspacedparallellinesperpendiculartothearrow.

InthefollowingexampleIintroducetwonewRfeatures:

Functionenvfitcanbecalledwithformulainterface.Formulahasaspecialcharactertilde(～),andtheleft-handsidegivestheordinationresults,andtheright-handsideliststheenvironmentalvariables.Inaddition,wemustde nethenameofthedatacon-tainingthe ttedvariables.

ThevariablesindataframesarenotvisibletoRsessionunlessthedataframeisattachedtothesession.Wemaynotwanttomakeallvariablesvisibletothesession,becausetheremaybesynonymousnames,andwemayusewrongvariableswiththesamenameinsomeanalyses.Wecanusefunctionwithwhichmakesthegivendataframevisibleonlytothefollowingcommand.

Nowwearereadyfortheexample.Wemakevector ttingforselectedvariablesandadd ttedsurfacesinthesameplot.

>ef<-envfit(vare.mds~Al+Ca,varechem)>plot(vare.mds,display="sites")>plot(ef)

>tmp<-with(varechem,ordisurf(vare.mds,Al,add=TRUE))

with(varechem,ordisurf(vare.mds,Ca,add=TRUE,col="green4"))

Functionordisurfreturnstheresultof ttedgam.Ifwesavethatresult,likewedidinthe rst twithAl,wecanuseitforfurtheranalyses,suchasstatisticaltestingandpredictionofnewvalues.Forinstance,fitted(ef)willgivetheactual ttedvaluesforsites.

3.3Factors

Classcentroidsareanaturalchoiceforfactorvariables,andR2canbeusedasagoodness-of- tstatistic.The“signi cance”canbetestedwithpermutationsjustlikeinvector tting.Variablescanbede nedasfactorsinR,andtheywillbetreatedaccordinglywithoutanyspecialtricks.

Asanexample,weshallinspectdunemeadowdatawhichhasseveralclassvariables.Functionenvfitalsoworkswithfactors:

>data(dune)

>data(dune.env)

>dune.ca<-cca(dune)

>ef<-envfit(dune.ca,dune.env,permutations=999)>

3ENVIRONMENTALINTERPRETATION***VECTORS

CA1CA2r2Pr(>r)

A10.99820.06060.310.043*---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1Pvaluesbasedon999permutations.***FACTORS:Centroids:CA1CA2Moisture1-0.75-0.14Moisture2-0.47-0.22Moisture40.18-0.73Moisture51.110.57ManagementBF-0.73-0.14ManagementHF-0.39-0.30ManagementNM0.651.44ManagementSF0.34-0.68UseHayfield-0.290.65UseHaypastu-0.07-0.56UsePasture0.520.05Manure00.651.44Manure1-0.46-0.17Manure2-0.59-0.36Manure30.52-0.32Manure4

-0.21

-0.88

Goodnessoffit:

r2Pr(>r)

Moisture0.410.005**Management0.440.001***Use0.180.079.Manure0.460.010**---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1Pvaluesbasedon999permutations.>plot(dune.ca,display="sites")>plot(ef)

Thenamesoffactorcentroidsareformedbycombiningthenameof

thefactorandthenameofthelevel.Nowtheaxesshowthecentroidsforthelevel,andtheR2valuesareforthewholefactor,justlikethesigni cancetest.Theplotlookscongested,andwemayusetricksof§2.6(p.12)tomakecleanerplots,butobviouslynotallfactorsarenecessaryininterpretation.

Packageveganhasseveralfunctionsforgraphicaldisplayoffactors.Functionordihulldrawsanenclosingconvexhullfortheitemsinaclass,ordispidercombinesitemstotheir(weighted)classcentroid,andordiellipsedrawsellipsesforclassstandarddeviations,standarder-rorsorcon denceareas.TheexampledisplaysalltheseforManagementtypeinthepreviousordinationandautomaticallylabelsthegroupsin

3.3Factors

ManagementNMManure0

AC1

UseHayfield18

Moisture5

1415A1

ManagementBFUsePasture16

Moisture1ManagementHF5Manure2

7Moisture2Manure1Manure38UseHaypastu

2ManagementSF1

Manure4Moisture449

123

13 1

2 1012

CA1

4CONSTRAINEDORDINATION

ordispidercommand:

AC1

2 101

CA1

>plot(dune.ca,display="sites",type="p")

>with(dune.env,ordiellipse(dune.ca,Management,kind="se",conf=0.95))>with(dune.env,ordispider(dune.ca,Management,col="blue",label=TRUE))>

with(dune.env,ordihull(dune.ca,Management,col="blue",lty=2))

Correspondenceanalysisisaweightedordinationmethod,andvegan

functionsenvfitandordisurfwilldoweighted tting,unlesstheuserspeci esequalweights.

Constrainedordination

Inunconstrainedordinationwe rst ndthemajorcompositionalvaria-tion,andthenrelatethisvariationtoobservedenvironmentalvariation.

Inconstrainedordinationwedonotwanttodisplayallorevenmostofthecompositionalvariation,butonlythevariationthatcanbeexplainedbytheusedenvironmentalvariables,orconstraints.Constrainedordina-tionisoftenknownas“canonical”ordination,butthisnameismisleading:thereisnothingparticularlycanonicalinthesemethods(seeyourfavoriteDictionaryfortheterm).Thenamewastakenintouse,becausethereisonespecialstatisticalmethod,canonicalcorrelations,buttheseindeedarecanonical:theyarecorrelationsbetweentwomatricesregardedtobesymmetricallydependentoneachother.Theconstrainedordinationisnon-symmetric:wehave“independent”variablesorconstraintsandwehave“dependent”variablesorthecommunity.Constrainedordinationratherisrelatedtomultivariatelinearmodels.

Theveganpackagehasthreeconstrainedordinationmethodswhichallareconstrainedversionsofbasicordinationmethods:

Constrainedanalysisofproximities(cap)infunctioncapscaleisrelatedtometricscaling(cmdscale).Itcanhandleanydissimilaritymeasuresandperformsalinearmapping.

Redundancyanalysis(rda)infunctionrdaisrelatedtoprincipalcomponentsanalysis.ItisbasedonEuclideandistancesandper-formslinearmapping.

Constrainedcorrespondenceanalysis(cca)infunctionccaisre-latedtocorrespondenceanalysis.ItisbasedonChi-squareddis-tancesandperformsweightedlinearmapping.

Wehavealreadyusedfunctionsrdaandccaforunconstrainedordination:theywillperformthebasicunconstrainedmethodasaspecialcaseifconstraintsarenotused.

Allthesethreeveganfunctionsareverysimilar.Thefollowingexam-plesmainlyusecca,butothermethodscanbeusedsimilarly.Actually,theresultsaresimilarlystructured,andtheyinheritpropertiesfromeachother.Forhistoricalreasons,ccaisthebasicmethod,andrdainheritspropertiesfromit.Functioncapscaleinheritsdirectlyfromrda,andthroughthisfromcca.Manyfunctions,arecommonwithallthesemeth-ods,andtherearespeci cfunctionsonlyifthemethoddeviatesfromits

4CONSTRAINEDORDINATION4.1Modelspeci cation

ancestor.Inveganversion2.0-6thefollowingclassfunctionsarede nedforthesemethods:

cca:add1,alias,anova,as.mlm,bstick,calibrate,coef,deviance,

drop1,eigenvals,extractAIC,fitted,goodness,model.frame,model.matrix,nobs,permutest,plot,points,predict,print,residuals,RsquareAdj,scores,screeplot,simulate,summary,text,tolerance,weights rda:as.mlm,biplot,coef,deviance,fitted,goodness,predict,RsquareAdj,scores,simulate,weights capscale:fitted,print,simulate.

Manyofthesemethodsareinternalfunctionsthatusersrarelyneed.

4.1Modelspeci cation

Therecommendedwayofde ningaconstrainedmodelistousemodelformula.Formulahasaspecialcharacter～,andonitsleft-handsidegivesthenameofthecommunitydata,andright-handgivestheequationforconstraints.Inaddition,youshouldgivethenameofthedatasetwhereto ndtheconstraints.This tsaccaforvarespecconstrainedbysoilAl,KandP:

>a<-cca(varespec~Al+P+K,varechem)>a

Call:cca(formula=varespec~Al+P+K,data=varechem)

InertiaProportionRank

Total2.0831.000Constrained0.6440.3093Unconstrained1.4390.69120

InertiaismeansquaredcontingencycoefficientEigenvaluesforconstrainedaxes:CCA1CCA2CCA30.3620.1700.113

Eigenvaluesforunconstrainedaxes:

CA1CA2CA3CA4CA5CA6CA7CA80.35000.22010.18510.15510.13510.10030.07730.0537(Showedonly8ofall20unconstrainedeigenvalues)

Theoutputissimilarasinunconstrainedordination.Nowthetotalinertiaisdecomposedintoconstrainedandunconstrainedcomponents.Therewerethreeconstraints,andtherankofconstrainedcomponentisthree.Therankofunconstrainedcomponentis20,whenitusedtobe23inthepreviousanalysis.Therankisthesameasthenumberofaxes:youhave3constrainedaxesand20unconstrainedaxes.Insomecases,theranksmaybelowerthanthenumberofconstraints:someoftheconstraintsaredependentoneachother,andtheyarealiasedintheanalysis,andaninformativemessageisprintedwiththeresult.

Itisverycommontocalculatetheproportionofconstrainedinertiafromthetotalinertia.However,totalinertiadoesnothaveaclearmean-ingincca,andthemeaningofthisproportionisjustasobscure.Inrda

4.1Modelspeci cation13

Cal.vul

Bet.pub

Bar.lyc75

216

Led.pal

Dic.fusCla.botPti.cilIch.eri

182

C23CVac.myrCla.criCla.arb0Des.fle20Cla.defDic.polEmp.nigCla.fimSte.spVac.uliCla.amaCla.uncCet.islPin.sylCla.cocCla.ranVac.vit

Cla.graDip.mon0

Ple.schKPol.junCla.corPoh.nutPol.pilCet.eri

Cla.chlCla.sp1127

Nep.arc

Pel.aphCla.steAl

Hyl.splCla.phy

Dic.sp

12Cla.cer

210Cet.nivP

3 2 1

012

CCA1

4 3 2

AC2

C 4A0

C 3C1 2 1

0 1 23

CCA1

610

ManagementBF11

Viclat

ManagementHF2

TripraAchmil

Plalan19

Brohor

RumaceAntodo

Trirep

Hyprad

LolperPoapra

BelperLeoautAC0Brarut

JunartCPoatriElyrep

ManagementNM

EmpnigPotpalAirpraSalrepJunbufElepal

RanflaSagproCalcus

9Alogen

Agrsto1

3 4

ManagementSF

ChealbCirarv20

2 13

2 10

CCA1

4CONSTRAINEDORDINATION

thiswouldbetheproportionofvariance(orcorrelation).Thismayhaveaclearermeaning,buteveninthiscasemostofthetotalinertiamayberandomnoise.Itmaybebettertoconcentrateonresultsinsteadoftheseproportions.

Basicplottingworksjustlikeearlier:

>plot(a)

havesimilarinterpretationas ttedvectors:thearrowpointstothedirec-tionofthegradient,anditslengthindicatesthestrengthofthevariable

inthisdimensionalityofsolution.Thevectorswillbeofunitlengthinfullranksolution,buttheyareprojectedtotheplaneusedintheplot.Thereisalsoaprimitive3Dplottingfunctionordiplot3d(whichneedsuserinteractionfor nalgraphs)thatshowsallarrowsinfulllength:

>ordiplot3d(a,type="h")

Withfunctionordirglyoucanalsoinspect3Ddynamicplotsthatcanbespinnedorzoomedintowithyourmouse.

Theformulainterfaceworkswithfactorvariablesaswell:

>a<-cca(dune~Management,dune.env)>plot(a)>a

Call:cca(formula=dune~Management,data=dune.env)InertiaProportionRank

Total2.1151.000Constrained0.6040.2853Unconstrained1.5110.71516

InertiaismeansquaredcontingencycoefficientEigenvaluesforconstrainedaxes:CCA1CCA2CCA30.3190.1820.103

Eigenvaluesforunconstrainedaxes:

CA1CA2CA3CA4CA5CA6CA70.447370.203000.163010.134570.129400.094940.07904

CA8CA9CA10CA11CA12CA13CA140.065260.050040.043210.038700.023850.017730.00917

CA15CA160.007960.00416

FactorvariableManagementhadfourlevels(BF,HF,NM,SF).Internally

Rexpressedthesefourlevelsasthreecontrasts(sometimescalled“dummyvariables”).Theappliedcontrastslooklikethis:

ManagementHFManagementNMManagementSF

BF000SF001HF100NM

010

Wedonotneedbutthreevariablestoexpressfourlevels:ifthereis

numberoneinacolumn,theobservationbelongstothatlevel,andif

4CONSTRAINEDORDINATIONthereisawholelineofzeros,theobservationmustbelongtotheomittedlevel,orthe rst.Thebasicplotfunctiondisplaysclasscentroidsinsteadofvectorsforfactors.

Inadditiontotheseordinaryfactors,Ralsoknowsorderedfactors.VariableMoistureindune.envisde nedasanorderedfour-levelfactor.Inthiscasethecontrastslookdi erent:

Moisture.LMoisture.QMoisture.C

1-0.67080.5-0.22362-0.2236-0.50.670840.2236-0.5-0.67085

0.67080.50.2236

Rusespolynomialcontrasts:thelineartermLisequaltotreatingMoistureasacontinuousvariable,andthequadraticQandcubicCtermsshownonlinearfeatures.Therewerefourdistinctlevels,andthenumberofcontrastsisoneless,justlikewithordinarycontrasts.Theordinationcon guration,eigenvaluesorrankdonotchangeifthefactorisunorderedorordered,butthepresentationofthefactorintheresultsmaychange:

>a<-cca(dune~Moisture,dune.env)>plot(a)

Nowplotshowsboththecentroidsoffactorlevelsandthecontrasts.

Ifwecouldchangetheorderedfactortoacontinuousvector,onlythelineare ectarrowwouldbeimportant.Iftheresponsetothevariableisnonlinear,thequadratic(andcubic)arrowswouldbelongaswell.

Ihaveexplainedonlythesimplestusageoftheformulainterface.Theformulaisverypowerfulinmodelspeci cation:youcantransformyourcontrastswithintheformula,youcande neinteractions,youcanusepolynomialcontrastsetc.However,modelswithinteractionsorpolyno-mialsmaybedi culttointerpret.

4.2Permutationtests

Thesigni canceofalltermstogethercanbeassessedusingpermutationtests:thedatarepermutedrandomlyandthemodelisre tted.Whenconstrainedinertiainpermutationsisnearlyalwayslowerthanobservedconstrainedinertia,wesaythatconstraintsaresigni cant.

Theeasiestwayofrunningpermutationtestsistousethemockanovafunctioninvegan:

>anova(a)

Permutationtestforccaunderreducedmodel

Model:cca(formula=dune~Moisture,data=dune.env)

DfChisqFN.PermPr(>F)

Model30.632.251990.005**Residual161.49---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1

TheModelreferstotheconstrainedcomponent,andResidualtotheunconstrainedcomponentoftheordination,Chisqisthecorresponding

4.2Permutationtests

Moisture4

213

Junbuf132

ACSagpro