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 ----

   # 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")

You are viewing proxied material from republic.circumlunar.space. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.