library(vegan)
data(dune)
data(dune.env)

#Atšķirības (dissimilarities)
d.euclidean<-dist(dune)
d.manhattan<-dist(dune,method="manhattan")
d<-d.bray<-vegdist(dune)
d.gower<-vegdist(dune,method="gower")
d.Jaccard<-2*d.bray/(1+d.bray)
d<-d.bray

#Klāsteranalīze un dendrogrammu attēli
klast1 <- hclust(d, method="single")
plot(klast1)
plot(klast1,hang=-1)

klast2 <- hclust(d, method="complete")
klast3 <- hclust(d, method="average")

par(mfrow=c(1,3))
plot(klast1,hang=-1)
plot(klast2,hang=-1)
plot(klast3,hang=-1)

#Nodala trīs klāsterus
plot(klast1,hang=-1)
rect.hclust(klast1,3)
plot(klast2,hang=-1)
rect.hclust(klast2,3)
plot(klast3,hang=-1)
rect.hclust(klast3,3)

#Piešķir klāsteru numurus katram parauglaukumam
k2<-cutree(klast2,3)
#Parauglaukumu skaits katrā klāsterī
table(k2)

k1<-cutree(klast1,3)
k3<-cutree(klast3,3)
table(k1)
table(k3)

#Parāda atšķirības klasifikācijas veidos
table(k2,k1)
table(k2,k3)

#Parāda veģetācijas datus 
dune[k2==1,]
dune[k2==2,]
dune[k2==3,]

#Parāda veģetācijas datus kompaktā veidā, sakārtotus pa klāsteriem
vegemite(dune,k2)

#k-means klāsteranalīze
km<-kmeans(d,3)
#Piešķir klāsteru numurus katram parauglaukumam
k4<-km$cluster

#Izveido failu ar analīzes rezultātiem
write.table(cbind(k1,k2,k3,k4),"Dune_klasteri.txt",sep="\t")

library(labdsv)
#Sugu sastopamība klāsteros
sastopamiba<-const(dune,k2)
write.table(sastopamiba,"Dune_klast_sastop.txt",sep="\t")

#Indikatorsugu analīze
library(indicspecies)
ind<-multipatt(dune,k2,func = "IndVal.g")
summary(ind)

#Nodala dažādu klāsteru skaitu
c2<-cutree(klast2,2)
c3<-cutree(klast2,3)
c4<-cutree(klast2,4)
c5<-cutree(klast2,5)

#Indikatorsugu analīze
ind2<-multipatt(dune,c2,func = "IndVal.g")
ind3<-multipatt(dune,c3,func = "IndVal.g")
ind4<-multipatt(dune,c4,func = "IndVal.g")
ind5<-multipatt(dune,c5,func = "IndVal.g")

#Indikatorsugu indikatorvērtības un būtiskums klāsteriem
s2<-ind2$sign
s3<-ind3$sign
s4<-ind4$sign
s5<-ind5$sign

#Atlasa tikai būtiskās indikatorsugas
s2a<-s2[s2$p.value<=0.05&!is.na(s2$p.value),]
s3a<-s3[s3$p.value<=0.05&!is.na(s3$p.value),]
s4a<-s4[s4$p.value<=0.05&!is.na(s4$p.value),]
s5a<-s5[s5$p.value<=0.05&!is.na(s5$p.value),]

#Indikatorsugu skaits
#Nk - klāsteru skaits
#N1 - tikai uz vienu klāsteri attiecināmo indikatorsugu skaits
#N2 - visu indikatorsugu skaits
#################################
Nk<-seq(from=2,to=5,by=1)
N1<-N2<-rep(0,4)
N2[1]<-nrow(s2a)
N2[2]<-nrow(s3a)
N2[3]<-nrow(s4a)
N2[4]<-nrow(s5a)
N1[1]<-nrow(s2a[rowSums(s2a[,1:2])==1,])
N1[2]<-nrow(s3a[rowSums(s3a[,1:3])==1,])
N1[3]<-nrow(s4a[rowSums(s4a[,1:4])==1,])
N1[4]<-nrow(s5a[rowSums(s5a[,1:5])==1,])

#Grafiks ar indikatorsugu skaitu atkarībā no klāsteru skaita
plot(c(2,5),c(0,20),type="n",ann=FALSE,axes=FALSE)
axis(side=1,cex=0.7,at=Nk,labels=Nk)
axis(side=2,cex=0.7)
mtext("Klāsteru skaits",side=1,cex=1,line=2)
mtext("Indikatorsugu skaits",side=2,cex=1,line=2)
lines(Nk,N1,col="red")
lines(Nk,N2,col="blue")
points(Nk,N1,pch=22,bg="red")
points(Nk,N2,pch=24,bg="blue")
legend(x="bottomright",legend=c("Viena klāstera indikatorsugas","Visas indikatorsugas"),
	lty=1,col=c("red","blue"),pch=c(22,24),pt.bg=c("red","blue"))
box()
	
#Dati uzdevumam
veg<-read.table("VegR.txt",header=TRUE,sep="\t")