# Football Crazy, Probability Mad: How much does it really cost to complete the Euro 2020 Sticker album?

Every four years football gives us the European Championship. 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 tournament. In doing so, they give us an intriguing mathematical challenge: how can we collect all 678 stickers as cheaply as possible?

This year there will be, as there always is, 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 678 stickers to collect
2) Each pack contains 5 stickers and costs £0.90
3) Each sticker is equally likely to occur in any pack
4) Individual stickers can be purchased on the Panini website for £0.28 each

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

You will see articles claiming that it costs £840 to complete the collection. This is wrong, as it doesn’t consider the possibility of buying individual missing stickers from the Panini website.

For those who just want to know how to complete the collection as cheaply as possible – buy 60 packs, and then buy all the missing stickers online. Total cost: £176. If you have a friend to swap duplicates with, you should each buy 85 packs, do your swaps, and then buy the remainder online. Total cost each: £164. If you have more than one friend in your swapping syndicate, then the price comes down again slightly, but by a smaller amount for each new friend you include.

If you care how to calculate this: read on!

To start – let’s think about the simplest way that someone can collect all the stickers. You could 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 678 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 678 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 10,000 times and have plotted the results below

I found that on average you would have to buy 934 packs to complete the collection! Which would cost a total of £840.60, 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 934 packs would give us 4670 stickers, of which only a meagre 15% 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 678 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 434 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 £800 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.28 each online, compared to the £0.18 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 934 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 60 packs, for a total cost of £175.93. A much more reasonable answer than 934 packs for a cost of £840

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/678. 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 677/678, 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 676/678. 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 could have arrived at the same answer

So is that the final word, shall we put some of our savings back in the bank and prepare to spend £176? 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!

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 blue 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). We have to buy more packs when playing with a partner than we do when alone – 85 packs compared to 60 when playing on our own. The key however is that the total cost to complete the collection is lower; £164. 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.28 each.

Adding more friends marginally lowers the cost still, but not by much. If you think about it, there might be extra swaps if you have two friends with you, but those swaps also have to be shared amongst more people. You do get a benefit in avoiding a few of those annoying cases where your friends’ swaps are of no use to you, but the benefit is small. To see more, check out my 2018 post on the World Cup sticker collection

https://goodscienceblog.wordpress.com/2018/03/12/how-much-does-it-costs-to-complete-the-world-cup-2018-sticker-album/

If you want to see how your progress compares with the expectation (to see how lucky you are), check out the excel spreadsheet below

# 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.

worldcupstickers_laurie

# St Petersburg Paradox – Code Efficiency

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.

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

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

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.

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.

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.

# RHINO HORN: LEGAL TRADE NOT AN OPTION – YET

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: http://www.sciencemag.org/content/339/6123/1038.summary

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

# FIFA World Cup – Caps and Goals: why some players get more than others

Football culminates in the FIFA World Cup, the premier showcase of global talent. We investigated factors influencing the number of international caps and goal efficiency of players selected for the 2014 World Cup, focussing on differences between continents. We found that goal efficiency is affected by whether players play in their home country, and age but these relationships differ markedly between continents. We also found that international caps are affected by continent, and the interaction between age and whether players play in their home country. We recommend further analysis be conducted on larger longitudinal data to determine if these relationships hold for the careers of individual players.

Introduction

With the FIFA World Cup now in full swing fans all over the world are debating who should start for their country, and betting on who will be top goalscorer. Here, we have taken a more analytical approach, using data on all thirty-two national squads to determine what really impacts number of caps and goals of international players, and whether these relationships differ between players from different continents. In the U.S sport as a whole has embraced statistics as a way of enriching the discussions of fans and boosting the success of clubs. European clubs are beginning to catch up, but the statistical revolution here has yet to hit full stride. Delving into data such as this could provide insight into a number of areas. An example comes from suspected age fraud in African players, with nations such as Nigeria unable to match their numerous successes at under-17 level in the senior game. Recent improvements in ageing testing at the under-17 world cup led to a large withdrawal of Nigerian players in 2009, and retrospective testing suggests as many as 35% of players were overage in the three previous tournaments. Rumours have surfaced that professional African players may be as many as ten years older than their official age. Using a large database it may be possible to use the relationship between age and goals for players of a known age to compare to the same relationship for players suspected to be older than they claim to determine their actual age. In recent years some of the traditional dogma of football has begun to be challenged, with work such as the book Soccernomics using data to test numerous assertions, such as racial discrimination in the U.K in the 80’s, and the true age at which players peak. The author’s analysis does an excellent job of highlighting the immense value of information such as this to clubs. Regarding age of peak performance, if a club knows which factors influence this, they can make tailored predictions of when each players will reach his peak, and possibly sell players for large fees when they start their statistical decline, but before other teams notice this decline. The two variables we have chosen to focus on is number of caps and number of goals, with a separate measure created for goals per cap to account for some player scoring lots of goals purely because they have played lots of games. Age is of course a factor we would expect to have a considerable impact on both number of caps and number of goals, but we will also be investigating how this relationship changes depending on which continent a player is from, his height, and whether he plays in his home country. With the top European leagues heavily overrepresented in the World Cup it is inevitable that many countries will rely on players that don’t play their club football in their home country. Debate continues in many countries on the effect this may have, and whether in some countries these players are disadvantaged due to playing football in possibly a different style to other members of the team. This factor may also interact with factors such as age in determining the number of caps a player gets. National coaches may be under pressure to pick young players who have performed well in their domestic league, rather than ageing players playing abroad. Our chief aim here is to assess whether any of these effects exist, and whether the relationships can be quantified to allow predictions to be made about the future caps and goals of players based on factors such as age, height and where they play their club football.

Methods We have used the database of all 736 players originally selected for the 23 player world cup squad of one of the 32 competing nations.

Statistical Analysis

Descriptive Statistics To begin we have graphed summary statistics on continental differences between key factors of the database, which has revealed some interesting differences.

What factors affect goal efficiency of players? To investigate the factors affecting goal efficiency we decide to look only at players listed as forwards, as a substantial number of players have scored very few international goals due to age or position, skewing the data. We created a general linear model (GLM) with goals per game (natural log transformed to normalise) as our response variable and explanatory variables age, continent, whether a player plays in his home country, and height.

What factors affect the number of caps a player gets? To investigate factors affecting the number of caps a player gets we constructed a quasipoisson GLM, using all players with at least 1 international cap. We fitted number of caps as the response variable alongside explanatory variables age, continent, height and whether a player plays in his home country.

Results Descriptive Statistics One interesting factor we investigated is height. At a national level there was considerable variation, with Ghana’s squad standing a diminutive 173.6cm tall on average, compared to the average German squad member at 185.2cm tall. Height is a trait with considerable global variation, impacted by numerous factors. We compared average height of a nation’s squad with average height of an adult male from that nation and found some striking differences.

Figure 1.0 – Height difference between average member of the 2014 World Cup squad of a nation and average adult male from that nation. Data shown for the five smallest and five largest differences.

African and Australian squads picked the youngest sides for the tournament (Figure 2.0) with Ghana, Nigeria and Australia all picking teams with an average age of under 26. Squads from the Americas were much older on average, with Argentina boasting the highest average age of 28.5 years. There was substantial goal efficiency variation within continents, translating to small differences in averages between the continents, with South America achieving the highest goal efficiency (Figure 3.0) alongside highest average age (Figure 2.0). Brazil, Mexico and Germany have the highest goal efficiency of all nations, with Greece, Costa Rica and Italy having the smallest.

There was much more variation in number of caps between continents, although on the whole continents with older players on average also unsurprisingly had players with more international experience on average. It is important to note at this stage that the continent of Australia is only represented by the nation of Australia, whereas there are 13 countries representing Europe so there are some issues in comparing these continents. The most experienced nations are Spain, Uruguay, Honduras and Germany all with an average of over 40 caps for each player. At the other end of the spectrum is Algeria, Australia and France with under 25 caps on average for each player.

Average number of goals was highly variable, and as a result it is difficult to suggest substantial differences between continents. At a national level, Germany and Spain we big leaders in average goals, with the average player from that nation having over 9 goals. Algeria, Colombia and Nigeria brought the squads with the fewest average goals.

What factors affect goal efficiency of players? Goal efficiency was significantly influenced by an interaction between home continent and player age (F5,140=2.4199, p=0.03871) and an interaction between continent and whether they played for their home country (F4,140=3.6123, p=0.007796)

What factors affect the number of caps a player gets? Number of caps was found to be significantly influenced by continent (Quasipoisson GLM, F1,627= 4.3008, p=0.0007445) and an interaction between player age and whether they played for their home country or not (F1,627 =10.579, p=0.001205). There was a trend for increasing number of caps with increasing height, but this relationship was not significant (F1,626 =3.3297,p=0.06852).

Does it matter if you play in your home country? In terms of number of international caps we found that for players playing their home country, number of caps increases slightly with age, but for players not playing in their home country number of caps massively increases with age. Again it is important to note that this data isn’t following the ageing of individual players, rather it is a snapshot in time of players of different ages. Nonetheless, it still shows that playing away from your home country can do wonders for your international career, and unsurprising trend that is probably driven by non-European nations due to the likes of Spain, England and Germany having players from their home country playing in the best leagues in the world. In this model number of caps was also significantly influenced by continent, unsurprising given that we know some continents brought more experienced teams than others (Figure 4). More interestingly there was a trend that was not quite significant of number of caps increasing with increasing height. More work on this is needed to confirm this relationship, particularly by controlling for position as goalkeepers and defenders may be taller and have longer international lifespan for reasons other than their height. Figure 1 shows the five lowest and five greatest differences in height between the nation and its football stars. Interestingly the Netherlands is the only country where footballers were smaller on average than males from that nation, although as arguably the tallest nation in the world that is somewhat unsurprising. At the other end of the spectrum both Ecuador and Nigeria named squads where players were on average over 15cm taller than the average adult male. With the heritability of human height at around 0.6-0.8 perhaps genetic factors are more important in becoming an international footballer in countries such as Nigerian, Ecuador and Iran compared to nations such as Spain, Ghana and Germany, where it takes much more that height to get spotted.

Do South America and Asian players get less efficient with age? Figure 6.0 shows that older players from South America and Asia show a decreased efficiency in terms of goals per game. It is important to first note that this is not data following the ageing of individual players, rather it is a cross section of efficiency of players in the squad as a whole. Furthermore this analysis was performed using only strikers, so cannot be explained by any effect of these continents tending to pick younger strikers and older defenders. Several possibilities remain. Firstly, young strikers may only be considered for selection if they have proved themselves efficient, whereas older players may still make the squad based on previous form, experience and leadership. Secondly players from South America and Asia may indeed peak at a young age. This could be easily tested using data from previous international or domestic competitions and investigating the change in performance of individual players with age. Perhaps the style of football utilized by South American teams favours younger players, or perhaps Africa has been unfairly scrutinized when investigating age fraud in young players. Further analysis is definitely warranted here.

Does it matter if you play in your home country? For goal efficiency it may depend on continent. We found a significant interaction between continent and home country affecting goal efficiency. In Europe and Asia playing in your home country led to decreased goals per game, whereas in the other continents it led to increased goals per game. This is in contrast to the general view that players in many continents benefit from playing abroad in the best European leagues, as they may actually have higher goal efficiency if they play in their home country. This could potentially be influenced by small sample sizes of some players, and a larger database is needed to see if this effect is real, and players perform better at international level if they play in their home country. In terms of number of international caps we found that for players playing their home country, number of caps increases slightly with age, but for players not playing in their home country number of caps massively increases with age. Again it is important to note that this data isn’t following the ageing of individual players, rather it is a snapshot in time of players of different ages. Nonetheless, it still shows that playing away from your home country can do wonders for your international career, and unsurprising trend that is probably driven by non-European nations due to the likes of Spain, England and Germany having players from their home country playing in the best leagues in the world. In this model number of caps was also significantly influenced by continent, unsurprising given that we know some continents brought more experienced teams than others (Figure 4). More interestingly there was a trend that was not quite significant of number of caps increasing with increasing height. More work on this is needed to confirm this relationship, particularly by controlling for position as goalkeepers and defenders may be taller and have longer international lifespan for reasons other than their height. In this report we have highlighted several intriguing statistical relationships, particularly related to how the factors affecting goal efficiency and number of international caps vary depending on which continent the players are from. These conclusions are based on a relatively small sample of the 736 players selected for the world cup, unevenly split up into continents. In the information age much larger datasets exist, on a much deeper array of football stats. This relatively untapped data source could enhance comparison of player ability and improve efficiency in the transfer market.