# Load libraries
library(vegan)
library(mvabund)
library(indicspecies)
library(Hmisc)

# Load data
data(tikus)
pre_post<-tikus$x$time!=81
community<-tikus$abund

# remove species with zero abundance
community<-community[,apply(community,2,sum)>0]

# indicator species analysis
indsp<-multipatt(community,cluster = pre_post,duleg = T,func = "IndVal.g")
indsp

# list only the significant indicators
indsp$sign[indsp$sign$p.value<0.05,]

# nmds
nmds<-metaMDS((community)^0.5,K=2,autotransform = F,trymax = 500)
# plot nmds
plot(nmds$points,xlab="Axis I",ylab="Axis II",type="n")
text(nmds$points,labels=tikus$x$time,col=c("red","black")[as.factor(pre_post)])
# polygon around pre samples
p<-nmds$points[pre_post==F,]
ch<-chull(p[,1],p[,2])
polygon(p[ch,1],p[ch,2],col = rgb(1,0,0,0.25))

# simper
sim<-simper(community,pre_post)
summary(sim)

# add the ten most incfuential species to the plot
top_ten<-sim$FALSE_TRUE$species[sim$FALSE_TRUE$ord][1:10]
top_ten_scores<-nmds$species[rownames(nmds$species) %in% top_ten,]
text(top_ten_scores,labels=rownames(top_ten_scores),col="blue")