Population Biology
How accurate is mark-recapture?

In the study of various species, it is often important for biologists to know, or at least estimate, the number of individuals in a given area (eg. the number of bison in Yellowstone State Park, or the number of trout in Crater Lake).  One common method is to actually capture some of the individuals and then release them back into the population with some sort of permanent mark or tag that indicates their capture.  Some time later, a second round of captures are made. Biologists can then estimate the total size of the population based on the number of newly captured individuals that already have tags.
	One basic mathematical assumption of this method is that individual creatures are caught at random, i.e. that no particular creature has a better chance of getting caught than another. In the following, we will use Mathematica to mimic a the random selection of individual creatures (using fish as an example) from the entire population in a particular habitat (using a pond as the example).

Let's create a virtual pond. The pond will have a certain number of fish in it, which will be represented by the parameter Popsize. For now, we will set the number of fish at a small value, say 23.

PopSize  =  23;

There are PopSize fish in the pond, some of which will eventually be tagged.  A good way to keep track of this is to use a Mathematica list, which can be created using the Table function.  The first part of the table function describes the values to include in the list, the second portion is itself a list that determines the length of the table. For example:

Table[ "Blind Mouse", {3}]

As all the fish in the pond are currently untagged, we will generate a list where each element tells the status of a particular fish, in this case "Untagged". 
We will refer to this list as our Population. Notice that there is one element for
each of the PopSize (= 23) fish.

PopSize =  23;
Population  =  Table["Untagged", {PopSize}]

At any time we can look at the status of a particular fish in the Population by giving its location in the list, using an index. For example, the fish associated with list element number one is, as expected, "Untagged".

Population[[1]]

It's time to go fishing! Here we will add another parameter, NumFirstCatch, denoting the number of fish we catch and tag with the first capture-and-release.

NumFirstCatch  =  7;

So seven fish now need to be labeled "Tagged".  We can do this by using
ReplacePart to change the list.  ReplacePart requires three inputs: the list, the new part, and the place in the list to be replaced.

ReplacePart[Population, "Tagged", 3]

ReplacePart actually creates a whole new list. To actually change the Population, we need to set it equal to the result of ReplacePart.

Population  =  ReplacePart[Population, "Tagged", 1]

We want to repeat this process some number of times equal to NumFirstCatch (in this case a total of seven times). We let Mathematica do the work for us by using a loop function. One common loop function takes the form Do[expr,{i,imin,imax}]. The loop uses a counter, in this case represented by the variable i, to count from imin  to imax, each time evaluating expr. So, to tag each of the first NumFirstCatch (= 7) fish in our population we use the counter i inside the ReplacePart statement, which is inside a
Do loop.

PopSize  =  23;
NumFirstCatch  =  7;
Population  =  
        Table["Untagged", {PopSize}];
Do[
Population  =  
          ReplacePart[Population, "Tagged", i], 
{i, 1, 
          NumFirstCatch}];

Notice that our Population now has seven "Tagged" fish.

Population

Time to go fishing again! This time we will plan to catch a few more (twelve) fish, by setting a variable NumSecondCatch.

NumSecondCatch  =  12;

This time, however, we have to be careful about how we catch the fish.  If we assume the first fish we catch is the first fish in the list, we will automatically catch a tagged fish.  The process assumes that the fish are caught at random, so we must use a random process, something equivalent to the roll of a die, to determine which fish we 'catch'.
	Mathematica can generate a wide variety of random numbers.  The basic
command Random[] gives an essentially random number from 0 to 1.

Random[]

With an integer specification, Mathematica will return either 0 or 1, with equal probability.

Random[Integer]

By adding a short list with minimum and maximum values, we can specify the random selection of one integer within a particular range.

Random[Integer, {5, 10}]

To choose a particular fish in our Population, we set an variable,
CaughtThisOne, equal to a random number from 1 to the size of the population (PopSize).

CaughtThisOne  =  Random[Integer, {1, PopSize}]

We can then use CaughtThisOne as an index to see the status of the particular fish we have caught.

Population[[CaughtThisOne]]

We could catch more fish (perhaps as many as NumSecondCatch?) and make a list of indices for fish we catch, by again using our Table function.

\(Table[
	
      CaughtThisOne\  = 
        Random[Integer, {1, 
            PopSize}], 
	{NumSecondCatch}]\)

Now we must be careful! Look closely at the SecondCatch list. Unfortunately, there is a good chance (97%) that you caught the same fish (i.e. picked the same index) more than once! (If not, you can evaluate the command above again...)
	Is this consistent with our physical system? Are we throwing the fish back in as we catch them so we could catch them twice in the same catch? Probably not, so we need to fix our computational model to match the real system.
	
	Notice that everytime we catch a fish the second time, we remove it from the lake (at least temporarily), leaving one fewer fish in the lake. To keep track of this we make two new lists, StillInLake and SecondCatch, The first is the same as Population to start, while the second is empty (as shown by an empty pair of brackets.

StillInLake  =  Population;
SecondCatch  =  {};

We can use the Length function to determine how many elements each list, and therefore how many fish we have caught or are StillInLake.

Length[StillInLake]

We again choose a particular fish, this time from those StillInLake. Notice that if the size of the list changes, so will the range of the Random numbers. We can then Delete the status of particular fish from those StillInLake, and Append it to the SecondCatch!

CaughtThisOne  =  
        Random[Integer, {1, 
            Length[StillInLake]}];
ThisFish  =  
        StillInLake[[CaughtThisOne]];
StillInLake  =  
        Delete[StillInLake, CaughtThisOne];
SecondCatch  =  Append[SecondCatch,  ThisFish];

StillInLake

SecondCatch

Now we simply Do this process over again until we have reached our
target NumSecondCatch. We put parentheses around the sequence to group it together inside a Do statement.

PopSize  =  23;
NumFirstCatch  =  7;
NumSecondCatch  =  12;
Population  =  
        Table["Untagged", {PopSize}];
Do[
Population  =  
          ReplacePart[Population, "Tagged", i], 
{i, 1, 
          NumFirstCatch}];

    
StillInLake  =  Population;
SecondCatch  =  {};
\(\(Do[
\((CaughtThisOne\  = 
	
            Random[Integer, {1, 
                Length[StillInLake]}]; 
ThisFish\  = \ 
            StillInLake[\([CaughtThisOne]\)]; 
StillInLake\  = \ Delete[StillInLake, CaughtThisOne]; 
SecondCatch\  = \ 
            Append[SecondCatch, \ ThisFish];)\), 
{i, 1, 
          NumSecondCatch}];\)\)

SecondCatch

The number of fish CaughtTwice can then be easily determined using
Count, which gives the number of elements in a list matching a certain value.

CaughtTwice  =  Count[SecondCatch, "Tagged"]

If the catch is truly random, probability says that the ratio of fish CaughtTwice to the NumSecondCatch should be the same as the ratio of the NumFirstCatch tagged fish to the 'unknown' total population. Using our simulation, we should be able to make an estimate of the total population, PopEst by solving for it in the ratio
	(CaughtTwice/NumSecondCatch) =  
        (NumFirstCatch/PopEst)

PopEst  =  N[NumFirstCatch*NumSecondCatch/CaughtTwice]

The N[] puts the output in number (not fractional) form. Notice that you can't do the division if CaughtTwice == 0, i.e. you can't make an estimate unless you recapture at least one tagged fish! We can use an If[condition, truepart, falsepart] statement to return the value zero instead of doing the division, and thus avoid a possible error.  The double equals (==) is used to make comparisons.

PopEst  =  
If[(CaughtTwice  ==  
            0),  
   0,  
   N[
          NumFirstCatch*
            NumSecondCatch/CaughtTwice]
   ]

How good is the estimate? We can calculate a precentage error by comparing to PopSize, which we already know. Abs[] is the absolute value function.

PctErr  =  N[Abs[PopSize  -  PopEst]/PopSize*100]

The error will range between 8.7 percent and 265 percent, although values from 8.7 to 30 percent are most likely. Notice, however, that we caught over half of the fish in the pond as part of the recapture, which, although really helpful for our accuracy, is probably not a reasonable sample size!

Putting it all together, we have a small set of commands with the parameters listed first, where we can easily find and adjust them, and where the final statement is a list of outputs, which include the estimate for the population, the actual population size, and the computed percentage error.  Let's see what we get with a population of 200 and catch and recatch both of 50:

PopSize  = 200;
NumFirstCatch  =  50;
NumSecondCatch  =  50;
Population  =  
        Table["Untagged", {PopSize}];
Do[
Population  =  
          ReplacePart[Population, "Tagged", i], 
{i, 1, 
          NumFirstCatch}];
StillInLake  =  Population;
SecondCatch  =  {};
\(\(Do[
\((CaughtThisOne\  = 
	
            Random[Integer, {1, 
                Length[StillInLake]}]; 
ThisFish\  = \ 
            StillInLake[\([CaughtThisOne]\)]; 
StillInLake\  = \ Delete[StillInLake, CaughtThisOne]; 
SecondCatch\  = \ 
            Append[SecondCatch, \ ThisFish];)\), 
{i, 1, 
          NumSecondCatch}];\)\)
CaughtTwice  =  
        Count[SecondCatch, "Tagged"];
PopEst  = 
        If[(CaughtTwice  ==  0),  
          0,  
N[
            NumFirstCatch*
              NumSecondCatch/CaughtTwice]] ;
PctErr  =  
        N[Abs[PopSize  -  PopEst]/
              PopSize*100];
{PopEst, PopSize, PctErr}

If we have a sequence of statements with clear inputs and outputs, it is often useful to define the group as a function. We will call our function GuessPop and require three inputs: the actual population size and the size of each catch. We use the trailing underline to designate a variable as an input parameter, and use the ': =" to tell Mathematica to 'evaluate later' (i.e. when the function is called with given inputs). We group the statements in the function with parentheses,

\(GuessPop[PopSize_, NumFirstCatch_, 
        NumSecondCatch_]\  := \ 
\((
Population\  = \ 
          Table["Untagged", {PopSize}]; 
Do[
Population\  = \ 
            ReplacePart[Population, "Tagged", i], 
{i, 
            1, NumFirstCatch}]; 
StillInLake\  = \ 
          Population; 
SecondCatch\  = \ {}; 
Do[
\((CaughtThisOne\  = 
	
              Random[Integer, {1, 
                  Length[StillInLake]}]; 
ThisFish\  = \ 
              StillInLake[\([CaughtThisOne]\)]; 
StillInLake\  = \ 
              Delete[StillInLake, 
                CaughtThisOne]; 
SecondCatch\  = \ 
              Append[SecondCatch, \ ThisFish];)\), 
{i, 1, 
            NumSecondCatch}]; 
CaughtTwice\  = \ 
          Count[SecondCatch, "Tagged"]; 
PopEst\  = 
          If[\((CaughtTwice\  == \ 0)\), \ 
            0, \ 
N[
              NumFirstCatch*
                NumSecondCatch/
                  CaughtTwice]]\ ; 
PctErr\  = \ 
          N[Abs[PopSize\  - \ PopEst]/
                PopSize*100]; 
{PopEst, PopSize, 
          PctErr}
)\)\)

Now we can put population and catch numbers directly into our new function!

GuessPop[23, 7, 12]

GuessPop[200, 50, 50]

So, how many fish do you need to tag and/or recatch to get a 'good' population guess? Because we now have a model, we can explore what happens relatively easily! Let's start with a medium sized population, say 150 fish, and catch and tag 20 percent, or 30 of those fish. If we don't catch any fish in the second catch, we will catch zero tagged fish, and predict zero fish in the pond, which is 100 percent error.

GuessPop[150, 30, 0]

On the other end, if we catch 150 fish in the second catch (all of them!), we will catch all 30 tagged fish, and our estimate will be perfect (0 percent error).

GuessPop[150, 30, 150]

To see what happens in between, we can simply look at results for second catches from 1 to 149. Just like the Do function, the Table function can also include a counter (often a single letter variable) that can be used to determine each value in the list that Table creates. For example, to generate an '11 times table', we could evaluate

Table[11*i, {i, 12}]

To create a list of population guesses for second catches from 1 to 150, we  can make a Table of the first part ([[1]]) of each of i guesses (remember, GuessPop returns three parts!), where i goes from 1 to PopSize (= 150). We give the list a name (Est2Catch) and we use the semicolon to make sure the list doesn't output (because it is 150 numbers long!). We can make a similar list of the percentage errors, which come from part three ([[3]]) and we label Err2Catch. Note that this might take some computational time!

PopSize  =  150;
Est2Catch  = 
Table[GuessPop[PopSize, 30, 
              i][[1]], {i, PopSize}];
Err2Catch  = 
Table[GuessPop[PopSize, 30, 
              i][[3]], {i, PopSize}];

To visualize our results, we can plot the data from each list using ListPlot. We store each graph by giving it a name. Each plot is a graph of the position in the list on the x-axis (which in this case corresponds to number in the second catch) and the value of list elements on the y-axis (corresponding here to the estimated population for the first graph and the percent error for the second).

PlotEst  =  
        ListPlot[Est2Catch, PlotRange -> {0, 2*PopSize}];

PlotErr  =  
        ListPlot[Err2Catch, PlotRange -> {0, 100}];

What happens to the accuracy of the guess as we catch more fish? At what point do you think you would 'guarantee' being within 20 percent of the actual population? We could make the visualization easier with some lines at the appropriate percentages, which we can create as a function using Plot.

EstHi20 =  
        Plot[(1 + 20/100)*PopSize, {x, 0, 
            PopSize}];
EstLo20  =  
        Plot[(1 - 20/100)*PopSize, {x, 0, 
            PopSize}];
Err20  =  Plot[20, {x, 0, PopSize}];

Show is useful to combine graphics pictures.

Show[PlotEst, EstHi20, EstLo20];
Show[PlotErr, Err20];

Let us try the same study, but this time we will vary the number of fish we tag, and keep constant (at 30) the number in the second catch (switching i and 30 in GuessPop).

PopSize  =  150;
Est1Catch  = 
Table[GuessPop[PopSize, i, 
              30][[1]], {i, PopSize}];
Err1Catch  = 
Table[GuessPop[PopSize, i, 
              30][[3]], {i, PopSize}];
PlotEst  =  
        ListPlot[Est1Catch, 
          PlotRange -> {0, 2*PopSize}];
PlotErr  =  
        ListPlot[Err1Catch, 
          PlotRange -> {0, 100}];
Show[PlotEst, EstHi20, EstLo20];
Show[PlotErr, Err20];

Does the order of the catch appear to make a significant difference in accuracy of the prediction?
	
	 Notice that each data point is a single guess for the population measurement, and can vary quite a bit due to the dependence on probability.
We might get a more consistent picture if we tried, say, 10 catches at every point and averaged the results.  To do so, we will need to define a Mean function, which takes a list (data), adds up the elements (Apply[Plus,...]) and divides by the number of elements (Length[data]). [Note this function is also available in the Mathematica DescriptiveStatistics package].

Mean[data_]  :=  
        N[Apply[Plus, data]/Length[data]];
Mean[{1, 2, 3, 4}]

To average 10 guesses, we make a Table with 10 values from GuessPop and apply Mean.

Mean[Table[GuessPop[150, 30, 5][[1]], {10}]]

If we use these average values in the catch tables instead, we find that the data points don't vary nearly as much.  Note, however, that the simulation takes a fairly long time, as we are making 150*10 = 1500 guesses, each of which requires a run of our GuessPop function! (Note: the following is a second catch simulation).

PopSize  =  150;
Est2Catch  = 
Table[
Mean[
            Table[GuessPop[PopSize, 30, i][[1]], {10}]], {i, 
            PopSize}];
Err2Catch  = 
Table[
Mean[
            Table[GuessPop[PopSize, 30, i][[3]], {10}]], {i, 
            PopSize}];
PlotEst  =  
        ListPlot[Est2Catch, 
          PlotRange -> {0, 2*PopSize}];
PlotErr  =  
        ListPlot[Err2Catch, 
          PlotRange -> {0, 100}];
Show[PlotEst, EstHi20, EstLo20];
Show[PlotErr, Err20];

Now what would you estimate is a good size second catch to guarantee 20 percent or less error? If you only caught 30 fish in the second catch (thus equal to the first catch), what error would you expect from your numerical guess? What do you get when you try a guess?

GuessPop[PopSize, 30, 30]

So, what do you think? Is capture and re-capture a good method of estimating population? Why or why not?

This notebook created by Dr. Garrett R. Love for the Shodor Education Foundation. Copyright 2001