Load required packages

library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-3
library(MASS)
library(vegan3d)
## Warning: package 'vegan3d' was built under R version 3.5.2

Ordination basic method

##Non-metric Multidimensional scaling
data(varespec)
vare.dis <- vegdist(varespec)
vare.mds0 <- isoMDS(vare.dis)
## initial  value 18.026495 
## iter   5 value 10.095483
## final  value 10.020469 
## converged
stressplot(vare.mds0, vare.dis)

ordiplot(vare.mds0, type = "t")
## species scores not available

vare.mds <- metaMDS(varespec, trace = FALSE)
vare.mds
## 
## Call:
## metaMDS(comm = varespec, trace = FALSE) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(sqrt(varespec)) 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.1825658 
## Stress type 1, weak ties
## Two convergent solutions found after 20 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'wisconsin(sqrt(varespec))'
metaMDS(comm = varespec, trace = FALSE)
## 
## Call:
## metaMDS(comm = varespec, trace = FALSE) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(sqrt(varespec)) 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.1843196 
## Stress type 1, weak ties
## Two convergent solutions found after 20 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'wisconsin(sqrt(varespec))'
plot(vare.mds, type = "t")

##Community dissimilarities
data(varechem)
rankindex(scale(varechem), varespec, c("euc","man","bray","jac","kul"))
##       euc       man      bray       jac       kul 
## 0.2396330 0.2735087 0.2837910 0.2837910 0.2839834
dis <- vegdist(decostand(varespec, "norm"), "euclid")
dis <- vegdist(decostand(varespec, "hell"), "euclidean")
d <- vegdist(varespec, "bray", binary = TRUE)
d <- designdist(varespec, "(A+B-2*J)/(A+B)")
d <- designdist(varespec, "(b+c)/(2*a+b+c)", abcd=TRUE)
##Comparing ordinations: Procrustes rotation
tmp <- wisconsin(sqrt(varespec))
dis <- vegdist(tmp)
vare.mds0 <- isoMDS(dis, trace = 0)
pro <- procrustes(vare.mds, vare.mds0)
pro 
## 
## Call:
## procrustes(X = vare.mds, Y = vare.mds0) 
## 
## Procrustes sum of squares:
## 0.1563
plot(pro)

plot(pro, kind = 2)

## Eigenvector methods
vare.pca <- rda(varespec)
vare.pca
## Call: rda(X = varespec)
## 
##               Inertia Rank
## Total            1826     
## Unconstrained    1826   23
## Inertia is variance 
## 
## Eigenvalues for unconstrained axes:
##   PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8 
## 983.0 464.3 132.3  73.9  48.4  37.0  25.7  19.7 
## (Showing 8 of 23 unconstrained eigenvalues)
plot(vare.pca)

sum(apply(varespec, 2, var))
## [1] 1825.659
biplot(vare.pca, scaling = -1)

vare.pca <- rda(varespec, scale = TRUE)
vare.pca
## Call: rda(X = varespec, scale = TRUE)
## 
##               Inertia Rank
## Total              44     
## Unconstrained      44   23
## Inertia is correlations 
## 
## Eigenvalues for unconstrained axes:
##   PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8 
## 8.898 4.756 4.264 3.732 2.964 2.884 2.727 2.180 
## (Showing 8 of 23 unconstrained eigenvalues)
plot(vare.pca, scaling = 3)

dim(varespec)
## [1] 24 44
vare.ca <- cca(varespec)
vare.ca
## Call: cca(X = varespec)
## 
##               Inertia Rank
## Total           2.083     
## Unconstrained   2.083   23
## Inertia is scaled Chi-square 
## 
## Eigenvalues for unconstrained axes:
##    CA1    CA2    CA3    CA4    CA5    CA6    CA7    CA8 
## 0.5249 0.3568 0.2344 0.1955 0.1776 0.1216 0.1155 0.0889 
## (Showing 8 of 23 unconstrained eigenvalues)
plot(vare.ca)

chisq.test(varespec/sum(varespec))
## Warning in chisq.test(varespec/sum(varespec)): Chi-squared approximation
## may be incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  varespec/sum(varespec)
## X-squared = 2.0832, df = 989, p-value = 1
plot(vare.ca, scaling = 1)

## Detrended correspondence analysis
vare.dca <- decorana(varespec)
vare.dca
## 
## Call:
## decorana(veg = varespec) 
## 
## Detrended correspondence analysis with 26 segments.
## Rescaling of axes with 4 iterations.
## 
##                   DCA1   DCA2    DCA3    DCA4
## Eigenvalues     0.5235 0.3253 0.20010 0.19176
## Decorana values 0.5249 0.1572 0.09669 0.06075
## Axis lengths    2.8161 2.2054 1.54650 1.64864
plot(vare.dca, display="sites")

## Ordination graphics
data(BCI)
mod <- decorana(BCI)
plot(mod)

names(BCI)[1:5]
## [1] "Abarema.macradenia"    "Vachellia.melanoceras" "Acalypha.diversifolia"
## [4] "Acalypha.macrostachya" "Adelia.triloba"
shnam <- make.cepnames(names(BCI))
shnam[1:5]
## [1] "Abarmacr" "Vachmela" "Acaldive" "Acalmacr" "Adeltril"
pl <- plot(mod, dis="sp")
identify(pl, "sp", labels=shnam)

## integer(0)
stems <- colSums(BCI)
plot(mod, dis="sp", type="n")
sel <- orditorp(mod, dis="sp", lab=shnam, priority=stems, pcol = "gray", pch="+")

plot(mod, dis="sp", type="n")
ordilabel(mod, dis="sp", lab=shnam, priority = stems)

sel[1:14]
## Abarmacr Vachmela Acaldive Acalmacr Adeltril Aegipana Alchcost Alchlati 
##    FALSE    FALSE    FALSE    FALSE    FALSE    FALSE    FALSE    FALSE 
## Alibedul Allopsil Alseblac Amaicory Anacexce Andiiner 
##    FALSE    FALSE    FALSE     TRUE    FALSE    FALSE

Environmental interpretation

## Vector fitting
data(varechem)
ef <- envfit(vare.mds, varechem, permu = 999)
ef
## 
## ***VECTORS
## 
##             NMDS1    NMDS2     r2 Pr(>r)    
## N        -0.05737 -0.99835 0.2536  0.059 .  
## P         0.61977  0.78478 0.1939  0.095 .  
## K         0.76652  0.64222 0.1809  0.113    
## Ca        0.68526  0.72830 0.4119  0.006 ** 
## Mg        0.63258  0.77450 0.4270  0.002 ** 
## S         0.19145  0.98150 0.1752  0.135    
## Al       -0.87156  0.49029 0.5269  0.001 ***
## Fe       -0.93596  0.35211 0.4450  0.002 ** 
## Mn        0.79872 -0.60171 0.5231  0.001 ***
## Zn        0.61757  0.78652 0.1879  0.106    
## Mo       -0.90311  0.42940 0.0609  0.503    
## Baresoil  0.92481 -0.38043 0.2508  0.048 *  
## Humdepth  0.93280 -0.36041 0.5201  0.002 ** 
## pH       -0.64795  0.76168 0.2308  0.055 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
plot(vare.mds, display = "sites")
plot(ef, p.max = 0.1)

ef <- envfit(vare.mds ~ Al + Ca, varechem)
plot(vare.mds, display = "sites")
plot(ef)
tmp <- with(varechem, ordisurf(vare.mds, Al, add = TRUE))
with(varechem, ordisurf(vare.mds, Ca, add = TRUE, col = "green4"))

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
## 
## Estimated degrees of freedom:
## 4.72  total = 5.72 
## 
## REML score: 156.6558
## Factors
data(dune)
data(dune.env)
dune.ca <- cca(dune)
ef <- envfit(dune.ca, dune.env, permutations = 999)
ef
## 
## ***VECTORS
## 
##         CA1      CA2     r2 Pr(>r)  
## A1 0.998160 0.060614 0.3104  0.048 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                  CA1     CA2
## Moisture1    -0.7484 -0.1423
## Moisture2    -0.4652 -0.2156
## Moisture4     0.1827 -0.7315
## Moisture5     1.1143  0.5708
## ManagementBF -0.7258 -0.1413
## ManagementHF -0.3867 -0.2960
## ManagementNM  0.6517  1.4405
## ManagementSF  0.3376 -0.6761
## UseHayfield  -0.2861  0.6488
## UseHaypastu  -0.0735 -0.5602
## UsePasture    0.5163  0.0508
## Manure0       0.6517  1.4405
## Manure1      -0.4639 -0.1738
## Manure2      -0.5872 -0.3600
## Manure3       0.5187 -0.3172
## Manure4      -0.2059 -0.8775
## 
## Goodness of fit:
##                r2 Pr(>r)   
## Moisture   0.4113  0.006 **
## Management 0.4441  0.003 **
## Use        0.1845  0.090 . 
## Manure     0.4552  0.008 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
plot(dune.ca, display = "sites")
plot(ef)

plot(dune.ca, display = "sites", type = "p")
with(dune.env, ordiellipse(dune.ca, Management, kind = "se", conf = 0.95))
with(dune.env, ordispider(dune.ca, Management, col = "blue", label= TRUE))
with(dune.env, ordihull(dune.ca, Management, col="blue", lty=2))

Constrained ordination

##Model specification
vare.cca <- cca(varespec ~ Al + P + K, varechem)
vare.cca
## Call: cca(formula = varespec ~ Al + P + K, data = varechem)
## 
##               Inertia Proportion Rank
## Total          2.0832     1.0000     
## Constrained    0.6441     0.3092    3
## Unconstrained  1.4391     0.6908   20
## Inertia is scaled Chi-square 
## 
## Eigenvalues for constrained axes:
##   CCA1   CCA2   CCA3 
## 0.3616 0.1700 0.1126 
## 
## Eigenvalues for unconstrained axes:
##    CA1    CA2    CA3    CA4    CA5    CA6    CA7    CA8 
## 0.3500 0.2201 0.1851 0.1551 0.1351 0.1003 0.0773 0.0537 
## (Showing 8 of 20 unconstrained eigenvalues)
plot(vare.cca)

ordiplot3d(vare.cca, type = "h")

dune.cca <- cca(dune ~ Management, dune.env)
plot(dune.cca)

dune.cca
## Call: cca(formula = dune ~ Management, data = dune.env)
## 
##               Inertia Proportion Rank
## Total          2.1153     1.0000     
## Constrained    0.6038     0.2855    3
## Unconstrained  1.5114     0.7145   16
## Inertia is scaled Chi-square 
## 
## Eigenvalues for constrained axes:
##   CCA1   CCA2   CCA3 
## 0.3186 0.1825 0.1027 
## 
## Eigenvalues for unconstrained axes:
##    CA1    CA2    CA3    CA4    CA5    CA6    CA7    CA8    CA9   CA10 
## 0.4474 0.2030 0.1630 0.1346 0.1294 0.0949 0.0790 0.0653 0.0500 0.0432 
##   CA11   CA12   CA13   CA14   CA15   CA16 
## 0.0387 0.0239 0.0177 0.0092 0.0080 0.0042
vare.cca <- cca(dune ~ Moisture, dune.env)
plot(vare.cca)

##  Permutation tests
anova(vare.cca)
## Permutation test for cca under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = dune ~ Moisture, data = dune.env)
##          Df ChiSquare      F Pr(>F)    
## Model     3   0.62831 2.2536  0.001 ***
## Residual 16   1.48695                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod <- cca(varespec ~ Al + P + K, varechem)
anova(mod, by = "term", step=200)
## Permutation test for cca under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = varespec ~ Al + P + K, data = varechem)
##          Df ChiSquare      F Pr(>F)    
## Al        1   0.29817 4.1440  0.001 ***
## P         1   0.18991 2.6393  0.005 ** 
## K         1   0.15605 2.1688  0.021 *  
## Residual 20   1.43906                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(mod, by = "margin", perm=500)
## Permutation test for cca under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = varespec ~ Al + P + K, data = varechem)
##          Df ChiSquare      F Pr(>F)    
## Al        1   0.31184 4.3340  0.001 ***
## P         1   0.16810 2.3362  0.017 *  
## K         1   0.15605 2.1688  0.023 *  
## Residual 20   1.43906                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(mod, by="axis", perm=1000)
## Permutation test for cca under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = varespec ~ Al + P + K, data = varechem)
##          Df ChiSquare      F Pr(>F)    
## CCA1      1   0.36156 5.0249  0.001 ***
## CCA2      1   0.16996 2.3621  0.039 *  
## CCA3      1   0.11262 1.5651  0.129    
## Residual 20   1.43906                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Session information

sessionInfo()
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 7 x64 (build 7601) Service Pack 1
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_Australia.1252  LC_CTYPE=English_Australia.1252   
## [3] LC_MONETARY=English_Australia.1252 LC_NUMERIC=C                      
## [5] LC_TIME=English_Australia.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] vegan3d_1.1-2   MASS_7.3-51.1   vegan_2.5-3     lattice_0.20-38
## [5] permute_0.9-4  
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.0              knitr_1.21             
##  [3] cluster_2.0.7-1         magrittr_1.5           
##  [5] scatterplot3d_0.3-41    xtable_1.8-3           
##  [7] R6_2.3.0                stringr_1.3.1          
##  [9] tools_3.5.1             webshot_0.5.1          
## [11] parallel_3.5.1          grid_3.5.1             
## [13] nlme_3.1-137            mgcv_1.8-26            
## [15] xfun_0.4                miniUI_0.1.1.1         
## [17] htmltools_0.3.6         crosstalk_1.0.0        
## [19] yaml_2.2.0              digest_0.6.18          
## [21] rgl_0.99.16             Matrix_1.2-15          
## [23] shiny_1.2.0             later_0.7.5            
## [25] htmlwidgets_1.3         promises_1.0.1         
## [27] manipulateWidget_0.10.0 mime_0.6               
## [29] evaluate_0.12           rmarkdown_1.11         
## [31] stringi_1.2.4           compiler_3.5.1         
## [33] jsonlite_1.6            httpuv_1.4.5.1