How to plot sites in different colours in a NMDS plot in R (vegan package)?

Question

I have data from grassland (sites 1-14) and quarry (sites 15-28) sites and want them to have different colouring in my NMDS plot that I created with the vegan package, e. g. site numbers 1-14 in red and 15-28 in blue.

I'm new to R and not familiar with editing of graphics.

vare.mds <- metaMDS(species, trace=FALSE)
plot(vare.mds, type = "n")
text(vare.mds, display = "sites")

Answer 1

Assuming that the samples are in the order you indicate (rows 1–14 are grassland, rows 15–28 are quarry) then the simplest way to do this in base graphics is to create a factor variable of the required length:

grp <- factor(rep(c('grassland', 'quarry'), each = 14))

Then create a vector of colours for the two groups

cols <- c('red','black')

Then you can colour the text created by the text() call by using the col argument, which is typically available for most base graphic plotting functions.

To see what will happen when you actually do the plot, what we are going to do is index the colour vector cols using grp the factor containing the group membership:

> cols[grp]
 [1] "red"   "red"   "red"   "red"   "red"   "red"   "red"   "red"   "red"  
[10] "red"   "red"   "red"   "red"   "red"   "black" "black" "black" "black"
[19] "black" "black" "black" "black" "black" "black" "black" "black" "black"
[28] "black"

grp is stored as a numeric vector whose elements index the level for each observation:

> as.numeric(grp)
 [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2

As a result, in cols[grp], grp becomes a numeric vector indexing the first or second colour respectively and as needed for each of your sites.

So, the text() will look like this:

text(vare.mds, display = "sites", col = cols[grp])

In you case, it would also be sufficient to simply call

text(vare.mds, display = "sites", col = rep(cols, each = 14))

because the groups are arranged in contiguous blocks of samples, but if this is not the case (and resorting is not necessarily a good thing if you forget to sort the same way associated data sets) then the method I describe above is beneficial.

To make a directly reproducible example, here's is the full thing worked through for the varespec data set provided with vegan

library('vegan')
data(varespec)

## varespec has on 24 observations so 1-12 will be grassland and 13-24 quarry
grp <- factor(rep(c('grassland', 'quarry'), each = 12))

## vector of colours
cols <- c('red', 'black')

set.seed(1)
ord <- metaMDS(varespec, trace=FALSE)
plot(ord, type = 'n')
text(ord, display = 'sites', col = cols[grp])
legend('bottomright', legend = tools::toTitleCase(levels(grp)),
       fill = cols, bty = 'n')

which produces