TITLE: An email about ordination and environmental fits
DATE: 2020-01-15
AUTHOR: John L. Godlee
====================================================================


I got an email from a colleague asking for my opinion on analysing
the environmental determinants of tree species diversity measured in
multiple plots across a province in their country. I did a load of
reading on this sort of thing about 6 months ago so I relied heavily
on my notes to write a response and also provide an R script with my
recommendations.

The original email, paraphrased slightly, said:

 Dear John

 I need your help with diversity analysis in R for my masters
 student. Our goal is to detect relations between environmental
 variables and diversity.

 All variables are numeric, we would like to know what kind of
 analysis you suggest?

 Ordinations? Regressions? PCA? NMDS?

 Can you guide Us?

They also sent a document with their overall research goal:

 Evaluate and relate by statistical analysis the effect of soil,
 climate and topography on the spatial distribution of diversity
 (Richness index Shannon index and Simpson index ) of vegetation in
 province.

y email response was:

 My preference for an ordination technique to analyse variation in
 species composition is NMDS. This is for two main reasons. The
 first is that NMDS deals well with instances where there are no
 shared species between two sites, which I imagine may happen a lot
 in your very large dataset, across a wide geographical area.
 Because Principle Co-ordinate Analysis (PCoA) is an eigenvalue
 analysis and must find a unique solution, when two or more sites
 share no common data, this unique solution cannot be found. The
 second reason is that in contrast to Principle Component Analysis
 (PCA), NMDS can properly represent variation in species abundances
 because it is not limited to using a euclidean distance method. I
 would recommend using the Bray-Curtis distance method.

 I don’t think it is wise to use all the environmental variables in
 the NMDS, especially because as you have shown, some of the
 variables are highly correlated. I think you should choose
 variables which you think have a valid ecological reason for
 influencing the species composition.

 In addition to investigating how species composition (species
 identity and relative abundance) varies with environmental
 variables, you may wish to investigate how species diversity
 (species richness and abundance evenness) varies with
 environmental variables. To do this, I would use a simple general
 linear model with multiple environmental variables as the
 predictor variables and the Shannon index as the response
 variable.

 I don’t have your dataset so I have provided in an attached R
 script the way I would perform this sort of analysis in R, using
 an example dataset provided in the {vegan} R package which I hope
 is similar to your data.

The R script I sent looked like this:

   # Environmental determinants of diversity
   # John Godlee ([email protected])
   # 2020-01-12

   # Preamble ----

   # Remove old crap
   rm(list = ls())
   # dev.off()

   # Set working directory

   # Packages
   library(vegan)  # diversity(), metaMDS(), data(dune)
   library(ggplot2)  # ggplot()
   library(dplyr)  # %>%, select()
   library(tidyr)  # gather()

   # Import data
   data(dune)
   ##' Site (rows) by species (columns) abundance matrix
   ##' 20 sites
   ##' 30 species

   env <- read.csv("env.csv")
   ##' Dataframe with environmental variables for each of the 20 sites
   ##' This is fake data I created for this example

   # Define a function to estimate the optimal number of dimensions for an NMDS
   NMDS.scree <- function(x) {
     plot(rep(1, 10), replicate(10, metaMDS(x, autotransform = F, k = 1)$stress),
       xlim = c(1, 10),ylim = c(0, 0.30),
       xlab = "# of Dimensions", ylab = "Stress", main = "NMDS stress plot")
     for (i in 1:10) {
       points(rep(i + 1,10),
         replicate(10, metaMDS(x, autotransform = F, k = i + 1)$stress)
       )
     }
   }

   # Run NMDS.scree
   NMDS.scree(dune)
   ##' The scree plot shows that 3 dimensions produces a stress value below 0.1,
   ##' which is indicates that with this many dimensions, the NMDS provides a good
   ##' representation of the difference in species composition.

   # Run NMDS with 3 dimensions
   dune_nmds <- metaMDS(dune, distance = "bray", try = 500, trymax = 500, k = 3,
     autotransform = FALSE)

   # Assess the fit of the NMDS with a stressplot (Shepard plot)
   stressplot(dune_nmds)
   ##' The stressplot shows a strong positive correlation
   ##' It plots the distances among objects in the ordination plot against the
   ##' original Bray-Curtis distances. A tight positive correlation between
   ##' original distance and transformed distance shows that the dimensionality
   ##' reduction was successful.

   # Extract final stress value of the NMDS
   nmdsstress <- dune_nmds$stress
   ##' Should be quoted in results, to show validity of results of NMDS

   # Extract site (plot) scores from NMDS analysis
   plot_scores <- as.data.frame(scores(dune_nmds))

   # Extract species scores from NMDS analysis
   species_scores <- as.data.frame(scores(dune_nmds, "species"))
   species_scores$species_binomial <- rownames(species_scores)

   # Fit environmental variables to NMDS
   dune_envfit <- envfit(dune_nmds, env[,
     c("BIO1", "BIO12", "BIO15", "BIO17", "elevation", "sand_coarse")],
   permutations = 999)
   ##' Note that these aren't my recommendations for the variables you should
   ##' use on your data, they are just commonly used variables.

   # Which environmental variables significantly affect species composition?
   dune_envfit

   # Create ggplot of NMDS output

   # Get arrow vectors
   dune_envfit_arrows <- data.frame(dune_envfit$vectors$arrows)
   dune_envfit_arrows$var <- rownames(dune_envfit_arrows)
   dune_envfit_arrows$r2 <- dune_envfit$vectors$r
   dune_envfit_arrows$p <- dune_envfit$vectors$pvals

   dune_nmds_plot <- ggplot() +
     geom_hline(aes(yintercept = 0), linetype = 2) +
     geom_vline(aes(xintercept = 0), linetype = 2) +
     geom_point(data = species_scores,
       aes(x = NMDS1, y = NMDS2), shape = 1, colour = "black") +
     geom_point(data = plot_scores,
       aes(x = NMDS1, y = NMDS2), shape = 2, colour = "red") +
     geom_segment(data = dune_envfit_arrows,
       aes(xend = NMDS1*r2*2, yend = NMDS2*r2*2, x = 0, y = 0),
       arrow = arrow(length = unit(0.05, "npc")),
       colour = "blue") +
     geom_text(data = dune_envfit_arrows,
       aes(x = NMDS1*r2*2, y = NMDS2*r2*2, label = var),
       colour = "blue") +
       coord_equal()

   # General linear model of Shannon index vs. environment

   # Calculate Shannon index from abundance matrix
   env$shannon <- diversity(dune)

   # Run linear model with multiple predictors
   dune_lm <- lm(shannon ~ BIO1 + BIO12 + BIO15 + BIO17 + elevation + sand_coarse,
     data = env)
   ##' You may find that the model should be more complex than this, possibly with
   ##' interaction terms, depending on the variables you choose. Also, you may
   ##' want to introduce a spatial auto-correlation term. Finally, you should check
   ##' that the predictor variables conform to the assumptions of a general
   ##' linear model.

   # Which environmental variables significantly (p<0.05) affect species diversity?
   summary(dune_lm)

   # Plot of linear relationship between shannon and chosen environmental variables
   env_gather <- env %>%
     dplyr::select(plotcode, shannon, BIO1, BIO12, BIO15, BIO17, elevation, sand_coarse) %>%
     gather(key = "var", value = "value", -plotcode, -shannon)

   dune_lm_plot <- ggplot(env_gather, aes(x = value, y = shannon)) +
     geom_point() +
     geom_smooth(method = "lm") +
     facet_wrap(~var, scales = "free_x")