I'm plotting a series of RDAs in R, each of which has 10+ environmental vectors. The environmental variables each fall in to one of 5 categories. I'd like the vector colors to reflect these categories. I did it in a somewhat ghetto fashion by producing the raw black and white plot, then tracing it in Powerpoint (super time consuming), but there are ALOT of inherent issues with this but it looks good at least. Here's that plot, note the vector colors.
I'd like to have my output be an R native plot. The plot I use to make the base plot is currently:
#download the data#
env<- read.csv("environmenta data.csv")
bio<- read.csv("bio data.csv")
#break out the groups of environmental vectors#
#these are purley for example#
group1<-as.matrix(subset(env,select=c(a,b,c))
group2<-as.matrix(subset(env,select=c(d,e,f,h))
group3<-as.matrix(subset(env,select=c(g))
group4<-as.matrix(subset(env,select=c(j,k,l,o))
group5<-as.matrix(subset(env,select=c(m,n,))
#run the RDA#
rda1<-rda(bio,env)
#Plot it#
plot(rda1,type="n",bty="n",main="",
xlab="XX% variance explained",
ylab="XX% variance explained",
col.main="black",col.lab="black", col.axis="white",
xaxt="n",yaxt="n")
#points(rda1,display="species",col="gray",pch=20)
#option to display species points
#text(rda2,display="species",col="gray")
#option for species labels
points(rda1,display="cn",col="black",lwd=2)#<<< guessing this is the key statement for my issue, but not sure
text(rda1,display="cn",col="black",cex=0.5)
I'm assuming the answer falls in this points--> col="..." command, but I'm not sure how to tell it to subset by groups. Any and all help would be appreciated.
There are a couple of issues here:
plot.cca
The latter doesn't really matter so much if you are only plotting the biplot scores, but a general solution to this kind of plotting requires that we consider that other scores may also be required to be draw.
Here's a reproducible example using data built into
library('vegan')
data(varespec, varechem)
Now we fit a silly ordination (don't do this at home folks)
ord <- rda(varespec ~ ., data = varechem)
Next, extract the biplot scores
bp <- scores(ord, display = 'bp')
and assume we have a variable that contains the groups we want to assign each biplot score to
f <- factor(sample(1:5, nrow(bp), replace = TRUE))
and the associated vector of colours we want to use
cols <- c('red','black','green','navy','purple')
and expand this vector of colours into one of length nrow(bp)
, i.e. one colour per biplot score/variable
cols <- cols[f]
Next start to plot, by preparing an empty plot region into which we can add the arrows
plot(ord, type = 'n')
In the above example, we leave set up the plot region so that species and site/sample scores are accommodated.
We now calculate a multiplier that allows the biplot scores to occupy a specified proportion (fill = 0.75
) of the plot region
mul <- ordiArrowMul(bp, fill = 0.75)
Finally, add the biplot arrows, coloured as per their group
arrows(0, 0, mul * bp[,1], mul * bp[,2],
length = 0.05, col = cols)
and add labels to the biplot arrows
labs <- rownames(bp)
text(ordiArrowTextXY(mul * bp, labs), labs, col = cols)
This produces
library(vegan)
#DATA
data(varespec)
data(varechem)
env = varechem
bio = varespec
rm(list = c("varechem", "varespec"))
#Assign colors based on groups for the columns of env
groups = rep(c("red", "blue"), length.out = NCOL(env))
rda1 = rda(X = bio, Y = env)
plot(rda1, type = "n")
points(rda1)
text(rda1, display = "bp", col = groups)
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With