Yesterday, I tweeted a photo of a ordination I plotted with
ggvegan, and thought I should show how I made it.
Ordinations can be plotted with base R (see
?plot.cca). This is fine for simple plots, but it is a lot of effort to make a complex plot. The
ggvegan package, written by Gavin Simpson, lets you use the power of
ggplot2 and can make complex plots easier to make.
Start by loading some packages.
ggvegan will need to be installed from GitHub if you don’t already have it.
#load packages library(tidyverse) library(vegan) #devtools::install_github("gavinsimpson/ggvegan") library(ggvegan)
I’m going to use the
pyrifos dataset, which is included in
vegan. The data are the abundances of aquatic invertebrates from an experiment in which twelve ditches were studies several times before and after insecticide treatment. We need to make a
tibble (or a
data.frame) of the predictors (dose, time and ditch ID).
#load data data("pyrifos") env <- tibble( week = rep(c(-4, -1, 0.1, 1, 2, 4, 8, 12, 15, 19, 24), each = 12), dose = factor(rep(c(0.1, 0, 0, 0.9, 0, 44, 6, 0.1, 44, 0.9, 0, 6), 11)), ditch = gl(12, 1, length=132))
We can run a principal components analysis with
pca <- rda(pyrifos)
Simple ordination plots can be made with the
autoplot function. This plot shows all the samples. This is fine, but we want to see dose and time information.
autoplot(pca, layers = "sites", arrows = FALSE)
At the moment, it is difficult to do much to this plot. But we can use the
fortify function from
ggvegan to extract the scores from the ordination object and then combine these with the predictor data before plotting with
ggplot. I’m going to colour the samples by dose, and draw a line through the samples in each ditch, with a point on the first value.
pca_fort <- fortify(pca, display = "sites") %>% bind_cols(env) ggplot(pca_fort, aes(x = PC1, y = PC2, colour = dose, group = ditch)) + geom_path() + #use geom_path not geom_line geom_point(aes(size = if_else(week == min(week), 1, NA_real_)), show.legend = FALSE) + scale_color_viridis_d() + scale_size(range = 2) + coord_equal()
Don’t forget to add
coord_equal() to ensure that the plot is correctly scaled.