[R-sig-eco] IndVal groups & Plotting an RDA with labels

Tue Jul 3 21:50:50 CEST 2012

Dear Caitlin Porter,

Your plotting code below is fooled up beyond all recognition. Most of your commands are complete mess. I suggest you start with simpler way and work out gradually changing things you want to change. It seems to me that you want to have something very simple and you do not need all that fuss.

Assuming that c1.rda is your result. You can start with

plot(c1.rda)

If it gave you only points and you wanted to have text, you can *add* this one argument (remember, you can recycle and edit the old commands in GUI):

plot(c1.rda, type = "t")

If you also want to change axis labels, then you can add arguments 'xlab' and 'ylab'.

If you don't get species or site labels, an obvious reason is that your data did not have those names. Either check and fix your data, or give those labels in 'labels' argument in text() command. I have a feeling that your data is badly in need of checking before the analysis. Do that first.

Vegan comes with some introductory documents. One of those is called "introduction" (and another is called "FAQ"). You can read those within vegan session using commands

vegandocs("intro")
vegandocs("FAQ")

and they both discuss plotting.

I give some comments on your code below:

On 03/07/2012, at 16:14 PM, Caitlin Porter wrote:
> 
> 
> 2.      My second question relates to plotting an RDA. I have been able to
> run an RDA in vegan package successfully but unable to plot it in a way
> that I can interpret. I need to label sites, species (response matrix) and
> environmental variables/PCA axes (explanatory matrix).
> 
> 
> 
> *#esqplot (MASS package) *
> 
> library(MASS)
> 
> #subset species and sites scores from the rda for first 10 RDA axes
> 
> sr <- scores(c1.rda, display = c("sites", "species"), choices =
> c(1,2,3,4,5,6), scaling = 2)
> 
> sites.only<- as.data.frame(sr$sites)
> 
> srsp <- as.data.frame(sr$species) # data frame with just the species in it
> 
> c1.site <- c1[,1] # object with just the site names from the original data
> set
> 
So your data had no names, but the names were the first variable. Please fix this and life is easier.

> cp.m <- merge(c1.site, sites.only, by=0, sort=FALSE) # merged site object
> with site names

Don't! You should not merge names as a variable: if you do not have names, give them (see command rownames()).
> 
> eqscplot(cp.m$RDA1, cp.m$RDA2, xlim=c(-1, 1), ylim=c(-1, 1), col="blue",
> xlab="RDA Axis 1", ylab="RDA Axis 2", cex=0.3) # defining the variables &
> limits of plot and what the symbols look like

Are you sure your scores are in the range you impose here? (from -1 to 1 in both axes)? I would be surprised, but this could happen. If you do not know the range, do not give it but trust R. 
> 
> text(cp.m$RDA1, cp.m$RDA2, labels=cp.m$x, col="black", cex=0.3) #adding
> names on plots

If this works, this adds names over points of the previous eqscplot(). Why would you like to do this=
> 
> text(srsp$RDA1, srsp$RDA2, labels=rownames(srsp),col= "red", cex=0.3) #adding
> names/species

> # Error in text.default(cp.m$RDA1, cp.m$RDA2, labels = cp.m$x, col =
> "black",  :   zero length 'labels'
> 
The error message claims that you had no rownames(srsp). Do you have rownames in your data? What does colnames(c1) return -- and also check what does rownames(c1) return. Do you have names? (vegan will invent rownames for scores if they were missing in your data, but you seem to use something else here, and the error message and command are inconsistent: error claims that labels were from cp.m$x (uh), but your commend specified rownames(srsp)).
> 
> 
> *#Borcard et al. 2011 plot function*
> 
> plot(c1.rda, main= "Triplot RDA - scaling 2 - wa scores")
> 
> spe.c1 <- scores(c1.rda, choices=1:6, scaling=2, display="sp")
> 
> arrows(0,0,spe.c1[,1], spe.c1[,4], length=0, lty=1, col="red")
> 
> 
The first command here specifies plotting of axis 1 and 2, but arrows ask for axes 1 and 4: make these consistent.

> 
> c1.rda.species <- as.data.frame(c1.rda$CCA$v)
> 
Above you asked species scores with scaling = 2, and here you extract unscaled species scores. This cannot work.

> c1.species1 <- c1.rda.species[abs(c1.rda.species$RDA1)>0.15 |
> abs(c1.rda.species$RDA2)>0.15,]
> 
> text(c1.species1$RDA1, c1.species1$RDA2, labels=rownames(c1.species1),
> cex=0.7, col="red")
> 
> 
> 
> *#plot.cca (Vegan package) *
> 
> plot(c1.rda, choices = c(1,2), display =c("sp", "wa", "cn"), scaling =2,
> type ="text", xlim= c(-1, 1), ylim = c(-1, 1), xlab = RDA Axis 1", ylab=
> "RDA Axis 2", cex=0.3)
> 
Why do you insist on setting xlim, ylim? Are your data really in this range? I doubt. 

You already got a comment on missing opening quote for xlab: this makes unbalanced quotes and gives you errors after errors.

> 
> 
> *# **Error: unexpected symbol in "plot(c1.rda, choices = c(1,2), display
> =c("sp", #"wa", "cn"), scaling =2, type ="text", xlim= c(-1, 1), ylim =
> c(-1, 1), xlab #= RDA Axis"*
> 
> 
> 
> ## S3 method for class 'cca':
> 
> text(c1.rda, display = "sites", labels, choices = c(1, 2), scaling = 2,
> 
>    arrow.mul, head.arrow = 0.05, select, const)
> 
> 
It looks like you have minimally edited copied & pasted vegan help page (here and all the following code). Don't do this but *follow* the documentation. If you give arguments like 'const', 'select', 'arrow.mul', you should give some value to those, or you get an error. Normally you do not need to give these arguments at all as the functions find defaults internally. You only must give them if you want to change those defaults (like 'select' some other cases than all), and then you must give the new values.

> 
> # would I put: scores(c1.rda, choices=1:2, scaling=2, display="sp") in for
> const
> 
This question makes no sense.

It seems that a good idea might be that you sit down and peacefully go through some simple introductory example (like that in vegandocs("intro")). Of course, there is PC-Ord...

Cheers, Jari Oksanen
> 
> 
> ## S3 method for class 'cca':
> 
> points(c1.rda, display = "sites", choices = c(1, 2), scaling = 2,
> 
>    arrow.mul, head.arrow = 0.05, select, const, ...)
> 
> 
> 
> ## S3 method for class 'cca':
> 
> scores(x, choices=c(1,2), display=c("sp","wa","cn"), scaling=2, ...)
> 
> 
> 
> ## S3 method for class 'rda':
> 
> scores(x, choices=c(1,2), display=c("sp","wa","cn"), scaling=2,
> 
>    const, ...)
> 
> 
> 
> ## S3 method for class 'cca':
> 
> summary(object, scaling = 2, axes = 6, display = c("sp", "wa",
> 
>    "lc", "bp", "cn"), digits = max(3, getOption("digits") - 3), ...)
> 
> 
> 
> ## S3 method for class 'summary.cca':
> 
> print(x, digits = x$digits, head = NA, tail = head, ...)
> 
> 
> 
> ## S3 method for class 'summary.cca':
> 
> head(x, n = 6, tail = 0, ...)
> 
> 
> 
> ## S3 method for class 'summary.cca':
> 
> tail(x, n = 6, head = 0, ...)
> 
> * *
> 
> * *
> -- 
> Caitlin Porter
> MSc candidate
> Biology Department
> Saint Mary's University - Halifax
> Office: S320
> (902) 719-4815
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-ecology mailing list
> R-sig-ecology at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

-- 
Jari Oksanen, Dept Biology, Univ Oulu, 90014 Finland
jari.oksanen at oulu.fi, Ph. +358 400 408593, http://cc.oulu.fi/~jarioksa