Biplot или построение частичного RDA с помощью ggplot2. что не работает со стрелками, вычисленными с помощью var_fit?

Частичное решение: как упоминалось в комментарии, чтобы имитировать поведение базы R, необходимо вычислить коэффициент умножения, чтобы правильно масштабировать стрелки. Это также было предложено в этом ответе из вопроса, ссылка на который есть в исходном посте.

Однако есть еще проблема, которую я не решил. Для базы R требуется коэффициент масштабирования, который может быть равен 1 или 2. Кажется, мое решение работает только с scaling = 2, как я показываю ниже. Я не уверен, почему scaling = 1 до сих пор не работает, но я веду расследование.

Обновлено: следуя комментарию ниже, я думаю, что решил проблему масштабирования. Кажется, теперь все работает:

# load libraries
library("ggplot2")
library("phyloseq")
library("vegan")
library("scales")

# set scaling factor
scaling_factor <- 1

graphics.off()

# load data
data(soilrep)

# add some random env vars
set.seed(0)
meta_data <- as(sample_data(soilrep), "data.frame")
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)

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

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

# 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")

# set levels
# as factor
meta_data$Treatment <- factor(meta_data$Treatment, levels=unique(meta_data$Treatment))
meta_data$warmed <- factor(meta_data$warmed, levels=unique(meta_data$warmed))

# base plotting
dev.new()

plot(prda, type = "none", scaling=scaling_factor)
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=scaling_factor))
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 = "Treatment"))

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

# 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 multiplying factors for the arrows
arrow_factor <- ordiArrowMul(var_fit)
# get arrows and scale them
arrows <- as.data.frame(vegan::scores(var_fit, display = "vectors"), scaling=scaling_factor) * arrow_factor
arrows$names <- rownames(arrows)

# select arrows which have meaningful p-value, if needed
# for this example, we keep a threshold of 0.99
threshold_p <- 0.99
signif_arrows <- names(which(var_fit$vectors$pvals<=threshold_p))
# subset
print_arrows <- arrows[signif_arrows, ]

# ggplotting
dev.new()

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

Обновлено:

дальнейшее дополнение, поскольку экспорт кажется сложным из-за некоторых других факторов, которые вводятся plot.envfit, как сообщается на странице справки. Проблема в том, что при построении графика с базой R есть параметр, т. е. arrow.mul, который, я думаю, должен быть установлен как коэффициент стрелки, используемый для умножения. Из help(envfit):

Длина стрелок для подобранных векторов регулируется автоматически. для физического размера графика, а длины стрелок не могут быть сравнивали по участкам. Для аналогичного масштабирования стрелок необходимо явно установите аргумент «arrow.mul» в команде «plot»; видеть «ordiArrowMul» и «ordiArrowTextXY».

Итак, это решение, которое даст один и тот же график, независимо от того, какой метод (базовый R или ggplot2) используется:

# load libraries
library("ggplot2")
library("phyloseq")
library("vegan")
library("scales")

# set scaling factor
scaling_factor <- 1

# load data
data(soilrep)

# add some random env vars
set.seed(0)
meta_data <- as(sample_data(soilrep), "data.frame")
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)

# export metadata
write.table(meta_data, "metadata_prints_to_screen.tsv", quote=F, row.names=F, sep = "\t")

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

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

# get multiplying factors for the arrows
arrow_factor <- ordiArrowMul(var_fit)

# 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")

# set levels
# as factor
meta_data$Treatment <- factor(meta_data$Treatment, levels=unique(meta_data$Treatment))
meta_data$warmed <- factor(meta_data$warmed, levels=unique(meta_data$warmed))

# new dev
dev.new()

# plot with base R
plot(prda, type = "none", scaling=scaling_factor)
with(meta_data, plot(var_fit, add=T, col = "black", p.max=0.99, arrow.mul=arrow_factor))
with(meta_data, points(prda, "sites", pch=shape[Treatment], bg=warmed_palette[warmed], col = "black", scaling=scaling_factor))
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 = "Treatment"))

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

# 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 arrows and scale them
arrows <- as.data.frame(vegan::scores(var_fit, display = "vectors"), scaling=scaling_factor) * arrow_factor
arrows$names <- rownames(arrows)

# select arrows which have meaningful p-value, if needed
# for this example, we keep a threshold of 0.99
threshold_p <- 0.99
signif_arrows <- names(which(var_fit$vectors$pvals<=threshold_p))
# subset
print_arrows <- arrows[signif_arrows, ]

# make the plot
export_ggplot <- ggplot(ggplot_rda) +
    geom_point(mapping = aes(x=RDA1, y=RDA2, fill=warmed, shape=treatment)) +
    coord_fixed() +
    scale_shape_manual(values=shape) +
    scale_fill_manual(values=warmed_palette) +
    geom_segment(data=print_arrows, aes(x = 0, xend = RDA1, y = 0, yend = RDA2), arrow = arrow(length = unit(0.25, "cm")), colour = "black") +
    geom_text(data=print_arrows, aes(x = RDA1, y = RDA2, label = names))+
    guides(fill=guide_legend(override.aes=list(shape=21))) +
    theme_bw()

# new dev
dev.new()
plot(export_ggplot)