rank lognormal
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