piLaboratory / sads

R package for fitting species abundance distributions
23 stars 4 forks source link

Include Gambin distribution #107

Open piklprado opened 8 years ago

piklprado commented 8 years ago

Gambin is definetly the most requested model for inclusion in the package. The gambin package seems to fit the model to data aggregated to octaves. I separated this issue from #7 because we have to check this and if so to figure out if a fit to the raw abundance vector is feasible.

piklprado commented 8 years ago

Carlos Barboza (arlosambarboza@gmail.com) got the following from Tom Mathews to make equivalent fits to gambim and poilog, as in http://onlinelibrary.wiley.com/doi/10.1111/ecog.00861/abstract:

logn <- function(data){
  data <- data[data>0]
  gambfit <- fitGambin(data)
  oct <- create_octaves(data)[,2]
  newdata <- rep(0:(length(oct)-1), time=oct)
  fitPLN <- poilogMLE(newdata, zTrunc=TRUE)
  logLikPLN <- fitPLN$logLval
  attr(logLikPLN, "df") <- 2
  attr(logLikPLN, "nobs") <- length(oct)
  class(logLikPLN) <- "logLik"
  pred <- dpoilog(0:(length(oct)-1), fitPLN$par[1], fitPLN$par[2])*length(newdata)
  ret <- gambfit
  ret$Alpha <- NULL
  ret$logLik <- logLikPLN
  ret$fitted.values <- pred
  ret$coefficients <- c("mu"=fitPLN$par[1], "sigma"= exp(fitPLN$par[2]))
  ret$MaxOctave <- length(ret$fitted.values)-1
  b <- BIC(ret)
  c <- AICc(ret)
  res=c(b,c)
  return(res)
}
logn(moths)
andrechalom commented 8 years ago

I'm reading the documentation on Gambin package to see what we can do, but it's hard to make comparisons as so many things are different. For instance, our "octav" method and their "create_octave" functions produce different results:

> octav(1:128)
Object of class "octav"
  octave upper Freq
1      0     1    1
2      1     2    1
3      2     4    2
4      3     8    4
5      4    16    8
6      5    32   16
7      6    64   32
8      7   128   64
9      8   256    0
> create_octaves(1:128)
  octave species
1      0       1
2      1       2
3      2       4
4      3       8
5      4      16
6      5      32
7      6      64
8      7       1

Also, the function above "converts" the poilog fit to a "binned" fit in order to compare them, but we would like to convert the gambin fit to a "unbinned" fit, in order to compare with all our other models.