[R-sig-eco] Fwd: [R] cca standard error species

[Jari Oksanen]

Paloma,

Firstly, I explain what the code snippet is supposed to do, and then I
wonder your specific error message.

On 3/05/10 14:41 PM, "prb" <eco.prb at gmail.com> wrote:

>  Dear all,
> 
> First of all I want to apologise for send the previous message to the wrong
> list. I am very grateful to Dr. J. Oksanen for his reply.
> 
> As J. Oksanen answered:
> 
> plot(vare.cca, scaling=2, type="n")
> text(vare.cca, dis="species", scaling=2, cex=0.6)
> with(varespec, ordiellipse(vare.cca, rep(1, 24), show= TRUE, display="lc",
> scaling=2, w= Cet.niv, col=2))
> 
The ordiellipse function adds so-called "species sd" to the graph. If
species have a Gaussian response to the ordination space, this "species sd"
gives an estimate of Gaussian width parameter. In the example above this
"species sd" is drawn for species Cet.niv (a lichen that used to be called
Cetraria nivalis until taxonomists changed all the names). This happens so
that ordiellipse() finds weighted variances and covariances for scores using
species abundances as weights ('w = Cet.niv'). For this we must tell that
all observations belong to one single group ('groups = rep(1, 24)' when the
data have 24 rows, and 'groups' vector will be 24 elements of ones). In
order to have these covariance ellipses to coincide with species scores we
must base the analysis on LC (linear combination) scores ('display = "lc"')
and appropriate scaling where species are weighted averages of sites
('scaling = 2'). This was behind all the fine details of arguments above
('show' is actually not necessary, and 'col' only sets the colour of the
ellipse).

Was this what you meant with "species sd"?

>> From this function I assume that the standard error will be supplied with
> ordiellipse function. I refered as "se" (standard error) as the error that
> will be calculated for species centroids displayed in the previous function
> with *dis="species"*, some similar to confidence intervals. In the papers
> that I read the mention it like "*Data represents means +-SE".*
> 
The usual practice is to show species sd, not their "se". The "se" can be
calculated, but it does not have the meaning of standard error of species
scores. Finding this needs much more complicated code and some unpublished
science (but can be approximated). Species sd refers to the width of species
response and not to the reliability or accuracy of species ordination
scores.

> I am deepening in the parameters of the function ordiellipse and I think
> that will be useful to define some arguments like kind:
> 
> with(varespec, ordiellipse(vare.cca, rep(1, 24), show= TRUE,
> display="lc",kind="se"
> scaling=2, w= Cet.niv, col=2))
> 
> However, when I applied this function to my particular data set, it returns
> an error message like:
> 
> Error in cov.wt(X, W) : weights must be non-negative and not all zero
> 
One thing that you must adapt is to use the correct length of 'groups'
vector (I should have been more explicit here). So the command is better
written as:

with(varespec, ordiellipse(vare.cca, groups = rep(1, nrow(varespec)), ...)

That is, you cannot use rep(1, 24) for any data, but the latter number (24)
should be equal to the number of rows in your data sets. If you didn't adapt
the command to your data dimensions, you may get the error message you
reported. This happens if the species has only zeros in the first 24 rows.
If you adopted the command to use the correct number of rows in your data
and still got the error message, then believe what R says: you should
non-negative and not all zero weights, where weights are the abundances for
the species you want to plot.

Cheers, Jari Oksanen

> I looked for an error like this, but I do not find any answer yet. I would
> appreciate any comment or references to this matter. Again, I apologize for
> my previous e-mail send to wrong list.
> 
> Best regards,
> 
> Paloma
> 
>  ---------- Forwarded message ----------
> prb <eco.prb <at> gmail.com> writes:
>> Dear all,
>> 
>> I realised a correspondence analysis with function cca() of vegan library.
>> Just like in Okansen (2010) in the example of R help:
>> 
>> library(vegan)
>> data(varespec)
>> data(varechem)
>> vare.cca<-cca(varespec~ Al + P + K, varechem)
>> 
>> With plot.cca() function I represented the species matrix in the next way:
>> plot(vare.cca,display="species")
>> 
>> Being similar to:
>> plot((c(-2,2)),(c(-2,2)), type="n", xlab="AXIS 1", ylab="AXIS 2")
>> points(vare.cca,display='species')
>> 
>> Now, I want to add to this graph a standard error in 'x' and 'y'
> direction.
>> How can I add it? I tried different options but I do not get it. It will
> be
>> just like the function ordiellipse of ellipse library and others like
>> ordihull of vegan library, but without add any factor level. The program
> CAP
>> of Anderson realice this kind of output and I am wondering if in R exists
>> some function implemented that deals with it. Some examples of these
>> graphics will be found in Quero et al., 2008 Basic and Applied Ecology 9:
>> 635-644; Maestre et al., 2009 Ecology Letters 12: 930 - 941.
>> 
> Paloma,
> 
> 
> This is a question on specific package, and I'm not sure that general R News
> is
> a suitable forum to answer this. I suggest you re-send this message to
> r-sig-ecol... at r-project.org (see special interest groups at
> 
> http://www.r-project.org/mail.html for more info). Although vegan indeed is
> a
> special, questions related to vegan are regularly discussed there and the
> subject may be tolerated at that forum (or we may assume as long as we are
> not
> told to go away). Another alternative is to use vegan fora at
> http://vegan.r-forge.r-project.org/ which are dedicated to vegan questions.
> 
> Then a brief preliminary aswer to your question: you should tell us what is
> "species standard error". I can guess what you mean, and if I guess
> correctly,
> this can be done. Here is one way how to do this for one species continuing
> your
> example:
> 
> plot(vare.cca, scaling=2, type="n")
> text(vare.cca, dis="species", scaling=2, cex=0.6)
> with(varespec, ordiellipse(vare.cca, rep(1, 24), show= TRUE, display="lc",
> scaling=2, w= Cet.niv, col=2))
> 
> Why this works and what this means is better explained in some more
> dedicated
> news group: there are many fine details. This is also based on the
> hypothesis
> that I guessed correctly what you wanted to do.
> 
> Cheers, Jari Okansen
> 
> [[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