From a37bca88f2b903a788b71dbc288056a90bfe7ffc Mon Sep 17 00:00:00 2001
From: "anouk.glad" <anouk.glad@irstea.fr>
Date: Tue, 18 Dec 2018 17:48:25 +0100
Subject: [PATCH] Revision

---
 Data/Regular_transect.dbf  | Bin 811 -> 625 bytes
 Data/Regular_transect.shp  | Bin 2036 -> 1508 bytes
 Data/Regular_transect.shx  | Bin 276 -> 228 bytes
 Data/transect_distance.tif | Bin 22765 -> 25574 bytes
 R_code.Rmd                 | 223 ++++++++++++++++++++++++++++++++-----
 5 files changed, 195 insertions(+), 28 deletions(-)

diff --git a/Data/Regular_transect.dbf b/Data/Regular_transect.dbf
index 85a249f0f8f55f985ea0fb4b25f4c4feace8ac74..96a8f3fef828e4741aefe82c77d3907df1f34d63 100644
GIT binary patch
delta 14
VcmZ3@_K}5^Ii6igU?Xb|6964O17iRH

literal 811
zcmaLTy$%6U6h+~aC?pcCLSsLIbN%<zXqZqaCknL*Ucj3pVnt-fEpEP<eXG@Jy2C+4
z++*~-f6gIdUGZ=8QjgC!>fLr<zaQDXI^IrAdpTd*uzO#e{3{h^Cry~0G-Y<ujM+(Z
zW+yF}owQ_j(u&zB$Di3rNi`=Wt(=s!a#GUDNl7awC9Ry4v~p6?%1KEpCl#tWsZh;H
Ig=*h=0KHshd;kCd

diff --git a/Data/Regular_transect.shp b/Data/Regular_transect.shp
index 74a60b2eacd9799e03c1fd9e181082cd00a2cf2b..8355e73eb862abc6d9933cd827a8929e8ae60e72 100644
GIT binary patch
literal 1508
zcmZ|P-Dv_*5QgF7s_{1_VgtdO1WXg!g}0VsDOO+&LK;~J2#Oa%KoCOqn0eP3_!L23
zXMEfRcT;-rC;2GO`({ejGkv{nx2tkC&uf?Qa-MrJ{=D^Zvb?=Y;{IWydi(Ln!p>G+
zEbM#gcm1xd*6#W{9wV(>W}FXg=aUb$b~pRz(#mCqOPkyB?(Ps>TDi<{X>;4&El<&<
zmCFp5HuuT9#W}jPa+%@M=61Y0{fREETxPhmxn1urF43ix%M6z`x98o}HM+ENnc>pr
z=4-G1pRIqRODmTdE^Y3>yPL-UwaPB7TxM#wGKb#XwW3Qaml-Z?$NTKvvJ+iexy*2B
nb4T7S2GOOJ%M6z`_r<%@QFLkLGQ*|Kef2K)SaxaUGE=*MdaPz~

literal 2036
zcmZ|QKQBXJ7>4oJX{ojx{e!`1A@~yg3f;<PG>XaSqZo}wqtR%ziA0(p5)q*yoKw$r
zJWq0Kn$+cf{d$^co6<%r$y;qQuBBAm)86+{etNq-D385&M}x9&dOOp8?2o38lKA`O
zO!4F9e-?IWl?OBXJN3N&tR1b~$z8lgTDfez9@<_{d8w5<ypJxeTsB<V+`4zq#?hsf
z%Z5vvJL}!-A-c43*>Gud=e)c76kS@mY`C<!4ey>mN0(ME8!m0`ymwnK(WRBkhD)2<
z^zQy^bZO<X;nL=o@4f2(?BXrDv~t;SX>;4&?R`X-RxTSZZSI114?d$yE0+zIHh0mx
zmtWDPmCJ@po11%gE&YCGah6stn>t&VOWr-qqDw264VSj-?RfX99$i|wY`C<!UGJ_p
zqDw264VN~z=iTFGbZO<X;nL>zy?fn`F0EWPT-w}a?{4MMrIpKuOPjmm-IGpqY2~uv
S(&nyum+xH8(#mC1XMX_CTEn>j

diff --git a/Data/Regular_transect.shx b/Data/Regular_transect.shx
index a0654096e4123f85dd0bcd9fead31b3c48893988..7b1fd0c16886af138b2489c76aa7f5ab69376d82 100644
GIT binary patch
delta 41
scmbQj^n`JOw0zMEW(EczVmPX>VMEd*ZAT!J;fG+PBZ%)LIMHeq00mhMO8@`>

delta 90
zcmaFDIE86~w0zeKW(EczVu=1<!nf$5wj+?q@Ix@t5yZbGG|_660RscmC!nAPP@No<
P_JY!7P<j!RJ_n)!16L8+

diff --git a/Data/transect_distance.tif b/Data/transect_distance.tif
index 83366a1faf1ae4cd7be11d7b633c6be442ff6a11..b1908305baacad9d6daea37a2857157caebdd134 100644
GIT binary patch
literal 25574
zcmeHPX;4#F6uwyqD2oLOh!`7D6ht5)5C|klG#Dd<En!gs4J1`$lO<vWDo7Pku_z4I
z1+^92Q4~Zhg1EF4sVE}if{IQR83a=mt;3+bdBho<PV2w^I&bFXeD}Wly)F05d+*+J
zKCyTP(nAQDB1EDf5=K<~L(<|vTqbL=H^vk#*1?#j#Tt$=iHfYSjlu=9QT66pOu^eM
zv~6BuJaJST{wORR)!r!A0BLw{$Ox_d&H&LcHo_Gr@+c?*V=s&=^9;}!O&g3)=8@52
zO+6K%>@qSk!RK-XuQ?b`!1Y~tKboZwxBWlX!}xfm_*@+D?XA~GNHPwg`z8n_m?Jc7
ziBJI@A^K#5DyAaDw?*ixJwlOA2=zN7l*K{F&;#@Nm@mM5A?DA+d|%8D!2Dp$55@e2
zm>-Gxa)fvyp}<d~Op`0*X>zZz3<gh}rc7WY$`h1)YpF~i6H8^{Inu=vvA<XnBC%$q
z$`e+^<84f9FDB1P)A+wwg7AN^oZVbM&+`%Ze{v`m)7g>5a%H=_a5TT4-GNjl6#9gI
zvI`E=k>k!~b6IRRraRN^<1TnkpWXrFz;fhpKV4ju_H3s6NJl_i?-B6PYODXn)OOqB
ze*&TBrV9wG!6d?C%&^mpVBE!tIx_WSp2n**u0=v0){yj0=7r({SzD$NcFPn!ZF!7d
zWgZFlaDBq+B16LBjRqrpver@vtB0tBwVOgS(D2Aem#`FMKv)toBrNt9!cx&V!fL#B
z)Z#7j9x3UK*z9@Q@?63;M%WP78R$S<ulE62xZu#czG8Y=inbj>szgoqQe=weNE!YZ
zG58Eou#~~#x-)IXAtKRC0@`H?<!Pz}yit`Hmzb2EnCTNAAG;zoDONE;aA#%)ro<-T
z-irAtN@ZetY|OIEOq7infow5Zm%PGq$@P__jTRJ!#Ur)z?|!J!&L^At!}X`Gy17kd
z#<kVau4aWoMP)B_nyXI9#j-DMS054CuAi$&Zl1VNG;6P*r)UT7mJPc^@Z<(}-3+Bc
z*=5_8V=fy>JhwDGGd@CJsFLlU|M=!v9geQzsBDQej@o*D_OWYys!b|$skHibuR6oU
z`Ym~9<iJhevbDJ>Igh$Sx^x=s<YiIi?`Ccpn`F}zJ>Z^in99mM9DabypA<SSKDQ(J
zY;S+uqID;hIizpaIYe5vRyo5`KhZ?6IVbw$L`R*)89%l!c%il*c(^^My`Qo0P{Y^a
z6RpvzUH$$qUzML4s2iql4?9!({7ub-a}A5qyBk*3)IZPa?7Yw^Kli@X>zeSOOCe?3
z2mbNLv~FZ%^2cMg+P?5@V}O45yZY9?=MSrr{4Dh^UYknm@phDc(7jZ#=>_$?laOWP
zRrz|dvEBFE9Phar$4ku42l%RoYTIUxu{hCYkUQ?z<GZaU`UJMWc_nBYY8pR#p1U_^
z+Wyru#~V3+DX{kNrW~HssBWn;Y`o!OKX$g$=INFVMJ`wRnpiiRt(qFC8=uWLm9@JF
zw^sdC>uDle)#p)E7vII*TKwnqfE$X!LI2>03vK~bid%?V?PKBTduphHD)24+O5|Vd
zqVx+2D)*dn-Lya_=tj!xdy+OYVW;G`)5YOe1KD(&*`aBk+ZP3-RY!W<a-~Y=g;bSW
zNrPXv><!`7_+4)GPaWRV8(KmyH*PB~KOGg)F-UuP*W4kWcWt|&%I!|j&7jog>L(W*
zuDB@+=HGv#YM$EFu`Y-GFopgqSNvARw|S6O@a+AA>%(zLEhQDR;%ZgHZHeK*9j?ox
zYEEaw$Ane1XGN~^xwql5rCG~1`Jf|LdiK=rqxLPVf}}!aZG&j@Q_lSt`5o+pq2$)_
zG<i<UYI^-o-8=bWW{*0`cZ1?xh-vAw+YteQJp+!qF;7-cJA0_FW0F<aou7U5cii<@
z-<-TWeRxN9f8Fu&V-HO`1|vdL4~pt~Gp@#*>AgL8@Oa101B>^??%+I9MQ@`mP%Q2G
zwqoDFuBTU{)#Ga{E0c%Jcjd<pdG#Dt-VK|4DL$>D?w6>+*7kG#8BO#xNtt_st~}z^
zTisHZu1P-=&OPwt+m-Km>3Jz0O<5I6&w;Re$#rM;)_l`7tEhG`YH8Jlp{&l+FP4<<
zQQ6JU_d8bTOD2=SAHg5NAHg5t`3TQPcs|1O5#CqfeHGqU;rkK3AL08Ez8_&c2gY+?
zJO{>CVSE+FS7Ce=<~v}%1Liwm{t@OMVg3>3A0Zw9;sGEY0OE5XJ_q7+AU+4;<se=T
z;^iQI72;PReih<ZA)XxK$swK`@;e~E1M)i{zXS5EAm0k|tsvhD@*g4p5%M1)9~ttI
zAs-p?kwG5-^Z`I00Q5RQuLJZtK(7Pzb3i`_^m9PZ3iPZ%&kFRcKwl2@<v?E!^d3R)
k5%eBG?-BG@L4OtWS3wUM^pHUh8T61rpB(hbTWt0I20(&^EdT%j

literal 22765
zcmeHPX;4#F6g~-&03wb?#Rb|Zt{5PgkPtSBki;kiVh9AmO4L}X5CQ}pHC0^3s93Cu
zXk99!C~X-Ptbikn3s8kYK!I9?mQpLNgJ>(&>L6`z9%%=s)8$vEKhB$Z?|gSTZ%OW%
z_wGICd~&%v*Z}|wfMx?Uf*9mLv*Jinrdx3^!8TSrj9^<Uw&;waF<=srZ;`^-A@$>|
z*oL$@Tl+jE*kwqcz#&*Qq`%WvN3iHUrvq%+<_NX~JCTZWW*dknIEdgPvm@AB`Vd@f
zro&Q8Jp&*!pAMtRz0{L^Gr<!{eFr&@Wt&0z{uldVa=l~ZUTT;CJq`e)od86m020Ok
zn4JO6lACIC1$e~<aGMISU<N>{JHR0ifXiHfZXUo0ApqYGATE&n2EhOyhX6E70e+SP
zyb%r{A%9Sk3SfUUK!7Y%92TxKYFB8D+8_ruJ3wyKCA0O~WL@ACRg^eNu8NX}sFsGy
z74mR(_!M@!HrbFw+PG7KxB*_4&VQ#R*8EA!+vi{CNyLiRE+yo9&f#%|0@2{{A1|zm
z3JsM+zjg{f-;>7^`igwK$<7n~eu@CE*M9-<@b(mlUf=k-**vZ<83mS26ikUf<Y4Kq
zju*e4Ew;2hV#QL@f@ze+j3q`hh{+JA>gbGOv&E_`rbUBaH__~h&C#Slx0Y=wneAq0
zE!*1_nQ3H%52qv(Mo^M8Iu7d7t<r{)cEF(IEH_(6iyj@Al(YdyO448iC5gRIk^v(r
zX|dWNTd*v6uw*wVGv`~&Qc7NlpGG}rqzCo9!Kb4%;*VbT<hbUiS^EKG=q<x*g-I5o
zVk-z@2eM_0RBUgNFL$O~Et7F5ILj2OHD)A}#ti*xeac$BNs^S5Xh=^<TrsHN%QZ!&
zB_@;6O7zooI{n(j<*Q64$OLxeoH{y_o-ih_*+84_%yg;l+UvVp;p(`pCNEN>h07BY
z#-CJv-%J-Vch}|@YmS!2GPCL@D@6fL^6`6R%lqVv#tGg#<hne2lhbCq<Kph3SrsSJ
znWt|ToL?N^kQGv>c2BQs-Ev%Tta?_0eb_HC^A1Lz9CK)DcD2p>?fIqtoNG;*?dq<Q
z^0D=)8>{a1#vBP%F%6~CADkTq=5!`sT~Ksqm(v1WdhIpuJcsNMuiW|vv;7>T%gQ&k
ze4YM0>BC(+&A0EkkMEzJ;=k)g5BG5I#(i}=_pe-2;px=fd*;?#pOnQv+frF`$>3i8
z<y!Zbi(2pdxvux=EG+7LSZeyENx1w{^HY<{&Fg!VkK2|7xJE?x->rSIrlP7dxBvTk
zzeZj|qaZ8m`K!6-=Gj)xGz^b^G?&v-Mq?|Ej=xaVNdK(;N^{nWzKgx>mmMn<>7A?~
zZ<|Zs29J#AcZC#(PG;12g~g2)DDz`R1r@2PU5cBg1~Md(g_Uo7cAm47K3|z*+kdTr
zv(EK;<f{)@H|8)ly*an1#?`&fa6jpuJZ;;aN!Obyg!G0~&ViX_`(7RJsXO1G@~HE0
zk@0JvbS(Fmrkk3)BrV3X1@jE1kB8MOFTCdv#?!I53RXv-kCc^E!jt4oP4#bxIyOvE
z$jbSV*2?ui@$%X;w}U~|EdwKiTYoYYL`3BHs3y0iJd2AvXWaIX6P>o<l(_W*r({Xp
z?GjOfxKu5cs}&EDWwHE$#Fdfbd}5`ya@&;e#N?K2T@^C+=E1bsPd)Z5S}YuMMz1N&
z=rpSD_-ycud0xA8$&%dO-jqhcfwIi6%GTH~RA!uU^Xn3J^m(=m)df2QFVycEBBCm_
zD?isXPTT&?!EBGM_gT|rt2=yG#+o%oN#XoO=cSXvB3fIYyxl%Q5SXds3)ba!-`$w1
z`A&V>oPBe`z@5C7wN)F855K6~xi|jq*&|2pzZe&CV6{N{?%7ontF9hST70^$;L*Kt
zr%tS1FDPH%)4#^~L`R?PvVr(=Nmq$J&Qn))D|X|;D=EBVJf-$tTk0ndCh1R>KMqi6
zYr9Qe9@U+k(R4HuG!rxvG!uNA;M)Y>Ciu~XA6@v-g>w^}o8a68r|>w1$0<Cnbm2-D
zu5{soFfIt=f-w3`&~Ji%6ZD&)dlB7>=w3uGJbK~L3y;An7_5T9Dj3s+F<ltbg`tla
z`iP;A7$J-i!Wbcpf%F(ikAd`<Z-V(Im~VplCYY9mX<3+-h1rXky@=V1nB<E|zL?~T
zx$u|^kGb$j8G)1$NEw05D#)yY%qmDggakxLK!iM9$kT;9UC7ggbaqH*hjezx`uP7~
JeXKuB{|zlG>kt3{

diff --git a/R_code.Rmd b/R_code.Rmd
index fec86ad..7e6c9a2 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)
 
-- 
GitLab