library(vegan)


## data jsou uložená v d.spe

summary(d.spe)

pca.1<-rda(decostand(d.spe, method="hellinger"))
pca.1 # summary pca analýzy
# Call: rda(X = decostand(d.spe, method = "hellinger"))
# 
#               Inertia Rank
# Total          0.5023     
# Unconstrained  0.5023   27
# Inertia is variance 
# 
# Eigenvalues for unconstrained axes:
#     PC1     PC2     PC3     PC4     PC5     PC6     PC7     PC8 
# 0.24913 0.06592 0.04615 0.03723 0.02125 0.01662 0.01477 0.01297 
# (Showing 8 of 27 unconstrained eigenvalues)

eigenvals(pca.1)/sum(eigenvals(pca.1)) # vysvětlená variabilita jednotlivých os
#       PC1       PC2       PC3       PC4       PC5       PC6       PC7       PC8 
# 0.4959363 0.1312153 0.0918782 0.0741173 0.0422980 0.0330895 0.0294018 0.0258249 
#       PC9      PC10      PC11      PC12      PC13      PC14      PC15      PC16 
# 0.0209811 0.0132689 0.0100335 0.0084766 0.0055073 0.0052003 0.0029957 0.0027615 
#      PC17      PC18      PC19      PC20      PC21      PC22      PC23      PC24 
# 0.0020652 0.0015557 0.0009454 0.0008487 0.0005598 0.0004356 0.0002752 0.0001744 
#      PC25      PC26      PC27 
# 0.0001062 0.0000294 0.0000183 

screeplot(pca.1, bstick=T) # zobrazí screeplot s broken stick modelem 
abline(h=mean(eigenvals(pca.1)), lty=2) # přidá do grafu čáru zobrazující průměrné
# eigenvalue přes všechny osy



## Ordinační diagram pro první a druhou osu
aa<-ordiplot(pca.1, display=c("sites", "species"), type='n', scaling=2, 
             xlab="PC1 (49.6%)", ylab="PC2 (13.1%)")## nakreslí prázdný diagram
text(aa, what="species",arrows=T, length=0.04, cex=0.7) # přidá do diagramu
# šipky druhů jejich popisky 
points(aa, what="sites", col="darkgrey", pch=16) # přidá do diagramu skóre lokalit


## Ordinační diagram pro první a třetí osu
aa<-ordiplot(pca.1, display=c("sites", "species"), type='n', scaling=2,
             choices=c(1,3), #volba os
             xlab="PC1 (49.6%)", ylab="PC3 (13.1%)")## nakreslí prázdný diagram
text(aa, what="species",arrows=T, length=0.04, cex=0.7) # přidá do diagramu
# šipky druhů jejich popisky 
points(aa, what="sites", col="grey", pch=16) # přidá do diagramu skóre lokalit

aa<-ordiplot(pca.1, display=c("sites", "species"), type='n', scaling=1,
             choices=c(1,3), #volba os
             xlab="PC1 (49.6%)", ylab="PC3 (13.1%)")## nakreslí prázdný diagram
text(aa, what="species",arrows=T, length=0.04, cex=0.7) # přidá do diagramu
# šipky druhů jejich popisky 
points(aa, what="sites", col="grey", pch=16) # přidá do diagramu skóre lokalit

pdf("ordiplot.pdf", 6,6)
aa<-ordiplot(pca.1, display=c("sites", "species"), type='n', scaling=2,
             choices=c(1,3), #volba os
             xlab="PC1 (49.6%)", ylab="PC3 (13.1%)",
             xlim=c(-1,1))## nakreslí prázdný diagram
text(aa, what="species",arrows=T, length=0.04, cex=0.7) # přidá do diagramu
# šipky druhů jejich popisky 
points(aa, what="sites", col="grey", pch=16) # přidá do diagramu skóre lokalit
dev.off()

### 8
bk<-read.delim("clipboard", row.names = 1)
summary(bk)
bk<-bk[,-2]
summary(bk)

pca.bk<-rda(bk, scale=T)
# pca.bk
# Call: rda(X = bk, scale = T)
# 
#               Inertia Rank
# Total              13     
# Unconstrained      13   13
# Inertia is correlations 
# 
# Eigenvalues for unconstrained axes:
#   PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8   PC9  PC10  PC11  PC12  PC13 
# 4.342 2.116 1.564 1.417 0.948 0.755 0.624 0.424 0.398 0.254 0.111 0.043 0.005 

screeplot(pca.bk, bstick = T)
abline(h=mean(eigenvals(pca.bk)), lty=2)

bb<-ordiplot(pca.bk, display=c("species", "sites"), scaling=3, 
             type="n")
text(bb, what="sp", arrows=T, length=0.05)
points(bb, what="si", pch="+", col="red")


# Pokus o analýzu se standardizací po vzorcích 
pca.bk<-rda(decostand(bk, method="normalize"), scale=T)
bb<-ordiplot(pca.bk, display=c("species", "sites"), scaling=3, 
             type="n")
text(bb, what="sp", arrows=T, length=0.05)
points(bb, what="si", pch="+", col="red")

# Jak vidíte, tak se ordinační diagram dost zádadně změnil a vede 
# k dost odlišným interpretacím. To je proto, že standardizace po
# vzorcích nedává u takovýchto dat smysl kvůli odlišným jednotkám a 
# měřítkům, tj. ze stejného důvodu proč vždy musíme provst standardizaci po 
# proměnných (druzích).