CODEWITHSUNDEEP
rvegan
i.b
i.b

Reputation: 337

how to plot only most abundant species in NMDS?

I need to plot an ordination plot showing only let s say the 20 most abundant species.

I tried to do the sum of the species colunm and then select only a certain sum value:

abu <- colSums(dune)
abu
sol <- metaMDS(dune)
sol
plot(sol, type="text", display="species", select = abu > 40)

I get this error: select is not a graphical parameter

I would expect to see only small number of species but it does not happen, how do you show only a small number of species in the NMDS plot?

Upvotes: 1

Views: 2157

Answers (1)

January
January

Reputation: 17090

This is not straightforward. You are getting an error because select is not a parameter for the plot. Unfortunately, the result of the analysis is not a data.frame that could be handled easily (e.g. with tidyverse), and even more unfortunately, the plot() function called is not your standard plot, but a method defined specifically for objects of this class. The authors of this method did not foresee your need, and therefore, we must make the plot manually. But to do that, we need to understand what is plotting and how.

Let us find out more about the object sol:

class(sol)
# [1] "metaMDS" "monoMDS"
methods(class="metaMDS")
# [1] goodness    nobs        plot        points      print       scores      sppscores<- text

Oh good, we have a plot method. After a moment of digging, we find it in the vegan package (not exported, so we need to access it via vegan:::plot.metaMDS). It appears to be a wrapper around a function called ordiplot. We edit the function with edit() to figure out what it is doing. Essentially, it boils down to the following (with loads of unnecessary code):

Y <- scores(sol, display="species")
plot(Y, type="n")
text(Y[,1], Y[,2], rownames(Y), col="red")

This is, more or less, your plot. Choosing the species to show is now trivial, but first we must make sure that rows of Y are in the same order as columns of dune:

all(colnames(dune) == rownames(Y))
Y.sel <- Y[colSums(dune) > 40, ]
plot(Y.sel[,1], Y.sel[,2], type="n", xlim=c(-.8, .8), ylim=c(-.4, .4))
text(Y.sel[,1], Y.sel[,2], rownames(Y.sel), col="red")

enter image description here

We can of course make a much nicer plot. For example, with ggplot (it is definitely possible to make a much nicer plot with base R as well). We could actually show the abundance of the plants using the size esthetics:

library(ggplot2)
library(ggrepel)

Y <- data.frame(Y)
Y$abundance <- colSums(dune)
Y$labels <- rownames(Y)

ggplot(Y, aes(x=NMDS1, y=NMDS2, size=abundance)) +
   geom_point() + geom_text_repel(aes(label=labels)) + 
   theme_minimal()

enter image description here

To filter the species by abundance, we now can do the following:

library(tidyverse)
Y %>% filter(abundance > 40) %>%
  ggplot(Y, aes(x=NMDS1, y=NMDS2, size=abundance)) +
    geom_point() + geom_text_repel(aes(label=labels)) + 
    theme_minimal()

Upvotes: 1

Related Questions