R/biomod2_classes_5.R

In biomod2: Ensemble Platform for Species Distribution Modeling

##' @name predict.em
##' @author Damien Georges
##' 
##' @title Functions to get predictions from \code{\link{biomod2_ensemble_model}} objects
##' 
##' @description This function allows the user to predict single models from 
##' \code{\link{biomod2_ensemble_model}} on (new) explanatory variables.
##' 
##' 
##' @param object a \code{\link{biomod2_ensemble_model}} object
##' @param newdata a \code{data.frame} or \code{\link[terra:rast]{SpatRaster}} object 
##' containing data for new predictions
##' @param \ldots (\emph{optional}) 
##' 
##' 
##' @seealso \code{\link{biomod2_ensemble_model}}
##' @family Toolbox functions
##' 
##' 
##' @importFrom terra app classify nlyr
##' @importFrom stats qt sd
##' 
NULL

#setGeneric("predict", def = function(object, ...) { standardGeneric("predict") })

# ---------------------------------------------------------------------------- #
# 9. biomod2_ensemble_model --------------------------------------------------
# ---------------------------------------------------------------------------- #

##' @name biomod2_ensemble_model
##' @aliases biomod2_ensemble_model-class
##' @aliases EMmean_biomod2_model-class
##' @aliases EMmedian_biomod2_model-class
##' @aliases EMcv_biomod2_model-class
##' @aliases EMci_biomod2_model-class
##' @aliases EMca_biomod2_model-class
##' @aliases EMwmean_biomod2_model-class
##' @aliases EMmode_biomod2_model-class
##' @aliases EMfreq_biomod2_model-class
##' @author Damien Georges
##' 
##' @title Ensemble model output object class (when running \code{BIOMOD_EnsembleModeling()})
##' 
##' @description Class created by \code{\link{BIOMOD_EnsembleModeling}}
##' 
##' 
##' @slot modeling.id a \code{character} corresponding to the name (ID) of the simulation set 
##' @slot model_name a \code{character} corresponding to the model name
##' @slot model_class a \code{character} corresponding to the model class
##' @slot model_options a \code{list} containing the model options
##' @slot model the corresponding model object
##' @slot scaling_model the corresponding scaled model object
##' @slot dir_name a \code{character} corresponding to the modeling folder
##' @slot resp_name a \code{character} corresponding to the species name
##' @slot expl_var_names a \code{vector} containing names of explanatory variables
##' @slot expl_var_type a \code{vector} containing classes of explanatory variables
##' @slot expl_var_range a \code{list} containing ranges of explanatory variables
##' @slot model_evaluation a \code{data.frame} containing the model evaluations
##' @slot model_variables_importance a \code{data.frame} containing the model variables importance
##' 
##' @param object a \code{\link{biomod2_ensemble_model}} object
##' 
##' @details 
##' 
##' \code{biomod2_model} is the basic object for \pkg{biomod2} ensemble species distribution models. 
##' \cr All listed classes below are derived from \code{biomod2_model}, and have a 
##' \code{model_class} slot specific value :
##' 
##' \itemize{
##'   \item \code{biomod2_ensemble_model} : \code{model_class} is \code{EM}
##'   \item \code{EMmean_biomod2_model} : \code{model_class} is \code{EMmean}
##'   \item \code{EMmedian_biomod2_model} : \code{model_class} is \code{EMmedian}
##'   \item \code{EMcv_biomod2_model} : \code{model_class} is \code{EMcv}
##'   \item \code{EMci_biomod2_model} : \code{model_class} is \code{EMci}
##'   \item \code{EMca_biomod2_model} : \code{model_class} is \code{EMca}
##'   \item \code{EMwmean_biomod2_model} : \code{model_class} is \code{EMwmean}
##'   \item \code{EMmode_biomod2_model} : \code{model_class} is \code{EMmode}
##'   \item \code{EMfreq_biomod2_model} : \code{model_class} is \code{EMfreq}
##' }
##' 
##' 
##' @seealso \code{\link{biomod2_model}}, \code{\link{BIOMOD_EnsembleModeling}}
##' @family Toolbox objects
##' 
##' 
##' @examples
##' 
##' showClass("biomod2_ensemble_model")
##' showClass("EMmean_biomod2_model")
##' showClass("EMmedian_biomod2_model")
##' showClass("EMcv_biomod2_model")
##' showClass("EMci_biomod2_model")
##' showClass("EMca_biomod2_model")
##' showClass("EMwmean_biomod2_model")
##' showClass("EMmode_biomod2_model")
##' showClass("EMfreq_biomod2_model")
##' 
##' 
NULL

##' @name biomod2_ensemble_model-class
##' @rdname biomod2_ensemble_model
##' @export
##' 

## 9.1 Class Definition ---------------------------------------------------------
setClass('biomod2_ensemble_model',
         representation(modeling.id = 'character'), ## maybe some additional args should be added here
         contains = 'biomod2_model',
         prototype = list(model_class = 'EM'),
         validity = function(object) { return(TRUE) })


## 9.2 Show method -------------------------------------------------------------
##' @rdname biomod2_ensemble_model
##' @importMethodsFrom methods show
##' @importFrom methods callNextMethod
##' @export
##' 

setMethod('show', signature('biomod2_ensemble_model'), function(object) { callNextMethod(object) })


## 9.3 predict2 method -------------------------------------------------------------

### biomod2_ensemble_model predict2.em doc -------------------------------------

##' @name predict2.em
##' @aliases predict2.biomod2_ensemble_model.SpatRaster
##' @aliases predict2.biomod2_ensemble_model.data.frame
##' @aliases predict2.EMmean_biomod2_model.SpatRaster
##' @aliases predict2.EMmean_biomod2_model.data.frame
##' @aliases predict2.EMmedian_biomod2_model.SpatRaster
##' @aliases predict2.EMmedian_biomod2_model.data.frame
##' @aliases predict2.EMcv_biomod2_model.SpatRaster
##' @aliases predict2.EMcv_biomod2_model.data.frame
##' @aliases predict2.EMci_biomod2_model.SpatRaster
##' @aliases predict2.EMci_biomod2_model.data.frame
##' @aliases predict2.EMca_biomod2_model.SpatRaster
##' @aliases predict2.EMca_biomod2_model.data.frame
##' @aliases predict2.EMwmean_biomod2_model.SpatRaster
##' @aliases predict2.EMwmean_biomod2_model.data.frame
##' @aliases predict2.EMmode_biomod2_model.SpatRaster
##' @aliases predict2.EMmode_biomod2_model.data.frame
##' @aliases predict2.EMfreq_biomod2_model.SpatRaster
##' @aliases predict2.EMfreq_biomod2_model.data.frame
##' @author Rémi Lemaire-Patin
##' 
##' @title Functions to get predictions from \code{\link{biomod2_ensemble_model}} objects
##' 
##' @description This function allows the user to predict single models from 
##' \code{\link{biomod2_ensemble_model}} on (new) explanatory variables.
##' 
##' 
##' @param object a \code{\link{biomod2_ensemble_model}} object
##' @param newdata a \code{data.frame} or \code{\link[terra:rast]{SpatRaster}} object 
##' containing data for new predictions
##' @param data_as_formal_predictions (\emph{optional, default} \code{FALSE}). A
##' \code{boolean} describing whether \code{newdata} is given as raw environmental 
##' data (\code{FALSE}) or as formal predictions of the individual models 
##' used to build the ensemble model (\code{TRUE}).
##' 
##' @param \ldots (\emph{optional}) 
##' @inheritParams predict2.bm
##' 
##' @seealso \code{\link{biomod2_ensemble_model}}
##' @family Toolbox functions
##' 
##' 
##' @importFrom terra rast app classify writeRaster
##' @keywords internal

NULL

##' 
##' @rdname predict2.em
##' 

setMethod('predict2', signature(object = 'biomod2_ensemble_model', newdata = "SpatRaster"),
          function(object, newdata, predfun, seedval = NULL, ...)
          {
            args <- list(...)
            data_as_formal_predictions <- args$data_as_formal_predictions
            if (is.null(data_as_formal_predictions)) {
              data_as_formal_predictions <- FALSE 
            }
            filename <- args$filename
            overwrite <- args$overwrite
            on_0_1000 <- args$on_0_1000
            na.rm <- args$na.rm
            
            if (is.null(na.rm)) { na.rm <- TRUE }
            if (is.null(overwrite)) { overwrite <- TRUE }
            if (is.null(on_0_1000)) { on_0_1000 <- FALSE }
            
            # additional arg retrieved for EMci
            sd_prediction <- args$sd_prediction
            mean_prediction <- args$mean_prediction
            side <- args$side
            data.type <- args$data.type
            # additional arg retrived for EMca
            thresh <- args$thresh
            # additional arg retrived for EMwmean
            penalization_scores <- args$penalization_scores
            mod.name <- args$mod.name
            
            
            if (data_as_formal_predictions) { 
              newdata <- subset(newdata, object@model)
            } else {
              newdata <- .get_formal_predictions(object, newdata, on_0_1000 = on_0_1000, seedval = seedval)
              #newdata <- apply(newdata, as.numeric) ## ? normalement pas utile ?
            }
            
            out <- predfun(newdata, 
                           on_0_1000 = on_0_1000,
                           mean_prediction = mean_prediction,
                           sd_prediction = sd_prediction, 
                           side = side,
                           data.type = data.type,
                           thresh = thresh, 
                           penalization_scores = penalization_scores,
                           mod.name = mod.name,
                           na.rm = na.rm)
            
            
            if (!is.null(out) & !is.null(filename)) {
              cat("\n\t\tWriting projection on hard drive...")
              if (on_0_1000 & !inherits(object, "EMcv_biomod2_model")) { ## projections are stored as positive integer
                writeRaster(out, filename = filename, overwrite = overwrite, datatype = "INT2S", NAflag = -9999)
              } else { ## keep default data format for saved raster
                writeRaster(out, filename = filename, overwrite = overwrite)
              }
              out <- rast(filename)
            }
            return(out)
          }
)

##' 
##' @rdname predict2.em
##' 

setMethod('predict2', signature(object = 'biomod2_ensemble_model', newdata = "data.frame"),
          function(object, newdata, predfun, seedval = NULL,  ...)
          {
            args <- list(...)
            data_as_formal_predictions <- args$data_as_formal_predictions
            if (is.null(data_as_formal_predictions)) {
              data_as_formal_predictions <- FALSE 
            }
            
            on_0_1000 <- args$on_0_1000
            if (is.null(on_0_1000)) { on_0_1000 <- FALSE }
            na.rm <- args$na.rm
            if (is.null(na.rm)) { na.rm <- TRUE }
            
            # additional arg retrieved for EMci
            sd_prediction <- args$sd_prediction
            mean_prediction <- args$mean_prediction
            side <- args$side
            data.type <- args$data.type
            # additional arg retrived for EMca
            thresh <- args$thresh
            # additional arg retrived for EMwmean
            penalization_scores <- args$penalization_scores
            
            if (data_as_formal_predictions) { 
              newdata <- newdata[ , object@model, drop = FALSE]
            } else {
              newdata <- .get_formal_predictions(object, newdata, on_0_1000 = on_0_1000, seedval = seedval)
              if (object@model_class %in% c("EMmode", "EMfreq")){
                newdata <- apply(newdata, 2, as.character)
              } else {
                newdata <- apply(newdata, 2, as.numeric)
              }
            }
            
            out <- predfun(newdata,
                           on_0_1000 = on_0_1000,
                           mean_prediction = mean_prediction,
                           sd_prediction = sd_prediction,
                           side = side,
                           data.type = data.type,
                           thresh = thresh, 
                           penalization_scores = penalization_scores,
                           na.rm = na.rm)
            return(out)
          }
)

# --------------------------------------------------------------------------- #
# 10.1 biomod2_ensemble_model subclass ---------------------------------------
# ---------------------------------------------------------------------------- #

### -------------------------------------------------------------------------- #
### 10.1 EMmean_biomod2_model ------------------------------------------------
### -------------------------------------------------------------------------- #
##' @name EMmean_biomod2_model-class
##' @rdname biomod2_ensemble_model
##' @export

setClass('EMmean_biomod2_model',
         representation(),
         contains = 'biomod2_ensemble_model',
         prototype = list(model_class = 'EMmean'),
         validity = function(object) { return(TRUE) })

##' 
##' @rdname predict2.em
##' 

setMethod('predict2', signature(object = 'EMmean_biomod2_model', newdata = "SpatRaster"),
          function(object, newdata, ...)
          {
            predfun <- function(newdata, on_0_1000, mod.name, na.rm, ...) {
              if (nlyr(newdata) == 1) {
                return(newdata)
              } else {
                return(
                  app(newdata, function(x) {
                    m <- mean(x, na.rm = na.rm)
                    if (on_0_1000) { m <- round(m) }
                    return(m)
                  }, wopt = list(names = mod.name))
                )
              }
            }
            # redirect to predict2.biomod2_ensemble_model.SpatRaster
            callNextMethod(object, newdata, predfun = predfun, ...)
          }
)

##' 
##' @rdname predict2.em
##' 

setMethod('predict2', signature(object = 'EMmean_biomod2_model', newdata = "data.frame"),
          function(object, newdata, ...)
          {
            predfun <- function(newdata, on_0_1000, na.rm, ...) {
              out <- rowMeans(newdata, na.rm = na.rm)
              if (on_0_1000) { out <- round(out) }
              out
            }
            # redirect to predict2.biomod2_ensemble_model.SpatRaster
            callNextMethod(object, newdata, predfun = predfun, ...)
          }
)


### -------------------------------------------------------------------------- #
### 10.2 EMmedian_biomod2_model ----------------------------------------------
### -------------------------------------------------------------------------- #
##' @name EMmedian_biomod2_model-class
##' @rdname biomod2_ensemble_model
##' @export

setClass('EMmedian_biomod2_model',
         representation(),
         contains = 'biomod2_ensemble_model',
         prototype = list(model_class = 'EMmedian'),
         validity = function(object) { return(TRUE) })

##' 
##' @rdname predict2.em
##' 

setMethod('predict2', signature(object = 'EMmedian_biomod2_model', newdata = "SpatRaster"),
          function(object, newdata, ...)
          {
            predfun <- function(newdata, on_0_1000, mod.name, na.rm, ...) {
              if (nlyr(newdata) == 1) {
                return(newdata)
              } else {
                return(
                  app(newdata, function(x) {
                    m <- median(x, na.rm = na.rm)
                    if (on_0_1000) { m <- round(m) }
                    return(m)
                  }, wopt = list(names = mod.name))
                )
              }
            }
            # redirect to predict2.biomod2_ensemble_model.SpatRaster
            callNextMethod(object, newdata, predfun = predfun, ...)
          }
)

##' 
##' @rdname predict2.em
##' 

setMethod('predict2', signature(object = 'EMmedian_biomod2_model', newdata = "data.frame"),
          function(object, newdata, ...)
          {
            predfun <- function(newdata, on_0_1000, na.rm, ...) {
              out <- apply(newdata, 1, median, na.rm = na.rm)
              if (on_0_1000) { out <- round(out) }
              out
            }
            # redirect to predict2.biomod2_ensemble_model.SpatRaster
            callNextMethod(object, newdata, predfun = predfun, ...)
          }
)

### -------------------------------------------------------------------------- #
### 10.3 EMcv_biomod2_model --------------------------------------------------
### -------------------------------------------------------------------------- #
##' @name EMcv_biomod2_model-class
##' @rdname biomod2_ensemble_model
##' @export

setClass('EMcv_biomod2_model',
         representation(),
         contains = 'biomod2_ensemble_model',
         prototype = list(model_class = 'EMcv'),
         validity = function(object) { return(TRUE) })

##' 
##' @rdname predict2.em
##' 

setMethod('predict2', signature(object = 'EMcv_biomod2_model', newdata = "SpatRaster"),
          function(object, newdata,  ...)
          {
            predfun <- function(newdata, on_0_1000, mod.name, na.rm, ...) {
              if (nlyr(newdata) <= 1) {
                stop(paste0("\n Model EMcv was not computed because only one single model was kept in ensemble modeling (", names(newdata), ")"))
              }
              out <- app(newdata, function(x) {
                sd(x, na.rm = na.rm) / mean(x, na.rm = na.rm) * 100
              }, wopt = list(names = mod.name))
              return(out)
            }
            # redirect to predict2.biomod2_ensemble_model.SpatRaster
            callNextMethod(object, newdata, predfun = predfun, ...)
          }
)

##' 
##' @rdname predict2.em
##' 

setMethod('predict2', signature(object = 'EMcv_biomod2_model', newdata = "data.frame"),
          function(object, newdata, ...)
          {
            predfun <- function(newdata, na.rm, ...) {
              if (ncol(newdata) <= 1) {
                stop(paste0("\n Model EMcv was not computed because only one single model was kept in ensemble modeling ("
                            , colnames(newdata), ")"))
              }
              out <- apply(newdata, 1,
                           function(x) {
                             ifelse(length(x) == 1, 0, 
                                    sd(x, na.rm = na.rm) / mean(x, na.rm = na.rm) * 100)
                           })
              return(out)
            } 
            # redirect to predict2.biomod2_ensemble_model.SpatRaster
            callNextMethod(object, newdata, predfun = predfun, ...)
          }
)

### -------------------------------------------------------------------------- #
### 10.4 EMci_biomod2_model --------------------------------------------------
### -------------------------------------------------------------------------- #
##' @name EMci_biomod2_model-class
##' @rdname biomod2_ensemble_model
##' @export

setClass('EMci_biomod2_model',
         representation(alpha = 'numeric', side = 'character'),
         contains = 'biomod2_ensemble_model',
         prototype = list(model_class = 'EMci', alpha = 0.05, side = 'superior'),
         validity = function(object) {
           if (!(object@side %in% c('inferior','superior'))) {
             stop("side arg should be 'inferior' or 'superior")
           } else { 
             return(TRUE) 
           }
         })

##' 
##' @rdname predict2.em
##' 

setMethod('predict2', signature(object = 'EMci_biomod2_model', newdata = "SpatRaster"),
          function(object, newdata, ...)
          {
            predfun <- function(newdata, on_0_1000, mod.name, na.rm, ...) {
              args <- list(...)
              mean_prediction <- args$mean_prediction
              sd_prediction <- args$sd_prediction
              side <- args$side
              data.type <- args$data.type
              
              if (is.null(mean_prediction)) { 
                mean_prediction <- app(newdata, mean, wopt = list(names = mod.name), na.rm = na.rm)
              }
              if (is.null(sd_prediction)) { 
                sd_prediction <- app(newdata, sd, wopt = list(names = mod.name), na.rm = na.rm) 
              }
              
              ci_prediction <-  switch(
                side,
                inferior = mean_prediction - (sd_prediction * (qt((1 - object@alpha / 2), df = length(object@model) + 1) / sqrt(length(object@model)))),
                superior = mean_prediction + (sd_prediction * (qt((1 - object@alpha / 2), df = length(object@model) + 1) / sqrt(length(object@model))))
              )
              
              if (data.type == "binary") {
                if (on_0_1000) {
                  ci_prediction <- classify(round(ci_prediction),
                                            matrix(c(-Inf, 0, 0, 1000, Inf, 1000),
                                                   nrow = 2, byrow = TRUE))
                } else {
                  ci_prediction <- classify(round(ci_prediction),
                                            matrix(c(-Inf, 0, 0, 1, Inf, 1),
                                                   nrow = 2, byrow = TRUE))
                }
              }
              ci_prediction
            } 
            # redirect to predict2.biomod2_ensemble_model.SpatRaster
            callNextMethod(object, newdata, predfun = predfun, side = object@side, data.type = object@model_type, ...)
          }
)

##' 
##' @rdname predict2.em
##' 

setMethod('predict2', signature(object = 'EMci_biomod2_model', newdata = "data.frame"),
          function(object, newdata, ...)
          {
            predfun <- function(newdata, na.rm, ...) {
              args <- list(...)
              mean_prediction <- args$mean_prediction
              sd_prediction <- args$sd_prediction
              side <- args$side
              data.type <- args$data.type
              on_0_1000 <- args$on_0_1000
              
              if (is.null(mean_prediction)) { 
                mean_prediction <- rowMeans(newdata, na.rm = na.rm)
              }
              if (is.null(sd_prediction)) { 
                sd_prediction <- apply(newdata, 1, sd, na.rm = na.rm)
              }
              
              ci_prediction <- switch(side,
                                      inferior = mean_prediction - (sd_prediction * (qt((1 - object@alpha / 2), df = length(object@model) + 1) / sqrt(length(object@model)))),
                                      superior = mean_prediction + (sd_prediction * (qt((1 - object@alpha / 2), df = length(object@model) + 1) / sqrt(length(object@model)))))
              
              # reclassify prediction to prevent from out of bounds prediction
              if (data.type == "binary") {
                if (on_0_1000) {
                  ci_prediction <- round(ci_prediction * 1000)
                  ci_prediction[ci_prediction > 1000] <- 1000
                  ci_prediction[ci_prediction < 0] <- 0
                } else {
                  ci_prediction[ci_prediction > 1] <- 1
                  ci_prediction[ci_prediction < 0] <- 0
                }
              }
              ci_prediction
            }
            # redirect to predict2.biomod2_ensemble_model.data.frame
            callNextMethod(object, newdata, predfun = predfun, side = object@side, data.type = object@model_type, ...)
          }
)


### -------------------------------------------------------------------------- #
### 10.5 EMca_biomod2_model --------------------------------------------------
### -------------------------------------------------------------------------- #
##' @name EMca_biomod2_model-class
##' @rdname biomod2_ensemble_model
##' @export

setClass('EMca_biomod2_model',
         representation(thresholds = 'numeric'),
         contains = 'biomod2_ensemble_model',
         prototype = list(model_class = 'EMca'),
         validity = function(object) { return(TRUE) })

##' 
##' @rdname predict2.em
##' 

setMethod('predict2', signature(object = 'EMca_biomod2_model', newdata = "SpatRaster"),
          function(object, newdata, data_as_formal_predictions = FALSE, ...)
          {
            args <- list(...)
            on_0_1000 <- args$on_0_1000
            if (is.null(on_0_1000)) { on_0_1000 <- FALSE }
            
            predfun <- function(newdata, on_0_1000, thresh, mod.name, na.rm, ...) {
              if (nlyr(newdata) == 1) {
                return(bm_BinaryTransformation(newdata, thresh))
              } else {
                return(
                  app(bm_BinaryTransformation(newdata, thresh), function(x) {
                    m <- mean(x, na.rm = na.rm)
                    if (on_0_1000) { m <- round(m * 1000) }
                    return(m)
                  }, wopt = list(names = mod.name))
                )
              }
            }
            if (on_0_1000) {
              thresh <- object@thresholds 
            } else {
              thresh <- object@thresholds / 1000 
            }
            # redirect to predict2.biomod2_ensemble_model.SpatRaster
            callNextMethod(object, newdata, predfun = predfun, thresh = thresh,
                           data_as_formal_predictions = data_as_formal_predictions, ...)
          }
)

##' 
##' @rdname predict2.em
##' 

setMethod('predict2', signature(object = 'EMca_biomod2_model', newdata = "data.frame"),
          function(object, newdata, data_as_formal_predictions = FALSE, ...)
          {
            args <- list(...)
            on_0_1000 <- args$on_0_1000
            if (is.null(on_0_1000)) { on_0_1000 <- FALSE }
            
            predfun <- function(newdata, na.rm, ...) {
              out <- rowMeans(bm_BinaryTransformation(newdata, thresh), na.rm = na.rm)
              if (on_0_1000) { out <- round(out * 1000) }
              out
            }
            if (on_0_1000) { 
              thresh <- object@thresholds
            } else {
              thresh <- object@thresholds / 1000
            }
            # redirect to predict2.biomod2_ensemble_model.data.frame
            callNextMethod(object, newdata, predfun = predfun,
                           data_as_formal_predictions = data_as_formal_predictions, ...)
          }
)


### -------------------------------------------------------------------------- #
### 10.6 EMwmean_biomod2_model -----------------------------------------------
### -------------------------------------------------------------------------- #
##' @name EMwmean_biomod2_model-class
##' @rdname biomod2_ensemble_model
##' @export

setClass('EMwmean_biomod2_model',
         representation(penalization_scores = 'numeric'),
         contains = 'biomod2_ensemble_model',
         prototype = list(model_class = 'EMwmean'),
         validity = function(object) { return(TRUE) })

##' 
##' @rdname predict2.em
##' 

setMethod('predict2', signature(object = 'EMwmean_biomod2_model', newdata = "SpatRaster"),
          function(object, newdata, data_as_formal_predictions = FALSE, ...)
          {
            if (ncol(newdata) < 1) {
              stop("Model EMwmean was not computed because no single model was kept in ensemble modeling")
            }
            predfun <- function(newdata, on_0_1000, penalization_scores, mod.name, ...) {
              if (nlyr(newdata) == 1) {
                return(newdata)
              } else {
                return(
                  app(newdata, function(x) {
                    wm <- sum(x * penalization_scores)
                    if (on_0_1000) { wm <- round(wm) }
                    return(wm)
                  }, wopt = list(names = mod.name))
                )
              }
            }
            # redirect to predict2.biomod2_ensemble_model.SpatRaster
            callNextMethod(object, newdata, predfun = predfun,
                           data_as_formal_predictions = data_as_formal_predictions,
                           penalization_scores = object@penalization_scores, ...)
          }
)

##' 
##' @rdname predict2.em
##' 

setMethod('predict2', signature(object = 'EMwmean_biomod2_model', newdata = "data.frame"),
          function(object, newdata, data_as_formal_predictions = FALSE, ...)
          {
            if (ncol(newdata) < 1) {
              stop("Model EMwmean was not computed because no single model was kept in ensemble modeling")
            }
            predfun <- function(newdata, on_0_1000, penalization_scores, na.rm, ...) {
              out <- apply(as.matrix(newdata), 1, function(x) penalization_scores*x)
              if (ncol(newdata) == 1) {
                out <- matrix(out, nrow = 1)
              }
              out <- as.vector(rowSums(t(out), na.rm = na.rm))
              if (on_0_1000) { out <- round(out) }
              out
            }
            # redirect to predict2.biomod2_ensemble_model.SpatRaster
            callNextMethod(object, newdata, predfun = predfun, 
                           penalization_scores = object@penalization_scores,
                           data_as_formal_predictions = data_as_formal_predictions,
                           ...)
          }
)

### -------------------------------------------------------------------------- #
### 10.7 EMmode_biomod2_model ------------------------------------------------
### -------------------------------------------------------------------------- #
##' @name EMmode_biomod2_model-class
##' @rdname biomod2_ensemble_model
##' @export

setClass('EMmode_biomod2_model',
         representation(),
         contains = 'biomod2_ensemble_model',
         prototype = list(model_class = 'EMmode'),
         validity = function(object) { return(TRUE) })

##' 
##' @rdname predict2.em
##' 

setMethod('predict2', signature(object = 'EMmode_biomod2_model', newdata = "SpatRaster"),
          function(object, newdata, ...)
          {
            predfun <- function(newdata, on_0_1000, mod.name, na.rm, ...) {
              if (nlyr(newdata) == 1) {
                return(newdata)
              } else {
                # return(
                #   app(newdata, .find_freq_mode, wopt = list(names = mod.name))
                # )
              
                levels <- unlist(unique(values(newdata[[1]], dataframe = T, na.rm = T)))## On assume qu'il les a tous....
                levels <- levels(levels)
                
                pred <- app(newdata, function(x){
                  m <- .find_mode(x)
                  m <- as.numeric(m)
                  return(m)
                }, wopt = list(names = mod.name))
                pred <- terra::subst(pred, 1:length(levels), levels)
                return(pred)
              }
            }
            # redirect to predict2.biomod2_ensemble_model.SpatRaster
            callNextMethod(object, newdata, predfun = predfun, ...)
          }
)

##' 
##' @rdname predict2.em
##' 

setMethod('predict2', signature(object = 'EMmode_biomod2_model', newdata = "data.frame"),
          function(object, newdata, ...)
          {
            predfun <- function(newdata, on_0_1000, na.rm, ...) {
              out <- unlist(apply(newdata, 1, .find_mode))
              out
            }
            # redirect to predict2.biomod2_ensemble_model.SpatRaster
            callNextMethod(object, newdata, predfun = predfun, ...)
          }
)


### -------------------------------------------------------------------------- #
### 10.8 EMfreq_biomod2_model ------------------------------------------------
### -------------------------------------------------------------------------- #
##' @name EMfreq_biomod2_model-class
##' @rdname biomod2_ensemble_model
##' @export

setClass('EMfreq_biomod2_model',
         representation(),
         contains = 'biomod2_ensemble_model',
         prototype = list(model_class = 'EMfreq'),
         validity = function(object) { return(TRUE) })

##' 
##' @rdname predict2.em
##' 

setMethod('predict2', signature(object = 'EMfreq_biomod2_model', newdata = "SpatRaster"),
          function(object, newdata, ...)
          {
            predfun <- function(newdata, on_0_1000, mod.name, na.rm, ...) {
              if (nlyr(newdata) == 1) {
                return(newdata)
              } else {
                return(
                  app(newdata, .find_freq_mode, wopt = list(names = mod.name))
                )
              }
            }
            # redirect to predict2.biomod2_ensemble_model.SpatRaster
            callNextMethod(object, newdata, predfun = predfun, ...)
          }
)

##' 
##' @rdname predict2.em
##' 

setMethod('predict2', signature(object = 'EMfreq_biomod2_model', newdata = "data.frame"),
          function(object, newdata, ...)
          {
            predfun <- function(newdata, on_0_1000, na.rm, ...) {
              out <- unlist(apply(newdata, 1, .find_freq_mode))
              out
            }
            # redirect to predict2.biomod2_ensemble_model.SpatRaster
            callNextMethod(object, newdata, predfun = predfun, ...)
          }
)