I am trying to plot a partial RDA and some environmental parameters on top of that.

I followed the idea from this post to plot the RDA with ggplot2. I also found this other post which, instead, allows for plotting with base R. I modified it, to add the arrows for environmental parameters.

I am aware of the existence of ggvegan but, even using fortify as suggested, i am unable to get the two plots, in particular the length of the arrows, to match. What am i doing wrong?

If i look at the content of arrows, it looks like ggplot2 is doing it’s job. So is it base R that is plotting arrows wrongly?

Here is the code:

library("ggplot2")
library("phyloseq")
library("vegan")
library("scales")

# load data
data(soilrep)

# add some random env vars
meta_data <- as.data.frame(sample_data(soilrep))
meta_data$env_var_one <- sample(100, nrow(meta_data), replace=T)
meta_data$env_var_two <- sample(100, nrow(meta_data), replace=T)
meta_data$env_var_three <- sample(100, nrow(meta_data), replace=T)

set.seed(1)

# compute partial RDA using formula
prda <- ordinate(soilrep, "RDA", formula=as.formula("soilrep~warmed+warmed"))

# compute fitting
var_fit <- envfit(prda, meta_data[, c("env_var_one", "env_var_two", "env_var_three")], perm=1000, na.rm=TRUE)

# define shapes for treatment legend
shape_legend <- c(15, 16, 17, 18)
# define shapes for treatment in plot
shape <- c(21, 22, 23, 24)
names(shape) <- unique(meta_data$Treatment)

# define colours for warmed
warmed_palette <- hue_pal()(2)[c(2, 1)]
names(warmed_palette) <- c("no", "yes")

# base plotting
dev.new()

par(mar=c(5, 5, 4, 2))
pl <- plot(prda, type="none", scaling=1)
with(meta_data, plot(var_fit, add=T, col="black", p.max=0.99))
with(meta_data, points(prda, "sites", pch=shape[Treatment], bg=warmed_palette[warmed], col="black", scaling=1))
with(meta_data, legend("topright", legend=levels(warmed), fill=warmed_palette, title="warmed"))
with(meta_data, legend("topleft", legend=levels(Treatment), pch=shape_legend, title="Time point"))

# create data frame for ggplot plotting
# get scores for sites
ggplot_rda <- as.data.frame(vegan::scores(prda, display="sites"))

# add some columns for plotting purposes
ggplot_rda$label <- rownames(ggplot_rda)
ggplot_rda$treatment <- meta_data$Treatment
ggplot_rda$warmed <- meta_data$warmed

# as factor
ggplot_rda$treatment <- factor(ggplot_rda$treatment, levels=unique(ggplot_rda$treatment))
ggplot_rda$warmed <- factor(ggplot_rda$warmed, levels=unique(ggplot_rda$warmed))

# get scores for arrows
arrows <- as.data.frame(vegan::scores(var_fit, display = "vectors"))
arrows$names <- rownames(arrows)

# ggplotting
dev.new()

ggplot(ggplot_rda) +
    geom_point(mapping = aes(x=RDA1, y=PC1, fill=warmed, shape=treatment)) +
    coord_fixed(ratio=1) +
    scale_shape_manual(values=shape) +
    scale_fill_manual(values=warmed_palette) +
    geom_segment(data=arrows, aes(x = 0, xend = RDA1, y = 0, yend = PC1), arrow = arrow(length = unit(0.25, "cm"))) +
    geom_text(data=arrows, aes(x = RDA1, y = PC1, label = names))+
    guides(fill=guide_legend(override.aes=list(shape=21)))