How to color code different sites in R using NMDS?

Question

I am doing an independent study looking at species community composition data for fish in the San Antonio River basin. I am using NMDS in R to meta analyze my results.

library(vegan)
test2<-data.frame(relative_abundance_fish_data[1:100,5:51])
test2.hel<-decostand(test2,method="hellinger")
set.seed()
nmds1<-metaMDS(test2.hel,autotransform=F)
ordiplot(nmds1) 

Then, I have my data divided up where each row is an independent data point. For example,

Site Code Site Date Year relative abundances of species... env variables...

Each site code is unique, but I only have 12 total sites. I would like to colorcode my sites based off of the SITES, not the SITE CODE. Is there a way to do this?

Please let me know if I need to clarily anything else!

I have seen answer threads that colorcode based on individual events (n=100)- I just want to colorcode among sites (n=12). In that same thread, there were an equal number of sites for each category, but that is not the case for my dataset. For example, I have 8 events for one site, and maybe only 6 for another.

Answer 1

It is often simpler to use a data set included with R that is similar to your data. For example the vegan package contains a data set varespec that may be similar to your data (?varespec to get a description of the dat). It contains 24 rows and 44 variables. It does not have groups, but we can add them:

library(vegan)
data(varespec)
dim(varespec)
# [1] 24 44
grps <- factor(c(rep("A", 6), rep("B", 6), rep("C", 6), rep("D", 6)))
grps
#  [1] A A A A A A B B B B B B C C C C C C D D D D D D
Levels: A B C D

The groups do not have to be contiguous if you have names for them, just convert variable to a factor. Now run the analysis:

var.hel <- decostand(varespec, method="hellinger")
nmds1 <- metaMDS(var.hel, autrotransform=FALSE)

We need to create integer values for the groups, A = 1, B = 2, etc. Then we need a vector of colors (or symbols) to use in the plot:

cdx <- as.integer(grps)
clr <- c("red", "blue", "green", "violet")

Finally we create the plot and add the points:

op <-ordiplot(nmds1, display="sites", type="none")
points(op, "sites", pch=16, col=clr[cdx])

This may not be exactly what you want, but it should provide a guide you can adapt to your data. Reading the manual page for ?ordiplot will help.