#####################################
#Load libraries
#####################################
library(vegan)
library(RColorBrewer)

env_data<-read.csv("environmental_data.csv", row.names='site', header=TRUE)

prin_comp<-rda(env_data[,2:14], scale=TRUE)
summary(prin_comp)
screeplot(prin_comp,bstick=TRUE)
ordiplot(prin_comp,type="text")
biplot(prin_comp,type = "points", axes=F, add=T)

# define x and y axis labels with percent variation explained for each axis
x_axis_text<-sprintf("Axis 1 (%.2f%%)",100*prin_comp$CA$eig[1]/prin_comp$CA$tot.chi)
y_axis_text<-sprintf("Axis 2 (%.2f%%)",100*prin_comp$CA$eig[2]/prin_comp$CA$tot.chi)
# Define a list of colors for symbols
point_colors<-c("red","blue","black","gold","tan","green","gray","purple")
point_colors<-brewer.pal(8,"Spectral")
# Get the row and column scores
pca_scores<-scores(prin_comp)
# Plot the points
plot(pca_scores$sites[,1],pca_scores$sites[,2],pch=21,bg=point_colors[env_data$drainage],xlab=x_axis_text,ylab=y_axis_text)
# Add arrows (vectors) for variables
arrows(0,0,pca_scores$species[,1],pca_scores$species[,2],lwd=1,length=0.1)
# Label the variables
text(pca_scores$species[,1],pca_scores$species[,2],rownames(pca_scores$species))
# Add a legend
legend("topleft",names(summary(env_data[,1])),text.col=point_colors)

# Calculate drainage centroids
drain_cent_axis1<-tapply(pca_scores$sites[,1],env_data$drainage,mean)
drain_cent_axis2<-tapply(pca_scores$sites[,2],env_data$drainage,mean)
# Plot drainage centroids as larger triangles
points(drain_cent_axis1,drain_cent_axis2,bg=point_colors,pch=24,cex=2)

source("sigpca2.R")
loading_sig<-sigpca2(env_data[,2:14],1000)

drain_sder_axis1<-tapply(pca_scores$sites[,1],env_data$drainage,sd)/(summary(env_data$drainage)^0.5)
drain_sder_axis2<-tapply(pca_scores$sites[,2],env_data$drainage,sd)/(summary(env_data$drainage)^0.5)
plot(pca_scores$sites[,1],pca_scores$sites[,2],type="n",xlab=x_axis_text,ylab=y_axis_text)
points(drain_cent_axis1,drain_cent_axis2,pch=21,cex=1.25,bg=point_colors)
arrows(drain_cent_axis1-drain_sder_axis1,drain_cent_axis2,drain_cent_axis1+drain_sder_axis1,drain_cent_axis2,length = 0.05,angle = 90,code=3)
arrows(drain_cent_axis1,drain_cent_axis2-drain_sder_axis2,drain_cent_axis1,drain_cent_axis2+drain_sder_axis2,length = 0.05,angle = 90,code=3)
# which loadings are significant 
filter<-loading_sig[,1]<0.001
arrows(0,0,pca_scores$species[filter,1],pca_scores$species[filter,2],col="red",length=0.1)
legend("topleft",names(summary(env_data[,1])),text.col=point_colors)
text(pca_scores$species[filter,1],pca_scores$species[filter,2],rownames(pca_scores$species)[filter])