#1 
library(readxl)
library(vegan)

pako.env<-read_excel("pako_opr.xlsx", sheet=2)
unique(pako.env$hab)


rda.1<-rda(decostand(pako, method="hell")~hab, data=pako.env)
rda.1
# Call: rda(formula = decostand(pako, method =
# "hell") ~ hab, data = pako.env)
# 
#               Inertia Proportion Rank
# Total          0.3778     1.0000     
# Constrained    0.1757     0.4651    4
# Unconstrained  0.2021     0.5349   22
# Inertia is variance 

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.17574 4.7831  1e-04 ***
# Residual 22  0.20208                  
# ---
# Signif. codes:  
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

RsquareAdj(rda.1)
# $r.squared
# [1] 0.4651425
# 
# $adj.r.squared
# [1] 0.3678956


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

gd<-goodness(rda.1, choices=c(1,2), summarize=T)
hist(gd)
sel<-gd>0.3

?plot.cca
aa<-ordiplot(rda.1, display=c("sp", "cn"), type="n", scaling="symmetric",
             correlation=T)
text(aa, what="sp", arrow=T, labels = labs, select=sel, cex=0.6)
text(aa, what="sp", arrow=T, labels = NULL, select=!sel, col="grey")
text(aa, what="cen", col=2)

#2
summary(pako.env)
dim(pako.env)
pako.en<-pako.env[,4:12]

library(psych)
pairs.panels(pako.en)
pako.en$velSD<-log(pako.en$velSD)
pako.en$shear_vel<-log(pako.en$shear_vel)
pako.en$roughness<-log(pako.en$roughness)

rda.0<-rda(decostand(pako, method="hell")~1, data=pako.en)
rda.f<-rda(decostand(pako, method="hell")~depth+velocity + 
             velSD + froude + Re + shear_vel + organic + 
    pom + roughness, data=pako.en)

as.formula(pako.en)
a1<-add1(rda.0, .~depth+velocity + velSD + froude + Re + shear_vel + organic + 
    pom + roughness, test="permutation")

p.adjust(a1$`Pr(>F)`)

rda.1<-update(rda.0, .~.+shear_vel)
a2<-add1(rda.1, .~depth+velocity + velSD + froude + Re + shear_vel + organic + 
    pom + roughness, test="permutation")

p.adjust(a2$`Pr(>F)`)

rda.2<-update(rda.1, .~.+velSD)
add1(rda.2, .~depth+velocity + velSD + froude + Re + shear_vel + organic + 
    pom + roughness, test="permutation")

rda.st.1<-ordistep(rda.0, rda.f, direction="forward")

anova(rda.f)
rda.st.2<-ordiR2step(rda.0, rda.f, direction="forward")


anova(rda.st.2)
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = decostand(pako, method = "hell") ~ shear_vel + froude, data = pako.en)
#          Df Variance      F Pr(>F)    
# Model     2  0.11977 5.5693  0.001 ***
# Residual 24  0.25806                  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

anova(rda.st.2, by="terms")
# Permutation test for rda under reduced model
# Terms added sequentially (first to last)
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = decostand(pako, method = "hell") ~ shear_vel + froude, data = pako.en)
#           Df Variance      F Pr(>F)    
# shear_vel  1 0.092088 8.5645  0.001 ***
# froude     1 0.027678 2.5741  0.037 *  
# Residual  24 0.258057                  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1



gd<-goodness(rda.st.2, choices=c(1,2), summarize=T)
hist(gd)
sel<-gd>0.3

aa<-ordiplot(rda.st.2, display=c("sp", "bp"), type="n", scaling="species",
             correlation=T)
text(aa, what="sp", arrow=T, labels = labs, select=sel, cex=0.6)
text(aa, what="sp", arrow=T, labels = NULL, select=!sel, col="grey")
text(aa, what="bip", col=2)

### Variation partitioning
vp.1<-varpart(decostand(pako, method="hell"), ~pako.env$hab, ~pako.env$shear_vel+pako.env$froude)
vp.1
# Partition of variance in RDA 
# 
# Call: varpart(Y = decostand(pako, method = "hell"), X =
# ~pako.env$hab, ~pako.env$shear_vel + pako.env$froude)
# 
# Explanatory tables:
# X1:  ~pako.env$hab
# X2:  ~pako.env$shear_vel + pako.env$froude 
# 
# No. of explanatory tables: 2 
# Total variation (SS): 9.8234 
#             Variance: 0.37782 
# No. of observations: 27 
# 
# Partition table:
#                      Df R.squared Adj.R.squared Testable
# [a+c] = X1            4   0.46514       0.36790     TRUE
# [b+c] = X2            2   0.26848       0.20752     TRUE
# [a+b+c] = X1+X2       6   0.52847       0.38701     TRUE
# Individual fractions                                    
# [a] = X1|X2           4                 0.17949     TRUE
# [b] = X2|X1           2                 0.01911     TRUE
# [c]                   0                 0.18841    FALSE
# [d] = Residuals                         0.61299    FALSE
# ---
# Use function ‘rda’ to test significance of fractions of interest

plot(vp.1)

anova(rda(decostand(pako, method="hell")~pako.env$hab+
            Condition(pako.env$shear_vel+pako.env$froude)))
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = decostand(pako, method = "hell") ~ pako.env$hab + Condition(pako.env$shear_vel + pako.env$froude))
#          Df Variance      F Pr(>F)    
# Model     4 0.098228 2.7568  0.001 ***
# Residual 20 0.178156                  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

anova(rda(decostand(pako, method="hell")~Condition(pako.env$hab)+
            pako.env$shear_vel+pako.env$froude))
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = decostand(pako, method = "hell") ~ Condition(pako.env$hab) + pako.env$shear_vel + pako.env$froude)
#          Df Variance     F Pr(>F)
# Model     2 0.023926 1.343  0.156
# Residual 20 0.178156

### 3
dca.1<-decorana(log1p(bk.spe))

rda.1<-rda(log1p(bk.spe)~bk.env$hemiparasite)
anova(rda.1)
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = log1p(bk.spe) ~ bk.env$hemiparasite)
#          Df Variance      F Pr(>F)    
# Model     3    9.450 1.9763  0.001 ***
# Residual 22   35.065                  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

names(bk.env)
bk.soil<-bk.env[,6:14]
pairs.panels(bk.soil)
bk.soil$PMehl3<-log(bk.soil$PMehl3)
bk.soil$Ca<-log(bk.soil$Ca)
bk.soil$Mg<-log(bk.soil$Mg)
bk.soil$K<-log(bk.soil$K)
bk.soil$Na<-log(bk.soil$Na)


rda.0<-rda(log1p(bk.spe)~1, data=bk.soil)
rda.f<-rda(log1p(bk.spe)~
  pHw +PMehl3 + LOI + N + Ca + Mg + K + Na + C.N,
  data=bk.soil)
anova(rda.f)
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = log1p(bk.spe) ~ pHw + PMehl3 + LOI + N + Ca + Mg + K + Na + C.N, data = bk.soil)
#          Df Variance      F Pr(>F)    
# Model     9   20.948 1.5802  0.001 ***
# Residual 16   23.567                  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


rda.st<-ordiR2step(rda.0, rda.f)
anova(rda.st, by="terms")
# Permutation test for rda under reduced model
# Terms added sequentially (first to last)
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = log1p(bk.spe) ~ LOI + pHw, data = bk.soil)
#          Df Variance      F Pr(>F)    
# LOI       1    6.740 4.4699  0.001 ***
# pHw       1    3.097 2.0538  0.009 ** 
# Residual 23   34.678                  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


vp.2<-varpart(log1p(bk.spe), ~bk.env$hemiparasite, 
              ~bk.soil$LOI+bk.soil$pHw)
# Partition of variance in RDA 
# 
# Call: varpart(Y = log1p(bk.spe), X = ~bk.env$hemiparasite,
# ~bk.soil$LOI + bk.soil$pHw)
# 
# Explanatory tables:
# X1:  ~bk.env$hemiparasite
# X2:  ~bk.soil$LOI + bk.soil$pHw 
# 
# No. of explanatory tables: 2 
# Total variation (SS): 1112.9 
#             Variance: 44.515 
# No. of observations: 26 
# 
# Partition table:
#                      Df R.squared Adj.R.squared Testable
# [a+c] = X1            3   0.21229       0.10487     TRUE
# [b+c] = X2            2   0.22096       0.15322     TRUE
# [a+b+c] = X1+X2       5   0.37521       0.21901     TRUE
# Individual fractions                                    
# [a] = X1|X2           3                 0.06579     TRUE
# [b] = X2|X1           2                 0.11414     TRUE
# [c]                   0                 0.03908    FALSE
# [d] = Residuals                         0.78099    FALSE
# ---
# Use function ‘rda’ to test significance of fractions of interest


plot(vp.2)

anova(rda((log1p(bk.spe)~bk.env$hemiparasite+ 
              Condition(bk.soil$LOI+bk.soil$pHw))))
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = log1p(bk.spe) ~ bk.env$hemiparasite + Condition(bk.soil$LOI + bk.soil$pHw))
#          Df Variance      F Pr(>F)   
# Model     3   6.8661 1.6458  0.002 **
# Residual 20  27.8124                 
# ---
# Signif. codes:  
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
anova(rda((log1p(bk.spe)~Condition(bk.env$hemiparasite)+ 
              bk.soil$LOI+bk.soil$pHw)))
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = log1p(bk.spe) ~ Condition(bk.env$hemiparasite) + bk.soil$LOI + bk.soil$pHw)
#          Df Variance      F Pr(>F)    
# Model     2   7.2524 2.6076  0.001 ***
# Residual 20  27.8124                  
# ---
# Signif. codes:  
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1