A simple illustration of an ecological lottery model in R
A colleague contacted me yesterday with some interesting modelling work based on the neutral community model of Stephen Hubbell.
Hubbell’s work fascinated me when it was first published. I have always assumed that it is easily misunderstood. If Hubbell is taken to imply that all species are equal, then the work appears to be eroding the basic subject matter of ecology. Ecologists tend to look for informative differences between species. However this is not exactly what the theory is based on. It is rather more subtle. The theory concerns the equality of individuals (in this sense it almost has political implications). Hubbell writes.
“The theory treats organisms in a community as essentially identical in their per capita probabilities of giving birth, dying, migrating, and speciating. This neutrality is defined at the individual level, not the species level. All that is required is that all individuals of every species obey exactly the same rules of ecological engagement. So, for example, if all individuals and species enjoy a frequency-dependent advantage in per capita birth rate when rare, and if this per capita advantage is exactly the same for each and every individual of a species of equivalent abundance then such a theory would qualify as a bona fide neutral theory by the present definition.”
I have illustrated this idea for students with a very small simulation in R . R code, like matlab or octave, is very terse and efficient so a couple of lines can do a fair amount of work. Lets assume a very simple reproductive rule applies to all individuals. The rule is that all have the same discrete generation time, in other words they all survive one time step of a simulation. At the end of the time step each individual produces two offspring. The lottery then applies. The offspring are randomly mixed and placed back on the finite space occupied by the previous generation. So half find a home and half “die”. The process is then repeated. This was exactly how I first implemented the model. I then re-implemented the idea in fewer lines by making the probability of selection depend on the proportion of individuals in each species, which amounts to the same thing.
This is all the code needed to run the model. (Try this if the quotation marks are not right lottery.doc)
nspecies<-8
gridsize<-20
time<-100
RUN<-function(){
X<<-matrix(0,nspecies,time)
palette(rainbow(nspecies))
mat<-sample(1:nspecies,gridsize*gridsize,replace=T)
X[,1]<-table(mat)
image(matrix(mat,gridsize,gridsize),col=sort(unique(mat)))
for (i in 2:time){
a<-X[,i-1]
mat<-sample(1:nspecies,gridsize*gridsize,prob=a,replace=T)
X[,i]<-table(factor(mat,levels=1:nspecies))
image(matrix(mat,gridsize,gridsize),col=sort(unique(mat)))
gc()
}
matplot(t(X),type=”l”,lwd=2,col=rainbow(nspecies))
}
#However to get a little tcltk interface add this
############################################
library(tcltk)
Run<-function(…){
nspecies<<-as.numeric(tclvalue(”nspecies”))
gridsize<<-as.numeric(tclvalue(”gridsize”))
time<<-as.numeric(tclvalue(”time”))
RUN()}
tt <-tktoplevel()
d1<-tklabel(tt,text=”Numero de species”)
d2<-tklabel(tt,text=”Tamaño del grid”)
d3<-tklabel(tt,text=”Tiempo”)
s1 <- tkscale(tt,from=0, to=20, variable=”nspecies”,
showvalue=TRUE, resolution=1, orient=”horiz”)
s2 <- tkscale(tt, from=2, to=100, variable=”gridsize”,
showvalue=TRUE, resolution=1, orient=”horiz”)
s3 <- tkscale(tt, from=10, to=500, variable=”time”,
showvalue=TRUE, resolution=10, orient=”horiz”)
but <- tkbutton(tt, text=”Run”,command=Run)
tkgrid(d1,s1)
tkgrid(d2,s2)
tkgrid(d3,s3)
tkgrid(but)
The model is to be played with (just download R from CRAN and paste all the code into the console), rather than watched but here are some example runs on You Tube,
The moral of the story is that under a neutral model some species start a random walk to extinction if co-existing in a small space with other species. However we cannot tell which will dominate and which become extinct from their characteristics, because all individuals, regardless of species are governed by exactly the same rules. The fate of species is determined by the cumulative luck or misfortune of the individuals that have been given the label. The larger the grid the more species that can be supported for a longer time, but extinction risk is always present if a long losing streak hits. The element of Hubbell’s theory that this model does not reproduce well is the log-normal (-ish) relative abundance distribution of the surviving species, although come to think of it I have not tried running the model with a really large number of species and a really large grid. R runs out of colours for the illustration, but this is not important.
Many interesting questions remain. Could a simulation model based on this sort of idea be used to predict extinction over a real landscape? Models of this sort seem far too abstract to provide real insight into the real world. But in some sense the process that Hubbell describes in a more formal mathematical way in his book probably does apply. But, how can these sort of ideas be communicated accurately to prevent them from provoking the sort of unproductive SLOSS debate that has often resulted from theory being pushed too far?
Rationality and the lottery
The BBC web site today contained what appears to me to be a misrepresentation of decision theory. The argument goes…
“Should you invest £2 a day or use it to buy lottery tickets?
Maths makes the decision obvious. Suppose you invest two quid every day at the reasonable rate of 10%. It will take you almost exactly 50 years to accumulate £1m. To earn this same £1m in the National Lottery, you would (on average) have to match five numbers and a bonus ball, at odds of 2,330,635-to-1.
If you spent two quid a day for 50 years you would total just over 36,500 tickets and would thus have only a 1-in-63 chance of making that million pounds. However, the available image of immediate wealth subverts this rationality.”
Is this right. Is it “obvious” as the author claims. No it is not. It is far from obvious.
The calculation of compound interest is correct, although banks do not normally compound interest on a daily basis and 10% is rather optimistic. You can check by simulating the arrangement as an R function using numerical integration.
f<-function(ndays=365*50,interest=0.1,value=2){
a<-numeric(ndays)
a[1]<-value
for (i in 2:(ndays)){
a[i]<-a[i-1]+value
a[i]<-a[i]+(interest/365)*a[i]}
a}
par(bg=grey(0.92))
plot(f(v=2),type=”l”,lwd=2,col=”red”,xlab=”Number of days”,ylab=”Accumulated value”)
grid(col=1)
The money in the bank grows healthily towards the one million target. So what is wrong with the argument? The author claims that the odds of winning a million on the lottery are 2,330,635:1. This is not a fair bet, but it is not such a bad one either. You have just under one chance in two million of winning the one million on offer. The expected value of your one pound ticket is the chances of winning (admittedly very small indeed) multiplied by the sum that would be won (and of course this is very large).
1/2330636* 1000000= 0.4290674
So the expected value of your ticket is about 43p. You have superficially wasted 57p.
The story about all the interest you would get by investing the money is a misleading red herring. If you took the conclusion of a 1:63 ratio between saving and gambling seriously it would persuade you not to buy a single lottery ticket even if the odds on winning were to become more favourable than one in a two million and bettered the value of the prize. Decisions between retaining a small sum with certainty and risking a big one do always involve subjective judgement, but few would not consider the lottery worth a shot if the prize of 1 million could be won at odds of (say) 200,000:1. The author of the article would (on this erroneous logic) still be convinced that it is better to put the money in the bank.
The formula for compound interest can be written as an R function in terms of the principal (p) the number of periods in a year that interest is paid (q), interest rate (i) and number of years(n)
f1<-function(p=1,i=0.1,q=365,n=1)p*(1+(i/q))^(n*q)
So using this function, lets think this all through calmly. If you were to win the lottery tomorrow and do the same with the money as you would have done with the two pounds you spent on the ticket, i.e. invest it at a compound interest of 10% you would be colossally wealthy in fifty years time. Using the same interest rate calculation that the author assumed you would have over 148,000,000.
f1(p=1000000,n=50)
[1] 148311560
On the other hand, if you were to win your million exactly fifty years from now you would just have your million at the end of the period. This would coincide with what you would have gained from saving.
So to reiterate, all wins before the final date are worth more than the saved money in the bank at the end of fifty years, The earlier you win the better. The only addition I have made to the authors’ own argument is to assume (quite fairly) that lottery winnings also gain compound interest. The comparison the author makes between the frugal saver and the lottery player is quite unfair. It uses only the absolute minimum that a lottery win would be worth as the baseline for comparison. The expected lottery winnings at the end of fifty years are quite clearly worth very much more than one million. In fact under this model, it is easy to show that the expected amount is exactly 0.4290674 times the money that would be in the bank if you had not played the lottery, providing comparable assumptions are made regarding the use of the money.
plot(f(v=2),type=”l”,lwd=2,col=”red”,xlab=”Number of days”,ylab=”Accumulated value”)
lines(f(v=2*0.4290674),type=”l”,lwd=2,col=”blue”)
The differences between the money paid and the in expected value (in purely monetary terms) doesn’t change.The ratio between the red line (saver) and the blue line (expected value from playing the lottery and investing the proceeds) stays the same. Lottery players have (on average) an expected value of around 43% that of the savers. They are worse off, but nowhere near as irrational as the article suggests..
But we can go a step further with the argument. As you think it through and apply common sense it gets better and better for the lottery player. You clearly wouldn’t ever dream of actually investing a million you won tomorrow in order to have megabucks in fifty years time. A small fortune is worth much more to you now than an unspendable fortune in the future. In fact, to you, it is almost certainly worth much more than 148 times its future value, given the positive, life enhancing, potential of a single million. After the first million the next 147 are increasingly irrelevant to your happiness. This could be written using a function that converts money into happiness. This is a curve that reaches some sort of asymptote. The absolute level of the asymptote varies between individuals, but the shape is fixed, even for Bill Gates.
At the same time the savings should be devalued by the probability of dying before they can be used, the bank suffering a fate worse than Northern Rock, a meteorite strike, or the consequences of catastrophic global warming among a multitude of other scenarios. Depending on just how much all these trade offs come in at (which again is a rather subjective matter), your lottery ticket could easily turn out to be worth more to you than the pound you paid for it.
The original article states ..
If you spent two quid a day for 50 years you would total just over 36,500 tickets and would thus have only a 1-in-63 chance of making that million pounds.
A 1-in-63 chance during a lifetime doesn’t sound so unlikely!
It can in fact be perfectly rational to buy a lottery ticket. Which is why so many rational people do.
End line: Why should this interest a forest ecologist? Because it might explain why sustainable forestry is so difficult!