Inclusive fitness errors and divorce in birds

singing-chaffinch-e1568134857283.png
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.

babytits.jpg

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.

Slide1
Slide2

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.

Slide3

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

Figure1_1

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.

v2_Figure2

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.

v2_Figure3

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

v2_Figure4

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

v2_Figure5

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

v2_Figure6

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?

v2_Figure7

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.

v2_Figure8

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.

v2_Figure9

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.

v2_Figure10

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.

worldcupstickers_laurie

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
 }
output
## [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)
library(tictoc)

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)
mean(winnings)
## [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

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

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

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

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
parse(text=(paste("Output[",i,"]*50",sep="")))
## expression(Output[1] * 50)

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

eval(parse(text=(paste("Output[",i,"]*50",sep=""))))
## [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){
belch3("loci",i,"<<-c(i/i)")
}

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){
  belch5("loci.output",i,"<<-matrix(rep(loci",i,"+runif(50),25),nrow=5,ncol=5)")
}

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){
eval(parse(text=(paste("loci.output",i,"<-matrix(rep(loci",i,"+runif(2),25),nrow=5,ncol=5)",sep=""))))
}

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

palette(rainbow(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){
  belch3("allele",i,"<<-matrix(runif(25,1,2),25,5)")
}

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

palette(rainbow(number.loci))  
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')")
}
  }

Simple Evolutionary Model Tutorial in R

Simple Evolutionary Model Tutorial

This tutorial will demonstrate how to construct a simple evolutionary model in R. More specifically, it will show how to set up a population where individuals interact and receive payoffs according to the well known Prisoner’s dilemma, and simulate that population evolving through many generations.

The PDF of this tutorial is available here

The aim of this tutorial is to cover the basics in running evolutionary simulations in R – something that is missing in many university courses and online tutorials. The code used is simple, but demonstrates the basic use of many key functions that are useful in all circumstances. The full code is posted at the end of the tutorial. As always, it is reccomended to use RStudio for ease of code editing, problem solving, plotting etc.

Scenario

The canonical example of the game involves members of a criminal gang, but I prefer to use a more real and recent example of doping in cycling. The key to any game is players, strategies, and payoffs. In this example two cyclists can choose to cooperate by following the doping rules or defect by cheating and taking performance enhancing drugs. The average payoffs are greater if all players cooperate, but rational individuals will always defect, as they receive a higher payoff from doing so, regardless of the choice of the other player. In the cycling example it is better to cheat if the other player is cooperating, as by doing so you can win more races. It is also better to cheat if the other player is cheating, as this is required to still have a chance of competing against the enhanced other player. The outcome of this is a suboptimal scenario for everyone.

eden-presentation

Obviously this is a simplification of the real scenario, but for more detail on the Prisoner’s dilemma in sports doping see: https://www.scientificamerican.com/article/the-doping-game-payoffs/ For anyone with a specific interest in game theory, the open course from Yale is particularly good http://oyc.yale.edu/economics/econ-159

The code

A useful first step in any code is to clear the memory using rm(list=ls()) and then set the number of generations, population size, and mutation rate.

# Clear Memory
rm(list=ls())

# Set number of generations
maxgen <- 20

# Set population size
N <- 1000

# Set the mutation rate
mutation.rate = 0.005

We can then set our payoffs, and visualise them in a dataframe. b is often used for ‘benefit’ of cooperation, and c for cost of cooperation.

# Set the payoffs
b <- 3
c <- 1

the paste function will be helpful for viewing our payoff matrix in the traditional form, with each box containing the payoffs for the two players, seperated by a comma. paste("Payoff =", b-c, sep = " ") will return the text “Payoff =” followed by whatever our value for b-c is, seperated by a space. This function becomes highly useful in more complicated settings. Here we are using it to make two vectors, x and y, which each contain the relevant payoffs for that box, seperated with a comma by the paste function. We can then name the columns C for cooperate and Dfor defect.

# visualise payoffs
x <- c(paste(b-c, b-c, sep = ", "), paste(-c, b, sep = ", "))
y <- c(paste(b, -c, sep = ", "), paste(0, 0, sep = ", "))
df <- data.frame(x,y)
colnames(df) <- c("C", "D")
rownames(df) <- colnames(df)
df
##       C     D
## C  2, 2 3, -1
## D -1, 3  0, 0

Now we can set our starting population, for which we will use all cooperators. Here we introduce the invaluable sample function. this works in the basic form of what to sample, how many, and whether to put the extracted value back in the mix. For example sample(c(1:50), 10, replace = TRUE) will sample from the numbers 1:50 ten times, and put each sample back into the mix, where it may be sampled again. In this example we are using the number 1 to mean a cooperator and -1 to mean a defecter. We can sample from a vector of just a single one (aslong as replace = TRUE) to generate a population of the right size such that all individuals will be 1 and therefore cooperators.

start <- c(1)
population <- sample(start, N, replace = TRUE)

If we wanted to start with an equal mix of cooperators and defectors we could use start <- c(1,-1) and sample from that, or if we wanted to use a certain proprtion we could use the ability of sample to make a biased sample. Here we will assign the probability 0.25 to cooperators, and 0.75 to defectors, and check the freqeuncy of our population.

start <- c(1, -1)
prob <- c(0.25, 0.75)
population <- sample(start, N, replace = TRUE, prob = prob)
length(population[population==1]) / length(population)
## [1] 0.236
length(population[population==-1]) / length(population)
## [1] 0.764

These last two lines illustrate some other simple and useful functions. length(population) will tell us how ‘long’ in individuals that vector is (or matrix row or column etc if we wanted), and can also be adapted to tell us how many of our population meet the condition of equalling one length(population[population==1]) which can be useful in extracting how frequent cooperators or defectors are in our population.

Next, we will set up a matrix to store our output from each generation. It is good practise to assign all the required rows and columns of our matrix at the start, add meaningful column names, and fill it with NA to help identify errors. output <- matrix(NA, maxgen, 4) will give us a matrix of NAs with ‘maxgen’ number of rows and 4 columns, in which we will store generation number, frequency of cooperators, freqeuncy of defectors, and average payoff.

output <- matrix(NA, maxgen, 4)
colnames(output) <- c("generation", "freq coop", "freq defect", "mean payoff")

We are now ready to open our loop. They key is to start with something (here it is our population) originally defined outside the loop, apply selection in terms of payoffs, mutations to bring in some randomness, and end up with a population at the end of the loop. Then we are ready to simply close the loop and feed back through to the start.

Simple Example Loop

A simple example loop could be constructed as follows;

# Begin Loop
for (z in 1:maxgen){
  # make two samples, to be matched as pairs of       individuals
  player1 <- sample(population,N,replace=FALSE)
  player2 <- sample(population,N,replace=FALSE)
  # make a vector for the payoff to cooperators
  coop.payoff <- 0
  # loop through the population, if player one and    player two both equal 1 (i.e. are cooperators)         assign a payoff to cooperation of b-c
  for(y in 1:N) {
   if(player1[y]==1 && player2[y]==1) {
      coop.payoff <- coop.payoff + (2 * (b - c))
   }
  }
  # After this loop, store the generation number and  payoff to cooperators in the output matrix
  output[z, 1] <- z
  output[z, 2] <- coop.payoff
}

Here, z will start at one, and go through the whole loop, and become two. This will continue until z has reached maxgen, our number of generations. Within each loop we are taking two vectors as samples of populations, which can be lined up so the first number of player 1 is partnered with the first number of player 2. We then use a second loop, going from 1 to our population size, where we use an if function to see if both players are 1 (are cooperators) and if so, we give cooperators a payoff. if(player1[y]==1 && player2[y]==1) {coop.payoff <- coop.payoff + 1} will add 1 to coop.payoff only if player1 and player2 at that place in the population (player1[2] will give the second number in player 1) both equal one.

This is obviously not a useful loop as we have not applied selection or mutation, essentially we are just seeing random variation in how many cooperators match up in each generation (and if we have a population of all 1s to begin with, this won’t change, so you can play around with the starting population)

This demonstrates the basic approach of having a loop of each generation, and nested within that having another loop that goes through each individual in the population and assigns payoffs. After the internal loop going through each individual in the payoff, the generation number (z here, as this is the letter chose for the loop of generations) and payoff are stored in the relevant part of the output matrix. The basic approach for doing this is to use something like output[z, 1] <- z which assigns row z column one to z. As we are looping through z generations, on each run through the loop output[z, 1] <- z will refer to the first row, then the second row, then the third row etc. if we wanted to view this, we could simple examine output[, 1] which will show us the first column of the matrix.

Actual Loop

As a loop has to be closed to run without error we will build up the code inside of the loop and then add the loop itself. To start we can take two vectors as samples of the population, which can be lined up so the first number of player 1 is partnered with the first number of player 2. We can then set the starting payoff as 0 for this generation.

# make pairs of individuals
  player1 <- sample(population,N,replace=FALSE)
  player2 <- sample(population,N,replace=FALSE)

# set payoff counters
  coop.payoff <- 0
  defect.payoff <- 0

Now we can prepare the internal loop where we run through the population assigning payoffs. This can be acheived with if functions for all of the possible options. These tend to run quite efficiently in R, and can be setout so if the condition is not met, the program simply moves on to the next if statement.

# Loop to calculate payoffs based on payoff matrix
   for(j in 1:N) {
    # If both players cooperate (=1) payoffs = b-c, b-c
    if(player1[j]==1 && player2[j]==1) {
      coop.payoff <- coop.payoff + (2 * (b - c))
        }
    # If player one cooperates and two defects (1, -1) payoffs = -c, b
    if(player1[j]==1 && player2[j]==-1) {
      coop.payoff <- coop.payoff - c
      defect.payoff <- defect.payoff + b
        }
     # If player two cooperates and one defects (-1, 1) payoffs = b, -c
     if(player1[j]==-1 && player2[j]==1) {
       coop.payoff <- coop.payoff - c
       defect.payoff <- defect.payoff + b
        }
     # If both players defect (-1, -1) payoffs = 0, 0
     if(player1[j]==-1 && player2[j]==-1) {
       coop.payoff <- coop.payoff + 0
       defect.payoff <- defect.payoff + 0
     }
   }

In this example we are using one vector for each strategy (coop.payoff) and (defect.payoff) that we are simply adding to or subtracting from based on the match with the other player. In this sense we are looking at payoffs to strategies rather than individuals. This is because it is much quicker to run the code this way and roughly equivalent. Sometimes it may be required to follow individual payoffs, which can be done by assigning payoff for each individual to a new vector. If this is something you are interested in, then something adapted from the code below could be useful.

player1.payoffs <- c(0)
player2.payoffs <- c(0)
for(j in 1:N) {
    # If both players cooperate (=1) payoffs = b-c, b-c
    if(player1[j]==1 && player2[j]==1) {
      player1.payoffs[j] <- b - c
      player2.payoffs[j] <- b - c
        }
}

Now that we have payoffs, we can apply selection to our population, generating a new population with the frequency of the two strategies based on their relative payoffs. In the below code we first assign any negative payoffs to 0 to allow meaningful relative payoffs to be calculated. Next, we sum the payoffs, and for strategy generate a new population based on relative payoffs. Note that we are also making our new population larger than we need to by multiplying by three. This is because we will soon be selecting a new population of 1000 for the next generation, and want to mimic some of the randomness of selection for the next generation of real scenarios, and to ensure we never end up with a decrease in population size. pop.offspring here uses the rep function, which repeats something a set number of times e.g. rep(1,5) repeats one, five times. pop.offspring therefore becomes a vector of 1s and -1s, with each strategy repeated a number of times according to their relative payoff. by sampling 1000 individuals from this new population, we have our population of selected individuals for the next generation.

# generate new population based on payoffs
  if(coop.payoff < 0){coop.payoff <- 0}
  SUM <- coop.payoff + defect.payoff
  pop.coop <- coop.payoff / SUM * N * 3
  pop.defect <- defect.payoff / SUM * N * 3
  pop.offspring <- c(rep(1,pop.coop), rep(-1,pop.defect))
  population <- sample(pop.offspring,N,replace=FALSE)
  
  # calculate frequency of cooperators and defectors in this population
  fc <- sum(population == 1) / (sum(population == 1) + sum(population == -1))
  fd <- sum(population == -1) / (sum(population == 1) + sum(population == -1))

We have also here calculated the freqeuncy of 1s and -1s in our new population, as this is what we aim to track over time.

The penultimate important step is to mutate our population. This can be done in a number of ways, but this way is the most intuitive. We define m as the number of mutations expected (mutation rate X population size). Then we use sample (it really is very useful) to pick m (here 5) numbers which represent the individuals to mutate (e.g. the 3rd, 454th, 667th, 787th and 990th). individuals at those ‘mutation places’ are then randomly given a value of -1 or 1. As such they may not actually mutate, which is something to remember when interpreting mutation rate effects.

# mutation
  m <- mutation.rate * N
  mutation.place <- sample(length(population), m, replace = FALSE)
  population[mutation.place] <- sample(c(-1, 1), 5, replace = TRUE)

In the final chunk of code we are simply recording generation number, frequency of cooperators, frequency of defectors, and average payoff in our output matrix. The last line will simply print the number of the current generation on the console. Helpful if you are running large numbers of generations and want to know how long you will be waiting!

  output[i, 1] <- i
  output[i, 2] <- fc
  output[i, 3] <- fd
  output[i, 4] <- SUM / N / 2
  print(paste("Generation =", i, sep = " "))

The loop can now be closed, so the full loop will be as follows;

# Begin Loop
for (i in 1:maxgen){
 
# make pairs of individuals
  player1 <- sample(population,N,replace=FALSE)
  player2 <- sample(population,N,replace=FALSE)

# set payoff counters
  coop.payoff <- 0
  defect.payoff <- 0
  
# Loop to calculate payoffs based on payoff matrix
   for(j in 1:N) {
    # If both players cooperate (=1) payoffs = b-c, b-c
    if(player1[j]==1 && player2[j]==1) {
      coop.payoff <- coop.payoff + (2 * (b - c))
        }
    # If player one cooperates and two defects (1, -1) payoffs = -c, b
    if(player1[j]==1 && player2[j]==-1) {
      coop.payoff <- coop.payoff - c
      defect.payoff <- defect.payoff + b
        }
     # If player two cooperates and one defects (-1, 1) payoffs = b, -c
     if(player1[j]==-1 && player2[j]==1) {
       coop.payoff <- coop.payoff - c
       defect.payoff <- defect.payoff + b
        }
     # If both players defect (-1, -1) payoffs = 0, 0
     if(player1[j]==-1 && player2[j]==-1) {
       coop.payoff <- coop.payoff + 0
       defect.payoff <- defect.payoff + 0
     }
   }
  
  # generate new population based on payoffs
  if(coop.payoff < 0){coop.payoff <- 0}
  SUM <- coop.payoff + defect.payoff
  pop.coop <- coop.payoff / SUM * N * 3
  pop.defect <- defect.payoff / SUM * N * 3
  pop.offspring <- c(rep(1,pop.coop), rep(-1,pop.defect))
  population <- sample(pop.offspring,N,replace=FALSE)
  
  # calculate frequency of cooperators and defectors in this population
  fc <- sum(population == 1) / (sum(population == 1) + sum(population == -1))
  fd <- sum(population == -1) / (sum(population == 1) + sum(population == -1))
  
  # mutation
  m <- mutation.rate * N
  mutation.place <- sample(length(population), m, replace = FALSE)
  population[mutation.place] <- sample(c(-1, 1), 5, replace = TRUE)
  
  output[i, 1] <- i
  output[i, 2] <- fc
  output[i, 3] <- fd
  output[i, 4] <- SUM / N / 2
  
  print(paste("Generation =", i, sep = " "))
 }

Plotting Figures

It is easy to visualise our results using simple plots. The online help for basic plotting is very comprehensive so I won’t bother to explain too much. Plot 1 is generation number against frequency of each strategy (cooperators in blue, defectors in red). Plot 2 is generation number against average payoff, with dotted lines for the individual options of b, b-c, 0, and -c

par(mfrow=c(2,1))
par(mar=c(4.5,4.1,1.1,1.1))
# blue = cooperators
plot(output[,1], output[,2], type = "l",lwd=2, col="blue", xlab = "generation", ylab = "frequency", ylim=c(0,1))
# red = defectors
points(output[,1], output[,3], type = "l",lwd=2, col="red")
plot(output[,1], output[,4], type = "l",lwd=2, col="black", xlab = "generation", ylab = "average payoff", ylim=c(-1,3))
abline(h = 0, lty = 3)
abline(h = (b), lty = 3, col= "green")
abline(h = (-c), lty = 3, col="red")
abline(h = (b-c), lty = 3, col="blue")

pdplot

As we can see, it doesn’t take many generations for defectors to invade the population, and everyone suffers as a result.

And now we have our model. A population evolving over time, with both selection and mutation, towards the inevitable ESS outcome of defectors reaching fixation.

The full code for this model can be found below. It is easily adaptable in some ways, the benefits and costs can be changed, but if we want to make the game totally different, for example into a snowdrift game, there are some alterations that can be made. These can be found below the full script

###### Simple Prisoner's Dilemma

# Laurie Belcher 07/11/2016

# Clear Memory
rm(list=ls())

# Set number of generations
maxgen <- 20

# Set population size
N <- 1000

# Set the mutation rate
mutation.rate = 0.005

# Set the payoffs
b <- 3
c <- 1

# visualise payoffs
paste("Payoff =", b-c, sep = " ")
x <- c(paste(b-c, b-c, sep = ", "), paste(-c, b, sep = ", "))
y <- c(paste(b, -c, sep = ", "), paste(0, 0, sep = ", "))
df <- data.frame(x,y)
colnames(df) <- c("C", "D")
rownames(df) <- colnames(df)
df

# set up starting population
# -1 = will defect
# 1 = will cooperate
start <- c(1, -1)
population <- sample(start, N, replace = TRUE)
prob <- c(1, 0)
population <- sample(start, N, replace = TRUE, prob = prob)
length(population[population==1]) / length(population)
length(population[population==-1]) / length(population)


# set up matrix to store output
output <- matrix(NA, maxgen, 4)
colnames(output) <- c("generation", "freq coop", "freq defect", "mean payoff")

# Begin Loop
for (i in 1:maxgen){
 
# make pairs of individuals
  player1 <- sample(population,N,replace=FALSE)
  player2 <- sample(population,N,replace=FALSE)

# set payoff counters
  coop.payoff <- 0
  defect.payoff <- 0
  
# Loop to calculate payoffs based on payoff matrix
   for(j in 1:N) {
    # If both players cooperate (=1) payoffs = b-c, b-c
    if(player1[j]==1 && player2[j]==1) {
      coop.payoff <- coop.payoff + (2 * (b - c))
        }
    # If player one cooperates and two defects (1, -1) payoffs = -c, b
    if(player1[j]==1 && player2[j]==-1) {
      coop.payoff <- coop.payoff - c
      defect.payoff <- defect.payoff + b
        }
     # If player two cooperates and one defects (-1, 1) payoffs = b, -c
     if(player1[j]==-1 && player2[j]==1) {
       coop.payoff <- coop.payoff - c
       defect.payoff <- defect.payoff + b
        }
     # If both players defect (-1, -1) payoffs = 0, 0
     if(player1[j]==-1 && player2[j]==-1) {
       coop.payoff <- coop.payoff + 0
       defect.payoff <- defect.payoff + 0
     }
   }
  
  # generate new population based on payoffs
  if(coop.payoff < 0){coop.payoff <- 0}
  SUM <- coop.payoff + defect.payoff
  pop.coop <- coop.payoff / SUM * N * 3
  pop.defect <- defect.payoff / SUM * N * 3
  pop.offspring <- c(rep(1,pop.coop), rep(-1,pop.defect))
  population <- sample(pop.offspring,N,replace=FALSE)
  
  # calculate frequency of cooperators and defectors in this population
  fc <- sum(population == 1) / (sum(population == 1) + sum(population == -1))
  fd <- sum(population == -1) / (sum(population == 1) + sum(population == -1))
  
  # mutation
  m <- mutation.rate * N
  mutation.place <- sample(length(population), m, replace = FALSE)
  population[mutation.place] <- sample(c(-1, 1), 5, replace = TRUE)
  
  output[i, 1] <- i
  output[i, 2] <- fc
  output[i, 3] <- fd
  output[i, 4] <- SUM / N / 2
  
  print(paste("Generation =", i, sep = " "))
 }
par(mfrow=c(2,1))
par(mar=c(4.5,4.1,1.1,1.1))
# blue = cooperators
plot(output[,1], output[,2], type = "l",lwd=2, col="blue", xlab = "generation", ylab = "frequency", ylim=c(0,1))
# red = defectors
points(output[,1], output[,3], type = "l",lwd=2, col="red")
plot(output[,1], output[,4], type = "l",lwd=2, col="black", xlab = "generation", ylab = "average payoff", ylim=c(-1,3))
abline(h = 0, lty = 3)
abline(h = (b), lty = 3, col= "green")
abline(h = (-c), lty = 3, col="red")
abline(h = (b-c), lty = 3, col="blue")

Making the code more adpatable to different games

Let’s make the payoff matrix into an actual matrix. With each of the four cells referring to the playoff to player 1 (who choses between row 1 and row 2)

payoffs <- matrix(0, 2, 2)
payoffs[1,1] <- b - c
payoffs[1,2] <- -c
payoffs[2,1] <- b
payoffs[2,2] <- 0
payoffs
##      [,1] [,2]
## [1,]    2   -1
## [2,]    3    0

This is the simple prisoner’s dilemma that we have been using for the whole tutorial. Now that we are using a matrix for payoffs instead of only b and c we need to change the payoffs loop to account for this

for(j in 1:N) {
    # If both players cooperate (=1) payoffs = b-c, b-c
    if(player1[j]==1 && player2[j]==1) {
      coop.payoff <- coop.payoff + (2 * payoffs[1,1])
    }
    # If player one cooperates and two defects (1, -1) payoffs = 0, b
    if(player1[j]==1 && player2[j]==-1) {
      coop.payoff <- coop.payoff + payoffs[1,2]
      defect.payoff <- defect.payoff + payoffs[2,1]
    }
    # If player two cooperates and one defects (-1, 1) payoffs = b, 0
    if(player1[j]==-1 && player2[j]==1) {
      coop.payoff <- coop.payoff + payoffs[1,2]
      defect.payoff <- defect.payoff + payoffs[2,1]
    }
    # If both players defect (-1, -1) payoffs = -1, -1
    if(player1[j]==-1 && player2[j]==-1) {
      defect.payoff <- defect.payoff + payoffs[2,2] + payoffs[2,2]
    }
  }

The advantage of this approach is that we can easily alter the game. For example, we can change the game to a snowdrift game. The snowdrift game differs from the prisoner’s dilemma in that there is a higher payoff to cooperating with a defector than from mutual defection. This is the snowdrift payoff matrix.

payoffs <- matrix(0, 2, 2)
payoffs[1,1] <- b - c
payoffs[1,2] <- 0
payoffs[2,1] <- b
payoffs[2,2] <- -c
payoffs
##      [,1] [,2]
## [1,]    2    0
## [2,]    3   -1

If we run our model with this payoff matrix, we find a different outcome

sdplot

A coexistence between cooperators and defectors occurs in this scenario (though we may wish to increase the number of generations to check that this is stable). There are several other ways that the payoffs can be structured and ways that this basic model can be extended.

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 theconversation.com 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.

telomeres

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

http://www.badscience.net/2014/12/my-bmj-editorial-on-academic-press-releases-routinely-misleading-the-public/

http://mbe.oxfordjournals.org/content/25/1/220.full

http://www.kcl.ac.uk/newsevents/news/newsrecords/2015/January/Exercise-allows-you-to-age-optimally.aspx

http://www.independent.co.uk/sport/cycling/the-secret-of-eternal-youth-skintight-lycra-and-a-bicycle-9959058.html

https://au.lifestyle.yahoo.com/prevention/a/25920895/study-cycling-slows-ageing-process/

WHEN CONSERVATION GETS CRUEL

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. http://www.guardian.co.uk/world/2013/feb/24/us-guam-snakes-mice-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