vegan_tutorial.Rmd

---
title: "Write-up of *vegan* Tutorial"
author: "Heili Lowman"
date: "6/4/2020"
output: html_document
---

# Vegan tutorial adapted from An Bui (EEMB, UCSB)

This is an R Markdown document in which I work through the *vegan* package tutorial for use in analyzing community data.

First, I will load the necessary packages and data:

```{r setup}

# Load libraries.
library(tidyverse) # Loads the tidyverse, including ggplot2.
library(forcats) # Loading these until I can get the tidyverse working...
library(vegan) # Loads the vegan package, which we will be exploring.
library(ggvegan) # Loads package to facilitate plotting from vegan in ggplot2.

# Load data.

birds <- read.csv("bta_snm_vbirds.csv") %>% # Bird community data.
  rename(site = X) %>%
  column_to_rownames("site") # Makes site numbers into row names.
# This data is sites (row) and bird species (column); all data entries are counts.

env <- read.csv("bta_snm_vfield.csv") # Environmental variables.
# This data is sites (X), stem density or # per hectare (all_stem_dens),  large stem basal area (big_stem_bas), percent canopy cover (can_cov_meas), mean canopy height (can_heig_mea), and ecological landtype (ELT).

site_type <- env %>% 
  select(X, ELT) %>% 
  rename(site = X, landtype = ELT) # Sets up metadata dataframe with new columns names for use later.

site_type$sitec <- as.character(site_type$site) # Create character column, otherwise below will not work.

```

# specnumber()

## How speciose are my communities?

specnumber() will tell you the number of species within each sample. You can then run an analysis of variance to ask if mean species richness is significantly different across sites.

```{r specnum}

sppr <- specnumber(birds) 

#spprdf <- as.data.frame(sppr)%>%
#  rowid_to_column(var='name') # Adds in a new column of sites for later use.

# To run an anova on the data, use a similar format as other models : response ~ dependent, data = environmental grouping

sppr_aov <- aov(sppr ~ landtype, data = site_type)
summary(sppr_aov)

```

There is no significant difference in species richness between ecological landtypes (p = 0.06). But, let's examine a plot of the data in any case.

```{r specplot}

sppr_df <- sppr %>% 
  enframe() %>% 
  full_join(site_type, by = c("name" = "sitec")) # "name" refers to the hidden column of characters numbering each of the sppr dataset values.

sppr_df$landtypef <- factor(sppr_df$landtype, levels = c( "do", "dm", "wm")) # Adds in new column for landtype as a leveled factor.

pal <- c("lightsalmon1", "tan4", "springgreen3")

plot_sppr <- ggplot(sppr_df, aes(x = landtypef, y = value, fill = landtypef)) +
  geom_boxplot() +
  scale_fill_manual(values = pal) +
  scale_x_discrete(labels = c("dry \n (n = 96)", "mix \n (n = 59)", "riparian \n (n = 55)")) +
  theme(legend.position = "none",
        plot.background = element_rect("white"),
        panel.background = element_rect("white"),
        panel.grid = element_line("grey90"),
        axis.line = element_line("gray25"),
        axis.text = element_text(size = 12, color = "gray25"),
        axis.title = element_text(color = "gray25"),
        legend.text = element_text(size = 12)) + 
  labs(x = "Ecological landtype",
       y = "Number of species per site",
       title = "Species richness")

plot_sppr

```

*Figure 1. Species Richness across Ecological Landtypes.*

Although number of species did not vary significantly by site, the median (50th percentile) richness in mixed (dry-mesic) and riparian (wet-mesic) sites appeared higher than at dry (dry-arid) sites and the interquartile range (25th to 75th percentile) at mixed and riparian sites was roughly similar. Furthermore, data collected at dry sites displayed a smaller interquartile range with more outliers at particularly speciose sites. It is worth noting, however, that nearly double the samples were collected from dry sites than at mixed or riparian sites.

# diversity()

## How diverse are my communities?

diversity() may calculate Shannon, Simpson, and Fisher's alpha. Here, we will specify Shannon, which takes into account species abundance and evenness. You can then run an analysis of variance to ask if mean species diversity is significantly different by landtype.

```{r diversity}

shannondiv <- diversity(birds) # Calculate Shannon diversity, the default for the index input.

head(shannondiv) # Display the first row of data.

sppdiv_aov <- aov(shannondiv ~ landtype, data = site_type) # Run the anova by landtype.

summary(sppdiv_aov) # Display anova results.

```

Again, there is no significant difference in bird species diversity (Shannon index) by ecological landtype (p = 0.06). But let's plot this out.

```{r divplot}

shandiv_df <- shannondiv %>% 
  # put all those calculations into a data frame
  enframe() %>% 
  # rename columns for ease of joining
  rename(site = name,
         shan_div = value)

div_plot_df <- shandiv_df %>% 
  # join with site_type
  full_join(site_type, ., by = c("sitec" = "site")) %>% 
  # group by landtype
  group_by(landtype) %>% 
  # calculate mean and standard error of diversity
  summarize(mean = round(mean(shan_div), 2),
            err = sd(shan_div)/sqrt(length(shan_div))) %>% 
  dplyr::mutate(label = "mean") %>% 
  unite("mean_label", label, mean, sep = " = ", remove = FALSE)

div_plot_df$landtypef <- factor(div_plot_df$landtype, levels = c( "do", "dm", "wm")) # Adds in new column for landtype as a leveled factor.

clean_background <- theme(plot.background = element_rect("white"),
        panel.background = element_rect("white"),
        panel.grid = element_line("white"),
        axis.line = element_line("gray25"),
        axis.text = element_text(size = 12, color = "gray25"),
        axis.title = element_text(color = "gray25"),
        legend.text = element_text(size = 12),
        legend.key = element_rect("white"))

plot_shandiv <- ggplot(div_plot_df, aes(x = landtypef, y = mean, fill = landtypef)) +
  geom_col(color = "black") +
  scale_fill_manual(values = pal) +
  geom_errorbar(aes(ymin = mean - err, ymax = mean + err), width = 0.5) +
  geom_text(aes(x = landtype, y = mean + err + 0.07, label = mean_label)) +
  scale_x_discrete(labels = c("dry \n (n = 96)", "mix \n (n = 59)", "riparian \n (n = 55)")) +
  scale_y_continuous(limits = c(0, 2.75), expand = c(0,0)) +
  clean_background + 
  theme(legend.position = "none") +
  labs(x = "Ecological landtype",
       y = "Mean Shannon diversity",
       title = "Shannon diversity")

plot_shandiv

```

*Figure 2. Species Diversity across Ecological Landtypes.*

There was no significant difference in species diversity by ecological landtype, and based on the above figure, this makes sense. The mean Shannon diversity indices for all three landtypes are within 0.1 of one another.

# adonis()

## How different are my communities in species composition?

adonis() may be used to run a permANOVA (**per**mutational **M**ultivariate **An**alysis **o**f **Va**riance) to assess differences in community composition based on dissimilarity. Here, the null hypothesis is that the centroids of you groups are equivalent for all groups.

```{r adonis}

bird_perm <- adonis(birds ~ ELT, data = env)

bird_perm

```

These results suggest that bird communities vary significantly by ecological landtype (p = 0.001). In order to represent this data in ordination space, there are a number of possible methods, but we will use PCA (**P**rinciple **C**omponents **A**nalysis). We'll use the rda() function in the vegan package to accomplish this.

```{r rdaplot}

birdPCA <- rda(birds) # Runs a redundancy analysis on the original birds dataset.

birdPCA # Provides results of this analysis.

PCAscores <- scores(birdPCA, display = "sites") %>% 
  as.data.frame() %>% 
  rownames_to_column("site") %>% 
  full_join(site_type, by = c("site" = "sitec")) # Creates dataset of analysis results.

PCAvect <- scores(birdPCA, display = "species") %>% 
  as.data.frame() # Creates dataset of vector values for each species.

PCAscores$landtypef <- factor(PCAscores$landtype, levels = c( "do", "dm", "wm")) # Adds in new column for landtype as a leveled factor.

plot_PCA <- ggplot() +
  geom_point(data = PCAscores, aes(x = PC1, y = PC2, color = landtypef)) +
  scale_color_manual(values = pal) +
  geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
  geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +
  geom_segment(data = PCAvect, aes(x = 0, y = 0, xend = PC1, yend = PC2), arrow = arrow(length = unit(0.2, "cm"))) +
  geom_text(data = PCAvect, aes(x = PC1, y = PC2, label = rownames(PCAvect))) +
  clean_background +
  labs(x = "PC1 (23.57%)",
       y = "PC2 (12.23%)",
       title = "Principal Components Analysis") 

plot_PCA

```

*Figure 3. Community Composition by Ecological Landtype Using ggplot2.*

Each point in the above figure is a site, colored by landtype, and each arrow is a bird species. The length of the arrow indicates the amount of variation in your communities explained by that particular variable (longer = larger), and the angle of the arrows relative to one another indicates the amount they correlate with one another (more obtuse = less correlated).

Another option is to plot this same figure using the ggvegan package.

```{r autoplot}

PCA_biplot <- autoplot(birdPCA) # Creates the plot.

PCA_biplot # Output of same, but less customized, figure 3.

```


*Figure 4. Community Composition across All Sites Sampled.*

And we can use the fortify() function to better customize and replicate Figure 3 within the ggvegan package.

```{r fortify}

PCA_fortify <- fortify(birdPCA) # This provides results in a dataframe.

PCA_fort_sites <- PCA_fortify %>% 
  filter(Score == "sites") %>% 
  full_join(., site_type, by = c("Label" = "sitec")) # This extracts the results of the PCA only for the sampling sites and joins it with the dataframe containing landtypes.

PCA_fort_species <- PCA_fortify %>% 
  filter(Score == "species") # This extracts the results for the PCA only for the species (so the vectors).

PCA_fort_sites$landtypef <- factor(PCA_fort_sites$landtype, levels = c( "do", "dm", "wm")) # Adds in new column for landtype as a leveled factor.

PCA_fortify_plot <- ggplot() +
  geom_point(data = PCA_fort_sites, aes(x = PC1, y = PC2, col = landtypef)) +
  geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
  geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +
  scale_color_manual(values = c("lightsalmon1", "tan4", "springgreen3")) +
  geom_segment(data = PCA_fort_species, aes(x = 0, xend = PC1, y = 0, yend = PC2), 
               arrow = arrow(length = unit(0.2, "cm"))) +
  geom_text(data = PCA_fort_species, aes(x = PC1, y = PC2, label = Label)) +
  clean_background +
  labs(x = "PC1 (23.57%)",
       y = "PC2 (12.23%)",
       title = "Principal Components Analysis - using fortify()") # This creates the plot.

PCA_fortify_plot # And this prints the plot to mimic figure 3.

```

*Figure 5. Community Composition by Ecological Landtype Using ggvegan.*

Now, yet another technique you may use is NMDS (*N*on-metric *M*ulti*d*imensional *S*caling), which allows you to collapse all species axes into 2 to plot in cartesian space and visualize the differences between samples and sites using metaMDS().

```{r metamds}

bird_NMDS <- metaMDS(birds) # Runs the NMDS analysis.

bird_NMDS # Presents summary of NMDS results.

stressplot(bird_NMDS) # For evaluating how well ordination represented data complexity.

```

*Figure 6. Stressplot of NMDS Results.*

The above figure shows how closely ordination along the y-axis represents dissimilarities along the  x-axis. Each point is a community, and their distance from the line represents the "stress," or how they are pulled from their original position to be represented in their ordination.

**Rule of thumb: High stress is bad, low stress is good.** NMDS stress will increase as sample size increases, but *generally* speaking values above 0.2 aren't great and suggest that the ordination didn't do a god job of representing community structure. (< 0.05 is excellent, < 0.1 is great, < 0.2 is good.) See An's original tutorial (https://rpubs.com/an-bui/vegan-cheat-sheet) for more exploration of this.

```{r nmdsplot}

plot(bird_NMDS) # Creates base R plot of results.

```

*Figure 7. Base R Plot of NMDS Results.*

```{r nmdsggplot}

plot_df <- scores(bird_NMDS, display = "sites") %>% 
  as.data.frame() %>% 
  rownames_to_column("site") %>% 
  full_join(site_type, by = c("site" = "sitec")) # Creates extracted dataset from NMDS analysis.

plot_df$landtypef <- factor(plot_df$landtype, levels = c( "do", "dm", "wm")) # Adds in new column for landtype as a leveled factor.

plot_nmds <- ggplot(plot_df, aes(x = NMDS1, y = NMDS2, color = landtypef, shape = landtypef)) +
  geom_point(size = 3, alpha = 0.8) +
  scale_color_manual(values = pal) +
  stat_ellipse(linetype = 2, size = 1) +
  clean_background +
  labs(title = "NMDS") # Creates plot with ellipses to denote 95% confidence intervals fitted into the spatial ordination.

plot_nmds

```

*Figure 8. Ggplot of NMDS Results.*

You can also create this same plot with vectors for each species added in to see how they contribute to relative differences in community composition.

```{r envfitplot}

fit <- envfit(bird_NMDS, birds, perm = 999) # Creates new dataset of NMDS results + species matrix.

fit_pvals <- fit$vectors$pvals %>% 
  as.data.frame() %>% 
  rownames_to_column("species") %>% 
  dplyr::rename("pvals" = ".") # Extracts p-value for each species.

fit_spp <- fit %>% 
  scores(., display = "vectors") %>% 
  as.data.frame() %>% 
  rownames_to_column("species") %>% 
  full_join(., fit_pvals, by = "species") %>% 
  filter(pvals == 0.001) # Extracts vector coordinates for each species, but then only keeps in species with p = 0.001.

nmds_plot_new <- ggplot(plot_df, aes(x = NMDS1, y = NMDS2)) +
  coord_fixed() +
  geom_point(aes(color = landtypef, shape = landtypef), size = 3, alpha = 0.8) +
  stat_ellipse(aes(color = landtype)) +
  scale_color_manual(values = pal) +
  geom_segment(data = fit_spp, aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2),
               arrow = arrow(length = unit(0.25, "cm")),
               col = "black") +
  geom_text(data = fit_spp, aes(label = species)) +
  clean_background # Creates new biplot.

nmds_plot_new

```

*Figure 9. NMDS Results with Species Vectors.*

# cca()

## How is community structure related to specific environmental variables?

All of the above examples are *unconstrained* ordination, which means that the ordination is done only considering species counts in the site * species matrix. However, *constrained* ordination is appropriate for when there is a matrix of communities (site * species) and another matrix of environmental characteristics (site * environment), which in this case is the env dataframe we've already loaded in.

*Constrained* ordination allows us to ask how environmental variables relate to community composition - here we will use CCA (**C**anonical **C**orresnpondence **A**nalysis) to examine the effects of canopy height, stems per hectare, and big stem diameter on bird community structure.

```{r cca}

birdCCA <- cca(birds ~ can_heig_mea + all_stem_den + big_stem_bas, data = env) # Runs the CCA analysis.

birdCCA # Reports the results.

ccaplot <- plot(birdCCA) # Plots the results in base R.

```

*Figure 10. Base R Plot of CCA Results.*

And we can again extract this data to create a more customized plot in ggplot.

```{r ggcca}

ccavectors <- as.matrix(scores(birdCCA, display = "bp", scaling = "species")*12.2) %>% 
  as.data.frame() # Pulls out environmental factors as vector data.

site_data <- scores(birdCCA, display = "sites") %>% 
  as.data.frame() %>% 
  rownames_to_column("site") %>% 
  full_join(site_type, by = c("site" = "sitec")) # Joins together CCA results with environmental dataframe.

site_data$landtypef <- factor(site_data$landtype, levels = c( "do", "dm", "wm")) # Adds in new column for landtype as a leveled factor.

species_data <- scores(birdCCA, display = "species") %>% 
  as.data.frame() # Pulls out bird community data by site.

plot_cca <- ggplot(site_data) +
  geom_point(aes(x = CCA1, y = CCA2, color = landtypef), shape = 19, size = 2, alpha = 0.8) +
  coord_fixed() +
  scale_color_manual(values = pal) +
  geom_segment(data = ccavectors, aes(x = 0, y = 0, xend = CCA1, yend = CCA2), arrow = arrow(length = unit(0.2, "cm"))) +
  geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
  geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +
  geom_point(data = species_data, aes(x = CCA1, y = CCA2), shape = 17, size = 2, color = "slateblue") +
  scale_x_continuous(limits = c(-12, 10)) +
  scale_y_continuous(limits = c(-3, 12)) +
  geom_text(data = ccavectors, aes(x = CCA1, y = CCA2, label = rownames(ccavectors)), nudge_x = 0.3, nudge_y = 0.3) +
  clean_background +
  labs(title = "Canonical Correspondence Analysis")

plot_cca

```

*Figure 11. Ggplot of CCA Results.*

NEED TO INSERT FINAL INTERPRETATION OF FIGURE 11 HERE.

End of R Markdown script.