eWalk, an R function to simulate two-dimensional evolutionary walks

[This article was first published on long time ago..., and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.



Hey there! 

 During last months, I was thinking with my friend Javi Fajardo (see his work here!) a little bit more about Brownian / OU process and their applications in biology. Remember that an OU model describes a stochastic Markov process that consist of two parts: first, a Brownian motion (random walk) model that represents a completely random process of variable evolution, and a second term that represents the attraction of an optimum value of the variable. For this, OU process can be considered as a variation of Brownian motion model. In the last posts we were talking about these kinds of models used to drive trait character evolution, specifically an ecological trait, such as optimum temperature or precipitation.

 
 MycolleagueJavi and I were talkingabout the applications of this models in two dimensions, representing for example space coordinates (latitude/longitude). They have been used, for example, to simulate predator/prey movements, home range movements, etc. We developed an R function similar to the previous eMotionfunction in order implement this kind of movement models: eWalk(from evolutiveWalk…). eWalkis similar to other movement-simulation functions in adehabitatLT package (Calenge, 2006). However, you can use eWalk to simulate a variety of processes as Brownian motion, OU processes or Lèvy flights (see below).

Here you have the (commented!) codeto load the function into R console:

   
   
eWalk <- function(sigma,lon,lat, generations,alpha_x,theta_x,alpha_y,theta_y, color,levy=FALSE, plot=TRUE) {
  route=data.frame(1:generations,1:generations)                             # a data.frame to store the route
  names(route)<-c("lon","lat")
  generation= generations                                               # generation counter
  
  while (generation > 0){                                                   # stops the simulation when generation counter is 0
    if (levy == TRUE){                                                      # only if levy=TRUE: Levy Flight!
      rndm=abs(rcauchy(1))                                                  # process to select if I should jump or not
      if (rndm > 15){                                                       # Levy jump probability threshold
        l=15                                                                # length of the jump
        x= lon + (l * (sigma * rnorm(1))) + (alpha_x * ( theta_x - lon))    # Brownian or OU process for longitude
        y= lat + (l * (sigma * rnorm(1))) + (alpha_y * ( theta_y - lat))    # Brownian or OU process for latitude
        loc=cbind(x,y)
        route[generation,]<- loc                                            # store new location
        generation= generation - 1                                        # advance to next generation
        lon<- x                                                        # go to the new position!
        lat<- y
        
      }
      else{                                                                 # if Levy jump is not selected, l=1 i.e. proceed with  
        l=1                                                                 #normal Brownian motion or OU process
        x= lon + (l * (sigma * rnorm(1))) + (alpha_x * ( theta_x - lon))
        y= lat + (l * (sigma * rnorm(1))) + (alpha_y * ( theta_y - lat))
        loc=cbind(x,y)
        route[generation,]<- loc
        generation= generation - 1  
        lon<- x      
        lat<- y
      }
    }
    else{                                                                   # if levy=FALSE, l=1, only Brownian motion or OU
      l=1
      x= lon + (l * (sigma * rnorm(1))) + (alpha_x * ( theta_x - lon))
      y= lat + (l * (sigma * rnorm(1))) + (alpha_y * ( theta_y - lat))
      loc=cbind(x,y)
      route[generation,]<- loc
      generation= generation - 1  
      lon<- x      
      lat<- y
    }
  }
  
  if (plot == TRUE) {                                                       # if plot=TRUE, plto the route!
    lines(route, col= color, lwd=3)  
  }
  return(route)                                                             # print the data.frame with each location during the walk!
  
}
>


 The needed parameters to run the function are similar to those for eMotion, but taking into account that now you have two variables (latitude/longitude, for example) at the same time.




eWalk(sigma,initial_lon,initial_lat,number of generations, lon_alpha, lon_optimum, y_alpha, y_optimum, color, levy=FALSE, plot=TRUE)



In addition, I have added a complement to simulate Lèvy flights, a non-gaussian process similar to Brownian motion but allowing for jumps during the motion. Lévy flights pattern has been found in searchingmovements of predators (Bartumeus et al., 2005), or have been used to simulatepunctuated equilibrium in traitsevolution (Gould and Eldredge, 1977; Landis et al., 2013). To perform Lèvy flight, eWalk selects a random number from Cauchy distribution and if it is higher than a threshold (15 by default), Brownian motion term is multiplied by a constant representing a jump (also 15 by default). Although this is not the formal expression of a Lèvy process, this is a simple way to simulate a negative exponential rate in for the length of each step.







Here you have some examples using eWalk and their code to reproduce them:

   
 par(mfrow=c(2,2))

a=eWalk(0.15,0,0,1000, 0, 0, 0, 0, color="blue", levy=FALSE, plot=FALSE)
plot(a,type="n", main="Brownian motion")
lines(a,lwd=3, col="blue")
points(a, cex=.3, col="black", pch=21)

a=eWalk(0.15,0,0,1000, 0, 0, 0, 0, color="red", levy=TRUE, plot=FALSE)
plot(a,type="n", main="Lêvy flight")
lines(a,lwd=3, col="red")
points(a, cex=.3, col="black", pch=21)

a=eWalk(0.15,0,0,1000, 0.01, 20, 0, 0, color="orange", levy=FALSE, plot=FALSE)
plot(a,type="n", main="OU one optimum")
lines(a,lwd=3, col="orange")
abline(v=20, col="red", lwd="2")
points(a, cex=.3, col="black", pch=21)

a=eWalk(0.15,0,0,1000, 0.01, 20, 0.01, 20, color="green", levy=FALSE, plot=FALSE)
plot(a,type="n", main="OU two optimum")
lines(a,lwd=3, col="green")
points(a, cex=.3, col="black", pch=21)
abline(v=20, col="red", lwd="2")
abline(h=20, col="red", lwd="2")




 References


-Bartumeus, F., da Luz, M. E., Viswanathan, G. M., & Catalan, J. (2005). Animal search strategies: a quantitative random‐walk analysis. Ecology, 86(11), 3078-3087.

-Calenge, C. (2006). The package “adehabitat” for the R software: a tool for the analysis of space and habitat use by animals. Ecological modelling, 197(3), 516-519.

-Gould, S. J., & Eldredge, N. (1977). Punctuated equilibria: the tempo and mode of evolution reconsidered. Paleobiology, 3(02), 115-151.

-Landis, M. J., Schraiber, J. G., & Liang, M. (2013). Phylogenetic analysis using Lévy processes: finding jumps in the evolution of continuous traits. Systematic biology, 62(2), 193-204.

To leave a comment for the author, please follow the link and comment on their blog: long time ago....

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)