[R-sig-eco] rank lognormal

[Jari Oksanen]

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