Nothing
.Biomod.Models.loop <- function(X,
modeling.id,
Model,
Options,
VarImport,
mod.eval.method,
SavePred,
xy=NULL,
scal.models = TRUE){
cat("\n\n-=-=-=- Run : ",X$name, '\n')
res.sp.run <- list()
for(i in 1:ncol(X$calibLines)){ # loop on RunEval
cat('\n\n-=-=-=--=-=-=-',paste(X$name,dimnames(X$calibLines)[[2]][i],sep=""),'\n')
res.sp.run[[dimnames(X$calibLines)[[2]][i]]] <- lapply(Model, .Biomod.Models,
Data = X$dataBM,
Options = Options,
calibLines = na.omit(X$calibLines[,i,]), ## transform 3D calibLines obj into a 1D vector
Yweights = na.omit(X$Yweights),
nam = paste(X$name,dimnames(X$calibLines)[[2]][i], sep=""),
VarImport = VarImport,
mod.eval.method = mod.eval.method,
evalData = X$evalDataBM,
SavePred = T,#SavePred,
xy = X$xy,
eval.xy = X$eval.xy,
scal.models = scal.models,
modeling.id = modeling.id)
names(res.sp.run[[dimnames(X$calibLines)[[2]][i]]]) <- Model
}
return(res.sp.run)
}
.Biomod.Models <- function (Model, Data, Options, calibLines, Yweights, nam, VarImport = 0,
mod.eval.method = c('ROC','TSS','KAPPA'), evalData = NULL,
SavePred = FALSE,
xy = NULL, eval.xy = NULL, scal.models = TRUE, modeling.id = ''){
################################################################################################
# 1. Print model running and getting model options =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
# check and get modified args if nececary
args <- .Biomod.Models.check(Model, Data, Options, calibLines, Yweights, mod.eval.method, evalData, scal.models)
if(is.null(args)){ # trouble in input data -> Not Run
return(0)
} else {
Data <- args$Data
Yweights <- args$Yweights
evalLines <- args$evalLines
Type <- args$Type
criteria <- args$criteria
Prev <- args$Prev
mod.eval.method <- args$mod.eval.method
evalData <- args$evalData
scal.models <- args$scal.models
resp_name <- args$resp_name
expl_var_names <- args$expl_var_names
compress.arg <- TRUE # ifelse(.Platform$OS.type == 'windows', 'gzip', 'xz')
}
categorial_var <- unlist(sapply(expl_var_names, function(x){if(is.factor(Data[,x])) return(x) else return(NULL)} ))
model_name <- paste(nam,'_',Model,sep="")
# defining the function outputs
ListOut <- list(evaluation = NULL,
var.import = NULL,
pred = NULL,
pred.eval = NULL,
calib.failure = NULL)
# CTA models creation =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
if (Model == "CTA") {
# converting cost argument
if(is.null(Options@CTA$cost)){
cost.tmp <- rep(1,(ncol(Data)-2))
} else{
cost.tmp <- Options@CTA$cost
}
if(Options@CTA$parms == 'default'){
model.sp <- try( rpart(makeFormula(colnames(Data)[1],
head(Data[,-c(1,ncol(Data)), drop=FALSE]),
'simple', 0),
data = Data[calibLines,],
weights = Yweights,
method = Options@CTA$method,
cost = cost.tmp,
control = eval(Options@CTA$control)) )
} else{
model.sp <- try( rpart(makeFormula(colnames(Data)[1],
head(Data[,-c(1,ncol(Data)), drop=FALSE]),
'simple', 0),
data = Data[calibLines,],
weights = Yweights,
method = Options@CTA$method,
parms = Options@CTA$parms,
cost = cost.tmp,
control = eval(Options@CTA$control)) )
}
if( !inherits(model.sp,"try-error") ){
# select best trees --------------- May be done otherway
tr <- as.data.frame(model.sp$cptable)
tr$xsum <- tr$xerror + tr$xstd
tr <- tr[tr$nsplit > 0, ]
Cp <- tr[tr$xsum == min(tr$xsum), "CP"]
model.sp <- prune(model.sp, cp = Cp[length(Cp)])
# creation of biomod2 model object
model.bm <- new("CTA_biomod2_model",
model = model.sp,
model_name = model_name,
model_class = 'CTA',
model_options = Options@CTA,
resp_name = resp_name,
expl_var_names = expl_var_names,
expl_var_type = get_var_type(Data[calibLines,expl_var_names,drop=F]),
expl_var_range = get_var_range(Data[calibLines,expl_var_names,drop=F]))
}
}
# end CTA models creation =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
# GAM models creation =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
if (Model == "GAM"){
# NOTE : To be able to take into account GAM options and weights we have to do a eval(parse(...))
# it's due to GAM implementation ( using of match.call() troubles)
# cat("\n\tUser defined control args building..")
# user.control.list <- Options@GAM$control
#
# if(Options@GAM$algo == 'GAM_gam'){
# default.control.list <- gam::gam.control()
# } else{
# default.control.list <- mgcv::gam.control()
# }
#
# control.list <- lapply(names(default.control.list), function(x){
# if(x %in% names(user.control.list)){
# return(user.control.list[[x]])
# } else {
# return(default.control.list[[x]])
# }
# })
# names(control.list) <- names(default.control.list)
### Old version
if(Options@GAM$algo == 'GAM_gam'){ ## gam package
# package loading
if(isNamespaceLoaded("mgcv")){
if(isNamespaceLoaded("caret")){unloadNamespace("caret")} ## need to unload caret before car
if(isNamespaceLoaded("car")){unloadNamespace("car")} ## need to unload car before mgcv
unloadNamespace("mgcv")
}
if(!isNamespaceLoaded("gam")){requireNamespace("gam", quietly = TRUE)}
cat('\n\t> GAM (gam) modelling...')
gamStart <- eval(parse(text=paste("gam::gam(",colnames(Data)[1] ,"~1 ," ,
" data = Data[calibLines,,drop=FALSE], family = ", Options@GAM$family$family,"(link = '",Options@GAM$family$link,"')",#eval(Options@GAM$family),
", weights = Yweights[calibLines])" ,sep="")))
model.sp <- try( gam::step.gam(gamStart, .scope(Data[1:3,-c(1,ncol(Data))], "gam::s", Options@GAM$k),
data = Data[calibLines,,drop=FALSE],
# keep = .functionkeep,
direction = "both",
trace= Options@GAM$control$trace,
control = Options@GAM$control))#eval(control.list)) )
} else { ## mgcv package
# package loading
# if( ("package:gam" %in% search()) ){ detach("package:gam", unload=TRUE)}
# if( ! ("package:mgcv" %in% search()) ){ require("mgcv",quietly=TRUE) }
if(isNamespaceLoaded("gam")){unloadNamespace("gam")}
if(!isNamespaceLoaded("mgcv")){requireNamespace("mgcv", quietly = TRUE)}
if(is.null(Options@GAM$myFormula)){
cat("\n\tAutomatic formula generation...")
gam.formula <- makeFormula(resp_name,head(Data[,expl_var_names,drop=FALSE]),Options@GAM$type, Options@GAM$interaction.level, k=Options@GAM$k)
} else{
gam.formula <- Options@GAM$myFormula
}
if (Options@GAM$algo == 'GAM_mgcv'){
cat('\n\t> GAM (mgcv) modelling...')
model.sp <- try(mgcv::gam(gam.formula,
data= Data[calibLines,,drop=FALSE],
family= Options@GAM$family,
weights = Yweights,
control = Options@GAM$control))
} else if (Options@GAM$algo == 'BAM_mgcv'){ ## big data.frame gam version
cat('\n\t> BAM (mgcv) modelling...')
model.sp <- try(mgcv::bam(gam.formula,
data=Data[calibLines,,drop=FALSE],
family=Options@GAM$family,
weights = Yweights,
control = Options@GAM$control))
}
}
if( !inherits(model.sp,"try-error") ){
model.bm <- new("GAM_biomod2_model",
model = model.sp,
model_name = model_name,
model_class = 'GAM',
model_subclass = Options@GAM$algo,
model_options = Options@GAM,
resp_name = resp_name,
expl_var_names = expl_var_names,
expl_var_type = get_var_type(Data[calibLines,expl_var_names,drop=F]),
expl_var_range = get_var_range(Data[calibLines,expl_var_names,drop=F]))
}
}
# end GAM models creation =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
# GBM models creation =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
if (Model == "GBM") {
model.sp <- try(gbm(formula = makeFormula(colnames(Data)[1],head(Data)[,expl_var_names,drop=FALSE], 'simple',0),
data = Data[calibLines,,drop=FALSE],
distribution = Options@GBM$distribution,
var.monotone = rep(0, length = ncol(Data)-2), # -2 because of removing of sp and weights
weights = Yweights,
interaction.depth = Options@GBM$interaction.depth,
n.minobsinnode = Options@GBM$n.minobsinnode,
shrinkage = Options@GBM$shrinkage,
bag.fraction = Options@GBM$bag.fraction,
train.fraction = Options@GBM$train.fraction,
n.trees = Options@GBM$n.trees,
verbose = Options@GBM$verbose,
#class.stratify.cv = Options@GBM$class.stratify.cv,
cv.folds = Options@GBM$cv.folds ))##, n.cores=1)) ## to prevent from parallel issues
if( !inherits(model.sp,"try-error") ){
best.iter <- try(gbm.perf(model.sp, method = Options@GBM$perf.method , plot.it = FALSE))
model.bm <- new("GBM_biomod2_model",
model = model.sp,
model_name = model_name,
model_class = 'GBM',
n.trees_optim = best.iter,
model_options = Options@GBM,
resp_name = resp_name,
expl_var_names = expl_var_names,
expl_var_type = get_var_type(Data[calibLines,expl_var_names,drop=F]),
expl_var_range = get_var_range(Data[calibLines,expl_var_names,drop=F]))
}
}
# end GBM models creation =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
# GLM models creation =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
if (Model == "GLM"){
## build the most complete model formula
if(is.null(Options@GLM$myFormula)){
glm.formula <- makeFormula(colnames(Data)[1],head(Data),Options@GLM$type, Options@GLM$interaction.level)
} else{
glm.formula <- Options@GLM$myFormula
}
if(Options@GLM$test != 'none'){
## make the model selection
glmStart <- glm(eval(parse(text=paste(colnames(Data)[1],"~1",sep=""))),
data = Data[calibLines,,drop=FALSE],
family = Options@GLM$family,
control = eval(Options@GLM$control),
weights = Yweights[calibLines],
mustart = rep(Options@GLM$mustart, sum(calibLines)),
model = TRUE)
## remove warnings
warn <- options('warn')
options(warn=-1)
model.sp <- try( stepAIC(glmStart,
glm.formula,
data = Data[calibLines,,drop=FALSE],
direction = "both", trace = FALSE,
k = criteria,
weights = Yweights[calibLines],
steps = 10000,
mustart = rep(Options@GLM$mustart, sum(calibLines))) )
## reexec warnings
options(warn)
} else {
## keep the total model
model.sp <- try( glm(glm.formula,
data = cbind(Data[calibLines,,drop=FALSE],matrix(Yweights[calibLines], ncol=1, dimnames=list(NULL, "Yweights"))),
family = Options@GLM$family,
control = eval(Options@GLM$control),
weights = Yweights,
# mustart = rep(Options@GLM$mustart, sum(calibLines)),
model = TRUE) )
}
if( !inherits(model.sp,"try-error") ){
# print the selected formula
cat("\n\tselected formula : ")
print(model.sp$formula, useSource=FALSE)
model.bm <- new("GLM_biomod2_model",
model = model.sp,
model_name = model_name,
model_class = 'GLM',
model_options = Options@GLM,
resp_name = resp_name,
expl_var_names = expl_var_names,
expl_var_type = get_var_type(Data[calibLines,expl_var_names,drop=F]),
expl_var_range = get_var_range(Data[calibLines,expl_var_names,drop=F]))
}
}
# end GLM models creation =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
# # MARS models creation -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
# if (Model == "MARS"){
# ## deal with nk argument
# ## if not defined, it will be setted up to default mars value i.e max(21, 2 * ncol(x) + 1)
# nk <- Options@MARS$nk
# if(is.null(nk)){
# nk <- max(21, 2 * length(expl_var_names) + 1)
# }
#
# model.sp <- try( mars(x = Data[calibLines,expl_var_names,drop=FALSE],
# y = Data[calibLines,1],
# degree = Options@MARS$degree,
# nk = nk,
# penalty = Options@MARS$penalty,
# thresh = Options@MARS$thresh,
# prune = Options@MARS$prune,
# w = Yweights[calibLines]) )
#
# if( !inherits(model.sp,"try-error") ){
#
# model.bm <- new("MARS_biomod2_model",
# model = model.sp,
# model_name = model_name,
# model_class = 'MARS',
# model_options = Options@MARS,
# resp_name = resp_name,
# expl_var_names = expl_var_names,
# expl_var_type = get_var_type(Data[calibLines,expl_var_names,drop=F]),
# expl_var_range = get_var_range(Data[calibLines,expl_var_names,drop=F]))
# }
# }
# # end MARS models creation -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
# MARS models creation -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
if (Model == "MARS"){
## build the most complete model formula
if(is.null(Options@MARS$myFormula)){
mars.formula <- makeFormula(colnames(Data)[1],head(Data)[, -ncol(Data), drop = FALSE],Options@MARS$type, Options@MARS$interaction.level)
} else{
mars.formula <- Options@MARS$myFormula
}
## deal with nk argument
## if not defined, it will be setted up to default mars value i.e max(21, 2 * ncol(x) + 1)
nk <- Options@MARS$nk
if(is.null(nk)){
# nk <- max(21, 2 * length(expl_var_names) + 1)
nk <- min(200, max(20, 2 * length(expl_var_names))) + 1
}
model.sp <- try(earth(formula = mars.formula,
data = Data[calibLines, , drop=FALSE],
weights = Yweights,
glm = list(family = binomial),
ncross = 0,
keepxy = FALSE,
# degree = Options@MARS$degree,
pmethod = Options@MARS$pmethod,
nprune = Options@MARS$nprune,
nk = nk,
penalty = Options@MARS$penalty,
thresh = Options@MARS$thresh))
if( !inherits(model.sp,"try-error") ){
model.bm <- new("MARS_biomod2_model",
model = model.sp,
model_name = model_name,
model_class = 'MARS',
model_options = Options@MARS,
resp_name = resp_name,
expl_var_names = expl_var_names,
expl_var_type = get_var_type(Data[calibLines,expl_var_names,drop=F]),
expl_var_range = get_var_range(Data[calibLines,expl_var_names,drop=F]))
}
}
# end MARS models creation -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
# FDA models creation =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
if (Model == "FDA") {
model.sp <- try( do.call(fda,
c( list( formula = makeFormula(colnames(Data)[1],head(Data)[,expl_var_names,drop=FALSE], 'simple',0),
data = Data[calibLines,,drop=FALSE],
method = eval(parse(text=call(Options@FDA$method))),
weights = Yweights[calibLines] ),
Options@FDA$add_args) ) )
if( !inherits(model.sp,"try-error") ){
model.bm <- new("FDA_biomod2_model",
model = model.sp,
model_name = model_name,
model_class = 'FDA',
model_options = Options@FDA,
resp_name = resp_name,
expl_var_names = expl_var_names,
expl_var_type = get_var_type(Data[calibLines,expl_var_names,drop=F]),
expl_var_range = get_var_range(Data[calibLines,expl_var_names,drop=F]))
}
}
# end FDA models creation =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
# ANN models creation =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
if (Model == "ANN") {
size = Options@ANN$size
decay = Options@ANN$decay
if(is.null(size) | is.null(decay) | length(size)>1 | length(decay)>1 ){
## define the size and decay to test
if(is.null(size)) size=c(2,4,6, 8)
if(is.null(decay)) decay=c(0.001, 0.01, 0.05, 0.1)
## do cross validation test to find the optimal values of size and decay parameters (prevent from overfitting)
CV_nnet = .CV.nnet(Input = Data[,expl_var_names,drop=FALSE],
Target = Data[calibLines,1],
size=size,
decay=decay,
maxit = Options@ANN$maxit,
nbCV = Options@ANN$NbCV,
W = Yweights[calibLines])
## get the optimised parameters values
decay <- CV_nnet[1, 2]
size <- CV_nnet[1,1]
}
model.sp <- try(nnet(formula = makeFormula(resp_name,head(Data[,expl_var_names,drop=FALSE]), 'simple',0),
data = Data[calibLines,,drop=FALSE],
size = size,
rang = Options@ANN$rang,
decay = decay,
weights=Yweights,
maxit = Options@ANN$maxit,
trace = FALSE))
if( !inherits(model.sp,"try-error") ){
model.bm <- new("ANN_biomod2_model",
model = model.sp,
model_name = model_name,
model_class = 'ANN',
model_options = Options@ANN,
resp_name = resp_name,
expl_var_names = expl_var_names,
expl_var_type = get_var_type(Data[calibLines,expl_var_names,drop=F]),
expl_var_range = get_var_range(Data[calibLines,expl_var_names,drop=F]))
}
}
# ANN models creation =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
# RF models creation =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
if (Model == "RF") {
if(Options@RF$do.classif){
# defining occurences as factor for doing classification and not regression in RF
dTmp <- Data[,1]
Data[,1] <- as.factor(Data[,1])
}
if(Options@RF$mtry == 'default'){
model.sp <- try(randomForest(formula = makeFormula(resp_name,head(Data), 'simple',0),
data = Data[calibLines,],
ntree = Options@RF$ntree,
#mtry = ifelse(Options@RF$ntree == 'default', round((ncol(Data)-1)/2), Options@RF$ntree ),
importance = FALSE,
norm.votes = TRUE,
strata = factor(c(0,1)),
nodesize = Options@RF$nodesize,
maxnodes = Options@RF$maxnodes) )
} else {
model.sp <- try(randomForest(formula = makeFormula(resp_name,head(Data), 'simple',0),
data = Data[calibLines,],
ntree = Options@RF$ntree,
mtry = Options@RF$mtry,
importance = FALSE,
norm.votes = TRUE,
strata = factor(c(0,1)),
nodesize = Options@RF$nodesize,
maxnodes = Options@RF$maxnodes) )
}
if(Options@RF$do.classif){
# canceling occurences class modifications
Data[,1] <- as.numeric(dTmp)
rm(dTmp)
}
if( !inherits(model.sp,"try-error") ){
model.bm <- new("RF_biomod2_model",
model = model.sp,
model_name = model_name,
model_class = 'RF',
model_options = Options@RF,
resp_name = resp_name,
expl_var_names = expl_var_names,
expl_var_type = get_var_type(Data[calibLines,expl_var_names,drop=F]),
expl_var_range = get_var_range(Data[calibLines,expl_var_names,drop=F]))
}
}
# end RF models creation -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
# SRE models creation =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
if (Model == "SRE"){
model.sp <- try(sre(Response = Data[calibLines,1],
Explanatory = Data[calibLines,expl_var_names,drop=FALSE],
NewData = NULL,
Quant = Options@SRE$quant,
return_extremcond=TRUE))
if( !inherits(model.sp,"try-error") ){
model.bm <- new("SRE_biomod2_model",
extremal_conditions = model.sp,
model_name = model_name,
model_class = 'SRE',
model_options = Options@SRE,
resp_name = resp_name,
expl_var_names = expl_var_names,
expl_var_type = get_var_type(Data[calibLines,expl_var_names,drop=F]),
expl_var_range = get_var_range(Data[calibLines,expl_var_names,drop=F]))
}
}
# end SRE models creation =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
# MAXENT.Phillips models creation -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
if (Model == "MAXENT.Phillips"){
MWD <- .Prepare.Maxent.WorkDir(Data, xy, calibLines, nam, VarImport = 0,
evalData, eval.xy, species.name=resp_name,
modeling.id = modeling.id,
background_data_dir = Options@MAXENT.Phillips$background_data_dir)
# run MaxEnt:
cat("\n Running Maxent...")
maxent.cmd <- paste0("java ",
ifelse(is.null(Options@MAXENT.Phillips$memory_allocated),"",paste("-mx",Options@MAXENT.Phillips$memory_allocated,"m",sep="")),
" -jar ", file.path(Options@MAXENT.Phillips$path_to_maxent.jar, "maxent.jar"),
" environmentallayers=\"", MWD$m_backgroundFile,
"\" samplesfile=\"", MWD$m_speciesFile,
"\" projectionlayers=\"", gsub(", ",",",toString(MWD$m_predictFile)),
"\" outputdirectory=\"", MWD$m_outdir, "\"",
" outputformat=logistic ",
# "jackknife maximumiterations=",Options@MAXENT.Phillips$maximumiterations,
ifelse(length(categorial_var),
paste(" togglelayertype=",categorial_var, collapse=" ",sep=""),
""),
" redoifexists",
" visible=", Options@MAXENT.Phillips$visible,
" linear=", Options@MAXENT.Phillips$linear,
" quadratic=", Options@MAXENT.Phillips$quadratic,
" product=", Options@MAXENT.Phillips$product,
" threshold=", Options@MAXENT.Phillips$threshold,
" hinge=", Options@MAXENT.Phillips$hinge,
" lq2lqptthreshold=", Options@MAXENT.Phillips$lq2lqptthreshold,
" l2lqthreshold=", Options@MAXENT.Phillips$l2lqthreshold,
" hingethreshold=", Options@MAXENT.Phillips$hingethreshold,
" beta_threshold=", Options@MAXENT.Phillips$beta_threshold,
" beta_categorical=", Options@MAXENT.Phillips$beta_categorical,
" beta_lqp=", Options@MAXENT.Phillips$beta_lqp,
" beta_hinge=", Options@MAXENT.Phillips$beta_hinge,
" betamultiplier=", Options@MAXENT.Phillips$betamultiplier,
" defaultprevalence=", Options@MAXENT.Phillips$defaultprevalence,
" autorun nowarnings notooltips noaddsamplestobackground")
system(command = maxent.cmd, wait = TRUE, intern = TRUE,
ignore.stdout = FALSE, ignore.stderr = FALSE)
model.bm <- new("MAXENT.Phillips_biomod2_model",
model_output_dir = MWD$m_outdir,
model_name = model_name,
model_class = 'MAXENT.Phillips',
model_options = Options@MAXENT.Phillips,
resp_name = resp_name,
expl_var_names = expl_var_names,
expl_var_type = get_var_type(Data[calibLines,expl_var_names,drop=F]),
expl_var_range = get_var_range(Data[calibLines,expl_var_names,drop=F]))
# for MAXENT.Phillips predicitons are calculated in the same time than models building to save time.
cat("\n Getting predictions...")
g.pred <- try(round(as.numeric(read.csv(MWD$m_outputFile)[,3]) * 1000))
# remove tmp dir
.Delete.Maxent.WorkDir(MWD)
}
# end MAXENT.Phillips models creation -=-=-=-=-=-=-=-=-=-=-=-=-= #
# MAXENT.Tsuruoka models creation -=-=-=-=-=-=-=-=-=-=-=-=-=-=-= #
if(Model == "MAXENT.Tsuruoka"){
model.sp <- try(maxent::maxent(feature_matrix = Data[calibLines, expl_var_names, drop = FALSE],
code_vector = as.factor(Data[calibLines, 1]),
l1_regularizer = Options@MAXENT.Tsuruoka$l1_regularizer,
l2_regularizer = Options@MAXENT.Tsuruoka$l2_regularizer,
use_sgd = Options@MAXENT.Tsuruoka$use_sgd,
set_heldout = Options@MAXENT.Tsuruoka$set_heldout,
verbose = Options@MAXENT.Tsuruoka$verbose))
if( !inherits(model.sp,"try-error") ){
model.bm <- new("MAXENT.Tsuruoka_biomod2_model",
model = model.sp,
model_name = model_name,
model_class = 'MAXENT.Tsuruoka',
model_options = Options@MAXENT.Tsuruoka,
resp_name = resp_name,
expl_var_names = expl_var_names,
expl_var_type = get_var_type(Data[calibLines,expl_var_names,drop=F]),
expl_var_range = get_var_range(Data[calibLines,expl_var_names,drop=F]))
}
}
# end of MAXENT.Tsuruoka models creation -=-=-=-=-=-=-=-=-=-=-=- #
# make prediction =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
if((Model != "MAXENT.Phillips")){
g.pred <- try(predict(model.bm, Data[, expl_var_names, drop = FALSE], on_0_1000 = TRUE))
}
# scale or not predictions =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= #
if(scal.models & !inherits(g.pred,'try-error')){
cat("\n\tModel scaling...")
# model.bm@scaling_model <- try(.scaling_model(g.pred/1000, Data[, 1])) ## without weigths
model.bm@scaling_model <- try( .scaling_model(g.pred/1000, Data[, 1], weights= Yweights )) ## with weights
g.pred <- try(predict(model.bm, Data[,expl_var_names,drop=FALSE], on_0_1000=TRUE))
}
# check predictions existance and stop execution if not ok -=-=- #
test_pred_ok <- TRUE
if (inherits(g.pred,"try-error")) { # model calibration or prdiction failed
test_pred_ok <- FALSE
cat("\n*** inherits(g.pred,'try-error')")
} else if (sum(!is.na(g.pred))<=1){ # only NA predicted
test_pred_ok <- FALSE
cat("\n*** only NA predicted")
} else if(length(unique(na.omit(g.pred))) <=1){ # single value predicted
test_pred_ok <- FALSE
cat("\n*** single value predicted")
}
if(test_pred_ok){
# keep the model name
ListOut$ModelName <- model_name
} else{
# keep the name of uncompleted modelisations
cat("\n ! Note : ", model_name, "failed!\n")
ListOut$calib.failure = model_name
return(ListOut) ## end of function.
}
# make prediction on evaluation data =-=-=-=-=-=-=-=-=-=-=-=-=-= #
if(!is.null(evalData)){
g.pred.eval <- try(predict(model.bm, evalData[,expl_var_names,drop=FALSE], on_0_1000=TRUE))
}
# save predictions -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
if(SavePred){
ListOut$pred <- g.pred
if(exists("g.pred.eval"))
ListOut$pred.eval <- g.pred.eval
}
# Model evaluation stuff =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= #
if(length(mod.eval.method) > 0){
cat("\n\tEvaluating Model stuff...")
## Check no NA in g.pred to avoid evaluation failures
na_cell_id <- which(is.na(g.pred))
if(length(na_cell_id)){
# g.pred.without.na <- g.pred[-na_cell_id]
evalLines <- evalLines[!(evalLines %in% na_cell_id)]
cat('\n\tNote : some NA occurs in predictions')
} #else {
# g.pred.without.na <- g.pred
# }
#Precision = max( (max(g.pred.without.na[evalLines]) - min(g.pred.without.na[evalLines]) ) / 50 , 1) # max 50 steps
cross.validation <- sapply(mod.eval.method,
Find.Optim.Stat,
Fit = g.pred[evalLines],
Obs = Data[evalLines,1])#,Precision = Precision)
rownames(cross.validation) <- c("Testing.data","Cutoff","Sensitivity", "Specificity")
if(exists('g.pred.eval')){
## Check no NA in g.pred to avoid evaluation failures
na_cell_id <- which(is.na(g.pred.eval))
if(length(na_cell_id)){
g.pred.eval.without.na <- g.pred.eval[-na_cell_id]
evalData <- evalData[-na_cell_id,]
cat('\n\tNote : some NA occurs in evaluation predictions')
} else {
g.pred.eval.without.na <- g.pred.eval
}
true.evaluation <- sapply(mod.eval.method,
function(x){
return( Find.Optim.Stat(Stat = x,
Fit = g.pred.eval.without.na,
Obs = evalData[,1],
Fixed.thresh = cross.validation["Cutoff",x]) )
})
cross.validation <- rbind(cross.validation["Testing.data",], true.evaluation)
rownames(cross.validation) <- c("Testing.data","Evaluating.data","Cutoff","Sensitivity", "Specificity")
}
ListOut$evaluation <- t(round(cross.validation,digits=3))
## store results
cross.validation <- t(round(cross.validation,digits=3))
ListOut$evaluation <- cross.validation
model.bm@model_evaluation <- cross.validation
## remove useless objects
rm(list=c('cross.validation') )
}
# End evaluation stuff =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= #
# Variables Importance -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
if (VarImport > 0){ # do Varimp stuff
cat("\n\tEvaluating Predictor Contributions...", "\n")
variables.importance <- variables_importance(model.bm, Data[, expl_var_names,drop=FALSE], nb_rand=VarImport)
model.bm@model_variables_importance <- variables.importance$mat
## we stored only the mean of variables importance run
ListOut$var.import <- round(rowMeans(variables.importance$mat, na.rm=T),digits=3)
## remove useless objects
rm(list=c('variables.importance') )
}
# if (VarImport > 0){ # do Varimp stuff
# # Create data.frame vhere corelation between predictions with and without permutation will be
# # stored
#
# cat("\n\tEvaluating Predictor Contributions...", "\n")
# VarImpTable <- matrix(data = 0, nrow = VarImport, ncol = length(expl_var_names))
# dimnames(VarImpTable) <- list(paste('rand', 1:VarImport, sep=""), expl_var_names)
#
# for(vari in expl_var_names){
# for (run in 1:VarImport) {
# ## create a new dataset with interest variable suffled
# TempDS <- Data[, expl_var_names,drop=FALSE]
# TempDS[, vari] <- sample(TempDS[, vari])
#
# if(Model != "MAXENT.Phillips"){
# ## make projection on suffled dataset
# shuffled.pred <- try(predict(model.bm, TempDS, on_0_1000=TRUE))
# } else{
# ## for MAXENT.Phillips, we have created all the permutation at model building step
# shuffled.pred <- try(round(as.numeric(read.csv(file.path(model.bm@model_output_dir, paste(nam, vari, run, "swd.csv", sep="_")))[,3])*1000) )
# ## scal suffled.pred if necessary
# if(length(getScalingModel(model.bm))){
# shuffled.pred <- try( round(.testnull(object = getScalingModel(model.bm), Prev = 0.5 , dat = data.frame(pred = shuffled.pred/1000) ) *1000) )
# # shuffled.pred <- round(as.numeric(predict(getScalingModel(model.bm), shuffled.pred/1000))*1000)
# }
# ## remove useless files on hard drive
# file.remove(list.files(path=model.bm@model_output_dir,
# pattern=paste(nam, vari, run, "swd", sep="_"),
# full.names=TRUE))
# }
#
# ## test if differences exist between the 2 vectors
# # check predictions existance and stop execution if not ok -=-=- #
# test_shuffled.pred_ok <- TRUE
# if (inherits(shuffled.pred,"try-error")) { # model calibration or prdiction failed
# test_shuffled.pred_ok <- FALSE
# } else if (sum(!is.na(shuffled.pred))<=1){ # only NA predicted
# test_shuffled.pred_ok <- FALSE
# } else if(length(unique(na.omit(shuffled.pred))) <=1){ # single value predicted
# test_shuffled.pred_ok <- FALSE
# } else if(length(shuffled.pred)!= length(g.pred)){
# test_shuffled.pred_ok <- FALSE
# }
#
# if(!test_shuffled.pred_ok){
# cat("\n ! Note : ", model_name, "variable importance for",vari,run,"failed!\n")
# VarImpTable[run,vari] <- 0
# } else{
# if(sum( g.pred != shuffled.pred, na.rm=T) == 0){
# VarImpTable[run,vari] <- 0
# } else {
# ## calculate correlation between vectors as proxy for variables importance
# VarImpTable[run,vari] <- 1 - max(round(cor(x=g.pred, y=shuffled.pred, use="pairwise.complete.obs", method="pearson"),digits=3),0,na.rm=T)
# }
# }
#
# }
# }
#
# ## store results
# model.bm@model_variables_importance <- VarImpTable
# ## we stored only the mean of variables importance run
# ListOut$var.import <- round(apply(VarImpTable, 2, mean, na.rm=T),digits=3)
# }
# End Variables Importance -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
# Model saving step =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
assign(x=paste( nam, Model, sep = "_"),
value= model.bm)
save(list=paste( nam, Model, sep = "_"),
file=file.path(resp_name, "models", modeling.id, paste( nam, Model, sep = "_")),
compress=compress.arg)
# End model saving step =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #
return(ListOut)
}
.Biomod.Models.check <- function(Model, Data, Options, calibLines, Yweights, mod.eval.method, evalData, scal.models, criteria=NULL, Prev=NULL){
# get species and expanatory variables names
resp_name <- colnames(Data)[1]
expl_var_names <- colnames(Data)[-1]
# replace Pseudo absences selected (NA) into true absences (0).. for model computing purpose
if(sum(is.na(Data[,1])))
Data[which(is.na(Data[,1])),1] <- 0
# Calib Lines checking
# & Test if there is absences AND presences in data given
if (sum(!calibLines)>0){ # data are splited into 2 set one for evaluation and an other for evaluation stuff
evalLines <- !calibLines
if(sum(Data[calibLines,1] == 0 ) == 0 || sum(Data[calibLines,1] == 0 ) == sum(calibLines) ||
sum(Data[evalLines,1] == 0) == 0 || sum(Data[evalLines,1] == 0) == sum(evalLines)){
warning(paste(colnames(Data)[1], " ", Model," was switch off because of no both
presences and absences data given",sep=""), immediate.=T)
return(NULL)
}
} else { # all values are taken for cali and valid stuff -----> Not so good but usefull for little data set
evalLines <- calibLines
if(sum(Data[,1] == 0 ) == 0 || sum(Data[,1] == 0 ) == nrow(Data)){
warning(paste(colnames(Data)[1], " ", Model," was switch off because of no both
presences and absences data given (full model)",sep=""), immediate.=T)
return(NULL)
}
}
# weights checking
if(is.null(Yweights)){
Yweights <- rep(1,nrow(Data))
}
if(Model %in% c('GBM', 'CTA', 'ANN', 'FDA', 'GAM', 'MARS')){ # this models required data and weights to be in a same datdaset
Data <- cbind(Data,Yweights)
}
# scaling parameter checking
# never scal SRE
if(Model == "SRE") scal.models <- FALSE
# always scal ANN, FDA
if(Model %in% c("ANN", "FDA") ) scal.models <- TRUE
# models options checking and printing
if (Model == "GLM"){
cat('\nModel=GLM')
if(!is.null(Options@GLM$myFormula)){
cat('\n\tformula = ', paste(Options@GLM$myFormula[2],Options@GLM$myFormula[1],Options@GLM$myFormula[3]))
} else{
cat(' (',Options@GLM$type,'with', ifelse(Options@GLM$interaction.level == 0, 'no interaction )', paste('order',Options@GLM$interaction.level,'interaction level )')))
}
if(Options@GLM$test == "AIC"){
criteria <- 2
cat("\n\tStepwise procedure using AIC criteria")
} else if(Options@GLM$test == "BIC"){
criteria <- log(ncol(Data))
cat("\n\tStepwise procedure using BIC criteria")
} else if(Options@GLM$test == "none"){
criteria <- 0
cat("\n\tNo stepwise procedure")
cat("\n\t! You might be confronted to models convergence issues !")
}
}
if (Model == "GBM") {
cat("\nModel=Generalised Boosting Regression \n")
cat("\t", Options@GBM$n.trees, "maximum different trees and ", Options@GBM$cv.folds,
" Fold Cross-Validation")
set.seed(456) # to be able to refind our trees MAY BE BAD
}
if (Model == "GAM") {
cat("\nModel=GAM")
cat("\n\t",Options@GAM$algo,"algorithm chosen")
# cat("\t", Options@GAM$spline, " Degrees of smoothing")
# if(is.null(Yweights)) Yweights <- rep(1,nrow(Data))
# Data <- cbind(Data,Yweights)
}
if (Model == "CTA") {
cat("\nModel=Classification tree \n")
cat("\t", Options@CTA$control$xval, "Fold Cross-Validation")
set.seed(123) # to be able to refind our trees MAY BE BAD
}
if (Model == "ANN") {
cat("\nModel=Artificial Neural Network \n")
cat("\t", Options@ANN$NbCV, "Fold Cross Validation + 3 Repetitions")
# cat("\tCalibration and evaluation phase: Nb of cross-validations: ",
# ncol(Ids), "\n")
set.seed(555) # to be able to refind our trees MAY BE BAD
}
if (Model == "SRE")
cat("\nModel=Surface Range Envelop")
if (Model == "FDA"){
cat("\nModel=Flexible Discriminant Analysis")
}
if (Model == "MARS"){
cat("\nModel=Multiple Adaptive Regression Splines")
if(!is.null(Options@MARS$myFormula)){
cat('\n\tformula = ', paste(Options@MARS$myFormula[2],Options@MARS$myFormula[1],Options@MARS$myFormula[3]))
} else{
cat(' (',Options@MARS$type,'with', ifelse(Options@MARS$interaction.level == 0, 'no interaction )', paste('order',Options@MARS$interaction.level,'interaction level )')))
}
cat("\n")
}
if (Model == "RF"){
cat("\nModel=Breiman and Cutler's random forests for classification and regression")
set.seed(71)
}
if(Model == 'MAXENT.Phillips'){
cat('\nModel=MAXENT.Phillips')
}
if(Model == 'MAXENT.Tsuruoka'){
cat('\nModel=MAXENT.Tsuruoka')
}
# if (Model == "GLM" | Model == "GAM")
# Prev <- sum(DataBIOMOD[, i + Biomod.material$NbVar])/nrow(DataBIOMOD)
## not exactly same as before
if (Model == "GLM" | Model == "GAM"){
Prev <- sum(Data[,1], na.rm=T)/length(Data[,1])
}
# Evaluation Check
available.eval.meth <- c('ROC','KAPPA','TSS','ACCURACY','BIAS','POD','FAR','POFD','SR','CSI',
'ETS','HK','HSS','OR','ORSS')
# if( Model %in% c('SRE') ) available.eval.meth <- available.eval.meth[which(available.eval.meth!='ROC')]
if(sum(!(mod.eval.method %in% available.eval.meth)) > 0 ){
warnings(paste(toString(mod.eval.method[!which(mod.eval.method %in% available.eval.meth)]),
' were switched off !', sep='' ), imediate = TRUE)
}
mod.eval.method <- mod.eval.method[which(mod.eval.method %in% available.eval.meth)]
return(list(Data=Data,
Yweights=Yweights,
evalLines=evalLines,
criteria=criteria,
Prev=Prev,
mod.eval.method=mod.eval.method,
evalData=evalData,
scal.models=scal.models,
resp_name=resp_name,
expl_var_names=expl_var_names))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.