comparison brt.R @ 3:001d7d101915 draft

planemo upload for repository https://github.com/galaxyecology/tools-ecology/tree/master/tools/Ecoregionalization_workflow commit e03df85746a3b61a382a5ee7e3357a8bf42a5097

2:1b9522671e90

#16/02/2023
## Analyse BRT data

### Clean environment 
rm(list = ls(all.names = TRUE))
options(warn=-1)

### load packages

library(dismo, warn.conflicts = FALSE)
library(gbm, warn.conflicts = FALSE)

#Carbo,Grav,Maxbearing,Maxmagnit,Meancurmag,Meansal,Meantheta,Mud,Prof,Rugosity,Sand,Seaice_prod,Sili,Slope,Standcurmag,Standsal,Standtheta

#Load functions

make.brt <- function(spe,data,pred.vars,env,nb_file){
   brt_step <- gbm.step(data= data, gbm.x = pred.vars, gbm.y = spe, family = "bernoulli", tree.complexity = 2, learning.rate = 0.0001,max.trees = 10000,plot.main = F)
   #plot
   if (is.null(brt_step)==FALSE){
     pdf(file = paste("BRT-",spe,".pdf"))
     gbm.plot(brt_step, write.title = T,show.contrib = T, y.label = "fitted function",plot.layout = c(3,3))
     dev.off()

    coord_quickmap()
  output_directory <- ggsave(paste("BRT-",spe,"_pred_plot.png"))
  
  #Write prediction in a file
  preds <- cbind(preds,spe)
  write.table(preds, paste(nb_file,"_brts_pred_ceamarc.tsv",sep=""), quote=FALSE, dec=".", row.names=F, col.names=T,append = T,sep="\t")
}

#### RUN BRT ####
nb_file = 0

for (file in data_files[[1]]) {
  species_data <- read.table(file, dec = dec_species, sep = "\t", header = TRUE, na.strings = "NA", colClasses = "numeric")
  nb_file = nb_file + 1
  `%!in%` <- Negate(`%in%`)
  sp = list()
  for (n in names(species_data)) {
    if (n %!in% names(env) && n != 'station'){
       sp = cbind(sp,n)
    }
  }     
  for (spe in sp){
   try(make.prediction.brt(make.brt(spe,species_data,pred.vars,env,nb_file)))
   }
}

cat("Here is the list of your abiotic parameters:\n")
cat(paste(pred.vars, collapse = ", "), "\n")

3:001d7d101915

#16/02/2023
## Analyse BRT data

### Clean environment 
rm(list = ls(all.names = TRUE))
options(warn=1)

### load packages

library(dismo, warn.conflicts = FALSE)
library(gbm, warn.conflicts = FALSE)

#Carbo,Grav,Maxbearing,Maxmagnit,Meancurmag,Meansal,Meantheta,Mud,Prof,Rugosity,Sand,Seaice_prod,Sili,Slope,Standcurmag,Standsal,Standtheta

#Load functions

make.brt <- function(spe,data,pred.vars,env,nb_file){
   cat(paste("  ", spe,":\n  -> optimising BRT model ",sep=""))
   lr <- 0.05
   no.trees <- 0
   while ( no.trees < 1000 & lr > 0.0005 ) {
    cat(".")
    try(brt_step <- gbm.step(data= data, gbm.x = pred.vars, gbm.y = spe, family = "bernoulli", tree.complexity = 2, learning.rate = lr,max.trees = 10000, plot.main = F))
    # if the gbm does not converge, the return object is null or of size 0
    if (!is.null(brt_step) ) {
      if (object.size(brt_step) > 0 ) {
        no.trees <- brt_step$gbm.call$best.trees
        print(no.trees)
      }
    } else {
      no.trees <- 0
      print(no.trees)
    }

    # decrease the learning rate
    lr <- lr / 2
    print(lr)
  }
   #plot
   if (is.null(brt_step)==FALSE){
     pdf(file = paste("BRT-",spe,".pdf"))
     gbm.plot(brt_step, write.title = T,show.contrib = T, y.label = "fitted function",plot.layout = c(3,3))
     dev.off()

    coord_quickmap()
  output_directory <- ggsave(paste("BRT-",spe,"_pred_plot.png"))
  
  #Write prediction in a file
  preds <- cbind(preds,spe)
  write.table(preds, paste(nb_file,"_brts_pred_ceamarc.tsv",sep=""), quote=FALSE, dec=".", row.names=F, col.names=!file.exists(paste(nb_file,"_brts_pred_ceamarc.tsv",sep="")),append = T,sep="\t")
}

#### RUN BRT ####
nb_file = 0

# Creating the %!in% operator
`%!in%` <- Negate(`%in%`)

# Data file browsing
for (file in data_files[[1]]) {
  
  # Reading the file
  species_data <- read.table(file, dec = dec_species, sep = "\t", header = TRUE, na.strings = "NA", colClasses = "numeric")
  nb_file = nb_file + 1
 
  # List to store species to predict
  sp = list()
  
  # Selection of columns that are not in 'env' and that are not coordinates or stations
  for (n in names(species_data)) {
    if (n %!in% names(env) && n != 'station' && n != 'decimalLatitude' && n != 'decimalLongitude' && n!='lat' && n!='long'){
       sp = c(sp,n) 
    }
  } 
  # Making predictions for each species 
  for (spe in sp){
   try(make.prediction.brt(make.brt(spe,species_data,pred.vars,env,nb_file)))
   }
}

#Display of abiotic parameters
cat("Here is the list of your abiotic parameters:\n")
cat(paste(pred.vars, collapse = ", "), "\n")