library(readxl)
library(cluster)
library(vegan)
library(indicspecies)
library(betapart)

jedno<-read_excel("PH_jednoleta_export_juice.xlsx", sheet=1)
vytr<-read_excel("PH_vytrvala_komplet_export_juice.xlsx", sheet=1)

jedno.rows<-jedno[,1]
jedno<-jedno[,-1]# Vymazani prvniho sloupce (kody vzorku) z dat

specnumber(jedno, MARGIN=2)# Frekvence druhu
table(specnumber(jedno, MARGIN=2))# Summary frekvenci

sum(jedno>0)# Celkovy pocet vyskytu

jedno.2<-jedno[colSums(jedno>0)>1]# Vytvoreni noveho data.frame bez  jednotlivych vyskytu
jedno.d<-vegdist(sqrt(jedno.2), method="bray")# matice nepodobnosti - Bray Curtis

### PAM

#Smycka pro testovani ruzneho k - uklada se Adjusted Rsq a prumerna sirka siluety
pam.k<-data.frame(k=NA, Rsq=NA, sil=NA)
for (k in 2:10){
pam.1<-pam(jedno.d, k=k)
pam.k[k-1,1]<-k

rs<-RsquareAdj(capscale(sqrt(jedno.2)~as.factor(pam.1$clustering)))
pam.k[k-1,2]<-rs$adj.r.squared

sil<-data.frame(silhouette(pam.1))
pam.k[k-1, 3]<-mean(sil$sil_width)
}

pam.k
plot(Rsq~k, data=pam.k)
plot(sil~k, data=pam.k)

pam.7<-pam(jedno.d, k=7)
plot(silhouette(pam.7))

palette(palette.colors(n=7, palette="R4"))

pco.1<-capscale(sqrt(jedno.2)~1, method="bray")
aa<-ordiplot(pco.1, type="n", display = c("si"), 
             scaling=1)
points(aa, what="si")
ordispider(aa, groups = pam.7$clustering, col=1:7)



specnumber(jedno, MARGIN=1)# Ukazuje velkou nerovnomernost druhove bohatosti
# Na takova data by se mohl hodit Simpsonuv index nepodobnosti, resp. jeho kvantitativni podoba - balancovane zmeny abundance
# Podrobnosti v Baselga, A. (2016). Partitioning abundancebased multiplesite dissimilarity into components: balanced variation in abundance and abundance gradients. 
# Methods in Ecology and Evolution, https://doi.org/10.1111/2041210X.12693 

?beta.pair.abund
jedno.d<-beta.pair.abund(sqrt(jedno.2))#Computes 
class(jedno.d[[1]])# Tohle obsahuje matici balanced change in abundance

# Testovani PAM pro ruzne k
pam.k<-data.frame(k=NA, Rsq=NA, sil=NA)
k<-2
for (k in 2:10){
pam.1<-pam(jedno.d[[1]], k=k)
pam.k[k-1,1]<-k

rs<-RsquareAdj(capscale(jedno.d[[1]]~as.factor(pam.1$clustering)))
pam.k[k-1,2]<-rs$adj.r.squared

sil<-data.frame(silhouette(pam.1$clustering, jedno.d[[1]]))
pam.k[k-1, 3]<-mean(sil$sil_width)
}
plot(Rsq~k, data=pam.k)
plot(sil~k, data=pam.k)

pam.4<-pam(jedno.d[[1]], k=4)
plot(silhouette(pam.4))

pco.1<-capscale(jedno.d[[1]]~1)
aa<-ordiplot(pco.1, type="n", display = c("si", "sp"), 
             scaling=1)
points(aa, what="si", col=pam.4$clustering, pch=16)
ordihull(aa, groups = pam.4$clustering, col=1:4)

#Summary phi-koeficientu pro jednotlive klastry
summary(multipatt(decostand(jedno.2, method="pa"), pam.4$clustering, 
          duleg=T,
          func="r.g"))

# Cluster analysis - hierarchicka
hc.1<-hclust(sqrt(jedno.d[[1]]), method="ward.D2") #Bezi na  matici balanced changes
plot(hc.1)

# Testovani orezu dendrogramu do rucnych poctu  skupin
hc.k<-data.frame(k=NA, Rsq=NA, sil=NA)
k<-2
for (k in 2:10){
ct.1<-cutree(hc.1, k=k)
hc.k[k-1,1]<-k

rs<-RsquareAdj(capscale(sqrt(jedno.d[[1]])~as.factor(ct.1)))
hc.k[k-1,2]<-rs$adj.r.squared

sil<-data.frame(silhouette(ct.1, sqrt(jedno.d[[1]])))
hc.k[k-1, 3]<-mean(sil$sil_width)
}

plot(Rsq~k, data=hc.k)
plot(sil~k, data=hc.k)

gr.5<-cutree(hc.1,k=5)
plot(silhouette(gr.5, jedno.d[[1]]))

pco.1<-capscale(sqrt(jedno.d[[1]])~1)
aa<-ordiplot(pco.1, type="n", display = c("si", "sp"), 
             scaling=1)
points(aa, what="si", col=gr.5, pch=16)
ordihull(aa, groups = gr.5, col=1:7)

summary(multipatt(decostand(jedno.2, method="pa"), gr.5, duleg=T, func="r.g"))

#Cviceni 8 - 1b na datech Petra Hubatky

pam.5<-pam(vegdist(log1p(jedno.2), dist="bray"), k=5)
cl.5s<-hclust(vegdist(log1p(jedno.2), dist="bray"), method="single")
cl.5c<-hclust(vegdist(log1p(jedno.2), dist="bray"), method="complete")
cl.5a<-hclust(vegdist(log1p(jedno.2), dist="bray"), method="average")

plot(cl.5s)
plot(cl.5c)
plot(cl.5a)

summary(multipatt(decostand(jedno.2, method="pa"), cutree(cl.5c, k=5),  duleg=T))

pco.1<-capscale(log1p(jedno.2)~1, method="bray", sqrt.dist=T)
aa<-ordiplot(pco.1, type="n", display = c("si"), 
             scaling=1)
points(aa, what="si", col=cutree(cl.5c, k=5), pch=16)
ordihull(aa, groups = cutree(cl.5c, k=5), col=1:7)

RsquareAdj(
capscale(log1p(jedno.2)~as.factor(cutree(cl.5c, k=5)), method="bray", sqrt.dist=T)
)