diff --git a/Data/Regular_transect.dbf b/Data/Regular_transect.dbf
index 85a249f0f8f55f985ea0fb4b25f4c4feace8ac74..96a8f3fef828e4741aefe82c77d3907df1f34d63 100644
Binary files a/Data/Regular_transect.dbf and b/Data/Regular_transect.dbf differ
diff --git a/Data/Regular_transect.shp b/Data/Regular_transect.shp
index 74a60b2eacd9799e03c1fd9e181082cd00a2cf2b..8355e73eb862abc6d9933cd827a8929e8ae60e72 100644
Binary files a/Data/Regular_transect.shp and b/Data/Regular_transect.shp differ
diff --git a/Data/Regular_transect.shx b/Data/Regular_transect.shx
index a0654096e4123f85dd0bcd9fead31b3c48893988..7b1fd0c16886af138b2489c76aa7f5ab69376d82 100644
Binary files a/Data/Regular_transect.shx and b/Data/Regular_transect.shx differ
diff --git a/Data/transect_distance.tif b/Data/transect_distance.tif
index 83366a1faf1ae4cd7be11d7b633c6be442ff6a11..b1908305baacad9d6daea37a2857157caebdd134 100644
Binary files a/Data/transect_distance.tif and b/Data/transect_distance.tif differ
diff --git a/R_code.Rmd b/R_code.Rmd
index fec86adc48894878adeb77737ab05992e9b6e2df..7e6c9a24af0b69aa4ebbafd19d6fa74a876aad5d 100644
--- a/R_code.Rmd
+++ b/R_code.Rmd
@@ -27,6 +27,9 @@ Rpackage_path<-"/home/anouk/R/x86_64-pc-linux-gnu-library/3.2/"
 ```
 
 ``` {r echo=TRUE, eval=TRUE, tidy=FALSE, message=FALSE, results='hide'}
+Rpackage_path<-"/home/anouk/R/x86_64-pc-linux-gnu-library/3.2/"
+.libPaths(Rpackage_path)
+
 library(shapefiles)
 library(logspline)
 library(rgdal)
@@ -54,6 +57,8 @@ library(grid)
 library(gtable)
 library(gridExtra)
 library(RStoolbox)
+library(DSsim)
+library(cvAUC)
 ##### linux paths ######
 data_path<-"/media/anouk/DATAPART1/These/Data/DATAGTJPROPRE/virtual_sp/code_v_sp/Spatial_sampling_bias/Data/"
 
@@ -616,7 +621,7 @@ plot(region, plot.units = "m")
 design <- make.design(transect.type = "line",
                       design.details = c("Parallel", "Systematic"),
                       region.obj = region,
-                      spacing = 300)
+                      spacing = 400)
 
 transects <- generate.transects(design, region = region)
 plot(region, plot.units = "m")
@@ -664,7 +669,7 @@ writeOGR(transect4, dsn = '.', layer = "Regular_transect", driver = "ESRI Shapef
 transect<-shapefile(paste0(data_path, 'Regular_transect.shp'))
 
 #Use raster of environmental variable in order set the right extent
- a<-raster(extent(new_raster3)) #choose a raster to set up extent
+ a<-raster(extent(new_env_raster)) #choose a raster to set up extent
  res(a)<-c(25,25)
  #Create a raster representing trails (same extent and resolution as environmental variables rasters)
  transectR <- rasterize(transect, a)
@@ -4568,19 +4573,26 @@ row_spec(4:6,  color = "black", background = "#fff")
 *Open environmental data*
 ```{r echo=TRUE, eval=TRUE, tidy=TRUE, message=FALSE}
 
-new_env_raster <- stack(
-                        paste0(data_path, "JURAraster_H.nb10_20relative_density.tif"),
-                        paste0(data_path, "JURAraster_H.nb20_30relative_density.tif"),
-                        paste0(data_path, "JURAraster_H.simpson.tif"),
-                        
+new_env_raster <- stack(paste0(data_path, "JURAraster_MV_1_8ha_Mean_H.nb10_20relative_density.tif"), 
+
+                        paste0(data_path, "JURAraster_MV_1_8ha_Mean_H.nb2_5ratio.tif"),
+
+                        paste0(data_path, "JURAraster_MV_1_8ha_Mean_H.simpson.tif"), 
+
+                        paste0(data_path, "JURAraster_MV_1_8ha_Mean_H.p25.tif"),
+
+                        paste0(data_path, "JURAraster_MV_1_8ha_sd_H.nb20_30relative_density.tif"),
+
+                        paste0(data_path, "JURAraster_MV_1_8ha_sd_H.nb2_5ratio.tif"),
+
 
                         quick=TRUE
 )
 
 
-names(new_env_raster) <- c("Canopy1020", "Canopy2030",  "Simpson")
+names(new_env_raster) <- c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")
 
-#crs(new_env_raster)<-"+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +vunits=m +no_defs"
+crs(new_env_raster)<-"+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +vunits=m +no_defs"
 
 new_env_raster<-normImage(new_env_raster)
 
@@ -4592,7 +4604,7 @@ new_env_raster_tot <- overlay(new_env_raster, rpoly, fun = function(x, y) {
   x[is.na(y[])] <- NA
   return(x)
 })
-names(new_env_raster_tot) <-c( "Canopy1020", "Canopy2030",  "Simpson" )
+names(new_env_raster_tot) <-c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")
 
 predictor_tot<-rasterToPoints(new_env_raster_tot)
 predictor_tot<-na.omit(predictor_tot)
@@ -4650,7 +4662,7 @@ new_env_raster_S <- overlay(new_env_raster, rpoly_S, fun = function(x, y) {
   x[is.na(y[])] <- NA
   return(x)
 })
-names(new_env_raster_S) <-c( "Canopy1020", "Canopy2030",  "Simpson" )
+names(new_env_raster_S) <-c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")
 
 
 # Non Sytematic
@@ -4663,7 +4675,7 @@ new_env_raster_R <- overlay(new_env_raster, rpoly_R, fun = function(x, y) {
   return(x)
 })
 
-names(new_env_raster_R) <-c( "Canopy1020", "Canopy2030",  "Simpson")
+names(new_env_raster_R) <-c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")
 ```
 *Select data (systemtic/Conventional design and Female/male)*
 ```{r echo=TRUE, eval=TRUE, tidy=TRUE, message=FALSE}
@@ -4745,6 +4757,92 @@ pres_pointsGDT<- data.frame(x=OBSGDTshp$coords.x1, y=OBSGDTshp$coords.x2, presen
 
 ```
 
+```{r echo=TRUE, eval=TRUE, tidy=TRUE, message=FALSE}
+#fonction for fold and AUC
+iid <- function(data, V, predictors_vector){
+
+.cvFolds <- function(Y, V){  #Create CV folds (stratify by outcome)
+Y0 <- split(sample(which(Y==0)), rep(1:V, length=length(which(Y==0))))
+Y1 <- split(sample(which(Y==1)), rep(1:V, length=length(which(Y==1))))
+
+folds <- vector("list", length=V)
+
+for (v in seq(V)) {folds[[v]] <- c(Y0[[v]], Y1[[v]])}
+
+return(folds)
+}
+  
+
+.doFit <- function(v, folds, data){  #Train/test glm for each fold
+fit <- maxnet(as.vector(data[-folds[[v]],"presence"]), data[-folds[[v]], predictors_vector], f = FORMULA, regmult = 1)
+pred <- predict(fit, newdata=data[folds[[v]], predictors_vector], type=c("exponential"), clamp=FALSE)
+return(pred)
+}
+
+set.seed(158)
+folds <- .cvFolds(Y=data$presence, V=V)  #Create folds
+
+predictions <- unlist(sapply(seq(V), .doFit, folds=folds, data=data))  #CV train/predict
+
+predictions[unlist(folds)] <- predictions  #Re-order pred va
+
+# Get CV AUC and confidence interval
+out <- ci.cvAUC(predictions=predictions, labels=data$presence, folds=folds, confidence=0.95)
+
+#Get models for output
+# list_fit<-list()
+# 
+# list_fit <- lapply(folds,  function(x) maxnet(as.vector(data[-x,"presence"]), data[-x, predictors_vector], f = FORMULA, regmult = 1) ) 
+# 
+# OUT<-list(out, list_fit)
+# return(OUT)
+
+return(out)
+}
+
+
+iid2 <- function(data, V, predictors_vector){
+
+.cvFolds <- function(Y, V){  #Create CV folds (stratify by outcome)
+Y0 <- split(sample(which(Y==0)), rep(1:V, length=length(which(Y==0))))
+Y1 <- split(sample(which(Y==1)), rep(1:V, length=length(which(Y==1))))
+
+folds <- vector("list", length=V)
+
+for (v in seq(V)) {folds[[v]] <- c(Y0[[v]], Y1[[v]])}
+
+return(folds)
+}
+  
+
+.doFit <- function(v, folds, data){  #Train/test glm for each fold
+fit <- maxnet(as.vector(data[-folds[[v]],"presence"]), data[-folds[[v]], predictors_vector], f = FORMULA, regmult = 1)
+data[, 10][data[, 10] != 0] <- 0
+pred <- predict(fit, newdata=data[folds[[v]], predictors_vector], type=c("exponential"), clamp=FALSE)
+return(pred)
+}
+
+set.seed(158)
+folds <- .cvFolds(Y=data$presence, V=V)  #Create folds
+
+predictions <- unlist(sapply(seq(V), .doFit, folds=folds, data=data))  #CV train/predict
+
+predictions[unlist(folds)] <- predictions  #Re-order pred va
+
+# Get CV AUC and confidence interval
+out <- ci.cvAUC(predictions=predictions, labels=data$presence, folds=folds, confidence=0.95)
+
+#Get models for output
+# list_fit<-list()
+# 
+# list_fit <- lapply(folds,  function(x) maxnet(as.vector(data[-x,"presence"]), data[-x, predictors_vector], f = FORMULA, regmult = 1) ) 
+# 
+# OUT<-list(out, list_fit)
+# return(OUT)
+
+return(out)
+}
+```
 
 ```{r echo=TRUE, eval=TRUE, tidy=TRUE, message=FALSE}
 
@@ -4773,11 +4871,16 @@ pres_bgGDT <- na.omit(pres_bgGDT)
 
 Pres_data<-as.vector(pres_bgGDT$presence)
 
-predictor_maxnet<-pres_bgGDT[, c( "Canopy1020", "Canopy2030",  "Simpson")]
+predictor_maxnet<-pres_bgGDT[, c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")]
 
+FORMULA<-maxnet.formula(Pres_data, predictor_maxnet, classes = "q")
 
 model_maxnet<-maxnet(Pres_data, predictor_maxnet, f = maxnet.formula(Pres_data, predictor_maxnet, classes = "q"), regmult = 1)
 
+var_names<-c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")
+
+out_F_1 <- iid(data=pres_bgGDT, V=10, var_names)
+
 #use it to predict results
 predictor<-rasterToPoints(new_env_raster_S)
 predictor<-na.omit(predictor)
@@ -4832,11 +4935,16 @@ pres_bgGDT <- na.omit(pres_bgGDT)
 
 Pres_data<-as.vector(pres_bgGDT$presence)
 
-predictor_maxnet<-pres_bgGDT[, c( "Canopy1020", "Canopy2030", "Simpson")]
+predictor_maxnet<-pres_bgGDT[, c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")]
 
+FORMULA<-maxnet.formula(Pres_data, predictor_maxnet, classes = "q")
 
 model_maxnet<-maxnet(Pres_data, predictor_maxnet, f = maxnet.formula(Pres_data, predictor_maxnet, classes = "q"), regmult = 1)
 
+var_names<-c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")
+
+out_F_2 <- iid(data=pres_bgGDT, V=10, var_names)
+
 #use it to predict results
 predictor<-rasterToPoints(new_env_raster_S)
 predictor<-na.omit(predictor)
@@ -4924,10 +5032,16 @@ pres_bgGDT <- na.omit(pres_bgGDT)
 
 Pres_data<-as.vector(pres_bgGDT$presence)
 
-predictor_maxnet<-pres_bgGDT[, c( "Canopy1020", "Canopy2030",   "Simpson", "SystTracks_distance5_5")]
+predictor_maxnet<-pres_bgGDT[, c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd", "SystTracks_distance5_5")]
+
+FORMULA<-maxnet.formula(Pres_data, predictor_maxnet, classes = "q")
 
 model_maxnet<-maxnet(Pres_data, predictor_maxnet, f = maxnet.formula(Pres_data, predictor_maxnet, classes = "q"), regmult = 1)
 
+var_names<-c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd", "SystTracks_distance5_5")
+
+out_F_3 <- iid2(data=pres_bgGDT, V=10, var_names)
+
 #use it to predict results
 #Set all values to 0
 raster_tracks_syst_1<-raster_tracks_syst
@@ -4943,7 +5057,7 @@ predictor$predicted<-predict(model_maxnet, predictor,type=c("exponential"), clam
 
 maps_predicted_F_3<-rasterFromXYZ(predictor[, c("x","y","predicted")], res=c(25,25))
 
-colnames(predictor_tot2) <- c("x", "y", "SystTracks_distance5_5", "Canopy1020", "Canopy2030", "Simpson")
+colnames(predictor_tot2) <- c("x", "y", "SystTracks_distance5_5", "Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")
 predictor_tot2$predicted_tot2<-predict(model_maxnet, predictor_tot2,type=c("exponential"), clamp=FALSE)
 
 maps_predicted_F_3b<-rasterFromXYZ(predictor_tot2[, c("x","y","predicted_tot2")], res=c(25,25))
@@ -4993,10 +5107,16 @@ pres_bgGDT <- na.omit(pres_bgGDT)
 
 Pres_data<-as.vector(pres_bgGDT$presence)
 
-predictor_maxnet<-pres_bgGDT[, c( "Canopy1020", "Canopy2030",  "Simpson")]
+predictor_maxnet<-pres_bgGDT[, c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")]
+
+FORMULA<-maxnet.formula(Pres_data, predictor_maxnet, classes = "q")
 
 model_maxnet<-maxnet(Pres_data, predictor_maxnet, f = maxnet.formula(Pres_data, predictor_maxnet, classes = "q"), regmult = 1)
 
+var_names<-c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")
+
+out_F_4 <- iid(data=pres_bgGDT, V=10, var_names)
+
 #use it to predict results
 predictor<-rasterToPoints(new_env_raster_R)
 predictor<-na.omit(predictor)
@@ -5047,11 +5167,16 @@ pres_bgGDT <- na.omit(pres_bgGDT)
 
 Pres_data<-as.vector(pres_bgGDT$presence)
 
-predictor_maxnet<-pres_bgGDT[, c( "Canopy1020", "Canopy2030",  "Simpson")]
+predictor_maxnet<-pres_bgGDT[, c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")]
 
+FORMULA<-maxnet.formula(Pres_data, predictor_maxnet, classes = "q")
 
 model_maxnet<-maxnet(Pres_data, predictor_maxnet, f = maxnet.formula(Pres_data, predictor_maxnet, classes = "q"), regmult = 1)
 
+var_names<-c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")
+
+out_F_5 <- iid(data=pres_bgGDT, V=10, var_names)
+
 #use it to predict results
 predictor<-rasterToPoints(new_env_raster_R)
 predictor<-na.omit(predictor)
@@ -5135,10 +5260,17 @@ pres_bgGDT <- na.omit(pres_bgGDT)
 
 Pres_data<-as.vector(pres_bgGDT$presence)
 
-predictor_maxnet<-pres_bgGDT[, c( "Canopy1020", "Canopy2030", "Simpson", "BiaisTracks_distance5_5")]
+
+predictor_maxnet<-pres_bgGDT[, c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd", "BiaisTracks_distance5_5")]
+
+FORMULA<-maxnet.formula(Pres_data, predictor_maxnet, classes = "q")
 
 model_maxnet<-maxnet(Pres_data, predictor_maxnet, f = maxnet.formula(Pres_data, predictor_maxnet, classes = "q"), regmult = 1)
 
+var_names<-c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd", "BiaisTracks_distance5_5")
+
+out_F_6 <- iid2(data=pres_bgGDT, V=10, var_names)
+
 #use it to predict results
 #Set all values to 1
 raster_tracks_bias_1<-raster_tracks_bias
@@ -5152,7 +5284,7 @@ predictor<-as.data.frame(predictor)
 
 predictor$predicted<-predict(model_maxnet, predictor,type=c("exponential"), clamp=FALSE)
 
-colnames(predictor_tot2) <- c("x", "y", "BiaisTracks_distance5_5", "Canopy1020", "Canopy2030", "Simpson")
+colnames(predictor_tot2) <- c("x", "y", "BiaisTracks_distance5_5", "Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")
 
 predictor_tot2$predicted_tot2<-predict(model_maxnet, predictor_tot2,type=c("exponential"), clamp=FALSE)
 
@@ -5211,10 +5343,16 @@ pres_bgGDT <- na.omit(pres_bgGDT)
 
 Pres_data<-as.vector(pres_bgGDT$presence)
 
-predictor_maxnet<-pres_bgGDT[, c( "Canopy1020", "Canopy2030", "Simpson")]
+predictor_maxnet<-pres_bgGDT[, c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")]
+
+FORMULA<-maxnet.formula(Pres_data, predictor_maxnet, classes = "q")
 
 model_maxnet<-maxnet(Pres_data, predictor_maxnet, f = maxnet.formula(Pres_data, predictor_maxnet, classes = "q"), regmult = 1)
 
+var_names<-c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")
+
+out_M_1 <- iid(data=pres_bgGDT, V=10, var_names)
+
 #use it to predict results
 predictor<-rasterToPoints(new_env_raster_S)
 predictor<-na.omit(predictor)
@@ -5261,11 +5399,16 @@ pres_bgGDT <- na.omit(pres_bgGDT)
 
 Pres_data<-as.vector(pres_bgGDT$presence)
 
-predictor_maxnet<-pres_bgGDT[, c( "Canopy1020", "Canopy2030",  "Simpson")]
+predictor_maxnet<-pres_bgGDT[, c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")]
 
+FORMULA<-maxnet.formula(Pres_data, predictor_maxnet, classes = "q")
 
 model_maxnet<-maxnet(Pres_data, predictor_maxnet, f = maxnet.formula(Pres_data, predictor_maxnet, classes = "q"), regmult = 1)
 
+var_names<-c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")
+
+out_M_2 <- iid(data=pres_bgGDT, V=10, var_names)
+
 #use it to predict results
 predictor<-rasterToPoints(new_env_raster_S)
 predictor<-na.omit(predictor)
@@ -5312,10 +5455,18 @@ pres_bgGDT <- na.omit(pres_bgGDT)
 
 Pres_data<-as.vector(pres_bgGDT$presence)
 
-predictor_maxnet<-pres_bgGDT[, c( "Canopy1020", "Canopy2030",  "Simpson", "SystTracks_distance5_5")]
+model_maxnet<-maxnet(Pres_data, predictor_maxnet, f = maxnet.formula(Pres_data, predictor_maxnet, classes = "q"), regmult = 1)
+
+predictor_maxnet<-pres_bgGDT[, c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd", "SystTracks_distance5_5")]
+
+FORMULA<-maxnet.formula(Pres_data, predictor_maxnet, classes = "q")
 
 model_maxnet<-maxnet(Pres_data, predictor_maxnet, f = maxnet.formula(Pres_data, predictor_maxnet, classes = "q"), regmult = 1)
 
+var_names<-c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd", "SystTracks_distance5_5")
+
+out_F_3 <- iid2(data=pres_bgGDT, V=10, var_names)
+
 predictor<-rasterToPoints(new_env_rastersyst3)
 predictor<-na.omit(predictor)
 predictor<-as.data.frame(predictor)
@@ -5323,7 +5474,7 @@ predictor<-as.data.frame(predictor)
 predictor$predicted<-predict(model_maxnet, predictor,type=c("exponential"), clamp=FALSE)
 
 #predictor_tot2<-predictor_tot2[, c(-7,-8)]
-colnames(predictor_tot2) <- c("x", "y", "SystTracks_distance5_5", "Canopy1020", "Canopy2030", "Simpson")
+colnames(predictor_tot2) <- c("x", "y", "SystTracks_distance5_5", "Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")
 
 predictor_tot2$predicted_tot2<-predict(model_maxnet, predictor_tot2,type=c("exponential"), clamp=FALSE)
 
@@ -5375,11 +5526,16 @@ pres_bgGDT <- na.omit(pres_bgGDT)
 
 Pres_data<-as.vector(pres_bgGDT$presence)
 
-predictor_maxnet<-pres_bgGDT[, c( "Canopy1020", "Canopy2030",  "Simpson")]
+predictor_maxnet<-pres_bgGDT[, c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")]
 
+FORMULA<-maxnet.formula(Pres_data, predictor_maxnet, classes = "q")
 
 model_maxnet<-maxnet(Pres_data, predictor_maxnet, f = maxnet.formula(Pres_data, predictor_maxnet, classes = "q"), regmult = 1)
 
+var_names<-c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")
+
+out_M_4 <- iid(data=pres_bgGDT, V=10, var_names)
+
 #use it to predict results
 predictor<-rasterToPoints(new_env_raster_R)
 predictor<-na.omit(predictor)
@@ -5427,11 +5583,16 @@ pres_bgGDT <- na.omit(pres_bgGDT)
 
 Pres_data<-as.vector(pres_bgGDT$presence)
 
-predictor_maxnet<-pres_bgGDT[, c( "Canopy1020", "Canopy2030", "Simpson")]
+predictor_maxnet<-pres_bgGDT[, c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")]
 
+FORMULA<-maxnet.formula(Pres_data, predictor_maxnet, classes = "q")
 
 model_maxnet<-maxnet(Pres_data, predictor_maxnet, f = maxnet.formula(Pres_data, predictor_maxnet, classes = "q"), regmult = 1)
 
+var_names<-c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")
+
+out_M_5 <- iid(data=pres_bgGDT, V=10, var_names)
+
 #use it to predict results
 predictor<-rasterToPoints(new_env_raster_R)
 predictor<-na.omit(predictor)
@@ -5476,17 +5637,23 @@ pres_bgGDT <- na.omit(pres_bgGDT)
 
 Pres_data<-as.vector(pres_bgGDT$presence)
 
-predictor_maxnet<-pres_bgGDT[, c( "Canopy1020", "Canopy2030",  "Simpson", "BiaisTracks_distance5_5")]
+predictor_maxnet<-pres_bgGDT[, c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd", "BiaisTracks_distance5_5")]
+
+FORMULA<-maxnet.formula(Pres_data, predictor_maxnet, classes = "q")
 
 model_maxnet<-maxnet(Pres_data, predictor_maxnet, f = maxnet.formula(Pres_data, predictor_maxnet, classes = "q"), regmult = 1)
 
+var_names<-c("Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd", "BiaisTracks_distance5_5")
+
+out_F_6 <- iid2(data=pres_bgGDT, V=10, var_names)
+
 predictor<-rasterToPoints(new_env_rasterbias3)
 predictor<-na.omit(predictor)
 predictor<-as.data.frame(predictor)
 
 predictor$predicted<-predict(model_maxnet, predictor,type=c("exponential"), clamp=FALSE)
 
-colnames(predictor_tot2) <- c("x", "y", "BiaisTracks_distance5_5", "Canopy1020", "Canopy2030", "Simpson")
+colnames(predictor_tot2) <- c("x", "y", "BiaisTracks_distance5_5", "Canopy1020","penetrationratio0205","Simpson", "Q25", "Canopy2030sd","penetrationratio0205sd")
 
 predictor_tot2$predicted_tot2<-predict(model_maxnet, predictor_tot2,type=c("exponential"), clamp=FALSE)