lunes, 23 de julio de 2018

Antipredator behavior with R (or why wildebeest should stay together)



Hi there!

 In 2010, when I was studying my Biology Degree at Universidad Complutense in Madrid, I fell in love with a documentary miniseries called Great Migrations (National Geographic). Their episodes talk about awesome migrations of animals around the globe. One of these chapters is “Great Migrations: Science Of Migrations”. It shows how scientists study the patterns and processes of animal migration. One of these researchers is John Fryxell, from Gelph University in Canada. John explains how mathematical models can help to understand movement patterns. Through simulations, he shows how wildebeests maintaining a clustered movement pattern can more effectively avoid to be hunted by a virtual lion (here is the Spanish version in YouTube). I have to confess that each time I see those images, I think that I would be completely fulfilled with this kind of jobs.





I modified a little bit my previous “go” function to reproduce the experiment shown by John. This is a small tribute to this documentary series and to all those people who inspire and motivate us and arouse our curiosity, as John.

You can find the code to perform these simulations with R below. To perform the gif of this post, you can follow this procedures. And if you want to be updated with my posts, find me on my new Twitter account @javi_ferlop

That is all, enjoy!




 # First, we prepare packages and raster layers to be used later. We also load  
 # the species class to define the hunter inital location and his environmental  
 # optimum (grass cover in this case), as in previous posts.  
   
 library (raster)  
 library(dismo)  
   
 gc <- raster(nrows = 100, ncols = 100, xmn = 0, xmx = 100, ymn = 0, ymx = 100)  
 gc[] <- runif(10000, -80, 180)  
 gc <- focal(gc, w = matrix(1, 5, 5), mean)  
 gc <- focal(gc, w = matrix(1, 5, 5), mean)  
   
 species <- setClass("species", slots = c(x = "numeric", y = "numeric",   
                      opt = "numeric"))  
   
 # Then we will load the hunt function. This is pretty similar to my previous  
 # "go" function. However, I added a "preys" attribute, a set of locations  
 # of available preys. After each movement of the hunter, the function checks  
 # if the hunter is in the same cell as any prey. If yes, it records this  
 # location and delete this prey from the available preys. The output is a 4  
 # column dataframe: the first two are the location of the hunter at each  
 # iteration. The last two show if there is any capture at this iteration.   
   
 hunt <- function (sp, env, n, sigma, theta_x, alpha_x, theta_y, alpha_y, preys) {  
  track <- data.frame() #create an empty data.frame to store each step coordinates  
  track[1,1] <- sp@x  #the first position will be the initial position of the specimen defined by the species class  
  track[1,2] <- sp@y  #the first position will be the initial position of the specimen defined by the species class  
  hunted <- data.frame()#create an empty data.frame to store hunted preys  
    
    
  for (step in 2:n) {  
     
   #update preys locations   
     
   prey_cells <- cellFromXY(env, as.matrix(preys + step))  
     
   #First, the function searches in the adjacent cells what have the value that best fits with the species requirements  
   lon_candidate<--9999  
   lat_candidate<--9999  
   #randomization/atraction routine  
   while ( is.na(extract(env, matrix(c(lon_candidate,lat_candidate),1,2)))) {  
    lon_candidate <- track[step-1,1]+ (sigma * rnorm(1)) + (alpha_x * ( theta_x - track[step-1,1]))  
    lat_candidate <- track[step-1,2]+ (sigma * rnorm(1)) + (alpha_y * ( theta_y - track[step-1,2]))  
   }  
     
   neig <- adjacent(env,   
            cellFromXY(env, matrix(c(lon_candidate,  
                        lat_candidate), 1,2)),   
            directions=8, pairs=FALSE )  
     
   options <- data.frame()  
   for (i in 1:length(neig)){  
    options[i,1]<-neig[i]  
    options[i,2]<- sp@opt - env[neig[i]]  
   }  
   option <- c(na.omit(options[abs(options$V2) == min(abs(na.omit(options$V2))), 1 ]),   
         na.omit(options[abs(options$V2) == min(abs(na.omit(options$V2))), 1 ]))  
     
   new_cell <- sample(option,1)  
   new_coords <- xyFromCell(env,new_cell)  
     
   track[step,1] <- new_coords[1]  
   track[step,2] <- new_coords[2]  
     
   if (new_cell %in% prey_cells) {  
    hunted[step,1] <- new_coords[1]  
    hunted[step,2] <- new_coords[2]  
    preys <- preys[-match(new_cell,prey_cells),]  
   }  
   else {  
    hunted[step,1] <- NA  
    hunted[step,2] <- NA  
   }  
  }  
  return(cbind(track,hunted))  
 }  
   
 # Later we define our predator, in this case, a lion! The inital possition will  
 # be at the middle of the space (50,50). We will also select the environmental  
 # optimum of our lion, 60% of grass cover (our lion hates bare ground xD)  
   
 lion <- species(x= 50, y =50, opt= 60)  
   
 # Now, we will set the locations of a bunch of wildebeests using "randomPoints"   
 # from "dismo" package. As our wildebeests will be running around, we will   
 # initially locate them at the bottom left.  
   
 r <- raster(nrows=500, ncols=500, xmn=-500, xmx=0, ymn=-500,ymx=0)  
 wb <- randomPoints(r, 10000)  
   
 # Here we have our wildebeests ready to run!  
   
 plot(wb)  
   
 # Let start the hunting!   
   
 hunting <- hunt(lion, gc, 600, 2, 0, 0, 0, 0, wb)  
   
 # In "hunting" object is stored each location of our lion during the 600   
 # iterations, as well as each time she catches a wildebeest. To see only the   
 # catches, you can use:  
   
 na.omit(hunting)  
   
 # Pretty awesome, right? But, what happens if we group our wildebeests in   
 # clusters? To do it, we will define another small function that creates   
 # clusters of points:  
   
 clus <- function(n,m1,m2){  
  df<-data.frame()  
  for (i in 1:n){  
     
   a <- rnorm(100,m1[i],1)  
   b <- rnorm(100,m2[i],1)  
   df <- rbind(df,data.frame(a,b))  
  }  
    
  return(df)  
 }  
   
 wb_clustered <- clus(100, runif(100,-500,0),runif(100,-500,0))  
   
 # Look at the wildebeests locations now!  
   
 plot(wb_clustered)  
   
 # Ok, let's go again!  
   
 hunting_clustered <- hunt(lion, gc, 600, 2, 0, 0, 0, 0, wb_clustered)  
   
 # Let's chek how many catches he have now:  
   
 na.omit(hunting_clustered)  
   
 # And comparing both strategies...  
   
 nrow(na.omit(hunting))  
 nrow(na.omit(hunting_clustered))  
   
 # ...looks like wildebeests should stay together!!  
 # You can run these simulations n times to see if the results are robust.  
   
 # Finally we can print each iteration to make a cool gif ;)  
   
 for (i in 1:600) {  
  id <- sprintf("%03d", i)  
  png(paste("lion1_",id,".png", sep=""), width=1300, height=690, units="px",  
    pointsize = 16)  
    
  par(mfrow=c(1,2))  
    
  plot(gc, xlim = c(6,94), ylim = c(6,94), main="random", legend=FALSE)  
  points(hunting[i,1:2], cex=1.5, col="red", pch = 0, lwd =3)  
  points(hunting[1:i,3:4], cex=1, col="blue", pch = 3, lwd =3)  
  points (wb +i, pch = 20)  
    
  plot(gc, xlim = c(6,94), ylim = c(6,94), main="clustered", legend=FALSE)  
  points(hunting_clustered[i,1:2], cex=1.5, col="red", pch = 0, lwd =3)  
  points(hunting_clustered[1:i,3:4], cex=1, col="blue", pch = 3, lwd =3)  
  points (wb_clustered +i, pch = 20)  
    
  dev.off()  
 }  
   
   

3 comentarios:

  1. Hey, really enjoyed the post, when i tried to run the code the gif part didn't work, my R didn't recognise 'gnu' or 'gnu_clustered' I assumed it was supposed to be 'wb' and 'wb_clustered'? And when I tried to find the gif it just saved them down as 'png' files. thanks for any help, Mark

    ResponderEliminar
    Respuestas
    1. Hey Mark! Yes, you are right, is a mistake! it should be wb and wb_clustered, sorry.

      And you are right, with this code you are just getting the png files... to build the GIF you should use ImageMagick or other software (there is an R package to do this also called magick). Check this other post and its comments:

      http://allthiswasfield.blogspot.com/2017/12/p-margin-bottom-0.html

      Eliminar
  2. I am a novice R user, so please excuse if the answer is very simple :)

    ResponderEliminar