Overview

Today:

• A capture–recapture problem
• Fitting a model “by eye” with simulation
• Where we are going this semester

Estimating a population size

Suppose you are a fisheries manager responsible for trout in a small lake.

In order to set limits on how many can be caught, you need to know:

How many trout are living in this lake?

You cannot do a census, but you can catch fish, mark them, and look for them in a new sample later. (Don’t click on the link if you’re squeamish, pictures of tagged fish.)

Marking

You begin by catching 60 trout.

You give each one a small numbered tag and release them back into the lake.

Recapturing and marking

Once a week for three weeks you return to the lake.

Each time, you catch 40 trout.

For every fish you catch:

• If it already has a tag, you record its ID.
• If it does not have a tag, you give it a new tag.

Then you release all 40 fish back into the lake (the number of marked fish in the lake increases over time).

How should we estimate the number of fish?

Suppose that we don’t really know any math or statistics, but we do have a computer and know how to use it.

We have some observed data, for instance, something like:

and we want to know which population sizes could plausibly have produced it.

How might we proceed?

Simulate fake datasets

For any proposed population size \(N\), we can simulate:

1. catching and tagging 60 fish,
2. returning a week later and catching 40,
3. marking the new fish,
4. repeating for three weeks.

Values of \(N\) that tend to give simulated data that look like the real data are plausible candidates, while values of \(N\) that tend to give simulated data that don’t look like the real data are implausible candidates.

The simulator

We might simulate the experiment something like this:

simulate_experiment <- function(N_pop, n_marked = 60, n_recaptured = 40) {
    # fish are represented internally by numbers 1,...,N_pop
    tag_id <- rep(NA_integer_, N_pop)
    # first visit: catch n_marked fish
    first_catch <- sample(1:N_pop, n_marked)
    # give them tag IDs 1,...,n_marked
    tag_id[first_catch] <- 1:n_marked
    next_tag_id <- n_marked + 1
    capture_history <- data.frame(
        week = 0,
        tag_id = 1:n_marked,
        already_tagged = FALSE
    )
    for (week in 1:3) {
        # catch n_recaptured fish this week
        caught <- sample(1:N_pop, n_recaptured)
        # which ones already have tags?
        was_tagged <- !is.na(tag_id[caught])
        # give new tag IDs to previously untagged fish
        new_fish <- caught[!was_tagged]
        if (length(new_fish) > 0) {
            new_ids <- next_tag_id:(next_tag_id + length(new_fish) - 1)
            tag_id[new_fish] <- new_ids
            next_tag_id <- next_tag_id + length(new_fish)
        }
        # record the IDs each fish has after being processed
        capture_history <- rbind(
            capture_history,
            data.frame(
                week = week,
                tag_id = tag_id[caught],
                already_tagged = was_tagged
            )
        )
    }
    return(capture_history)
}

For example:

fish_data <- simulate_experiment(N_pop = 200)
fish_data |> head()

##   week tag_id already_tagged
## 1    0      1          FALSE
## 2    0      2          FALSE
## 3    0      3          FALSE
## 4    0      4          FALSE
## 5    0      5          FALSE
## 6    0      6          FALSE

And we can summarise how many marked fish were captured each week.

fish_data <- simulate_experiment(N_pop = 200)
## how many recaptures we get each week
fish_data |> group_by(week) |> summarise(sum(already_tagged))

## # A tibble: 4 × 2
##    week `sum(already_tagged)`
##   <dbl>                 <int>
## 1     0                     0
## 2     1                    11
## 3     2                    15
## 4     3                    23

Run it again and you’ll get something a little different.

Let’s see whether this works

I generated an observed dataset (observed_data). A summary, showing the number of captured fish which had already been marked each week, is below.

# the observed_data object was generated at the top of the script in a code block that I don't show on the slides
# number recaught each week
observed_data |> group_by(week) |> summarise(sum(already_tagged))

## # A tibble: 4 × 2
##    week `sum(already_tagged)`
##   <dbl>                 <int>
## 1     0                     0
## 2     1                    24
## 3     2                    29
## 4     3                    26

Let’s use the simulator to decide what value of N we think generated them.

What did we actually do while we were experimenting?

• Look for simulations whose average was close to the observed average?
• Look at all three numbers at once?
• Try to reproduce the observed values exactly?
• Run many simulations at each value of N?
• Mostly change N until things “looked about right”?

Compare to a standard estimator

You don’t have to know about this estimator, it’s not going to be important for the class.

It is moment-based and apparently fairly standard in this domain.

schnabel_from_fish_data <- function(fish_data) {
    weeks <- sort(unique(fish_data$week))
    catch <- fish_data |> group_by(week) |> summarise(catch = n()) |> pull(catch)
    recaps <- fish_data |> group_by(week) |> summarise(recaps = sum(already_tagged)) |> pull(recaps)
    newmarks <- catch - recaps
    fishmethods::schnabel(
        catch = catch,
        recaps = recaps,
        newmarks = newmarks
    )
}
schnabel_from_fish_data(observed_data)

##                             N        invSE      LCI      UCI CI_Distribution
## Schnabel             112.9114 0.0009964344 83.14210 175.8890               t
## Schumacher-Eschmeyer 114.7646 0.0007979571 82.32609 189.3885               t

So did we get the right answer?

Look at the data more carefully.

# in a simulated dataset, how many fish are caught once, twice, three times, etc?
simulate_experiment(N_pop = 125)$tag_id |> table() |> table()

## 
##  1  2  3  4 
## 43 36 19  2

# in the observed data, how many fish are caught once, twice, three times, four times?
observed_data$tag_id |> table() |> table()

## 
##  1  2  3  4 
## 63 11 13 14

Number of fish caught 1, 2, 3, 4 times in simulated (black) and observed (red) datasets.

What’s going on?

I actually simulated from a model with 200 individuals where 15% of the population is easy to catch and the remainder is hard to catch.

Code is in the un-displayed code block at the top of the Rmd file.

What did we just do?

We had

• some observed data,
• a computer model that could generate fake data for any proposed \(N\),
• and a parameter we wanted to learn about.

Without knowing much statistics, we were able to get a reasonable estimate of \(N\) and a reasonable quantification of the associated uncertainty.

The basic idea was:

Parameter values are plausible if they generate data that look like the data we observed.

This idea gets us a long way.

Can we get a computer to do the same thing?

What we did by eye was something like

\[ \text{try }N \longrightarrow \text{simulate data} \longrightarrow \text{compare to observed data} \longrightarrow \text{try another }N. \]

Could we specify this procedure precisely enough that a computer could do it?

If we can, then we can also ask statistical questions about the procedure:

• Does it recover the right parameter?
• How variable is the estimate?
• When does it fail?

To do so, we will need to pin down exactly what we are comparing, how close various aspects of the data need to be, and so on.

If we can train a computer to estimate parameters, can we then quantify uncertainty in those estimates?

A single estimate like \(\hat N = 113\) does not tell us how much the data actually tell us about \(N\).

We would also like to know:

Which other values of \(N\) are reasonably compatible with the data?

We will look into how simulation can give us

• standard errors
• confidence intervals
• distributions over plausible parameter values

Can we diagnose whether the model is any good?

The fish example showed that we can match some aspects of the data and still get a poor estimate due to model mis-specification.

We will use simulation to diagnose model quality.

Where are we going?

Generative models
Defining generative models, forward simulation, prior predictive checks

↓

Make fitting-by-eye systematic
Simulated method of moments, indirect inference + associated uncertainty quantification, simulation-based optimization

↓

Bayesian version of simulation-based inference Conditioning as rejection sampling, approximate Bayesian computation

↓

Non-parametric generative models and their uses Bootstrap, permutation tests

↓

Try to be more efficient Importance sampling, ABC-SMC, synthetic likelihood, simulated likelihood,
neural posterior/density estimation, neural ratio estimation