martes, 9 de mayo de 2017

niceOverPlot, or when the number of dimensions does matter

  Hi there!

    Over the last few months, my lab-mate Irene Villa (see more of her work here!)
and I have been discussing ecological niche overlap. The niche concept dates back to ideas first proposed by ornithologist J. Grinnell (1917). Later on, G.E. Hutchinson
(1957) defined the ecological niche of a species as the n-dimensional hyper-volume of ecological factors that define a space where the species can exist indefinitely. These ecological factors can be either abiotic (fundamental niche) or biotic variables, and the intersection of both defines the realized niche, or potential ecological range where the species can survive. The last necessary element for the existence of a species in a particular place (actual presence or distribution) is related to the “movement” ability or dispersal capacity of the species to arrive to those places where ecological factors are suitable. This theory is nicely summarized in the BAM diagram (Biotic-Abiotic-Movement; Soberón, 2007).

    In the last decades, with the development of cartographic information about climatic variables (mainly temperature and precipitation) and mathematical algorithms, many research efforts have been directed towards modeling the ecological niche of species from presence records and bioclimatic information (Elith & Leathwick, 2009). As these techniques were developed, many researchers focused on the comparison of ecological niches of two species, and the degree of overlap between them. A wide range of literature has been published about this topic (Peterson, 2011) but, in the last years, a general consensus has been reached about the methods to measure this overlap. Briefly, the method consists of two steps: 1) building an environmental space using bioclimatic values in the presence records and background points with any multivariate technique (usually Principal Component Analysis), and 2) representing both species niches in this space via kernel density, comparing them using D or I indexes and applying similarity or identity tests to assess significance (Warren et al., 2008).

    Following this general approach, the R package ecospat (Di Cola et al., 2016) provides an easy way to perform a Principal Component Analysis and to compare two species niches using D overlapping index, identity tests, etc. However, Irene and I realized that, when using this approach, in many cases for which niche overlap was not obtained, this was due to the lack of overlap for just one of the two PC axes, while a complete overlap was retrieved for the other axis. Taking into account that both axes do not usually contribute equally to explain environmental space (the first axis usually explains much more variability than the following ones), we worked on a function using the ggplot2 R package to create a nice plot for visual exploration of the Principal Component scores in both axes. This plot allows us to examine the factors that are or are not producing niche overlap, since it shows three plots in one: each axis separately and the joint distribution.

   Here is the function code (too large to print here). You can download and run it in an R console. The following is an example with two wall lizards from a Mediterranean area. You can download ecological values for the 19 bioclimatic variables at presence points of these two species here.

 #######   EXAMPLE   ######  
 # Read wall lizards ecological data  
 # Note: This data represent the values for the 19 bioclimatic variables at occurrence   
 # points for two Mediterranean wall lizards, Iberian wall lizard (Podarcis hispanicus,   
 # mainly in Iberian Peninsula) and Lilford's wall lizard (Podarcis lilfordi, from the   
 # Balearic Islands). Data were downloaded from GBIF* and Worldclim (  
 # The first 1125 rows are the bioclimatic values for P.hispanicus presences, while the  
 # next 44 rows represent bioclimatic values for P.lilifordi. The following rows represent   
 # available bioclimatic conditions (background) obtained from 10000 random points   
 # generated in a buffer of 400 km approximatively around presence points. Some inaccuracies   
 # (such as 2 Lilford's wall lizard presences in Iberian Peninsula) were not removed (this is   
 # just an example!!).  
 # *Data searches:  
 # Podarcis hispanicus: (30th April 2017) GBIF Occurrence Download  
 # Podarcis lilfordi: (30th April 2017) GBIF Occurrence Download  
 # From this data, we perform a Principal Component Analysis using "dudi.pca" function from "ecospat" package.  
 # We will choose 2 axes for the environmental space representation.  
 pca_Ph_Pl <-dudi.pca(na.omit(data_Ph_Pl, center = T, scale = T, scannf = F, nf = 2))  
 # Now, we can use directly this result with the niceOverPlot function to represent a   
 # 2D environmental space in a central plot, and each environmental gradient represented   
 # by each axis at top and right. We must provide to the function with the number of presences  
 # of each species in the same order as in the input data (nº of presences  
 # for Sp1 and nº of presences for Sp2)  
 niceOverPlot(pca_Ph_Pl, n1=1125 , n2= 44)  

Results from niceOverPlot (Blue area: P. lilifordi; Pink area: P. hispanicus)


  This plot could help us to interpret the results obtained from niche overlap analysis for these two species. In this example, we can see that, in 2D environmental space, the bioclimatic preferences of the two wall lizard species do not overlap. But you can notice that this lack of overlap is due to Axis2 of the PCA, which represents around 18% of contribution in environmental space. Now, we can focus on the bioclimatic variables that are involved in this axis to understand which factors are preventing the overlap.
In upcoming posts, we will discuss some alternative ways to address the D overlap index to take into account cases like our example.

Have fun!

Podarcis lilifordi, from the Balearic Islands


-Soberón, J. (2007). Grinnellian and Eltonian niches and geographic distributions of species. Ecology letters, 10(12), 1115-1123.

-Elith, J., & Leathwick, J. R. (2009). Species distribution models: ecological explanation and prediction across space and time. Annual review of ecology, evolution, and systematics, 40, 677-697.

-Peterson, A. T. (2011). Ecological niche conservatism: A time‐structured review of evidence. Journal of Biogeography, 38(5), 817-827.

-Di Cola, V., Broennimann, O., Petitpierre, B., Breiner, F. T., D'Amen, M., Randin, C., ... & Pellissier, L. (2016). ecospat: an R package to support spatial analyses and modeling of species niches and distributions. Ecography.

domingo, 4 de diciembre de 2016

rWind R package released!

Hi there! 

 Let me introduce you rWind, an R package with several tools for downloading, editing and converting wind data from Global Forecast System ( in other formats as raster for GIS! Wind data is a powerful source of information that could be used for many purposes in biology and other sciences: from the design of air pathways for airplanes to the study of the dispersion routes of plants or bird migrations. Making more accessible this kind of data to scientist and other users is the objective of ERDDAP (, a web service to dive into a lot of weather and oceanographic data-bases and download it easily. 

 I was using specifically one of the ERDDAP data-bases to get wind direction and speed from satellite data, the NOAA/NCEP Global Forecast System (GFS) Atmospheric Model ( At first, I was following this wonderful post from Conor Delaney ( to download and fix the data to be used as a GIS layer. However, I needed soon to download and modify a lot of wind data, so I started to write some R functions to automate the different tasks. Finally, I decided to put all together into an R package and upload it to CRAN repository to make it available for other users that could be interested in this kind of data. Here I give you a reference manual and an R code with a brief tutorial to get familiar with the utilities of the rWind package!

If you have any doubt or you want to report a bug or make any suggestion, please, comment the post or write me: or 

 Enjoy it!

Javier Fernández-López (2016). rWind: Download, Edit and Transform Wind Data from GFS. R package version 0.1.3.

 # Download and install "rWind" package from CRAN:  
 # You should install also "raster" package if you do not have it   
 # "rWind" is a package with several tools for downloading, editing and transforming wind data from Global Forecast   
 # System (GFS, see <>) of the USA's  
 # National Weather Service (NWS, see <>).  
 # > Javier Fernández-López (2016). rWind: Download, Edit and Transform Wind Data from GFS. R package version 0.1.3.  
 # >  
 # First, we can download a wind dataset of a specified date from GFS using wind.dl function  
 # help(wind.dl)  
 # Download wind for Spain region at 2015, February 12, 00:00  
 # help(wind.dl)  
 # By default, this function generates an R object with downloaded data. You can store it...  
 # or download a CVS file into your work directory with the data using type="csv" argument:  
 wind.dl(2015,2,12,0,-10,5,35,45, type="csv")  
 # If you inspect inside wind_data object, you can see that data are organized in a weird way, with  
 # to rows as headers, a column with date and time, longitude data expressed in 0/360 notation and wind  
 # data defined by the two vector components U and V. You can transform these data in a much more nice format
 # using "" function:  
 # Now, data are organized by latitude, with -180/180 and U and V vector components are transformed  
 # into direction and speed. You can export the data.frame as an CVS file to be used with a GIS software  
 write.csv(wind_data, "wind_data.csv")  
 # Once you have data organized by latitude and you have direction and speed information fields,  
 # you can use it to create a raster layer with wind2raster function to be used by GIS software or to be plotted   
 # in R, for example.  
 # As raster layer can only store one information field, you should choose between direction (type="dir")  
 # or speed (type="speed").  
 r_dir <- wind2raster(wind_data, type="dir")  
 r_speed <- wind2raster(wind_data, type="speed")   
 # Now, you can use rworldmap package to plot countries contours with your direction and speed data!  
 newmap <- getMap(resolution = "low")  
 plot(r_dir, main="direction")  
 lines(newmap, lwd=4)  
 plot(r_speed, main="speed")  
 lines(newmap, lwd=4)  
 # Additionally, you can use arrowDir and Arrowhead (from "shape" package) functions to plot wind direction  
 # over a raster graph:  
 alpha<- arrowDir(wind_data)  
 plot(r_speed, main="wind direction (arrows) and speed (colours)")  
 lines(newmap, lwd=4)  
 Arrowhead(wind_data$lon, wind_data$lat, angle=alpha, arr.length = 0.12, arr.type="curved")  
 # If you want a time series of wind data, you can download it by using a for-in loop:  
 # First, you should create an empty list where you will store all the data  
 wind_serie<- list()  
 # Then, you can use a wind.dl inside a for-in loop to download and store wind data of   
 # the first 5 days of February 2015 at 00:00 in Europe region. It could take a while...  
 for (d in 1:5){  
 # Finally, you can use wind.mean function to calculate wind average   
 r_average_dir<-wind2raster(wind_average, type="dir")  
 r_average_speed<-wind2raster(wind_average, type="speed")  
 plot(r_average_dir, main="direction average")  
 lines(newmap, lwd=1)  
 plot(r_average_speed, main="speed average")  
 lines(newmap, lwd=1)  

miércoles, 2 de noviembre de 2016

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

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.

 My colleague Javi and I were talking about 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 eMotion function in order implement this kind of movement models: eWalk (from evolutiveWalk...). eWalk is 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!) code to 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
  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
        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))
        route[generation,]<- loc
        generation= generation - 1  
        lon<- x      
        lat<- y
    else{                                                                   # if levy=FALSE, l=1, only Brownian motion or OU
      x= lon + (l * (sigma * rnorm(1))) + (alpha_x * ( theta_x - lon))
      y= lat + (l * (sigma * rnorm(1))) + (alpha_y * ( theta_y - lat))
      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 searching movements of predators (Bartumeus et al., 2005), or have been used to simulate punctuated equilibrium in traits evolution (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:


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")


-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.

domingo, 10 de julio de 2016

eMotion, an R function for evolutionary motion of character values

Hi all!

 Some time ago, I was studying an article from Tamara Münkemüller et al. 2015: "Phylogenetic niche conservatism – common pitfalls and ways forward". It´s a nice study where the ability of current tools in order to identify the evolutionary model that drives the changes in a trait (in this case, an ecological suitability) is discussed.
 They simulate the evolution of a trait under several evolutionary processes and then, try to apply the tools that are commonly used to test what kind of process is behind the evolution of the trait, allowing us to infer several patterns or processes as phylogenetic niche conservatism. Their results are really interesting, and I recommend you to take a look to the whole research!

However, I would like to focus in one of the main pictures of the article:

 Looks really nice! As a beginner, it was very difficult to me to understand this awesome figure. From my biological background, to study all this mathematical stuff requires a lot of time that I should invert in laboratory or field work... but it´s also very important to know the underlying statistical background of all the analysis that I would like to apply to my data in the future! So I decided to try to replicate some of those boxes by myself using R language. Many times, it´s the best way to learn!

 The figure shows trait evolution over time under Ornstein–Uhlenbeck 1 parameter in the left column, Ornstein–Uhlenbeck evolutionary at middle column, and Ornstein–Uhlenbeck extended (multiple optimum) at right column, with different selection strength in each row. This means that the first row is actually a Brownian Motion process, since alpha (selection strength) is 0.

 The graphic is the result of applying the OU process equation:

x(t+dt)= x(t)+ alpha(theta(t)-x(t))dt + (sigma dt E)

where , x(t) is the trait value at time=t, alpha is the selection strength with which a specific trait value (optimum) is wanted, theta is the optimal value of the trait in time=t, sigma is Brownian motion rate (the amount of random trait variability per time) and E is a normal distribution(0,1).

 Following this equation, I developed a simple and intuitive R function that simulate an OU evolutionary process for a hypothetical character in two species: eMotion (from "evolutionary motion")

 You can provide each parameter of the OU equation for the two species, and then you should define the number of generations (time) and replicates (number of simulations per species).

 Here you have the code to load the function into R console:

 eMotion <- function (sigma, theta, alpha, sigma2, theta2, alpha2, generations, n) {  
   plot(seq(1:generations),seq(1:generations), ylim=c(-200,200), type="n",   
       xlab="generations", ylab="trait value", main=paste("Sp1- s=", sigma,  
       " th=", theta, " a=", alpha,"  Sp2- s=", sigma2, " th=", theta2,   
       " a=", alpha2, sep=""))  
   simulations= list(1:n)  
   for (j in 1:n){  
     for (i in 1:generations) {   
       a= a + (sigma * rnorm(1)) + (alpha * ( theta - a))  
   for(k in 1:n){  
     lines(bb, col="red")  
   simulations= list(1:n)  
   for (j in 1:n){  
     for (i in 1:generations) {   
       a= a + (sigma2 * rnorm(1)) + (alpha2 * ( theta2 - a))  
   for(k in 1:n){  
     lines(bb, col="blue")  

 Once that you have loaded the function, you can use it to create simple plots to observe the effect of each parameter of OU equation in resultant graphic. For instance, here we describe an OU process of a trait in two species. The species 1 (red) follows an OU process of character evolution, with a sigma value (Brownian motion rate) of 2 and an optimum in 100, while species 2 (blue) follows an OU model with a sigma value of 0.5 and an optimum in -100. In both cases, alpha is 0.002.

eMotion(sigma, theta, alpha, sigma2, theta2, alpha2, generations, n)

 eMotion(2, 100, 0.002, 0.5,-100, 0.002, 10000, 20)  

The function automatically generates this simple plot, where you can observe that species 1 (red), is looking for a trait value around 100 (optimum), while species 2 (blue) is trying to get an optimum near to -100. Although both species have the same alpha (selection strength), the red one has a Brownian rate (sigma) much higher than the blue one. That explains why species 1 (red) has a higher range of variation than species 2 (blue), where all values are always near to the optimum (-100).

 You can use the function to explore other options in parameter values, which could be useful in order to teach about evolutionary processes or to understand the Ornstein–Uhlenbeck equation. Here you have one more example, enjoy!

 eMotion(1, 0, 0, 1,0, 0.01, 10000, 20)  

lunes, 27 de junio de 2016

nichEvolve documentation and code

Here you have some documentation about the nichEvolve model and a link to download the NetLogo code. I will write soon some practical experiments to test with the model... enjoy!
Link:  nichEvolve.nlogo

*(Don't you have NetLogo installed in your computer? NetLogo is a nice tool to implement simple agent-based models, and very useful to introduce yourself in the world of dynamic models!) Download from here: NetLogo


Here is a simple approach to explain how different models of trait evolution can work over one or two different species. Two available models of trait evolution are available in the general framework of Ornstein-Uhlenbeck process, with the possibility of adapting it into a simpler Brownian Motion process by changing the settings.
The model could be usefull as a virtual laboratory to teach about this concepts.

author: Javier Fernandez-Lopez


The model represent the evolution of a single ecological trait that could be defined as "ecological suitability", and could be understood for instance, as a specific requirement of temperature or humidity. The program allows an investigator to observe the effect of changes in each parameter of the Ornstein-Uhlenbeck equation

x(t+dt)= x(t)+ alpha(theta(t)-x(t))dt + (sigma dt E)

which define the evolution through time of a trait, where x(t) is the trait value at time=t, sigma is Brownian motion rate (the amount of random variability of the trait per time), theta is the optimal value of the trait in time=t, and alpha is the selection strengh with which a specific trait value is wanted.

 By changing equation parameters, the researcher can check the evolution of the trait value over time in the graphic box at the right side of the model. The "optimum" boxes show as well the theoretical value of the trait value. If speciation switch is "off", the model represents just one species. When speciation is turned "on", the species is split (red and blue) and the trait evolution occur independently. Now, theta parameter (OU-optimum) could differ for each species.

 Left side of the model shows a espatialy-splicit representation of the model. 200 mobile agents of each colour are shown over a grid that simules different ecological conditions. Remember that if speciation is "off", both colours represent the same species with identical ecological requirements. When the simulation starts, each agent try to find a patch around him that best match with his ecological requirements (theorical optimun niche values). The "realized" boxes show the mean of the values of the environmental conditions that the mobile agents have found by tracking the surronding environment. You can also compare the theorical optimum with the realized niche of mobile agents.

 Following KISS philosophy (Keep It Simple, Stupid!) we have not included some processes as extinction when theoretical and realized niches doesn't match, or reproduction of mobile agents.

 Finally, we have included as an experiment the option of "climate-change". When this switch is "on", it only affects to the spatially explicit representation of the model. The environment value of each patch increases by time, changing the niche availability for the mobile agents and affecting (or not) their distribution. The background process of trait evolution is not affected by "climate-change".


Click "Setup" to initialize the model. Then, you can use "Step" to see how the model runs step by step, or use "Go" button to run the model in a continuous way.

1) Start with a simple Brownian Motion model for a single species (speciation "off"). For this, alpha slide should be 0, theta's values are not important, since alpha (optimum selection strengh) is zero. You could test for changes in graphical box when sigma values are too large or too small.

2) Once that you know the effect of the Brownian motion rate (sigma), you could turn on speciation switch, to split the species. Now you can observe how both species are evolving independently, and, if you are lucky, you could see how their distributions differs from each other in the 2D display.

3) Now you can explore OU models increasing the value of alpha. You can try first using the same theta value for both species and changing it latter. You can also observe how interact Brownian motion rate (sigma) with optimum slection strengh (alpha)

4) Finally you could check the behaivour of mobile agents when climate-change is "on"


 By observing the patterns drawn in graphical box by different combinations of each parameter values, you could think about how difficult it is to guess the actual background process that drives a specific trait evolution if we only know the final state of the values.  


Fernandez-Lopez, J. 2016. nichEvolve: a niche evolution model