#####################################
#Load libraries
#####################################
library(BiodiversityR)
library(vegan)

#####################################
#Load Data
#####################################
# Load species data
community<-read.csv(file="black_creek_community.csv", row.names='samples', sep=",", header=TRUE)

# Load environmental data
envdata<-read.csv(file="black_creek_environmental.csv", row.names='samples', sep=",", header=TRUE)

# Drop rare species
community<-community[,apply(community>0,2,sum)>2]
community<-decostand(community, method="hellinger")
# Define a reduced and full models for AIC model selection
reduced_cca<-cca(log(community+1) ~ 1, envdata)
full_cca<-cca(log(community+1) ~ . ,envdata)
# model selection
model_fit <- ordiR2step(reduced_cca, scope=formula(full_cca))

# This final model has two variables, compare the constrained intertial to the full model with all variables
summary(model_fit)
vif.cca(model_fit)
# Now we have a simpler model that is explaining just as much variation

# Perform an overal permutive ANOVA of the CCA
anova(model_fit)
anova(model_fit, by = "terms", perm = 1000)
anova(model_fit, by = "axis", perm = 1000)

# define x and y axis labels with percent variation explained for each axis
x_axis_text<-sprintf("Axis 1 (%.2f%%)",100*model_fit$CCA$eig[1]/model_fit$CA$tot.chi)
y_axis_text<-sprintf("Axis 2 (%.2f%%)",100*model_fit$CCA$eig[2]/model_fit$CA$tot.chi)
# Get the row and column scores
cca_scores<-scores(model_fit)
constraints<-scores(model_fit, display="bp")

# Get the 15 most abundant species
abundance_cutoff<-sort(apply(community,2,sum),decreasing=T)[15]
plot_species<-apply(community,2,sum)>abundance_cutoff
species_matrix<-subset(cca_scores$species,plot_species)
sp_names<-make.cepnames(row.names(species_matrix))

# Triplot
plot(cca_scores$sites[,1],cca_scores$sites[,2],pch=16,col="black",cex=0.8,xlab=x_axis_text,ylab=y_axis_text)
# light dotted lines for x and y axes
abline(h=0,lty="dotted",col="gray")
abline(v=0,lty="dotted",col="gray")
text(species_matrix[,1],species_matrix[,2],sp_names,col="blue")
# Add arrows (vectors) for variables
arrows(0,0,constraints[,1],constraints[,2],lwd=1,length=0.1, col="red")
# Label the variables - the *1.15 tries to prevent the labels from overlapping arrow tips
text(constraints[,1]*1.15,constraints[,2]*1.15,rownames(constraints),col="red")

# Biplot with species and vectors
plot(species_matrix[,1],species_matrix[,2],col="black",cex=0.1,xlab=x_axis_text,ylab=y_axis_text)
abline(h=0,lty="dotted",col="gray")
abline(v=0,lty="dotted",col="gray")
text(species_matrix[,1],species_matrix[,2],sp_names,col="blue")
arrows(0,0,constraints[,1]*0.85,constraints[,2]*0.85,lwd=1,length=0.2, col="red")
text(constraints[,1],constraints[,2],rownames(constraints),col="red")

# Same, but use ordipointlabel to minimize overlap among points
plot(species_matrix[,1],species_matrix[,2],col="black",cex=0.1,xlab=x_axis_text,ylab=y_axis_text)
abline(h=0,lty="dotted",col="gray")
abline(v=0,lty="dotted",col="gray")
ordipointlabel(model_fit,add=T,display="species",select=plot_species,col="blue", font=c(3,1))
arrows(0,0,constraints[,1]*0.85,constraints[,2]*0.85,lwd=1,length=0.2, col="red")
text(constraints[,1],constraints[,2],rownames(constraints),col="red")