Add arrows in RDA in R

Question

I am relatively new to R and I am trying to get my head around how to do ordination techniques in R, so that I don't need to use other software. I am trying to get a PCA with environmental factors in the place of species. As I have sites which differ qualitatively (in terms of land use) I wanted to be able to show that difference in the final plot (with different colours). Therefore, I used the method a la Gavin Simpson with the package vegan. So far so good. Here is also the code that I used for that:

with(fish, status)
scl <- -1 ## scaling = -1
colvec <- c("red2", "mediumblue")
plot(pond.pca, type = "n", scaling = scl)
with(fish, points(pond.pca, display = "sites", col = colvec[status], scaling = scl, pch = 21, bg = colvec[status]))
head(with(fish, colvec[status]))
text(pond.pca, display = "species", scaling = scl, cex = 0.8, col = "darkcyan")
with(fish, legend("topright", legend = levels(status), bty = "n", col = colvec, pch = 21, pt.bg = colvec))

The problem arises when I try to put arrows for my environmental variables in the ordination plot. If I use biplot and other functions like ordiplot etc. I ll not be able to keep the different colours for my two types of sites, therefore I don't want to use those. If I use the command here:

plot(envfit(pond.pca, PondEnv38, scaling=-1), add=TRUE, col="black")

I get nice arrows, only the are not aligned (and in some cases are completely opposite) with the environmental variables that I ve given with the code before (line 5). I tried to change the scaling but they just cannot align.

Does anyone know how to deal with that problem?

Any tips would be useful.

Answer 1

It is not clear what you are doing wrong as you don't provide a reproducible example of the problem and I am having difficulty following your description of what is wrong. Here is a fully worked out example for you to follow that does what you seem to being trying to do.

data(varespec)
data(varechem)

ord <- rda(varespec)

set.seed(1)
(fit <- envfit(ord, varechem, perm = 999))

## make up a fake `status`
status <- factor(rep(c("Class1","Class2"), times = nrow(varespec) / 2))

> head(status)
[1] Class1 Class2 Class1 Class2 Class1 Class2

Now plot

layout(matrix(1:2, ncol = 2))
## auto version
plot(fit, add = FALSE)

## manual version with extra things
colvec <-  c("red","green")
scl <- -1
plot(ord, type = "n", scaling = scl)
points(ord, display = "sites", col = colvec[status], pch = (1:2)[status])
points(ord, display = "species", pch = "+")
plot(fit, add = TRUE, col = "black")
layout(1)

Which gives

And all the arrows seem to be pointing as they would if you plotted the envfit object directly.