R tips from around the web...

The great thing of R, is that the number of available resources on the web is increasing dramatically. If you cannot afford expensive books, or you are looking straight questions or some tutorials you will find almost anything you need. Here's my short list of best place to learn and explore R other than the CRAN website:

Searching, Exploring and Procrastinating

R Seek: You probably use google while searching for specific commands. But this search engine is much more powerful. It looks for R related websites, and you can also filter to "functions", "blogs", "books" and etc. Awesome!

R Graph Gallery : This is unfortunately not updated anymore, but you can find fancy graphs and codes which can inspire you for something even better.

R Graphical Manual: Visual learning is always the best. This is probably the principle of this website, which allows you to search for functions and packages from the standpoint of their graphical outputs.

R Reference Card : This is actually from CRAN, but it's an extremely handy reference card that you can print and put on your office wall to show people how geeky you are. Most of them are daily uses functions that you'll probably learn soon anyway, but it's extremely useful if you've just got started.

R Bloggers: New research method, elaborated codes, inspiring plots, novel ideas on everything related to R cannot be found in books. The cutting edge is in the blogosphere, where hundreds of peer R users are perhaps tackling the same issue you are working one. This websites put together almost 200 bloggers from around the world. Just put in you're RSS feed, and you'll get about 10 daily posts with new ideas. Inspiring.


Learning and Tutorials

Quick-R. Nice and intuitive website on R, covering wide range of topics from analysis to plotting

Quantitative Archaeology Wiki . This is still under development but has some sections on R specifically designed for archaeologists.

There are of course plenty of other tutorials around the web. If you want use GRASS and R together you can find a nice introductory tutorial on geostats here. A number of academics also share their courses online so you can sneak in and learn some great stuff like this one. My favourite course isRichard McElreath's  Statistical Thinking in Evolutionary Ecology Course. This is simply mind-opening.

Running an R-based ABM in parallel on a Multicore Desktop

I've been running for the last couple of months a lot of simulations written in R on legion cluster  here at UCL. But I still do many things on our quad-core server here at the institute, and parallelising the simulations is extremely handy.

So here's my recipe:

1)Firstly, the main ingredient,  you need an R function for your model where you have a set of parameters to sweep, something like:

Rsim(InitialPopulation=100,K=5000,ReproductiveRate=0.2,TransmissionRate=0.4,MutationRate=0.1)
{
#Put Something amazing and cool here
# ...


return(list(pop=pop,fitness=fitness,diversity=diversity))
}

The function can return everything, in my case this will be a list with three vectors, say $pop, $fitness and $diversity.

2)For the actual parallelisation, you need the following two packages: foreach and doMC and their dependencies. After loading both packages you specify the number of cores you want to use, as follows:

library(foreach)
library(doMC)

cores=4  
registerDoMC(cores)


3)Now, at this point I create a sweepSpace, which is basically a data.frame with all the combinations of parameters that I want to explore, repeated n times, with n being the number of runs. So in my case I want to sweep the MutationRate, the TransmissionRate and ReproductiveRate for 100 runs.


runs=50
mutationSweep=c(0.0001,0.001,0.01,0.1)
reproductionSweep=c(0.1,0.2,0.3)
transmissionSweep=c(0,0.5,1.0)


sweepSpace=expand.grid(runs=1:runs,m= mutationSweep,r= reproductionSweep,z= transmissionSweep)


4)Now we are ready to run the simulation with the following code


res=foreach(a=sweepSpace$runs,b=sweepSpace$m,c=sweepSpace$r,d=weepSpace$z)%dopar%{Rsim(InitialPopulation=100,K=5000,ReproductiveRate=c,TransmissionRate=d, TransmissionRate =b)}


The foreach() function will read each row of sweepSpace and assign the values of each column  to the letters, and then use these for the Rsim() function. The resulting object is a list with a length equal to the number of rows of sweepSpace (thus in this case 500), containing the results of Rsim() in each of them. The parallelisation allows  the simultaneous but independent run for all combinations and the storage of the results which follows the order in the sweepSpace. So if you have a specific parameter combination which takes longer to compute, the foreach() function allows to carry on, in the meanwhile, all the other combinations and put the results in the correct place.

5)Get and plot your results!

The sweepSpace can also be used as map to go through the resulting list.
For instance, if you want to create a function which allows the plotting of all the population dynamics on a single plot, you can write something as the following code:

plotSimPop=function(res,sweepSpace, ReproductiveRate, TransmissionRate, TransmissionRate)
{
#First we want to retrieve the index of list having the specific parameters that we've requested:


tmp=res[which(sweepSpace$r==ReproductiveRate&sweepSpace$m== MutationRate&sweepSpace$z==TransmissionRate)]


#we then retrieve the number of runs, and create a matrix with our results, with row corresponding to the number of runs, and columns to the number of timesteps (here set to 1000):


nruns=length(tmp)
resMat=matrix(0,nruns,1000)


#We then loop through the results collecting and storing the population dynamics in our matrix


for (x in 1:nruns)
{
resMat[x,]=res[[x]]$pop
}


#And we can now plot the results, with the average dynamics in red:
average=apply(resMat,2,mean,na.rm=TRUE)
plot(resMat[1,],type="n",xlim=c(1,1000),ylim=c(min(resMat),max(resMat),xlab="Time",ylab="Population")
for (x in 1:nruns)
{
lines(resMat[x,],col="grey")
}
lines(average,col="red")
}

apply() function and ABM in R

I know know...I've been away again...
We (myself and Mark Lake) are presenting a paper at the CECD conference and we have still some to stuff to finish...so I'm really, really busy... I'll post asap a much more detailed post on the conference and on our paper, but before that I just wanted to share an useful link I found this morning which would have been handy a couple of month ago.
As I said, I'm writing my ABM in R this time. There are many good reasons (but also bad reasons) for this which I'm gonna write on another post (yes I keep promising...). Having said that, R is terribly slow. Yes  you can write things in C and call them inside your function, but the main reason is that it is terribly bad on looping. And an ABM involves a lot of looping. Then I realised that many people avoid using  loops in R, and instead use the series of "apply" function. These are however hard to grasp, but this blog explains it very elegantly. I wish I had read this long ago...Anyway, when you master the "apply" family, you can also play around with the mcapply() function of the multicore package which parallels the apply function through your cores, fastening alot your simulation!

R Plotting tips

Yes...I haven't been updating the blog for a very long while...But I'm really busy with many things right know.... Our paper for the "Cultural Evolution in Spatially Structured Populations" have been accepted so we are currently working on that, and I have also stuffs for the PhD and a couple of papers I'm working on...busy busy busy.....
I've also started writing a couple of ABM using R, which sounds crazy at first (also at second, and third) but it has some nice things which I'll write about extensively in a future post.
But for know, I just wanted to start a series of very small posts (mainly for archaeologists) of small tips, which are astonishing simple concepts which however takes a couple of hours of googling and forum foraging...
For instance, have you ever tried to plot a time series of BC or BP dates? Suppose you have a sequence of count per century as follows:

data<-c(789,100,923,444,224,192,83,45,32,21,19,22,23,42,120)

plotting this as a timeseries is very simple:

plot(data,type="l")

and then you realise that you want something meaningful on the x-axis and you write the follow

dates<-c(3500,3400,3300,3200,3100,3000,2900,2800,2700,2600,2500,2400,2300,2200,2100)

 perhaps, if you know a bit of R you'll choose the more elegant nested function

dates<-sort(seq(2100,3500,100),decreasing=TRUE)

In any case you'll try to plot this as follow:

plot(x=dates,y=data,type="l")

and you'll find out that R ignored the ordering of the vector dates, and it even reversed your time-series.

My practical solution was to use negative values on the plot, and then delete the "-" with gimp or something (yes I should really be ashamed of myself).
Well for the small portion of people who had the same problem and here's the solution

plot(x=dates,y=data,type="l",xlim=c(max(dates),min(dates)))

Basically you can tell to the plot function that the range of values for the x axis is from the greatest value (the oldest date in our case, thus the largest number) to the smallest value. R will simply then read the values of dates in the correct order and plot the TS in the way it should look like.
Easy.