library(maptools) library(dismo) library(rgbif) library(XML) library(ggplot2) library(HH) library(rgeos) library(raster) library(biomod2) library(dplyr) library(ggdendro) # set work directory # setwd() varimp <- read.table(file="VarImp.csv",header=T, sep=';', fill=TRUE, check.names=TRUE, stringsAsFactors=FALSE) head(varimp) varimp.max <- tail(sort(varimp$MaxResVol),10) varimp.max <- dplyr::filter(varimp,MaxResVol %in% varimp.max) varimp.max head(variables) presences <- variables[,46:50] #Change numbers if necessary[,46:122] head(presences) suma <- presences%>% summarise_all(sum) suma <- suma%>% select_if(~max(.)>20) suma species <- colnames(suma) species rm(presences,suma) for (i in species) { i2 <- sub("_",".",i) variables$BassinIDmajority <- as.factor(variables$BassinIDmajority) #Manual preselection SRE background_data <- variables sre_i <- sre(Response = background_data[,i], Explanatory = background_data[,names(variables.tabla2)[names(variables.tabla2)!="BassinIDmajority"]], NewData = background_data[,names(variables.tabla2)[names(variables.tabla2)!="BassinIDmajority"]], Quant = 0.025, return_extremcond = F) #table(background_data[,i]) #background_data_unsuit <- variables$MGRS_Id[sre_i==0] # Merge preselected Absences and real presences for model building biomod_modelbuilding_df <- rbind(variables[variables[,i]>0,c(i,"MGRS_Id",names(variables.tabla2))], variables[sre_i==0,c(i,"MGRS_Id",names(variables.tabla2))]) biomod_modelbuilding_df[biomod_modelbuilding_df[,i]==0,i] <- NA summary(biomod_modelbuilding_df) # Create dataframe for model projections biomod_modelprojection_df <- variables[,c(i,"MGRS_Id",names(variables.tabla2))] nrow(biomod_modelbuilding_df) projection_data_ls<- list(biomod_modelprojection_df) names(projection_data_ls) <- "current" scenarios <- c(-0.15,-0.10,-0.05,0.05,0.10,0.15) for (k in scenarios) { scen_name <- paste("scen_",as.character(k),sep="") biomod_modelprojection_df_k <- biomod_modelprojection_df biomod_modelprojection_df_k$MaxResVol <- biomod_modelprojection_df_k$MaxResVol*(1+k) biomod_modelprojection_df_k$LocResCap <- biomod_modelprojection_df_k$LocResCap*(1+k) projection_data_ls[[scen_name]] <- biomod_modelprojection_df_k } head(projection_data_ls[[2]]) class(biomod_modelbuilding_df[,i]) # 1. Formatting Data BiomodData_i <- BIOMOD_FormatingData(resp.var = as.numeric(biomod_modelbuilding_df[,i]), expl.var = biomod_modelbuilding_df[,names(variables.tabla2)], resp.name = i, PA.nb.rep = 3, PA.nb.absences = 1000, #PA.strategy = 'sre', PA.strategy = 'random', #PA.sre.quant = 0.025, na.rm = TRUE, resp.xy = cbind(biomod_modelbuilding_df$MGRS_Id,biomod_modelbuilding_df$MGRS_Id)) # 2. Defining Models Options using default options. BiomodOption <- BIOMOD_ModelingOptions(GAM=list(k=3)) # 3. Run Biomod Models set.seed(999) BiomodModelOut_sp_i <- BIOMOD_Modeling( BiomodData_i, #models = c("GBM"), models = c("GLM","GBM","RF","MAXENT.Phillips"), models.options = BiomodOption, NbRunEval=10, #n cross-validation DataSplit=70, VarImport=0, models.eval.meth = c('TSS','ROC','KAPPA'), do.full.models=FALSE, modeling.id="Test_Model") #write.csv(get_variables_importance(BiomodModelOut_sp_i), # file = paste("H:/Carlos/Doctorat/SDM/SDM/Results/",i2,"/varimp_models_",i2,".csv",sep=""), # row.names=TRUE) save(BiomodModelOut_sp_i, file=paste("H:/Carlos/Doctorat/SDM/SDM/Results/",i2,"/BiomodModelOut_sp_i.out",sep="")) # 4. Projection Models scenario_runs <- c("current",paste("scen_",as.character(scenarios),sep="")) for(j in scenario_runs){ # Individual models projections onenvironmental conditions print(paste("This is: ",j,sep="")) biomod_modelprojection_df_j <- projection_data_ls[[j]] BiomodProjection_sp_i_j <- BIOMOD_Projection( modeling.output = BiomodModelOut_sp_i, new.env = biomod_modelprojection_df_j[,names(variables.tabla2)], xy.new.env = cbind(biomod_modelprojection_df_j$MGRS_Id,biomod_modelprojection_df_j$MGRS_Id), proj.name = j, binary.meth = c('TSS','ROC'), selected.models = "all", compress = TRUE, build.clamping.mask = TRUE) } # 5. Doing Ensemble Modelling BiomodEM_sp_i <- BIOMOD_EnsembleModeling( modeling.output = BiomodModelOut_sp_i, chosen.models = 'all', em.by = 'all', eval.metric = c('ROC'), eval.metric.quality.threshold = c(0.7), models.eval.meth = c('TSS','ROC','KAPPA'), prob.mean = TRUE, prob.cv = TRUE, prob.ci = TRUE, prob.ci.alpha = 0.05, prob.median = TRUE, committee.averaging = TRUE, prob.mean.weight = TRUE, prob.mean.weight.decay = 'proportional', VarImport = 1) write.csv(get_variables_importance(BiomodEM_sp_i), file = paste("H:/Carlos/Doctorat/SDM/SDM/Results/",i2,"/varimp_",i2,".csv",sep=""), row.names=TRUE) save(BiomodEM_sp_i, file=paste("H:/Carlos/Doctorat/SDM/SDM/Results/",i2,"/BiomodEM_sp_i.out",sep="")) # 6. Projecting the ensemble models scenario_runs <- c("current",paste("scen_",as.character(scenarios),sep="")) for(j in scenario_runs){ proj_output_j <- get(load(paste("H:/Carlos/Doctorat/SDM/SDM/Results/",i2,"/proj_",j,"/",i2,".",j,".projection.out",sep=""))) BIOMOD_EnsembleForecasting( projection.output = proj_output_j, EM.output = BiomodEM_sp_i, proj.name=paste(j,"_ensemble_all",sep="")) } # 7. Output for(j in scenario_runs){ BIOMOD_ens_Output <- get(load(paste("H:/Carlos/Doctorat/SDM/SDM/Results/",i2,"/proj_",j,"_ensemble_all","/proj_",j,"_ensemble_all_",i2,"_ensemble.RData",sep=""))) #colnames(BIOMOD_ens_Output) #summary(BIOMOD_ens_Output[,"Sander.lucioperca_EMwmeanByROC_mergedAlgo_mergedRun_mergedData"]) biomod_modelprojection_df_j <- projection_data_ls[[j]] length(biomod_modelprojection_df_j$MGRS_Id) BIOMOD_ens_Output <- cbind(data.frame(biomod_modelprojection_df_j$MGRS_Id),BIOMOD_ens_Output) BIOMOD_ens_Output$biomod_modelprojection_df_j.MGRS_Id <- as.character(BIOMOD_ens_Output$biomod_modelprojection_df_j.MGRS_Id) BIOMOD_ens_Output[,paste(i2,"_EMciByROC_mergedAlgo_mergedRun_mergedData",sep="")] <- BIOMOD_ens_Output[,paste(i2,"_EMciSupByROC_mergedAlgo_mergedRun_mergedData",sep="")]-BIOMOD_ens_Output[,paste(i2,"_EMciInfByROC_mergedAlgo_mergedRun_mergedData",sep="")] summary(BIOMOD_ens_Output) write.csv(BIOMOD_ens_Output, file = paste("H:/Carlos/Doctorat/SDM/SDM/Results/",i2,"/",i2,"_",j,".csv",sep=""), row.names=FALSE) write.table((data.frame("String","Real","Real","Real","Real","Real","Real","Real","Real")), file = paste("H:/Carlos/Doctorat/SDM/SDM/Results/",i2,"/",i2,"_",j,".csvt",sep=""), row.names=FALSE,col.names = FALSE,sep=",") } # 8. Model evaluation scores models_evaluation_i <- get(load(paste("H:/Carlos/Doctorat/SDM/SDM/Results/",i2,"/",i2,".Test_Model.models.out",sep=""))) PresAbs_summary_i <- get(load(paste("H:/Carlos/Doctorat/SDM/SDM/Results/",i2,"/.BIOMOD_DATA/Test_Model/formated.input.data",sep=""))) AUC_data_i <- models_evaluation_i@models.evaluation@val["ROC","Testing.data",,,] mean_AUC_i <- mean(AUC_data_i) sd_AUC_i <- sd(AUC_data_i) TSS_i <- mean(models_evaluation_i@models.evaluation@val["TSS","Testing.data",,,]) Specificity_i <- mean(models_evaluation_i@models.evaluation@val["ROC","Specificity",,,]) Sensitivity_i <- mean(models_evaluation_i@models.evaluation@val["ROC","Sensitivity",,,]) n_Pres_i <- sum(PresAbs_summary_i@data.species==1,na.rm=T) n_PseudoAbs_i <- sum(is.na(PresAbs_summary_i@data.species)) #eval[i, 1] <- i #eval[i, 2] <-mean_AUC_i #eval[i, 3] <-sd_AUC_i #eval[i, 4] <-TSS_i #eval[i, 5] <-Specificity_i #eval[i, 6] <-Sensitivity_i #eval[i, 7] <-n_Pres_i #eval[i, 8] <-n_PseudoAbs_i as.data.frame(cbind(mean_AUC_i, sd_AUC_i, TSS_i, Specificity_i, Sensitivity_i, n_Pres_i, n_PseudoAbs_i)) write.table(as.data.frame(cbind(c("mean_AUC_i","sd_AUC_i","TSS_i","Specificity_i","Sensitivity_i","n_Pres_i","n_PseudoAbs_i"), c(mean_AUC_i,sd_AUC_i,TSS_i,Specificity_i,Sensitivity_i,n_Pres_i,n_PseudoAbs_i))), file = paste("H:/Carlos/Doctorat/SDM/SDM/Results/",i2,"/",i2,"_EVAL.csvt",sep=""), row.names=FALSE,col.names = FALSE,sep=",") # 8. Ensemble Model evaluation scores ens_models_evaluation_i <- get(load(paste("H:/Carlos/Doctorat/SDM/SDM/Results/",i2,"/",i2,".Test_Modelensemble.models.out",sep=""))) ens_models_evaluation_i@em.computed } #write.csv(eval, # file = paste("H:/Carlos/Doctorat/SDM/SDM/Results/Evaluation.csv",sep=""), # row.names=TRUE)