Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Include Gambin distribution #107

Open
piklprado opened this issue Dec 4, 2015 · 2 comments
Open

Include Gambin distribution #107

piklprado opened this issue Dec 4, 2015 · 2 comments

Comments

@piklprado
Copy link
Contributor

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
Copy link
Contributor Author

Carlos Barboza ([email protected]) 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 andrechalom added this to the sads 1.0.0 milestone Dec 9, 2015
@andrechalom
Copy link
Member

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

2 participants