Error when fitting Fisher's log series

[Pakillo]

Hi,

First of all, thanks for the great package!

We are using it to fit SADs to a number of communities but find errors sometimes. A reproducible example:

abundance <- c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 4L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 3L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 3L, 1L, 4L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 
1L, 1L, 1L, 3L, 3L, 1L, 1L, 2L, 1L, 3L, 1L, 1L, 2L, 2L, 1L, 2L, 
2L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 4L, 1L, 2L, 
1L, 2L, 5L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 1L, 1L, 1L, 2L, 1L, 1L, 
1L, 4L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 7L, 7L, 1L, 1L, 
1L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 3L, 1L, 1L, 1L, 1L, 5L, 1L, 1L, 1L, 4L, 1L, 1L, 1L, 
1L, 1L, 1L, 3L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 2L, 1L, 1L, 1L, 5L)

fitsad(abundance, sad = "ls")

Gives error:
Error in uniroot(f1, interval = c(1/N, N)) :
f() values at end points not of opposite sign

Is that expected for some communities in which abundance data don't fit this model, or is there something else going on?

Many thanks in advance

[andrechalom]

Hello!

This is expected sometimes, as the numeric optimization is a somewhat delicate process, and we have tried to provide initial guesses and optimzation parameter that work well for larger published datasets.

If you have a good guess of the model parameters, it helps to provide it as the start.value argument. In your case, the abundances are all very close to 1, so the alpha parameter is very high. Your model is fitted with:

fitsad(abundance, sad = "ls", start.value=200)

Unfortunately, there is no general rule for what to do when the model fitting fails.

[Pakillo]

Thanks for your response!

Interestingly, if we assume that the function is decreasing (which seems reasonable for Fisher's log-series, if not all SADs), we can specify that to uniroot (using extendIt = "downX", see https://stat.ethz.ch/R-manual/R-devel/library/stats/html/uniroot.html), and this avoids the error.

For example, taking the code from your fitls function:

x <- c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 4L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 3L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 3L, 1L, 4L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 
1L, 1L, 1L, 3L, 3L, 1L, 1L, 2L, 1L, 3L, 1L, 1L, 2L, 2L, 1L, 2L, 
2L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 4L, 1L, 2L, 
1L, 2L, 5L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 1L, 1L, 1L, 2L, 1L, 1L, 
1L, 4L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 7L, 7L, 1L, 1L, 
1L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 3L, 1L, 1L, 1L, 1L, 5L, 1L, 1L, 1L, 4L, 1L, 1L, 1L, 
1L, 1L, 1L, 3L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 2L, 1L, 1L, 1L, 5L)

 S <- length(x)
 N <- sum(x)
 if (missing(start.value)){
    f1 <- function(a) {
      S + a*log((a/(a + N)))
 }
 sol <- uniroot(f1, interval = c(1/N, N), extendInt = "downX")
alfa <- sol$root

alfa = 324.533

Do you see that as a good solution? Or there are caveats?

[piklprado]

Your solution reached a higher likelihood value, nice! I am flagging this issue as an 'enhancement'. We wiil check if we can inlude this improvement in the next version. Thanks.