# This script calculates forward selection on data from Carpathian wetlands (Hajek et al.)
library (vegan)

# Loading data
vasc <- read.delim ('http://www.davidzeleny.net/anadat-r/data-download/vasc_plants.txt', row.names = 1)
chem <- read.delim ('http://www.davidzeleny.net/anadat-r/data-download/chemistry.txt', row.names = 1)

decorana (vasc)
# Call:
#   decorana(veg = vasc) 
# 
# Detrended correspondence analysis with 26 segments.
# Rescaling of axes with 4 iterations.
# 
# DCA1   DCA2    DCA3    DCA4
# Eigenvalues     0.4172 0.1657 0.11176 0.09538
# Decorana values 0.4244 0.1657 0.09826 0.08197
# Axis lengths    3.0884 2.5454 1.40922 1.64015

vasc.hell <- decostand (vasc, method = "hell")
# chem2 <- chem[,c(1,2,3,4,5,6,7,8,9,10,11,12,13,14)]
# chem2 <- chem[,1:14]
chem2 <- chem[,-15]
names (chem2)

# Run the global RDA model
rda.all <- rda (vasc.hell ~ ., data = chem2)
rda.all

# variation explained by individual axes
total.inertia <- rda.all$tot.chi
eig.constrained <- rda.all$CCA$eig
eig.unconstr <- rda.all$CA$eig

# percentage variation explained by axes
expl.constr <- eig.constrained/total.inertia
expl.unconstr <- eig.unconstr/total.inertia

# drawing this variation
barplot (c(expl.constr, expl.unconstr))

rda.0 <- rda (vasc.hell ~ 1, data = chem2)
rda.0

expl.pca <- rda.0$CA$eig/total.inertia
expl.pca

# Calculates forward selection:
fw <- ordistep (rda.0, scope = formula (rda.all), direction = 'forward')
fw

# Draw the results of RDA with variables selected by forward selection:
windows ()
plot (fw)
anova (fw)