#1-2
library(readxl)
bru<-read_excel("bruslarky.xlsx")
summary(bru)

#3
library(vegan)
bru.pca<-rda(bru[,4:11], scale=T)

#4
bru.pca
# Call: rda(X = bru[, 4:11], scale = T)
# 
#               Inertia Rank
# Total               8     
# Unconstrained       8    8
# Inertia is correlations 
# 
# Eigenvalues for unconstrained axes:
#   PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8 
# 7.639 0.251 0.056 0.021 0.014 0.008 0.007 0.005 

screeplot(bru.pca)

#5
par(mar=c(5,4,2,5))
aa<-ordiplot(bru.pca, type="n", display=c("si", "sp"))
points(aa, what="si", cex=0.2, pch=16)
text(aa, what="sp", arrow=T)

length(unique(bru$instar))
# ordihull(bru.pca, groups=bru$instar)
# ?palette
palette(palette.colors(6, "Okabe-Ito", alpha=0.5))
ordihull(bru.pca, groups=bru$instar,
         col=1:6,
         draw="polygon", lty=0)
legend(x=4.0, y=0, legend=1:6, fill=1:6, 
       title="Instar", xpd=T)

#6
par(mar=c(5,4,2,8))
bru.pca.6<-rda(bru[bru$instar==6,4:11], scale=T)
aa<-ordiplot(bru.pca.6, display=c("si", "sp"),
             type="n")
# points(aa, what="si", cex=0.2, pch=16)
text(aa, what="sp", arrow=T)
ordispider(bru.pca.6, groups=bru$species[
bru$instar==6], col=1:6, lwd=2)
legend(x=3.5, y=0, legend=unique(bru$species), col=1:6, 
      lwd=2, title="Species", xpd=T)
#7
pako<-read_excel("pako_opr.xlsx")
#8
summary(pako)
pako<-pako[,-1]
#9
pako[is.na(pako)]<-0
# 10
pako<-pako[,colSums(pako>0)>1]

#11
dca.pako<-decorana(log1p(pako))
dca.pako
# decorana(veg = log1p(pako)) 
# 
# Detrended correspondence analysis with 26 segments.
# Rescaling of axes with 4 iterations.
# Total inertia (scaled Chi-square): 0.9145 
# 
#                        DCA1    DCA2    DCA3    DCA4
# Eigenvalues          0.2515 0.11365 0.05767 0.03509
# Additive Eigenvalues 0.2515 0.11361 0.05727 0.03414
# Decorana values      0.2620 0.08353 0.04597 0.02029
# Axis lengths         1.9168 1.38527 1.07346 0.79291

labs<-make.cepnames(names=names(pako), seconditem = T)
substr(labs,5,5)<-toupper(substr(labs,5,5))

hist(weights(dca.pako, display="sp"))
sel<-weights(dca.pako, display="sp")>20



aa<-ordiplot(dca.pako, type = "n")
points(aa,what="si")
text(aa,what="sp", labels=labs, cex=0.6,
     select=sel)
points(aa,what="sp", cex=0.6,
     select=!sel, pch=16, col="gray")

#13
pako.nmds<-metaMDS(log1p(pako))
pako.nmds$stress
stressplot(pako.nmds)

pako.nmds.3<-metaMDS(log1p(pako), k=3)
pako.nmds.3$stress
stressplot(pako.nmds.3)

pako.nmds.4<-metaMDS(log1p(pako), k=4)
pako.nmds.4$stress
stressplot(pako.nmds.4)

sel<-colSums(pako>0)>10
table(sel)

aa<-ordiplot(pako.nmds.3, display=c("si", "sp"),
            type="n")
points(aa, what="si")
text(aa, what="sp", labels=labs, cex=0.6, 
     select=sel)