#####################################
#Load libraries
#####################################
library(BiodiversityR)
library(labdsv)

#####################################
#Load Data
#####################################
sample_data<-read.csv(file="sample2.csv", row.names=1, header=TRUE)
# The vear and creek variables in the file will serve as factors
year<-as.factor(sample_data$year)
creek<-sample_data$creek
# We need to have the community data in its own matrix
community<-subset(sample_data,select=-c(year,creek))

#log transform community data if you want to use it, note that untransformed data is used below
community_proportion <- disttransform(community, method='profiles')

#####################################
#Run ANOSIM
#####################################
#ANOSIM testing for differences among years using a bray-curtis distance matrix
community_dist<-vegdist(community_proportion, method='bray')
year_ano<-anosim(community_dist,year, permutations=5000)
summary(year_ano)
#plot 
plot(year_ano)
# values in boxplots are mean ranks within and between groups
tapply(year_ano$dis.rank,year_ano$class.vec,mean)
# Can calculate R with these values
mean_ranks<-tapply(year_ano$dis.rank,year_ano$class.vec,mean)
between<-mean_ranks[1]
within<-mean(mean_ranks[2:4])
(between-within)/((36*35)/4)
#plot a frequency distribution ofthe permuted R values
hist(year_ano$perm)

#ANOSIM testing for differences among creeks using a bray-curtis distance matrix
community_dist<-vegdist(community_proportion, method='bray')
creek_ano<-anosim(community_dist,creek, permutations=5000)
summary(creek_ano)
#plot 
plot(creek_ano)
#plot a frequency distribution ofthe permuted R values
hist(creek_ano$perm)

#ANOSIM testing for differences among years using a qualitative distance matrix
pa_dist<- dsvdis(community, index="steinhaus")
pa_year_ano<-anosim(pa_dist,year, permutations=5000)
summary(pa_year_ano)
#plot 
plot(pa_year_ano)
#plot a frequency distribution ofthe permuted R values
hist(year_ano$perm)


#####################################
#Run Multi Response Permutation Procedure (MRPP)
#####################################
#MRPP using bray-curtis distance and testing for differences among years
community_dist<-vegdist(community_proportion, method='bray')
year_mrpp<-mrpp(community_dist,year,permutations=5000)
year_mrpp
#historgram of permuted delta scores
hist(year_mrpp$boot.deltas)


#MRPP using bray-curtis distance and testing for differences among creeks
community_dist<-vegdist(community, method='bray')
creek_mrpp<-mrpp(community_dist,creek,permutations=5000)
creek_mrpp
#historgram of permuted delta scores
hist(creek_mrpp$boot.deltas)

#####################################
# multivariate ANOVA for distance matrix
#####################################

permanova<-adonis(community_proportion ~ creek*year, method="bray", permutations=10000)
print(permanova)
#histogram for assessment of the first factor
hist(permanova$f.perms[,1])
#histogram for assessment of the second factor
hist(permanova$f.perms[,2])
#species coefficient scores (??)
print(permanova$coefficients)
# Look for species that have large differences in coefficients among factor levels.
# For example, Species7 has very large coefficient scores among years
plot(year,community$Species7)
# In contrast, Species6 has very similar coefficients for years
plot(year,community$Species6)


# The multivariate ANOVA for distance matrix with Euclidean distance and one variable should be identical to a univariate ANOVA
# community[,1] - uses just one variable (first species) from the community data
permanova<-adonis(community[,1] ~ creek*year, method="euclidean", permutations=5000)
permanova

# Function lm for a standard linear model
# Note the similarity between this and the adonis approach
model<-lm(community[,1]~creek*year)
anova(model)

#####################################
# Analysis of Molecular Variation
#####################################
library(ade4)
# packjage comes with the original Excoffier et al (1992) human dataset
data(humDNAm)
amovahum <- amova(humDNAm$samples, sqrt(humDNAm$distances), humDNAm$structures)
amovahum
randtest(amovahum,nrepet=1000)

library(pegas)
# note that when you load this, R tells you function "amova" from ade4 is being masked
# this means you are now using this version of amova. 
library(ape) # package that contains sample genetic data
data(woodmouse) # load sequence dataset for woodmouse
# calculate genetic distance with dist.dan function
d <- dist.dna(woodmouse)
# set up a factor separating the 15 mouse samples into populations A and B
pops <- factor(c(rep("A", 7), rep("B", 8)))
amova(d ~ pops, nperm = 1000)

# note that the results are basically identical to what adonis function returns
adonis(d~pops, permutations=1000)