setwd() veg <- read.csv("Abundance data_vegetation.csv", header=T, sep=",") veg[1:28,1:5] com <- veg[,-c(1:3)] com[1:5,1:5] ## nmds using presence absence data com1 <- sign(com) ### change into presence/absence data library(vegan) ord_veg <- metaMDS(com1,'jaccard',autotransform=F, k=2) ord1 <- metaMDS(com1,'jaccard',autotransform=F, k=3) ord_veg <- metaMDS(com1,'jaccard',autotransform=F,k=4) ord_veg <- metaMDS(com1,'jaccard',autotransform=F,k=5) tiff(filename="NMDS_vegetation.tiff", res = 600, width=3000, height=3000, compression = "lzw") par(mfrow=c(1,1)) par(mar=c(6,6,2,2)) par(pty='s') groupz <- sort(unique(veg$HAB)) ordiplot(ord1, disp="sites", type="n",cex.axis=1.5,cex.lab=2,xlim = c(-1,0.5),ylim=c(-0.5,0.5)) #ordiplot(ord1, disp="sites", type="text",cex.axis=1.5,cex.lab=2,xlim = c(-1,0.5),ylim=c(-0.5,0.5)) for(i in seq(groupz)) { ordiellipse(ord1, veg$HAB, kind="se", conf=0.95, draw="polygon", col=i, show.groups=groupz[i]) } points(ord1, 'sites', col = as.numeric(veg$HAB),pch=16,cex=1.5) ordisurf(ord1,veg$BA, add=TRUE,col='black') dev.off() ##data merging with decomposition data data3 <- veg[c(2)] nmds_val <- ord1$points nmds_val <- data.frame(nmds_val) names(nmds_val) <- c("veg_MDS1","veg_MDS2","veg_MDS3") data4 <- data.frame(data3, nmds_val) ##the data4 were merged with "Data_landscape experiment" ########## ##Soil data has not been made available. #In case of necessity, please contact Marleen de Blécourt at "" load('soil.Rdata') head(soil) soil <- soil[order(soil$FOR_TYPE),] names(soil)[16]<-paste("pH") names(soil)[21]<-paste("Elevation") names(soil)[20]<-paste("Slope") sl_pca <- soil[,6:21] names(sl_pca) plot(sl_pca) sl_pca <- sl_pca[,c(1:3,6:12,15:16)] sd<-as.data.frame(scale(sl_pca)) head(sd) pca2 <- prcomp(sd) summary(pca2) biplot(pca2,xlabs = rep (".", dim(pca2$x)[1] )) # find the scores on axes required site.scr2 <- pca2$x# x spp.scr2 <- pca2$rotation [,c(1,2)]# y ## generate x,y lims xlims2 <- range(site.scr2[,1], spp.scr2[,1]) ylims2 <- range(site.scr2[,2], spp.scr2[,2]) ## plotting biplot with forest types with different colour points tiff(filename = "soil.pca.for_type11.tiff", res = 600, width=4000, height=4000, compression = 'lzw') biplot(pca2, choices=1:2, cex = 1.2, col="black", scale=1,xlabs = rep ("", dim(pca2$x)[1] )) points((pca2$x[1:50,1])*1.5, (pca2$x[1:50,2])*1.5, col = "black", pch = 20,cex=1.5) points((pca2$x[51:109,1])*1.5, (pca2$x[51:109,2])*1.5, col = "red", pch = 20,cex=1.5) points((pca2$x[110:139,1])*1.5, (pca2$x[110:139,2])*1.5, col = "green", pch = 20,cex=1.5) dev.off() ####the PC1 and PC2 values at subplot levels were merged #with decomposition data to get "Data_environmental factors.csv"