# 1
# Slozeni spolecenstva
pako<-read.csv2("pako_opr.csv", row.names = 1)
summary(pako)
pako[is.na(pako)]<-0

#Prediktory
pako.env<-read.delim("pako_env.txt", row.names=1)

#RDA

rda.1<-rda(decostand(pako, method="hell")~hab, data=pako.env)

rda.1
#               Inertia Proportion Rank
# Total          0.3783     1.0000     
# Constrained    0.1757     0.4644    4
# Unconstrained  0.2026     0.5356   22

RsquareAdj(rda.1)
# $r.squared
# [1] 0.4644437
# 
# $adj.r.squared
# [1] 0.3670698

anova(rda.1, permutations=how(nperm=9999))
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 9999
# 
# Model: rda(formula = decostand(pako, method = "hell") ~ hab, data = pako.env)
#          Df Variance      F Pr(>F)    
# Model     4  0.17572 4.7697  1e-04 ***
# Residual 22  0.20262  
# ---
# Signif. codes:  
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

aa<-ordiplot(rda.1, type="n", display = c("sp","cn"), 
             correlation=T)
text(aa, what="sp", arrows = T, length=0.03)
text(aa, what="cen", col="red", arrows=F)



#2
names(pako.env)

rda.0<-rda(decostand(pako, method="hell")~1, data=pako.env)
tests<-add1(rda.0,.~+depth+velocity+velSD+froude+Re+shear_vel+organic+pom+roughness,
     test="perm", permutations = how(nperm=999))
#           Df     AIC      F Pr(>F)    
# <none>       -25.262                  
# depth      1 -24.829 1.4943  0.165    
# velocity   1 -30.438 7.6114  0.001 ***
# velSD      1 -28.324 5.1560  0.001 ***
# froude     1 -30.714 7.9467  0.001 ***
# Re         1 -27.828 4.6072  0.001 ***
# shear_vel  1 -29.352 6.3259  0.001 ***
# organic    1 -24.679 1.3476  0.201    
# pom        1 -28.530 5.3863  0.001 ***
# roughness  1 -24.173 0.8585  0.519 

p.adjust(tests$`Pr(>F)`, method="holm")

rda.1<-update(rda.0, .~.+froude)
tests<-add1(rda.1,.~+depth+velocity+velSD+froude+Re+shear_vel+organic+pom+roughness,
     test="perm", permutations = how(nperm=999))

rda.full<-rda(decostand(pako, method="hell")~.,
              data=pako.env[,-c(1:2)])

rda.b<-ordistep(rda.full, direction="back")
summary(rda.b)
anova(rda.b, by="margin")
# Model: rda(formula = decostand(pako, method = "hell") ~ velocity + velSD + froude, data = pako.env[, -c(1:2)])
#          Df Variance      F Pr(>F)  
# velocity  1 0.023595 2.2049  0.048 *
# velSD     1 0.024616 2.3002  0.032 *
# froude    1 0.021392 1.9990  0.065 .
# Residual 23 0.246134               


aa<-ordiplot(rda.b, type="n", display = c("sp","bp"), 
             correlation=T)
text(aa, what="sp", arrows = T, length=0.03)
text(aa, what="bip", col="red", arrows=T, length=0.05)

#### 3
bk.env<-read.csv2("env.csv", row.names = 1)
bk.spe<-read.csv2("veg.csv", row.names = 1)

cca.full<-cca(log(bk.spe+1)~., data=bk.env[,6:14])
anova(cca.full)

cca.0<-cca(log(bk.spe+1)~1, data=bk.env[,6:14])

cca.step<-ordistep(cca.0, cca.full, direction = "forward")
anova(cca.step, by="terms")

anova(cca.step)
# Model: cca(formula = log(bk.spe + 1) ~ LOI + pHw + K, data = bk.env[, 6:14])
#          Df ChiSquare      F Pr(>F)    
# Model     3   0.69312 1.9445  0.001 ***
# Residual 22   2.61403                  
# ---
# Signif. codes:  
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

RsquareAdj(cca.step)
# $r.squared
# [1] 0.209583
# 
# $adj.r.squared
# [1] 0.1016426

pcca.1<-update(cca.step, .~.+Condition(bk.env$hemiparasite))
RsquareAdj(pcca.1)
# $r.squared
# [1] 0.1742598
# 
# $adj.r.squared
# [1] 0.08106232

anova(pcca.1)
#          Df ChiSquare      F Pr(>F)    
# Model     3    0.5763 1.7004  0.001 ***
# Residual 19    2.1465 

pcca.2<-cca(log(bk.spe+1)~hemiparasite+Condition(LOI + pHw
+ K), data=bk.env)
RsquareAdj(pcca.2)
# $r.squared
# [1] 0.1413581
# 
# $adj.r.squared
# [1] 0.04411461

anova(pcca.2)
#          Df ChiSquare      F Pr(>F)    
# Model     3   0.46749 1.3793  0.001 ***
# Residual 19   2.14654

pcca.3<-cca(log(bk.spe+1)~hemiparasite+LOI + pHw
+ K, data=bk.env)
RsquareAdj(pcca.3)
# $r.squared
# [1] 0.3509411
# 
# $adj.r.squared
# [1] 0.1442775

anova(pcca.3)
#          Df ChiSquare      F Pr(>F)    
# Model     6    1.1606 1.7122  0.001 ***
# Residual 19    2.1465  

aa<-ordiplot(pcca.3, display = c("sp", "cn", "bp"), 
             type="n", correlation=T)
points(aa, what="sp")
text(aa, what="cen", col="red", labels=unique(bk.env$hemiparasite))
text(aa, what="bip", col="red", arrows = T, length=0.05)