Plotting a subset of envfit results onto an ordination

Question

I'm working on a figure for a publication where we are looking at a combination of plant coverage and other environmental data on differing communities. I am trying to make a multi-panel figure, with panels that display all the envfit results, one that displays only the plants, and one that displays only the other enviro. Because of the complexity of the figure, it's actually a little easier to construct in the base plot function than in ggvegan.

My challenge is figuring out how to subset the results of the envfit analysis object for the different panels. A simplified example would be:

library(vegan)
data("mite")
data("mite.env")
set.seed(55)
nmds<-metaMDS(mite)
set.seed(55)
ef<-envfit(nmds, mite.env, permu=999)
plot(ef, p.max = .05)

which produces this figure

For sake of the example, does anyone have suggestions on a way I could create two separate figures, one with only the WatrCont vector and one with only the SubsDens vector? I'm sure there is a way to pull specific results out of the ef object, but my coding is not savvy enough to understand how. Additionally, is there a way to have the jumble of text at the center not overlap, similar to jitter in ggplot?

Thank y'all for all of your help!

Answer 1

I would suggest extracting the data from nmds and ef and using ggplot to add the required elements to your plots.

Here is an example:

library(vegan)
library(ggplot2)
data("mite")
data("mite.env")
set.seed(55)
nmds<-metaMDS(mite)
set.seed(55)
ef<-envfit(nmds, mite.env, permu=999)

# Get the NMDS scores
nmds_values <- as.data.frame(scores(nmds))

# Get the coordinates of the vectors produced for continuous predictors in your envfit
vector_coordinates <- as.data.frame(scores(ef, "vectors")) * ordiArrowMul(ef)

# Plot the vectors separately
ggplot(nmds_values,
       aes(x=NMDS1, y = NMDS2)) +
  geom_point() +
  geom_segment(aes(x=0, y=0, xend=NMDS1, yend=NMDS2),
               vector_coordinates[1,]) +
  geom_text(aes(x=NMDS1,y=NMDS2),
            vector_coordinates[1,],
            label=row.names(vector_coordinates[1,]))

ggplot(nmds_values,
       aes(x=NMDS1, y = NMDS2)) +
  geom_point() +
  geom_segment(aes(x=0, y=0, xend=NMDS1, yend=NMDS2),
               vector_coordinates[2,]) +
  geom_text(aes(x=NMDS1,y=NMDS2),
            vector_coordinates[2,],
            label=row.names(vector_coordinates[2,]))

You can play around with the colours, size of the different elements as you see fit. Coordinates for categorical predictors can be extracted in a similar manner.