Resource selection function

#Import Fish database
load("Fish72.Rdata")
fish<-Fish72; rm(Fish72)

#Import the raster layer of the lake to plot the fish locations on the map
load("Bathymetry_Ledoux.Rdata")

#Look at the fish data frame structure
head(fish)

#plot the raster
image(ledoux_depth,col=palette,ylim=c(780,1150),main="The 596 locations of the Fish recorded in 4h")

#add the fish locations on the map
points(fish$X,fish$Y,pch=16,cex=0.7)

#plot the raster
image(ledoux_depth,col=palette,ylim=c(780,1150),main="Step comparison")

#add the observed steps on the map (length plus angle)
lines(fish$X[c(1:300,400,500,590)],fish$Y[c(1:300,400,500,590)],lty=4)
lines(fish$X[c(300,400)],fish$Y[c(300,400)],lwd=2)
lines(fish$X[c(400,500)],fish$Y[c(400,500)],lwd=2)


#add simulated steps on the map
lines(c(fish$X[300],800),c(fish$Y[300],1050),col="red",lwd=1.5)
lines(c(fish$X[300],840),c(fish$Y[300],970),col="red",lwd=1.5)
lines(c(fish$X[300],950),c(fish$Y[300],1025),col="red",lwd=1.5)

lines(c(fish$X[400],750),c(fish$Y[400],1050),col="forestgreen",lwd=1.5)
lines(c(fish$X[400],700),c(fish$Y[400],1070),col="forestgreen",lwd=1.5)
lines(c(fish$X[400],790),c(fish$Y[400],1055),col="forestgreen",lwd=1.5)

# The function as.ltraj enable us to transform a matrix with time and locations into a succession of steps, with length of the step, the angle between a step and the previous one, the duration of the step etc.

#NB: when we have several fish or periods sampled, this function allows to deal with it by specifying the name of each individuals and the burst (i.e. a succesion of relocation for a specific animal at a specific date). It will result in a list containing all the calculated trajectories.

track_fish=as.ltraj(fish[,c("X","Y")],fish$posix_est,id=fish$Name,
                     burst=fish$burst,slsp="missing",proj4string = crs("+proj=aeqd +lat_0=46.802927 +lon_0=-73.277817 +x_0=1000 +y_0=1000 +datum=WGS84 +units=m +no_defs"))

#lets see what we have. Here, we have one fish tracked during four hours, we have just one burst so one element in the list.  
head(track_fish[[1]])

#The function redisltraj enables us to homogenize each step to one-minute steps by specifying the step duration in seconds. It will just recalculate the locations of each steps by interpolation.

fixed_track <- redisltraj(track_fish, 60, type="time")

head(fixed_track[[1]])

#plot the background bathymetry map, with a little zoom in the area it used
image(ledoux_depth,col=palette,ylim=c(800,1100),xlim=c(700,1200),main="Raw location vs. Fixed track")

#add the locations of the fish
points(fish$X,fish$Y,pch=16,cex=0.7)

#add the calculated trajectory of the tracked fish
lines(fixed_track[[1]]$x,fixed_track[[1]]$y,col="salmon",lwd=1.5)

# Set the minimum depth the fish is able to go
lim_depth=min(fish$Depth,na.rm=T)
print(lim_depth)

# Create a raster containing only areas of the lake with minimum depth.

# To do so, I inelegantly transform the bathymetry raster into a data frame.  
map_temp=as.data.frame(ledoux_depth, row.names=NULL, optional=FALSE, xy=TRUE,
                     na.rm=FALSE, long=FALSE)

# Then I select only the locations of the lake where the lake is deep enough to welcome the observed fish
map_temp=map_temp[which(map_temp$layer>=lim_depth),]

# I create an empty raster with the dimension of the lake original raster
 e <- extent(c(min(map_temp$x),max(map_temp$x),min(map_temp$y),max(map_temp$y)))
  ncellx=(e[2]-e[1])
  ncelly=(e[4]-e[3])
  empty_r<- raster(e, ncol=ncellx, nrow=ncelly)

# Finally I integrate the selected locations into the empty raster
r <- rasterize(map_temp[, 1:2], empty_r, map_temp[,3], fun=mean)

# I built a constraint function which must have SpatialPixelsDataFrame object so we convert the raster
map <- as(r, "SpatialPixelsDataFrame")
crs(map)<-NA

# Here the area in the lake we will let our simulated trajectories go, depending on the depth distribution of the observed fish
image(ledoux_depth,col="black",ylim=c(780,1150),xlim=c(300,1600))
par(new=T)
image(r,ylim=c(780,1150),col=palette,xlim=c(300,1600))

#constraint function which suppress the location if coordinates of a simulated trajectory (x) fell outside of the constrained area (par). This constraint function is one of the arguments integrated to the function simulating trajectories.

consfun <- function(x, par){
  coordinates(x) <- x[,1:2]
  ov <- over(x, geometry(par))
  return(all(!is.na(ov)))
}

  mo <- NMs.randomCRW(na.omit(fixed_track), rangles=TRUE, rdist=TRUE,
                      constraint.func=consfun,
                      constraint.par = map, nrep=9,fixedStart=TRUE)

 Simulated_tracks=testNM(mo)

# Transform the list of trajectories into a data frame, adding a column "no" for the trajectory id; 1 being the observed one, and 2 to 10, the 9 simulations

no=1
RSF_database=cbind(fixed_track[[1]][,c(1,2)],no)

for (i in 1:9){
  no=i+1
  RSF_database=rbind(RSF_database,cbind(Simulated_tracks[[1]][[i]][,c(1,2)],no))
}

# A color vector to have the observed trajectory in color and simulated ones in grey.
couleurs=c("salmon",gray.colors(9, start = 0.3, end = 0.9, gamma = 2.2, rev = FALSE))

# plot the map
image(ledoux_depth,col=palette,ylim=c(780,1150),main="Simulated and observed trajectories")

#add all trajectory lines
for (i in 10:1){
  temp=subset(RSF_database,no==i)
  lines(temp$x,temp$y,col=couleurs[i])
}

distri_depth=fish$Depth
hist(distri_depth,main="",xlab="Depth (m)")

##Lets pass quickly on the recalculation of real depths for the observed trajectory

fixed_track[[1]]$depth=NA

#for loop to assign a depth to each relocation in the observed trajectory
for (i in 1:nrow(fixed_track[[1]])){

  # the time at the relocation i
  time_stamp=as.numeric(substr(fixed_track[[1]]$date[i],12,13))*3600+
                  as.numeric(substr(fixed_track[[1]]$date[i],15,16))*60+
                  as.numeric(substr(fixed_track[[1]]$date[i],18,19))

  # Gathering time for all depth measurements
  time_table=as.numeric(substr(fish[,"posix_est"],12,13))*3600+
                  as.numeric(substr(fish[,"posix_est"],15,16))*60+
                  as.numeric(substr(fish[,"posix_est"],18,19))

  # Selecting the depth where the time of measurement was the closest to the time of the relocation i              
  ind_timing1=which(abs(time_stamp-time_table)==min(abs(time_stamp-time_table),na.rm=T))
          print(c(i,length(ind_timing1)))

  # If both previous and next closest depth measurement are equally distant, we average the depth between the two to have the depth at relocation i
   if(length(ind_timing1)==2){
                fixed_track[[1]]$depth[i]=mean(fish[,"Depth"][ind_timing1])
   }

  # If we have only one closest depth measurement
   if(length(ind_timing1)==1){

  # If the closest depth measurement is after the time of relocation i, we take this measurement and the one before and then make a linear interpolation to retrieve the depth at relocation i
    if(time_stamp-time_table[ind_timing1]<0){
      ind_timing2=ind_timing1-1
      time=time_table[ind_timing1]-time_table[ind_timing2]
      depth_var=fish[ind_timing1,"Depth"]-fish[ind_timing2,"Depth"]
      depth_inter=(time_stamp-time_table[ind_timing2])*(depth_var/time)+fish[ind_timing2,"Depth"]
      fixed_track[[1]]$depth[i]=depth_inter
    }

  # If the closest depth measurement is before the time of relocation i, we take this measurement and the one after and then make a linear interpolation to retrieve the depth at relocation i
    if(time_stamp-time_table[ind_timing1]>0){
      ind_timing2=ind_timing1+1
      time=time_table[ind_timing2]-time_table[ind_timing1]
      depth_var=fish[ind_timing2,"Depth"]-fish[ind_timing1,"Depth"]
      depth_inter=(time_stamp-time_table[ind_timing1])*(depth_var/time)+fish[ind_timing1,"Depth"]
      fixed_track[[1]]$depth[i]=depth_inter

    }
  # If the closest depth measurement is at the exact time of relocation i, we directly take this measurement to assign the depth at relocation i
    if(time_stamp-time_table[ind_timing1]==0){
      fixed_track[[1]]$depth[i]=fish[ind_timing1,"Depth"]
    }
  }
}

  # Create a column to inform the depth of the water column at each simulated location
RSF_database$max_depth=NA

  # Use the extract function of the raster package to retrieve the water column depth for each simulated locations
RSF_database$max_depth=extract(ledoux_depth,RSF_database[,c(1:2)])

  # Empty vector
depth_random=NA
  #for loop where we paste each depth randomized sampling with the vector
  for (i in which(RSF_database$no!=1)){
    depth_random=c(depth_random,sample(distri_depth[which(distri_depth<=RSF_database$max_depth[i])],1))
  }
  #Then, in the database regrouping all observed and simulated trajectories, we add the depth column

RSF_database$depth=c(fixed_track[[1]]$depth,depth_random[-1])

#Lets see how our depth distribution behaved

ggplot(RSF_database,aes(x=depth))+
  geom_histogram(position="dodge",bins=15)+
  facet_wrap(~no)

  # Create a vector of depth breaks
depth_vector=seq(-0.5,16.5,1)

  # Load the zooplankton database
load("large_Z.Rdata")

  # Create the column with future zooplankton concentrations
RSF_database$LZ=NA

  # Create empty raster to create zooplankton raster to be extracted for each depth layer.
e <- extent(ledoux_depth)
  #Wished horizontal resolution : 5m per 5m
ncellx=(e[2]-e[1])/5
ncelly=(e[4]-e[3])/5
r<- raster(e, ncol=ncellx, nrow=ncelly)

  # For loop to assign zooplankton concentrations at both simulated and observed locations for each depth layer
for (prof in 1:(length(depth_vector)-1)){

  # Subset of the trajectories database with only locations in the right depth layer
  temp=RSF_database[which(RSF_database$depth>=depth_vector[prof]&RSF_database$depth<depth_vector[prof+1]),]

  # If we have any locations here
    if(nrow(temp)>0){

  # Select the zooplankton measurements for the corresponding depth layer
      tempZ=Large_zooplankton[which(Large_zooplankton$Depth==(depth_vector[prof]+0.5)),c("X","Y","Predicted")]

  # Transform it into a raster  
      temp_rast_LZ <- rasterize(tempZ[, 1:2], r, tempZ[,3], fun=mean)

  # Extract the zooplankton concentrations for all the locations in this depth layer
      RSF_database$LZ[which(RSF_database$depth>=depth_vector[prof]&RSF_database$depth<depth_vector[prof+1])]=extract(temp_rast_LZ,RSF_database[which(RSF_database$depth>=depth_vector[prof]&RSF_database$depth<depth_vector[prof+1]),c(1,2)])
    }
}

head(RSF_database)

RSF_database$weight=ifelse(RSF_database$no==1,1,5000)

  #To conditional regression, we need a TRUE / FALSE variable. Here, only the no 1 trajectory is TRUE since it is the observed one
RSF_database$case=RSF_database$no==1

  # Scale the variables
RSF_database$LZ=scale(RSF_database$LZ)
RSF_database$max_depth=scale(RSF_database$max_depth)

  # Write the regression model
RSF_model= glm(case ~ max_depth + LZ,data = RSF_database,weight=weight, family = binomial(link = "logit"))

summary (RSF_model)

RSS_zoo=exp(RSF_model$coefficients[3])
paste0("RSS = ",round(RSS_zoo,3))

  # Here a database containing regression estimates calculated 10 times with various number of simulated trajectories for our fish
N_Simul=read.csv("S4_base.csv")
N_Simul_Z=N_Simul[which(N_Simul$param=="Large zooplankton"),]


ggplot(N_Simul_Z, aes(y=coeff,group=N_random_tracks))+
  geom_boxplot(outlier.shape = NA)+
  theme_bw()+
  xlab("Number of simulated trajectories")+
  scale_x_continuous(breaks=c(seq(-0.4,0.4,0.8/6)),
                   labels=as.character(seq(0,300,50)))+
  geom_vline(xintercept=0,lty=2)+
  theme(axis.title.x = element_text(size=14),axis.text.x = element_text(size=14),axis.text.y = element_text(size=14))+
  theme(legend.position = "None")+
  ylab(paste0("RSF estimate"))

  # Import database where we calculated 500 times the regression based one simulated trajectory and 150 other simulated trajectory
H0_500=read.csv("HO_500.csv")

    # Import database where we calculated a various number of times the regression based one simulated trajectory and 150 other simulated trajectory
H0_data= read.csv("S5DB.csv",row.names=1)

ggplot(H0_500, aes(x=RSF_estimate))+
  geom_histogram(fill="gray80",colour="gray50",bins=18,aes(y=(..count../sum(..count..)) * 100))+
  theme_bw()+
  geom_vline(xintercept = quantile(H0_500$RSF_estimate,c(0.025,0.975)),lty=2)+
  xlab("RSF estimate")+
  theme(axis.title.y = element_text(size=14),axis.title.x = element_text(size=14),axis.text.x = element_text(size=14),axis.text.y = element_text(size=14))+
  theme(legend.position = "None")+
  ylab(paste0("Percentage of occurrence (%)"))

  ggplot(H0_data, aes(y=inf,x=N_H0_RSF))+
  geom_smooth(col="gray50")+


  geom_smooth(aes(y=sup,x=N_H0_RSF),col="gray50")+
  geom_point()+
  geom_point(aes(y=sup,x=N_H0_RSF))+
  theme_bw()+
  geom_vline(xintercept = 500,lty=2)+
  xlab("Number of RSF iterations")+
  theme(axis.title.y = element_text(size=14),axis.title.x = element_text(size=14),axis.text.x = element_text(size=14),axis.text.y = element_text(size=14))+
  theme(legend.position = "None")+
  ylab(paste0("Boundaries of the 5% confidence interval of RSF estimates"))