I am sure there is an easy answer for this. I have been using this code for years and the other day at work I went to make an ordination, and it stopped working. :( Granted, it has been at least a year since I ran this particular code so maybe the package has changed with an update. Also, there is probably a better way to do this, but this is code I got from someone and it has always worked.
I will give the example with an older data set of mine, since I am running into the same issue with this old data set as the one I have at work. Pretty basic, ecological data with otu frequencies and site characteristics (this is a small look, there are three treatments and three species, 46 otus). This is similar to the dune or other included data sets in R:
For the ordination, I want the shapes to be the different plant species and the colors to be the different treatments. Here is the code I have always used in the past and had it work ever time:
data1<-read.csv("treatment_by_species.csv")
#data frame with fungal otus
fungal_measure<-data.frame(data1[,3:48],row.names=NULL)
#Lables for graph
trtlabels<-data1[,c(1,1)]
splabels<-data1[,c(1,2)]
#distance matrix
dm=bcdist(fungal_measure)
meta_out=metaMDS(dm,k=2,trymax=9999,utotransform=F)
meta_out
windows(title="meta_out")
stressplot(meta_out)
#graph
color1=c("black","grey","grey96")
windows()
plot(meta_out$points,type="n",xlim=c(-0.6,0.6),ylim=c(-0.6,0.6))
points(meta_out$points,cex=1.5,pch=c(21,22,24)[splabels$Species],
bg=color1[trtlabels$Treatment])
legend("topleft",bty="n",pch=c(21,22,24,15,15,15),
col=c("black","black","black","black","grey","grey96"),
legend=c("Aspen","Pine","Spruce","FFM","Peat","Subsoil"))
So what happens is the blank graph will come up just fine, but the points won't graph at all and there is no error message. I am guessing it has something to do with how I am coding the labels. If someone could point to a post or easier way to do this, it would be most appreciated. Always a bummer when things stop working.
