######################### # Example: artificial fish data ######################### library(vegan) fishes <-read.table("DataAndR/artificialFishes.txt",h=T) fishes # response variables Resp <- fishes[,c("Sp1","Sp2","Sp3","Sp4","Sp5","Sp6")] # transform coral sand and other into factors fishes$Coral <- as.factor(fishes$Coral) fishes$Sand <- as.factor(fishes$Sand) fishes$Other <- as.factor(fishes$Other) # all substrates as one variable Substrate <- as.factor(c(rep("Sand",3),rep(c("Other","Coral"),3),"Other")) #rda myrda <- rda(Resp~Depth+Coral+Sand,data=fishes) pdf(file="Figures/fishesBip2.pdf") plot(myrda, scaling=2) dev.off() pdf(file="Figures/fishesBip1.pdf") plot(myrda, scaling=1) dev.off() # rda with Substrate myrda2 <- rda(Resp~Depth+Substrate,data=fishes) pdf(file="Figures/fishesSubsBip2.pdf") plot(myrda2,scaling=2) dev.off() ######################### # Example: height shoe weight ######################### shsw <-read.table("DataAndR/HeightShoeWeight.txt",header=TRUE) attach(shsw) head(shsw) hsw <- shsw[,2:4] head(hsw) fm.col <- character() fm.col[sex==0] <- "blue" fm.col[sex==1] <- "red" sqrt( length(sex)-1 ) # = 15 ######################### # Example: Height shoe weight RDA ######################### library(vegan) Expl <- hsw[,c(1,2,3)] Resp <- hsw[,c(1,2,3)] myrda <- rda(Resp~shoe,scale=TRUE,data=Expl) plot(myrda,scaling=1,type="n",main="Distance triplot, scaling=1") segments(x0=0,y0=0, x1=scores(myrda, display="species", scaling=1)[,1], y1=scores(myrda, display="species", scaling=1)[,2]) text(myrda, display="sp", scaling=1, col=2) text(myrda, display="bp", scaling=1, row.names(scores(myrda, display="bp")), col=3) points(myrda, display=c("lc"), scaling=1, pch=1,col=fm.col) plot(myrda,scaling=2,type="n",main="Correlation triplot, scaling=2") segments(x0=0,y0=0, x1=scores(myrda, display="species", scaling=2)[,1], y1=scores(myrda, display="species", scaling=2)[,2]) text(myrda, display="sp", scaling=2, col=2) text(myrda, display="bp", scaling=2, row.names(scores(myrda, display="bp")), col=3) points(myrda, display=c("lc"), scaling=2, pch=1,col=fm.col) # without shoe as Response variable Resp <- hsw[,c(1,3)] myrda <- rda(Resp~shoe,scale=TRUE,data=Expl) plot(myrda,scaling=1,type="n",main="Distance triplot, scaling=1") segments(x0=0,y0=0, x1=scores(myrda, display="species", scaling=1)[,1], y1=scores(myrda, display="species", scaling=1)[,2]) text(myrda, display="sp", scaling=1, col=2) text(myrda, display="bp", scaling=1, row.names(scores(myrda, display="bp")), col=3) points(myrda, display=c("lc"), scaling=1, pch=1,col=fm.col) plot(myrda,scaling=2,type="n",main="Correlation triplot, scaling=2") segments(x0=0,y0=0, x1=scores(myrda, display="species", scaling=2)[,1], y1=scores(myrda, display="species", scaling=2)[,2]) text(myrda, display="sp", scaling=2, col=2) text(myrda, display="bp", scaling=2, row.names(scores(myrda, display="bp")), col=3) points(myrda, display=c("lc"), scaling=2, pch=1,col=fm.col) ######################### # Example: RIKZ ######################### # Variables are explained on the slides and # in Chapter 27 of Zuur: Analysing Ecological Data library(vegan) RIKZ <- read.table("DataAndR/RIKZGroups.txt", header = TRUE) Species <- RIKZ[,2:5] #Data were square root transformed Species.sq <- sqrt(Species) I1 <- rowSums(Species) #Could be used to drop sites with a total #of 0. Just in case you want to apply one #of these db-RDA transformation ExplVar <- RIKZ[, c("angle1","exposure","salinity", "temperature","NAP","penetrability", "grainsize","humus","chalk", "sorting1")] RIKZ_RDA<-rda(Species.sq, ExplVar, scale=T) plot(RIKZ_RDA,scaling=2,main="Correlation biplot") # Correlation biplot, scaling=2 plot(RIKZ_RDA, scaling=2,main="Correlation",type="n") segments(x0=0,y0=0, x1=scores(RIKZ_RDA, display="species", scaling=2)[,1], y1=scores(RIKZ_RDA, display="species", scaling=2)[,2]) text(RIKZ_RDA, display="sp", scaling=2, col=2) text(RIKZ_RDA, display="bp", scaling=2, row.names(scores(RIKZ_RDA, display="bp")), col=3) text(RIKZ_RDA, display=c("sites"), scaling=2,labels=rownames(Species.sq)) cor(Species.sq,ExplVar) # Distance biplot, scaling=1 plot(RIKZ_RDA, scaling=1,main="Distance biplot",type="n") segments(x0=0,y0=0, x1=scores(RIKZ_RDA, display="species", scaling=1)[,1], y1=scores(RIKZ_RDA, display="species", scaling=1)[,2]) text(RIKZ_RDA, display="sp", scaling=1, col=2) text(RIKZ_RDA, display="bp", scaling=1, row.names(scores(RIKZ_RDA, display="bp")), col=3) text(RIKZ_RDA, display=c("sites"), scaling=1,labels=rownames(Species.sq)) # Proportion of variance explained RIKZ_RDA # 0.3751 = total variance explained by explanatory variables summary(RIKZ_RDA) # order of importance # need rda expressed with the formula option RIKZ_RDA2 <-rda(Species.sq~angle1+exposure+salinity, data=ExplVar, scale=T) anova.cca(RIKZ_RDA2) # whole model anova.cca(RIKZ_RDA2, by="margin") # marginal test for each expanatory variable