Mercurial > repos > ecology > ecoregion_brt_analysis

comparison brt.R @ 0:7d77be8fab92 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyecology/tools-ecology/tree/master/tools/Ecoregionalization_workflow commit e03df85746a3b61a382a5ee7e3357a8bf42a5097
author ecology
date Wed, 11 Sep 2024 09:18:37 +0000
parents
children d63d74194100
comparison
equal deleted inserted replaced
-1:000000000000 0:7d77be8fab92
1 #16/02/2023
2 ## Analyse BRT data
3
4 ### Clean environment
5 rm(list = ls(all.names = TRUE))
6 options(warn=1)
7
8 ### load packages
9
10 library(dismo, warn.conflicts = FALSE)
11 library(gbm, warn.conflicts = FALSE)
12 library(ggplot2, warn.conflicts = FALSE)
13
14
15 #load arguments
16 args = commandArgs(trailingOnly=TRUE)
17 if (length(args)==0)
18 {
19 stop("This tool needs at least one argument")
20 }else{
21 enviro <- args[1]
22 species_files <- args[2]
23 abio_para <- args[3]
24 dec_env <- args[8]
25 dec_species <- args[9]
26 }
27
28 ### load data
29
30 env = read.table(enviro, dec = dec_env, header = TRUE, sep="\t", na.strings = "-9999")
31 pred_vars = strsplit(abio_para, ",")[[1]]
32 data_files = strsplit(species_files,",")
33
34 pred.vars <- character(length(pred_vars))
35
36 for (i in seq_along(pred_vars)) {
37 pred_var_col <- as.numeric(pred_vars[i])
38 pred.vars[i] <- names(env)[pred_var_col]}
39
40 #environemental parameters
41 #Carbo,Grav,Maxbearing,Maxmagnit,Meancurmag,Meansal,Meantheta,Mud,Prof,Rugosity,Sand,Seaice_prod,Sili,Slope,Standcurmag,Standsal,Standtheta
42
43 #Load functions
44
45 make.brt <- function(spe,data,pred.vars,env,nb_file){
46 cat(paste(" ", spe,":\n -> optimising BRT model ",sep=""))
47 lr <- 0.05
48 no.trees <- 0
49 while ( no.trees < 1000 & lr > 0.0005 ) {
50 cat(".")
51 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))
52 # if the gbm does not converge, the return object is null or of size 0
53 if (!is.null(brt_step) ) {
54 if (object.size(brt_step) > 0 ) {
55 no.trees <- brt_step$gbm.call$best.trees
56 print(no.trees)
57 }
58 } else {
59 no.trees <- 0
60 print(no.trees)
61 }
62
63 # decrease the learning rate
64 lr <- lr / 2
65 print(lr)
66 }
67 #plot
68 if (is.null(brt_step)==FALSE){
69 pdf(file = paste("BRT-",spe,".pdf"))
70 gbm.plot(brt_step, write.title = T,show.contrib = T, y.label = "fitted function",plot.layout = c(3,3))
71 dev.off()
72 #total deviance explained as (Leathwick et al., 2006)
73 total_deviance <- brt_step$self.statistics$mean.null
74 cross_validated_residual_deviance <- brt_step$cv.statistics$deviance.mean
75 total_deviance_explained <- (total_deviance - cross_validated_residual_deviance)/total_deviance
76 #Validation file
77 valid = cbind(spe,brt_step$cv.statistics$discrimination.mean,brt_step$gbm.call$tree.complexity,total_deviance_explained)
78 write.table(valid, paste(nb_file,"_brts_validation_ceamarc.tsv",sep=""), quote=FALSE, dec=".",sep="\t" ,row.names=F, col.names=F,append = T)}
79
80 return(brt_step)
81 }
82
83 make.prediction.brt <- function(brt_step){
84 #predictions
85 preds <- predict.gbm(brt_step,env,n.trees=brt_step$gbm.call$best.trees, type="response")
86 preds <- as.data.frame(cbind(env$lat,env$long,preds))
87 colnames(preds) <- c("lat","long","Prediction.index")
88 #carto
89 ggplot()+
90 geom_raster(data = preds , aes(x = long, y = lat, fill = Prediction.index))+
91 geom_raster(data = preds , aes(x = long, y = lat, alpha = Prediction.index))+
92 scale_alpha(range = c(0,1), guide = "none")+
93 scale_fill_viridis_c(
94 alpha = 1,
95 begin = 0,
96 end = 1,
97 direction = -1,
98 option = "D",
99 values = NULL,
100 space = "Lab",
101 na.value = "grey50",
102 guide = "colourbar",
103 aesthetics = "fill")+
104 xlab("Longitude") + ylab("Latitude")+ ggtitle(paste(spe,"Plot of BRT predictions"))+
105 theme(plot.title = element_text(size = 10))+
106 theme(axis.title.y = element_text(size = 10))+
107 theme(axis.title.x = element_text(size = 10))+
108 theme(axis.text.y = element_text(size = 10))+
109 theme(axis.text.x = element_text(size = 10))+
110 theme(legend.text = element_text(size = 10))+
111 theme(legend.title = element_text(size = 10))+
112 coord_quickmap()
113 output_directory <- ggsave(paste("BRT-",spe,"_pred_plot.png"))
114
115 #Write prediction in a file
116 preds <- cbind(preds,spe)
117 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")
118 }
119
120 #### RUN BRT ####
121 nb_file = 0
122
123 # Creating the %!in% operator
124 `%!in%` <- Negate(`%in%`)
125
126 # Data file browsing
127 for (file in data_files[[1]]) {
128
129 # Reading the file
130 species_data <- read.table(file, dec = dec_species, sep = "\t", header = TRUE, na.strings = "NA", colClasses = "numeric")
131 nb_file = nb_file + 1
132
133 # List to store species to predict
134 sp = list()
135
136 # Selection of columns that are not in 'env' and that are not coordinates or stations
137 for (n in names(species_data)) {
138 if (n %!in% names(env) && n != 'station' && n != 'decimalLatitude' && n != 'decimalLongitude' && n!='lat' && n!='long'){
139 sp = c(sp,n)
140 }
141 }
142 # Making predictions for each species
143 for (spe in sp){
144 try(make.prediction.brt(make.brt(spe,species_data,pred.vars,env,nb_file)))
145 }
146 }
147
148 #Display of abiotic parameters
149 cat("Here is the list of your abiotic parameters:\n")
150 cat(paste(pred.vars, collapse = ", "), "\n")
151
152