1

I used the Vegan package for RDA and got its plot but I want to customize the plot and make a better one for publication. The plot(rda) of vegan gave a plot and the arrows had normal size (the graph with blue arrows) but when I used the following code to customize the plot, the arrows got shorter and overlapped (plot with red arrows). I played around with the code to add "length" and change the arrow's length but it did not work. I also tried to change xlim and ylim in "plot" to change the plot size but that made everything even more shrunk. A similar question has been asked before but there was no helpful answer for it. Any input would be appreciated.

Original plot with Vigan

Customized plot with short and overlapped arrows

The code used for customizing the plot which results in short arrows:


perc <- round(100*(summary(my_rda2)$cont$importance[2, 1:2]), 2)

sc_si <- scores(my_rda2, display="sites", choices=c(1,2), scaling=1)
sc_sp <- scores(my_rda2, display="species", choices=c(1,2), scaling=1)
sc_bp <- scores(my_rda2, display="bp", choices=c(1, 2), scaling=1)

# Set up a blank plot with scaling, axes, and labels
plot(my_rda2,
     scaling = 1, # set scaling type 
     type = "none", # this excludes the plotting of any points from the results
     frame = FALSE,
     # set axis limits
     xlim = c(-1,1), 
     ylim = c(-1,1),
     # label the plot (title, and axes)
     main = "Triplot RDA - scaling 1",
     xlab = paste0("RDA1 (", perc[1], "%)"), 
     ylab = paste0("RDA2 (", perc[2], "%)") 
)
# add points for site scores
points(sc_si, 
       pch = 21, # set shape (here, circle with a fill colour)
       col = "black", # outline colour
       bg = "steelblue", # fill colour
       cex = 1.2) # size
# add points for species scores
points(sc_sp, 
       pch = 22, # set shape (here, square with a fill colour)
       col = "black",
       bg = "#f2bd33", 
       cex = 1.2)
# add text labels for species abbreviations
text(sc_sp + c(0.03, 0.05), # adjust text coordinates to avoid overlap with points 
     labels = rownames(sc_sp), 
     col = "black", 
     font = 2, # bold
     cex = 0.4)
# add arrows for effects of the expanatory variables
arrows(0,0, # start them from (0,0)
       sc_bp[,1], sc_bp[,2], # end them at the score value
       col = "red",
       lwd = 3)

# add text labels for arrows
text(x = sc_bp[,1] -0.05, # adjust text coordinate to avoid overlap with arrow tip
     y = sc_bp[,2] - 0.03, 
     labels = rownames(sc_bp), 
     col = "red", 
     cex = 1, 
     font = 1)

I expected that the customized plot looks similar to the original one but the arrows are shorter and overlapped. I tried to change the arguments of the "plot" and "arrow" but they did not work.

Gavin Simpson
  • 170,508
  • 25
  • 396
  • 453
Fereshteh Shahoveisi
  • 11
  • 1
  • 1
    What code did you use to generate the default? If you literally just used the `plot()` method, note that the default scaling is not `scaling = 1`. – Gavin Simpson Dec 02 '22 at 08:51
  • Thank you for your response. I will comment the code I used. So, I changed the scaling to 1, that will make the plot() to generate the scaling 1 plot. Right? – Fereshteh Shahoveisi Dec 03 '22 at 17:13
  • my_rda2 <- rda(spe_hellinger ~ Mg + M3 + K + ENR, data = env) RsquareAdj(my_rda)$adj.r.squared # Using all explanatory variables RsquareAdj(my_rda2)$adj.r.squared # Using 4 explanatory variables sqrt(vif.cca(my_rda)) # Multicollinearity before variables selection sqrt(vif.cca(my_rda2)) # Multicollinearity after variables selection anova.cca(my_rda2, by = "axis") # Test which axis are significant anova.cca(my_rda2, by = "terms") # Test which terms are significant my_rda$CA$eig[my_rda$CA$eig > mean(my_rda$CA$eig)] #residual variation plot(my_rda2) – Fereshteh Shahoveisi Dec 03 '22 at 17:15
  • The key point is that the default plot does not draw arrows to coordinate (scores) endpoints, but it scales arrows to fill the current plot area. See help(ordiArrowMul) for documentation. BTW, what do you think scaling = 1 does or "should do"? – Jari Oksanen Dec 04 '22 at 12:28
  • You can *edit* your question to include this additional info by the way. Anyway, I'd missed that you used `arrows()` then `text()`. If you'd used `text(my_rda2, display = "bp")` instead you'd have gotten the correct multiplier that @JariOksanen mentioned. See `?text.cca` and `?points.cca` for lower level helpers, esp the former to see how we compute the multiplier, and `?ordArrowMul` for the detail – Gavin Simpson Dec 05 '22 at 07:35

0 Answers0