Deterministic-Continuous models vs. Stochastic-Discrete models

Recently, I've been in several workshops/discussions where several people were using continuos numerical models rather than ABM. A very common question from these guys are along the line of "why we should use ABM when we have mathematical models?". The argument is a known controversy (see McElreath and Boyd 2007 for a typical critique of ABM from a mathematical standpoint). Generally ABM is considered weak as it is computationally far more demanding, cannot be solved analytically and requires the exploration of a phase space (which sometime is not done, or not done properly) which is often too large. ABM users defend themselves by pointing out the importance of stochastic components, and the difficulty to represent complex behaviour which cannot be represented (easily at least) with a set of differential equations. My position stays somewhere in the middle, as I know how sometimes the temptation is to go for ABM due to my mathematical deficiency as a modeller. Ideally I try to invest on both types of modelling, and not to be epistemologically confined in one or another. Having said that one aspect which gives an advantage to ABM is the capacity to model dynamics with small population sizes.
While trying to verify the ABM of my doctoral project, I found a stunning example of this. In order to test if my ABM was working properly, I've created a mathematical equation of the same model, to see if their behaviour were consistent. Of course there were small differences, the mathematical formula was not discrete and non-stochastic. Thus for instance, if the population reproductive rate is 0.3, then the new population size will be:

old population size + old population size * 0.3

This means that if the original population size was 7, then the new population size would be 9.1. In the stochastic discrete model the different is that the population size is always an integer, and that the reproduction is basically modelled through a random pick from a binomial distribution as follows:

old population + Bin(old population, 0.3)

where the Bin() stay for a random pick of a Binomial distribution, with number of trials = old population, and success rate 0.3. In this case the most typical outcome will be the same (9.1), but there will be chances for smaller or larger number (you can visualise this in R, with hist(7+rbinom(100,size=7,prob=0.3);Notice that the mean of such frequency distribution will be roughly 9.1).


The first figure shows an example which compares a single run of the ABM (solid line) with the expectation derived from the numerical equation. Although the exact values are slightly different, the qualitative long-term outcome of the two cases are identical.

The second figure shows the same model with different parameter settings. The deterministic/numerical simulation expects an equilibrium around 4 individuals (actually 3.824 persons) while the discrete/stochastic model shows a second increase in population size at around timestep 170. The model actually continues such dynamic (not shown here) with other instances of growth and decline.

The third figure shows the very same model, expect that this time one of the parameter were set in a way that the population size was 25 times larger. Although the timing between the deterministic and stochastic model is slightly different, the overall shape, along with the long term equilibrium appears to be essentially the same.

So what is the point of all these figures? If you are dealing with low population sizes, the fact that your model is deterministic and continuous matters, and, depending on the model, matters a lot.  This is not just a problem of dealing with continuous numbers (e.g. 1.5 persons) but also because the consequences of small fluctuations are far more bigger in smaller population (the greatest real-world example is probably the effects of genetic drift). If you are dealing with large populations, chance events will be averaged out before it effects determines any qualitative change in the system. If you are dealing with smaller population, chance events might radically condition the output. This happens because some basin of attractions are so small that can be captured only when we deal with non-discrete population sizes. For instance the first figure had an equilibrium of 3.824. In a discrete model we have either 3 or 4 individuals, and both are likely to be outside of such basin. If you have a model which is scaled by x25 (as in the third figure) we will have an equilibrium of ca 95.59 individuals, and the discrete model can get closer to this (95 or 96), allowing to reach (and stay) within the basin of attraction. If the process we are trying to model is occurring at smaller population sizes, continuous model becomes unrealistic as decimal values are no longer an acceptable starting assumption. Of course you could tweak the mathematics and reason in discrete terms, but this becomes complicated and requires additional assumptions on how to convert continuous outputs into discrete values. In that case ABM is a robust alternative which can definitely give some contribution despite the known drawbacks.

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.

Multiple Y-axis in a R plot

I often have to to plot multiple time-series with different scale of values for comparative purposes, and although placing them in different rows are useful, placing on a same graph is still useful sometimes...
I searched a bit about this, and found some nice suggestions for 2 Y-axis, but haven't found anything for more 2+. So here's my solution:



#Create Dataset

time<-seq(7000,3400,-200)
pop<-c(200,400,450,500,300,100,400,700,830,1200,400,350,200,700,370,800,200,100,120)
grp<-c(2,5,8,3,2,2,4,7,9,4,4,2,2,7,5,12,5,4,4)
med<-c(1.2,1.3,1.2,0.9,2.1,1.4,2.9,3.4,2.1,1.1,1.2,1.5,1.2,0.9,0.5,3.3,2.2,1.1,1.2)


#Define Margins. The trick is to use give as much space possible on the left margin (second value)
par(mar=c(5, 12, 4, 4) + 0.1)

#Plot the first time series. Notice that you don't have to draw the axis nor the labels

plot(time, pop, axes=F, ylim=c(0,max(pop)), xlab="", ylab="",type="l",col="black", main="",xlim=c(7000,3400))
points(time,pop,pch=20,col="black")
axis(2, ylim=c(0,max(pop)),col="black",lwd=2)
mtext(2,text="Population",line=2)

#Plot the second time series. The command par(new=T) is handy here. If you just need to plot two timeseries, you could also use the right vertical axis as well. In that case you have to substitute "2" with "4" in the functions axis() and mtext(). Notice that in both functions lines is increased so that the new axis and its label is placed to the left of the first one. You don't need to increase the value if you use the right vertical axis.

par(new=T)
plot(time, med, axes=F, ylim=c(0,max(med)), xlab="", ylab="", 
type="l",lty=2, main="",xlim=c(7000,3400),lwd=2)
axis(2, ylim=c(0,max(med)),lwd=2,line=3.5)
points(time, med,pch=20)
mtext(2,text="Median Group Size",line=5.5)


#Plot the third time series. Again the line parameter are both further increased.


par(new=T)
plot(time, grp, axes=F, ylim=c(0,max(grp)), xlab="", ylab="", 
type="l",lty=3, main="",xlim=c(7000,3400),lwd=2)
axis(2, ylim=c(0,max(grp)),lwd=2,line=7)

points(time, grp,pch=20)
mtext(2,text="Number of Groups",line=9)

#We can now draw the X-axis, which is of course shared by all the three time-series.

axis(1,pretty(range(time),10))
mtext("cal BP",side=1,col="black",line=2)

#And then plot the legend.

legend(x=7000,y=12,legend=c("Population","Median Group Size","Number of Groups"),lty=c(1,2,3))

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

Le cadavre exquis boira le vin nouveau

I'm currently reading "Fooled by Randomness" by N.Talib, and I'm really enjoying a chapter dedicated on randomness in humanities... So let me start with two quotes:

" 'Reality is part of the dialectic of consciousness' says Derrida; however, according to Scuglia[1] , it is not so much reality that is part of the dialectic of consciousness, but rather the absurdity, and hence the futility, of reality."

"Sound is the change in the specific condition of segregation of the material parts, and in the negation of this condition; merely an abstract or an ideal ideality, as it were, of that specification [...]"

Very exciting. Now you have to guess: one of the two quotes is from Hegel, and the other one has been generated from a computer using the Postmodernism generator, created by Andrew C. Bulhak using the Dada Engine. (you can find the whole article I've created here, and if you want to make your own paper just click here).
Now, it's all fun, and if you know some philosophy you've probably guessed right (Hegel is the second one, also because it's very unlikely that he've quoted Derrida which was born 101 years after his death...). The point is, with another random combination (without Derrida being quoted) could someone have fooled you? To be honest I'm not a philosopher nor a postmodernist and I'll be concerned if you fool me with a scientific paper on spatial analysis rather than this (it would be nice to actually have a positivist generator to see how it would work...). I simply don't care much, but actually Alan Sokal, a physicist  of New York University, tried to submit a paper on Social Text by randomly creating (not with a computer this time) meaningless sentences... Well you can read his accepted and published paper here....Later in the same year, Sokal explains his point in another paper, (this time generated non-randomly I think):

"For some years I've been troubled by an apparent decline in the standards of intellectual rigor in certain precincts of the American academic humanities [...]. I decided to try a modest (though admittedly uncontrolled) experiment: Would a leading North American journal of cultural studies -- whose editorial collective includes such luminaries as Fredric Jameson and Andrew Ross -- publish an article liberally salted with nonsense if (a) it sounded good and (b) it flattered the editors' ideological preconceptions?"


and then states:

"Social Text's acceptance of my article exemplifies the intellectual arrogance of Theory -- meaning postmodernist literarytheory -- carried to its logical extreme. No wonder they didn't bother to consult a physicist. If all is discourse and ``text,'' then knowledge of the real world is superfluous; even physics becomes just another branch of Cultural Studies. If, moreover, all is rhetoric and ``language games,'' then internal logical consistency is superfluous too: a patina of theoretical sophistication serves equally well. Incomprehensibility becomes a virtue; allusions, metaphors and puns substitute for evidence and logic. My own article is, if anything, an extremely modest example of this well-established genre." (bold mine)

Scary, but perhaps true also in some archaeology and anthropology? While asking this question I've just read on Carl Lipo's blog how the executive committee of the  American Anthropological Association  are proposing to rephrase  their mission statement removing the word "science" (see here).

UPDATES: Simulpast, STDM and paper

Lot's of things happened the last few weeks so I need to update you on couple of things...
I've been in Barcelona last week and spent a wonderful time with the folks of the Spanish Research Council, where they allowed me to present two papers related to my PhD. They are just going to start a very exciting project called simulpast from early next year:
I quote from their website:

"The aim of the project is to develop an innovative and interdisciplinary methodological framework to model and simulate ancient societies and their relationship with environmental transformations. The propject will include 11 Research Groups of 7 different Institutions with more than 60 researchers from several fields (archaeology, anthropology, computer science, environmental studies, physics, mathematics and sociology). The leader institution is the IMF-CSIC in Barcelona."

I've never heard of any archaeological project centred on computational modelling having such broad range of case studies. This is a great opportunity and I'm really looking forward on the project outcomes. And I guess this can also be a great leap forward in terms of terms of standardisation and communicability of models. Good Luck and Thanks for the Tapas!!!

***

In the mean time I'm keeping myself (and Mark Lake) busy, as we are working for a paper which will go deeper on some of the topics we've explored for the CECD conference this September. We'll mainly focus on cultural transmission models of fitness-enhancing traits (2 and n-traits) with frequency dependency of the fitness and different types of Carrying Capacity (shared and independent) looking at short term dynamics and long term equilibrium of adoption rate and trait diversity. Stay tuned for more info!!!

***

The last but not the least! I'm quite excited, since there will be a International Symposium of Spatio-Temporal Analysis and Data Mining hosted here at UCL on July!!! This is a great chance to see many advanced techniques in spatio-temporal analysis and simulation which might give us some new perspectives in archaeology!!!

Inference From Confirming Evidence

Every once in a while I read papers in Archaeology claiming about there allegedly scientific methods failing to choose the right type of hypothesis testing for their models. While reading Taleb's "The Black Swan" I come across a very nice  psychological experiment conducted 50 years ago by Wason (1960).
The basic idea is that you have a data-set with a specific pattern, and in order to explain the underlying rule you conduct a series of experiments to test and propose a model. Of course in archaeology you cannot strictly do experiments all the time, but you can look for other data which will support or not support your model.
Now the experiment is based on a simple sequence of number, and one should simply "discover" the underlying rule. The player can propose another sequence of numbers and the experimenter will tell you whether such sequence can been generated from the same algorithm of the original data-set or not.

The sequence is:

2-4-6

Now, most people will most likely propose something like following sequence:

8-10-12

which will basically test the model increase by two. Now if the experimenter will tell you "Yes", you'll be probably quite happy about that, and probably you will write a paper for the Journal of Integer Sequences, which believe it or not, actually exists) with something lines:
"Our experiment confirmed our hypothesis of the increase by two rule". Very few will claim that this model is wrong and will be part of the scientific knowledge of your field (I'm sure that none of the editors of Journal of Integer Sequences will accept your paper, nor will fall on this trap...)
Now the big problem here is that the algorithm with which 2-4-6 has been generated was numbers in ascending order. Thus actually 5-6-7 would also have been accepted by the experimenter. However, it's vert likely that most people uses these experiments to confirm their model rather trying to falsify it. Testing 5-6-7 would have allowed a re-evaluation of the originally proposed rule and might have lead to the right answer.
The question now is: what sequence of number are you proposing in your archaeological research?


References:

• Taleb, N.N.,2007, The Black Swan: The Impact of the Highly Improbable, Random House, 
• Wason, P.C., 1960, On the failure to eliminate hypotheses in a conceptual task, The Quarterly Journal of Experimental Psychology, 12: 3, 129-140.