Inclusive fitness errors and divorce in birds

In Darwin’s theory of evolution by natural selection the organism was king – natural selection would act to favour the traits that helped an organism to survive and reproduce. Whilst this idea explained so much of the adaptation and apparent design that we observe in the natural world, it also presented a problem. Why do organisms such as workers in social insect colonies sacrifice their own reproduction to help the queen reproduce and raise her offspring? The solution to this came from the work of Bill Hamilton, later clarified and promoted by Richard Dawkins in The Selfish Gene. Hamilton showed how the seemingly paradoxical altruistic behaviour of organisms (like workers helping the queen) made sense if we considered the genes involved. The great insight was to show how natural selection can favour an organism to sacrifice its own reproduction if it can provide benefits to related individuals, with whom the organism has a predictable chance of sharing genes. In short, a gene for altruism can spread even if it causes the individual carrying it to sacrifice its reproduction, as long as those effects are outweighed by the benefits the altruistic act provides to a relative, modified by the degree to which the altruistic individual shares genes with the beneficiary. This logic forms the basis of the famous ‘Hamilton’s rule’, which shows that a gene for a social behaviour such as altruism will spread if rb-c>0, where c and b are the costs and benefits of the social behaviour on the actor and recipient respectively, and r is the relatedness of the actor to the recipient. As a simple example, imagine that I can pay a cost c=1 to provide a benefit b=3 to another individual. This behaviour will be favoured if rb-c>0, or 3r-1>0 in this case, which simplifies to r>1/3. In this case, the behaviour would be favoured if the recipient was a brother (r=0.5) but not a half-brother (r=0.25).

Hamilton’s ‘inclusive fitness’ theory that is captured simply by Hamilton’s rule has been extremely useful in explaining so much of the complex cooperation (and conflict) that we see in nature, partly because of the simple logic of Hamilton’s rule, and the basic idea of individuals valuing other individuals according to the degree of relatedness between them (e.g. Haldane’s idea of jumping into a river and sacrificing yourself if you can save 2 siblings, or 4 half-siblings, or 8 cousins etc.). However, the simple verbal logic can sometimes lead you astray. A recent series of papers by Caro (2016), Bebbington & Kingma (2017), and Levin (2019) provide a great example of this, which I will discuss here.


Imagine a species of bird which is monogamous and produces one offspring per year, which both parents work to feed at the nest. The chick signals its desire for food, and could demand more food than it really needs to provide the maximum benefit to itself, or could demand less. By demanding less food than it could, the chick is paying a cost, but the consequences of doing so are to increase the chance of the parents surviving to breed the following year – producing a sibling (r=0.5) which shares genes with our altruistic chick. As such, we could imagine a scenario where the costs (to the chick) and benefits (of allowing the parents to produce a chick the following year) are such that this behaviour is favoured by natural selection.


With the basic scenario now set-up, we can ask interesting questions about what would happen if we changed certain aspects of our scenario to fit in with what we know about the kinds of behaviours we see in birds in the real world. Imagine first that our chick pays the same cost to be prudent with demands for food, and provides the same benefit of survival odds to the parents, but that the female has mated with multiple male (as is the very common in birds). The consequence of this multiple-mating would be that the offspring produced might be only a half-sibling (r=0.25) of our altruistic chick. Consequently, the altruism in the chick may no longer be favoured.

So far, our verbal arguments of the logic of Hamilton’s rule have served us well, and everyone would agree with the analysis shown here for these scenarios. However, what if we make a different change to our assumption about how the birds pair up. Imagine now that the scenario is the same, expect that our bird parents now ‘divorce’ after they have raised our altruistic chick, and (thanks to the altruism of our chick not demanding too much food) both survive to reproduce the following season with new partners. Would this effect the chance of altruism being favoured?

We might answer this problem by considering how the altruistic act of our chick is now helping the mother to survive and reproduce the following year – producing a half-sibling of our chick with her new partner (r=0.25). The chick’s altruism also helps the father survive and reproduce – producing a half sibling with his new partner (r=0.25). Helping its parents produce two extra half-siblings (in the case of divorce) is equivalent to helping the parents produce one extra sibling (if they stay together), so we would argue that the divorce has made no difference to whether the trait is favoured. At first glance it is hard to fault this explanation: it sounds logical, and it sounds like it has correctly interpreted the effects of producing siblings versus half siblings. When I was a student just becoming interested in social evolution, I’m sure I would have completely accepted this logic. However, this answer is wrong*, and I think it’s a really good example of where the simple logic of Hamilton’s rule can mislead us.


To explain why the idea that divorce wouldn’t make a difference in whether altruism is favoured is wrong, we can go back to what Hamilton himself pointed out about his own theory. The critical point is that inclusive fitness isn’t about calculating all of the offspring produced by relatives, rather it is about those which are caused by the behaviour. This is important because what we are aiming to do is calculate whether a rare mutant allele for altruism will spread in our population, which will depend on the effects that are specifically caused by that allele, not on the overall fitness of the individual’s relatives. In order to make this calculation, we have to be very careful about what effects can actually be ascribed as being due to the behaviour of the allele. Going back to our scenario, we can illustrate this. Without divorce, our altruistic chick is giving a benefit (let’s call it B) to each parent, and those benefits are combining (i.e. B+B=2B) to help those parents produce one extra offspring (a sibling of our chick) the following year. With divorce, our altruistic chick is giving a benefit (B) to each parent, and each parent is apparently able to use that to produce one extra offspring (a half-sibling) the following year. The error is that without the divorce we are saying that the parent only needs B to produce an offspring the following year, whereas with divorce we said they needed 2B to produce an offspring. Of course, the extra B comes from the new partner of the mother and father – but those effects are not caused by the altruism of our chick, and so aren’t part of the fitness effect of the trait. Instead, they are part of the fitness effect of some other trait, perhaps a different (competing) allele in a different chick. As such, we can only attribute half of the production of the half-sibling as an effect of the altruism gene. To count the extra B as part of the fitness effect of altruism would be to commit a so-called ‘double counting’ error, where a fitness effect has been attributed to multiple different genes. This error becomes obvious in formal mathematical treatments of inclusive fitness theory, but can often be obscured in seemingly sound logic of verbal arguments based on the simple idea of Hamilton’s rule. I think that the example of the consequences of ‘divorce’ on selection for altruism provide a compelling example of this, and one which could be very instructive to people learning about social evolution

*It is important to note that the argument is wrong here for the example we are considering – which is very much focused on birds, for which parental care is hugely important and constitutes a major part of the investment required to successfully raise offspring. If there was no parental care, and the important step in producing offspring was simply surviving to produce a successful gamete, then divorce wouldn’t make a difference. This is because the benefit provided to the parent by the offspring (B) would be sufficient to produce the half-sibling, which could be fully attributed as an effect of the altruism gene in our focal chick. In this way, the answer to the question of whether divorce matters for altruism depends very much on the biology of the system in question – another important aspect that should always be considered. However, in the case of birds it is very clear that the limiting factor is the investment of the parents in feeding the chick, not simply the production of gametes.

Football Crazy, Probability Mad: How much does it really cost to complete the World Cup 2018 Sticker album?

Every four years football gives us the World Cup; the premier showcase of global talent. At the same time, the Italian company founded by the Panini brothers gives us something far greater. Panini produce a sticker album with hundreds of stickers to collect, representing the players, stadiums and nations at the World Cup. In doing so, they give us an intriguing mathematical challenge: how can we collect all 682 stickers as cheaply as possible?

Since the release of this tournament’s sticker album there have been several articles claiming to calculate the cost of completing the collection. Most of these articles make the same basic mistakes, and are hugely misleading. Always remember to check the assumptions and logical basis of any of these analyses. This post is based on the facts and assumptions listed below

1) There are 682 stickers to collect
2) Each pack contains 5 stickers and costs £0.80
3) Each sticker is equally likely to occur in any pack
4) Individual stickers can be purchased on the Panini website for £0.22 each

Point 4 is ignored by almost everyone who attempts to calculate the cost of completing the collection – but it makes a huge difference.

With this in mind, let’s get started. What is the simplest way that someone can collect all the stickers? Keep buying individual packs until you have them all! But how many packs do you need to buy? We can answer this through simulations. We use random numbers to randomly pick 5 from the 682 available stickers to create one pack of stickers. Then we count how many of those stickers are unique – not duplicates. We then randomly pick another 5 stickers to make another pack, and add those to our simulated collection. We can continue like this until our individual has collected all 682 stickers. The beauty of doing this on a computer is that we can then replicate this a large number of times, to work out how long it would take to complete the collection on average. I did this 100,000 times and have plotted the results below


I found that on average you would have to buy 968 packs to complete the collection! Which would cost a total of £774.40, more than most would be willing to spend on a football sticker album.

Why do we need to buy so many packs to complete the collection? Buying 968 packs would give us 4840 stickers, of which only a meagre 14% would be unique.

The best way of seeing why we have to buy so many packs is another simulation. But this time, instead of simulating the continual buying of packs until the 6820 total is reached, we will simulate the buying of a number of packs equal to n, where n is any number between 1 and 1000. We can then record the number of unique stickers we expect to get for each number of packs we buy. This is repeated 10000 times for each number of packs, to get a good estimate of the average number of unique stickers obtained from each number of packs bought.


As you can see, we get a nice curve. And this makes good logical sense. At the start, for each new pack we buy its quite likely that all the stickers we get will be new to us, and contribute to our collection. Towards the end of the curve when we have already got most of the collection, each new pack is likely to contain only stickers that we already have. Note that we already have most of the stickers when we have bought 500 packs, but on average it took an extra 468 packs purchased to complete the collection in our simulations from earlier!

At this point do we give up hope, or start saving up more than £700 and preparing to empty the local Sainsburys of its supply of stickers? Well not quite yet. At any point we can go online and order individual stickers, the ones we need to complete the collection. Why don’t we just do this from the beginning, you may ask. Well aside from that being blatant cheating, individual stickers cost £0.22 each online, compared to the £0.16 they cost when you buy them as part of a pack. However, it seems clear that with this information we no longer have to buy 968 packs. As we have already noted, by 500 packs we already have most of the collection.

So with this new information, the key question becomes ‘at what point should I stop buying packs, and buy individual stickers instead?’. We can easily calculate this by taking our figure above, for the number of unique items we get for each number of packs bought, and converting this into a final cost, but adding the cost of buying that number of packs to the cost of buying the number of stickers required to complete the collection.


The y-axis of this graph is cost, so we can look for the point where this is lowest to find the optimal number of packs, and the cost of buying that number of packs. The answer is 44 packs, for a total cost of £143.88. A much more reasonable answer than 968 packs for a cost of £774.40

So far we have been using simulations to answer this question, but can we answer the same question with simple maths? Intuitively we know that we are dealing with simple probabilities, the probability of getting any sticker by picking one at random is 1/682. Can we however use this to calculate how many unique stickers we get for each number of packs bought? Yes we can.

This kind of situation is known in maths as the ‘Coupon collectors problem’. Let’s use p to represent the probability of getting a new unique sticker to add to our collection for each new sticker we buy. The chance of getting a new unique sticker from our first one is 1, since we haven’t seen any other stickers so there is 0 chance for a duplicate. The chance of getting a new unique sticker from our second pick is 681/682, since there is one option that we could pick that might be a duplicate. Similarly, for the third pick the chance of adding to our collection is 680/682. This can be formalised to say that the probability of getting a new unique sticker on any draw i = 1 + (n-1)/n + (n-2)/n + … and so on until we have a number of terms equal to i. Can we use this simple maths to work out the time it takes to reach any defined number of unique cards? Yes we can! time and probability can be simple inverses of each other, if defined carefully. For example, if I said that the chance of rolling an even number on a dice is 1/2, we would expect to take 2 rolls to see an even. So we can use our probability equation to calculate the time t to collect any number of unique items = 1 + n/(n-1) + n/(n-2) … and so on. We must be careful to include the fact that we don’t buy individual items, we buy stickers in packs of 5, so we have to divide this time by 5 to convert it to number of packs. Using this we can recreate the mathematical version of the simulations we did earlier, and find it to be indistinguishable


Next, we can see if this gives a different answer to the ultimate question of how many packs we should buy


Here, we can see that this solution is exactly the same as our simulations! We should still be buying 44 packs and then buying the rest of the stickers online. Our simulations, and our simple maths, are in agreement.

So is that the final word, shall we put some of our savings back in the bank and prepare to spend £144? Well actually that’s not the end of it either.

We are a social species, and we love to trade. So I could bring in a partner to collect stickers with me. And we would both benefit because I could give her my swaps for free, and receive her swaps for free!

Before you continue, have a quick think and make a prediction about what effect this will have on the optimal number of packs we each need to buy to complete the collection. We will assume that both partners buy the same number of packs, and give all of their duplicates to the other player. Compared to the optimal number of packs for a player on their own, do you think we should buy more packs when playing with a partner, or less? Why do you think this?

For me, the easiest way to think about the effect of another player is that it gives us a loads of extra packs for free. For any number of packs we buy, we inevitably get a certain number of duplicates. Our partner has the same. Those duplicates are then packaged up and given to the other player as free packs. This has the effect of shifting the curve for the number of unique cards we get for any number of packs bought


In red here we can see how the black curve we had before shifts when we take into account all the free cards we get from our partner. So what is the effect on the optimal number of packs we should buy?


In red here we can see the new cost function, how much it costs to complete the collection for any number of packs bought (before doing swaps and buying the rest individually). Contrary to my original expectation we have to buy more packs when playing with a partner than we do when alone! 74 packs compared to 44 when playing on our own. The key however is that the total cost to complete the collection is lower; £136.20. Buying more packs each gives us more swaps to give to the other player, allowing us to obtain a larger collection before we have to buy the comparatively expensive individual stickers online for £0.22 each.

This calculation was based on our maths based calculations, but I’m a big fan of checking these using simulations, where we can be confident that we are recreating the effect we want to observe, and check if we have made any big maths mistakes. The graph below shows the average cost to complete the collection based on the number of packs bought when there are two players cooperating.


This simulation actually has some element of randomness, in this example 133 packs was on average the optimal number of packs to buy, but repeating this leads to slightly different answers. The important thing is the shape of the curve, which is lowest at around 74 packs, much higher than the 44 packs for a single player and in full agreement with our mathematical calculation.

So there we have it, you should grab a friend and each buy 74 packs, then give all your duplicates to your partner, take their duplicates, and buy the rest of the stickers online. Incidentally, you should expect to have to buy around 349 stickers online if you want to complete the collection for as little money as possible (these can be purchased on the Panini website in batches of up to 50). You may wish to compare your progress as you go against the mathematical expectation. By pure luck you could end up doing better than expected, and benefit from stopping buying packs much sooner to save more money.

The remaining question is what happens if you have not one partner but two or three, or even ten. I have used the same principle of swapping and simulated this scenario, with the result shown below.


Here we find that with 3 players we actually don’t end up buying any more or less packs than we do with 2 players. But why is this? shouldn’t adding more members to our syndicate of sticker collectors make a larger difference than this? Again we have to check our assumptions, and how realistic they are. Up until now we have instructed our simulated to players to use a simple swapping system, whereby each player simply passes all of their duplicates to the next player. We could visualize this by imagining a 10 player system where all players sit in a circle, and pass all their swaps to the player to their left. You can easily see the problem with this system. Many of the swaps that you get won’t be helpful to you, because you will already have that sticker. Players can be much savvier with their swapping, trading individual swaps to any of the other players for specific individual stickers that they still need. In this final section I have named this system ‘strategic swapping’ and enacted it by having all players pool all of their swaps together, and take turns strategically picking one sticker at a time that they need to further their collection.

I have done this for 2, 3, and 4 players and shown the main result below.


In black you can see the number of unique stickers achieved for each number of packs bought when collecting alone. In red, blue and green we see 2, 3, and four players respectively. As we add more players to our syndicate of strategic swapping collectors, we increase the number of unique stickers we get for any number of sticker packs we each buy. However, for each extra person we add, the benefit of adding that player is less than the benefit of adding the previous player.

How does this impact the amount of money needed to complete the collection? For two players we needed to purchase 74 packs and then swap and buy the rest online. For three players we need to purchase 82 packs, and for four players it’s 91 packs. Switching from two to four players brings the total cost down from about £136 to £130. A modest, but worthwhile saving.

As a final note, I have attached here (see the link below) a spreadsheet that all players can use to track their progress. For each pack purchased simply add in the total number of unique stickers you now own (without any swapping) and you can see how you compare to the mathematical expectation. You can see that my progress for the first few packs is right in line with that we expect. May the odds be ever in your favour.


St Petersburg Paradox – Code Efficiency

The Paradox

The paradox involves a casino game which works as follows. A (fair) coin is flipped, and if it lands on heads, the game ends here and the player wins £2. If the coin however lands on tails, the game proceeds to round two. In round two the coin is flipped again, and if it lands on heads, the game ends with a payoff of £4 to the players, and if it lands on tails the game proceeds to round three. The potential payoff to the player doubles in each round, whilst the probability of this happening halves (due to the increasingly unlikeliness of flipping many consecutive tails). In the original paradox this can continue for an infinite number of rounds.

##      Round Payoff Probability
## [1,]     1      2     0.50000
## [2,]     2      4     0.25000
## [3,]     3      8     0.12500
## [4,]     4     16     0.06250
## [5,]     5     32     0.03125

The paradox arises when you ask how much should you be willing to pay for a chance to enter this game? You have half a chance of winning £2, a quarter of a chance of winning £4, and eight of a chance of winning £8 and so on until infinity. £2? £20? £20000?

Expected Value

To answer this, we use the expected value of the game. To illustrate this we can use the simple example of rolling dice. Let’s say the dice is not fair, so there is half a chance of rolling a six, and a 10% chance of rolling each of the remaining five numbers. Let’s then say that the payoff for rolling a six is £20, and the payoff for rolling any other number is 0. So 50% of the time we get £20 (for rolling a six), and 50% of the time we nothing (for rolling any other number). The expected value of this game is the amount we win if each event occurs, weighted by the chance of each of these events occuring. So here, our expected value is 0.5 x £20 + 0 x £1 , which equals £10. This is how much we expect to win (on average) by playing this game.

St Petersburg Expected Value

Going back to the St Petersburg Paradox, our expected value is again the amount we win from each event, multiplied by the chance of each event occuring. So for the first few events (shown in the table) this is (2 x 0.5) + (4 x 0.5) + (8 x 0.125) + (16 x 0.0625) + (32 x 0.3125). We then notice that each of these values = 1, because the probability of winning any amount is 1 / the amount. As such, our expected value, the amount we expect to win on avergae, is 1 + 1 + 1 + 1 + 1. If the game goes on until infinity, we can see that the expected value is 1 + 1 + 1… all the way to infinity, so we should be willing to pay an infinite amount of money to play this game! This is clearly a ridiculous answer, but what would be a more reasonable answer? The 18th century French mathematician Comte de Buffon approached the problem by paying a child to play 2048 trials of the game, but we can do better by making a computer play the game.

The Finite game

The game we are going to play is a slight variation of the St Petersburg paradox in that it is finite meaning the casino has only a certain amount of money to give you. We will call this W. This means that there is a maximum number of rounds, which we will call N which we can calculate as the rounded down (or floor) of log2 W. For example, if W = 20, N = floor(log2 20) = 4. In this instance, 50 doesn’t fall on the normal payoff scale (it is between 16 and 32) so the rules of the game state that if the game reaches the round where the casino has to pay out more than the size of their bank (W) they simply pay out their entire bank.

The code

So what is the intuitive way of approaching this problem? It is to do exactly what Buffon did ~250 years ago, but with a computer simulating trials of the game, rather than paying a child to manually flip a coin. To start with we can set our two parameters, which is the number of simulations sims and the size of our casino’s bank W, and calculate the maximum number of rounds N

  sims <- 1000000      # Number of simulations
  W <- 1000000          # size of bank
  N <- floor(log(W, 2)) # maximum number of rounds

First we can make a simple vector to store the outcome of each trial

output <- c()

Now we can make the loop that will go through each round number, work out what the outcome of the flip of the coin in that round is, and either record the payoff corresponding to that round and stop (if we flip heads) or proceed straight to the next round (if we flip tails). We can use a random number between 0 and 1 to simulate a coin flip (if it is > 0.5, we call that heads) and the payoff is simply 2i where i is the round we are currently in (i.e., if we flip heads in round 3, the payoff is 23 = 8)

 for(i in 1:(N+1)){
    # Flip coin, if heads store payoff in output matrix
    if(runif(1) >= 0.5){output <- 2^i} # payoff = 2^round
    if(i == (N+1)){output <- W}        # Or whole bank 
                                  # if we reach the last round
    if(output >=2){break} # if we win, stop flipping
## [1] 2

This will give us a payoff for one trial of the game. In order to make many attempts of the game, we need an outer loop to run through this procedure many times. For this we will need a matrix to record the output of each trial

output <- matrix(0, sims, 1)

We can then make our loop. Note that we have altered the inner loop slightly to record in a certain cell of the matrix, as we will be generating a new payoff for each trial. For a more detailed explanations of using loops, see here.

for(j in 1:sims){ # Outer loop j for many trials 
  for(i in 1:(N+1)){ # Inner loop i for each round
    if(runif(1) >= 0.5){output[j, 1] <- 2^i} # payoff = 2^round
    if(i == (N+1)){output[j, 1] <- W} # Or whole bank
                                    # if we reach final round
        if(output[j, 1] >=2){break} # if we win, stop flippinh
mean(output[, 1]) # output the mean payoff of all trials
## [1] 20.49615

We can time how long this takes using the tictoc package which is really simple to use. First use the following code to install (if required) and load the package

list.of.packages <- c("tictoc")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

To use it simply write tic() where you wish to start timing, and toc() at the end. Now we can see how long it takes to run 1 million trials of the game, when the maximum payoff is £1000000.

tic()   # start timing here
output <- matrix(0, sims, 1)
for(j in 1:sims){ # Outer loop j for many trials 
  for(i in 1:(N+1)){ # Inner loop i for each round
    if(runif(1) >= 0.5){output[j, 1] <- 2^i} # payoff = 2^round
    if(i == (N+1)){output[j, 1] <- W} # Or whole bank
                                    # if we reach final round
        if(output[j, 1] >=2){break} # if we win, stop flippinh
mean(output[, 1]) # output the mean payoff of all trials
toc()   # stop timing here

Probably ~15 seconds, not bad! But we can do so much better. We have approached the problem like a human, not a computer. We’ve logically programmed the steps that need to happen, in the order that they need to happen, but much of this is unneccesary. We are wastefully creating many random numbers to decide whether we should progress to the next round, and using a loop to go through every round of the game. A computer can do things far more effectively. Let’s start again by redefining our starting parameters

sims <- 1000000      # Number of simulations
  W <- 1000000          # size of bank
  N <- floor(log(W, 2)) # maximum number of rounds

Now let’s ask; what are the potential payoffs of this game? well we know from earlier that the formula for the payoff of each round is simply 2iwhere i is the round we are currently in. We can easily generate this list of possible payoffs, 2^seq(1:N) will generate this. We must also remember that the final possible payoff is simply the size of the casinos bank, W, which we can add to our list by using a vector

payoffs <- c(2^seq(1:N),W)

We can now sample from these payoffs to determine the payoff of each trial, without needing to simulate up to N random coin flips. We could do this using the sample function, for example like this

winnings <- sample(payoffs, sims, replace = TRUE)
## [1] 102837.7

Here we are generating a sample of the potential payoffs, doing so for sims number of trials, and replacing back each payoff into the available pool once we have sampled it.

The problem of course is that it is not equally likely that we win each of these payoffs, so we can’t just randomly sample them. However, each payoff has an associated probability, which we can use to make a biased sample of the payoffs that is exactly equivalent to playing the game many times, without having to waste time generating random coin flips; we can simply generate a single number, the payoff, for each round.

The biased sample is generated by simply attaching probabilities to our sample, a simple example could be as follows:

payoffs <- c(0, 1) # vector of possible payoffs
probs <- c(0.2, 0.8) # vector of associated probabilities
# we are making the probability of winning 0 = 0.2
# and the probability of winning 1 = 0.8
sample(payoffs, 10, prob = probs, replace = TRUE) # 10 samples 
##  [1] 1 1 1 0 1 1 1 1 1 1

We can use this for the St Petersburg paradox; we already know what the vector of possible payoffs is payoffs <- c(2^seq(1:N),W) and the probabilities are simply 1 / each payoff probabilities <- c(1/(2^seq(1:(N+1)))). Putting these into the sample function, we can then make our biased sample, with the size of the sample equal to the number of trials of the game we want to run, and then simply calculate the mean of that sample. We will time this code with tictoc to see how quick it is comapred to the previous method

  payoffs <- c(2^seq(1:N),W)
  probabilities <- c(1/(2^seq(1:(N+1))))
  winnings <- sample(payoffs, sims, prob = probabilities,
                      replace = TRUE)

It is several hundred times quicker! And we can even improve a little on that, by running the code on a single line

mean(sample(c(2^seq(1:N),W), sims, prob = c(1/(2^seq(1:(N+1)))), replace = TRUE))

This approach is so much faster, with the key being to think about how you can recreate a ‘real’ process like flipping a coin many times, with a mathematical equivalency, like taking a sample of possible outcomes, biased by their likelihood. Although this approach would be silly if we were to manually play the game; it would be like the French mathematician taking a million pieces of parchment, and writing the payoffs on a proportion of them equal to the likeliood of them occuring, it works so much faster when coding. If you try bumping up the number of trials you will see that the time difference between the two approaches is even greater.

A very useful function: eval parse text paste

A very useful function: eval parse text paste

When creating evolutions models and simulations in R, code can quickly become messy and slow. Functions that require fewer lines of code can speed up your program, and make it easier to follow what is happening. There are also often variables, such as number of alleles or loci, that need to be manipulated easily, without the requirement of changing multiple lines of code. This is where eval parse text paste becomes useful. The four stages of this can be understood by considering each of them individually, in reverse order.

paste simply converts something into text, and can combine several pieces of text together. For example. In the code below, we can paste the generation number of a program i next to the word generation. ,sep="" means there is no gap between the two pieces we are pasting together. We could make this to something else such as an underscore or a space if we wished.

i <- 1
paste("Generation ",i,sep="")
## [1] "Generation 1"

text takes this pasted expression as the input for parse, which converts it into an expression. For example take the code below. We are making a vector Output to store something, and creating an expression that can multiply the result for our value of i (generation) and multiply it by 50

i <- 1
Output <- c()
Output[i] <- i+1
## expression(Output[1] * 50)

eval then evaluates this expression, meaning it calculates the result of the expression, for that value of i

## [1] 100

We could use this to give us a Result for each generation i using the Output. The basic format is eval(parse(text=(paste("something",i,"something else",sep="")))).

Output <- c()
Results <- c()
for(i in 1:10){
  Output[i] <- i+1
  Results[i] <- eval(parse(text=(paste("Output[",i,"]*i",sep=""))))

This is a basic application that could be completed with other methods, but there are much more useful features of eval parse text paste.

At this stage I find it helpful to create a function to contain eval parse text paste. The code at the start of each line using this expression is quite long, making for some long lines of code that can be hard to read. Making functions in R is really easy. I won’t go into detail, but here is a good link for the basics of making your own simple fuctions

This function here is called belch3 because it is an eval expression with 3 pieces.

belch3 <- function(x, y, z) {
   eval(parse(text=(paste(x, y, z,sep=""))))

You can make more functions for arguments with 4 or 5 pieces, which are sometimes necessary

belch4 <- function(x, y, z, a) {
  eval(parse(text=(paste(x, y, z, a, sep=""))))

belch5 <- function(x, y, z, a, b) {
  eval(parse(text=(paste(x, y, z, a, b, sep=""))))

We can now start to explore a useful application of these functions, in creating the right number of variables for a scenario. Imagine you are creating an evolutionary model and need a vector for each locus in your model, but each time you run the model you will use a different number of loci. Manually creating the right number of variables is tedious, and unneccessary with this function.

number.loci <- 3
for(i in 1:number.loci){

We now have the correct number of variables produced each time, simply by changing the variable number.loci. Note that when creating variables we use <<- instead of just <-. This is to ensure that the variable is available to us, rather than just contained within the function. See this link if you want to find out more.

An example of using one of our longer functions is if we want to make another variable that uses information from one that we have just made. For example, here we are making a matrix for each loci, using information from that locus and random numbers (for no meaningful reason, simply to show the kinds of things that could be done)

for(i in 1:number.loci){

This function is extremely useful for generating and storing variables where you are varying something such as allele number or loci number, where each of these variants require their own vectors/matrices that would normally have to be generated manually, with many more lines of code. If you don’t want to use the function, you can just use the full eval argument on each line

for(i in 1:5){

but this is quite a bit longer, and can be more difficult to read, so I prefer to use functions.

There are many other uses for this function, with another key one being for plotting graphs.

Here, we are making a seperate graph for each locus, plotting the loci.output matrix. We don’t have to make each graph individually. This function will plot the correct number of graphs from number.loci

for(i in 1:number.loci){
  belch5("plot(seq(1:25),loci.output",i,", col = i, type='l', xlab = 'Generation', ylab = 'Loci Output', main='Locus ",i,"')")

We could also plot more than one line on each of these graphs, for example for different alleles at each loci. For this purpose I am going to use the eval function to generate a matrix for each allele, with a column in that allele for each loci and a row for each generation.

number.allele <- 4
for(i in 1:number.allele){

We can now plot the lines for each allele onto the graph for each loci. All in a very short section of code, that is fully flexible to changes in allele number and loci number

for(i in 1:number.loci){
  belch5("plot(seq(1:25),loci.output",i,", col = i, type='l',lwd=3,ylim=c(1,2), xlab = 'Generation', ylab = 'Loci Output', main='Locus ",i,"')")
  for(j in 1:number.allele){
    belch3("points(seq(1:25),allele",j,"[,i], col=j, type='l')")

Does cycling really slow ageing?

In the science section of several media sites this week was a story about ageing and exercise. The headline in The Independent was ‘The secret of eternal youth: skin-tight Lycra and a bicycle’ and the story also appeared on Yahoo headlined ‘Study: Cycling slows ageing process’. A game I frequently like to play when viewing science in the media is look at the headline and imagine a rough plan of what would be required to test that statement or questions. Frequently the media have answered a question that the researchers themselves haven’t even asked. In this instance I imagined a study where a series of people had been followed throughout their lives, and cycling had been isolated as a factor that slowed the anti-ageing process, without being confounded by other factors that cyclists are likely to share such as a general active lifestyle. Cycling would also have to be compared with other forms of exercise to justify its solitary elevation as an anti-ageing tool in the media articles.

The original research paper that forms this story is an early online article from the Journal of Physiology ‘An investigation into the relationship between age and physiological function in highly active older adults’. The main thrust of the paper is about finding ‘biomarkers’ for ageing. This is basically identifying physiological factors that can be used to accurately age an individual. The main purpose of such a pursuit is to test possible interventions to slow the effects of ageing. Having a baseline where individuals could be accurately aged would allow any anti-ageing effects to be quantified, and compared against each other for the magnitude of the effects. The identification of biomarkers for ageing has proved problematic so far, as many factors affect the physiological decline with age we see. Genetics and lifestyle are two particularly important factors that result in candidate biomarkers being inaccurate predictors of age. The authors here argue that lifestyle is particularly important as a confounding factor, with a sedentary lifestyle having large impacts on physiological function that stop any biomarker being able to accurately age individuals, as in sedentary individuals physiological decline occurs due to that lifestyle, as well as due to time.

To resolve this issue the authors decided to use a subset of people, who lead active lives. The reasoning for this is that in theory these individuals have not damaged their physiology through inactivity, so it should be possible to predict certain factors based on their inevitable decline with age. Elderly cyclists were chosen as a rather arbitrary group, as activity levels are high. Cardiovascular, respiratory and cognitive factors were amongst the candidate biomarkers, alongside factors such as bone density. The results were that, although some factors showed a trend towards a decline with age (as biomarkers are predicted to do) variation was high, and no factor was able to accurately predict age. VO2 max (maximum oxygen consumption) was the best predictor. These results show that sedentary lifestyles aren’t the full explanation for the lack of effectiveness of biomarkers in predicting age. Ageing is a highly individualistic process, with declines not happening consistently across individuals, even when removing the effects of a sedentary lifestyle. There are still genetic differences that this paper can’t account for, and the study is also cross-sectional which causes problems (see explanation later) but this is still a highly interesting result, that casts doubt on the likelihood of biomarkers being able to accurately pinpoint age, and opens up several avenues for future studies.

The paper is very good, and certainly worthy of attention, but its transformation into a justification of cycling as a means of holding back father time is an inaccurate representation of what the paper was all about. The paper isn’t about the merits of cycling to slow ageing, it is well known that an active lifestyle can have positive effects on physiological function, and no attempt has been made to compare cycling to other forms of exercise. The authors themselves state ‘In the absence of clear evidence defining the amount of exercise necessary to combat the negative effects of inactivity, we pragmatically set standards for acceptance into the study’. The main argument of the paper is about ageing still being highly varied and individualistic, even when removing the effects of exercise levels. So, at what stage did this research morph into a praise of cycling as a way of turning back time? The general method of scientific papers reaching the press is through universities releasing press releases, which then form the skeleton of articles in the media. Recent research in health has begun to look at the roll of press releases in reflecting research, and influencing media articles. An article in the BMJ found that over a third of press releases on health research from 20 top UK universities contained exaggerated advice, claims of causation for correlational studies, or unjustifiable applications of animal studies to human health. Tellingly this then passed on to news articles. When exaggeration was present 58-86% of articles featured the same exaggeration, compared to 10-17% of articles containing exaggeration where it wasn’t present in the press release. An editorial by Ben Goldacre makes several suggestions to counter the problem of misleading press releases, including accountability of authors (both press and academic), transparency and feedback.

Miscommunication of science is common in the media, and universities themselves need to take more responsibility to ensure that they aren’t causing this. News outlets respond to the authority of universities, and at present many universities clearly aren’t standing up to that responsibility by ensuring their communication of science is accurate and clear. Many researchers now use social media to communicate their work, and websites such as allow readers to receive news direct from the authors, but the mainstream media still uses press releases to write their stories, and academics and universities need to take greater responsibility. The title of the press release for this article? ‘Exercise allows you to age optimally’ followed by a picture of a bicycle.

A quick side note – cross sectional vs longitudinal studies

The study in question is cross sectional, meaning that when looking at a ‘decline’ with age, we are actually looking at a difference between people of one age compared to people of another age. This is different to longitudinal studies where individuals would be followed throughout the time period, so any decline would be an actual decline within an individual. To illustrate the issues with cross sectional studies consider the following example.


This is hypothetical data describing a real pattern found in a bird species, Leach’s storm petrel. The y axis displays telomere length. Telomeres are sections of DNA that form a protective cap on the end of chromosomes. Each time a cell divides the telomeres get shorter, until they reach a critical point where cell division can no longer occur. In this example it appears that Leach’s storm petrel’s telomeres are getting longer with age. Telomerase is an enzyme which can repair telomeres, as this result would suggest that the action of telomerase and the lengthening of telomeres is important for lifespan in this species. However, this data is cross sectional – meaning that age 1 on the graph is completely different individuals to age 8. An alternative explanation for the results is that as this is telomeres might not be getting longer, there might be large variation in telomere length in youngsters, and something about long telomeres to start with makes you more likely to live to old age, making the average telomere length longer for old birds, and giving rise to the apparent lengthening of telomeres. This hypothesis would suggest that it is having long telomeres to start with that is important, not having ones that get longer. Only a longitudinal study could fully detect what was happening.

Links and references


From the archives; I originally wrote this in March, 2013

This week has seen news of a new plan to control the highly invasive brown tree snake on the island of Guam. The plan involves toxic mice, and has met with fierce opposition from animal rights groups such as PETA.

First it is important to know a bit of background about Guam, and its large snake problem. Guam is a small island around 2000km north of Papua New Guinea, in the western Pacific. As a U.S territory its economy benefits from the location of an American military base, but its main source of income and jobs is the tourism industry.

For the past 70 years or so the people of Guam have shared their island not only with their rich native fauna, but also with an invader of great consequence, the brown tree snake (Boiga irregularis). There are well documented ecological factors that make some species more susceptible to extinction from invasives than others, and unfortunately island species tend to fit these criteria. Since its arrival the arboreal snake has freely expanded its population to its now huge density.

The impacts of this come in two interlinked flavours. Firstly, impacts on native wildlife; the snake has caused the local extinction of almost all native birds and reptiles, including the endemic Guam Flycatcher (now extinct) and the Guam Rail (now being bred in captivity for reintroductions elsewhere). 9 of 11 forest bird species previously found on Guam no longer occur there. As always with invasives, the impacts don’t stop there. Through impacts on the birds (which are important pollinators) native plants are also in decline. Secondly, impacts occur through conflict with people. Numerous power outages are directly caused by the snakes, along with a decline in tourism over and above that expected from purely the negative effects of the global recession.

So what are the U.S Department of Agriculture doing about this problem, which is having undisputed and hugely negative effects on the people of Guam and their wildlife? This is where we revisit the news that has been widely reported this week. The USDA are planning to use helicopters to drop dead mice laced with the chemical acetaminophen and attached to small parachutes, in order to kill the snakes and reduce the population. Now, firstly it pays to get to the root of the story, ideally a publication or press release from the organisation in question. A quick search led me to see that this story isn’t new at all. National Geographic reported in 2010 that the “toxic mice” were being airdropped in Guam. And the USDA have published reports detailing their efforts in controlling the snakes, and have been using acetaminophen since 2001. The primary concern of groups opposing the plans (ie PETA) is that the drug induces coma and death in the snakes, and takes up to 60 hours to have its affect. Other concerns are that airdrops are simply too random, and effects on non-target species are likely to be large.

These are all legitimate concerns, and animal rights groups play an important role in ensuring we aren’t exploiting species, but In this case I believe people need to do more research before they flat out refute something. Gone are the days when nations would take drastic biological decisions without a seconds thought (eg the disastrous use of myxomatosis to control Australia rabbits). The USDA conducts active research into plans to control invasives. A report from 2007 tested various types of parachute for the dead mice to reduce accidental death of non-target species – mainly coconut and hermit crabs. Far from the random and gung-ho approach the media have enjoyed portraying this is, it is a lot more calculated than that.

In terms of effectiveness this tactic appears to have a high potential, the fact that the brown tree snake will eat prey it has not consumed itself is fairly rare in snakes and is a weakness that conservationists need to exploit, in much the same way medical research looks for weaknesses in viruses. The cruelty issue is a tough one, without getting into a debate of the pain feeling abilities of reptiles, it is a cruel thing to do to any animal, and the snakes didn’t ask to be stowaways and transported to Guam. Having said that, I have come to the conclusion that it is a necessary evil. The USDA currently removes around 7000 snakes a year, out of an estimated population of over 1 million. That statistic alone shows that current efforts just aren’t going to be enough to ever make a serious change. There are also concerns that despite the presence of expert snake detecting dogs at all airports and ports on Guam the snakes could spread to other vulnerable islands. The opportunity is here to make a real dent in the snake population, and improve the livelihood of the people of Guam also. Focussing purely on the cruelty is not looking at the real picture, of the knock on effects on other species and people that this species is having. Sometimes conservation has to be cruel.

brown tree snake


Below is a response to a paper published in Science March 2013 about legalizing the trade in Rhino horn.

Original paper can be found here:

Response to Biggs et al (2013) Science 339

Biggs et al put forth the case the now is the time to legalise the trade in rhino horn, as the best option we have left. Whilst the overall argument of proponents of this idea is convincing and provocative, Biggs argument was much reduced, and featured several key omissions. An article in a journal such as this is sure to reach a huge audience, but the opportunity to present a compelling argument for legalisation of a trade that is inherently abhorrent to many people was somewhat missed.

When mentioning current destinations for horn, the authors mentioned Yemeni daggers and traditional Chinese medicine, and yet failed to mention Vietnam. Demand for horn in Yemen is actually decreasing whilst that in Vietnam appears to be driving the ever increasing poaching levels. It is a crucial error to miss out Vietnam, where a politician is said to have started the now widespread idea that rhino horn can cure cancer. The “cure” is now also being marketed and extended to maladies such as hangovers. Failure to mention this growing demand centre is a worrying sign. It is also of great importance when it comes to predicting the implications of legalising the trade. A lack of case studies prevent anyone from truly knowing what will happen if the trade is legalised, and yet this doesn’t stop the authors making spurious claims. The potential of a huge increase in demand is answered by the authors asserting that with stockpiles and harvesting a growing natural population future demand can be met. Whilst current stockpiles are estimated to be sufficient for ~20 years of current demand, we have no knowledge of how demand will change, and claims that natural harvesting can meet an unknown demand are inherently dubious. One key point that the authors neglected is the potential that legalising the trade could give legitimacy to Vietnamese claims that rhino horn is a broad cure for common illness. This could potentially hugely increase demand, beyond that sustained by natural harvesting. The authors back up their idea that demand may not increase by refuting claims that one off sales of elephant ivory to China in 2008 increased demand and/or poaching. Whilst concrete empirical evidence is indeed lacking, it is foolish to think that demand and poaching of elephants has decreased. Another issue with this, is that even if poaching wasn’t to increase, and we don’t know if that is true, demand is the key variable of interest.

Criminal syndicates who transport horn from Africa to Asia often do so under the protection of fraudulent CITES paperwork. One could easily argue that a complete legalisation of the trade will make the illegal trader easier, and not halt it. This may result in illegal trade still finding its way into the market. There is also the current international issue that trade and use of rhino horn is illegal in all destination countries, so any decision of African countries to legalise trade has to occur concurrently with Asian legalisation, something they may be reluctant to agree to. At this point, criminal syndicates in Asia would still benefit from a legal trade.

Biggs et al briefly mention a ‘safeguard’ against legalisation going badly wrong, namely halting the legal trade. A counter argument to this is that the damage done in this scenario would be huge. The potential Asian market for horn is 2 billion people, and the more people exposed to horn, the more who are likely to demand it. For this reason it can be argued that legal trade is largely irreversible without huge negative effects, and therefore needs to be thoroughly and meticulously planned, not rushed into.

I am inclined to agree that, in time, a legal trade is going to be the only card left to play. Killing poachers isn’t an effective conservation action, as the incentives present a seemingly limitless supply of potential poachers. Whilst acknowledging that time isn’t something available in abundance to Rhino, the framework for a legal trade simply isn’t yet in place. Under the values of the precautionary principle we need evidence that legalisation won’t cause an explosion in demand. Populations of Black Rhino are still increasing, and until we have more information on the impacts of a legal trade, and plans of how precisely it will work, the focus is still rightly on preventing poaching, and, of utmost importance, education and awareness in Asia.

© Brent Stirton / Veolia Environnement Wildlife Photographer of the Year 2012

© Brent Stirton / Veolia Environnement Wildlife Photographer of the Year 2012