#####################################
#Load libraries
#####################################
library(BiodiversityR)
library(vegan)
library(MASS)
library(car)
library(ade4)
#####################################
#Load Data
#####################################

community<-read.csv(file="example.csv", row.names='Sites', sep=",", header=TRUE)

# Simple corespondence analysis
ca<-cca(community)
summary(ca)
ca$rowsum
ca$colsum

ordiplot(ca, type="t")

# copy the matrix, make one value negative
community_with_negatives<-community
community_with_negatives[3,1]<- -0.7071
# try CA with negative value in matrix
cca(community_with_negatives)

# DCA on the same dataset
dca<-decorana(community)
summary(dca)
# 2D plot with text labels
ordiplot(dca, type="t")

# Perform a DCA without automatic downweighting of rare species and without detrending
# The decorana function uses a slightly different approach than cca (reciprocal averaging, axes labeled RA)
# However, note the similarity of eigenvalues and plots of CA and the untrended DCA
dca<-decorana(community, iweigh=0, ira=1)
summary(dca)
# 2D plot with text labels
ordiplot(dca, type="t")

# Discriminant function example
# load spaeth dataset
spaeth<-read.csv("spaeth.csv",header=T,row.names=1)
# save factors
spaeth_factors<-spaeth[,1:3]
# take factors out of dataset
spaeth<-spaeth[,4:43]
# eliminate species occuring less than two times
spaeth<-spaeth[,apply(spaeth>0,2,sum)>2]
# CA of log+1 transformed data
ca<-cca(log(spaeth+1))
# Make a nice plot
point_colors<-c("red","blue","black","gold")
ax1_label<-sprintf("Axis 1 %.2f%%", ca$CA$eig[1]/sum(ca$CA$eig)*100)
ax2_label<-sprintf("Axis 2 %.2f%%", ca$CA$eig[2]/sum(ca$CA$eig)*100)
plot(ca$CA$u[,1],ca$CA$u[,2],pch=21,bg=point_colors[spaeth_factors$creek],ylab=ax2_label,xlab=ax1_label)
# Plot the species...only plot the 10 most abundant species
abundance_cutoff<-sort(apply(spaeth,2,sum),decreasing=T)[10]
plot_species<-apply(spaeth,2,sum)>abundance_cutoff
species_matrix<-subset(ca$CA$v,plot_species)
# There is a vegan function to make 8 letter species abbreviatioins, use it to shorten names
sp_names<-make.cepnames(row.names(species_matrix))
text(species_matrix,sp_names,col="black")
# Add arrows (vectors) for variables
arrows(0,0,species_matrix[,1],species_matrix[,2],lwd=1,length=0.2)
# Add a legend
creeks<-names(summary(spaeth_factors$creek))
legend("topleft",creeks,text.col=point_colors)

# Put centroids and 75% confidence ellipses around the four creeks
# dataEllips function in the car package does this
for(c in 1:length(creeks))
{
  dataEllipse(subset(ca$CA$u[,1:2],spaeth_factors$creek==creeks[c]),add=T, levels=0.75, col=point_colors[c], center.pch=17)
}

dfa<-discrimin.coa(log(spaeth+1),spaeth_factors$creek,scannf=F,nf=6)
# Make a nice plot
point_colors<-c("red","blue","black","gold")
ax1_label<-sprintf("Axis 1 %.2f%%", dfa$eig[1]/sum(dfa$eig)*100)
ax2_label<-sprintf("Axis 2 %.2f%%", dfa$eig[2]/sum(dfa$eig)*100)
plot(dfa$li[,1],dfa$li[,2],pch=21,bg=point_colors[spaeth_factors$creek],ylab=ax2_label,xlab=ax1_label)
# Plot the species...only plot the 10 most abundant species
abundance_cutoff<-sort(apply(spaeth,2,sum),decreasing=T)[10]
plot_species<-apply(spaeth,2,sum)>abundance_cutoff
species_matrix<-subset(dfa$fa,plot_species)
# There is a vegan function to make 8 letter species abbreviatioins, use it to shorten names
sp_names<-make.cepnames(row.names(species_matrix))
text(species_matrix,sp_names,col="black")
# Add arrows (vectors) for variables
arrows(0,0,species_matrix[,1],species_matrix[,2],lwd=1,length=0.2)
# Add a legend
creeks<-names(summary(spaeth_factors$creek))
legend("topleft",creeks,text.col=point_colors)

# Put centroids and 75% confidence ellipses around the four creeks
# dataEllips function in the car package does this
for(c in 1:length(creeks))
{
  dataEllipse(as.matrix(subset(dfa$li[,1:2],spaeth_factors$creek==creeks[c])),add=T, levels=0.75, col=point_colors[c], center.pch=17)
}