#### Prima ordinace
library(vegan)

#### 1. Jednoducha prima ordinace
# druhova data jsou v pako.spe
# promenne prostredi jsou v pako. env

rda.1<-rda(decostand(log(pako.spe+1), method="hellinger")~hab, data=pako.env)
rda.1 #zobrazeni eigenvalues. Je videt, ze typ biotopu
# vysvetluje hodne variability ve srovnani se zbytkovymi PCA osami

# Call: rda(formula = decostand(log(pako.spe + 1), method = "hellinger") ~
# hab, data = pako.env)
# 
#               Inertia Proportion Rank
# Total          0.2862     1.0000     
# Constrained    0.1135     0.3966    4
# Unconstrained  0.1727     0.6034   22
# Inertia is variance 
# 
# Eigenvalues for constrained axes:
#    RDA1    RDA2    RDA3    RDA4 
# 0.07363 0.02002 0.01294 0.00690 
# 
# Eigenvalues for unconstrained axes:
#     PC1     PC2     PC3     PC4     PC5     PC6     PC7     PC8 
# 0.03743 0.01976 0.01748 0.01402 0.01112 0.00884 0.00789 0.00721 
# (Showing 8 of 22 unconstrained eigenvalues)

anova(rda.1) # test vyznamnosti modelu
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = decostand(log(pako.spe + 1), method = "hellinger") ~ hab, data = pako.env)
#          Df Variance      F Pr(>F)    
# Model     4  0.11349 3.6148  0.001 ***
# Residual 22  0.17268                  

RsquareAdj(rda.1) #Vysvetlena variabilita modelu
# $r.squared ### Normalni R2 modelu
# [1] 0.396584
# 
# $adj.r.squared    ### Upravene R2 modelu
# [1] 0.286872

#Ordinacni diagram
aa<-ordiplot(rda.1, display=c("cn", "sp"), type="n") #Vytvotreni prazdneho diagramu,
# s rozmery pro prediktor a druhy
# Prediktor je v tomto pripade faktor. Zobrazujeme proto jako centroidy, ne jako sipky.

# vygenerovani zkratek druhovych jmen
sh.labels<-make.cepnames(names(pako.spe))
substr(sh.labels,5,5)<-toupper(substr(sh.labels,5,5))

text(aa, what="cen", arrows=F, length=0.1, col=2)# Pridani centroidu promenne prostredi
# !!!! Pozor, parametr what musi byt "centroids" nebo zkratka, treba "cen". 
# what = "cn" nefunguje.

text(aa, what="sp", labels=sh.labels, arrows=T, length=0.05, cex=0.7)# pridani
#druhovych sipek

# 2. Postupny vyber modelu 
rda.20<-rda(decostand(log(pako.spe+1), method="hellinger")~1, data=pako.env[,-c(1,2)])
rda.2<-rda(decostand(log(pako.spe+1), method="hellinger")~., data=pako.env[,-c(1,2)])
anova(rda.2)

add1(rda.20, scope=rda.2, test="permutation", permutations=how(nperm=999))
#           Df     AIC      F Pr(>F)    
# <none>       -32.801                  
# depth      1 -32.170 1.3009  0.184    
# velocity   1 -37.018 6.4731  0.001 ***
# velSD      1 -34.543 3.7166  0.001 ***
# froude     1 -37.156 6.6352  0.001 ***
# Re         1 -35.162 4.3830  0.001 ***
# shear_vel  1 -35.563 4.8219  0.001 ***
# organic    1 -32.930 2.0518  0.042 *  
# pom        1 -36.157 5.4863  0.001 ***
# roughness  1 -31.721 0.8666  0.534    
# ---

### Manualni forward selection prediktoru
rda.21<-update(rda.20, .~.+froude)
add1(rda.21, scope=rda.2, test="permutation", permutations=how(nperm=999))
#           Df     AIC      F Pr(>F)  
# <none>       -37.156                
# depth      1 -36.349 1.0843  0.306  
# velocity   1 -36.812 1.5174  0.074 .
# velSD      1 -36.445 1.1734  0.224  
# Re         1 -36.991 1.6872  0.047 *
# shear_vel  1 -36.103 0.8560  0.632  
# organic    1 -36.611 1.3281  0.139  
# pom        1 -36.710 1.4214  0.118  
# roughness  1 -36.474 1.2000  0.237  
# ---

rda.21<-update(rda.21, .~.+Re)##Nejlepe (a prukazne) vychazi Re
add1(rda.21, scope=rda.2, test="permutation", permutations=how(nperm=999))
anova(rda.21, by="terms")

## Automaticka selekce pomoci upraveho R2 dle Blanchet et al 2008.
rda.step<-ordiR2step(rda.21, scope=rda.2)

##Automaticka backward selekce
rda.step<-ordistep(rda.2, scope=rda.20, direction = "backward")

#Ordinacni diagram
aa<-ordiplot(rda.21, display=c("bp", "sp"), type="n") #Vytvotreni prazdneho diagramu,
# s rozmery pro prediktor a druhy
# Prediktor je v tomto pripade ciselna promenna. Zobrazujeme proto jako sipky.

# vygenerovani zkratek druhovych jmen
sh.labels<-make.cepnames(names(pako.spe))
substr(sh.labels,5,5)<-toupper(substr(sh.labels,5,5))

text(aa, what="bip", arrows=T, length=0.1, col=2)# Pridani sipky promennych prostredi - 
# !!!! Pozor, parametr what musi byt "biplot" nebo zkratka, treba "bip". 
# what = "bp" nefunguje. 

text(aa, what="sp", labels=sh.labels, arrows=T, length=0.05, cex=0.7)# pridani
#druhovych sipek


# 3. Rozklad variability (Variation partitioning)
# Druhova data jsou v spe
# Promenne prostredi v env

#3.1 Vyber vyznamnych pudnich parametru

rda.s0<-rda(decostand(log(spe+1), method="hellinger")~1, data=env)
rda.sf<-rda(decostand(log(spe+1), method="hellinger")~+pHw + PMehl3 +
     LOI + N + Ca + Mg + K + Na+C.N, data=env)
rda.soil<-ordiR2step(rda.s0, rda.sf)

#3.2 Rozklad variance
vp.1<-varpart(log(spe+1), ~+hemiparasite, ~+pHw + LOI, data=env, transfo = "hel")
vp.1
# Partition of variance in RDA 
# 
# Call: varpart(Y = log(spe + 1), X = ~+hemiparasite, ~+pHw + LOI, data = env,
# transfo = "hel")
# Species transformation:  hellinger
# Explanatory tables:
# X1:  ~+hemiparasite
# X2:  ~+pHw + LOI 
# 
# No. of explanatory tables: 2 
# Total variation (SS): 13.881 
#             Variance: 0.55523 
# No. of observations: 26 
# 
# Partition table:
#                      Df R.squared Adj.R.squared Testable
# [a+b] = X1            3   0.19474       0.08493     TRUE
# [b+c] = X2            2   0.19841       0.12871     TRUE
# [a+b+c] = X1+X2       5   0.34177       0.17722     TRUE
# Individual fractions                                    
# [a] = X1|X2           3                 0.04851     TRUE
# [b]                   0                 0.03642    FALSE
# [c] = X2|X1           2                 0.09228     TRUE
# [d] = Residuals                         0.82278    FALSE
# ---
# Use function ‘rda’ to test significance of fractions of interest
# Nakresleni mnozin s podily vysvetlene variability
plot(vp.1)

## Testovani jednotlivych komponent pomoci parcialnich RDA
# Parcialni efekt druhu poloparazita

rda.hemi<-rda(spe~+hemiparasite+Condition(pHw + LOI), data=env)
anova(rda.hemi)
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = spe ~ +hemiparasite + Condition(pHw + LOI), data = env)
#          Df Variance      F Pr(>F)  
# Model     3   217.19 1.5651  0.016 *
# Residual 20   925.16    

rda.soil<-rda(spe~Condition(hemiparasite)+pHw + LOI, data=env)
anova(rda.soil)
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = spe ~ Condition(hemiparasite) + pHw + LOI, data = env)
#          Df Variance      F Pr(>F)    
# Model     2   252.40 2.7282  0.001 ***
# Residual 20   925.16                  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1   

rda.all<-rda(spe~+hemiparasite+pHw+ LOI,
              data=env)
anova(rda.all)
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = spe ~ +hemiparasite + pHw + LOI, data = env)
#          Df Variance      F Pr(>F)    
# Model     5   518.19 2.2404  0.001 ***
# Residual 20   925.16                  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

