Search code examples
rggplot2veganmds

plot vectors on MDS using ggplot

I need to plot a MDS with vectors showing changes in abundance for my species. I need the plot with just vectors

This is my data with abundance for each species and code

library(vegan)

A <- c(54.67, 37.67, 19.33, 0, 6, 8, 84.67, 0,0,0,0,0,0,0)
B <- c(3.67, 10.33, 32.67, 5.33, 20.33, 5.33, 4.67, 3, 4, 0.01, 0.1, 0, 5, 0)
C <- c(10, 1.67, 2.67, 1.67, 11.33, 1.33, 1, 2, 2.77, 0, 0.02, 1,3,0)
D <- c(1,10.33, 2.33, 28.33, 29.33, 4.33, 21, 6.97, 4.47, 0, 0.16, 11, 4,0)
df <- cbind(A, B, C, D)
row.names(df) <- c('B_2016', 'Emb_2016', 'Fes_2016', 'Ku_2016', 'Ra_2016', 'Ud_2016',
                   'Ve_2016', 'Ba_2017', 'Emb_2017', 'Fes_2017', 'Ku_2017', 'Ra_2017', 
                   'Ud_2017', 'Ve_2017')

mds <- metaMDS(df, distance='bray')

I am using these codes to create a dataframe 

mdspoints <- data.frame(scores(mds))
mdsvectors <- data.frame(mds$species)

and this is the code I am using to graph

g <- ggplot(data = mds, aes(MDS1, MDS2)) + 
  geom_segment(data = mdsvectors, aes(x=0, xend=MDS1, y=0, yend=MDS2),
               arrow = arrow(length = unit(0.5, "cm")),
               colour="grey", inherit_aes = FALSE) + 
  geom_text(data=mdspoints, aes(x=MDS1, y=MDS2, label=species), size=5)

but I cannot graph anything and get an error (Error: ggplot2 doesn't know how to deal with data of class metaMDS/monoMDS).

I would like something like this

Thank you

Solution

  • I'm not sure exactly what you are aiming for, based on your code, but here are some things you may want to take note.

    Point 1: Don't put anything in the top level ggplot() call, unless you want the subsequent layers to inherit it.

    Instead of:

    g <- ggplot(data = mds, aes(MDS1, MDS2)) +

    Use:

    g <- ggplot() +

    You've already created data frames mdspoints & mdsvectors, and none of your geom layers requires anything from mds. You really don't need it here. But since it is there, ggplot will check it.

    Had mds been a data frame, it would pass ggplot's check, & be ignored since the subsequent layers don't need it. However, it is metaMDS / monoMDS object, causing ggplot to throw the error you saw.

    Point 2: Check if your data frames are what you expected.

    Your code includes the following line:

      geom_text(data=mdspoints, aes(x=MDS1, y=MDS2, label=species), size=5)

    This tells ggplot that for plotting labels, it should look at the mdspoints data frame, and search for variables named MDS1 / MDS2 / species.

    This is what was actually created from mdspoints <- data.frame(scores(mds)):

    > mdspoints
             NMDS1         NMDS2
B_2016   -141.6526 -6.290613e-01
Emb_2016 -141.8424 -3.280861e-01
Fes_2016 -142.1144 -4.456856e-01
Ku_2016  -141.8335  3.674413e-01
Ra_2016  -141.8977  2.283486e-02
Ud_2016  -141.8824 -1.480702e-01
Ve_2016  -141.5302 -3.732303e-01
Ba_2017  -141.9265  2.233302e-01
Emb_2017 -141.9695  1.210940e-01
Fes_2017 -140.6462  1.430899e-01
Ku_2017  -141.8616  2.216499e-01
Ra_2017  -141.7638  7.116520e-01
Ud_2017  -142.0109  1.130730e-01
Ve_2017  1842.9317 -3.167902e-05

    So, NMDS1 / NMDS2 instead of MDS1 / MDS2, and no column name for "species". Are the row names corresponding to species? I'm not sure as I don't use the vegan package myself, but a quick look at the help file for its scores() function uncovers the following:

    ## Default S3 method:
scores(x, choices, display=c("sites", "species"), ...)

    Which suggests that this may be the scores for sites, rather than species. If that understanding is correct, you'll want to specify "species" when you create mdspoints, and manually create the species column from row names:

    mdspoints <- data.frame(scores(mds, "species"))
mdspoints$species <- row.names(mdspoints)

    Result

    Here's what the plot could look like:

    ggplot() + 
  geom_segment(data = mdsvectors, aes(x=0, xend=MDS1, y=0, yend=MDS2),
               arrow = arrow(length = unit(0.5, "cm")),
               colour="grey") +
  geom_text(data=mdspoints, aes(x=NMDS1, y=NMDS2, label=species), size=5) +
  labs(x = "NMDS1", y = "NMDS2") + # add axis labels
  theme_classic()                  # use a white theme for better contrast

    plot