guiblanchet / HMSC

Hierarchical modelling of species community
40 stars 13 forks source link

Error in Rsquared() for poisson model #13

Closed laurajanegraham closed 6 years ago

laurajanegraham commented 6 years ago

If the model is of family = poisson or family = overPoisson, the call to Rsquared() will result in the following error:

Error in Rsquared(model, averageSp = FALSE) : object 'R2' not found

I realise from the documentation now that it is not intended to calculate R2 for these models (and I think it also makes sense not to include this).

A more informative error message here would be really helpful though, for example:

Rsquared <- function (hmsc, newdata = NULL, averageSp = TRUE) 
{
  if (is.null(newdata)) {
    Y <- hmsc$data$Y
  }
  if (!is.null(newdata)) {
    Y <- newdata$Y
  }
  nsp <- ncol(Y)
  Ypred <- predict(hmsc, newdata = newdata)
  if (any(class(hmsc) == "probit")) {
    Y0 <- Y == 0
    Y1 <- Y == 1
    R2 <- numeric()
    for (i in 1:nsp) {
      R2[i] <- mean(Ypred[Y1[, i], i]) - mean(Ypred[Y0[, 
                                                       i], i])
    }
  }
  if (any(class(hmsc) == "gaussian")) {
    ssY <- colSums(scale(Y, scale = FALSE)^2)
    ssRes <- colSums((Y - predict(hmsc))^2)
    R2 <- 1 - ssRes/ssY
  }
  if (any(class(hmsc) %in% c("poisson", "overPoisson"))) {
    stop("R2 may only be calculated for probit and gaussian models!")
  }

  if (averageSp) {
    R2 <- mean(R2)
  }
  return(R2)
}
guiblanchet commented 6 years ago

Thanks for your suggestion! I just add your suggestion to the code of the R2squared function.

laurajanegraham commented 6 years ago

Thanks :) Great R package by the way.