1

I am a novice in using R for multivariate analysis . I am trying to get a RDA plot depicting the relationship between my species abundance and environmental data. I have 6 environmental variables. But when I obtain the plot, I am able to see only only two vectors representing two variables alone. The commands I have used are below.

data <- read.csv("all_data.csv",h=T); 
library(vegan)
sp1 <- data[,c("Sample","Acidobacteria","Actinobacteria","Aquificae","Bacteroidetes")];
env1 <- data[,c("Nitrogen","TOC","Phosphate","Sand","Silt","Clay")];
myrda <- rda(sp1,env1)
plot(myrda,scaling=2)

Someone please help me out with this. I wish to see all the 6 environmental parameters in my RDA plot.

thelatemail
  • 91,185
  • 12
  • 128
  • 188
Vijay
  • 11
  • 1
  • 3
  • Please make the code in your question reproducible? – Roman Luštrik Apr 02 '14 at 06:37
  • I have given the exact code what I had tried. Please let me know if more information is required to be added in. – Vijay Apr 02 '14 at 06:41
  • Eh, one time I don't link to the post... Here are a few examples of how to improve a question. http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example – Roman Luštrik Apr 02 '14 at 08:37

1 Answers1

0

Here is an example using vegan's example data varespec and varechem. The plot of the rda model automatically displays all 14 environmental variables:

library(vegan)
data(varespec)
data(varechem)
myrda <- rda(varespec, varechem)
myrda
colnames(varechem) # 14 variables
plot(myrda,scaling=2) # 14 vectors shown

enter image description here

Maybe double check that your data.frames correctly contain variable names so thet the plot knows where to grab labels. I would also make sure that your data splitting is working correctly - I don't think that your method will always work. Here is a possible alternative that should:

sp.incl <- match(c("Sample","Acidobacteria","Actinobacteria","Aquificae","Bacteroidetes"), colnames(data))
sp1 <- data[,sp.incl]

env.incl <- match(c("Nitrogen","TOC","Phosphate","Sand","Silt","Clay"), colnames(data))
env1 <- data[,env.incl]
Marc in the box
  • 11,769
  • 4
  • 47
  • 97
  • I tried the above piece of code, but ended up with the same plot with two vectors alone. The message I could see from myrda is as below. Call: rda(X = phyla, Y = envdata) Inertia Proportion Rank Total 69.53 1.00 Constrained 69.53 1.00 2 Unconstrained 0.00 0.00 0 Inertia is variance Some constraints were aliased because they were collinear (redundant) Eigenvalues for constrained axes: RDA1 RDA2 61.680 7.852 Kindly advise. – Vijay Apr 02 '14 at 09:12
  • You tried with the example vegan data or your own? You may find some sound advice [here](https://stat.ethz.ch/pipermail/r-sig-ecology/2008-June/000195.html), where a similar problem is addressed. The author of vegan, Jari Oksanen, seems to suggest that this is a problem that might be addressed in future versions of vegan. – Marc in the box Apr 02 '14 at 09:27
  • Thanks for that link. Is there any other package of R that can used to get RDA plots? – Vijay Apr 02 '14 at 09:30
  • You can use the base plot functions. A guide from Gavin Simpson can be found [here](http://www.fromthebottomoftheheap.net/2012/04/11/customising-vegans-ordination-plots/). – Marc in the box Apr 02 '14 at 09:40
  • If you look at the output you gave above, you will see you have line `"Constrained 69.53 1.00 2"` which tells you that your constrained ("explained") inertia is 69.53 and all (`1.00`) of that is explained by two (`2`) axes. The next line confirms that there indeed is no residual variance. You have two-dimensions and that's all. However, modern version of `rda` in **vegan** would display all variables unless they are copies or combinations of each other. Please check your environmental data (and version of **vegan**: should be >2). – Jari Oksanen Apr 04 '14 at 09:20