tag:blogger.com,1999:blog-90452406177308555742024-03-20T08:07:11.437-07:00long time ago...Ecological modeling, phylogeography and R tools to analyze biological processes and patternsJavier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.comBlogger20125tag:blogger.com,1999:blog-9045240617730855574.post-47429629738489366382022-11-26T14:06:00.000-08:002022-11-26T14:06:05.286-08:00Getting coordinates (lon,lat) from road name and kilometer in Spain<p> </p><p style="text-align: justify;"> Hi there!</p><p style="text-align: justify;"></p><p style="text-align: justify;"> Early this year we published a research about getting animal abundance estimates from wildlife - vehicles collisions in Ecography journal:</p><p style="text-align: justify;"><a href="https://onlinelibrary.wiley.com/doi/full/10.1111/ecog.06113">Can we model distribution of population abundance from wildlife–vehicles collision data?</a> <br /></p><p> </p><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEifDu421ZNT_bG2ZYkTX2ZCygyrEuw4pJsHjxfstRFxBv0Hh_FmxTs4Ls-N5f-jYIH8KnHncOgwGnPTr4RLlcxpay9fYKQNvchnRANKnihWxzE6ev15f-F-rmuDGEHn-OVU3qj3INsnwxjBRsRoCLI4KU0J_8fSRyS0xpY6CJphAJWSrnJ7HxG2uKurYg/s1248/ecog12885-fig-0002-m.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="1070" data-original-width="1248" height="461" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEifDu421ZNT_bG2ZYkTX2ZCygyrEuw4pJsHjxfstRFxBv0Hh_FmxTs4Ls-N5f-jYIH8KnHncOgwGnPTr4RLlcxpay9fYKQNvchnRANKnihWxzE6ev15f-F-rmuDGEHn-OVU3qj3INsnwxjBRsRoCLI4KU0J_8fSRyS0xpY6CJphAJWSrnJ7HxG2uKurYg/w538-h461/ecog12885-fig-0002-m.jpg" width="538" /></a></div><br /><p></p><p style="text-align: justify;">The first thing we had to deal with was to obtain geographical coordinates (latitude and longitude) from a database with the name of the roads and the kilometer point. At the end, I wrote a small chunk of code to do that by using the <a href="http://centrodedescargas.cnig.es/CentroDescargas/buscar.do?filtro.codFamilia=REDTR">Transport Network from the Spanish Geographic Institute</a>. I first downloaded all the kilometer point datasets (by provinces) and I merged them all in a single file. You can find this file <a href="https://github.com/jabiologo/roads/blob/main/puntosKilometricosSpain.csv">here</a> (but it's old and current datasets should be updated, so I recommend you to download it again from the original sources). Then I used R language to create a routine to "search" the coordinates that better matched my points.<br /></p><p style="text-align: justify;">I recently updated my code, warping it in a function and I created a toy-dataset with 5 points to show you how it works (see below, comments in Spanish).<br /></p><p style="text-align: justify;"> <br /></p><pre style="background: rgb(240, 240, 240); border: 1px dashed rgb(204, 204, 204); color: black; font-family: arial; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; overflow-wrap: normal; word-wrap: normal;"> library(dplyr)
library(RCurl)
# Guardamos las URLs de los datos almacenados en mi Github (https://github.com/jabiologo/)
roadData <- "https://github.com/jabiologo/roads/blob/main/roadData.RData?raw=true"
# Cargamos la base de datos del IGM con los puntos kilométricos
# Cargamos la base de datos con nuestros nombres de carreteras y kilómetros
load(url(roadData))
# Base de datos del IGN
head(km)
# Mis puntos para los cuales quier las coordenadas
# Nótese que si se quieren utilizar unos datos propios con este script,
# estos deben tener el mismo formato que el objeto "puntos"
head(puntos)</code></pre><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj4dQCu3yRp89WhdboVSy5NMBF2U6vivF6YFM2Tv9FP2Db29TPX32LQckio-LujvKRlPOWoqH_CPR7sQe14kKNLRKxVwcspkcatebJvgfl5GM2JKO1Hj8T8Mhb-4Mh1-EgZaWn3SMGgcmYvLAUUiAhSu1c_qKhs9L1Nu11beS5ul6l25FhW1LumAUyJ8w/s1207/Captura%20de%20pantalla%20de%202022-11-26%2022-35-37.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="203" data-original-width="1207" height="91" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj4dQCu3yRp89WhdboVSy5NMBF2U6vivF6YFM2Tv9FP2Db29TPX32LQckio-LujvKRlPOWoqH_CPR7sQe14kKNLRKxVwcspkcatebJvgfl5GM2JKO1Hj8T8Mhb-4Mh1-EgZaWn3SMGgcmYvLAUUiAhSu1c_qKhs9L1Nu11beS5ul6l25FhW1LumAUyJ8w/w536-h91/Captura%20de%20pantalla%20de%202022-11-26%2022-35-37.png" width="536" /></a></div><p></p><p></p><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg4rYhlCA2tJS9_51iZe2tbfkro7wsIBfNLupS_hGXeHRC6lDgkTzzpQPMpgYVT_KUidcaLCnvy4lwaDJtaXOVk1s2IHdltCN6cb0Mtqf2Ed0lGYFSzcRUtc3qwPHGU4c-HPX6tuXTtfeO9JlLJecQKKCR2w0ZT_5nbn_vUWmoTPv073ki4ZTJfx2ac2A/s182/Captura%20de%20pantalla%20de%202022-11-26%2022-35-43.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="176" data-original-width="182" height="94" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg4rYhlCA2tJS9_51iZe2tbfkro7wsIBfNLupS_hGXeHRC6lDgkTzzpQPMpgYVT_KUidcaLCnvy4lwaDJtaXOVk1s2IHdltCN6cb0Mtqf2Ed0lGYFSzcRUtc3qwPHGU4c-HPX6tuXTtfeO9JlLJecQKKCR2w0ZT_5nbn_vUWmoTPv073ki4ZTJfx2ac2A/w97-h94/Captura%20de%20pantalla%20de%202022-11-26%2022-35-43.png" width="97" /></a> <br /></div>
<pre style="background: rgb(240, 240, 240); border: 1px dashed rgb(204, 204, 204); color: black; font-family: arial; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; overflow-wrap: normal; word-wrap: normal;"> # Cargamos la función
pkm2xy <- function(events, database){
# Almacenaremos todas las carreteras que no se encuentren en nuestra base de datos
noMatch <- events[!(events[,1] %in% database$Nombre),]
# Filtramos nuestros datos
events <- events[events[,1] %in% database$Nombre,]
# Creamos una serie de columas para almacenar las coordenadas y otra info
events$X <- NA
events$Y <- NA
events$error <- NA
events$RoadAssigned <- NA
events$PkmAssigned <- NA
# Mediante un bucle iremos obteniendo nuestras coordenadas
for(i in 1:nrow(events)){
# Seleccionamos de nuestra base de datos aquel punto kilométrico que se
# encuentre más cerca de nuestro kilómetro
sel <- database %>% filter(Nombre == events[i,1]) %>%
filter(abs(numero - events[i,2]) == min(abs(numero - events[i,2])))
# Almacenamos las coordenadas
events[i,3:4] <- sel[1,1:2]
# Calculamos el error entre el punto kilométrico obtenido y el deseado
events[i,5] <- abs(sel$numero[1] - events[i,2])
# Almacenamos el nombre de la carretera y el punto kilométrico que hemos
# seleccionado en nuestra base de datos del IGN
events[i,6] <- sel$Nombre[1]
events[i,7] <- sel$numero[1]
print(i)
}
# Retornamos una lista con dos elementos:
# Aquellos puntos para los cuales se han conseguido las coordenadas
# Aquellos puntos para los que no se ha encontrado la carretera
return(list(events,noMatch))
}
# Corremos nuestra función con el fichero de ejemplo
pkm2xy(puntos,km)
</code></pre><p> </p><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgValgfXf7RY6tQkrhwGERhkltGahJ2OX3Oto2PsVk_HxxZ0g5ccSWx8LYvIvnrxprOhHd-q9M6U_pckY35muQPOWVPaBoB3gtNA-f32o7glRUhTBcsUluGFafYvkXaFyYe5ZF1nkwxXou4r5qbT74766VQ0ZHPdF4e0yeTlhN5qpMJUscppRiK01wSnw/s649/Captura%20de%20pantalla%20de%202022-11-26%2022-41-35.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="266" data-original-width="649" height="207" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgValgfXf7RY6tQkrhwGERhkltGahJ2OX3Oto2PsVk_HxxZ0g5ccSWx8LYvIvnrxprOhHd-q9M6U_pckY35muQPOWVPaBoB3gtNA-f32o7glRUhTBcsUluGFafYvkXaFyYe5ZF1nkwxXou4r5qbT74766VQ0ZHPdF4e0yeTlhN5qpMJUscppRiK01wSnw/w508-h207/Captura%20de%20pantalla%20de%202022-11-26%2022-41-35.png" width="508" /></a></div><br /><p></p><p style="text-align: justify;">As you can see, there was a road (CR-P-7111) that was not present in our database, so we couldn't find the coordinates for this point. However, we obtained the coordinates for the rest points. The function also show you the error in kilometers between the point and the latitude/longitude obtained. I recommend to perform a random test to visually check that everything was right, that is, if your dataset is large, select a number of points and manually check with Google Maps or similar if the function worked well.</p><p style="text-align: justify;">You can find the code and the data also in my <a href="https://github.com/jabiologo/roads/">GitHub repository</a>.<br /></p><p style="text-align: justify;"> </p><p style="text-align: justify;">Hope it's useful! </p><p style="text-align: justify;"> Best<br /></p><p style="text-align: justify;"> <u><b>References</b></u></p><div class="gs_citr" style="text-align: justify;" tabindex="0">Fernández‐López, J., Blanco‐Aguiar, J.
A., Vicente, J., & Acevedo, P. (2022). Can we model distribution of
population abundance from wildlife–vehicles collision data?. <i>Ecography</i>, <i>2022</i>(5), e06113.</div><p><br /></p>Javier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.com0tag:blogger.com,1999:blog-9045240617730855574.post-8464471701069963102021-06-28T22:12:00.002-07:002021-06-28T22:12:38.653-07:00Missing Migrants, tracking human deaths along migratory routes<p> Hi there,</p><p>
</p><p style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
Several years ago I received the book “<a href="https://www.msf.es/actualidad/mediterraneo/nunca-hubieran-sido">Como si nunca hubieran sido</a>”
as a Christmas present. In this poem, <a href="http://t.co/lAeixVvJ2X?amp=1">Javier Gallego</a> and <a href="https://juangallegopintura.jimdofree.com/">Juan Gallego</a>
talk about the real drama suffered by immigrants trying to reach
Europe through the Mediterranean route. The current pandemic
situation has <span lang="en-US">overshadowed</span> many problems
that affect millions of people each year. <a href="https://missingmigrants.iom.int/">Missing Migrants Project</a> tracks incidents involving migrants, including refugees and
asylum-seekers, who have died or gone missing in the process of
migration towards an international destination. It shows the
situation of many people who are forced to emigrate even during the
pandemic time. And it also helps me to contextualize my first world
problems. </p><p style="line-height: 100%; margin-bottom: 0cm; text-align: justify;"> </p><p style="line-height: 100%; margin-bottom: 0cm; text-align: justify;"></p><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgkJOBfzZPycx_83hhPWZMkMFDoc7sIpBOpiBUx-rNDfRH55zcmmwOJ_BI1xsexVhh3DO4caD8isLomsyiRiSFddE4VRX5vldpEenmDJ5Pntm74SirPavgJY4tVvAsmiFnre49l_IsiVsra/s900/hmi.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="597" data-original-width="900" height="265" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgkJOBfzZPycx_83hhPWZMkMFDoc7sIpBOpiBUx-rNDfRH55zcmmwOJ_BI1xsexVhh3DO4caD8isLomsyiRiSFddE4VRX5vldpEenmDJ5Pntm74SirPavgJY4tVvAsmiFnre49l_IsiVsra/w400-h265/hmi.png" width="400" /></a></div><br /><p></p><p style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">Here you have a simple analyses about the human migration
fatalities during the last years in the Mediterranean corridor using
R. </p><p style="line-height: 100%; margin-bottom: 0cm; text-align: justify;"> </p>
<pre style="background: rgb(240, 240, 240) none repeat scroll 0% 0%; border: 1px dashed rgb(204, 204, 204); color: black; font-family: arial; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; overflow-wrap: normal; word-wrap: normal;"> library(dplyr)
library(lubridate)
library(purrr)
library(incidence)
library(ggmap)
library(scales)
# Let's download and read the data from the Missing Migrants Project web site
url <- "https://missingmigrants.iom.int/global-figures/all/csv"
mm <- read.csv("/home/javifl/Descargas/MissingMigrants-Global-2021-06-27T21-15-07.csv")
# Filter the data for Mediterranean region
mm <- mm %>% filter(Region == "Mediterranean") %>% filter(Migration.Route != "")
mm$lon <- as.numeric(unlist(map(strsplit(mm$Location.Coordinates, ", "), 2)))
mm$lat <- as.numeric(unlist(map(strsplit(mm$Location.Coordinates, ", "), 1)))
mm$date <- mdy(mm$Reported.Date)
# plot the incidence of fatalities
mmInc <- incidence(mm$date, interval = 30, groups = mm$Migration.Route)
plot(mmInc, labels_week = F, stack = TRUE, border = "grey") +
scale_x_date(labels = date_format("%B %Y"))
</code></pre>
<p style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
</p>
<p style="line-height: 100%; margin-bottom: 0cm; text-align: justify;"> </p><p style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">We can plot the
monthly incidence of accidents involving migrants using the R package
<a href="https://cran.r-project.org/web/packages/incidence/vignettes/overview.html"><i>incidence</i></a>. Here you can see how
during the first months of 2020 the migration fatalities decreased,
probably due to the pandemic impact. However, human migration
incidents increased during the second half of 2020 despite of
pandemic restrictions in West and Central Mediterranean routes.</p><p style="line-height: 100%; margin-bottom: 0cm; text-align: justify;"> </p>
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjJYuWwBLBXNXNaKVWzTUBJ4n5nnPLur9YfgJmk8HRqHfLW-npcQfl7syEotjMOGMQ_pAQUBZ_O8I_Qec8eceNc7acMwrWdGTXFHcnuFUSR9_MagLMwsQ7E7X_4LjROPr-c7_EPvZXT6Wyt/s900/hist_missing.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="597" data-original-width="900" height="265" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjJYuWwBLBXNXNaKVWzTUBJ4n5nnPLur9YfgJmk8HRqHfLW-npcQfl7syEotjMOGMQ_pAQUBZ_O8I_Qec8eceNc7acMwrWdGTXFHcnuFUSR9_MagLMwsQ7E7X_4LjROPr-c7_EPvZXT6Wyt/w400-h265/hist_missing.png" width="400" /></a></div>
<p style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">Finally, we can plot monthly incidents using <i><a href="https://www.nceas.ucsb.edu/sites/default/files/2020-04/ggmapCheatsheet.pdf">ggmap</a></i> R package.</p>
<pre style="background: rgb(240, 240, 240) none repeat scroll 0% 0%; border: 1px dashed rgb(204, 204, 204); color: black; font-family: arial; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; overflow-wrap: normal; word-wrap: normal;"> months <- seq(ymd("2014-01-1"),ymd(max(mm$date)), by = '1 month')
myLocation<-c(min(mm$lon)-2.5, min(mm$lat)-2, max(mm$lon)+1, max(mm$lat)+1)
myMap<-get_map(location=myLocation, source="stamen", maptype="watercolor", crop=FALSE)
for (i in 1:length(months)-1){
int <- interval(months[i], months[i+1])
newmm <- mm[mm$date %within% int,]
id <- sprintf("%02d", i)
if (nrow(newmm) > 0){
png(paste("mm",id,".png", sep=""), width=900, height=610, units="px",
pointsize=18)
print(ggmap(myMap) +
geom_point(aes(x=lon, y=lat, size = 2), data=newmm, col="orange", alpha=0.5) +
#geom_density2d(aes(x=lon, y=lat), data=newmm, size = 0.3) +
#stat_density2d(data = newmm,
# aes(x = lon, y = lat, fill = ..level.., alpha = ..level..), size = 0.01,
# bins = 16, geom = "polygon") + scale_fill_gradient(low = "green", high = "red") +
scale_alpha(range = c(0, 0.2), guide = FALSE) +
theme(legend.position = "none", plot.title = element_text(size = 25)) +
ggtitle(paste(year(months[i]), month(months[i], label = T, abbr = F), sep=" - ")))
dev.off()
}
print(id)
}
</code></pre>
<p><style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent }</style></p><p></p><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEia_VY2FhjRrH6omHpt1hiyaLH-vqQKw25JgDmYYnv81eoR1gmIdlEbLS1P0xNwhP0Pw-5NkUVF0ux7xAvQ2kbDwDelY6OjC2X-Z-WM529dTg6O0g9ijzFMetjcQ5lI1t3CT2E3JeyXKghN/s900/mm.gif" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="597" data-original-width="900" height="265" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEia_VY2FhjRrH6omHpt1hiyaLH-vqQKw25JgDmYYnv81eoR1gmIdlEbLS1P0xNwhP0Pw-5NkUVF0ux7xAvQ2kbDwDelY6OjC2X-Z-WM529dTg6O0g9ijzFMetjcQ5lI1t3CT2E3JeyXKghN/w400-h265/mm.gif" width="400" /></a></div><br /><p><a href="https://www.youtube.com/watch?v=4E1vIbuWJ0U">Stay safe!</a></p><p>"No soy un número ni parte de una cifra </p><p>aunque se paga por igual la misma tarifa..."<br /></p><p><br /></p>Javier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.com6tag:blogger.com,1999:blog-9045240617730855574.post-8396801939679938212021-01-22T22:07:00.005-08:002021-01-24T06:14:44.760-08:00Are betting houses randomly distributed in Madrid? Point pattern analysis in R<p align="justify"><span style="font-size: small;">
<span style="font-size: medium;">Hi there! <br /></span></span><span style="font-size: medium;">
<span><br /> Last December 13th there was a demonstration in Carabanchel, my home neighborhood in Madrid (Spain) against the proliferation of betting houses (small casinos where you can get cheap drinks or coffees while you spend your money betting or gambling). Many neighborhood associations usually complain against this kind of places arguing they bait young people with very low prices, who can develop gambling addiction in the future. Moreover some organizations ensure that the proliferation of this kind of places is focused in working class neighborhoods with low incomes and high unemployed rates (see this <a href="https://stopcasasdeapuestas.com/static/docs/informe-locales-apuestas-fravm.pdf" target="_blank">great report</a> for more information about that - in Spanish). So, how could we know if betting houses are randomly distributed in Madrid?
<br /><br /></span> <span><t>Imagine that you have a 25x25 (625 cells) square grid and you want to scatter 2500 points. If the points are randomly distributed, <b>they will have the same probability of occupy whatever cell</b>. That is, the first point will have a probability of 1/625 of being in cell number 1, the same probability (1/625) probability of being at cell number 2, and so on. Thus, if we have 2500 points in 625 cells, each cell would have 2500/625 = 4 points in average. This is the meaning of the <b>lambda</b> from a Poisson distribution, <b>the mean or the most probable value</b>.
</t></span></span></p><div class="separator" style="clear: both;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjqdGFUUh1uPhrgdM3Xo1cho2apdaN-Gkn92jrDm3EJ0_Begbr8PN9ZSboRDktbKoc7E_YNtR2l8b2-nw-UPh56aDEWcWde81YmlMLBN1v0J-f8TuNAUvn8FqRsMpgaBC7uyeaTZaCkc2_S/s1532/fig1.png" style="display: block; padding: 1em 0px; text-align: center;"><img alt="" border="0" data-original-height="570" data-original-width="1532" height="209" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjqdGFUUh1uPhrgdM3Xo1cho2apdaN-Gkn92jrDm3EJ0_Begbr8PN9ZSboRDktbKoc7E_YNtR2l8b2-nw-UPh56aDEWcWde81YmlMLBN1v0J-f8TuNAUvn8FqRsMpgaBC7uyeaTZaCkc2_S/w561-h209/fig1.png" width="561" /></a>
<p align="justify"><span style="font-size: small;"> <span style="font-size: medium;">Imagine now that we have the same points in the same grid, but now they are clustered in 3 groups. In this dummy example, we can easily guess that they are not randomly distributed by just observing the point pattern. But how could we demonstrate mathematically/numerically they are not randomly distributed?
</span>
</span></p><div class="separator" style="clear: both;"><span style="font-size: small;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjqjJbTWQKix76sWwTR9WO94-ME0jpoiu0V_7ufxQy1UHqCD9ffK-liOVgpVjBiXAdRx7E0sTmFYaKY4Xwt7ueDFOP0sG9qEep_N8JUFVbzWhmaX5qRXXPViat_R8pcs3kjlWBHfvYMaqWX/s1532/fig2.png" style="display: block; padding: 1em 0px; text-align: center;"><img alt="" border="0" data-original-height="570" data-original-width="1532" height="205" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjqjJbTWQKix76sWwTR9WO94-ME0jpoiu0V_7ufxQy1UHqCD9ffK-liOVgpVjBiXAdRx7E0sTmFYaKY4Xwt7ueDFOP0sG9qEep_N8JUFVbzWhmaX5qRXXPViat_R8pcs3kjlWBHfvYMaqWX/w552-h205/fig2.png" width="552" /></a></span></div><span style="font-size: small;"><span style="font-size: medium;">
</span><p align="justify"><span style="font-size: medium;"> Well, given we know the shape of the histogram (distribution) that should have a completely random distribution (in our case a <b>Poisson(λ=4)</b>), an easy and straightforward way to demonstrate mathematically that the second point patter is not random at all would be to compare both distributions using for example a Kolmogorov–Smirnov test in R.<br /><br /></span><span style="font-size: medium;"> Lets do the same with betting houses in Madrid city (Spain). Betting places locations can be downloaded from the city government website, as well as other data such as neighborhood limits or family income per neighborhood.I prepared a .zip with the needed files which can be downloaded from <a href="https://drive.google.com/file/d/1ZyrgVzj1yuraFE4-15CAANJVMrfQE231/view?usp=sharing">here</a>.</span></p><p align="justify"><span style="font-size: medium;">There are currently 409 betting houses. We will divide the city in a grid of 18x18 cells of 1 km size as follows</span></p><p style="line-height: 100%; margin-bottom: 0cm;"></p><p style="line-height: 100%; margin-bottom: 0cm;"><span lang="en-US"> </span></p>
<pre style="background: rgb(240, 240, 240) none repeat scroll 0% 0%; border: 1px dashed rgb(204, 204, 204); color: black; font-family: arial; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; overflow-wrap: normal; word-wrap: normal;"> # Load required packages
library(rgdal)
library(raster)
library(dismo)
library(GISTools)
# Load data
nb <- readOGR("/home/javifl/gambling/SHP_ETRS89/barrios_income_pop.shp")
bh <- read.csv("/home/javifl/gambling/betting_house.csv")
# Create the grid
bh.sp <- SpatialPointsDataFrame(coords=bh[,1:2], data=bh)
ext <- extent(c(min(bh$x),min(bh$x)+18000,min(bh$y),min(bh$y)+18000))
landscape <- raster(ncols=18,nrows=18,ext=ext, vals=1:324)
grid <- rasterToPolygons(landscape)
# Plot betting houses in Madrid
par(mar=c(0.1,0.1,1,0.1))
plot(nb, main = "City of Madrid")
plot(grid, lwd = 0.4, add = TRUE)
points(bh.sp, pch=16, cex=0.7, col="black")
</code></pre>
</span></div><p> </p><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhvUolQqQcLEJmsmZvka73LXS3KfpPtqVZCFOAO-6L3oiqFvnb2hpHpizyMkjKPHv859mllA6WaeJYKL8Hk0cFLZEzCtW3DyOtZqeCJMAV02DhpaVqZlHq43WeeuSKWG6z3v4jj1_KLfRLL/s667/fig3.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="667" data-original-width="615" height="548" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhvUolQqQcLEJmsmZvka73LXS3KfpPtqVZCFOAO-6L3oiqFvnb2hpHpizyMkjKPHv859mllA6WaeJYKL8Hk0cFLZEzCtW3DyOtZqeCJMAV02DhpaVqZlHq43WeeuSKWG6z3v4jj1_KLfRLL/w506-h548/fig3.png" width="506" /></a></div><p style="text-align: left;"><br /></p><p style="text-align: justify;"><span style="font-size: medium;"> Cool. Since Madrid is an heterogeneous city, we will select for our analysis only those cells in which there is one or more betting houses. Then, we'll spread random points in the selected squares to compare both histograms/distributions: random points VS betting houses. </span><br /></p><p>
</p><pre style="background: rgb(240, 240, 240) none repeat scroll 0% 0%; border: 1px dashed rgb(204, 204, 204); color: black; font-family: arial; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; overflow-wrap: normal; word-wrap: normal;"> # Select cells with betting houses
grid$bh <- poly.counts(bh.sp, grid)
sel.sqr <- (grid[grid$bh > 0,])
landscape[grid$bh == 0] <- NA
# Plot
par(mar=c(0.1,0.1,0.1,0.1))
plot(sel.sqr)
lines(nb, lwd = 0.2)
points(bh.sp, pch=16, cex=0.7, col="black")
# Add random points in selected cells
rp <- randomPoints(disaggregate(landscape, 100), nrow(bh))
points(rp, pch=16, cex=0.7, col = "red" )
rp.sp <- SpatialPointsDataFrame(coords=rp[,1:2], data=data.frame(rp))
</code></pre><p>
</p><br /><p></p><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhJEy_QQHQAQaAwH2QSjsB4VZ7OayTO7g9nWgUNaDsSqLRVm1KcIA7qF2yRSk036PGoQyolJ71I8ZHHA3ebMHgyvT0wovEZ6oqVRPVF1iJAZ15b3o7aedi_iix-F7rbO-q-yHzuDewKTTwR/s667/fig4.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="667" data-original-width="615" height="498" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhJEy_QQHQAQaAwH2QSjsB4VZ7OayTO7g9nWgUNaDsSqLRVm1KcIA7qF2yRSk036PGoQyolJ71I8ZHHA3ebMHgyvT0wovEZ6oqVRPVF1iJAZ15b3o7aedi_iix-F7rbO-q-yHzuDewKTTwR/w459-h498/fig4.png" width="459" /></a></div><br /><p></p><p></p><p style="text-align: justify;"> <span style="font-size: medium;">Well done! We have selected 118 cells with at least one betting house and we have scattered 409 random points on them. So, now the last task is to count the number of betting houses and the number of random points at each cell. Then, we could compare both histograms/ distributions, and decide if they present a similar pattern.</span><br /></p><p>
</p><pre style="background: rgb(240, 240, 240) none repeat scroll 0% 0%; border: 1px dashed rgb(204, 204, 204); color: black; font-family: arial; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; overflow-wrap: normal; word-wrap: normal;"> par(mfrow = c(1,2))
hist(poly.counts(rp.sp, sel.sqr), right=F, main = "random points",
xlab="number of random points by cell")
hist(poly.counts(bh.sp, sel.sqr), right=F, main = "betting houses",
xlab="number of betting houses by cell")
</code></pre>
<br /><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgTm5XTjANK9IN3WH3J0Zqp-nqYmcm9XPKw33RMWoeGvB5-dgzmPFPVoB1ok0kH8zoo0iFr_uKwYm1r0JOKm8ACx7m35W6E1WGWlqPxRytk2uDW16QX1Z0jPGTpC4AxJT5BSfr0_JVUnPTx/s1255/fig5.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="595" data-original-width="1255" height="260" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgTm5XTjANK9IN3WH3J0Zqp-nqYmcm9XPKw33RMWoeGvB5-dgzmPFPVoB1ok0kH8zoo0iFr_uKwYm1r0JOKm8ACx7m35W6E1WGWlqPxRytk2uDW16QX1Z0jPGTpC4AxJT5BSfr0_JVUnPTx/w547-h260/fig5.png" width="547" /></a></div><div style="text-align: justify;"><br /></div><p style="text-align: justify;"><span style="font-size: medium;">As you can see, under a random point pattern the most frequent value is between 3 and 4. Since we have 409 points randomly distributed over 108 cells, we expected to have 409/108 = 3.78 points at each cell. It looks pretty good!</span></p><p style="text-align: justify;"><span style="font-size: medium;">However, we can see that the most frequent value in the betting house histogram is between 0 and 2... This is because there are many empty cells, while there are some of them with high number of betting houses (a clustered pattern). We could compare both distributions mathematically using a </span><span style="font-size: small;"><span style="font-size: medium;">Kolmogorov–Smirnov test.<br /></span>
</span></p><pre style="background: rgb(240, 240, 240) none repeat scroll 0% 0%; border: 1px dashed rgb(204, 204, 204); color: black; font-family: arial; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><span style="font-size: small;"><code style="color: black; overflow-wrap: normal; word-wrap: normal;"> #Kolmogorov–Smirnov test
ks.test(poly.counts(rp.sp, sel.sqr), poly.counts(bh.sp, sel.sqr))
</code></span></pre><span style="font-size: small;">
</span><p></p><p><span style="font-size: medium;"> </span></p><p style="text-align: justify;"><span style="font-size: medium;">Bonus Track: so, what variable or variables drive the distribution of betting houses in Madrid? Here you have a clue... but take the data and explore it yourself!</span></p><p><span style="font-size: medium;"> <br /></span></p><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEghDabHDk4ZBmGV72aKm1PYVKXuKJ-LZ6AHF9ZApmQLJpCAdPMMWTKwY7zXFgbBdssnCz4Z4pwajGaJ7MEmD0MKZ_tlCPCQPASNBwOUE2_X2NiCNb-4uxouA_9KBfAu7heivpO6IhHOGtY9/s790/Rplot.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="702" data-original-width="790" height="388" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEghDabHDk4ZBmGV72aKm1PYVKXuKJ-LZ6AHF9ZApmQLJpCAdPMMWTKwY7zXFgbBdssnCz4Z4pwajGaJ7MEmD0MKZ_tlCPCQPASNBwOUE2_X2NiCNb-4uxouA_9KBfAu7heivpO6IhHOGtY9/w438-h388/Rplot.png" width="438" /></a></div><p></p><p><span style="font-size: medium;"> </span></p><p><span style="font-size: medium;">Stay safe!</span></p><p><span style="font-size: medium;"> </span></p><p><span style="font-size: medium;"> </span></p><div class="separator" style="clear: both; text-align: center;"><span style="font-size: medium;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgdaADDK7l-sv11Sh_JBz1nrI5u7PW0TVwAD9lUqco1erwdbC4g5gCoJctlWizfi2k3S9wopRkuARuE1w326LLn6iDV7wYPyLasMvAvh6dARwj1tU2oFWyhTXlpldVzb6iDChqy-e6S9Shn/s2048/IMG_20201213_131006_237.jpg" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="1536" data-original-width="2048" height="373" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgdaADDK7l-sv11Sh_JBz1nrI5u7PW0TVwAD9lUqco1erwdbC4g5gCoJctlWizfi2k3S9wopRkuARuE1w326LLn6iDV7wYPyLasMvAvh6dARwj1tU2oFWyhTXlpldVzb6iDChqy-e6S9Shn/w497-h373/IMG_20201213_131006_237.jpg" width="497" /></a></span></div><span style="font-size: medium;"><br /></span><p></p><p><span style="font-size: medium;"> </span></p><p><span style="font-size: medium;"> </span></p><p></p><p></p>Javier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.com7Carabanchel, 28047 Madrid, España40.38783 -3.7448740.384561318360596 -3.7491615344238283 40.391098681639406 -3.740578465576172tag:blogger.com,1999:blog-9045240617730855574.post-12880100254938860342020-07-09T23:38:00.000-07:002020-07-09T23:38:06.262-07:00N-mixture models to estimate animal abundance with R and unmarked<div style="text-align: justify;">
<br /></div>
<div style="text-align: justify;">
Hi there! </div>
<div style="text-align: justify;">
<br /></div>
<div style="text-align: justify;">
During the pandemic lockdown I was studying how N-mixture models to estimate animal abundance works. As a result, <a href="https://jabiologo.github.io/web/tutorials/nmixture.html" target="_blank">I wrote this small tutorial in R</a> (Spanish) that helped me to better understand such models. Any correction or suggestion can be made by posting here or at jflopez.bio@gmail.com.<br /><br />Hope you enjoy it!</div>
<br />
<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjTzxIBNtqVCcV123E-Kc8Dtq62sCn_PHBrXNZ_FdR08_us13tdUHZeEUJgDDwROkrW3pNrVzVuzl3k-5mVnfkbnYJSfxz380FdD1V_Q6ONPiikWiclR8vP60eS_qmSsq4eZJoaPr7lNW8C/s1600/post.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="828" data-original-width="827" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjTzxIBNtqVCcV123E-Kc8Dtq62sCn_PHBrXNZ_FdR08_us13tdUHZeEUJgDDwROkrW3pNrVzVuzl3k-5mVnfkbnYJSfxz380FdD1V_Q6ONPiikWiclR8vP60eS_qmSsq4eZJoaPr7lNW8C/s400/post.png" width="399" /></a></div>
<br />Javier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.com2tag:blogger.com,1999:blog-9045240617730855574.post-12602646664257500202020-04-03T17:16:00.000-07:002020-04-03T17:16:07.386-07:00Another "flatten the COVID-19 curve" simulation... in R<br />
Hi there,<br />
<br />
This is the best meme I've found during these days...<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjqIMzJz3EVfxv9tj7x7lYmAcGmx-PKdSbVwB0j5QtqDaDZFD_YkgpXUdF63lBdO1CWgts6IyD3-BxHK05go0vMKwLDJ78T3ZLFRhCtwCUtG1J2-9WzaFGskROTTTqHSNaBKcqy1GJwzmHk/s1600/89761491_2926059660770912_698880553033662464_n.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="500" data-original-width="750" height="266" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjqIMzJz3EVfxv9tj7x7lYmAcGmx-PKdSbVwB0j5QtqDaDZFD_YkgpXUdF63lBdO1CWgts6IyD3-BxHK05go0vMKwLDJ78T3ZLFRhCtwCUtG1J2-9WzaFGskROTTTqHSNaBKcqy1GJwzmHk/s400/89761491_2926059660770912_698880553033662464_n.jpg" width="400" /></a></div>
<br />
<br />
<br />
Well, here it is my "BUT" contribution. Some weeks ago The Washington Post published this <a href="https://www.washingtonpost.com/graphics/2020/world/corona-simulator/" target="_blank">simulations</a> about how "social distancing" could help to "flat the curve" of COVID-19 infections. I fell in love with these simulations because their simplicity and explanatory power. Indeed, you can use the pretty similar principles to simulate predator <a href="http://allthiswasfield.blogspot.com/2018/07/antipredator-behavior-with-r-or-why.html" target="_blank">hunt behavior</a> and other movement patterns... I wrote some R code to have a tiny version of them...<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgLGsixFbGK2o0OgDSbkLWfF9aW70R4dVCOA6zqzYNdHOrI0OzoGjpZdcumg8zwxvdE8WvQfF6SzVIRsVCTYvMSL7vahkg4Qz3BOUr4wd1Vo9OZBsEhOU72sqn2cFddHIhNC3v6T3GBiZEP/s1600/corona_469.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="800" data-original-width="780" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgLGsixFbGK2o0OgDSbkLWfF9aW70R4dVCOA6zqzYNdHOrI0OzoGjpZdcumg8zwxvdE8WvQfF6SzVIRsVCTYvMSL7vahkg4Qz3BOUr4wd1Vo9OZBsEhOU72sqn2cFddHIhNC3v6T3GBiZEP/s400/corona_469.png" width="390" /></a></div>
<br />
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> library(raster)
library(dismo)
r <- raster(nrows = 100, ncols = 100, xmn = 0, xmx = 100, ymn = 0, ymx = 100)
r[] <- 0
steps <- 500
n <- 100
locations <- data.frame(matrix(ncol = n, nrow = steps))
pp <- randomPoints(r,n)
cell <- cellFromXY(r, pp)
infected <- 1:10
pob <- 1:n
ratio <- data.frame(1:steps, NA)
pob_inf <- list()
for (j in 1:steps){
print(j)
locations[j,] <- cell
for(i in 1:n){
cell[i] <- adjacent(r, cell[i], 8)[round(runif(1,0.51,nrow(adjacent(r, cell[i], 8))+0.49),0),2]
}
new_inf <- cell %in% cell[infected]
infected <- pob[new_inf]
ratio[j,2] <- length(infected)
pob_inf[[j]] <- infected
}
locations2 <- data.frame(matrix(ncol = n, nrow = steps))
cell2 <- cellFromXY(r, pp)
infected2 <- 1:10
pob2 <- 1:n
ratio2 <- data.frame(1:steps, NA)
pob_inf2 <- list()
for (j in 1:steps){
print(j)
locations2[j,] <- cell2
for(i in 1:25){
cell2[i] <- adjacent(r, cell2[i], 8)[round(runif(1,0.51,nrow(adjacent(r, cell2[i], 8))+0.49),0),2]
}
new_inf2 <- cell2 %in% cell2[infected2]
infected2 <- pob2[new_inf2]
ratio2[j,2] <- length(infected2)
pob_inf2[[j]] <- infected2
}
</code></pre>
<br />
Let's make some plots to <a href="http://allthiswasfield.blogspot.com/2017/12/p-margin-bottom-0.html" target="_blank">put them together in a GIF</a> and better visualize the results...<br />
<br />
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> num <- seq(1,500,4)
for (p in 1:125){
id <- sprintf("%03d", num[p])
png(paste("corona_",id,".png", sep=""), width=780, height=800, units="px", pointsize = 15)
layout(matrix(c(1,1,2,3),2,2,byrow = TRUE))
plot(ratio[1:num[p],],pch=19, xlim = c(0,500), ylim = c(0,100), ylab = "nº of infected",
xlab = "time", col = "red", cex.axis = 1.4, cex.lab = 1.4, main = "Infected curves", cex.main = 2)
points(ratio2[1:num[p],],pch=19,col="blue")
legend("topright", legend = c("free movement", "restricted movement"),lwd = c(4,4), col = c("red","blue"), cex = 1.5 )
plot(r, col = "white", legend = FALSE, axes = FALSE, main = "free movement", cex.main = 1.8)
points(xyFromCell(r,as.numeric(locations[num[p],])), cex = 1.2, col = "grey40", pch = 18)
points(xyFromCell(r,as.numeric(locations[num[p],]))[pob_inf[[num[p]]],], pch = 18, col = "red", cex = 1.4)
plot(r, col = "white", legend = FALSE, axes = FALSE, main = "restricted movement", cex.main = 1.8)
points(xyFromCell(r,as.numeric(locations2[num[p],])), cex = 1.2, col = "grey40", pch = 18)
points(xyFromCell(r,as.numeric(locations2[num[p],]))[pob_inf2[[num[p]]],], pch = 18, col = "blue", cex = 1.4)
dev.off()
}
</code></pre>
<br />
<br />
Done! <br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgmaIyoyFq-bsRqhbSuL9XbRwCddhK1uDhEnrUvXITdpNt89k1bEFtJX3ih4dYG9yfaMqSptfYr7dWWeNs5DS5L2FiQfb1GNNyAdKAAi0SipkAk34kDddCQiXklvdK5e1s3q6jC0ZWtYTzL/s1600/corona.gif" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="800" data-original-width="780" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgmaIyoyFq-bsRqhbSuL9XbRwCddhK1uDhEnrUvXITdpNt89k1bEFtJX3ih4dYG9yfaMqSptfYr7dWWeNs5DS5L2FiQfb1GNNyAdKAAi0SipkAk34kDddCQiXklvdK5e1s3q6jC0ZWtYTzL/s400/corona.gif" width="390" /></a></div>
<br />
Stay safe!<br /><br />
<br />
<br />
<br />
<br />
<br />
<br />Javier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.com2tag:blogger.com,1999:blog-9045240617730855574.post-86159576169076836322020-02-18T14:53:00.001-08:002020-02-18T14:53:42.570-08:00rWind is working again!<br />
Yep, after several months fallen, <a href="https://cran.r-project.org/web/packages/rWind/index.html" target="_blank">rWind R package</a> is working again! I'm sorry, but I'm too busy lately with my PhD dissertation and I have not all the time that I'd like to improve rWind 😞. The problem was due to the change of the URL to GFS server, and it was solved thanks to a <a href="https://github.com/jabiologo/rWind/pull/14" target="_blank">pull request in GitHub</a>... Thank you very much! The new version on CRAN (1.1.5) comes with some new functions to download sea current data from the <a href="https://coastwatch.pfeg.noaa.gov/erddap/info/jplOscar_LonPM180/index.html" target="_blank">OSCAR</a> database! New examples with sea and wind currents coming soon... also <span class="css-901oao css-16my406 r-1qd0xha r-ad9z0x r-bcqeeo r-qvutc0">on Twitter @javi_ferlop</span><br />
<br />
Enjoy!<br />
<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhGGb6gho8XjkawkEitH7VEwU8cwgshn7QvW780v6BMGvlFJ3oSX2vzTOLVAfrmc35gtsXLm0tDaVVWwdaLtDemaDXMwEeaze38pSVfvZedlh03iI6kw1XIUsXGIgb7LL1ETqjVjfEY0KLW/s1600/currents.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="489" data-original-width="776" height="251" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhGGb6gho8XjkawkEitH7VEwU8cwgshn7QvW780v6BMGvlFJ3oSX2vzTOLVAfrmc35gtsXLm0tDaVVWwdaLtDemaDXMwEeaze38pSVfvZedlh03iI6kw1XIUsXGIgb7LL1ETqjVjfEY0KLW/s400/currents.png" width="400" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
</div>
<br />Javier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.com0tag:blogger.com,1999:blog-9045240617730855574.post-69730461801374404132019-04-14T14:36:00.002-07:002019-04-14T14:36:49.435-07:00Understanding Bayesian Inference with a simple example in R!<br />
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
Hi there!</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
Last summer, the
Royal Botanical Garden (Madrid, Spain) hosted the first edition of
MadPhylo, a workshop about Bayesian Inference in phylogeny using
RevBayes. It was a pleasure for me to be part of the organization staff with
John Huelsenbeck, Brian Moore, Sebastian Hoena, Mike May, Isabel
Sanmartin and Tamara Villaverde. Next edition of Madphylo <a href="https://www.madphylo.com/" target="_blank">will be held June 10, 2019 to June 19, 2019at the Real Jardín Botánico de Madrid</a>. If you are interested in
Bayesian Inference and phylogeny just can't miss it! You'll learn the RevBayes language, a programming language to perform phylogeny (and other) analyses under a Bayesian framework! Apply <a href="https://www.madphylo.com/apply" target="_blank">here</a>!</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhga6kA7jBiXWh2xtZLESUISQdNfZsMwouvGvr2U4KydX3ZXeSG3EzC8nE2icYOCQ1tkNFXdBheJil-vcaYDTj4rReDFkDOIpOA2e6LMuOfOnpHDUNr0MF6rEC1Hp6dmm15P6HOmu20IM-i/s1600/pic_1.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="450" data-original-width="1023" height="175" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhga6kA7jBiXWh2xtZLESUISQdNfZsMwouvGvr2U4KydX3ZXeSG3EzC8nE2icYOCQ1tkNFXdBheJil-vcaYDTj4rReDFkDOIpOA2e6LMuOfOnpHDUNr0MF6rEC1Hp6dmm15P6HOmu20IM-i/s400/pic_1.png" width="400" /></a></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
The first days were focused to explain how we can use the Bayesian
framework to estimate the parameters of a model. I'm not an expert in Bayesian Inference at all, but in this post I'll try to reproduce one of the first Madphylo <a href="https://revbayes.github.io/tutorials/mcmc_binomial/" target="_blank">tutorials</a> in R language. As a simple example,
we'll use a coin flipping experiment. We can model this experiment with
a binomial distribution:</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
binomial(<i>n</i>, <i>p</i>)</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
where <i>n</i> is
the number of tosses and the parameter <i>p</i> is the probability to
obtain a head. Given the number of tosses and the number of heads (<i>h</i>),
we can use Bayesian Inference to compute the probability to obtain a
head (<i>p</i>).</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
If we use Bayes’ theorem, we have that the probability of a
specific value of <i>p</i> given the number of heads (<i>h</i>) is:</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style> </div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style> </div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEghjSeQbhCjrYfNCvrFq1j3QPuItFGuSqXJjpgleYi2zScRsb0E1D1HxMQKrYmOec4468Q2vXu1vnU7OdgsiFLE4upstWCC2cyz16bGktiIFK_VyzjvT6B_yX2U3LXzGJsmWm-p4dC9VFTA/s1600/bayes1.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="216" data-original-width="721" height="59" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEghjSeQbhCjrYfNCvrFq1j3QPuItFGuSqXJjpgleYi2zScRsb0E1D1HxMQKrYmOec4468Q2vXu1vnU7OdgsiFLE4upstWCC2cyz16bGktiIFK_VyzjvT6B_yX2U3LXzGJsmWm-p4dC9VFTA/s200/bayes1.png" width="200" /></a></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style>
</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
where Pr(<i>h</i>|<i>p</i>) is the likelihood of the data given a
value of <i>p</i>, Pr(<i>p</i>) is the prior probability for <i>p</i>
and Pr(<i>h</i>) is the marginal likelihood for the data. The prior
probability of our parameter represent our believe about what values
can take the parameter. If we think our parameter can take whatever
value with the same probability, we use an uniform (flat) prior.</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
We can use a Markov Chain Monte Carlo (MCMC) to introduce many
different values of <i>p</i> and get the posterior probability of
these values. To do this, we can use the Metropolis-Hastings MCMC
algorithm, with the next steps (simplified):</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
- Step 1) Set an initial value for <i>p.</i></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> p <- runif(1, 0, 1)
</code></pre>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
- Step 2) Propose a new value of <i>p</i>, called <i>p</i>'.</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style> </div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> p_prime <- p + runif(1, -0.05, 0.05)
</code></pre>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
- Step 3) Compute the acceptance probability of this new value for the
parameter. We have to check if the new value improves the
posterior probability given our data. This can be seen as the ratio:
Pr(<i>p'</i>|<i>h</i>) / Pr(<i>p</i>|<i>h</i>).</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style> </div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhlje1SvRaMn40jCaaKRu7E52JPl9zKsgSQ0pRBXlIvnaRB0osEQOYw-hNs0DIhwqJfzeWP30_h8EUvh7wFepXR8ah5rnhBeYdHIT3CQYQt-0EvaZb6gv4WMO22hJ1nGSRy7kdF6MCrEiwg/s1600/bayes2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="421" data-original-width="1102" height="76" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhlje1SvRaMn40jCaaKRu7E52JPl9zKsgSQ0pRBXlIvnaRB0osEQOYw-hNs0DIhwqJfzeWP30_h8EUvh7wFepXR8ah5rnhBeYdHIT3CQYQt-0EvaZb6gv4WMO22hJ1nGSRy7kdF6MCrEiwg/s200/bayes2.png" width="200" /></a></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
It can be simplify as:</div>
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style><br />
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh3Vsa-w-LNYsoVpTU5_7JBqX0hM4bCWnWD7Wfb9g8craI_6Izgyvyc23_PKvwGKmKQLLIhftZmnJAIOcaJiJ10Nsyd568RtGzwJxNj9xEWoK26wT2hwbPX9KJsztSnpKps6lMebh4YTHP-/s1600/bayes3.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="267" data-original-width="831" height="63" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh3Vsa-w-LNYsoVpTU5_7JBqX0hM4bCWnWD7Wfb9g8craI_6Izgyvyc23_PKvwGKmKQLLIhftZmnJAIOcaJiJ10Nsyd568RtGzwJxNj9xEWoK26wT2hwbPX9KJsztSnpKps6lMebh4YTHP-/s200/bayes3.png" width="200" /></a></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
The advantage of this method is that we avoid to compute the marginal
likelihood, that is often difficult to obtain with more complex
models. Let’s stop here a little bit to explain each term of this
equation.</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
Since our main model
is a binomial model (coin toss), the likelihood function Pr(<i>h</i>|<i>p</i>)
can be defined in R language as: </div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> # help(dbinom)
likelihood <- function(h, n, p){
lh <- dbinom(h, n, p)
lh
}
</code></pre>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
Pr(<i>p</i>) is our prior probability, that is, our believe about
what values can take <i>p</i><span style="font-style: normal;">. The
only thing that we know is that it must be a value between 0 and 1,
since it is a probability. If we think that all values have the same
probability, we can define a flat prior using the </span><i>beta</i><span style="font-style: normal;">
distribution. A <i>beta</i> distribution with parameters </span><span style="font-family: "liberation" serif , serif;"><span style="font-style: normal;">β</span></span><span style="font-style: normal;">(1,1)
is a flat distribution between 0 and 1 (you can learn more about <i>beta</i>
distribution <a href="https://en.wikipedia.org/wiki/Beta_distribution" target="_blank">here</a>). </span><span style="font-style: normal;">Then, in
R we can obtain Pr(<i>p</i>) under a </span><span style="font-family: "liberation" serif , serif;"><span style="font-style: normal;">β</span></span><span style="font-style: normal;">(1,1)
</span><span style="font-style: normal;">as: </span></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style> </div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> # help(dbeta)
dbeta(p, 1, 1)
</code></pre>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
Now, the acceptance probability (R, see equations in Step 3) will be the minimum value: 1 or the ratio of
posterior probabilities given the different <i>p</i>.
We express this equation in R language as follows:</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style> </div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> R <- likelihood(h,n,p_prime)/likelihood(h,n,p) * (dbeta(p_prime,1,1)/dbeta(p,1,1))
</code></pre>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
- Step 4) Next, we generate a uniform random number between 0 and 1. If this
number is < R, we will accept the new value for <i>p</i> (<i>p</i>') and we
update the value of <i>p</i> = <i>p</i>'. If not, the change is rejected.</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> if (R > 1) {R <- 1}
random <- runif (1, 0, 1)
if (random < R) {
p <- p_prime
}
</code></pre>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
- Step 5) Now we record the current value of <i>p</i>.</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> posterior[i,1] <- log(likelihood(h, n, p))
posterior[i,2] <- p
</code></pre>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
Finally, we should repeat this loop many times to obtain a good estimate of <i>p</i><span style="font-style: normal;">.
This can be easily done in R using a </span><i>For</i><span style="font-style: normal;">
loop (check the full code below).</span></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style> </div>
<div style="line-height: 100%; margin-bottom: 0cm;">
</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
Following the <a href="https://revbayes.github.io/tutorials/mcmc_binomial/" target="_blank">example</a> studied in MadPhylo with RevBayes, I wrote some
code to reproduce this example with R. We toss a coin 100 times, and
we obtain 73 heads. Is my coin biased? Let's check what is the
probability to obtain a head given the data:</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> # Set the numer of tosses.
n <- 100
# Set the number of heads obtained.
h <- 73
# Define our likelihood function.
# Since our model is a binomial model, we can use:
likelihood <- function(h,n,p){
lh <- dbinom(h,n,p)
lh
}
# Set the starting value of p
p <- runif(1,0,1)
# Create an empty data.frame to store the accepted p values for each iteration.
# Remember: "the posterior probability is just an updated version of the prior"
posterior <- data.frame()
# Set the lenght of the loop (Marcov Chain, number of iterations).
nrep <- 5000
# Start the loop (MCMC)
for (i in 1:nrep) {
# Obtain a new proposal value for p
p_prime <- p + runif(1, -0.05,0.05)
# Avoid values out of the range 0 - 1
if (p_prime < 0) {p_prime <- abs(p_prime)}
if (p_prime > 1) {p_prime <- 2 - p_prime}
# Compute the acceptance proability using our likelihood function and the
# beta(1,1) distribution as our prior probability.
R <- likelihood(h,n,p_prime)/likelihood(h,n,p) * (dbeta(p_prime,1,1)/dbeta(p,1,1))
# Accept or reject the new value of p
if (R > 1) {R <- 1}
random <- runif (1,0,1)
if (random < R) {
p <- p_prime
}
# Store the likelihood of the accepted p and its value
posterior[i,1] <- log(likelihood(h, n, p))
posterior[i,2] <- p
print(i)
}
</code></pre>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
Let's plot some results...<br />
<br />
<br /></div>
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> par(mfrow= c(1,2))
prior <- rbeta(5000, 1,1)
plot(1:5000 ,posterior$V2, cex=0, xlab = "generations", ylab = "p",
main = "trace of MCMC\n accepted values of parameter p\n prior = beta(1,1) generations = 5000")
lines(1:5000, posterior$V2, cex=0)
abline(h=mean(posterior$V2), col="red")
plot(density(posterior$V2), xlim = c(min(min(prior),min((posterior$V2))), max(max(prior),max((posterior$V2)))),
ylim = c(0, max(max(density(prior)$y),max((density(posterior$V2)$y)))), main= "prior VS posterior\n prior= beta(1,1)",
lwd=3, col="red")
lines(density(prior), lwd=3, lty=2, col="blue")
legend("topleft", legend=c("prior density","posterior density"),
col=c("blue","red"), lty=c(3,1), lwd=c(3,3), cex = 1)
</code></pre>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi9U1IEMCYzMk9SaHO-WrFJRdO38QNX7ejUE5Ph9VuSGp7phAtC9WuZ-0eAGojE1U9uVzzTY7vZylKvgH2zK0pkUnihQc0qOU-z3cXhYuV1DKI9EOaMdoRdrejE2XERs4vtf4zBdXTep322/s1600/both.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="588" data-original-width="1323" height="177" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi9U1IEMCYzMk9SaHO-WrFJRdO38QNX7ejUE5Ph9VuSGp7phAtC9WuZ-0eAGojE1U9uVzzTY7vZylKvgH2zK0pkUnihQc0qOU-z3cXhYuV1DKI9EOaMdoRdrejE2XERs4vtf4zBdXTep322/s400/both.png" width="400" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
Well, we can see that the probability to obtain a head given our data is around 0.7, so our coin must be a fake! </div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
You can play with the code and explorewith a different number of tosses, or the effect of a different prior for <i>p</i>...</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEguiyav_dIf8XGtNKJ-nbjf_CLhHLfJ_oTh-Z9L6mvqChWhtTN3RDctIPtKFdBpz2Ezspa7yjO_OMYx7xGWafkuUEBCJPyjojgjSuX2Us6e2qPpPs-u4el21iOWUobYZKOHETe2NJ7OrRHS/s1600/plot_4.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="805" data-original-width="1319" height="243" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEguiyav_dIf8XGtNKJ-nbjf_CLhHLfJ_oTh-Z9L6mvqChWhtTN3RDctIPtKFdBpz2Ezspa7yjO_OMYx7xGWafkuUEBCJPyjojgjSuX2Us6e2qPpPs-u4el21iOWUobYZKOHETe2NJ7OrRHS/s400/plot_4.png" width="400" /></a></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
If you want to learn more about Bayesian Inference, I recommend you these YouTube <a href="https://www.youtube.com/watch?v=3OJEae7Qb_o&t=4s" target="_blank">videos</a> by <a href="http://www.sumsar.net/about.html" target="_blank">Rasmus Bååth</a>, or this wonderful <a href="https://github.com/rmcelreath/statrethinking_winter2019" target="_blank">book/course</a> by <a href="https://xcelab.net/rm/" target="_blank">Richard McElreath</a></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
That's all!<br />
<br />
Enjoy!</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style></div>
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style><br />
<div style="line-height: 100%; margin-bottom: 0cm;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style></div>
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style>Javier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.com2tag:blogger.com,1999:blog-9045240617730855574.post-89217386810789248992018-11-25T20:31:00.001-08:002020-12-10T02:07:47.309-08:00Plotting wind highways using rWind<script type='text/javascript' src='https://d1bxh8uas1mnw7.cloudfront.net/assets/embed.js'></script>
<br />
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<br /></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
Hi there!</div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<br /></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
Our manuscript about rWind R package has been recently accepted for
publication in <a href="https://onlinelibrary.wiley.com/doi/pdf/10.1111/ecog.03730" target="_blank">Ecography</a>! As you know, rWind is a tool used to
download and manage wind data, with some utilities that make easy to
include wind information in ecological or evolutionary analyses (<a href="https://www.tandfonline.com/doi/abs/10.1080/15230406.2017.1403376" target="_blank">or others!</a>).</div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<br /></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
Though there are several examples in the Supporting material of the
publication, and others in the <a href="https://cran.r-project.org/web/packages/rWind/vignettes/Shortest_paths.html" target="_blank">vignette</a> included in the CRAN version,
here you have a full tutorial with the main utilities of rWind...
enjoy!</div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<br /></div>
<h3 lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<u><b>Computing wind connectivity of Azores Archipelago</b></u></h3>
<div style="text-align: justify;">
<u><b>
</b></u></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<u><b>
</b></u></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
</div>
<div align="justify" lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
<br />
In this example, we will compute wind connectivity of Azores Islands
with different points of Europe, Africa and North America during June
2018. I will divide each task of this tutorial in different sections
for clarification.</div>
<div align="justify" lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgrU6dQmfw1DVr-Eit1S6nxP-5LpxQl-Z3JCqYrcZDTpS85YTcTa0tCpG_HrbUScG7_pPGA6mdDVGb6eFkfmZaIcppsdsPkqPMuDGl7__XuhA7tIfntHDnZAbQoPlzvinUR3_riEAwjx3rV/s1600/paths.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="725" data-original-width="1101" height="260" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgrU6dQmfw1DVr-Eit1S6nxP-5LpxQl-Z3JCqYrcZDTpS85YTcTa0tCpG_HrbUScG7_pPGA6mdDVGb6eFkfmZaIcppsdsPkqPMuDGl7__XuhA7tIfntHDnZAbQoPlzvinUR3_riEAwjx3rV/s400/paths.png" width="400" /></a></div>
</div>
<h4 align="justify" lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
</h4>
<h4 align="justify" lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
Downloading wind time series using lubridate and wind.dl_2()</h4>
<div align="justify" lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div align="justify" lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
In the last version of rWind (1.0.4), we can use lubridate to create a list of dates which will be used latter with wind.dl_2 to download wind time series. The object obtained with wind.dl_2 is an "rWind_series" "list". This object can be easily used with other rWind functions as we'll see later.</div>
<div align="justify" lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
</div>
<div align="justify" lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<pre style="background: rgb(240, 240, 240) none repeat scroll 0% 0%; border: 1px dashed rgb(204, 204, 204); color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; overflow-wrap: normal; word-wrap: normal;"> # Make sure that you have the last version of rWind
devtools::install_github("jabiologo/rWind")
library(rWind)
library(lubridate)
# Here, we use ymd_hms from lubridate package to create a sequence of dates
dt <- seq(ymd_hms(paste(2018,6,1,12,00,00, sep="-")),
ymd_hms(paste(2018,6,30,12,00,00, sep="-")),by="1 days")
# Now we can use wind.dl_2 with this sequence of dates. Have into account
# that it could take a while, they are 30 datasets and it's a big area.
ww <- wind.dl_2(dt,-85,5,20,60)
</code></pre>
<div align="justify" lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
</div>
<div align="justify" lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div align="justify" lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
<br />
<h4 align="justify" lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
Obtaining maximum wind speed during the time series with dplyr</h4>
<div align="justify" lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
</div>
<div align="justify" lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
<br />
We can use the rWind function called "tidy" to play with the "rWind_series" "list" object to obtain, for example, the
maximum speed reported for each cell in the study area during the
time series.<br />
<br />
<pre style="background: rgb(240, 240, 240) none repeat scroll 0% 0%; border: 1px dashed rgb(204, 204, 204); color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; overflow-wrap: normal; word-wrap: normal;"> # After tidy our wind data, we can use dplyr utilities to easily
# obtain several metrics from our wind data. Notice that wind average
# is a special case, and it will be discussed later
t_ww <- tidy(ww)
library(dplyr)
g_ww <- t_ww %>% group_by(lat, lon)
# Now, we select only the maximum speed values for each coordinate
max_ww <- g_ww %>% summarise(speed = max(speed))
maxw <- cbind(max_ww$lon, max_ww$lat, max_ww$speed)
head(maxw)
# We can also check the maximum speed value in the whole study area
max(maxw[,3])
# Finally, we can transform this data into a raster layer
rmax <- rasterFromXYZ(maxw)
acol <- colorRampPalette(c("white", "blue", "darkblue"))
plot(rmax, col=acol(1000), main= "Maximum wind speed reported",
xlab="Longitude", ylab="Lattitude")
library(rworldmap)
lines(getMap(resolution = "high"), lwd=2)
</code></pre>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgzIei0Ug58cVGeg7MIEdG650fPd_YzG9G7ZCLlB7XEUekI7Nb19jaybhHthNr-wl8ovYpOH0fkzCHy3SXYwiibfXtVx-ZouFP1hQxKrsXAMOdmWnUnBvMp_6DSnUyO5RR2R9vfkTlh5urz/s1600/max_speed.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="725" data-original-width="1101" height="262" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgzIei0Ug58cVGeg7MIEdG650fPd_YzG9G7ZCLlB7XEUekI7Nb19jaybhHthNr-wl8ovYpOH0fkzCHy3SXYwiibfXtVx-ZouFP1hQxKrsXAMOdmWnUnBvMp_6DSnUyO5RR2R9vfkTlh5urz/s400/max_speed.png" width="400" /></a></div>
<br />
<br />
<h4 align="justify" lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
</h4>
<h4 align="justify" lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
Computing averages of wind speeds and directions </h4>
<div align="justify" lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div align="justify" style="line-height: 100%; margin-bottom: 0cm;">
<span lang="en-US">Since
wind data is a vectorial magnitude given in two components (U and V),
to compute speed and direction averages we should use both vectors. To make it easier, rWind includes wind.mean function.
Notice that wind.mean is specially </span><span lang="en-US">developed
to work with wind.dl_2 outputs, so it will require "rWind_series"
object as input.</span></div>
</div>
</div>
<div align="justify" lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
<pre style="background: rgb(240, 240, 240) none repeat scroll 0% 0%; border: 1px dashed rgb(204, 204, 204); color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; overflow-wrap: normal; word-wrap: normal;"> # We can easily use wind.mean() function
w_mean <- wind.mean(ww)
# And then use wind2raster directly to create a raster layer
r_mean <- wind2raster(w_mean)
# We can plot a subset of this raster around Azores Islands.
# Using "arrowDir" we can include arrows for wind direction with
# "Arrowhead" function from "shape" R package
plot(r_mean$wind.speed, main= "Wind speed and direction average",
col=acol(1000), xlim=c(-38,-18), ylim=c(33,43))
alpha <- arrowDir(w_mean)
library(shape)
Arrowhead(w_mean$lon, w_mean$lat, angle=alpha,
arr.length = 0.2, arr.type="curved")
lines(getMap(resolution = "low"), lwd=4)
</code></pre>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg7hUZPRl2cWdr8Dju53pv2WhuyWj8dP_AbfQWpsu2HMYR8i24yUVZK2iNSbaNdGQA7TJrHWCOSI1M7ObL_J4WIRMQAhfn3V28IhBU2ZJPne1Vh1r-gtSNnByTZ3lhsuAUBddeVEl12NdBs/s1600/azores.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="725" data-original-width="1021" height="283" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg7hUZPRl2cWdr8Dju53pv2WhuyWj8dP_AbfQWpsu2HMYR8i24yUVZK2iNSbaNdGQA7TJrHWCOSI1M7ObL_J4WIRMQAhfn3V28IhBU2ZJPne1Vh1r-gtSNnByTZ3lhsuAUBddeVEl12NdBs/s400/azores.png" width="400" /></a></div>
<br />
<h4 align="justify" style="line-height: 100%; margin-bottom: 0cm;">
<span lang="en-US">W</span><span lang="en-US">ind connectivity
between Azores and mainland</span></h4>
<div align="justify" style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div align="justify" style="line-height: 100%; margin-bottom: 0cm;">
<span lang="en-US">To
complete this tutorial, we will use the downloaded wind data to
compute wind connectivity between several
points in mainland </span><span lang="en-US"><span lang="en-US">and Azores Archipelago</span>. This kind of measures can be useful in order to
establish theories about island colonization by plants or other
organisms that could be influenced by wind currents. We
will also use this connectivity to compute least cost paths between
two locations, in order to virtually plot a wind highway connecting
mainland with Azores islands.</span></div>
<div align="justify" style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<pre style="background: rgb(240, 240, 240) none repeat scroll 0% 0%; border: 1px dashed rgb(204, 204, 204); color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; overflow-wrap: normal; word-wrap: normal;"> # First, we define our locations. 7 locations in the East
# and 7 in the West. The last one is a point in the Azores islands
loc <- matrix(c(-7, -4, -1, -9, -6,-10,-15,
-60, -56,-63,-76,-76,-81,-80, -27,
55, 50, 45, 40, 35, 30, 25,
55, 50, 45, 40, 35, 30, 25, 39), 15, 2)
colnames(loc) <- c("lon", "lat")
rownames(loc) <- c(seq(1,14,1), "Azores")
# Check how "loc" is looking
tail(loc)
# You can plot them
#plot(r_mean$wind.speed, col=acol(1000))
#lines(getMap(resolution = "low"), lwd=2)
#points(loc, col="red4", pch = 19, cex = 2)
</code></pre>
</div>
<div align="justify" lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
<br />
<br />
<div align="justify" style="line-height: 100%; margin-bottom: 0cm;">
<span lang="en-US">Now, we'll execute the next </span><span lang="en-US">tasks</span><span lang="en-US">:</span></div>
<div align="justify" style="line-height: 100%; margin-bottom: 0cm;">
</div>
<ol>
<li><span lang="en-US">We'll
use wind2raster to obtain raster l</span><span lang="en-US">a</span><span lang="en-US">yers
of wind direction and speed of each day in the time series. </span></li>
<li><span lang="en-US">Then,
for each day, we</span><span lang="en-US">'</span><span lang="en-US">ll
compute the conductance matrix, </span><span lang="en-US">a matrix
with connectivity values between every single cell of the entire
study area.</span></li>
<li><span lang="en-US">Later,
with each conductance matrix, we'll compute the cost </span><span lang="en-US">between
</span><span lang="en-US">our</span><span lang="en-US"> locations</span><span lang="en-US">
and </span><span lang="en-US">will </span><span lang="en-US">store
them in a cost_list object.</span></li>
<li><span lang="en-US">Finally,
for each day, i</span><span lang="en-US">f the Cost is not Inf</span><span lang="en-US">inite</span><span lang="en-US">,
</span><span lang="en-US">we'll </span><span lang="en-US">get the
shortest path between </span><span lang="en-US">two </span><span lang="en-US">selected
locations, and </span><span lang="en-US">we'll</span><span lang="en-US">
store </span><span lang="en-US">it</span><span lang="en-US"> into a
“paths” object.</span></li>
</ol>
<br />
<pre style="background: rgb(240, 240, 240) none repeat scroll 0% 0%; border: 1px dashed rgb(204, 204, 204); color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; overflow-wrap: normal; word-wrap: normal;"> # 1) Raster layers from wind time series
layers <- wind2raster(ww)
# 2) Conductance from wind direction and speed
Conductance <- flow.dispersion(layers, type="passive",
output="transitionLayer")
# 3) and 4) Using a loop to compute cost between
# locations and least cost paths between location 11
# and Azores islands (as an example).
cost_list <- array(NA_real_, dim=c(15,15,30))
paths <- list(1:30)
library(gdistance)
for (i in 1:30){
cost_list[,,i] <- costDistance(Conductance[[i]],loc)
if (costDistance(Conductance[[i]],loc[11,] , loc[15,]) != Inf){
paths[[i]] <- shortestPath(Conductance[[i]], loc[11,], loc[15,],
output="SpatialLines")
print(i)
}
}
</code></pre>
<br />
<br />
<div align="justify" style="line-height: 100%; margin-bottom: 0cm;">
<span lang="en-US">Now, we will manage a little bit the results to
plot them in a cool way.</span></div>
<div align="justify" style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div align="justify" style="line-height: 100%; margin-bottom: 0cm;">
<span lang="en-US"><br /></span></div>
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style><br />
<pre style="background: rgb(240, 240, 240) none repeat scroll 0% 0%; border: 1px dashed rgb(204, 204, 204); color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; overflow-wrap: normal; word-wrap: normal;"> # Transform costs in connectivity
connectivity <- 1/cost_list
# Obtain connectivity averages
conn_avg <- apply(connectivity, c(1, 2), mean, na.rm = TRUE)
# Here we filter only connectivity from mainland to Azores
# As connectivity has arbitrary magnitude and units, we will
# multiply connectivity values 1000 times to avoid decimal zeros.
conn_azores <- conn_avg[,15]
conn_azores <- conn_azores *1000
# Finally, we can plot connectivity values relative to the width
# of the connections between mainland and Azores. You can power
# the values to increase the differences between them.
plot(r_mean$wind.speed, main= "Connectivity from mainland to Azores",
col=acol(1000))
lines(getMap(resolution = "low"), lwd=1)
for (i in 1:14){
lines(rbind(loc[i,], loc[15,]), lwd = conn_azores[i]^3, col="red4")
}
points(loc, col="red", pch = 19, cex = 1.5)
text( loc+2, labels=c(round(conn_azores[1:14], 2),"Azores"),
col="white", cex= 1.5, font=2)
</code></pre>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEieg78HENaz80sqwDo6_AEyZS9oaBdbR3MOHyJwVW1p8Gen1nfNQ_Ljy3opgwCYfx84FxT7ct3QW32rF3y3HKr6C9NkqcYzOSjWOxHtA7tC5uFkTv6Q7UvWxyrVBJuUdBY8xuekjTAAJdy4/s1600/connectivity.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="667" data-original-width="1021" height="261" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEieg78HENaz80sqwDo6_AEyZS9oaBdbR3MOHyJwVW1p8Gen1nfNQ_Ljy3opgwCYfx84FxT7ct3QW32rF3y3HKr6C9NkqcYzOSjWOxHtA7tC5uFkTv6Q7UvWxyrVBJuUdBY8xuekjTAAJdy4/s400/connectivity.png" width="400" /></a></div>
Last but not least, we can use the shortest paths stored for each day to plot a density heatmap. This can help us to visualize wind corridors across our study area (our wind highways!).<br />
<br />
<pre style="background: rgb(240, 240, 240) none repeat scroll 0% 0%; border: 1px dashed rgb(204, 204, 204); color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; overflow-wrap: normal; word-wrap: normal;"> # We'll use spatstat library to convert the shortest
# paths lines into an psp object.
library(spatstat)
paths_clean <- paths[!sapply(paths, is.null)]
paths_merged <- paths_clean[[1]]
for (h in 2:length(paths_clean)) {
paths_merged <- rbind(paths_merged,paths_clean[[h]])
}
paths_psp <- as(paths_merged, "psp")
# Now, we can obtain a kernel density of the paths
# and convert it into a raster layer that can be plotted.
# You can play with sigma value to modify the smoothing.
lines_kernel <- density(paths_psp, sigma=1, dimyx=c(350,410))
kernel <- raster(lines_kernel)
#kernel[kernel<(maxValue(kernel)*0.1)]<-NA
ext <- extent(r_mean$wind.speed)
kernel <- extend(kernel, ext)
kernel[is.na(kernel)] <- minValue(kernel)
acol <- colorRampPalette(c("grey40","darkblue","red2",
"orange","yellow","white"))
plot(kernel, col=acol(1000), zlim=c(-0.01,5),
main = "Least cost paths density", xlab="Longitude",
ylab="Lattitude")
lines(getMap(resolution = "high"), lwd=2)
</code></pre>
<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjWCWI6lG_u_nhOfLC3LJAW_quAO2zpx-0UL9m3c_BWqAIiT1PInUmivdQfQrn-Ga7rY5Iv7H6HgO4pI8NmstrErsfwu2FVUh2_sGMz0KdxlUwKopAD4zYoutScq-DpNUoEBoOn4fSSf7rh/s1600/paths.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="725" data-original-width="1101" height="262" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjWCWI6lG_u_nhOfLC3LJAW_quAO2zpx-0UL9m3c_BWqAIiT1PInUmivdQfQrn-Ga7rY5Iv7H6HgO4pI8NmstrErsfwu2FVUh2_sGMz0KdxlUwKopAD4zYoutScq-DpNUoEBoOn4fSSf7rh/s400/paths.png" width="400" /></a></div>
<br />
<br />
<br />
Isn't it cool?<br />
<br />
That's all for now!<br />
<br />
<br />
<div class='altmetric-embed' data-badge-type='donut' data-doi="110.1111/ecog.03730"></div>
<br />
<br /></div>
<div align="justify" lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style></div>
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; background: transparent none repeat scroll 0% 0%; }</style>Javier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.com7tag:blogger.com,1999:blog-9045240617730855574.post-71507856553401555652018-07-23T12:40:00.001-07:002018-07-25T12:16:07.773-07:00Antipredator behavior with R (or why wildebeest should stay together)<div class="separator" style="clear: both; text-align: center;">
</div>
<br />
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
Hi there! </div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
In 2010, when I was
studying my Biology Degree at Universidad Complutense in Madrid, I
fell in love with a documentary miniseries called <a href="https://en.wikipedia.org/wiki/Great_Migrations" target="_blank">Great Migrations</a>
(<a href="http://www.nationalgeographic.com.au/tv/great-migrations/" target="_blank">National Geographic</a>). 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 <a href="https://www.uoguelph.ca/ib/fryxell" target="_blank">John Fryxell</a>, 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 (<a href="https://www.youtube.com/watch?v=Fa0rxb_LbfM" target="_blank">here</a> 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.</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjeWJbS_X2Lsbym2Ku_75G2JtIB0V31UNsrqlxcOtykw8z2uC6l3f6ybMV2DQsA2AGImHOKdhGPWFEdNzE4Pkl_-zCWDfm-sjSsudiAzUpYov28NdluvISWcAelwYiSqw_vC3X6sVTwzb9-/s1600/random_cluster.gif" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="690" data-original-width="1300" height="210" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjeWJbS_X2Lsbym2Ku_75G2JtIB0V31UNsrqlxcOtykw8z2uC6l3f6ybMV2DQsA2AGImHOKdhGPWFEdNzE4Pkl_-zCWDfm-sjSsudiAzUpYov28NdluvISWcAelwYiSqw_vC3X6sVTwzb9-/s400/random_cluster.gif" width="400" /></a></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
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.</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
You can find the code to perform these simulations with R below. To perform the gif
of this post, you can follow <a href="http://allthiswasfield.blogspot.com/2017/12/p-margin-bottom-0.html" target="_blank">this procedures</a>. And if you want to be
updated with my posts, find me on my new Twitter account <a href="https://twitter.com/javi_ferlop" target="_blank">@javi_ferlop</a></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
That is all, enjoy!</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; }a:link { }</style>
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> # 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()
}
</code></pre>
Javier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.com3tag:blogger.com,1999:blog-9045240617730855574.post-43866926604749368722018-07-04T16:13:00.000-07:002018-07-04T16:16:23.280-07:00New rWind release on CRAN! (v1.0.2) <br />
<br />
<div style="line-height: 100%; margin-bottom: 0cm;">
Hi there!</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
Just a few lines to
inform you about the new release of rWind R package (v1.0.2). This
version have several new features. Here you have an example of one of
them, the function wind.dl_2 to download time series of wind data. In
this example we create a bunch of PNG files to be converted later
into a GIF of wind speed of east Asia region.</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
Enjoy!</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh7ZPt_e2LSnLyauCwI-otcHeAq7yJZTht1CuOpuW0HYL45e1SQUr8jnEiAyY4obhsk-udtZ_wgryiiySsSBRDDQxgqlUld1Ub6knrHmIZID5OhDpRXv2-rjrjJLXNMSIfXb3CA2L45Akcp/s1600/wind.gif" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="600" data-original-width="1000" height="240" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh7ZPt_e2LSnLyauCwI-otcHeAq7yJZTht1CuOpuW0HYL45e1SQUr8jnEiAyY4obhsk-udtZ_wgryiiySsSBRDDQxgqlUld1Ub6knrHmIZID5OhDpRXv2-rjrjJLXNMSIfXb3CA2L45Akcp/s400/wind.gif" width="400" /></a></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 115%; }</style>
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> # We will create a gif of wind speed for east Asia during the end of June (2018)
# First, we should load the packages that we will use. Use install.packages() to
# download and install the new version of rWind (v1.0.2)
install.package("rWind")
library(rWind)
library(fields)
library(shape)
library(rworldmap)
library(lubridate)
# Now, we use lubridate package to create a sequence of dates/times (each three
# hours)
dt <- seq(ymd_hms(paste(2018,6,25,00,00,00, sep="-")),
ymd_hms(paste(2018,7,4,21,00,00, sep="-")),by="3 hours")
# Now, we use the new function wind.dl_2 to download the whole time series of
# wind data. We use the "dt" object created with lubridate to provide the input
# to wind.dl_2. Since it's a large area and many days, it could take a while...
wind_series <- wind.dl_2(dt,90,150,5,40)
# Next, we can use wind2raster function from rWind package directly over the
# output from wind.dl_2, it has been adapted to work with this lists of wind
# data.
wind_series_layer <- wind2raster(wind_series)
# Finally, we will create a bunch of PNG files to be converted later in a GIF
id<-0
for (i in 1:72) {
id <- sprintf("%03d", i)
png(paste("asia",id,".png", sep=""), width=1000, height=600, units="px",
pointsize=18)
image.plot(wind_series_layer[[i]]$wind.speed, col=bpy.colors(1000),
zlim=c(0,18), main =wind_series[[i]]$time[1])
lines(getMap(resolution = "low"), lwd=3)
dev.off()
}
</code></pre>
Javier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.com4tag:blogger.com,1999:blog-9045240617730855574.post-80531123512014361832018-05-01T12:33:00.000-07:002018-05-01T12:33:16.458-07:00Simulating animal movements and habitat use<br />
Hi there!<br />
<br />
I was training some ways to simulate animal (or other organisms) movements having into account habitat suitability. To do this, I used my previous eWalk model as the underlying process to simulate random or directional walks. This model is based on Brownian / Ornstein–Uhlenbeck process. You can find <a href="http://allthiswasfield.blogspot.com.es/2016/11/ewolk-r-function-to-simulate-two.html" target="_blank">more about eWalk model here</a>!<br />
<br />
Today, I will add one more element to this movement simulations. In this case, we will have into account the habitat or environmental preferences of the simulated species, to perform a simulation like this:<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhpQiylC2glBGF3vqVBepXn3jPdwL1oHSdrFJPIhKmoQrmssS8IMiZ6knskNZD5Cnc5rkjWkWZjvqBCJAiaN8YddGyU2HZV3sjJRnvaefXFFfOZ5dZYTyNQlmktV4c5NaRbqNi9Bsbmahmn/s1600/paths.gif" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="831" data-original-width="800" height="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhpQiylC2glBGF3vqVBepXn3jPdwL1oHSdrFJPIhKmoQrmssS8IMiZ6knskNZD5Cnc5rkjWkWZjvqBCJAiaN8YddGyU2HZV3sjJRnvaefXFFfOZ5dZYTyNQlmktV4c5NaRbqNi9Bsbmahmn/s320/paths.gif" width="308" /></a></div>
<br />
<br />
First, we will create a raster layer as a random environmental variable, for example tree cover.<br />
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">
library (raster)
library (dismo)
tc <- raster(nrows=100, ncols=100, xmn=0, xmx=100, ymn=0,ymx=100)
tc[] <- runif(10000, -80, 180)
tc <- focal(tc, w=matrix(1, 5, 5), mean)
tc <- focal(tc, w=matrix(1, 5, 5), mean)
plot(tc)
</code></pre>
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi71tkRkDeLPpA67op-jGdFLFfjFXfgMmtFE8_KYKuPDortXIfHhXyYJmBUBsr-EpDDamOLAKjTcxlqNFP69vaf6znWpbbp9jfBtzhECdGwKanAinca0Mp3MozS-hQE0ERL0B5rLxzkAxn2/s1600/tc.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="831" data-original-width="800" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi71tkRkDeLPpA67op-jGdFLFfjFXfgMmtFE8_KYKuPDortXIfHhXyYJmBUBsr-EpDDamOLAKjTcxlqNFP69vaf6znWpbbp9jfBtzhECdGwKanAinca0Mp3MozS-hQE0ERL0B5rLxzkAxn2/s400/tc.png" width="383" /></a></div>
Second, we will define the species class. The species will be defined by their position (coordinates), and their optimum for the environmental variable.<br />
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">
species <- setClass("species", slots=c(x="numeric", y="numeric", opt="numeric"))
</code></pre>
<br />
<div class="separator" style="clear: both; text-align: center;">
</div>
Here we will define the Red deer species as a specimen in the coordinates (50,50) and an optimum of 80 for the environmental variable (tree cover). In the same way, we will define the Egyptian mongoose as a specimen in the coordinates (50,50) and an optimum of 30 for the tree cover variable.<br />
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">
Red_deer <- species(x= 50, y =50, opt= 90)
Egyptian_mongoose <- species(x= 50, y =50, opt= 30)
</code></pre>
<br />
Now, we will load the "<i>go</i>" function (I do not have a name for it yet). It require a species (<i>sp</i>), a raster layer with any environmental variable (<i>env</i>), number of iterations (<i>n</i>), a Brownian motion parameter (that is, how random is the movement of your species), a geographical optimum (the wanted destination of your species <i>theta_x</i> and <i>theta_y</i>), and the attraction strength or "interest" of the species to get this position (<i>alpha_x</i> and <i>alpha_y</i>). The syntaxis should be something like this:<br />
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">
path <- go (sp, env, n, sigma, theta_x, alpha_x, theta_y, alpha_y)
</code></pre>
<br />
Here is the function to load (I will comment the function in a future post):<br />
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">
go <- function (sp, env, n, sigma, theta_x, alpha_x, theta_y, alpha_y) {
track <- data.frame()
track[1,1] <- sp@x
track[1,2] <- sp@y
for (step in 2:n) {
neig <- adjacent(env,
cellFromXY(env, matrix(c(track[step-1,1],
track[step-1,2]), 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(options[abs(na.omit(options$V2)) == min(abs(na.omit(options$V2))), 1 ],
options[abs(na.omit(options$V2)) == min(abs(na.omit(options$V2))), 1 ])
new_cell <- sample(option,1)
new_coords <- xyFromCell(env,new_cell)
lon_candidate<--9999
lat_candidate<--9999
while ( is.na(extract(env, matrix(c(lon_candidate,lat_candidate),1,2)))) {
lon_candidate <- new_coords[1]+ (sigma * rnorm(1)) + (alpha_x * ( theta_x - new_coords[1]))
lat_candidate <- new_coords[2]+ (sigma * rnorm(1)) + (alpha_y * ( theta_y - new_coords[2]))
}
track[step,1] <- lon_candidate
track[step,2] <- lat_candidate
}
return(track)
}
</code></pre>
<br />
<br />
Well, now we can perform a simple experiment with our two specimens. We will simulate random movement of these two species having into account their environmental optimums. The "<i>go</i>" function will return us the track or the path followed by each specimen (coordinates by each step).<br />
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">
deer_simul <- go (Red_deer, tc, 100, 2, 90, 0, 90, 0)
mongoose_simul <- go (Egyptian_mongoose, tc, 100, 2, 90, 0, 90, 0)
</code></pre>
<br />
We can plot the paths...<br />
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">
plot(tc)
lines(deer_simul, lwd=1.5, col="red")
points(deer_simul, cex=0.3, col="red")
lines(mongoose_simul, lwd=1.5, col="blue")
points(mongoose_simul, cex=0.3, col="blue")
legend("topleft", legend=c("deer","mongoose"), col=c("red","blue"),
lty=c(1,1), lwd=c(2,2))
</code></pre>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh4MtpCAGWTwyz9hmBka_r6H3VVfHB-TjoTzqv33QrmJ3g8RvEKrEcieK4HNWyyv3dwsTNcW3ku7UEqDLjzHZIT3eIhWqlcKLhEVQAYLgl1hetmzI9kpQtWBHVBedMeRlsLyTe_DtrtP07t/s1600/paths.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="831" data-original-width="800" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh4MtpCAGWTwyz9hmBka_r6H3VVfHB-TjoTzqv33QrmJ3g8RvEKrEcieK4HNWyyv3dwsTNcW3ku7UEqDLjzHZIT3eIhWqlcKLhEVQAYLgl1hetmzI9kpQtWBHVBedMeRlsLyTe_DtrtP07t/s400/paths.png" width="385" /></a></div>
To test if each species is actually "searching" their environmental optimum, we can extract the environmental values by step for each species and plot their density distributions.<br />
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">
plot(density(extract(tc, deer_simul)),lwd=3, col="red", xlim=c(20,80),
ylim=c(0,max(c(density(extract(tc, deer_simul))$y,
density(extract(tc, mongoose_simul))$y))),
main="locations density distribution", xlab="tree cover")
lines(density(extract(tc, mongoose_simul)),lwd=3, col="blue")
legend("topleft", legend=c("deer","mongoose"), col=c("red","blue"),
lty=c(1,1), lwd=c(3,3))
</code></pre>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh-jxbLXQMlRqzo4D-VmxpF85mkvZ3z-GwZB5Kci7RBBAz6SOaAjtuM-AlUqEmVIL40zBL-fW-umxjynANW0rIoNU3QcZQpI1KCfTCWMimt0L5D19v7BLZ1pl1LHtTIc7pOFmgeLPk_ZMhB/s1600/Rplot.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="831" data-original-width="854" height="388" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh-jxbLXQMlRqzo4D-VmxpF85mkvZ3z-GwZB5Kci7RBBAz6SOaAjtuM-AlUqEmVIL40zBL-fW-umxjynANW0rIoNU3QcZQpI1KCfTCWMimt0L5D19v7BLZ1pl1LHtTIc7pOFmgeLPk_ZMhB/s400/Rplot.png" width="400" /></a></div>
So, we can see that the deer is actually using patches with a higher value of tree cover than the mongoose... our simulation worked! You can use the code of <a href="http://allthiswasfield.blogspot.com.es/2017/12/p-margin-bottom-0.html" target="_blank">this post</a> to perform a GIF like the one above.<br /><br />
That's all I have to say about this for now... In the next posts we will simulate more animal movements and migrations!<br />
<br />
<br />
PD: Let me show you a nice song about the <a href="https://www.youtube.com/watch?v=pQKUgoPZEks" target="_blank">mongoose</a>...<br />
<br />
<br />
<br />Javier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.com2tag:blogger.com,1999:blog-9045240617730855574.post-82642042937386953382018-03-14T18:05:00.002-07:002018-03-14T18:14:44.491-07:00Nice ggplot with sad data: something happens with women in science<br />
<div style="text-align: justify;">
</div>
<div style="text-align: justify;">
Last March 8, millions of women in more than 170 countries around the world joined street protests calling for "a society free of sexist oppression, exploitation and violence". Spanish strike was one of the <a href="http://www.bbc.com/news/world-europe-43324406?ocid=socialflow_facebook&ns_mchannel=social&ns_campaign=bbcnews&ns_source=facebook" target="_blank">most numerous</a>, where around 5.3 million women joined the strike. Today, I will use some data (<a href="https://drive.google.com/open?id=10jJyjcaqtCgBvvc5ed2HnsI3nOUMDPw3" target="_blank">download</a>) regarding the scientific career in Spain (<a href="http://www.csic.es/mujeres-y-ciencia" target="_blank">2016</a>) to make a plot using ggplot2 and discuss about it. Here you have the graph and the code:</div>
<div style="text-align: justify;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhdg2NHIp0fpjguIWA_wxlVdmE5vg7wXgDxeskuvMJWe70Rw1qWZI8CIx7ewDsDKkJdbcr6yJcj8SLpS1NxgkfM_lBZofwfBlQYpmEWZ75wlxexc2Z9gKao7zL70Mi5d7m0tt8G4ghouMgn/s1600/women.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="669" data-original-width="1034" height="258" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhdg2NHIp0fpjguIWA_wxlVdmE5vg7wXgDxeskuvMJWe70Rw1qWZI8CIx7ewDsDKkJdbcr6yJcj8SLpS1NxgkfM_lBZofwfBlQYpmEWZ75wlxexc2Z9gKao7zL70Mi5d7m0tt8G4ghouMgn/s400/women.png" width="400" /></a></div>
<div style="text-align: justify;">
<br /></div>
<div style="text-align: justify;">
</div>
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> library(ggplot2)
library(scales)
dat <- read.csv("wis.csv")
dat$pos <-
factor(dat$pos, levels = c("pre", "post", "ryc", "ct", "ic" , "pi"))
ggplot(dat, aes(x = pos, y = per)) +
geom_line(size = 3, aes(group = sex, color = sex)) +
geom_point() +
geom_text(aes(label = round(per * 100, 1)), vjust = -1, size = 4) +
labs(
title = "Women in Science \n Spanish National Research Council",
x = "research stage",
y = "sex ratio",
color = " \n",
subtitle = "(2016)"
) +
scale_color_manual(labels = c("women", "men"),
values = c("purple", "orange")) +
scale_x_discrete(
breaks = c("pre", "post", "ryc", "ct", "ic" , "pi"),
labels = c(
"Predoctoral",
"Postdoctoral",
"Ramón y Cajal",
"Científico\nTitular",
"Investigador\nCientífico",
"Profesor de\nInvestigación"
)
) +
scale_y_continuous(labels = percent) +
theme(
text = element_text(size = 15),
legend.text = element_text(size = 15),
axis.text.x = element_text(face = "bold", size = 11, angle = 30),
axis.title = element_text(size = 14, face = "bold"),
plot.title = element_text(hjust = 0.5, face = "bold", size = 20),
plot.subtitle = element_text(hjust = 0.5, size = 19)
)
</code></pre>
<br />
<br />
<div style="text-align: justify;">
In X axis, we can see the stages at the Spanish National Research Council (known in Spanish as CSIC) positions, from young PhD candidates (left - “Predoctoral”) to the highest level, Research Professor (right - “Profesor de Investigación”). It’s easy to understand: in Spain, there are more women starting a scientific career than men, but only few of them get to the highest positions in CSIC.<br /><br />We usually call this the “scissors graph”, since the two curves soon cross each other changing the trend shown at the beginning. For the last stages this imbalance between males and females can be due to the age (mean age of Research Professor in CSIC is around 58 years old, and things have changed a lot in Spain since 30 years ago). However, the graph shows that for the earlier stages this imbalance is still being a problem: just from PhD candidate stage to first postdoctoral position, women lose around 15% positions in favor of men. <a href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2247379/" target="_blank">And this is not only happening in Spain</a>. Only active politics and special measures orientated to these first stages could change the trend for the highest positions in the future!<br /><br />If you want to know a little bit more about women in science in Spain, here you have a <a href="https://drive.google.com/open?id=14IZi1MTgkpWK_op4LDb2s7HLWKTcI1tc" target="_blank">poster</a> made by some friends from the Ecology Department in Universidad de Alcalá.<br /></div>
Javier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.com3tag:blogger.com,1999:blog-9045240617730855574.post-24559026351456549692017-12-09T10:18:00.000-08:002017-12-09T11:44:03.961-08:00Brownian Motion GIF with R and ImageMagick<style type="text/css">p { margin-bottom: 0.25cm; line-height: 120%; }</style>
<br />
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 120%; }</style>
<br />
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
Hi
there!</div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<br /></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
Last
Monday we celebrated a “Scientific Marathon” at Royal Botanic
Garden in Madrid, a kind of mini-conference to talk about our
research. I was talking about the relation between fungal spore size and
environmental variables such as temperature and precipitation. To
make my presentation more friendly, I created a GIF to explain the
Brownian Motion model. In evolutionary biology, we can use this model
to simulate the random variation of a continuous trait through time.
Under this model, we can notice how closer species tend to maintain
closer trait values due to shared evolutionary history. You have a
lot of information about Brownian Motion models in evolutionary biology
<a href="https://lukejharmon.github.io/pcm/chapter3_bmintro/">everywhere</a>!</div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
Here I
will show you how I built a GIF to explain Brownian Motion in my talk
using R and ImageMagick.</div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<br /></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<br /></div>
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> # First, we simulate continuous trait evolution by adding in each iteration
# a random number from a normal distribution with mean equal to 0 and standard
# deviation equal to 1. We simulate a total of 4 processes, to obtain at first
# two species and a specieation event at the middle of the simulation, obtaining
# a total of 3 species at the end.
df1<- data.frame(0,0)
names(df1)<- c("Y","X")
y<-0
for (g in 1:750){
df1[g,2] <- g
df1[g,1] <- y
y <- y + rnorm(1,0,1)
}
#plot(df1$X,df1$Y, ylim=c(-100,100), xlim=c(0,1500), cex=0)
#lines(df1$X,df1$Y, col="red")
df2<- data.frame(0,0)
names(df2)<- c("Y","X")
y<-0
for (g in 1:1500){
df2[g,2] <- g
df2[g,1] <- y
y <- y + rnorm(1,0,1)
}
#lines(df2$X,df2$Y, col="blue")
df3<- data.frame(750,df1[750,1])
names(df3)<- c("Y","X")
y<-df1[750,1]
for (g in 750:1500){
df3[g-749,2] <- g
df3[g-749,1] <- y
y <- y + rnorm(1,0,1)
}
#lines(df3$X,df3$Y, col="green")
df4<- data.frame(750,df1[750,1])
names(df4)<- c("Y","X")
y<-df1[750,1]
for (g in 750:1500){
df4[g-749,2] <- g
df4[g-749,1] <- y
y <- y + rnorm(1,0,1)
}
#lines(df4$X,df4$Y, col="orange")
</code></pre>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgUq6sMYTboJGOApSuaDQpjnpZNtRHoaONkK4UH0ALYM7nBCGeKKUzavV7_4HidJo0iJX3M5yPV7cMIJxLLEJg04-BG7G5tenSlFMO9OM6kD5ieJK_PQLvkRfZlKaYyigKDxreM7-l_9wqP/s1600/bm1498.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="570" data-original-width="900" height="252" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgUq6sMYTboJGOApSuaDQpjnpZNtRHoaONkK4UH0ALYM7nBCGeKKUzavV7_4HidJo0iJX3M5yPV7cMIJxLLEJg04-BG7G5tenSlFMO9OM6kD5ieJK_PQLvkRfZlKaYyigKDxreM7-l_9wqP/s400/bm1498.png" width="400" /></a></div>
<br />
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> # Now, we have to plot each simmulation lapse and store them in our computer.
# I added some code to make lighter the gif (plotting just odd generations) and
# to add a label at the speciation time. Note that, since Brownan Model is a
# stocasthic process, my simulation will be different from yours.
# You should adjust labels or repeat the simulation process if you don't
# like the shape of your plot.
parp<-rep(0:1, times=7, each= 15)
parp<- c(parp, rep(0, 600))
for (q in 1:750){
if ( q %% 2 == 1) {
id <- sprintf("%04d", q+749)
png(paste("bm",id,".png", sep=""), width=900, height=570, units="px",
pointsize=18)
par(omd = c(.05, 1, .05, 1))
plot(df1$X,df1$Y, ylim=c(-70,70), xlim=c(0,1500), cex=0,
main=paste("Brownian motion model \n generation=", 749 + q) ,
xlab="generations", ylab="trait value", font.lab=2, cex.lab=1.5 )
lines(df1$X,df1$Y, col="red", lwd=4)
lines(df2$X[1:(q+749)],df2$Y[1:(q+749)], col="blue", lwd=4)
lines(df3$X[1:q],df3$Y[1:q], col="green", lwd=4)
lines(df4$X[1:q],df4$Y[1:q], col="orange", lwd=4)
if (parp[q]==0)
text(750, 65,labels="speciation event", cex= 1.5, col="black", font=2)
if (parp[q]==0)
arrows(750, 60, 750, 35, length = 0.20, angle = 30, lwd = 3)
dev.off()
}
}
</code></pre>
<br />
Now, you just have to use ImageMagick to put all the PNG files together in a GIF using a command like this in a terminal:<br />
<br />
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> convert -delay 10 *.png bm.gif
</code></pre>
<br />
Et voilà!<br />
<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEglV0F5VWR80Lmpgi_b4YvUf3jK-QxTdOT_4TRW6WHV4isUnd_bu7iub0VnXDZ0pLO2dwX7BLXQWlYXfFuglgK0QzLRAJJHrxj7KpvpuAqrSEuGVqOY4R3NOL48rz5LdmFhYt-ul3XUUaot/s1600/bm_fast2.gif" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="570" data-original-width="900" height="251" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEglV0F5VWR80Lmpgi_b4YvUf3jK-QxTdOT_4TRW6WHV4isUnd_bu7iub0VnXDZ0pLO2dwX7BLXQWlYXfFuglgK0QzLRAJJHrxj7KpvpuAqrSEuGVqOY4R3NOL48rz5LdmFhYt-ul3XUUaot/s400/bm_fast2.gif" width="400" /></a></div>
<br />Javier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.com4tag:blogger.com,1999:blog-9045240617730855574.post-58378685109059054502017-12-05T08:54:00.000-08:002017-12-05T13:34:39.841-08:00Computing wind average in an area using rWind<style type="text/css">p { margin-bottom: 0.25cm; line-height: 120%; }</style>
<br />
<div style="line-height: 100%; margin-bottom: 0cm;">
Hi all!</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
A researcher asked
me last week about how to compute wind average for an
area using rWind. I wrote a simple function to do this using dplyr
library (following the advice of my friend <a href="https://www.researchgate.net/profile/Javier_Fajardo4">Javier Fajardo</a>). The function will be added to
rWind package as soon as possible. Meanwhile, you can test the
results... enjoy!</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhEY4kjnGxa1srutDLPmUCOppg6tQi5IIs-y_CQ-dkl_Bcu6G8l0o_AZls3lGCC_vQ1rI_adKRlTkNxb_Ubo4LkmZ4K3d8MwY3-ut0HpdST0gXp59trzHdMEOW3Ikw6KGNjkJkD24yAvZQ2/s1600/Rplot04.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="831" data-original-width="796" height="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhEY4kjnGxa1srutDLPmUCOppg6tQi5IIs-y_CQ-dkl_Bcu6G8l0o_AZls3lGCC_vQ1rI_adKRlTkNxb_Ubo4LkmZ4K3d8MwY3-ut0HpdST0gXp59trzHdMEOW3Ikw6KGNjkJkD24yAvZQ2/s320/Rplot04.png" width="306" /></a></div>
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> # First, charge the new function
library(dplyr)
wind.region <- function (X){
X[,3] <- X[,3] %% 360
X[X[,3]>=180,3] <- X[X[,3]>=180,3] - 360
avg<-summarise_all(X[,-1], .funs = mean)
wind_region <- cbind(X[1,1],avg)
return(wind_region)
}
</code></pre>
<br />
Once you have
charged the function, let’s do an example...<br />
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> # Get some wind data and convert it into a raster to be plotted later
library(rWind)
library(raster)
wind_data<-wind.dl(2015,2,12,0,-10,5,35,45)
wind_fitted_data <- wind.fit(wind_data)
r_speed <- wind2raster(wind_fitted_data, type="speed")
</code></pre>
<br />
Now, you can use the new function to obtain wind average in the study area:<br />
<br />
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> myMean <- wind.region(wind_data)
myMean
</code></pre>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjIm8W3mq9VX_0ixucmFiep1IzcJQzVzSV5GtXOdKjt-OFcI2WHQtDP5dIKTIehVtA575kyM44aFy8CXZrWCcle0S8_AmL_DQ1sZuElSWOUo5wAil5E-Vx9f3VNjNOwsf5WFz00saMLVx9D/s1600/myMean.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="47" data-original-width="733" height="25" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjIm8W3mq9VX_0ixucmFiep1IzcJQzVzSV5GtXOdKjt-OFcI2WHQtDP5dIKTIehVtA575kyM44aFy8CXZrWCcle0S8_AmL_DQ1sZuElSWOUo5wAil5E-Vx9f3VNjNOwsf5WFz00saMLVx9D/s400/myMean.png" width="400" /></a></div>
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> # Now, you can use wind.fit to get wind speed and direction.
myMean_fitted <- wind.fit(myMean)
myMean_fitted
</code></pre>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiu-y8dtecpMDkkoaRekfbh8q-GduhQAn7eLVVXEdyey4hwWHIxtwr_Gs0XGv25zwfUstfYJo1rUf2GB4qbVqf-dwkyxtCJqajtAwvkM5tqHs85eEzGvqh6SCm8Xqucn29sqluFoyucabTw/s1600/myMean2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="50" data-original-width="207" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiu-y8dtecpMDkkoaRekfbh8q-GduhQAn7eLVVXEdyey4hwWHIxtwr_Gs0XGv25zwfUstfYJo1rUf2GB4qbVqf-dwkyxtCJqajtAwvkM5tqHs85eEzGvqh6SCm8Xqucn29sqluFoyucabTw/s1600/myMean2.png" /></a></div>
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> # Finally, let's plot the results!
library(rworldmap)
library(shape)
plot(r_speed)
lines(getMap(resolution = "low"), lwd=4)
alpha <- arrowDir(myMean_fitted)
Arrowhead(myMean_fitted$lon, myMean_fitted$lat, angle=alpha,
arr.length = 2, arr.type="curved")
text(myMean_fitted$lon+1, myMean_fitted$lat+1,
paste(round(myMean_fitted$speed,2), "m/s"), cex = 2)
</code></pre>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEji8rI8n0iXko25r6juWX3Aeysh-sS_hMmIHMygvguE2SV-nNZlQU99AJwtcnle6LX0fIt6ITS95wQkxS47H3s4FEwlOIKhr-hRvArUp0s3xhyphenhyphenobqf5AZNsRkvRrLtMT1q8beuIbjeGYMYR/s1600/Rplot04.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="831" data-original-width="796" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEji8rI8n0iXko25r6juWX3Aeysh-sS_hMmIHMygvguE2SV-nNZlQU99AJwtcnle6LX0fIt6ITS95wQkxS47H3s4FEwlOIKhr-hRvArUp0s3xhyphenhyphenobqf5AZNsRkvRrLtMT1q8beuIbjeGYMYR/s400/Rplot04.png" width="382" /></a></div>
<br />Javier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.com2tag:blogger.com,1999:blog-9045240617730855574.post-57598369853900613422017-05-09T15:44:00.000-07:002017-05-10T05:33:25.616-07:00niceOverPlot, or when the number of dimensions does matter<style type="text/css">p { margin-bottom: 0.25cm; direction: ltr; color: rgb(0, 0, 10); line-height: 120%; text-align: left; }p.western { font-family: "Liberation Serif",serif; font-size: 12pt; }p.cjk { font-family: "Noto Sans CJK SC Regular"; font-size: 12pt; }p.ctl { font-family: "FreeSans"; font-size: 12pt; }a:link { }</style>
<br />
<div align="justify" class="western" style="line-height: 100%; margin-bottom: 0cm; orphans: 0; page-break-after: avoid; widows: 0;">
<span style="font-family: inherit;"><span style="font-size: small;"><span style="font-family: inherit;">Hi there!<br /><br /> Over the last few months, my lab-mate Irene Villa <span style="font-family: inherit;">(see more of her work <a href="https://www.researchgate.net/profile/Irene_Villa-Machio" target="_blank">here!</a>) </span></span> 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 </span></span><span style="font-size: small;"><span style="font-family: inherit;">(1957)</span></span><span style="font-size: small;"><span style="font-family: inherit;"> 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 <a href="https://www.researchgate.net/publication/5990885_Grinnellian_and_Eltonian_niches_and_geographic_distributions_of_species/figures" target="_blank">BAM diagram</a> (Biotic-Abiotic-Movement; Soberón, 2007).<br /><br /> 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).<br /><br /> 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.</span></span></div>
<span style="font-family: inherit;"><span style="font-size: small;"><br /></span></span>
<span style="font-family: inherit;"><span style="font-size: small;"><br /></span></span>
<span style="font-family: inherit;"><span style="font-size: small;"> Here is the <a href="https://drive.google.com/open?id=0BzwZKMfGd8jvVXRvTUUtbU1HNEE" target="_blank">function code</a> (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 <a href="https://drive.google.com/open?id=0BzwZKMfGd8jvUF9sMkg5bGZpRWc" target="_blank">here</a>.</span></span><br />
<div align="justify" class="western" style="line-height: 100%; margin-bottom: 0cm; orphans: 0; page-break-after: avoid; widows: 0;">
<span style="font-family: inherit;"><span style="font-size: small;"><span style="font-family: inherit;"><br /></span></span></span>
<span style="font-family: inherit;"><span style="font-size: small;"><span style="font-family: inherit;"><br /></span></span></span></div>
<div align="justify" class="western" style="line-height: 100%; margin-bottom: 0cm; orphans: 0; page-break-after: avoid; widows: 0;">
<pre style="background: rgb(240, 240, 240) none repeat scroll 0% 0%; border: 1px dashed rgb(204, 204, 204); color: black; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><span style="font-family: inherit;"><span style="font-size: small;"><span style="font-family: inherit;"><code style="color: black; overflow-wrap: normal;">
##########################
####### EXAMPLE ######
##########################
library(ecospat)
library(ggplot2)
library(grid)
library(gridExtra)
library(gtable)
library(RColorBrewer)
# Read wall lizards ecological data
data_Ph_Pl<-read.csv("./data_Ph_Pl.csv")
# 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 (http://worldclim.org/).
# 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: GBIF.org (30th April 2017) GBIF Occurrence Download http://doi.org/10.15468/dl.yx4fbg
# Podarcis lilfordi: GBIF.org (30th April 2017) GBIF Occurrence Download http://doi.org/10.15468/dl.qqmzq3
# 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))
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)
</code></span></span></span></pre>
</div>
<div style="text-align: justify;">
<span style="font-family: inherit;"><span style="font-size: small;"><span style="font-family: inherit;"><br /></span></span></span>
</div>
<div style="text-align: justify;">
<span style="font-family: inherit;"><span style="font-size: small;"><span style="font-family: inherit;">Results from niceOverPlot (Blue area: <i>P. lilifordi</i>; Pink area: <i>P. hispanicus</i>)</span></span></span></div>
<div style="text-align: justify;">
<span style="font-family: inherit;"><span style="font-size: small;"><span style="font-family: inherit;"><br /></span></span></span>
</div>
<div class="separator" style="clear: both; text-align: justify;">
<span style="font-family: inherit;"><span style="font-size: small;"><span style="font-family: inherit;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEixBGXc4lSpYWgz6Y7uYsKoh86ALhgPKQ2GzKggQWH61z6kEKY8biOn3q6c8OEXeH32RuTSQZGfePTDfpssf7oSpHFhyphenhyphenYJCPAJFYMBH2QBkTIeS9jenSQSQPwuZ3PeYA_dzCuPr0J2E2oOC/s1600/niceOverlap_podarcis.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="640" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEixBGXc4lSpYWgz6Y7uYsKoh86ALhgPKQ2GzKggQWH61z6kEKY8biOn3q6c8OEXeH32RuTSQZGfePTDfpssf7oSpHFhyphenhyphenYJCPAJFYMBH2QBkTIeS9jenSQSQPwuZ3PeYA_dzCuPr0J2E2oOC/s640/niceOverlap_podarcis.png" width="640" /></a></span></span></span></div>
<div style="text-align: justify;">
<span style="font-family: inherit;"><span style="font-size: small;"><span style="font-family: inherit;"> </span></span></span></div>
<div style="text-align: justify;">
<span style="font-family: inherit;"><span style="font-size: small;"><span style="font-family: inherit;"><br /></span></span><span style="font-size: small;">
<span style="font-size: x-small;"><span style="font-family: inherit;"> <span style="font-size: small;">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.</span></span></span></span></span></div>
<div style="text-align: justify;">
In upcoming posts, we will discuss some alternative ways to address the D overlap index to take into account cases like our example.</div>
<div style="text-align: justify;">
<br />
<span style="font-family: inherit;">Have fun!</span></div>
<div style="text-align: justify;">
<span style="font-family: inherit;"><span style="font-size: small;"><br /></span></span></div>
<div class="separator" style="clear: both; text-align: center;">
<span style="font-family: inherit;"><span style="font-size: small;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjx6yhHW-_Xdn0t6Yv-gbO8u1kdcfPYC0kMfhasz_8pDNvJPgKEHhZCn43FFTUwIF3F-cyBIqul5HwvuHh9vFohzrjkuOfrTOWy1gtCCbyNDVyaxWFR5w7oSREo9zV1mN0_Ku95lzoajq9s/s1600/Pl.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="220" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjx6yhHW-_Xdn0t6Yv-gbO8u1kdcfPYC0kMfhasz_8pDNvJPgKEHhZCn43FFTUwIF3F-cyBIqul5HwvuHh9vFohzrjkuOfrTOWy1gtCCbyNDVyaxWFR5w7oSREo9zV1mN0_Ku95lzoajq9s/s400/Pl.png" width="400" /></a></span></span></div>
<div style="text-align: center;">
<i>Podarcis lilifordi,</i> from the Balearic Islands<br />
<br />
<div style="text-align: left;">
<b>References</b><br />
<br />
<br />
-Soberón, J. (2007). Grinnellian and Eltonian niches and geographic distributions of species. Ecology letters, 10(12), 1115-1123.<br />
<br />
-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.<br />
<br />
-Peterson, A. T. (2011). Ecological niche conservatism: A time‐structured review of evidence. Journal of Biogeography, 38(5), 817-827.<br />
<br />
-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.</div>
</div>
Javier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.com14tag:blogger.com,1999:blog-9045240617730855574.post-56433136695027142232016-12-04T10:40:00.000-08:002016-12-04T13:41:46.296-08:00rWind R package released!<div style="text-align: justify;">
Hi there! </div>
<div style="text-align: justify;">
<br /></div>
<div style="text-align: justify;">
Let me introduce you <a href="https://cran.r-project.org/web/packages/rWind/" target="_blank"><b>rWind</b></a>, an <b>R package</b> with several tools for <b>downloading</b>, <b>editing</b> and <b>converting</b> <b>wind data</b> from Global Forecast System (<a href="https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/global-forcast-system-gfs">https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/global-forcast-system-gfs</a>) in other formats as <b>raster for GIS</b>! 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 href="http://coastwatch.pfeg.noaa.gov/erddap/index.html">http://coastwatch.pfeg.noaa.gov/erddap/index.html</a>), a web service to dive into a lot of weather and oceanographic data-bases and download it easily. </div>
<div style="text-align: justify;">
<br /></div>
<div style="text-align: justify;">
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 (<a href="http://oos.soest.hawaii.edu/erddap/info/NCEP_Global_Best/index.html">http://oos.soest.hawaii.edu/erddap/info/NCEP_Global_Best/index.html</a>). At first, I was following this wonderful post from Conor Delaney (<a href="http://www.digital-geography.com/cloud-gis-getting-weather-data/#.WERgamd1DCL">http://www.digital-geography.com/cloud-gis-getting-weather-data/#.WERgamd1DCL</a>) 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 <a href="https://drive.google.com/open?id=0BzwZKMfGd8jvandyUG1aMjdxR0U" target="_blank">reference manual</a> and an R code with a brief <a href="https://drive.google.com/open?id=0BzwZKMfGd8jvaUdRSGtOajMyYW8" target="_blank">tutorial</a> to get familiar with the utilities of the rWind package!<br />
<br />
If you have any doubt or you want to report a bug or make any suggestion, please, comment the post or write me: jflopez@rjb.csic.es or jflopez.bio@gmail.com. </div>
<div style="text-align: justify;">
<br /></div>
<div style="text-align: justify;">
Enjoy it!<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjdQ8tMP0_rEMghc8-tpcbSuCOlE61KQJUJEgvzr24NaeX9rFMHalKY5VkU4aHntL-6kHu_QWkWB0yGHuVc0txspbYgR5SKUawPsBVYw26D_JJF2tgTxj4jWNIVTJqxxlXFf04ptDMKO05j/s1600/blog4.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="287" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjdQ8tMP0_rEMghc8-tpcbSuCOlE61KQJUJEgvzr24NaeX9rFMHalKY5VkU4aHntL-6kHu_QWkWB0yGHuVc0txspbYgR5SKUawPsBVYw26D_JJF2tgTxj4jWNIVTJqxxlXFf04ptDMKO05j/s320/blog4.png" width="320" /></a></div>
<br />
<br />
<u>Javier Fernández-López</u> (2016). rWind: Download, Edit and Transform Wind Data from GFS. R package version 0.1.3. https://CRAN.R-project.org/package=rWind</div>
<div style="text-align: justify;">
<br /></div>
<div style="text-align: justify;">
<br /></div>
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">
# Download and install "rWind" package from CRAN:
install.packages("rWind")
# You should install also "raster" package if you do not have it
library(rWind)
library(raster)
packageDescription("rWind")
help(package="rWind")
# "rWind" is a package with several tools for downloading, editing and transforming wind data from Global Forecast
# System (GFS, see <https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/global-forcast-system-gfs>) of the USA's
# National Weather Service (NWS, see <http://www.weather.gov/>).
citation("rWind")
# > Javier Fernández-López (2016). rWind: Download, Edit and Transform Wind Data from GFS. R package version 0.1.3.
# > https://CRAN.R-project.org/package=rWind
# 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)
wind.dl(2015,2,12,0,-10,5,35,45)
# By default, this function generates an R object with downloaded data. You can store it...
wind_data<-wind.dl(2015,2,12,0,-10,5,35,45)
head(wind_data)
</code></pre>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgP_Q0Ubrs1c_75g7FVvIFxKApjgU2HdpjH2vrNlrZPwoLHXa_DVtGr51DYqVTdd9P5t4GQZf5225zJmEaIewzPiBSovp4_kqV5UYMK6b70Rc_rBSfZMqsGRbW1R-gTeDqYc1eMQKjhWdc1/s1600/blog_rWInd1.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="78" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgP_Q0Ubrs1c_75g7FVvIFxKApjgU2HdpjH2vrNlrZPwoLHXa_DVtGr51DYqVTdd9P5t4GQZf5225zJmEaIewzPiBSovp4_kqV5UYMK6b70Rc_rBSfZMqsGRbW1R-gTeDqYc1eMQKjhWdc1/s400/blog_rWInd1.png" width="400" /></a></div>
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">
# or download </code><code style="color: black; word-wrap: normal;"><code style="color: black; word-wrap: normal;">a CVS file </code>into your work directory with the data using type="csv" argument:
getwd()
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 "wind.fit" function:
#help(wind.fit)
wind_data<-wind.fit(wind_data)
head(wind_data)
</code></pre>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg1jvc3vyDnvfpi_3qIhcoPffwwEfBgIEtX3G4RhSaj4RcMX_GNUKVzRwoOf4yeqfawu1rddNYQMpGABEkLuplbOARgSbc3Yt71bLkhL-oTiObSk8qPUoVDnAP6hQGRoc6vEL7E2RI1jEpH/s1600/Captura+de+pantalla+de+2016-12-04+20-08-08.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="171" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg1jvc3vyDnvfpi_3qIhcoPffwwEfBgIEtX3G4RhSaj4RcMX_GNUKVzRwoOf4yeqfawu1rddNYQMpGABEkLuplbOARgSbc3Yt71bLkhL-oTiObSk8qPUoVDnAP6hQGRoc6vEL7E2RI1jEpH/s320/Captura+de+pantalla+de+2016-12-04+20-08-08.png" width="320" /></a></div>
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> # 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!
#install.packages("rworldmap")
library(rworldmap)
newmap <- getMap(resolution = "low")
par(mfrow=c(1,2))
plot(r_dir, main="direction")
lines(newmap, lwd=4)
plot(r_speed, main="speed")
lines(newmap, lwd=4)
</code></pre>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiZBHcWrVbF36uL7vThxc9NMFkf-yH6UrhoVSMVJHXIAto4BnpXlwIsXN-keGdralvEUUSyjCRTpcvPdJTkTZTxDQlp12weLvmVupIiRpT4Go8rla18u19KH5X6g_-FZBe351KHXgaGDkOj/s1600/rWind_1.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="150" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiZBHcWrVbF36uL7vThxc9NMFkf-yH6UrhoVSMVJHXIAto4BnpXlwIsXN-keGdralvEUUSyjCRTpcvPdJTkTZTxDQlp12weLvmVupIiRpT4Go8rla18u19KH5X6g_-FZBe351KHXgaGDkOj/s400/rWind_1.png" width="400" /></a></div>
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">
# Additionally, you can use arrowDir and Arrowhead (from "shape" package) functions to plot wind direction
# over a raster graph:
#install.packages("shape")
library(shape)
dev.off()
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")
</code></pre>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg7mUm-BkONXky34YwHOvqJGw0RpMnT00fsqbN8jAVSOsUUqVgjkCSpcHOfxW1QBwxgVggwQbvfbYchN1NKld5BHtb_gpJHbOOmOUJa-NyWNCyCMJhdCsX-e5fsw1vksEOyyvI5r0fSUWJu/s1600/blog4.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="287" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg7mUm-BkONXky34YwHOvqJGw0RpMnT00fsqbN8jAVSOsUUqVgjkCSpcHOfxW1QBwxgVggwQbvfbYchN1NKld5BHtb_gpJHbOOmOUJa-NyWNCyCMJhdCsX-e5fsw1vksEOyyvI5r0fSUWJu/s320/blog4.png" width="320" /></a></div>
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">
# 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){
w<-wind.dl(2015,2,d,0,-10,30,35,70)
wind_serie[[d]]<-w
}
wind_serie
# Finally, you can use wind.mean function to calculate wind average
wind_average<-wind.mean(wind_serie)
wind_average<-wind.fit(wind_average)
r_average_dir<-wind2raster(wind_average, type="dir")
r_average_speed<-wind2raster(wind_average, type="speed")
par(mfrow=c(1,2))
plot(r_average_dir, main="direction average")
lines(newmap, lwd=1)
plot(r_average_speed, main="speed average")
lines(newmap, lwd=1)
</code></pre>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiuPFKsTfHQ54R-8UBvLUyEMffYnnHq3dlGG9rKb4oQ-J-QlFFZgWxZLS0T487p-nk7wUpDZ4ayXO7qC0h2olyM5y7WzU_z7O-A39DT5P3OKIaZ_VPr3N_cu9HdxUJq0F1nALucbHK8wFWi/s1600/rwind_plot2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="281" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiuPFKsTfHQ54R-8UBvLUyEMffYnnHq3dlGG9rKb4oQ-J-QlFFZgWxZLS0T487p-nk7wUpDZ4ayXO7qC0h2olyM5y7WzU_z7O-A39DT5P3OKIaZ_VPr3N_cu9HdxUJq0F1nALucbHK8wFWi/s400/rwind_plot2.png" width="400" /></a></div>
Javier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.com20tag:blogger.com,1999:blog-9045240617730855574.post-82458545381860667972016-11-02T07:24:00.001-07:002016-11-02T09:34:34.453-07:00eWalk, an R function to simulate two-dimensional evolutionary walks<style type="text/css">p { margin-bottom: 0.25cm; line-height: 120%; }</style>
<br />
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<span style="font-size: small;">Hey
there! </span></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<br /></div>
<div style="text-align: justify;">
<span style="font-size: small;">
</span></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<span style="font-size: small;"> During
last months, I was thinking with my friend Javi Fajardo (see his work
<a href="https://www.researchgate.net/profile/Javier_Fajardo4" target="_blank">here</a>!) 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.</span></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh2OQdFya04KeTmPD3DHafYracIk2iIw92L_4XvZdLCA14u_LR5B3NvviwaImkUzd-w7RhwxLArLaNsvGyI1qcIfKppHuWUCXlxR-ViC4Bd0RfEonh-pC2mwqe4rZRf_ytzoYAB9p8JvKLJ/s1600/Sin+t%25C3%25ADtulo_2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="164" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh2OQdFya04KeTmPD3DHafYracIk2iIw92L_4XvZdLCA14u_LR5B3NvviwaImkUzd-w7RhwxLArLaNsvGyI1qcIfKppHuWUCXlxR-ViC4Bd0RfEonh-pC2mwqe4rZRf_ytzoYAB9p8JvKLJ/s400/Sin+t%25C3%25ADtulo_2.png" width="400" /></a></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 120%; }</style>
</div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<span style="font-family: inherit;"><span style="font-style: normal;"> <span style="font-size: small;">M</span></span><span style="font-size: small;"><span style="font-style: normal;">y</span></span><span style="font-size: small;"><i>
</i></span><span style="font-size: small;"><span style="font-style: normal;">colleague</span></span><span style="font-size: small;"><span style="font-style: normal;">
Javi and I were </span></span><span style="font-size: small;"><span style="font-style: normal;">talking</span></span><span style="font-size: small;"><span style="font-style: normal;">
about the applications of this models in two dimensions, representing
for example space coordinates (latitude/longitude). </span></span><span style="font-size: small;"><span style="font-style: normal;">They
have been used, for example, to simulate predator/prey movements,
home range movements, etc. </span></span><span style="font-size: small;"><span style="font-style: normal;">We
developed an R function similar to the previous </span></span><span style="font-size: small;"><i>eMotion</i></span><span style="font-size: small;"><span style="font-style: normal;">
function in order implement this kind of movement models: </span></span><span style="font-size: small;"><i><b>eWalk</b>
(</i></span><span style="font-size: small;"><span style="font-style: normal;">from
</span></span><span style="font-size: small;"><i>evolutiveWalk...</i></span><span style="font-size: small;"><span style="font-style: normal;">). </span></span><span style="font-size: small;"><i>eWalk</i></span><span style="font-size: small;"><span style="font-style: normal;">
is similar to other movement-simulation functions in </span></span><span style="font-size: small;"><i>adehabitatLT
</i></span><span style="font-size: small;"><span style="font-style: normal;">package
(Calenge, 20</span></span><span style="font-size: small;"><span style="font-style: normal;">0</span></span><span style="font-size: small;"><span style="font-style: normal;">6).
However, you can use </span></span><span style="font-size: small;"><i>eWalk
</i></span><span style="font-size: small;"><span style="font-style: normal;">to
simulate a variety of processes as Brownian motion, OU processes or
Lèvy flights (see below).</span></span></span></div>
<div style="text-align: justify;">
<span style="font-size: small;">
</span></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<span style="font-size: small;"><br /></span>
</div>
<div style="text-align: justify;">
<span style="font-size: small;">
</span></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<span style="font-family: inherit;"><span style="font-size: small;"><span style="font-style: normal;">Here
you have the </span></span><span style="font-size: small;"><span style="font-style: normal;">(</span></span><span style="font-size: small;"><span style="font-style: normal;">commented</span></span><span style="font-size: small;"><span style="font-style: normal;">!)
</span></span><span style="font-size: small;"><span style="font-style: normal;">code</span></span><span style="font-size: small;"><span style="font-style: normal;">
</span></span><span style="font-size: small;"><span style="font-style: normal;">to
load the function into R console:</span></span></span></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<br /></div>
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">
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!
}
></-></-></code></pre>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<br /></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<span style="font-family: inherit;"><span style="font-size: small;"><span style="font-style: normal;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 120%; }</style>
</span></span></span></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
<span style="font-size: small;"><span style="font-family: inherit;"><span style="font-style: normal;"> The
needed parameters to run the function are similar to those for
</span><i>eMotion</i><span style="font-style: normal;">,
but taking into account that now you have two variables
(latitude/longitude, for example) at the same time.</span></span></span></div>
<span style="font-size: small;"><span style="font-family: inherit;">
</span></span>
<br />
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
<span style="font-size: small;"><br /></span></div>
<span style="font-size: small;">
</span>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhxlaR7sipQpSsWHSwGB_uO3Mr3hv-mXvqqHb7kwpNyDBzVFEY-EvC7bo3z58q88ApVwsT1VI6HDbn5BHzX5jl6zuc-PsdpdQG458obo269pdbTpqWoD-UyMpVgpXmmyj96HiS9sY60Zotl/s1600/Rplot.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhxlaR7sipQpSsWHSwGB_uO3Mr3hv-mXvqqHb7kwpNyDBzVFEY-EvC7bo3z58q88ApVwsT1VI6HDbn5BHzX5jl6zuc-PsdpdQG458obo269pdbTpqWoD-UyMpVgpXmmyj96HiS9sY60Zotl/s400/Rplot.png" width="380" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<span style="font-size: small;"></span></div>
<span style="font-size: small;">
</span>
<br />
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 120%; }</style><span style="font-size: small;">
</span></div>
<span style="font-size: small;">
</span>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
<span style="font-family: inherit;"><span style="font-style: normal;">eWalk(sigma,initial_lon,initial_lat,number
of generations, lon_alpha, lon_optimum, y_alpha, y_optimum, color,
levy=FALSE, plot=TRUE)</span></span></div>
<div style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<span style="font-family: inherit;"><span style="font-size: small;"><span style="font-style: normal;"></span></span></span></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm;">
<span style="font-size: small;">
</span></div>
<span style="font-family: inherit;"><span style="font-size: small;"><span style="font-style: normal;"></span></span></span><span style="font-family: inherit;"><span style="font-size: small;"><span style="font-style: normal;"></span></span></span><span style="font-family: inherit;"><span lang="en-US"><span style="font-style: normal;"></span></span></span><br />
<div style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgCpR8i8l4yUlkkC0P0N57QSvDmzPC0kAN_QlMGiFucn55-5pkK1arrLmiB95L9ZrO7ct4kuGIb-zuk8pkrFSvDOsnezK2RzBBPZWERoO_04weHGYrTanbh3ZfbhxPa_MmhAI6WDFaZe4xy/s1600/cauchy.png" imageanchor="1" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" height="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgCpR8i8l4yUlkkC0P0N57QSvDmzPC0kAN_QlMGiFucn55-5pkK1arrLmiB95L9ZrO7ct4kuGIb-zuk8pkrFSvDOsnezK2RzBBPZWERoO_04weHGYrTanbh3ZfbhxPa_MmhAI6WDFaZe4xy/s320/cauchy.png" width="304" /></a></div>
<span style="font-family: inherit;"><span lang="en-US"><span style="font-style: normal;">I</span></span><span lang="en-US"><span style="font-style: normal;">n
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
</span></span><span lang="en-US"><span style="font-style: normal;">searching</span></span><span lang="en-US"><span style="font-style: normal;">
movements of predators </span></span><span lang="en-US"><span style="font-style: normal;">(Bartumeus
et al., 2005)</span></span><span lang="en-US"><span style="font-style: normal;">,
or have been used to simulate</span></span><span lang="en-US"><span style="font-style: normal;"><span style="font-weight: normal;">
punctuated equilibrium in trait</span></span></span><span lang="en-US"><span style="font-style: normal;"><span style="font-weight: normal;">s</span></span></span><span lang="en-US"><span style="font-style: normal;"><span style="font-weight: normal;">
evolution (Gould </span></span></span><span lang="en-US"><span style="font-style: normal;"><span style="font-weight: normal;">and
Eldredge, 1977; </span></span></span><span lang="en-US"><span style="font-style: normal;"><span style="font-weight: normal;">Landis
et al., 201</span></span></span><span lang="en-US"><span style="font-style: normal;"><span style="font-weight: normal;">3</span></span></span><span lang="en-US"><span style="font-style: normal;"><span style="font-weight: normal;">).
</span></span></span><span lang="en-US"><span style="font-style: normal;"><span style="font-weight: normal;">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). </span></span></span><span lang="en-US"><span style="font-style: normal;"><span style="font-weight: normal;">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.</span></span></span></span></div>
<div style="text-align: justify;">
<span style="font-family: inherit;">
</span></div>
<div style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<span style="font-family: inherit;"><br /></span>
</div>
<div style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<span style="font-family: inherit;"><br /></span>
</div>
<div style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<span style="font-family: inherit;"><br /></span>
</div>
<div style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<span style="font-family: inherit;"><br /></span>
</div>
<div style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<span style="font-family: inherit;"><br /></span>
</div>
<div style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<span style="font-family: inherit;"><br /></span>
</div>
<div style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<span style="font-family: inherit;"><br /></span>
</div>
<div style="text-align: justify;">
<span style="font-family: inherit;">
</span></div>
<div style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<span style="font-family: inherit;"><span lang="en-US"><span style="font-style: normal;"><span style="font-weight: normal;">Here
you have some examples using eWalk and their code to reproduce them:</span></span></span></span></div>
<div style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">
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")
</code></pre>
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjkQgFtLBuxCjJYJIG6JUZgx-2gQaIAE4OLpoZ60PiHoxXCRYiRjKN1iiAcDnhhdn6CZEN8mhQUJSUWl_8Jzw6jkz57zDRAObPwyj4rr5GlxVVtYrzW9-xIaZvKUgBohJW0e6OID7t5dHQ1/s1600/4plots2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjkQgFtLBuxCjJYJIG6JUZgx-2gQaIAE4OLpoZ60PiHoxXCRYiRjKN1iiAcDnhhdn6CZEN8mhQUJSUWl_8Jzw6jkz57zDRAObPwyj4rr5GlxVVtYrzW9-xIaZvKUgBohJW0e6OID7t5dHQ1/s400/4plots2.png" width="380" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<br />
<div style="text-align: justify;">
<span style="font-family: inherit;">
</span></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<b><span style="font-family: inherit;"><span style="font-style: normal;"> References</span></span></b></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<span style="font-family: inherit;"><span style="font-style: normal;"><br /></span></span></div>
<div lang="en-US" style="line-height: 100%; margin-bottom: 0cm; text-align: justify;">
<br />
<style type="text/css">p { margin-bottom: 0.25cm; line-height: 120%; }</style>
</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
-Bartumeus, F., da
Luz, M. E., Viswanathan, G. M., & Catalan, J. (2005). Animal
search strategies: a quantitative random‐walk analysis. <i>Ecology</i>,
<i>86</i>(11), 3078-3087.</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
-Calenge, C. (2006).
The package “adehabitat” for the R software: a tool for the
analysis of space and habitat use by animals. <i>Ecological
modelling</i>, <i>197</i>(3), 516-519.</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
-Gould, S. J., &
Eldredge, N. (1977). Punctuated equilibria: the tempo and mode of
evolution reconsidered. <i>Paleobiology</i>, <i>3</i>(02), 115-151.</div>
<div style="line-height: 100%; margin-bottom: 0cm;">
<br /></div>
<div style="line-height: 100%; margin-bottom: 0cm;">
-Landis, M. J.,
Schraiber, J. G., & Liang, M. (2013). Phylogenetic analysis using
Lévy processes: finding jumps in the evolution of continuous traits.
<i>Systematic biology</i>, <i>62</i>(2), 193-204.</div>
Javier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.com2tag:blogger.com,1999:blog-9045240617730855574.post-71992522976269314432016-07-10T17:52:00.000-07:002016-07-15T15:53:47.260-07:00eMotion, an R function for evolutionary motion of character valuesHi all!<br />
<br />
Some time ago, I was studying an article from Tamara Münkemüller et <i>al</i>. 2015: "<b>Phylogenetic niche conservatism – common pitfalls and ways forward</b>". 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.<br />
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!<br />
<br />
However, I would like to focus in one of the main pictures of the article:<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgNKBAYQGvZez6dCseCkhE1vNymqTDsHJE47O8L2BTsjtJv_2PYyF-v54lIZtQ5GkT69N5on-nr4w2Mh6hWiQffo3zPq4JMRhisjFmvrDTo2KhLG7SeRgH6Aw5AdmtiU384243RYOUSrpaY/s1600/tamara.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="385" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgNKBAYQGvZez6dCseCkhE1vNymqTDsHJE47O8L2BTsjtJv_2PYyF-v54lIZtQ5GkT69N5on-nr4w2Mh6hWiQffo3zPq4JMRhisjFmvrDTo2KhLG7SeRgH6Aw5AdmtiU384243RYOUSrpaY/s400/tamara.png" width="400" /></a></div>
<br />
<br />
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!<br />
<br />
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.<br />
<br />
The graphic is the result of applying the OU process equation:<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi5hfwhVq96jBcsrNR23xrmzdAHLmK0tl05aelzv4Puc7sEZ26_uJgfp_pGLGSA4YJ_wMRRIrtj_0zI-dSO4cWTeyFCnYsyFzZWOre-iumw_fpHzQ4iO_tw7bdK1o6KKjA9Pa9-3jzAc9Qo/s1600/tamara_formula.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="35" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi5hfwhVq96jBcsrNR23xrmzdAHLmK0tl05aelzv4Puc7sEZ26_uJgfp_pGLGSA4YJ_wMRRIrtj_0zI-dSO4cWTeyFCnYsyFzZWOre-iumw_fpHzQ4iO_tw7bdK1o6KKjA9Pa9-3jzAc9Qo/s400/tamara_formula.png" width="400" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<span style="background-color: white; font-family: "open sans"; font-size: 13.2px; line-height: 18.48px;">x(t+dt)= x(t)+ alpha(theta(t)-x(t))dt + (sigma dt E)</span></div>
<br />
where , <span style="background-color: white; font-family: "open sans"; line-height: 18.48px;">x(t) is the trait value at time=t, </span><span style="background-color: white; font-family: "open sans"; line-height: 18.48px;">alpha is the selection strength with which a specific trait value (optimum) is wanted, </span><span style="background-color: white; font-family: "open sans"; line-height: 18.48px;">theta is the optimal value of the trait in time=t, </span><span style="background-color: white; font-family: "open sans"; line-height: 18.48px;">sigma is Brownian motion rate (the amount of random trait variability per time) and </span>E is a normal distribution(0,1).<br />
<br />
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")<br />
<br />
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). <br />
<br />
Here you have the code to load the function into R console:<br />
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> 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=""))
trait_a=numeric()
simulations= list(1:n)
for (j in 1:n){
a=0
for (i in 1:generations) {
a= a + (sigma * rnorm(1)) + (alpha * ( theta - a))
trait_a[i]=a
}
aa=data.frame(trait_a)
simulations[j]=aa
}
for(k in 1:n){
bb=data.frame(simulations[k])
lines(bb, col="red")
}
#######
trait_a=numeric()
simulations= list(1:n)
for (j in 1:n){
a=0
for (i in 1:generations) {
a= a + (sigma2 * rnorm(1)) + (alpha2 * ( theta2 - a))
trait_a[i]=a
}
aa=data.frame(trait_a)
simulations[j]=aa
}
for(k in 1:n){
bb=data.frame(simulations[k])
lines(bb, col="blue")
}
}
</code></pre>
<br />
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.<br />
<br />
eMotion(sigma, theta, alpha, sigma2, theta2, alpha2, generations, n)<br />
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> eMotion(2, 100, 0.002, 0.5,-100, 0.002, 10000, 20) </code></pre>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjRX_8KAnC53a_Yt6TUdCgxb9DiYQXC-KqySQXrWRuSeFJn95ZQ5bV3sKLvgxi7i4Eg3Er2II6JIYmWXbGQi9-lwUWZRPUWOdiRltcU9UhIov_kZM-4d2Oxt9FFlscGK0h5EVXhgdpy7SaO/s1600/eMotion.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjRX_8KAnC53a_Yt6TUdCgxb9DiYQXC-KqySQXrWRuSeFJn95ZQ5bV3sKLvgxi7i4Eg3Er2II6JIYmWXbGQi9-lwUWZRPUWOdiRltcU9UhIov_kZM-4d2Oxt9FFlscGK0h5EVXhgdpy7SaO/s400/eMotion.png" width="400" /></a></div>
<br />
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).<br />
<br />
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!<br />
<br />
<pre style="background: #f0f0f0; border: 1px dashed #cccccc; color: black; font-family: "arial"; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> eMotion(1, 0, 0, 1,0, 0.01, 10000, 20) </code></pre>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgOrqDO4dEEdcE5j88y2ULPhliGd3hVKryhlb4eAmgP6A0qQ-zjCA8ed5Crl5IjeXj9235SExaRARsFCe874VxWW2Q5x7dsCXq8iVQFB3Pb46z_J4-nAl6000Wl8c_MXGDsQzEnq7J7ZAbG/s1600/eMotion.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgOrqDO4dEEdcE5j88y2ULPhliGd3hVKryhlb4eAmgP6A0qQ-zjCA8ed5Crl5IjeXj9235SExaRARsFCe874VxWW2Q5x7dsCXq8iVQFB3Pb46z_J4-nAl6000Wl8c_MXGDsQzEnq7J7ZAbG/s400/eMotion.png" width="400" /></a></div>
Javier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.com0tag:blogger.com,1999:blog-9045240617730855574.post-32376906625630116252016-06-27T03:02:00.002-07:002016-07-15T16:03:42.244-07:00nichEvolve documentation and code<br />
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!<br />
Link: <a href="https://drive.google.com/open?id=0BzwZKMfGd8jvdVF3N2Z3NHJVbmc" target="_blank">nichEvolve.nlogo</a><br />
<br />
*(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: <a href="https://ccl.northwestern.edu/netlogo/" target="_blank">NetLogo</a><br />
<br />
<b style="background-color: #b4a7d6;">WHAT IS IT?</b><br />
<br />
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.<br />
The model could be usefull as a virtual laboratory to teach about this concepts.<br />
<br />
author: Javier Fernandez-Lopez <span style="background-color: white; font-family: "arial" , sans-serif; font-size: 13px; white-space: nowrap;"><span style="color: blue;">jflopez.bio@gmail.com</span></span><br />
<br />
<b style="background-color: #b4a7d6;">HOW IT WORKS</b><br />
<br />
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<br />
<br />
x(t+dt)= x(t)+ alpha(theta(t)-x(t))dt + (sigma dt E)<br />
<br />
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.<br />
<br />
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.<br />
<br />
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.<br />
<br />
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.<br />
<br />
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".<br />
<br />
<span style="background-color: #b4a7d6;"><b>HOW TO USE IT</b></span><br />
<br />
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.<br />
<br />
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.<br />
<br />
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.<br />
<br />
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)<br />
<br />
4) Finally you could check the behaivour of mobile agents when climate-change is "on"<br />
<br />
<br />
<span style="background-color: #b4a7d6;"><b>THINGS TO NOTICE</b></span><br />
<br />
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. <br />
<br />
<br />
<b style="background-color: #b4a7d6;">HOW TO CITE</b><br />
<br />
Fernandez-Lopez, J. 2016. nichEvolve: a niche evolution model<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhhF24FXTQPENq3nKX2PzJp9EunbPnhS_Q2qAC2FNI2DeZbkeFXeE9n5I4_RKHMqneQUcZSepB1_RsIrSkqDEY1lEAeN-NqTgh7-F_rC5EkgjUZ5BRHmiMiWoF8RWx_2nR0TSNKItr9vbQq/s1600/nichEvolve+interface2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="263" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhhF24FXTQPENq3nKX2PzJp9EunbPnhS_Q2qAC2FNI2DeZbkeFXeE9n5I4_RKHMqneQUcZSepB1_RsIrSkqDEY1lEAeN-NqTgh7-F_rC5EkgjUZ5BRHmiMiWoF8RWx_2nR0TSNKItr9vbQq/s400/nichEvolve+interface2.png" width="400" /></a></div>
<br />
<div class="separator" style="clear: both; text-align: center;">
</div>
<br />Javier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.com1tag:blogger.com,1999:blog-9045240617730855574.post-28443098501335651852016-06-26T01:01:00.000-07:002016-07-15T16:05:20.322-07:00nichEvolve, a niche evolution model with NetLogo<br />
During the last year I was reading some articles about niche conservatism/evolution. There are a lot of theories and methodologies to deal with this kind of questions, but sometimes it's hard to understand all this mathematical and theoretical concepts about models of evolution: Brownian motion models, Ornstein-Uhlenbeck processes, etc. Following KISS philosophy (Keep It Simple, Stupid!) I wrote some code using the wonderful tool NetLogo to develop a very simple model that helps to understand some simple processes related with ecological traits evolution...<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh5yQWsZrmRw8zwIG45wnsFAnf8pfslVDo_SCTA1NzHtBabQ3cwh4QuxBblPBUrGkUPbejWfADBB-qqZcM6saCbUzIlpBWWTMb4nf3YCSiYt18Y6xNOu-WV3LSMLJP0tVnFiU2hBQnGO_xs/s1600/niche_graph.png" imageanchor="1" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" height="147" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh5yQWsZrmRw8zwIG45wnsFAnf8pfslVDo_SCTA1NzHtBabQ3cwh4QuxBblPBUrGkUPbejWfADBB-qqZcM6saCbUzIlpBWWTMb4nf3YCSiYt18Y6xNOu-WV3LSMLJP0tVnFiU2hBQnGO_xs/s400/niche_graph.png" width="400" /></a></div>
The model implements two usual processes used to describe the evolution of trait values: Brownian motion and Ornstein-Uhlenbeck. While Brownian motion process represents random evolution of a trait value, depending only on Brownian rate (amount of change allowed per time), Ornstein-Uhlenbeck process means directional evolution towards an optimum with a specific selection strength. In the model you can explore diverse combinations of each parameter and observe the results also in a spatial explicit display with mobile agents that track environmental conditions trying to match their theoretical requirements.<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiYMbJECbJgS-hraWHX0RNSCTwPvO-rr1QH2D1i9afBwhZL67d0px18vMDWw0aNGJspZWe2QKjAkk6woSgQw3IU2hbatntWLuSlqXdxFu4SebxHKFjmDAWPOlvWGdOQdQKpMtxVrTiV57Bu/s1600/niche_evolve+interface.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="162" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiYMbJECbJgS-hraWHX0RNSCTwPvO-rr1QH2D1i9afBwhZL67d0px18vMDWw0aNGJspZWe2QKjAkk6woSgQw3IU2hbatntWLuSlqXdxFu4SebxHKFjmDAWPOlvWGdOQdQKpMtxVrTiV57Bu/s400/niche_evolve+interface.png" width="400" /></a></div>
<br />
Code will be uploaded soon, so you could modify it and report any suggestion or bug.<br />
<br />
Cheers!<br />
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<br />
<br />Javier Fernández-Lópezhttp://www.blogger.com/profile/04195699007682768268noreply@blogger.com2