[R-sig-eco] rank lognormal
Mario José
mariognu-eco at yahoo.com.br
Wed May 30 20:44:29 CEST 2012
Hi Oksanen,
thank you for help. I had not thought about estimated parameters. I need search how to do this.
Converting fitted values to rank abundance from same models is very difficult. I can't found how to converting fitted logserie, for instance, to rank abundance. In May (1975) text, he explain how to do and in Stevens (2009 - Primer of ecology with R) there are a implementation in R. But it result in values to plot and not a list of values ranked. I try too implement May recommendation, but I cant do this.
The recommendation is use rank abundance to compare different fitted models with data. My intention is to use ranks of models to test my data. I can not use AIC approach because I am using a third model (Tokeshi) that not have implementations for this.
Thank you again.
Regards,
Mario
----- Mensagem original -----
> De: Jari Oksanen <jari.oksanen at oulu.fi>
> Para: Mario José <mariognu-eco at yahoo.com.br>; "r-sig-ecology at r-project.org" <r-sig-ecology at r-project.org>
> Cc:
> Enviadas: Quarta-feira, 30 de Maio de 2012 10:45
> Assunto: RE: [R-sig-eco] rank lognormal
>
> Mario,
>
> There are two important differences in your code and in Wilson's text:
>
> 1) \Phi^{-1} (or your fi^-1) is standard normal, i.e, in R qnorm(y, 0, 1)
> instead of your qnorm(y, mulog, sdlog). The mulog and sdlog come back to the
> play with your multiplication and addition in the last line.
>
> 2) the term (S - i + 0.5)/S is the same as reverse ppoints function in R. Now
> you use reverse as the last term of the multiplication, but non-reversed in your
> fi() function. In fact f(p[i]) * p[i] should do the trick, as these two p[i]
> terms must be the same -- one cannot be reversed and another non-reversed.
>
> A working function with minimal changes is:
>
> rank.lognormal <-
> function(x){
> S <- length(x);
> xlog <- log(x);
> p <- rev(ppoints(xlog));
> mulog <- mean(xlog);
> sdlog <- sd(xlog);
> fi <- function(y){ qnorm(y) };
> sapply(1:S, FUN=function(i){ exp(1) ^ (mulog + sdlog * fi(p[i]) * (S-i+0.5)/S)
> });
> }
>
> Then some R things.
>
> 1) do not use ";" to end the line
> 2) exp(1)^(x) == exp(x) -- use the latter form
> 3) you have unnecessary sapply: you can directly work with vectors without
> accessing individual vector elements.
>
> Actually, you write the function as a one-liner, and with some elaboration you
> can streamline it otherwise, but a direct re-implementation of Bastow's
> formulation and a re-incarnation of your original code is:
>
> rank.lognormal <-
> function (x)
> {
> p <- rev(ppoints(x))
> logx <- log(x)
> exp(mean(logx) + sd(logx) * qnorm(p) * p)
> }
>
> This only returns the fitted values, like your original code. You need some
> changes to also return the estimated parameters.
>
> The fitted totals are now closer to observed totals, but not identical.
>
> Cheers, Jari Oksanen
> ________________________________________
> From: r-sig-ecology-bounces at r-project.org [r-sig-ecology-bounces at r-project.org]
> on behalf of Mario José [mariognu-eco at yahoo.com.br]
> Sent: 25 May 2012 16:13
> To: r-sig-ecology at r-project.org
> Subject: [R-sig-eco] rank lognormal
>
> Hi all,
>
> I trying develop a code to implement a recommendation of Wilson et al (1991). In
> they text formula:
>
> "...
>
> LnA' = fitted mean ln abundance; sigma = fitted standar deviation of ln
> abundance.
>
> In ranked-abundance terms:
>
> lnAi = lnA' + sigma * fi^-1 * (S-i+0.5)/S
>
> Where
> S = number of species; Ai = abundance of ith out of S species; fi^-1 =
> inverse cumulative distribution function of a norma distribution, i.e.
> the ln abundance at which the area under the normal curve is the value
> indicated.
> ..."
>
> I developed this code:
>
> rank.lognormal<-function(x){
> S <- length(x);
> xlog <- log(x);
> p <- ppoints(xlog);
> mulog <- mean(xlog);
> sdlog <- sd(xlog);
> fi <- function(y){ qnorm(y, mulog, sdlog) };
> sapply(1:S, FUN=function(i){ exp(1) ^ (mulog + sdlog * fi(p[i]) * (S-i+0.5)/S)
> });
> }
>
> bci<-c(25, 24, 22, 21, 18, 17, 15, 14, 14, 13, 13, 12, 12, 11, 11,
> 10, 9, 8, 7, 6, 6, 6, 6, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4,
> 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
> 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1);
>
> sum(rank.lognormal(bci))
> plot(log(rank.lognormal(bci)))
>
> But the sum of fitted rank don't match with the sum of rank abundance. The
> curve don't like lognormal too.
>
> I know rad.lognormal of vegan package, but I'd like follow Wilson
> recommendation and understand how this procedure works.
>
> Thank you.
>
> Best,
>
> Mario
>
> Wilson et al 1991 Methods for fitting dominance/diversity curves - Journal of
> Vegetation Science.
> http://onlinelibrary.wiley.com/doi/10.2307/3235896/abstract
>
> ...................................................
>
> Master student in Ecology
>
> Institute of Biology, Dept. Plant Biology, Ecology Lab.
> State University of Campinas - UNICAMP
> Campinas, São Paulo, Brazil
>
>
> _______________________________________________
> R-sig-ecology mailing list
> R-sig-ecology at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>
More information about the R-sig-ecology
mailing list