##############################################################

###UNCONSTRAINED ORDINATION ANALYSIS--------------------------

library(vegan)

##PCA with enviro data

PCA <- rda(env.tab.s, scale = F)

PCA


#Total inertia

sum (apply (env.tab.s, 2, var))


#Variance explained by PC1

4.597/15


#plot PCA results

biplot(PCA, display = 'species', choices=c(1,2))

biplot(PCA, display = 'species', choices=c(2,3))


##selection of PCA axes

#above-average eigenvalue

mean(PCA$CA$eig)

barplot(PCA$CA$eig)

lines(c(0,15), c(1,1), col="red")


#PCA scores

PCA.sc <- scores(PCA, choices = c(1,2), display = "sites")

head(PCA.sc)


#map first PCA axis

library(berryFunctions)

colPoints(coord$X, coord$Y, PCA.sc[,1], add=F)

colPoints(coord$X, coord$Y, PCA.sc[,2], add=F)


###PCoA----------------------------------------------------------

#claculate PCoA

pcoa <- cmdscale(beta.bc, k=3)

head(pcoa)


#plot ordination diagram

plot(pcoa[,1], pcoa[,2])


#add partition from clustering analysis

cl.ward #<- cutree(clust.ward, 4)

plot(coord, pch=21, bg=my.cols[cl.ward])


plot(pcoa[,1:2], pch=21, bg=my.cols[cl.ward])


#plot 2nd and 3rd axis

plot(pcoa[,2:3], pch=21, bg=my.cols[cl.ward])


#fit enviro. variables

ef <- envfit(pcoa[,1:2], env.tab.s)

ef

plot(pcoa[,1:2], pch=21, bg=my.cols[cl.ward])

plot(ef)


#fit to 2nd and 3rd axis

ef2 <- envfit(pcoa[,2:3], env.tab.s)

ef2

plot(pcoa[,2:3], pch=21, bg=my.cols[cl.ward])

plot(ef2)


#plot surface

plot(pcoa[,1:2], pch=21, col="black", bg=my.cols[cl.ward], cex=1.3)

ordisurf(pcoa[, 1:2], forest$Diversity, col="black", add=T)


###NMDS-----------------------------------------------------------

##2D

nmds <- metaMDS(beta.b, k=3, trymax = 30, previous.best=F)

nmds

stressplot(nmds)

ordiplot(nmds)


#export scores

nmds.sc <- as.data.frame(scores(nmds))


#add partition from clustering object

plot(nmds.sc[,1:2], pch=21, col="black", bg=my.cols[cl.ward], cex=1.3)


#plot 3D

library("vegan3d")

library("rgl")


ordirgl (nmds, col=my.cols[cl.ward])

rgl.bg(color = "white")