# Load libraries
library(BiodiversityR)

# load data
spaeth = read.csv('spaeth.csv', row.names = 'sites')

# Remove the first three columns (factors) and save a community
community = spaeth[,-(1:3)]

# The year, creek, and season variables in the file will serve as factors
year<-as.factor(spaeth$year)
creek = spaeth[,'creek']
season = spaeth[,'season']

# Rarefy dataset to 50 individuals per sample (rrarefy)
rarefied_community<-rrarefy(community, sample=50)

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

# Perform ANOSIM to test for community differences among years
community_dist<-vegdist(community_proportion, method='bray')
year_ano<-anosim(community_dist, year, permutations=5000)
summary(year_ano)

# Perform ANOSIM to test for community differences among creeks
creek_ano<-anosim(community_dist, creek, permutations=5000)
summary(creek_ano)

# Perform MRP to test for community differences among years
year_mrpp<-mrpp(community_dist,year,permutations=5000)
year_mrpp

# Perform MRP to test for community differences among creeks
creek_mrpp<-mrpp(community_dist,creek,permutations=5000)
creek_mrpp

# Perform PERMANOVA to test for community differences among years
permanova<-adonis(community_proportion ~ year*season, method="bray", permutations=10000, strata=creek)
print(permanova)